Geometrical Response Basis Sets generated using CP2K Atom Code

svn-origin-rev: 17568
This commit is contained in:
Jürg Hutter 2016-11-25 16:05:28 +00:00
parent 10bff1e61c
commit 368a85fda5
11 changed files with 1549 additions and 55 deletions

View file

@ -42,7 +42,7 @@ MODULE ai_onecenter
PUBLIC :: sg_overlap, sg_kinetic, sg_nuclear, sg_erf, sg_gpot, &
sg_proj_ol, sg_coulomb, sg_exchange, sg_kinnuc, &
sg_erfc
sg_erfc, sg_conf
PUBLIC :: sto_overlap, sto_kinetic, sto_nuclear
CONTAINS

View file

@ -242,7 +242,6 @@ CONTAINS
! wavefunction history
NULLIFY (history%dmat, history%hmat)
CALL atom_history_init(history, atom%optimization, fmat%op)
is_converged = .FALSE.
iter = 0
DO !SCF Loop

View file

@ -9,6 +9,7 @@ MODULE atom_energy
USE atom_electronic_structure, ONLY: calculate_atom
USE atom_fit, ONLY: atom_fit_density,&
atom_fit_kgpot
USE atom_grb, ONLY: atom_grb_construction
USE atom_operators, ONLY: atom_int_release,&
atom_int_setup,&
atom_ppint_release,&
@ -433,6 +434,13 @@ CONTAINS
END DO
END DO
! generate a geometrical response basis (GRB)
iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS", extension=".log")
IF (iw > 0) THEN
CALL atom_grb_construction(atom_info, atom_section, iw)
END IF
CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
! Analyze basis sets
iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ANALYZE_BASIS", extension=".log")
IF (iw > 0) THEN
@ -590,7 +598,7 @@ CONTAINS
amat(i, i) = amat(i, i)-1._dp
END DO
IF (MAXVAL(ABS(amat)) > 1.e-12) THEN
WRITE (iw, '(A,G20.10)') " Orthogonality error ", MAXVAL(ABS(amat))
WRITE (iw, '(/,A,G20.10)') " Orthogonality error ", MAXVAL(ABS(amat))
END IF
DEALLOCATE (amat)

View file

