Line data Source code
1 : MODULE m_atom2
2 : use m_juDFT
3 : ! *************************************************************
4 : ! fully relativistic atomic program based on the subroutines
5 : ! differ, outint and inwint by d.d.koelling
6 : ! erich wimmer august 1981
7 : ! modified for use in start-up generator. m.w. april 1982
8 : ! modified by adding stabilizing well. m. weinert may 1990
9 : ! *************************************************************
10 : CONTAINS
11 84 : SUBROUTINE atom2(&
12 : & atoms, xcpot, input, ntyp, jrc, rnot1,&
13 : & qdel,&
14 84 : & rhoss, nst, lnum, eig, vbar,l_valence)
15 :
16 : USE m_intgr, ONLY: intgr1, intgr0
17 : USE m_constants
18 : USE m_potl0
19 : USE m_stpot1
20 : ! USE m_setcor
21 : USE m_differ
22 : USE m_types
23 : IMPLICIT NONE
24 : ! ..
25 : ! .. Scalar Arguments ..
26 :
27 : TYPE(t_atoms), INTENT(IN) :: atoms
28 : CLASS(t_xcpot), INTENT(IN) :: xcpot
29 : TYPE(t_input), INTENT(IN) :: input
30 : INTEGER, INTENT(IN) :: jrc, ntyp
31 : REAL, INTENT(IN) :: rnot1, qdel
32 : REAL, INTENT(OUT) :: rhoss(:, :) !(mshd,input%jspins),
33 : REAL, INTENT(OUT) :: eig(29, input%jspins), vbar(input%jspins)
34 : INTEGER, INTENT(OUT) :: nst, lnum(29)
35 : LOGICAL,OPTIONAL,INTENT(IN)::l_valence
36 : ! ..
37 : ! .. Local Scalars ..
38 : REAL c, d, delrv, dist, distol, e, fisr, fj, fl, fn, h,&
39 : & p, p1, pmax, pmin, r, r3, rn, rnot, z, zero, bmu_l, rho
40 : INTEGER i, inr0, it, itmax, k, l, n, ispin, kk, ierr, msh_l,isp
41 : LOGICAL conv, lastit, l_start
42 : ! ..
43 : ! .. Local Arrays ..
44 84 : REAL a(jrc), b(jrc), dens(jrc), occ(29, input%jspins)
45 84 : REAL rad(jrc), rev(29, input%jspins), ahelp(jrc), ain(jrc),&
46 84 : & rh(jrc), vr(jrc), f(0:3),&
47 84 : & vr1(jrc, input%jspins), vr2(jrc, input%jspins), vx(atoms%msh, input%jspins), vxc(atoms%msh, input%jspins)
48 : INTEGER kappa(29), nprnc(29)
49 : ! ..
50 : ! ..
51 : ! .. Data statements ..
52 : !----> distol set from 1.0e-6 to 1.0e-3
53 : DATA zero, distol/0.0e0, 1.0e-3/
54 : ! ..
55 84 : c = c_light(1.0)
56 125212 : vxc(:, :) = 0.0
57 125212 : vx(:, :) = 0.0
58 : !
59 84 : WRITE (oUnit, FMT=8000)
60 : 8000 FORMAT(' subroutine atom2 entered')
61 84 : z = atoms%zatom(ntyp)
62 84 : n = jrc
63 84 : rnot = rnot1
64 84 : itmax = 100
65 84 : pmin = 0.01e0
66 84 : pmax = 0.2e0
67 84 : h = atoms%dx(ntyp)
68 84 : d = exp(h)
69 84 : r = rnot
70 72758 : DO i = 1, n
71 72674 : rad(i) = r
72 72758 : r = r*d
73 : enddo
74 84 : rn = rad(n)
75 84 : bmu_l = atoms%bmu(ntyp)
76 : !IF (bmu_l > 0.001 .AND. atoms%numStatesProvided(ntyp) .NE. 0) CALL &
77 : ! judft_warn("You specified both: inital moment and occupation numbers.", &
78 : ! hint="The inital moment will be ignored, set magMom=0.0", calledby="atom2.f90")
79 : !CALL setcor(ntyp, input%jspins, atoms, input, bmu_l, nst, kappa, nprnc, occ)
80 84 : CALL atoms%econf(ntyp)%get_core(nst,nprnc,kappa,occ,l_valence)
81 :
82 84 : if (input%jspins == 2) THEN
83 664 : bmu_l=sign(1.,sum(occ(:nst,1)-occ(:nst,2)))
84 664 : if (any(bmu_l*(occ(:nst,1)-occ(:nst,2))<-1.E-5)) call judft_warn("Inconsistent polarization of starting density")
85 54 : if (bmu_l<0) THEN
86 11 : DO i=1,nst
87 10 : bmu_l=occ(i,1)
88 10 : occ(i,1)=occ(i,2)
89 11 : occ(i,2)=bmu_l
90 : ENDDO
91 : bmu_l=-1
92 : ENDIF
93 : ENDIF
94 : !
95 : !---> for electric field case (sigma.ne.0), add the extra charge
96 : !---> to the uppermost level; ignore the possible problem that
97 : !---> the occupations may not be between 0 and 2
98 84 : IF (input%jspins == 1) THEN
99 30 : occ(nst, 1) = occ(nst, 1) + qdel
100 : ELSE
101 54 : occ(nst, 1) = occ(nst, 1) + qdel/2.
102 54 : occ(nst, input%jspins) = occ(nst, input%jspins) + qdel/2.
103 : ENDIF
104 : !
105 : CALL stpot1(&
106 : & jrc, n, z, rad,&
107 84 : & vr1)
108 72758 : DO i = 1, n
109 72758 : vr1(i, input%jspins) = vr1(i, 1)
110 : ENDDO
111 : !
112 : ! start iterating
113 : !
114 84 : lastit = .false.
115 84 : conv = .true.
116 84 : delrv = 0.100
117 84 : inr0 = log(5.0/rnot)/h + 1.5
118 1977 : DO it = 1, itmax
119 5270 : DO ispin = 1, input%jspins
120 : !
121 : !----> load potential
122 2889249 : DO i = 1, n
123 2889249 : vr(i) = vr1(i, ispin)
124 : ENDDO
125 : !----> adding stabilizing well: m. weinert
126 304605 : DO i = inr0, n
127 304605 : vr(i) = vr(i) + rad(i)*delrv*(rad(i) - rad(inr0))
128 : ENDDO
129 : !----> note that vr contains r*v(r) in hartree units
130 2889249 : DO i = 1, n
131 2889249 : rhoss(i, ispin) = zero
132 : ENDDO
133 3293 : if (lastit) THEN
134 138 : inquire (file="startcharges", exist=l_start)
135 138 : if (l_start) then
136 0 : OPEN (61, file="startcharges")
137 : DO WHILE (.true.)
138 0 : read (61, *, end=888, err=888) i, rho
139 0 : if (i == z) then
140 0 : occ(nst, 1) = occ(nst, 1) + rho
141 0 : goto 888
142 : endif
143 : enddo
144 : 888 continue
145 0 : close (61)
146 : endif
147 : endif
148 41514 : DO 90 k = 1, nst
149 36244 : fn = nprnc(k)
150 36244 : fj = iabs(kappa(k)) - 0.5e0
151 36244 : fl = fj + 0.5e0*isign(1, kappa(k))
152 36244 : e = -2*(z/(fn + fl))**2
153 36244 : ierr = -1
154 36244 : msh_l = jrc
155 72488 : DO WHILE (ierr .NE. 0)
156 : CALL differ(&
157 : & fn, fl, fj, c, z, h, rnot, rn, d, msh_l, vr,&
158 : & e,&
159 36244 : & a, b, ierr)!keep
160 36244 : msh_l = msh_l - 1
161 36244 : IF (jrc - msh_l > 100) CALL juDFT_error(&
162 0 : & "atom2", calledby="atom2")
163 : ENDDO
164 72488 : DO i = msh_l + 1, jrc
165 36244 : a(i) = a(msh_l)
166 72488 : b(i) = b(msh_l)
167 : ENDDO
168 :
169 33010176 : DO i = 1, n
170 33010176 : rh(i) = occ(k, ispin)*(a(i)**2 + b(i)**2)
171 : ENDDO
172 : !+ldau
173 36244 : IF (lastit) THEN ! calculate slater interals
174 1511 : l = int(fl)
175 : ! write(*,*) nprnc(k),l
176 4341 : DO kk = 0, 2*l, 2 ! F0 for s, F0 + F2 for p etc.
177 2830 : r = rnot
178 2595799 : DO i = 1, n
179 2592969 : ain(i) = a(i)**2*r**(-kk - 1) ! prepare inner integrand
180 2595799 : r = r*d
181 : ENDDO
182 : CALL intgr1(ain, rnot, h, n, & ! integrate&
183 2830 : &ahelp)
184 2830 : r = rnot
185 2592969 : DO i = 1, n - 1
186 2590139 : ain(i) = a(i)**2*r**kk*(ahelp(n) - ahelp(i))
187 2592969 : r = r*d
188 : ENDDO
189 : CALL intgr0(ain, rnot, h, n - 1, & ! integrate 2nd r&
190 4341 : &f(kk/2))
191 :
192 : ENDDO
193 : ! write(*,*) (hartree_to_ev_const*2*f(kk),kk=0,l)
194 : ENDIF
195 : !-ldau
196 36244 : eig(k, ispin) = e
197 : !----> calculate <r>
198 33010176 : DO i = 1, n
199 33010176 : a(i) = (a(i)**2 + b(i)**2)*rad(i)
200 : ENDDO
201 36244 : CALL intgr1(a, rnot, h, n, b)
202 36244 : rev(k, ispin) = b(n)
203 33010176 : DO i = 1, n
204 33010176 : rhoss(i, ispin) = rhoss(i, ispin) + rh(i)
205 : ENDDO
206 3293 : 90 ENDDO
207 : ENDDO
208 : !
209 : ! solve poisson's equation
210 : !
211 1715360 : DO i = 1, n
212 1715360 : dens(i) = rhoss(i, 1)
213 : ENDDO
214 1977 : IF (input%jspins == 2) THEN
215 1173889 : DO i = 1, n
216 1173889 : dens(i) = dens(i) + rhoss(i, input%jspins)
217 : ENDDO
218 : ENDIF
219 1977 : CALL intgr1(dens, rnot, h, n, a)
220 1715360 : DO i = 1, n
221 1715360 : rh(i) = dens(i)/rad(i)
222 : ENDDO
223 1977 : CALL intgr1(rh, rnot, h, n, b)
224 1977 : fisr = b(n)
225 1715360 : DO i = 1, n
226 1715360 : vr(i) = (a(i) + rad(i)*(fisr - b(i)) - z)
227 : ENDDO
228 : !+ta
229 5270 : DO ispin = 1, input%jspins
230 2891226 : DO i = 1, n
231 2889249 : rhoss(i, ispin) = rhoss(i, ispin)/(fpi_const*rad(i)**2)
232 : ENDDO
233 : ENDDO
234 1977 : IF (xcpot%needs_grad()) THEN
235 1150 : CALL potl0(xcpot, input%jspins, atoms%dx(ntyp), rad, rhoss, vxc)
236 : ELSE
237 827 : CALL xcpot%get_vxc(input%jspins, rhoss, vxc, vx)
238 : ENDIF
239 5270 : DO ispin = 1, input%jspins
240 2891226 : DO i = 1, n
241 2889249 : vr2(i, ispin) = vr(i) + vxc(i, ispin)*rad(i)
242 : ENDDO
243 : ENDDO
244 : !-ta
245 : ! determine distance of potentials
246 : !
247 1977 : r3 = rn**3
248 1977 : dist = 0.0
249 5270 : DO ispin = 1, input%jspins
250 2889249 : DO i = 1, n
251 2889249 : a(i) = (vr2(i, ispin) - vr1(i, ispin))**2
252 : ENDDO
253 3293 : CALL intgr1(a, rnot, h, n, b)
254 5270 : dist = dist + sqrt((3.0e0/r3)*b(n))
255 : ENDDO
256 1977 : IF (lastit) GO TO 190
257 1893 : IF (dist < distol) lastit = .true.
258 : ! mix new input potential
259 1893 : p1 = 1.0e0/dist
260 1893 : p = min(pmax, p1)
261 1893 : p = max(p, pmin)
262 1893 : WRITE (oUnit, FMT=8060) it, dist, p
263 1893 : p1 = 1.0e0 - p
264 5048 : DO ispin = 1, input%jspins
265 2770185 : DO i = 1, n
266 2768292 : vr1(i, ispin) = p1*vr1(i, ispin) + p*vr2(i, ispin)
267 : ENDDO
268 : ENDDO
269 : ENDDO
270 : !
271 : ! output
272 : !
273 0 : WRITE (oUnit, FMT=8030) dist
274 84 : conv = .false.
275 : ! list eigenvalues
276 84 : 190 IF (conv) WRITE (oUnit, FMT=8040) it, dist
277 222 : DO isp = 1, input%jspins
278 138 : ispin=merge(isp,3-isp,(bmu_l>0).or.(input%jspins<2))
279 138 : WRITE (oUnit, '(a8,i2)') 'spin No.',ispin
280 1649 : DO k = 1, nst
281 1511 : fj = iabs(kappa(k)) - 0.5e0
282 1511 : l = fj + 0.5e0*isign(1, kappa(k)) + 0.01e0
283 1511 : lnum(k) = l
284 1511 : WRITE (oUnit, FMT=8050) nprnc(k), kappa(k), l, fj,&
285 3160 : & occ(k, isp), eig(k, isp), rev(k, isp)
286 : ENDDO
287 : !
288 : !---> guess enpara if it doesn't exist, using floating energy parameters
289 : !
290 138 : i = atoms%jri(ntyp) - (log(4.0)/atoms%dx(ntyp) + 1.51)
291 138 : vbar(isp) = vr1(i, isp)/(rnot*exp(atoms%dx(ntyp)*(i - 1)))
292 222 : WRITE (oUnit, '(/,'' reference energy = '',2f12.6,/)') vbar(isp)
293 : ENDDO
294 :
295 : 8030 FORMAT(/, /, /, ' $$$ error: not converged, dist=', f10.6,/)
296 : 8040 FORMAT(/, /, 3x, 'converged in', i4, ' iterations to a distance of',&
297 : & e12.5, ' har', /, /, 3x, 'n kappa l j ', 5x,&
298 : & 'occ. eigenvalue (har) <r> ',/)
299 : 8050 FORMAT(3x, i1, i5, i5, f6.1, 2(3x, f7.2, 1x, 2f12.6))
300 : 8060 FORMAT('it,dist,p=', i4, 2f12.5)
301 :
302 84 : IF (input%jspins>1.and.bmu_l<0) THEN
303 11 : DO i=1,nst
304 10 : bmu_l=eig(i,1)
305 10 : eig(i,1)=eig(i,2)
306 11 : eig(i,2)=bmu_l
307 : ENDDO
308 966 : DO i=1,size(rhoss,1)
309 965 : bmu_l=rhoss(i,1)
310 965 : rhoss(i,1)=rhoss(i,2)
311 966 : rhoss(i,2)=bmu_l
312 : ENDDO
313 : ENDIF
314 :
315 84 : END SUBROUTINE atom2
316 : END MODULE m_atom2
|
