36 #include "./base/base_uses.f90"
41 LOGICAL :: has_nonlocal = .false.
42 INTEGER :: n_nonlocal = 0
44 LOGICAL,
DIMENSION(0:5) :: is_nonlocal = .false.
45 REAL(kind=
dp),
DIMENSION(:),
POINTER :: a_nonlocal => null()
46 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: h_nonlocal => null()
47 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: c_nonlocal => null()
48 LOGICAL :: has_local = .false.
49 INTEGER :: n_local = 0
50 REAL(kind=
dp) :: zval = 0.0_dp
51 REAL(kind=
dp) :: ac_local = 0.0_dp
52 REAL(kind=
dp),
DIMENSION(:),
POINTER :: a_local => null()
53 REAL(kind=
dp),
DIMENSION(:),
POINTER :: c_local => null()
54 LOGICAL :: has_nlcc = .false.
56 REAL(kind=
dp),
DIMENSION(:),
POINTER :: a_nlcc => null()
57 REAL(kind=
dp),
DIMENSION(:),
POINTER :: c_nlcc => null()
64 CHARACTER(len=*),
PARAMETER,
PRIVATE :: modulen =
'atom_sgp'
81 TYPE(atom_ecppot_type),
OPTIONAL :: ecp_pot
83 TYPE(gto_basis_set_type),
OPTIONAL,
POINTER :: orb_basis
84 REAL(kind=
dp),
DIMENSION(6) :: error
87 INTEGER,
DIMENSION(3) :: mloc
88 LOGICAL :: is_ecp, is_upf
89 REAL(kind=
dp) :: errcc, rcut
90 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cgauss, cutpots, cutpotu
91 TYPE(atom_basis_type) :: basis
92 TYPE(opmat_type),
POINTER :: core, hnl, score, shnl
95 IF (
PRESENT(orb_basis))
THEN
102 IF (
PRESENT(ecp_pot)) is_ecp = .true.
104 IF (
PRESENT(upf_pot)) is_upf = .true.
105 cpassert(.NOT. (is_ecp .AND. is_upf))
110 ALLOCATE (cutpotu(n))
111 rcut = maxval(upf_pot%r)
112 CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
114 ALLOCATE (cutpots(n))
115 CALL cutpot(cutpots, basis%grid%rad, rcut, 2.5_dp)
118 ALLOCATE (cutpots(n))
124 CALL ecp_sgp_constr(ecp_pot, sgp_pot, basis)
126 CALL upf_sgp_constr(upf_pot, sgp_pot, basis)
134 NULLIFY (score, shnl)
139 CALL ecpints(hnl%op, basis, ecp_pot)
141 CALL upfints(core%op, hnl%op, basis, upf_pot, cutpotu, sgp_pot%ac_local)
146 CALL sgpints(score%op, shnl%op, basis, sgp_pot, cutpots)
149 IF (sgp_pot%has_local)
THEN
150 n = min(3, ubound(core%op, 3))
151 error(1) = maxval(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
152 mloc = maxloc(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
153 error(4) = error(1)/max(abs(score%op(mloc(1), mloc(2), mloc(3))), 1.e-6_dp)
155 IF (sgp_pot%has_nonlocal)
THEN
156 n = min(3, ubound(hnl%op, 3))
157 error(2) = maxval(abs(hnl%op(:, :, 0:n) - shnl%op(:, :, 0:n)))
158 mloc = maxloc(abs(hnl%op(:, :, 0:n) - shnl%op(:, :, 0:n)))
159 error(5) = error(1)/max(abs(hnl%op(mloc(1), mloc(2), mloc(3))), 1.e-6_dp)
161 IF (sgp_pot%has_nlcc)
THEN
166 DO i = 1, sgp_pot%n_nlcc
167 cgauss(:) = cgauss(:) + sgp_pot%c_nlcc(i)*exp(-sgp_pot%a_nlcc(i)*upf_pot%r(:)**2)
169 errcc = sum((cgauss(:) - upf_pot%rho_nlcc(:))**2*upf_pot%r(:)**2*upf_pot%rab(:))
170 errcc = sqrt(errcc/real(n, kind=
dp))
202 TYPE(atom_p_type),
DIMENSION(:, :),
POINTER :: atom_info
203 TYPE(section_vals_type),
POINTER :: input_section
204 INTEGER,
INTENT(IN) :: iw
206 INTEGER :: i, n, ppot_type
207 LOGICAL :: do_transform, explicit, is_ecp, is_upf
208 REAL(kind=
dp) :: errcc, rcut
209 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cgauss, cutpots, cutpotu
210 TYPE(atom_ecppot_type),
POINTER :: ecp_pot
212 TYPE(atom_type),
POINTER :: atom_ref
214 TYPE(opmat_type),
POINTER :: core, hnl, score, shnl
217 IF (.NOT. explicit)
RETURN
219 IF (iw > 0)
WRITE (iw,
'(/," ",79("*"),/,T24,A,/," ",79("*"))')
"SEPARABLE GAUSSIAN PSEUDOPOTENTIAL"
221 atom_ref => atom_info(1, 1)%atom
223 ppot_type = atom_ref%potential%ppot_type
224 SELECT CASE (ppot_type)
226 IF (iw > 0)
WRITE (iw,
'(" GTH Pseudopotential is already in SGP form. ")')
227 do_transform = .false.
229 do_transform = .true.
232 ecp_pot => atom_ref%potential%ecp_pot
234 do_transform = .true.
237 upf_pot => atom_ref%potential%upf_pot
239 IF (iw > 0)
WRITE (iw,
'(" No Pseudopotential available for transformation. ")')
240 do_transform = .false.
246 IF (do_transform)
THEN
248 CALL ecp_sgp_constr(ecp_pot, sgp_pot, atom_ref%basis)
250 CALL upf_sgp_constr(upf_pot, sgp_pot, atom_ref%basis)
257 IF (do_transform)
THEN
261 NULLIFY (score, shnl)
268 ALLOCATE (cutpotu(n))
269 rcut = maxval(upf_pot%r)
270 CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
271 n = atom_ref%basis%grid%nr
272 ALLOCATE (cutpots(n))
273 CALL cutpot(cutpots, atom_ref%basis%grid%rad, rcut, 2.5_dp)
275 n = atom_ref%basis%grid%nr
276 ALLOCATE (cutpots(n))
281 CALL ecpints(hnl%op, atom_ref%basis, ecp_pot)
283 CALL upfints(core%op, hnl%op, atom_ref%basis, upf_pot, cutpotu, sgp_pot%ac_local)
288 CALL sgpints(score%op, shnl%op, atom_ref%basis, sgp_pot, cutpots)
290 IF (sgp_pot%has_local)
THEN
291 n = min(3, ubound(core%op, 3))
292 errcc = maxval(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
294 WRITE (iw,
'(" Local part of pseudopotential")')
295 WRITE (iw,
'(" Number of basis functions ",T77,i4)') sgp_pot%n_local
296 WRITE (iw,
'(" Max. abs. error of matrix elements ",T65,f16.8)') errcc
299 IF (sgp_pot%has_nonlocal)
THEN
301 WRITE (iw,
'(" Nonlocal part of pseudopotential")')
302 WRITE (iw,
'(" Max. l-quantum number",T77,i4)') sgp_pot%lmax
303 WRITE (iw,
'(" Number of projector basis functions ",T77,i4)') sgp_pot%n_nonlocal
304 DO i = 0, sgp_pot%lmax
305 errcc = maxval(abs(hnl%op(:, :, i) - shnl%op(:, :, i)))
306 WRITE (iw,
'(" Max. abs. error of matrix elements: l=",i2,T69,f12.8)') i, errcc
310 IF (sgp_pot%has_nlcc)
THEN
315 DO i = 1, sgp_pot%n_nlcc
316 cgauss(:) = cgauss(:) + sgp_pot%c_nlcc(i)*exp(-sgp_pot%a_nlcc(i)*upf_pot%r(:)**2)
318 errcc = sum((cgauss(:) - upf_pot%rho_nlcc(:))**2*upf_pot%r(:)**2*upf_pot%rab(:))
319 errcc = sqrt(errcc/real(n, kind=
dp))
325 WRITE (iw,
'(" Non-linear core correction: core density")')
326 WRITE (iw,
'(" Number of basis functions ",T77,i4)') sgp_pot%n_nlcc
327 WRITE (iw,
'(" RMS error of core density ",T69,f12.8)') errcc
346 IF (iw > 0)
WRITE (iw,
'(" ",79("*"))')
356 SUBROUTINE ecp_sgp_constr(ecp_pot, sgp_pot, basis)
358 TYPE(atom_ecppot_type) :: ecp_pot
360 TYPE(atom_basis_type) :: basis
362 INTEGER :: i, ia, ir, j, k, l, n, na, nl, nr
363 REAL(kind=
dp) :: alpha, amx, cmx
364 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: al, cl, cpot, pgauss
365 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: cmat, qmat, score, sinv, smat, tmat
366 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rad
368 sgp_pot%has_nlcc = .false.
370 sgp_pot%has_local = .false.
374 sgp_pot%has_nonlocal = .true.
379 sgp_pot%n_nonlocal = nl
380 sgp_pot%lmax = ecp_pot%lmax
381 ALLOCATE (sgp_pot%a_nonlocal(nl))
382 ALLOCATE (sgp_pot%h_nonlocal(nl, 0:ecp_pot%lmax))
383 ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:ecp_pot%lmax))
385 ALLOCATE (al(nl), cl(nl))
386 ALLOCATE (smat(nl, nl), sinv(nl, nl))
387 ALLOCATE (tmat(nl, nl), cmat(nl, nl))
389 amx = maxval(ecp_pot%bpot)
391 cmx = cmx**(1._dp/real(nl - 1,
dp))
394 al(ir) = amx*cmx**(ir - 1)
397 sgp_pot%a_nonlocal(1:nl) = al(1:nl)
400 rad => basis%grid%rad
401 ALLOCATE (cpot(nr), pgauss(nr))
402 DO l = 0, ecp_pot%lmax
404 ALLOCATE (score(na, na), qmat(na, nl))
406 DO k = 1, ecp_pot%npot(l)
407 n = ecp_pot%nrpot(k, l)
408 alpha = ecp_pot%bpot(k, l)
409 cpot(:) = cpot + ecp_pot%apot(k, l)*rad**(n - 2)*exp(-alpha*rad**2)
413 score(i, j) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
414 score(j, i) = score(i, j)
419 pgauss(:) = exp(-al(i)*rad(:)**2)*rad(:)**l
421 qmat(ia, i) = integrate_grid(basis%bf(:, ia, l), pgauss(:), basis%grid)
426 tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
427 smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
429 cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
430 cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
434 CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
436 sgp_pot%h_nonlocal(1:nl, l) = cl(1:nl)
437 sgp_pot%c_nonlocal(1:nl, 1:nl, l) = cmat(1:nl, 1:nl)
438 sgp_pot%is_nonlocal(l) = .true.
440 DEALLOCATE (score, qmat)
442 DEALLOCATE (cpot, pgauss)
443 DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
445 END SUBROUTINE ecp_sgp_constr
453 SUBROUTINE upf_sgp_constr(upf_pot, sgp_pot, basis)
457 TYPE(atom_basis_type) :: basis
459 CHARACTER(len=4) :: ptype
460 INTEGER :: ia, ib, ipa, ipb, ir, la, lb, na, nl, &
463 REAL(kind=
dp) :: cpa, cpb, eee, ei, errcc, errloc, rc, &
465 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: al, ccharge, cgauss, cl, pgauss, pproa, &
466 pprob, tv, vgauss, vloc, ww
467 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: cmat, qmat, score, sinv, smat, tmat
468 TYPE(atom_basis_type) :: gbasis
471 IF (upf_pot%is_ultrasoft .OR. upf_pot%is_paw .OR. upf_pot%is_coulomb)
THEN
472 sgp_pot%has_nonlocal = .false.
473 sgp_pot%n_nonlocal = 0
474 sgp_pot%has_local = .false.
476 sgp_pot%has_nlcc = .false.
485 ww(:) = upf_pot%r(:)**2*upf_pot%rab(:)
488 sgp_pot%has_local = .true.
490 ALLOCATE (vloc(nr), vgauss(nr))
494 CALL erffit(sgp_pot%ac_local, upf_pot%vlocal, upf_pot%r, zval)
497 IF (upf_pot%r(ir) < 1.e-12_dp)
THEN
498 vgauss(ir) = -sqrt(2.0_dp)*zval/
rootpi/sgp_pot%ac_local
500 rc = upf_pot%r(ir)/sgp_pot%ac_local/sqrt(2.0_dp)
501 vgauss(ir) = -zval/upf_pot%r(ir)*erf(rc)
504 vloc(:) = upf_pot%vlocal(:) - vgauss(:)
509 ALLOCATE (al(nl), cl(nl))
514 ostate%rhoend = 1.e-12_dp
515 ostate%rhobeg = 5.e-2_dp
521 IF (ostate%state == 2)
THEN
523 al(ir) = x(1)*x(2)**(ir - 1)
525 CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, ostate%f)
527 IF (ostate%state == -1)
EXIT
533 al(ir) = x(1)*x(2)**(ir - 1)
535 CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, errloc)
537 ALLOCATE (sgp_pot%a_local(nl), sgp_pot%c_local(nl))
539 sgp_pot%a_local(1:nl) = al(1:nl)
540 sgp_pot%c_local(1:nl) = cl(1:nl)
541 DEALLOCATE (vloc, vgauss)
545 ptype = adjustl(trim(upf_pot%pseudo_type))
546 IF (ptype(1:2) ==
"NC" .OR. ptype(1:2) ==
"US")
THEN
548 ELSE IF (ptype(1:2) ==
"SL")
THEN
551 cpabort(
"Pseudopotential type: ["//adjustl(trim(ptype))//
"] not known")
555 IF (upf_pot%l_max < 0)
THEN
556 sgp_pot%n_nonlocal = 0
558 sgp_pot%has_nonlocal = .false.
562 sgp_pot%has_nonlocal = .true.
568 sgp_pot%n_nonlocal = nl
569 sgp_pot%lmax = upf_pot%l_max
570 ALLOCATE (sgp_pot%a_nonlocal(nl))
571 ALLOCATE (sgp_pot%h_nonlocal(nl, 0:upf_pot%l_max))
572 ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:upf_pot%l_max))
574 ALLOCATE (al(nl), cl(nl))
575 ALLOCATE (smat(nl, nl), sinv(nl, nl))
576 ALLOCATE (tmat(nl, nl), cmat(nl, nl))
579 al(ir) = 10.0_dp*0.60_dp**(ir - 1)
582 sgp_pot%a_nonlocal(1:nl) = al(1:nl)
585 ALLOCATE (pgauss(nr), vloc(nr))
586 DO la = 0, upf_pot%l_max
587 IF (la == upf_pot%l_local) cycle
588 sgp_pot%is_nonlocal(la) = .true.
590 ALLOCATE (score(na, na), qmat(na, nl))
592 vloc(:) = upf_pot%vsemi(:, la + 1) - upf_pot%vlocal(:)
595 score(ia, ib) = sum(vloc(:)*gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*ww(:))
596 score(ib, ia) = score(ia, ib)
601 pgauss(:) = exp(-al(ir)*upf_pot%r(:)**2)*upf_pot%r(:)**la
602 eee =
rootpi/(2._dp**(la + 2)*
dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
603 pgauss(:) = pgauss(:)/sqrt(eee)
605 qmat(ia, ir) = sum(gbasis%bf(:, ia, la)*pgauss(:)*ww)
609 tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
610 smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
612 cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
613 cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
614 CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
618 ei =
rootpi/(2._dp**(la + 2)*
dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
619 cmat(ir, 1:nl) = cmat(ir, 1:nl)/sqrt(ei)
621 sgp_pot%h_nonlocal(1:nl, la) = cl(1:nl)
622 sgp_pot%c_nonlocal(1:nl, 1:nl, la) = cmat(1:nl, 1:nl)
623 sgp_pot%is_nonlocal(la) = .true.
624 DEALLOCATE (score, qmat)
627 sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/sqrt(
fourpi)
629 DEALLOCATE (pgauss, vloc)
630 DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
632 sgp_pot%has_nonlocal = .true.
634 ALLOCATE (pproa(nr), pprob(nr), pgauss(nr))
635 np = upf_pot%number_of_proj
637 ALLOCATE (al(nl), cl(nl))
638 ALLOCATE (smat(nl, nl), sinv(nl, nl))
639 ALLOCATE (tmat(nl, nl), cmat(nl, nl))
643 al(ir) = 10.0_dp*0.60_dp**(ir - 1)
646 sgp_pot%lmax = maxval(upf_pot%lbeta(:))
647 sgp_pot%n_nonlocal = nl
648 ALLOCATE (sgp_pot%a_nonlocal(nl))
649 ALLOCATE (sgp_pot%h_nonlocal(nl, 0:sgp_pot%lmax))
650 ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:sgp_pot%lmax))
652 sgp_pot%a_nonlocal(1:nl) = al(1:nl)
655 DO la = 0, sgp_pot%lmax
656 sgp_pot%is_nonlocal(la) = .true.
658 ALLOCATE (score(na, na), qmat(na, nl))
662 lb = upf_pot%lbeta(ipa)
664 pproa(:) = upf_pot%beta(:, ipa)
666 lb = upf_pot%lbeta(ipb)
668 pprob(:) = upf_pot%beta(:, ipb)
669 eee = upf_pot%dion(ipa, ipb)
671 cpa = sum(pproa(:)*gbasis%bf(:, ia, la)*ww(:))
673 cpb = sum(pprob(:)*gbasis%bf(:, ib, la)*ww(:))
674 score(ia, ib) = score(ia, ib) + cpa*eee*cpb
675 score(ib, ia) = score(ia, ib)
682 pgauss(:) = exp(-al(ir)*upf_pot%r(:)**2)*upf_pot%r(:)**la
683 eee =
rootpi/(2._dp**(la + 2)*
dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
684 pgauss(:) = pgauss(:)/sqrt(eee)
686 qmat(ia, ir) = sum(gbasis%bf(:, ia, la)*pgauss(:)*ww)
690 tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
691 smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
693 cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
694 cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
695 CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
699 ei =
rootpi/(2._dp**(la + 2)*
dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
700 cmat(ir, 1:nl) = cmat(ir, 1:nl)/sqrt(ei)
702 sgp_pot%h_nonlocal(1:nl, la) = cl(1:nl)
703 sgp_pot%c_nonlocal(1:nl, 1:nl, la) = cmat(1:nl, 1:nl)
704 sgp_pot%is_nonlocal(la) = .true.
705 DEALLOCATE (score, qmat)
708 sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/sqrt(
fourpi)
710 DEALLOCATE (pgauss, pproa, pprob)
711 DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
715 IF (upf_pot%core_correction)
THEN
716 sgp_pot%has_nlcc = .true.
718 sgp_pot%has_nlcc = .false.
723 IF (sgp_pot%has_nlcc)
THEN
724 ALLOCATE (ccharge(nr), cgauss(nr))
725 ccharge(:) = upf_pot%rho_nlcc(:)
727 ALLOCATE (al(nl), cl(nl), tv(nl))
728 ALLOCATE (smat(nl, nl), sinv(nl, nl))
732 al(ir) = 10.0_dp*0.6_dp**(ir - 1)
738 CALL sg_overlap(smat(1:nl, 1:nl), 0, al(1:nl), al(1:nl))
740 cgauss(:) = exp(-al(ir)*upf_pot%r(:)**2)
741 tv(ir) = sum(cgauss*ccharge*ww)
744 cl(1:nl) = matmul(sinv(1:nl, 1:nl), tv(1:nl))
747 cgauss(:) = cgauss(:) + cl(ir)*exp(-al(ir)*upf_pot%r(:)**2)
749 errcc = sum((cgauss - ccharge)**2*ww)
750 ALLOCATE (sgp_pot%a_local(nl), sgp_pot%c_local(nl))
752 sgp_pot%a_nlcc(1:nl) = al(1:nl)
753 sgp_pot%c_nlcc(1:nl) = cl(1:nl)
754 DEALLOCATE (ccharge, cgauss)
755 DEALLOCATE (al, cl, tv, smat, sinv)
760 END SUBROUTINE upf_sgp_constr
770 IF (
ASSOCIATED(sgp_pot%a_nonlocal))
DEALLOCATE (sgp_pot%a_nonlocal)
771 IF (
ASSOCIATED(sgp_pot%h_nonlocal))
DEALLOCATE (sgp_pot%h_nonlocal)
772 IF (
ASSOCIATED(sgp_pot%c_nonlocal))
DEALLOCATE (sgp_pot%c_nonlocal)
774 IF (
ASSOCIATED(sgp_pot%a_local))
DEALLOCATE (sgp_pot%a_local)
775 IF (
ASSOCIATED(sgp_pot%c_local))
DEALLOCATE (sgp_pot%c_local)
777 IF (
ASSOCIATED(sgp_pot%a_nlcc))
DEALLOCATE (sgp_pot%a_nlcc)
778 IF (
ASSOCIATED(sgp_pot%c_nlcc))
DEALLOCATE (sgp_pot%c_nlcc)
791 SUBROUTINE upfints(core, hnl, basis, upf_pot, cutpotu, ac_local)
792 REAL(kind=
dp),
DIMENSION(:, :, 0:) :: core, hnl
793 TYPE(atom_basis_type),
INTENT(INOUT) :: basis
795 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: cutpotu
796 REAL(kind=
dp),
INTENT(IN) :: ac_local
798 CHARACTER(len=4) :: ptype
799 INTEGER :: i, j, k1, k2, la, lb, m, n
800 REAL(kind=
dp) :: rc, zval
801 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: spot
802 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: spmat
803 TYPE(atom_basis_type) :: gbasis
812 spot(:) = upf_pot%vlocal(:)
815 IF (upf_pot%r(i) < 1.e-12_dp)
THEN
818 rc = upf_pot%r(i)/ac_local/
sqrt2
819 spot(i) = spot(i) + zval/upf_pot%r(i)*erf(rc)
822 spot(:) = spot(:)*cutpotu(:)
828 ptype = adjustl(trim(upf_pot%pseudo_type))
829 IF (ptype(1:2) ==
"NC" .OR. ptype(1:2) ==
"US")
THEN
831 n = maxval(gbasis%nbas(:))
832 m = upf_pot%number_of_proj
833 ALLOCATE (spmat(n, m))
836 la = upf_pot%lbeta(i)
837 DO j = 1, gbasis%nbas(la)
838 spmat(j, i) = integrate_grid(upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
842 la = upf_pot%lbeta(i)
844 lb = upf_pot%lbeta(j)
846 DO k1 = 1, gbasis%nbas(la)
847 DO k2 = 1, gbasis%nbas(la)
848 hnl(k1, k2, la) = hnl(k1, k2, la) + spmat(k1, i)*upf_pot%dion(i, j)*spmat(k2, j)
855 ELSE IF (ptype(1:2) ==
"SL")
THEN
857 DO la = 0, upf_pot%l_max
858 IF (la == upf_pot%l_local) cycle
859 m =
SIZE(upf_pot%vsemi(:, la + 1))
861 spot(:) = upf_pot%vsemi(:, la + 1) - upf_pot%vlocal(:)
862 spot(:) = spot(:)*cutpotu(:)
866 hnl(i, j, la) = hnl(i, j, la) + &
867 integrate_grid(spot(:), &
868 gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
869 hnl(j, i, la) = hnl(i, j, la)
875 cpabort(
"Pseudopotential type: ["//adjustl(trim(ptype))//
"] not known")
881 END SUBROUTINE upfints
889 SUBROUTINE ecpints(hnl, basis, ecp_pot)
890 REAL(kind=
dp),
DIMENSION(:, :, 0:) :: hnl
891 TYPE(atom_basis_type),
INTENT(INOUT) :: basis
892 TYPE(atom_ecppot_type) :: ecp_pot
894 INTEGER :: i, j, k, l, m, n
895 REAL(kind=
dp) :: alpha
896 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cpot
897 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rad
899 rad => basis%grid%rad
905 DO l = 0, ecp_pot%lmax
907 DO k = 1, ecp_pot%npot(l)
908 n = ecp_pot%nrpot(k, l)
909 alpha = ecp_pot%bpot(k, l)
910 cpot(:) = cpot(:) + ecp_pot%apot(k, l)*rad(:)**(n - 2)*exp(-alpha*rad(:)**2)
912 DO i = 1, basis%nbas(l)
913 DO j = i, basis%nbas(l)
914 hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
915 hnl(j, i, l) = hnl(i, j, l)
921 END SUBROUTINE ecpints
931 SUBROUTINE sgpints(core, hnl, basis, sgp_pot, cutpots)
932 REAL(kind=
dp),
DIMENSION(:, :, 0:) :: core, hnl
933 TYPE(atom_basis_type),
INTENT(INOUT) :: basis
935 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: cutpots
937 INTEGER :: i, ia, j, l, m, n, na
938 REAL(kind=
dp) :: a, zval
939 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cpot, pgauss
940 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: qmat
941 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rad
943 rad => basis%grid%rad
948 IF (sgp_pot%has_local)
THEN
952 DO i = 1, sgp_pot%n_local
953 cpot(:) = cpot(:) + sgp_pot%c_local(i)*exp(-sgp_pot%a_local(i)*rad(:)**2)
955 cpot(:) = cpot(:)*cutpots(:)
961 IF (sgp_pot%has_nonlocal)
THEN
963 ALLOCATE (pgauss(1:m))
964 n = sgp_pot%n_nonlocal
966 DO l = 0, sgp_pot%lmax
967 cpassert(l <= ubound(basis%nbas, 1))
968 IF (.NOT. sgp_pot%is_nonlocal(l)) cycle
971 ALLOCATE (qmat(na, n))
973 a = sgp_pot%a_nonlocal(i)
974 pgauss(:) = cutpots(:)*exp(-a*rad(:)**2)*rad(:)**l
976 qmat(ia, i) = integrate_grid(basis%bf(:, ia, l), pgauss(:), basis%grid)
979 qmat(1:na, 1:n) = sqrt(
fourpi)*matmul(qmat(1:na, 1:n), sgp_pot%c_nonlocal(1:n, 1:n, l))
983 hnl(i, j, l) = hnl(i, j, l) + qmat(i, ia)*qmat(j, ia)*sgp_pot%h_nonlocal(ia, l)
985 hnl(j, i, l) = hnl(i, j, l)
993 END SUBROUTINE sgpints
1002 SUBROUTINE erffit(ac, vlocal, r, z)
1003 REAL(kind=
dp),
INTENT(OUT) :: ac
1004 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: vlocal, r
1005 REAL(kind=
dp),
INTENT(IN) :: z
1007 REAL(kind=
dp),
PARAMETER :: rcut = 1.4_dp
1009 INTEGER :: i, j, m, m1
1010 REAL(kind=
dp) :: a1, a2, an, e2, en, rc
1011 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: epot, rval, vpot
1014 ALLOCATE (epot(m), vpot(m), rval(m))
1015 cpassert(
SIZE(vlocal) == m)
1016 IF (r(1) > r(m))
THEN
1018 vpot(m - i + 1) = vlocal(i)
1019 rval(m - i + 1) = r(i)
1022 vpot(1:m) = vlocal(1:m)
1027 IF (rval(i) > rcut)
THEN
1038 an = a1 + i*0.025_dp
1039 rc = 1._dp/(an*sqrt(2.0_dp))
1041 epot(j) = vpot(j) + z/rval(j)*erf(rval(j)*rc)
1043 en = sum(abs(epot(m1:m)*rval(m1:m)**2))
1051 DEALLOCATE (epot, vpot, rval)
1053 END SUBROUTINE erffit
1068 SUBROUTINE pplocal_error(nl, al, cl, vloc, vgauss, gbasis, rad, ww, method, errloc)
1070 REAL(kind=
dp),
DIMENSION(:) :: al, cl, vloc, vgauss
1071 TYPE(atom_basis_type) :: gbasis
1072 REAL(kind=
dp),
DIMENSION(:) :: rad, ww
1073 INTEGER,
INTENT(IN) :: method
1074 REAL(kind=
dp) :: errloc
1076 INTEGER :: ia, ib, ir, ix, la, na
1077 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tv
1078 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: rmat, sinv, smat
1079 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: gmat
1082 IF (method == 1)
THEN
1083 ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
1085 vgauss(:) = exp(-al(ir)*rad(:)**2)
1087 smat(ir, ix) = sum(vgauss(:)*exp(-al(ix)*rad(:)**2)*ww(:))
1089 tv(ir) = sum(vloc(:)*vgauss(:)*ww(:))
1092 cl(1:nl) = matmul(sinv(1:nl, 1:nl), tv(1:nl))
1095 ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
1099 DO la = 0, min(ubound(gbasis%nbas, 1), 3)
1100 na = gbasis%nbas(la)
1101 ALLOCATE (rmat(na, na), gmat(na, na, nl))
1106 rmat(ia, ib) = sum(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vloc(:)*ww(:))
1107 rmat(ib, ia) = rmat(ia, ib)
1111 vgauss(:) = exp(-al(ir)*rad(:)**2)
1114 gmat(ia, ib, ir) = sum(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vgauss(:)*ww(:))
1115 gmat(ib, ia, ir) = gmat(ia, ib, ir)
1120 tv(ir) = tv(ir) + accurate_dot_product(rmat, gmat(:, :, ir))
1122 smat(ir, ix) = smat(ir, ix) + accurate_dot_product(gmat(:, :, ix), gmat(:, :, ir))
1123 smat(ix, ir) = smat(ir, ix)
1126 DEALLOCATE (rmat, gmat)
1130 cl(1:nl) = matmul(sinv(1:nl, 1:nl), tv(1:nl))
1135 vgauss(:) = vgauss(:) + cl(ir)*exp(-al(ir)*rad(:)**2)
1137 errloc = sum((vgauss - vloc)**2*ww)
1139 DEALLOCATE (tv, smat, sinv)
1141 END SUBROUTINE pplocal_error
1150 SUBROUTINE cutpot(pot, r, rcut, rsmooth)
1151 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: pot
1152 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: r
1153 REAL(kind=
dp),
INTENT(IN) :: rcut, rsmooth
1156 REAL(kind=
dp) :: rab, rx, x
1159 cpassert(n <=
SIZE(r))
1164 IF (rab > rcut)
THEN
1166 ELSE IF (rab > rcut - rsmooth)
THEN
1167 rx = rab - (rcut - rsmooth)
1169 pot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
1173 END SUBROUTINE cutpot
subroutine, public sg_overlap(smat, l, pa, pb)
...
subroutine, public set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid, cp2k_norm)
...
subroutine, public atom_sgp_release(sgp_pot)
...
subroutine, public sgp_construction(sgp_pot, ecp_pot, upf_pot, orb_basis, error)
...
subroutine, public atom_sgp_construction(atom_info, input_section, iw)
...
Define the atom type and its sub types.
subroutine, public release_opmat(opmat)
...
subroutine, public release_atom_basis(basis)
...
subroutine, public init_atom_basis_default_pp(basis)
...
subroutine, public create_opmat(opmat, n, lmax)
...
subroutine, public atom_basis_gridrep(basis, gbasis, r, rab)
...
Routines that process Quantum Espresso UPF files.
Some basic routines for atomic calculations.
subroutine, public numpot_matrix(imat, cpot, basis, derivatives)
Calculate a potential matrix by integrating the potential on an atomic radial grid.
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
real(kind=dp), parameter, public rootpi
real(kind=dp), parameter, public fourpi
real(kind=dp), parameter, public sqrt2
Collection of simple mathematical functions and subroutines.
subroutine, public get_pseudo_inverse_diag(a, a_pinverse, rskip)
returns the pseudoinverse of a real, symmetric and positive definite matrix using diagonalization.
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
subroutine, public powell_optimize(n, x, optstate)
...