@ -313,6 +313,7 @@ CONTAINS
END IF
basis%bf = 0._dp
basis%dbf = 0._dp
basis%ddbf = 0._dp
nr = basis%grid%nr
DO l = 0, 3
DO i = 1, basis%nbas(l)
@ -322,6 +323,8 @@ CONTAINS
ear = EXP(-al*basis%grid%rad(k)**2)
basis%bf(k, i, l) = rk**l*ear
basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l-1)-2._dp*al*rk**(l+1))*ear
basis%ddbf(k, i, l) = (REAL(l*(l-1), dp)*rk**(l-2)- &
2._dp*al*REAL(2*l+1, dp)*rk**(l)+4._dp*al*rk**(l+2))*ear
END DO
END DO
END DO
@ -334,6 +337,7 @@ CONTAINS
END DO
basis%bf = 0._dp
basis%dbf = 0._dp
basis%ddbf = 0._dp
nr = basis%grid%nr
DO l = 0, 3
DO i = 1, basis%nbas(l)
@ -342,9 +346,11 @@ CONTAINS
pf = (2._dp*al)**nl*SQRT(2._dp*al/fac(2*nl))
DO k = 1, nr
rk = basis%grid%rad(k)
ear = EXP(-al*basis%grid%rad(k))
basis%bf(k, i, l) = pf*rk**(nl-1)*ear
basis%dbf(k, i, l) = pf*(REAL(nl-1, dp)/rk-al)*rk**(nl-1)*ear
ear = rk**(nl-1)*EXP(-al*rk)
basis%bf(k, i, l) = pf*ear
basis%dbf(k, i, l) = pf*(REAL(nl-1, dp)/rk-al)*ear
basis%ddbf(k, i, l) = pf*(REAL((nl-2)*(nl-1), dp)/rk/rk &
-al*REAL(2*(nl-1), dp)/rk+al*al)*ear
END DO
END DO
END DO
@ -398,6 +404,7 @@ CONTAINS
END IF
basis%bf = 0._dp
basis%dbf = 0._dp
basis%ddbf = 0._dp
nr = basis%grid%nr
DO l = 0, 3
DO i = 1, basis%nbas(l)
@ -407,6 +414,8 @@ CONTAINS
ear = EXP(-al*basis%grid%rad(k)**2)
basis%bf(k, i, l) = rk**l*ear
basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l-1)-2._dp*al*rk**(l+1))*ear
basis%ddbf(k, i, l) = (REAL(l*(l-1), dp)*rk**(l-2)- &
2._dp*al*REAL(2*l+1, dp)*rk**(l)+4._dp*al*rk**(l+2))*ear
END DO
END DO
END DO
@ -419,6 +428,7 @@ CONTAINS
END DO
basis%bf = 0._dp
basis%dbf = 0._dp
basis%ddbf = 0._dp
nr = basis%grid%nr
DO l = 0, 3
DO i = 1, basis%nbas(l)
@ -427,9 +437,11 @@ CONTAINS
pf = (2._dp*al)**nl*SQRT(2._dp*al/fac(2*nl))
DO k = 1, nr
rk = basis%grid%rad(k)
ear = EXP(-al*basis%grid%rad(k))
basis%bf(k, i, l) = pf*rk**(nl-1)*ear
basis%dbf(k, i, l) = pf*(REAL(nl-1, dp)/rk-al)*rk**(nl-1)*ear
ear = rk**(nl-1)*EXP(-al*rk)
basis%bf(k, i, l) = pf*ear
basis%dbf(k, i, l) = pf*(REAL(nl-1, dp)/rk-al)*ear
basis%ddbf(k, i, l) = pf*(REAL((nl-2)*(nl-1), dp)/rk/rk &
-al*REAL(2*(nl-1), dp)/rk+al*al)*ear
END DO
END DO
END DO
@ -476,15 +488,15 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'atom_fit_pseudo', &
routineP = moduleN//':'//routineN
CHARACTER(len=2) :: pc1, pc2
CHARACTER(len=2) :: pc1, pc2, pct
INTEGER :: i, i1, i2, iw, j, j1, j2, k, l, m, n, &
n10, np, nre, nreinit, ntarget
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: xtob
INTEGER, DIMENSION(0:3) :: ncore
LOGICAL :: explicit
REAL(KIND=dp) :: afun, charge, de, deig, drho, dx, eig, fopt, oc, pchg, peig, pv, rcm, rcov, &
rmax, t_ener, t_semi, t_valence, t_virt, w_ener, w_node, w_psir0, w_semi, w_valence, &
w_virt, wtot
rmax, semicore_level, t_ener, t_psir0, t_semi, t_valence, t_virt, w_ener, w_node, &
w_psir0, w_semi, w_valence, w_virt, wtot
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x, xi
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: wem
REAL(KIND=dp), ALLOCATABLE, &
@ -506,6 +518,7 @@ CONTAINS
CALL section_vals_val_get(powell_section, "WEIGHT_POT_NODE", r_val=w_node)
CALL section_vals_val_get(powell_section, "WEIGHT_DELTA_ENERGY", r_val=w_ener)
CALL section_vals_val_get(powell_section, "TARGET_PSIR0", r_val=t_psir0)
CALL section_vals_val_get(powell_section, "WEIGHT_PSIR0", r_val=w_psir0)
CALL section_vals_val_get(powell_section, "RCOV_MULTIPLICATION", r_val=rcm)
@ -514,6 +527,8 @@ CONTAINS
CALL section_vals_val_get(powell_section, "TARGET_POT_SEMICORE", r_val=t_semi)
CALL section_vals_val_get(powell_section, "TARGET_DELTA_ENERGY", r_val=t_ener)
CALL section_vals_val_get(powell_section, "SEMICORE_LEVEL", r_val=semicore_level)
n = SIZE(atom_info, 1)
m = SIZE(atom_info, 2)
ALLOCATE (wem(n, m))
@ -598,12 +613,14 @@ CONTAINS
atom%orbitals%wrefene(k, l, 1) = w_virt
atom%orbitals%wrefchg(k, l, 1) = w_virt/100._dp
atom%orbitals%crefene(k, l, 1) = t_virt
atom%orbitals%reftype(k, l, 1) = "U1"
ntarget = ntarget+2
wtot = wtot+atom%weight*(w_virt+w_virt/100._dp)
ELSE
atom%orbitals%wrefene(k, l, 1) = w_virt/100._dp
atom%orbitals%wrefchg(k, l, 1) = 0._dp
atom%orbitals%crefene(k, l, 1) = t_virt*10._dp
atom%orbitals%reftype(k, l, 1) = "U2"
ntarget = ntarget+1
wtot = wtot+atom%weight*w_virt/100._dp
END IF
@ -611,23 +628,27 @@ CONTAINS
atom%orbitals%wrefene(k, l, 1) = w_semi
atom%orbitals%wrefchg(k, l, 1) = w_semi/100._dp
atom%orbitals%crefene(k, l, 1) = t_semi
atom%orbitals%reftype(k, l, 1) = "SC"
ntarget = ntarget+2
wtot = wtot+atom%weight*(w_semi+w_semi/100._dp)
ELSE
IF (ABS(atom%state%occupation(l, k)-REAL(4*l+2, KIND=dp)) < 0.01_dp .AND. &
ABS(atom%orbitals%refene(k, l, 1)*evolt) > 30._dp) THEN
ABS(atom%orbitals%refene(k, l, 1)) > semicore_level) THEN
! full shell semicore
atom%orbitals%wrefene(k, l, 1) = w_semi
atom%orbitals%wrefchg(k, l, 1) = w_semi/100._dp
atom%orbitals%crefene(k, l, 1) = t_semi
atom%orbitals%reftype(k, l, 1) = "SC"
wtot = wtot+atom%weight*(w_semi+w_semi/100._dp)
ELSE
atom%orbitals%wrefene(k, l, 1) = w_valence
atom%orbitals%wrefchg(k, l, 1) = w_valence/100._dp
atom%orbitals%crefene(k, l, 1) = t_valence
atom%orbitals%reftype(k, l, 1) = "VA"
wtot = wtot+atom%weight*(w_valence+w_valence/100._dp)
END IF
IF (l == 0) THEN
atom%orbitals%tpsir0(k, 1) = t_psir0
atom%orbitals%wpsir0(k, 1) = w_psir0
wtot = wtot+atom%weight*w_psir0
END IF
@ -670,12 +691,14 @@ CONTAINS
atom%orbitals%wrefene(k, l, 1:2) = w_virt
atom%orbitals%wrefchg(k, l, 1:2) = w_virt/100._dp
atom%orbitals%crefene(k, l, 1:2) = t_virt
atom%orbitals%reftype(k, l, 1:2) = "U1"
ntarget = ntarget+4
wtot = wtot+atom%weight*2._dp*(w_virt+w_virt/100._dp)
ELSE
atom%orbitals%wrefene(k, l, 1:2) = w_virt/100._dp
atom%orbitals%wrefchg(k, l, 1:2) = 0._dp
atom%orbitals%crefene(k, l, 1:2) = t_virt*10.0_dp
atom%orbitals%reftype(k, l, 1:2) = "U2"
wtot = wtot+atom%weight*2._dp*w_virt/100._dp
ntarget = ntarget+2
END IF
@ -683,23 +706,27 @@ CONTAINS
atom%orbitals%wrefene(k, l, 1:2) = w_semi
atom%orbitals%wrefchg(k, l, 1:2) = w_semi/100._dp
atom%orbitals%crefene(k, l, 1:2) = t_semi
atom%orbitals%reftype(k, l, 1:2) = "SC"
ntarget = ntarget+4
wtot = wtot+atom%weight*2._dp*(w_semi+w_semi/100._dp)
ELSE
IF (ABS(atom%state%occupation(l, k)-REAL(2*l+1, KIND=dp)) < 0.01_dp .AND. &
ABS(atom%orbitals%refene(k, l, 1)*evolt) > 30._dp) THEN
ABS(atom%orbitals%refene(k, l, 1)) > semicore_level) THEN
atom%orbitals%wrefene(k, l, 1:2) = w_semi
atom%orbitals%wrefchg(k, l, 1:2) = w_semi/100._dp
atom%orbitals%crefene(k, l, 1:2) = t_semi
atom%orbitals%reftype(k, l, 1:2) = "SC"
wtot = wtot+atom%weight*2._dp*(w_semi+w_semi/100._dp)
ELSE
atom%orbitals%wrefene(k, l, 1:2) = w_valence
atom%orbitals%wrefchg(k, l, 1:2) = w_valence/100._dp
atom%orbitals%crefene(k, l, 1:2) = t_valence
atom%orbitals%reftype(k, l, 1:2) = "VA"
wtot = wtot+atom%weight*2._dp*(w_valence+w_valence/100._dp)
END IF
ntarget = ntarget+4
IF (l == 0) THEN
atom%orbitals%tpsir0(k, 1:2) = t_psir0
atom%orbitals%wpsir0(k, 1:2) = w_psir0
wtot = wtot+atom%weight*2._dp*w_psir0
END IF
@ -779,8 +806,9 @@ CONTAINS
ELSE
WRITE (pc2, "(I2)") MIN(NINT(pchg), 99)
END IF
WRITE (iunit, '(I5,I5,F14.2,F21.10,F14.6,"[",A2,"]",F13.6,"[",A2,"]")') &
l, k, oc, eig*evolt, deig*evolt, pc1, drho, pc2
pct = atom%orbitals%reftype(k, l, 1)
WRITE (iunit, '(I5,I5,F14.2,F21.10,A3,F11.6,"[",A2,"]",F13.6,"[",A2,"]")') &
l, k, oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2
END DO
END IF
END DO
@ -809,8 +837,9 @@ CONTAINS
ELSE
WRITE (pc2, "(I2)") MIN(NINT(pchg), 99)
END IF
WRITE (iunit, '(I5,I5,2X,A5,F11.2,F19.10,F13.6,"[",A2,"]",F12.6,"[",A2,"]")') &
l, k, "alpha", oc, eig*evolt, deig*evolt, pc1, drho, pc2
pct = atom%orbitals%reftype(k, l, 1)
WRITE (iunit, '(I5,I5,2X,A5,F11.2,F19.10,A3,F10.6,"[",A2,"]",F12.6,"[",A2,"]")') &
l, k, "alpha", oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2
oc = atom%state%occb(l, k)
eig = atom%orbitals%enerb(k, l)
deig = eig-atom%orbitals%refene(k, l, 2)
@ -829,8 +858,9 @@ CONTAINS
ELSE
WRITE (pc2, "(I2)") MIN(NINT(pchg), 99)
END IF
WRITE (iunit, '(I5,I5,2X,A5,F11.2,F19.10,F13.6,"[",A2,"]",F12.6,"[",A2,"]")') &
l, k, " beta", oc, eig*evolt, deig*evolt, pc1, drho, pc2
pct = atom%orbitals%reftype(k, l, 2)
WRITE (iunit, '(I5,I5,2X,A5,F11.2,F19.10,A3,F10.6,"[",A2,"]",F12.6,"[",A2,"]")') &
l, k, " beta", oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2
END DO
END IF
END DO
@ -950,16 +980,27 @@ CONTAINS
ELSE
WRITE (pc2, "(I2)") MIN(NINT(pchg), 99)
END IF
WRITE (iunit, '(I5,I5,F14.2,F21.10,F14.6,"[",A2,"]",F13.6,"[",A2,"]")') &
l, k, oc, eig*evolt, deig*evolt, pc1, drho, pc2
pct = atom%orbitals%reftype(k, l, 1)
WRITE (iunit, '(I5,I5,F14.2,F21.10,A3,F11.6,"[",A2,"]",F13.6,"[",A2,"]")') &
l, k, oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2
END DO
END IF
END DO
np = atom%state%maxn_calc(0)
DO k = 1, np
CALL atom_wfnr0(pv, atom%orbitals%wfn(:, k, 0), atom%basis)
pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun*100._dp
WRITE (iunit, '(" s-states"," N=",I5,T40,"Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') k, pv, NINT(pchg)
IF (ABS(atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp) THEN
IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 1))) THEN
pv = 0.0_dp
ELSE
pv = 10._dp*(ABS(pv)-ABS(atom%orbitals%tpsir0(k, 1)))
END IF
pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun
ELSE
pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun*100._dp
END IF
WRITE (iunit, '(" s-states"," N=",I5,T40,"Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') &
k, pv, NINT(pchg)
END DO
ELSE
!spin polarization
@ -986,8 +1027,9 @@ CONTAINS
ELSE
WRITE (pc2, "(I2)") MIN(NINT(pchg), 99)
END IF
WRITE (iunit, '(I5,I5,A,F11.2,F20.10,F13.6,"[",A2,"]",F11.6,"[",A2,"]")') &
l, k, " alpha", oc, eig*evolt, deig*evolt, pc1, drho, pc2
pct = atom%orbitals%reftype(k, l, 1)
WRITE (iunit, '(I5,I5,A,F11.2,F20.10,A3,F10.6,"[",A2,"]",F11.6,"[",A2,"]")') &
l, k, " alpha", oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2
oc = atom%state%occb(l, k)
eig = atom%orbitals%enerb(k, l)
deig = eig-atom%orbitals%refene(k, l, 2)
@ -1006,19 +1048,39 @@ CONTAINS
ELSE
WRITE (pc2, "(I2)") MIN(NINT(pchg), 99)
END IF
WRITE (iunit, '(I5,I5,A,F11.2,F20.10,F13.6,"[",A2,"]",F11.6,"[",A2,"]")') &
l, k, " beta", oc, eig*evolt, deig*evolt, pc1, drho, pc2
pct = atom%orbitals%reftype(k, l, 2)
WRITE (iunit, '(I5,I5,A,F11.2,F20.10,A3,F10.6,"[",A2,"]",F11.6,"[",A2,"]")') &
l, k, " beta", oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2
END DO
END IF
END DO
np = atom%state%maxn_calc(0)
DO k = 1, np
CALL atom_wfnr0(pv, atom%orbitals%wfna(:, k, 0), atom%basis)
pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun*100._dp
IF (ABS(atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp) THEN
IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 1))) THEN
pv = 0.0_dp
ELSE
pv = 10._dp*(ABS(pv)-ABS(atom%orbitals%tpsir0(k, 1)))
END IF
pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun
ELSE
pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun*100._dp
END IF
WRITE (iunit, '(" s-states"," N=",I5,T35,"Alpha Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') &
k, pv, NINT(pchg)
!
CALL atom_wfnr0(pv, atom%orbitals%wfnb(:, k, 0), atom%basis)
pchg = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv/afun*100._dp
IF (ABS(atom%orbitals%tpsir0(k, 2)) > 1.e-14_dp) THEN
IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 2))) THEN
pv = 0.0_dp
ELSE
pv = 10._dp*(ABS(pv)-ABS(atom%orbitals%tpsir0(k, 2)))
END IF
pchg = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv/afun
ELSE
pchg = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv/afun*100._dp
END IF
WRITE (iunit, '(" s-states"," N=",I5,T36,"Beta Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') &
k, pv, NINT(pchg)
END DO
@ -1320,7 +1382,16 @@ CONTAINS
IF (l == 0) THEN
IF (atom%orbitals%wpsir0(k, 1) > 0._dp) THEN
CALL atom_wfnr0(pv, atom%orbitals%wfn(:, k, l), atom%basis)
pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv
IF (ABS(atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp) THEN
IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 1))) THEN
pv = 0.0_dp
ELSE
pv = 10._dp*(ABS(pv)-ABS(atom%orbitals%tpsir0(k, 1)))
END IF
pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv
ELSE
pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv*100._dp
END IF
afun = afun+pv
END IF
END IF
@ -1373,12 +1444,30 @@ CONTAINS
IF (l == 0) THEN
IF (atom%orbitals%wpsir0(k, 1) > 0._dp) THEN
CALL atom_wfnr0(pv, atom%orbitals%wfna(:, k, l), atom%basis)
pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv
IF (ABS(atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp) THEN
IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 1))) THEN
pv = 0.0_dp
ELSE
pv = 10._dp*(ABS(pv)-ABS(atom%orbitals%tpsir0(k, 1)))
END IF
pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv
ELSE
pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv*100._dp
END IF
afun = afun+pv
END IF
IF (atom%orbitals%wpsir0(k, 2) > 0._dp) THEN
CALL atom_wfnr0(pv, atom%orbitals%wfnb(:, k, l), atom%basis)
pv = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv
IF (ABS(atom%orbitals%tpsir0(k, 2)) > 1.e-14_dp) THEN
IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 2))) THEN
pv = 0.0_dp
ELSE
pv = 10._dp*(ABS(pv)-ABS(atom%orbitals%tpsir0(k, 2)))
END IF
pv = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv
ELSE
pv = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv*100._dp
END IF
afun = afun+pv
END IF
END IF

