Line data Source code
1 : MODULE m_core
2 :
3 : CONTAINS
4 :
5 2 : SUBROUTINE core(mrad,vt,bt,zz,stval,dx,nlshell,nqntab,lqntab,jtop,ectab,rhochr,rhospn)
6 :
7 : USE m_constants
8 : USE m_felim
9 : USE m_findlim
10 : USE m_cnodes
11 : USE m_coredir
12 : USE m_ccsdnt
13 : USE m_rsimp
14 : USE m_cfnorm
15 : USE m_coreerr
16 : USE m_corehff
17 : USE m_nwrfst
18 : USE m_rinvgj
19 :
20 : IMPLICIT NONE
21 :
22 : INTEGER, INTENT (IN) :: mrad,jtop,nlshell
23 : REAL, INTENT (IN) :: dx,stval,zz
24 : INTEGER, INTENT (IN) :: lqntab(15),nqntab(15)
25 : REAL, INTENT (IN) :: bt(mrad),vt(mrad)
26 : REAL, INTENT (OUT) :: rhochr(mrad),rhospn(mrad)
27 : REAL, INTENT (INOUT) :: ectab(100)
28 :
29 : REAL :: bhf,bsum,btot,dvstep,ec,elim,qcor,spcor,tolvar,vz,xmj,ec_sv
30 : INTEGER :: i,ic,ic1,ic2,ie,iflag,ii,ilshell,imin,ir,ish,istart,iter
31 : INTEGER :: itermax,iv,j,jv,kap1,kap2,l,lll,muem05,n,ncor,nmatch,node
32 : INTEGER :: nqn,nrmax,nsh,nsol,nval,nvar,nzero,s,t
33 : LOGICAL :: ferro
34 :
35 2 : REAL :: bb(mrad),bhff(15),dedv(4,4),dv(4),dvde(4,4),err(4),errnew(4)
36 2 : REAL :: fc(2,2,mrad),fck(2,2,mrad),gc(2,2,mrad),gck(2,2,mrad)
37 2 : REAL :: piw(2,2),pow(2,2),qiw(2,2),qow(2,2),rc(mrad),rc2(mrad)
38 2 : REAL :: rchr(mrad),rspn(mrad),var(4),varnew(4),vv(mrad)
39 : INTEGER kap(2)
40 : CHARACTER txtl(0:3)
41 : CHARACTER*3 txtk(4)
42 :
43 : DATA txtl/'s','p','d','f'/
44 : DATA txtk/'1/2','3/2','5/2','7/2'/
45 : ! DATA nqntab/1,2,2,3,3,3,4,4,4,4,5,5,5,6,6/
46 : ! DATA lqntab/0,0,1,0,1,2,0,1,2,3,0,1,2,0,1/
47 2 : tolvar = 1.0e-12
48 2 : dvstep = 0.010
49 2 : itermax = 150
50 :
51 2 : nrmax = mrad
52 1348 : DO n = 1,nrmax
53 1346 : rc(n) = exp(stval+ (n-1)*dx)
54 1348 : rc2(n) = rc(n)**2
55 : END DO
56 :
57 1348 : DO n = 1,nrmax
58 1346 : rhochr(n) = 0.00
59 1348 : rhospn(n) = 0.00
60 : END DO
61 :
62 2 : bsum = 0.0
63 1348 : DO n = 1,jtop
64 1346 : vv(n) = vt(n)
65 1346 : bb(n) = bt(n)
66 1348 : bsum = bsum + abs(bb(n))
67 : END DO
68 :
69 2 : ncor = 0
70 12 : DO ilshell = 1,nlshell
71 12 : ncor = ncor + 2* (2*lqntab(ilshell)+1)
72 : END DO
73 2 : nval = int(zz) - ncor
74 :
75 2 : IF (bsum.GT.1.0E-8) THEN
76 2 : ferro = .true.
77 : ELSE
78 0 : ferro = .false.
79 0 : WRITE (oUnit,FMT='('' PARAMAGNETIC CASE'',/,'' *****************'')')
80 : END IF
81 :
82 2 : WRITE (oUnit,FMT='(/,'' ATOMIC NUMBER : '',F6.2,/,'' CORE ELECTRONS : '',I5,/,'' VAL. ELECTRONS : '',I5)') zz,ncor,nval
83 2 : WRITE (oUnit,FMT='( /,'' MESH RC(1) : '',F12.6,//10x,''MUE KAP ITER ENERGY '',''(Hart)'', " HF (KG) ")') rc(1)
84 :
85 :
86 : ! ---------------------------------------
87 : ! INITIALIZE QUANTUM NUMBERS NQN AND L
88 : ! ---------------------------------------
89 :
90 2 : btot = 0.0
91 2 : ic = 0
92 :
93 12 : DO ilshell = 1,nlshell ! nl - loop
94 10 : nqn = nqntab(ilshell)
95 10 : l = lqntab(ilshell)
96 10 : nsh = 2* (2*l+1)
97 10 : bhff(ilshell) = 0.0
98 10 : ish = 0
99 :
100 38 : DO muem05 = -l - 1,+l ! \mu - loop
101 28 : xmj = muem05 + 0.5
102 28 : kap1 = -l - 1
103 28 : kap2 = l
104 28 : kap(1) = kap1
105 28 : kap(2) = kap2
106 :
107 28 : lll = l* (l+1)
108 28 : IF (abs(xmj).GT.l) THEN
109 20 : nsol = 1
110 : ELSE
111 8 : nsol = 2
112 : END IF
113 :
114 28 : nvar = 2*nsol
115 :
116 74 : DO s = 1,nsol ! s - loop : solutions for each (nl,\mu)
117 108 : DO i = 1,2
118 252 : DO j = 1,2
119 97128 : DO n = 1,mrad
120 96912 : gc(i,j,n) = 0.0
121 96912 : gck(i,j,n) = 0.0
122 96912 : fc(i,j,n) = 0.0
123 97056 : fck(i,j,n) = 0.0
124 : END DO
125 : END DO
126 : END DO
127 :
128 36 : ic = ic + 1
129 36 : ish = ish + 1
130 36 : t = 3 - s
131 36 : ec = ectab(ic)
132 36 : write(*,*) ic,ec
133 :
134 36 : IF (ic > 1) THEN
135 34 : IF (ectab(ic-1) > 0.1) THEN
136 0 : vv(:) = vv(:) + 2*ec_sv
137 0 : bb(:) = bb(:) + 2*ec_sv
138 : ENDIF
139 : ENDIF
140 36 : IF (ec.GT.0.1) THEN
141 0 : vv(:) = vv(:) - 2*ectab(ic)
142 0 : bb(:) = bb(:) - 2*ectab(ic)
143 0 : ec_sv = ectab(ic)
144 0 : ec = -ectab(ic)
145 : ELSE
146 : ec_sv = 0.0
147 : ENDIF
148 36 : istart = 1
149 :
150 36 : IF (ish.GT.1) GO TO 30
151 :
152 : ! energy limits : elim and emax
153 10 : CALL felim(mrad,lll,zz,nqn,vv,rc,elim)
154 10 : write(*,*) 'elim',elim
155 :
156 10 : 20 IF (ec.LE.elim) ec = elim*0.70
157 :
158 : ! turning point and physical infinity : nmatch and nzero
159 10 : CALL findlim(mrad,lll,ec,vv,rc,nmatch,nzero)
160 36 : write(*,*) 'nmatch,nzero',nmatch,nzero
161 : 30 CONTINUE
162 :
163 : ! nodes for large component : if IFLAG=0 number of nodes is OK
164 36 : iflag = 0
165 36 : CALL cnodes(mrad,iflag,s,ec,l,xmj,nqn,vv,bb,rc,dx,nmatch,nzero,gc,fc,pow,qow,piw,qiw,node)
166 36 : IF (iflag.NE.0) GO TO 20
167 :
168 : ! setup of Newton-Raphson algorithm
169 36 : CALL nwrfst(mrad,nsol,s,t,nmatch,nzero,ferro,ec,rc,pow,piw,gc,err,var,dv,varnew,errnew)
170 :
171 36 : CALL coreerr(err,var,s,nsol,pow,qow,piw,qiw)
172 140 : DO iv = 1,nvar
173 140 : dv(iv) = var(iv)
174 : END DO
175 :
176 36 : iter = 0
177 : ! iterations start
178 148 : 50 iter = iter + 1
179 :
180 : !----------------------------------
181 : ! CHECK WHETHER NUMBER OF NODES O.K.
182 : !----------------------------------
183 148 : IF (iter.GT.1) THEN
184 112 : node = 0
185 56045 : DO n = 2, (nmatch-1)
186 56045 : IF (gc(s,s,n)*gc(s,s,n-1).LT.0.0) node = node + 1
187 : END DO
188 :
189 112 : IF (node.NE. (nqn-l-1)) THEN
190 0 : IF (node.GT. (nqn-l-1)) THEN
191 0 : ec = 1.2*ec
192 0 : write (*,*) 'up',node,ec
193 : ELSE
194 0 : ec = 0.8*ec
195 0 : write (*,*) 'do',node,ec
196 : END IF
197 0 : istart = istart + 1
198 0 : IF (istart.LT.20) GO TO 20
199 : END IF
200 : END IF
201 :
202 : ! - atomic energy search -
203 450 : DO iv = 2,nvar
204 1368 : DO jv = 1,nvar
205 1368 : varnew(jv) = var(jv)
206 : END DO
207 302 : varnew(iv) = var(iv) + dv(iv)*dvstep
208 :
209 302 : IF (abs(var(iv)).EQ.0.00) THEN
210 0 : IF (ferro) THEN
211 0 : WRITE (oUnit,FMT='(A,I3,A)') ' VAR ',iv,' = 0 ??????!!!!!'
212 : END IF
213 : ELSE
214 302 : IF (abs(dv(iv)/var(iv)).LT.tolvar) varnew(iv) = var(iv)*(1.00+sign(dvstep*tolvar,dv(iv)))
215 : END IF
216 :
217 302 : CALL coreerr(errnew,varnew,s,nsol,pow,qow,piw,qiw)
218 :
219 1516 : DO ie = 1,nvar
220 1368 : IF ((errnew(ie)-err(ie)).EQ.0.00) THEN
221 0 : dedv(ie,iv) = 0.00
222 0 : IF ((ie.EQ.iv) .AND. .NOT.ferro) THEN
223 0 : dedv(ie,iv) = 1.00
224 : END IF
225 : ELSE
226 1066 : dedv(ie,iv) = (errnew(ie)-err(ie)) / (varnew(iv)-var(iv))
227 : END IF
228 : END DO
229 : END DO
230 :
231 598 : DO jv = 1,nvar
232 598 : varnew(jv) = var(jv)
233 : END DO
234 148 : varnew(1) = var(1) + dv(1)*dvstep
235 148 : IF (abs(dv(1)/var(1)).LT.tolvar) varnew(1) = var(1) * (1.00+sign(dvstep*tolvar,dv(1)))
236 :
237 148 : CALL coredir(mrad,varnew(1),l,xmj,1,vv,bb,rc,dx,nmatch,nzero,gc,fc,pow,qow,piw,qiw)
238 148 : CALL coredir(mrad,varnew(1),l,xmj,2,vv,bb,rc,dx,nmatch,nzero,gc,fc,pow,qow,piw,qiw)
239 :
240 148 : CALL coreerr(errnew,varnew,s,nsol,pow,qow,piw,qiw)
241 :
242 598 : DO ie = 1,nvar
243 598 : dedv(ie,1) = (errnew(ie)-err(ie))/ (varnew(1)-var(1))
244 : END DO
245 :
246 148 : CALL rinvgj(dvde,dedv,4,nvar)
247 :
248 598 : DO iv = 1,nvar
249 450 : dv(iv) = 0.00
250 1966 : DO ie = 1,nvar
251 1966 : dv(iv) = dv(iv) + dvde(iv,ie)*err(ie)
252 : END DO
253 598 : var(iv) = var(iv) - dv(iv)
254 : END DO
255 :
256 148 : IF (var(1).GT.0.00) THEN
257 0 : WRITE (oUnit,FMT=*) ' WARNING FROM <CORE> E=',var(1)
258 0 : var(1) = -0.20
259 : END IF
260 :
261 148 : CALL coredir(mrad,var(1),l,xmj,1,vv,bb,rc,dx,nmatch,nzero,gc,fc,pow,qow,piw,qiw)
262 148 : CALL coredir(mrad,var(1),l,xmj,2,vv,bb,rc,dx,nmatch,nzero,gc,fc,pow,qow,piw,qiw)
263 :
264 148 : CALL coreerr(err,var,s,nsol,pow,qow,piw,qiw)
265 :
266 148 : ec = var(1)
267 :
268 : ! check relative change in parameters
269 : ! parameters 3 and 4 = 0 for paramagnetic case !
270 148 : IF (iter.LT.itermax) THEN
271 148 : IF ((nsol.EQ.2) .AND. .NOT.ferro) THEN
272 0 : DO iv = 3,4
273 0 : err(iv) = 0.00
274 0 : errnew(iv) = 0.00
275 0 : var(iv) = 0.00
276 0 : varnew(iv) = 0.00
277 0 : dv(iv) = 0.00
278 : END DO
279 : END IF
280 252 : DO iv = 1,nvar
281 252 : IF (abs(var(iv)).EQ.0.00) THEN
282 0 : IF (ferro) THEN
283 0 : WRITE (oUnit,FMT='(A,I3,A)') ' VAR ',iv,' = 0 ??????!!!!!'
284 : END IF
285 : ELSE
286 216 : IF (abs(dv(iv)/var(iv)).GT.tolvar) GO TO 50
287 : END IF
288 : END DO
289 : ELSE
290 : WRITE (oUnit, FMT='('' ITERATION NOT CONVERGED AFTER'',I3,'' STEPS !'',/,'' PARAMETERS:'',4E18.10,/,'' LAST CORR.:'',4E18.10,/,'' LAST ERROR:'',4E18.10)') &
291 0 : itermax, (var(iv),iv=1,4),(dv(iv),iv=1,4), (err(ie),ie=1,4)
292 : END IF
293 :
294 : ! end of iterations
295 36 : CALL cfnorm(mrad,s,t,nsol,nmatch,jtop,var,gc,fc,rc,rc2,dx,gck,fck)
296 :
297 : !--------------------------
298 : ! CALCULATE CHARGE DENSITY
299 : !--------------------------
300 36 : CALL ccsdnt(mrad,s,jtop,nsol,l,xmj,kap1,kap2,gck,fck,rc2,rchr,rspn)
301 24264 : DO ir = 1,jtop
302 24228 : rhochr(ir) = rhochr(ir) + rchr(ir)
303 24264 : rhospn(ir) = rhospn(ir) + rspn(ir)
304 : END DO
305 :
306 : !-----------------
307 : ! HYPERFINE FIELD
308 : !-----------------
309 36 : CALL corehff(mrad,kap1,kap2,xmj,s,nsol,bhf,gck,fck,rc,dx,jtop)
310 36 : bhff(ilshell) = bhff(ilshell) + bhf
311 :
312 : ! final info
313 36 : ectab(ic) = ec + 2*ec_sv
314 :
315 84 : IF (ish.LT.nsh) THEN
316 26 : WRITE (oUnit,FMT=8000) nqn,txtl(l),txtk(iabs(kap(s))),(2*muem05+1),kap(s),iter,ec/2.
317 : ELSE
318 10 : WRITE (oUnit,FMT=8000) nqn,txtl(l),txtk(iabs(kap(s))),(2*muem05+1),kap(s),iter,ec/2.
319 :
320 : !----------------------------
321 : ! CHECK CONSISTENCY OF RESULTS
322 : !----------------------------
323 10 : IF (l.NE.0) THEN
324 4 : ic1 = ic - nsh + 1
325 4 : ic2 = ic
326 4 : IF (ectab(ic2).GE.ectab(ic1)) THEN
327 : imin = ic1
328 : vz = +1.00
329 : ELSE
330 0 : imin = ic2
331 0 : vz = -1.00
332 : END IF
333 4 : iflag = 0
334 4 : ii = 1
335 4 : DO i = ic1 + 1,ic2,2
336 12 : IF (vz* (ectab(i)-ectab(i-ii)).LT.0.0) iflag = 1
337 12 : ii = 2
338 : END DO
339 4 : IF (ectab(ic1+2).GT.ectab(imin)) iflag = 1
340 8 : DO i = ic1 + 4,ic2 - 1,2
341 4 : IF (ectab(i).GT.ectab(imin)) iflag = 1
342 8 : IF (vz* (ectab(i)-ectab(i-ii)).GT.0.0) iflag = 1
343 : END DO
344 :
345 4 : IF (ferro .AND. (iflag.EQ.1)) THEN
346 0 : WRITE (oUnit,FMT='('' >>> CHECK DATA '', '' E(KAP,MJ) SHOULD BE MONOTONOUS AND '', '' E(+L,MJ) < E(-L-1,MJ) '',//)')
347 : END IF
348 : END IF
349 : END IF
350 :
351 : END DO ! s - loop end
352 :
353 : END DO ! \mu - loop end
354 :
355 10 : WRITE (oUnit,FMT='(I4,A1,20X,F16.3)') nqn,txtl(l),bhff(ilshell)
356 12 : btot = bhff(ilshell) + btot
357 :
358 : END DO ! nl - loop end
359 :
360 2 : WRITE (oUnit,FMT='(''HFTOT:'',20X,F16.3)') btot
361 :
362 2 : qcor = rsimp(mrad,rhochr,rc,jtop,dx)
363 2 : WRITE (oUnit,FMT=8020) 'charge',qcor
364 2 : spcor = rsimp(mrad,rhospn,rc,jtop,dx)
365 2 : WRITE (oUnit,FMT=8020) 'spin',spcor
366 :
367 : ! that's all
368 :
369 :
370 : 8000 FORMAT (I4,A1,A3,I3,'/2',2I4,1X,F14.7,F16.3,:,F16.3,/)
371 : 8010 FORMAT (/,' NQN=',I2,' L=',I2,' KAP=',I2,' MJ=',I2,'/2 IC=',I3,' ISH=',I2,/,' E(',A5,')=',F15.5,/,' NMATCH =',I5,' R=',F10.5,/,' NZERO =',I5,' R=',F10.5,/,' NODES =',I5,' RAT=',E11.4)
372 : 8020 FORMAT (' integrated ',A,':',F12.8)
373 :
374 2 : END SUBROUTINE core
375 : END MODULE m_core
|
