Line data Source code
1 : MODULE m_cored
2 : CONTAINS
3 0 : SUBROUTINE cored(input, jspin, atoms, rho, sphhar, l_CoreDenPresent, vr, qint, rhc, tec, seig, EnergyDen)
4 : ! *******************************************************
5 : ! ***** set up the core densities for compounds. *****
6 : ! ***** d.d.koelling *****
7 : ! *******************************************************
8 : USE m_juDFT
9 : USE m_intgr, ONLY : intgr3,intgr0,intgr1
10 : USE m_constants
11 : !USE m_setcor
12 : USE m_differ
13 : USE m_types
14 : USE m_cdn_io
15 : USE m_xmlOutput
16 : IMPLICIT NONE
17 :
18 : TYPE(t_input),INTENT(IN) :: input
19 : TYPE(t_sphhar),INTENT(IN) :: sphhar
20 : TYPE(t_atoms),INTENT(IN) :: atoms
21 : !
22 : ! .. Scalar Arguments ..
23 : INTEGER, INTENT (IN) :: jspin
24 : LOGICAL, INTENT (IN) :: l_CoreDenPresent
25 : REAL, INTENT (OUT) :: seig
26 : ! ..
27 : ! .. Array Arguments ..
28 : REAL, INTENT(IN) :: vr(atoms%jmtd,atoms%ntype)
29 : REAL, INTENT(INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
30 : REAL, INTENT(INOUT) :: rhc(atoms%msh,atoms%ntype,input%jspins)
31 : REAL, INTENT(INOUT) :: qint(atoms%ntype,input%jspins)
32 : REAL, INTENT(INOUT) :: tec(atoms%ntype,input%jspins)
33 : REAL, INTENT(INOUT), OPTIONAL :: EnergyDen(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
34 : ! ..
35 : ! .. Local Scalars ..
36 : REAL eig,fj,fl,fn,qOutside,rad,rhos,rhs,sea,sume,t2
37 : REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight, aux_weight
38 : INTEGER i,j,jatom,korb,n,ncmsh,nm,nm1,nst ,l,ierr
39 : ! ..
40 : ! .. Local Arrays ..
41 :
42 527 : REAL rhcs(atoms%msh),rhoc(atoms%msh),rhoss(atoms%msh),vrd(atoms%msh),f(0:3)
43 527 : REAL rhcs_aux(atoms%msh), rhoss_aux(atoms%msh) !> quantities for energy density calculations
44 1441 : REAL occ(maxval(atoms%econf%num_states)),a(atoms%msh),b(atoms%msh),ain(atoms%msh),ahelp(atoms%msh)
45 1441 : REAL occ_h(maxval(atoms%econf%num_states),2)
46 2355 : INTEGER kappa(maxval(atoms%econf%num_states)),nprnc(maxval(atoms%econf%num_states))
47 : CHARACTER(LEN=20) :: attributes(6)
48 : REAL stateEnergies(29)
49 : ! ..
50 :
51 527 : c = c_light(1.0)
52 527 : seig = 0.
53 : !
54 527 : IF (input%frcor.and. l_CoreDenPresent) THEN
55 0 : DO n = 1,atoms%ntype
56 0 : rnot = atoms%rmsh(1,n) ; dxx = atoms%dx(n)
57 0 : ncmsh = NINT( LOG( (atoms%rmt(n)+10.0)/rnot ) / dxx + 1 )
58 0 : ncmsh = MIN( ncmsh, atoms%msh )
59 : ! ---> update spherical charge density
60 0 : DO i = 1,atoms%jri(n)
61 0 : rhoc(i) = rhc(i,n,jspin)
62 0 : rho(i,0,n,jspin) = rho(i,0,n,jspin) + rhoc(i)/sfp_const
63 : ENDDO
64 : ! ---> for total energy calculations, determine the sum of the
65 : ! ---> eigenvalues by requiring that the core kinetic energy
66 : ! ---> remains constant.
67 0 : DO i = 1,atoms%jri(n)
68 0 : rhoc(i) = rhoc(i)*vr(i,n)/atoms%rmsh(i,n)
69 : ENDDO
70 0 : nm = atoms%jri(n)
71 0 : CALL intgr3(rhoc,atoms%rmsh(1,n),atoms%dx(n),nm,rhos)
72 0 : sea = tec(n,jspin) + rhos
73 0 : WRITE (oUnit,FMT=8030) n,jspin,tec(n,jspin),sea
74 0 : seig = seig + atoms%neq(n)*sea
75 : ENDDO
76 : RETURN
77 : END IF
78 :
79 : ! ---> set up densities
80 1441 : DO jatom = 1,atoms%ntype
81 914 : sume = 0.
82 914 : z = atoms%zatom(jatom)
83 : ! rn = rmt(jatom)
84 914 : dxx = atoms%dx(jatom)
85 914 : bmu = 0.0
86 : !CALL setcor(jatom,input%jspins,atoms,input,bmu,nst,kappa,nprnc,occ_h)
87 914 : CALL atoms%econf(jatom)%get_core(nst,nprnc,kappa,occ_h)
88 :
89 :
90 5912 : occ(1:nst) = occ_h(1:nst,jspin)
91 :
92 914 : rnot = atoms%rmsh(1,jatom)
93 914 : d = EXP(atoms%dx(jatom))
94 914 : ncmsh = NINT( LOG( (atoms%rmt(jatom)+10.0)/rnot ) / dxx + 1 )
95 914 : ncmsh = MIN( ncmsh, atoms%msh )
96 914 : rn = rnot* (d** (ncmsh-1))
97 914 : WRITE (oUnit,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
98 625082 : DO j = 1,atoms%jri(jatom)
99 624168 : rhoss(j) = 0.0
100 624168 : if(present(EnergyDen)) rhoss_aux(j) = 0.0
101 625082 : vrd(j) = vr(j,jatom)
102 : ENDDO
103 : !
104 914 : IF (input%l_core_confpot) THEN
105 : !---> linear extension of the potential with slope t1 / a.u.
106 914 : t1=0.125
107 : t1 = MAX( (vrd(atoms%jri(jatom)) - vrd(atoms%jri(jatom)-1)*d)*&
108 914 : d / (atoms%rmt(jatom)**2 * (d-1) ) , t1)
109 914 : t2=vrd(atoms%jri(jatom))/atoms%rmt(jatom)-atoms%rmt(jatom)*t1
110 914 : rr = atoms%rmt(jatom)
111 : ELSE
112 0 : t2 = vrd(atoms%jri(jatom)) / ( atoms%jri(jatom) - ncmsh )
113 : ENDIF
114 914 : IF ( atoms%jri(jatom) < ncmsh) THEN
115 90554 : DO i = atoms%jri(jatom) + 1,ncmsh
116 89640 : rhoss(i) = 0.
117 89640 : if(present(EnergyDen)) rhoss_aux(i) = 0.0
118 90554 : IF (input%l_core_confpot) THEN
119 89640 : rr = d*rr
120 89640 : vrd(i) = rr*( t2 + rr*t1 )
121 : ! vrd(i) = 2*vrd(jri(jatom)) - rr*( t2 + rr*t1 )
122 : ELSE
123 0 : vrd(i) = vrd(atoms%jri(jatom)) + t2* (i-atoms%jri(jatom))
124 : ENDIF
125 : !
126 : ENDDO
127 : END IF
128 :
129 914 : nst = atoms%econf(jatom)%num_core_states ! for lda+U
130 :
131 914 : IF (input%gw==1 .OR. input%gw==3)&
132 0 : & WRITE(15) nst,atoms%rmsh(1:atoms%jri(jatom),jatom)
133 :
134 914 : stateEnergies = 0.0
135 5912 : DO korb = 1,nst
136 5912 : IF (occ(korb) /= 0.0) THEN
137 4998 : fn = nprnc(korb)
138 4998 : fj = iabs(kappa(korb)) - .5e0
139 :
140 : ! weight = 2*fj + 1.e0
141 : ! IF (bmu > 99.) weight = occ(korb)
142 4998 : weight = 2*occ(korb)
143 :
144 4998 : fl = fj + (.5e0)*isign(1,kappa(korb))
145 :
146 4998 : eig = -2* (z/ (fn+fl))**2
147 :
148 4998 : CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, eig, a,b,ierr)
149 4998 : stateEnergies(korb) = eig
150 4998 : WRITE (oUnit,FMT=8010) fn,fl,fj,eig,weight
151 :
152 4998 : IF (ierr/=0) CALL juDFT_error("error in core-level routine" ,calledby ="cored")
153 4998 : IF (input%gw==1 .OR. input%gw==3) WRITE (15) NINT(fl),weight,eig,&
154 0 : a(1:atoms%jri(jatom)),b(1:atoms%jri(jatom))
155 :
156 4998 : sume = sume + weight*eig/input%jspins
157 4076324 : DO j = 1,ncmsh
158 4071326 : rhcs(j) = weight* (a(j)**2+b(j)**2)
159 4076324 : rhoss(j) = rhoss(j) + rhcs(j)
160 : ENDDO
161 :
162 4998 : IF(present(EnergyDen)) THEN
163 : !rhoss_aux = rhoss
164 0 : DO j = 1,ncmsh
165 : ! for energy density we want to multiply the weights
166 : ! with the eigenenergies
167 0 : rhoss_aux(j) = rhoss_aux(j) + (rhcs(j) * eig)
168 : ENDDO
169 : ENDIF
170 : ENDIF
171 : ENDDO
172 :
173 : ! ---->update spherical charge density rho with the core density.
174 : ! ---->for spin-polarized (jspins=2), take only half the density
175 914 : nm = atoms%jri(jatom)
176 625082 : DO j = 1,nm
177 624168 : rhoc(j) = rhoss(j)/input%jspins
178 625082 : rho(j,0,jatom,jspin) = rho(j,0,jatom,jspin) + rhoc(j)/sfp_const
179 : ENDDO
180 :
181 914 : IF(present(EnergyDen)) then
182 0 : DO j = 1,nm
183 : EnergyDen(j,0,jatom,jspin) = EnergyDen(j,0,jatom,jspin) &
184 0 : + rhoss_aux(j) /(input%jspins * sfp_const)
185 : ENDDO
186 : ENDIF
187 :
188 714722 : rhc(1:ncmsh,jatom,jspin) = rhoss(1:ncmsh) / input%jspins
189 61962 : rhc(ncmsh+1:atoms%msh,jatom,jspin) = 0.0
190 :
191 914 : seig = seig + atoms%neq(jatom)*sume
192 625082 : DO i = 1,nm
193 625082 : rhoc(i) = rhoc(i)*vr(i,jatom)/atoms%rmsh(i,jatom)
194 : ENDDO
195 914 : CALL intgr3(rhoc,atoms%rmsh(1,jatom),atoms%dx(jatom),nm,rhs)
196 914 : tec(jatom,jspin) = sume - rhs
197 914 : WRITE (oUnit,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
198 :
199 : ! ---> simpson integration
200 914 : rad = atoms%rmt(jatom)
201 : ! qOutside is the charge outside a single MT sphere of the considered atom type
202 914 : qOutside = rad*rhoss(nm)/2.
203 914 : DO nm1 = nm + 1,ncmsh - 1,2
204 44533 : rad = d*rad
205 44533 : qOutside = qOutside + 2*rad*rhoss(nm1)
206 44533 : rad = d*rad
207 44533 : qOutside = qOutside + rad*rhoss(nm1+1)
208 : ENDDO
209 914 : qOutside = 2*qOutside*dxx/3
210 : !+sb
211 914 : WRITE (oUnit,FMT=8020) qOutside/input%jspins
212 : !-sb
213 914 : qint(jatom,jspin) = qOutside*atoms%neq(jatom)
214 6398 : attributes = ''
215 914 : WRITE(attributes(1),'(i0)') jatom
216 914 : WRITE(attributes(2),'(i0)') NINT(z)
217 914 : WRITE(attributes(3),'(i0)') jspin
218 914 : WRITE(attributes(4),'(f18.10)') tec(jatom,jspin)
219 914 : WRITE(attributes(5),'(f18.10)') sume
220 914 : WRITE(attributes(6),'(f9.6)') qOutside/input%jspins
221 : CALL openXMLElementForm('coreStates',(/'atomType ','atomicNumber ','spin ','kinEnergy ',&
222 : 'eigValSum ','lostElectrons'/),&
223 6398 : attributes,RESHAPE((/8,12,4,9,9,13,6,3,1,18,18,9/),(/6,2/)))
224 5912 : DO korb = 1, atoms%econf(jatom)%num_core_states
225 4998 : fj = iabs(kappa(korb)) - .5e0
226 : ! weight = 2*fj + 1.e0
227 : ! IF (bmu > 99.) weight = occ(korb)
228 4998 : weight = occ(korb)
229 4998 : fl = fj + (.5e0)*isign(1,kappa(korb))
230 34986 : attributes = ''
231 4998 : WRITE(attributes(1),'(i0)') nprnc(korb)
232 4998 : WRITE(attributes(2),'(i0)') NINT(fl)
233 4998 : WRITE(attributes(3),'(f4.1)') fj
234 4998 : WRITE(attributes(4),'(f20.10)') stateEnergies(korb)
235 4998 : WRITE(attributes(5),'(f15.10)') weight
236 : CALL writeXMLElementForm('state',(/'n ','l ','j ','energy','weight'/),&
237 30902 : attributes(1:5),RESHAPE((/1,1,1,6,6,2,2,4,20,15/),(/5,2/)))
238 : END DO
239 2355 : CALL closeXMLElement('coreStates')
240 : ENDDO
241 :
242 : RETURN
243 :
244 : 8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f9.6,5x,&
245 : & 'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,&
246 : & 'weight')
247 : 8010 FORMAT (12x,2f5.0,f6.1,f10.4,f12.4)
248 : 8020 FORMAT (f20.8,' electrons lost from core.')
249 : 8030 FORMAT (10x,'atom type',i5,' (spin',i2,') ',/,10x,&
250 : & 'kinetic energy=',e20.12,5x,'sum of the eigenvalues=',&
251 : & e20.12)
252 527 : END SUBROUTINE cored
253 : END MODULE m_cored
|