1213
src/atom_grb.F Normal file

File diff suppressed because it is too large Load diff

View file

@ -189,8 +189,9 @@ MODULE atom_types
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: refene, refchg, refnod
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wrefene, wrefchg, wrefnod
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: crefene, crefchg, crefnod
REAL(KIND=dp), DIMENSION(:, :), POINTER :: wpsir0
REAL(KIND=dp), DIMENSION(:, :), POINTER :: wpsir0, tpsir0
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rcmax
CHARACTER(LEN=2), DIMENSION(:, :, :), POINTER :: reftype
END TYPE atom_orbitals
!> \brief Operator matrices
@ -1011,8 +1012,11 @@ CONTAINS
orbs%crefnod = 0._dp
ALLOCATE (orbs%rcmax(mo, 0:3, 2))
orbs%rcmax = 0._dp
ALLOCATE (orbs%wpsir0(mo, 2))
ALLOCATE (orbs%wpsir0(mo, 2), orbs%tpsir0(mo, 2))
orbs%wpsir0 = 0._dp
orbs%tpsir0 = 0._dp
ALLOCATE (orbs%reftype(mo, 0:3, 2))
orbs%reftype = "XX"
END SUBROUTINE create_atom_orbs
@ -1070,6 +1074,12 @@ CONTAINS
IF (ASSOCIATED(orbs%wpsir0)) THEN
DEALLOCATE (orbs%wpsir0)
END IF
IF (ASSOCIATED(orbs%tpsir0)) THEN
DEALLOCATE (orbs%tpsir0)
END IF
IF (ASSOCIATED(orbs%reftype)) THEN
DEALLOCATE (orbs%reftype)
END IF
DEALLOCATE (orbs)

