1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!  This module constains the following active subroutines:
3!   InitOrbit, DestroyOrbit, CopyOrbit, InitFC, DestroyFC,
4!      InitPot, DestroyPot, CopyPot, InitSCF, CopySCF
5!
6!  Note that all energies including tau are in Rydberg units
7!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
8
9#if defined HAVE_CONFIG_H
10#include "config.h"
11#endif
12
13MODULE atomdata
14
15  IMPLICIT NONE
16
17  TYPE OrbitInfo
18     CHARACTER(132) :: exctype
19     INTEGER :: nps, npp, npd ,npf, npg, norbit
20     INTEGER, POINTER :: np(:) => null()
21     INTEGER, POINTER :: l(:) => null()
22     INTEGER, POINTER :: kappa(:) => null()
23     REAL(8), POINTER :: eig(:) => null()
24     REAL(8), POINTER :: occ(:) => null()
25     REAL(8), POINTER :: wfn(:,:) => null()
26     REAL(8), POINTER :: lwfn(:,:) => null()
27     REAL(8), POINTER :: otau(:,:) => null() ! kinetic energy density for orbital
28     REAL(8), POINTER :: lqp(:,:) => null()  ! only used for HF
29     REAL(8), POINTER :: X(:,:) => null()    ! identical to HF%SumY(:,:)
30     LOGICAL, POINTER :: iscore(:) => null()
31     REAL(8),POINTER :: den(:) => null() ! accumulated over states
32     REAL(8),POINTER :: tau(:) => null() ! accumulated over states
33     REAL(8),POINTER :: deltatau(:) => null() !tau-tauW   (tauW==Weizsacker)
34                                      ! note this is evaluated with 1/4pir^2
35  END TYPE OrbitInfo
36
37  TYPE FCinfo
38     REAL(8), POINTER :: coreden(:) => null()
39     REAL(8), POINTER :: valeden(:) => null()
40     REAL(8), POINTER :: coretau(:) => null()
41     REAL(8), POINTER :: valetau(:) => null()
42     REAL(8) :: zvale,zcore
43  END TYPE FCinfo
44
45  TYPE PotentialInfo
46     CHARACTER(2) :: sym
47     INTEGER :: nz     !  nz is nuclear charge
48     REAL(8) :: zz        !  zz=nz is nuclear charge
49     REAL(8) :: q,v0,v0p  !  q is total electron charge
50     !  v0,v0p are potential value and deriv at r=0
51     REAL(8) :: Nv0,Nv0p    !  finite nucleus value and deriv at 0
52     REAL(8), POINTER :: rv(:)  => null()
53     REAL(8), POINTER :: rvn(:) => null()
54     REAL(8), POINTER :: rvh(:) => null()
55     REAL(8), POINTER :: rvx(:) => null()
56     !  rv(n) is  veff * r
57     !  rvh is hartree potential for den
58     !  rvn is nuclear potential
59     !  rvx is exchange-correlation potential
60     LOGICAL :: needvtau
61     REAL(8), POINTER :: vtau(:) => null() !for meta-gga
62     INTEGER :: finitenucleusmodel
63     ! Based on models 2, 3, 4, 5 discussed by Dirk Anrae ,
64     !   Physics Reports 336 (2000) 413-525
65     !    default is 0 for previous Gaussian model
66     !    for finitenucleusmodel<0, finite nucleus is false
67  END TYPE PotentialInfo
68
69  TYPE SCFInfo
70     INTEGER :: iter
71     REAL(8) :: delta,eone,ekin,estatic,ecoul,eexc,oepcs,etot
72     REAL(8) :: valekin,valecoul,valeexc,corekin,evale ! used in frozencore only
73  END TYPE SCFInfo
74
75  LOGICAL :: frozencorecalculation
76  LOGICAL :: setupfrozencore
77  LOGICAL :: scalarrelativistic
78  LOGICAL :: diracrelativistic
79  LOGICAL :: BDsolve
80  LOGICAL :: finitenucleus
81  LOGICAL :: gaussianshapefunction,besselshapefunction
82  LOGICAL :: ColleSalvetti
83  LOGICAL :: HFpostprocess
84  LOGICAL :: localizedcoreexchange
85  LOGICAL :: exploremode
86  INTEGER :: nlogderiv
87  REAL(8) :: minlogderiv,maxlogderiv
88
89
90CONTAINS
91
92!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
93!   Subroutine InitOrbit  -- used in CopyOrbit
94!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
95
96  SUBROUTINE InitOrbit(Orbit,norbit,n,exctype)
97    TYPE (OrbitInfo), INTENT(INOUT) :: Orbit
98    INTEGER, INTENT(IN) :: n,norbit
99    CHARACTER(*),INTENT(IN) :: exctype
100    INTEGER :: ok
101    CALL DestroyOrbit(Orbit)
102    Orbit%norbit=norbit;Orbit%exctype=TRIM(exctype)
103    Orbit%nps=0;Orbit%npp=0;Orbit%npd=0;Orbit%npf=0;Orbit%npg=0
104    ALLOCATE(Orbit%np(norbit),Orbit%l(norbit),Orbit%eig(norbit),&
105&            Orbit%occ(norbit),Orbit%iscore(norbit),&
106&            stat=ok)
107    IF (ok/=0) STOP 'Error in allocation of nl, l, occ...'
108    Orbit%iscore=.false.
109    Orbit%np=0;Orbit%l=0
110    Orbit%eig=0.d0;Orbit%occ=0.d0
111    ALLOCATE(Orbit%wfn(n,norbit),Orbit%otau(n,norbit), &
112&             Orbit%den(n),Orbit%tau(n),Orbit%deltatau(n),stat=ok)
113    IF (ok/=0) STOP 'Error in allocation of wfn, den...'
114    Orbit%wfn=0.d0;Orbit%den=0.d0;Orbit%tau=0.d0;Orbit%otau=0.d0
115    Orbit%deltatau=0.d0
116    If (diracrelativistic) then
117       ALLOCATE(Orbit%lwfn(n,norbit),Orbit%kappa(norbit),stat=ok)
118       IF (ok/=0) STOP 'Error in allocation of lwfn,kappa'
119       Orbit%lwfn=0.d0
120       Orbit%kappa=0.d0
121    else
122       NULLIFY(Orbit%lwfn,Orbit%kappa)
123    Endif
124    If (exctype == "HF".or.exctype == "EXXKLI") then
125       ALLOCATE(Orbit%lqp(norbit,norbit),Orbit%X(n,norbit),stat=ok)
126       IF (ok/=0) STOP 'Error in allocation of lqp, X...'
127    ELSE
128       NULLIFY(Orbit%lqp,Orbit%X)
129    ENDIF
130  END SUBROUTINE InitOrbit
131
132  SUBROUTINE DestroyOrbit(Orbit)
133    TYPE (OrbitInfo), INTENT(INOUT) :: Orbit
134    IF (ASSOCIATED(Orbit%np)) DEALLOCATE(Orbit%np)
135    IF (ASSOCIATED(Orbit%l)) DEALLOCATE(Orbit%l)
136    IF (ASSOCIATED(Orbit%kappa)) DEALLOCATE(Orbit%kappa)
137    IF (ASSOCIATED(Orbit%iscore)) DEALLOCATE(Orbit%iscore)
138    IF (ASSOCIATED(Orbit%eig)) DEALLOCATE(Orbit%eig)
139    IF (ASSOCIATED(Orbit%occ)) DEALLOCATE(Orbit%occ)
140    IF (ASSOCIATED(Orbit%wfn)) DEALLOCATE(Orbit%wfn)
141    IF (ASSOCIATED(Orbit%otau)) DEALLOCATE(Orbit%otau)
142    IF (ASSOCIATED(Orbit%lwfn)) DEALLOCATE(Orbit%lwfn)
143    IF (ASSOCIATED(Orbit%den)) DEALLOCATE(Orbit%den)
144    IF (ASSOCIATED(Orbit%tau)) DEALLOCATE(Orbit%tau)
145    IF (ASSOCIATED(Orbit%deltatau)) DEALLOCATE(Orbit%deltatau)
146    IF (ASSOCIATED(Orbit%lqp)) DEALLOCATE(Orbit%lqp)
147    IF (ASSOCIATED(Orbit%X)) DEALLOCATE(Orbit%X)
148  END SUBROUTINE DestroyOrbit
149
150!!!!!!!!!!!!!!!!!!!!!!!!!
151!  CopyOrbit(source,copy)
152!!!!!!!!!!!!!!!!!!!!!!!!!
153  SUBROUTINE CopyOrbit(SOrbit,COrbit)
154    TYPE(OrbitInfo),INTENT(INOUT)::SOrbit
155    TYPE(OrbitInfo),INTENT(INOUT)::COrbit
156    INTEGER::n
157    n=SIZE(SOrbit%den,1)
158    CALL InitOrbit(COrbit,SOrbit%norbit,n,SOrbit%exctype)
159    COrbit%nps=SOrbit%nps
160    COrbit%npp=SOrbit%npp
161    COrbit%npd=SOrbit%npd
162    COrbit%npf=SOrbit%npf
163    COrbit%npg=SOrbit%npg
164    COrbit%np(1:SOrbit%norbit)=SOrbit%np(1:SOrbit%norbit)
165    COrbit%l(1:SOrbit%norbit)=SOrbit%l(1:SOrbit%norbit)
166    COrbit%eig(1:SOrbit%norbit)=SOrbit%eig(1:SOrbit%norbit)
167    COrbit%occ(1:SOrbit%norbit)=SOrbit%occ(1:SOrbit%norbit)
168    COrbit%wfn(:,1:SOrbit%norbit)=SOrbit%wfn(:,1:SOrbit%norbit)
169    COrbit%otau(:,1:SOrbit%norbit)=SOrbit%otau(:,1:SOrbit%norbit)
170    COrbit%iscore(1:SOrbit%norbit)=SOrbit%iscore(1:SOrbit%norbit)
171    COrbit%den=SOrbit%den
172    COrbit%tau=SOrbit%tau
173    COrbit%deltatau=SOrbit%deltatau
174    If (diracrelativistic) then
175       COrbit%lwfn(:,1:SOrbit%norbit)=SOrbit%lwfn(:,1:SOrbit%norbit)
176       COrbit%kappa(1:SOrbit%norbit)=SOrbit%kappa(1:SOrbit%norbit)
177    Endif
178
179    IF (SOrbit%exctype == "HF".or.SOrbit%exctype == "EXXKLI") then
180        COrbit%X(:,1:SOrbit%norbit)=SOrbit%X(:,1:SOrbit%norbit)
181        COrbit%lqp(1:SOrbit%norbit,1:SOrbit%norbit)= &
182&                SOrbit%lqp(1:SOrbit%norbit,1:SOrbit%norbit)
183    ENDIF
184  END SUBROUTINE CopyOrbit
185
186  SUBROUTINE InitFC(FC,n)
187    INTEGER, INTENT(IN) :: n
188    TYPE (FCInfo), INTENT(INOUT) :: FC
189    INTEGER :: ok
190    CALL DestroyFC(FC)
191    FC%zvale=0.d0;FC%zcore=0.d0
192    ALLOCATE(FC%coreden(n),FC%valeden(n),stat=ok)
193    IF (ok/=0) STOP 'Error in allocation of coreden, valeden,...'
194    FC%coreden=0.d0;FC%valeden=0.d0
195    ALLOCATE(FC%coretau(n),FC%valetau(n),stat=ok)
196    IF (ok/=0) STOP 'Error in allocation of coretau, valetau,...'
197    FC%coretau=0.d0;FC%valetau=0.d0
198  END SUBROUTINE InitFC
199
200  SUBROUTINE DestroyFC(FC)
201    TYPE (FCInfo), INTENT(INOUT) :: FC
202    IF (ASSOCIATED(FC%valeden)) DEALLOCATE(FC%valeden)
203    IF (ASSOCIATED(FC%coreden)) DEALLOCATE(FC%coreden)
204    IF (ASSOCIATED(FC%valetau)) DEALLOCATE(FC%valetau)
205    IF (ASSOCIATED(FC%coretau)) DEALLOCATE(FC%coretau)
206  END SUBROUTINE DestroyFC
207
208!!!!!!!!!!!!!!!!!!!!!!!!!
209!  CopyFC(source,copy)
210!!!!!!!!!!!!!!!!!!!!!!!!!
211  SUBROUTINE CopyFC(SFC,CFC)
212    TYPE(FCInfo),INTENT(IN) :: SFC
213    TYPE(FCInfo),INTENT(INOUT) :: CFC
214    INTEGER :: n
215    n=SIZE(SFC%coreden,1)
216    CALL InitFC(CFC,n)
217    CFC%zvale=SFC%zvale
218    CFC%zcore=SFC%zcore
219    CFC%coreden(1:n)=SFC%coreden(1:n)
220    CFC%valeden(1:n)=SFC%valeden(1:n)
221    CFC%coretau(1:n)=SFC%coretau(1:n)
222    CFC%valetau(1:n)=SFC%valetau(1:n)
223  END SUBROUTINE CopyFC
224
225  SUBROUTINE InitPot(Pot,n)
226    INTEGER, INTENT(IN) :: n
227    TYPE (PotentialInfo), INTENT(INOUT) :: Pot
228    INTEGER :: ok
229    CALL DestroyPot(Pot)
230!   Pot%sym="";Pot%nz=0;Pot%zz=0.d0;Pot%q=0.d0;Pot%v0=0.d0;Pot%v0p=0.d0
231    ALLOCATE(Pot%rv(n),Pot%rvn(n),Pot%rvh(n),Pot%rvx(n),Pot%vtau(n),stat=ok)
232    IF (ok/=0) STOP 'Error in allocation of Pot%rv, Pot%rvh...'
233    Pot%rv=0.d0;Pot%rvn=0.d0;Pot%rvh=0.d0;Pot%rvx=0.d0;Pot%vtau=0.d0
234    Pot%needvtau=.false.
235  END SUBROUTINE InitPot
236
237  SUBROUTINE DestroyPot(Pot)
238    TYPE (PotentialInfo), INTENT(INOUT) :: Pot
239    IF (ASSOCIATED(Pot%rv)) DEALLOCATE(Pot%rv)
240    IF (ASSOCIATED(Pot%rvn)) DEALLOCATE(Pot%rvn)
241    IF (ASSOCIATED(Pot%rvh)) DEALLOCATE(Pot%rvh)
242    IF (ASSOCIATED(Pot%rvx)) DEALLOCATE(Pot%rvx)
243    IF (ASSOCIATED(Pot%vtau)) DEALLOCATE(Pot%vtau)
244  END SUBROUTINE DestroyPot
245
246!!!!!!!!!!!!!!!!!!!!!!!!!
247!  CopyPot(source,copy)
248!!!!!!!!!!!!!!!!!!!!!!!!!
249  SUBROUTINE CopyPot(SPot,CPot)
250    TYPE(PotentialInfo),INTENT(IN) :: SPot
251    TYPE(PotentialInfo),INTENT(INOUT) :: CPot
252    INTEGER :: n
253    n=SIZE(SPot%rv,1)
254    CALL InitPot(CPot,n)
255    CPot%nz=SPot%nz
256    CPot%zz=SPot%zz
257    CPot%sym=SPot%sym
258    CPot%q=SPot%q
259    CPot%v0=SPot%v0
260    CPot%v0p=SPot%v0p
261    CPot%finitenucleusmodel=SPot%finitenucleusmodel
262    CPot%Nv0=SPot%Nv0
263    CPot%Nv0p=SPot%Nv0p
264    CPot%rv(1:n)=SPot%rv(1:n)
265    CPot%rvn(1:n)=SPot%rvn(1:n)
266    CPot%rvh(1:n)=SPot%rvh(1:n)
267    CPot%rvx(1:n)=SPot%rvx(1:n)
268    CPot%vtau(1:n)=SPot%vtau(1:n)
269    CPot%needvtau=SPot%needvtau
270  END SUBROUTINE CopyPot
271
272  SUBROUTINE InitSCF(SCF)
273    TYPE(SCFInfo),INTENT(INOUT)::SCF
274    SCF%iter=0
275    SCF%delta=0.d0;SCF%eone=0.d0;SCF%ekin=0.d0;SCF%estatic=0.d0
276    SCF%ecoul=0.d0;SCF%eexc=0.d0;SCF%oepcs=0.d0;SCF%etot=0.d0
277    SCF%valekin=0.d0;SCF%valecoul=0.d0;SCF%valeexc=0.d0
278    SCF%corekin=0.d0;SCF%evale=0.d0
279  END SUBROUTINE InitSCF
280
281!!!!!!!!!!!!!!!!!!!!!!!!!
282!  CopySCF(source,copy)
283!!!!!!!!!!!!!!!!!!!!!!!!!
284  SUBROUTINE CopySCF(SSCF,CSCF)
285    TYPE(SCFInfo),INTENT(IN)::SSCF
286    TYPE(SCFInfo),INTENT(INOUT)::CSCF
287    CALL InitSCF(CSCF)
288    CSCF%iter=SSCF%iter
289    CSCF%delta=SSCF%delta
290    CSCF%eone=SSCF%eone
291    CSCF%ekin=SSCF%ekin
292    CSCF%ecoul=SSCF%ecoul
293    CSCF%estatic=SSCF%estatic
294    CSCF%eexc=SSCF%eexc
295    CSCF%etot=SSCF%etot
296    CSCF%valekin=SSCF%valekin
297    CSCF%valecoul=SSCF%valecoul
298    CSCF%valeexc=SSCF%valeexc
299    CSCF%corekin=SSCF%corekin
300    CSCF%evale=SSCF%evale
301  END SUBROUTINE CopySCF
302
303END MODULE atomdata
304