View file

@ -55,7 +55,7 @@ MODULE atom_utils
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_utils'
PUBLIC :: atom_condnumber, atom_completeness
PUBLIC :: atom_condnumber, atom_completeness, atom_basis_condnum
PUBLIC :: atom_set_occupation, get_maxl_occ, get_maxn_occ
PUBLIC :: atom_denmat, atom_density, atom_core_density
PUBLIC :: integrate_grid, atom_trace, atom_solve
@ -2014,7 +2014,7 @@ CONTAINS
CALL init_spherical_harmonics(3, 0)
cnum = 0.0_dp
DO i = 1, 9
ci = 2.0_dp*(0.8_dp+i*0.1_dp)
ci = 2.0_dp*(0.85_dp+i*0.05_dp)
rad(i) = crad*ci
CALL atom_basis_condnum(basis, rad(i), cnum(i))
WRITE (iw, '(A,F15.3,T50,A,F14.4)') " Lattice constant:", &
@ -2169,7 +2169,7 @@ CONTAINS
!> \param cnum ...
! **************************************************************************************************
SUBROUTINE atom_basis_condnum(basis, rad, cnum)
TYPE(atom_basis_type), POINTER :: basis
TYPE(atom_basis_type) :: basis
REAL(KIND=dp), INTENT(IN) :: rad
REAL(KIND=dp), INTENT(OUT) :: cnum

View file

@ -304,7 +304,7 @@ CONTAINS
! dump grid density and xcpot (xc energy?)
IF (.FALSE.) THEN
CALL open_file(unit_number=unit_nr, file_name="atom.dat")
CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atom.dat")
DO i = 1, SIZE(xcpot(:, 1, 1))
WRITE (unit_nr, *) atom%basis%grid%rad(i), rho(i, 1), xcpot(i, 1, 1)
ENDDO
@ -315,12 +315,26 @@ CONTAINS
deriv => xc_dset_get_derivative(deriv_set, "(rho)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 0)
IF (.FALSE.) THEN
CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atompot.dat")
DO i = 1, SIZE(xcpot(:, 1, 1))
WRITE (unit_nr, *) atom%basis%grid%rad(i), rho(i, 1), xcpot(i, 1, 1)
ENDDO
CALL close_file(unit_nr)
ENDIF
DEALLOCATE (rho)
END IF
IF (needs%norm_drho) THEN
deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 1)
IF (.FALSE.) THEN
CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atomdpot.dat")
DO i = 1, SIZE(xcpot(:, 1, 1))
WRITE (unit_nr, *) atom%basis%grid%rad(i), drho(i, 1), xcpot(i, 1, 1)
ENDDO
CALL close_file(unit_nr)
ENDIF
DEALLOCATE (drho)
END IF
IF (needs%tau) THEN

View file

@ -205,30 +205,35 @@ CONTAINS
NULLIFY (print_key, keyword)
! Print key section
CALL cp_print_key_section_create(print_key, "PROGRAM_BANNER", &
description="Controls the printing of the banner of the ATOM program", &
print_level=silent_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "METHOD_INFO", &
description="Controls the printing of method information", &
print_level=medium_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "BASIS_SET", &
description="Controls the printing of the basis sets", &
print_level=high_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "POTENTIAL", &
description="Controls the printing of the potentials", &
print_level=high_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create( &
print_key, "FIT_DENSITY", &
description="Fit the total electronic density to a linear combination of Gaussian functions", &
@ -242,8 +247,10 @@ CONTAINS
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "FIT_KGPOT", &
description="Fit an approximation to the non-additive kinetic energy potential used in KG", &
description="Fit an approximation to the non-additive"// &
" kinetic energy potential used in KG", &
print_level=high_print_level, filename="__STD_OUT__")
CALL keyword_create(keyword, name="NUM_GAUSSIAN", &
description="Number of Gaussian terms for the fit", &
@ -257,10 +264,10 @@ CONTAINS
default_i_val=4)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "RESPONSE_BASIS", &
description="Calculate a response basis set contraction scheme", &
print_level=high_print_level, filename="__STD_OUT__")
@ -279,18 +286,114 @@ CONTAINS
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "GEOMETRICAL_RESPONSE_BASIS", &
description="Calculate a response basis set based on a set of geometrical exponents", &
print_level=high_print_level, filename="__STD_OUT__")
!
CALL keyword_create(keyword, name="DELTA_CHARGE", &
description="Variation of charge used in finite difference calculation", &
usage="DELTA_CHARGE real ", type_of_var=real_t, &
default_r_val=0.05_dp)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="DERIVATIVES", &
description="Number of wavefunction derivatives to calculate", &
usage="DERIVATIVES integer ", type_of_var=integer_t, &
default_i_val=3)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="QUADRATURE", &
description="Algorithm to construct the atomic radial grids", &
usage="QUADRATURE (GC_SIMPLE|GC_TRANSFORMED|GC_LOG)", &
enum_c_vals=s2a("GC_SIMPLE", "GC_TRANSFORMED", "GC_LOG"), &
enum_i_vals=(/do_gapw_gcs, do_gapw_gct, do_gapw_log/), &
enum_desc=s2a("Gauss-Chebyshev quadrature", &
"Transformed Gauss-Chebyshev quadrature", &
"Logarithmic transformed Gauss-Chebyshev quadrature"), &
default_i_val=do_gapw_log)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="GRID_POINTS", &
description="Number of radial grid points", &
usage="GRID_POINTS integer", &
default_i_val=400)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="NUM_GTO_CORE", &
description="Number of Gaussian type functions for s, p, d, ... "// &
"for the main body of the basis", &
usage="NUM_GTO 6 ", n_var=1, type_of_var=integer_t, &
default_i_val=-1)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="NUM_GTO_EXTENDED", &
description="Number of Gaussian type functions for s, p, d, ... "// &
"for the extension set", &
usage="NUM_GTO 4 ", n_var=1, type_of_var=integer_t, &
default_i_val=-1)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="NUM_GTO_POLARIZATION", &
description="Number of Gaussian type functions for the polarization set", &
usage="NUM_GTO 4 ", n_var=1, type_of_var=integer_t, &
default_i_val=-1)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="EXTENSION_BASIS", &
description="Number of basis functions for s, p, d, ... "// &
"for the extension set", &
usage="EXTENSION_BASIS 4 3 2 1 ", n_var=-1, type_of_var=integer_t, &
default_i_val=-1)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="GEOMETRICAL_FACTOR", &
description="Geometrical basis: factor C in a*C^k (initial value for optimization)", &
usage="GEOMETRICAL_FACTOR real", &
default_r_val=2.3_dp)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="GEO_START_VALUE", &
description="Geometrical basis: starting value a in a*C^k (initial value for optimization)", &
usage="GEO_START_VALUE real", &
default_r_val=0.06_dp)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="CONFINEMENT", &
description="Onset value of barrier confinement potential [Bohr]", &
usage="CONFINEMENT real", &
default_r_val=8.00_dp)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="NAME_BODY", &
description="Specifies the body of the basis set name ", &
usage="NAME_BODY <char>", &
type_of_var=char_t, default_c_val="GRB", n_var=-1)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
!
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "SCF_INFO", &
description="Controls the printing of SCF information", &
print_level=medium_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "ORBITALS", &
description="Controls the printing of the optimized orbitals information", &
print_level=high_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "ANALYZE_BASIS", &
description="Calculates some basis set analysis data", &
print_level=high_print_level, filename="__STD_OUT__")
@ -307,24 +410,28 @@ CONTAINS
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "FIT_PSEUDO", &
description="Controls the printing of FIT PSEUDO task", &
print_level=medium_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "FIT_BASIS", &
description="Controls the printing of FIT BASIS task", &
print_level=medium_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "UPF_FILE", &
description="Write GTH pseudopotential in UPF format", &
print_level=high_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! Print key section
CALL cp_print_key_section_create(print_key, "ENERGIES_MINUS_KINETIC", &
description="Print out the total energy and orbital energies without "// &
"the kinetic energy component. Useful for atomic calculations used"// &
@ -333,7 +440,7 @@ CONTAINS
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! ADMM Analysis
! Print key section: ADMM Analysis
CALL cp_print_key_section_create(print_key, "ADMM", &
description="Analysis of ADMM approximation to exact exchange", &
print_level=high_print_level, filename="__STD_OUT__")
@ -345,16 +452,6 @@ CONTAINS
CALL atom_basis_section(subsection)
CALL section_add_subsection(print_key, subsection)
CALL section_release(subsection)
! CALL keyword_create(keyword, name="OVERLAP_CONDITION_NUMBER", &
! description="Condition number of the basis set overlap matrix calculated for a cubic crystal", &
! usage="OVERLAP_CONDITION_NUMBER <logical>", type_of_var=logical_t, default_l_val=.FALSE.)
! CALL section_add_keyword(print_key, keyword)
! CALL keyword_release(keyword)
! CALL keyword_create(keyword, name="COMPLETENESS", &
! description="Calculate a completeness estimate for the basis set.", &
! usage="COMPLETENESS <logical>", type_of_var=logical_t, default_l_val=.FALSE.)
! CALL section_add_keyword(print_key, keyword)
! CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
@ -1093,6 +1190,14 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="TARGET_PSIR0", &
description="Minimum value for the wavefunctions at r=0 (only occupied states)"// &
" Value=0 means keeping wfn(r=0)=0", &
usage="TARGET_PSIR0 0.50", &
default_r_val=0._dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="WEIGHT_PSIR0", &
description="Weight for the wavefunctions at r=0 (only occupied states)", &
usage="WEIGHT_PSIR0 0.01", &
@ -1107,6 +1212,13 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="SEMICORE_LEVEL", &
description="Energy at which to consider a full shell as semicore", &
usage="SEMICORE_LEVEL 1.0", &
default_r_val=1._dp, unit_str="hartree")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_powell_section
! **************************************************************************************************

View file

@ -0,0 +1,48 @@
&GLOBAL
PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
ELEMENT C
RUN_TYPE ENERGY
ELECTRON_CONFIGURATION CORE 2s2 2p2
CORE [He]
MAX_ANGULAR_MOMENTUM 3
&METHOD
METHOD_TYPE KOHN-SHAM
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END METHOD
&PRINT
&GEOMETRICAL_RESPONSE_BASIS
NUM_GTO_CORE 6
NUM_GTO_EXTENDED 4
NUM_GTO_POLARIZATION 4
EXTENSION_BASIS 2 1
CONFINEMENT 8.0
NAME_BODY GRB-q4
&END
&END
&POTENTIAL
PSEUDO_TYPE GTH
&GTH_POTENTIAL
2 2 0 0
0.33842037096913 2 -8.80041712690571 1.33249221738633
1
0.30268766321861 1 9.60946227253083
&END
CONFINEMENT_TYPE BARRIER
CONFINEMENT 200. 4.0 12.0
&END POTENTIAL
&POWELL
ACCURACY 1.e-14
STEP_SIZE 0.20
&END
&END ATOM

View file

@ -1,12 +1,12 @@
C-HF-q4.inp 56 1.e-06 0.0006073242
H-B3LYP-q1.inp 56 1.e-07 0.0350327794
O-B97-q6.inp 56 1.e-06 0.0002771652
O-PBE0-q6.inp 56 1e-06 0.0867371679
O-PBE0-q6.inp 56 1e-06 0.0867738353
O_KG.inp 61 5e-04 0.9385713054E-01
C_nlcc_1.inp 35 1.e-12 -6.841350380260
C_nlcc_2.inp 35 1.e-12 -6.886838336582
C_nlcc_3.inp 61 5.e-03 0.0483032344
C_nlcc_4.inp 61 2.e-05 0.0976869423
C_nlcc_3.inp 61 5.e-03 0.0487762073
C_nlcc_4.inp 61 2.e-05 0.0982881676
upf1.inp 35 1.e-12 -5.352864657647
upf2.inp 35 1.e-12 -7.123642611762
ecp1.inp 35 1.e-12 -5.362932801634
@ -14,4 +14,5 @@ ecp2.inp 35 1.e-12
C_basis1.inp 35 1.e-12 -37.747721651890
C_basis2.inp 35 1.e-12 -37.728835366125
C_basis3.inp 35 1.e-12 -13.923782735900
C_grb.inp 35 1.e-12 -5.360947119808
#EOF