Line data Source code
1 : MODULE m_hybrid_core
2 :
3 : ! read core radial wavefunctions from corebas
4 : ! corebas is written in cored.F
5 : ! (core basis functions can be read in once during an iteration)
6 :
7 : CONTAINS
8 16 : SUBROUTINE corewf(atoms, jsp, input,&
9 16 : vr, lmaxcd, maxindxc, fmpi, lmaxc, nindxc, core1, core2, eig_c)
10 :
11 : USE m_juDFT
12 : USE m_types
13 : USE m_constants
14 :
15 : IMPLICIT NONE
16 :
17 : TYPE(t_mpi), INTENT(IN) :: fmpi
18 : TYPE(t_input), INTENT(IN) :: input
19 : TYPE(t_atoms), INTENT(IN) :: atoms
20 :
21 : ! - scalars -
22 : INTEGER, INTENT(IN) :: jsp
23 : INTEGER, INTENT(IN) :: lmaxcd
24 : INTEGER, INTENT(INOUT) :: maxindxc
25 :
26 : ! -arrays -
27 : INTEGER, INTENT(INOUT) :: lmaxc(:)
28 : INTEGER, INTENT(INOUT) :: nindxc(0:lmaxcd, atoms%ntype)
29 :
30 : REAL, INTENT(IN) :: vr(:, :, :)!(atoms%jmtd,atoms%ntypd,input%jspins)
31 : REAL, INTENT(INOUT) :: core1(:, :, 0:, :) !(atoms%jmtd,maxindxc,0:lmaxcd,atoms%ntype)
32 : REAL, INTENT(INOUT) :: core2(:, :, 0:, :) !(jmtd,maxindxc,0:lmaxcd,ntype)
33 :
34 : REAL, INTENT(INOUT) :: eig_c(maxindxc, 0:lmaxcd, atoms%ntype)
35 :
36 : ! - local scalars -
37 : INTEGER :: ncstd
38 : INTEGER :: ok, itype, i, l
39 :
40 : REAL :: weight1, weight2
41 : ! - local arrays -
42 16 : INTEGER, ALLOCATABLE :: nindxcr(:, :)
43 :
44 16 : REAL, ALLOCATABLE :: core1r(:, :, :, :), core2r(:, :, :, :)
45 16 : REAL, ALLOCATABLE :: eig_cr(:, :, :)
46 :
47 16 : call timestart("corewf")
48 40 : ncstd = maxval(atoms%econf%num_core_states)
49 64 : allocate(nindxcr(0:ncstd, atoms%ntype), stat=ok)
50 :
51 : ! generate relativistic core wave functions( ->core1r,core2r )
52 : CALL calcorewf( input, jsp, atoms,&
53 : ncstd, vr,&
54 16 : lmaxc, nindxcr, core1r, core2r, eig_cr, fmpi)
55 :
56 88 : nindxc = 0
57 :
58 : ! average over core states that only differ in j
59 : ! core functions with l-qn equal 0 doesnot change during the average process
60 :
61 40 : nindxc(0, :) = nindxcr(0, :)
62 40 : DO itype = 1, atoms%ntype
63 : core1(:, :nindxc(0, itype), 0, itype)&
64 44616 : = core1r(:, 0, :nindxc(0, itype), itype)
65 : core2(:, :nindxc(0, itype), 0, itype)&
66 44616 : = core2r(:, 0, :nindxc(0, itype), itype)
67 : eig_c(:nindxc(0, itype), 0, itype)&
68 96 : = eig_cr(0, :nindxc(0, itype), itype)
69 : END DO
70 :
71 40 : DO itype = 1, atoms%ntype
72 64 : DO l = 1, lmaxc(itype)
73 24 : weight1 = 2*(l - 0.5) + 1
74 24 : weight2 = 2*(l + 0.5) + 1
75 48 : IF (modulo(nindxcr(l, itype), 2) == 0) THEN
76 24 : DO i = 1, nindxcr(l, itype), 2
77 32 : nindxc(l, itype) = nindxc(l, itype) + 1
78 : core1(:atoms%jri(itype), nindxc(l, itype), l, itype) =&
79 : (weight1*core1r(:atoms%jri(itype), l, i, itype) +&
80 : weight2*core1r(:atoms%jri(itype), l, i + 1, itype))&
81 24928 : /(weight1 + weight2)
82 : core2(:atoms%jri(itype), nindxc(l, itype), l, itype) =&
83 : (weight1*core2r(:atoms%jri(itype), l, i, itype) + &
84 : weight2*core2r(:atoms%jri(itype), l, i + 1, itype))&
85 24928 : /(weight1 + weight2)
86 :
87 : eig_c(nindxc(l, itype), l, itype) =&
88 : (weight1*eig_cr(l, i, itype) +&
89 : weight2*eig_cr(l, i + 1, itype))&
90 32 : /(weight1 + weight2)
91 : END DO
92 : ELSE
93 0 : DO i = 1, nindxcr(l, itype) - 1, 2
94 0 : nindxc(l, itype) = nindxc(l, itype) + 1
95 : core1(:atoms%jri(itype), nindxc(l, itype), l, itype) =&
96 : (weight1*core1r(:atoms%jri(itype), l, i, itype) +&
97 : weight2*core1r(:atoms%jri(itype), l, i + 1, itype))&
98 0 : /(weight1 + weight2)
99 : core2(:atoms%jri(itype), nindxc(l, itype), l, itype) =&
100 : (weight1*core2r(:atoms%jri(itype), l, i, itype) + &
101 : weight2*core2r(:atoms%jri(itype), l, i + 1, itype))&
102 0 : /(weight1 + weight2)
103 :
104 : eig_c(nindxc(l, itype), l, itype) =&
105 : (weight1*eig_cr(l, i, itype) +&
106 : weight2*eig_cr(l, i + 1, itype))&
107 0 : /(weight1 + weight2)
108 : END DO
109 0 : nindxc(l, itype) = nindxc(l, itype) + 1
110 : core1(:atoms%jri(itype), nindxc(l, itype), l, itype)&
111 0 : = core1r(:atoms%jri(itype), l, nindxcr(l, itype), itype)
112 : core2(:atoms%jri(itype), nindxc(l, itype), l, itype)&
113 0 : = core2r(:atoms%jri(itype), l, nindxcr(l, itype), itype)
114 : eig_c(nindxc(l, itype), l, itype)&
115 0 : = eig_cr(l, nindxcr(l, itype), itype)
116 : END IF
117 :
118 : END DO
119 : END DO
120 :
121 16 : deallocate(nindxcr, core1r, core2r, eig_cr)
122 :
123 88 : IF (maxindxc /= maxval(nindxc))&
124 0 : call judft_error('corewf: counting error nindxc')
125 16 : call timestop("corewf")
126 16 : END SUBROUTINE corewf
127 :
128 16 : SUBROUTINE calcorewf( input, jspin, atoms,&
129 16 : ncstd, vr,&
130 16 : lmaxc, nindxcr, core1, core2, eig_c, fmpi)
131 :
132 : USE m_intgr, ONLY: intgr3, intgr0, intgr1
133 : USE m_constants
134 : USE m_differ
135 : USE m_types
136 : IMPLICIT NONE
137 :
138 : TYPE(t_mpi), INTENT(IN) :: fmpi
139 : TYPE(t_input), INTENT(IN) :: input
140 : TYPE(t_atoms), INTENT(IN) :: atoms
141 :
142 : ! - scalars -
143 : INTEGER, INTENT(IN) :: ncstd
144 : INTEGER, INTENT(IN) :: jspin
145 : INTEGER, INTENT(INOUT):: lmaxc(:)
146 :
147 : ! - arrays -
148 : INTEGER, INTENT(INOUT):: nindxcr(0:ncstd, atoms%ntype)
149 : REAL, INTENT(IN) :: vr(:, :, :)!(atoms%jmtd,atoms%ntypd,input%jspins)
150 : REAL, ALLOCATABLE :: core1(:, :, :, :), core2(:, :, :, :)
151 : REAL, ALLOCATABLE :: eig_c(:, :, :)
152 :
153 : ! - local scalars -
154 : INTEGER :: i, j, itype, korb, ncmsh, nst, ierr
155 : REAL :: e, fj, fl, fn, t2, c, bmu, weight
156 : REAL :: d, dxx, rn, rnot, z, t1, rr
157 : LOGICAL, SAVE :: first = .true.
158 :
159 : ! - local arrays -
160 : INTEGER :: kappa(29), nprnc(29)
161 16 : REAL :: vrd(atoms%msh)
162 16 : REAL :: occ(29), occ_h(29, 2), a(atoms%msh), b(atoms%msh)
163 : REAL, ALLOCATABLE, SAVE:: vr0(:, :, :)
164 :
165 : ! - intrinsic functions -
166 : INTRINSIC exp, iabs, isign
167 16 : call timestart("calcorewf")
168 :
169 16 : c = c_light(1.0)
170 :
171 16 : IF (first) THEN
172 30 : allocate(vr0(atoms%jmtd, atoms%ntype, input%jspins))
173 : END IF
174 :
175 16 : IF (input%frcor) THEN
176 0 : IF (first) THEN
177 0 : vr0 = vr
178 0 : first = .false.
179 : END IF
180 : ELSE
181 25368 : vr0 = vr
182 16 : first = .false.
183 : END IF
184 :
185 : ! this loop determines the dimensions
186 :
187 208 : lmaxc = 0; nindxcr = 0
188 40 : DO itype = 1, atoms%ntype
189 24 : z = atoms%zatom(itype)
190 24 : dxx = atoms%dx(itype)
191 24 : bmu = 0.0
192 24 : call atoms%econf(itype)%get_core(nst,nprnc,kappa,occ_h)
193 :
194 : IF ((bmu > 99.)) THEN
195 : occ(1:nst) = input%jspins*occ_h(1:nst, jspin)
196 : ELSE
197 144 : occ(1:nst) = occ_h(1:nst, 1)
198 : END IF
199 24 : rnot = atoms%rmsh(1, itype)
200 24 : d = exp(atoms%dx(itype))
201 24 : ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1)
202 24 : ncmsh = min(ncmsh, atoms%msh)
203 24 : rn = rnot*(d**(ncmsh - 1))
204 :
205 24 : nst = atoms%econf(itype)%num_core_states
206 :
207 160 : DO korb = 1, nst
208 144 : IF (occ(korb) > 0) THEN
209 120 : fn = nprnc(korb)
210 120 : fj = iabs(kappa(korb)) - 0.5
211 120 : weight = 2*fj + 1.
212 : IF (bmu > 99.) weight = occ(korb)
213 120 : fl = fj + (0.5)*isign(1, kappa(korb))
214 120 : e = -2*(z/(fn + fl))**2
215 :
216 120 : nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1
217 120 : lmaxc(itype) = max(lmaxc(itype), NINT(fl))
218 : ENDIF
219 : END DO
220 :
221 : END DO
222 :
223 288 : allocate(core1(atoms%jmtd, 0:maxval(lmaxc), maxval(nindxcr), atoms%ntype))
224 288 : allocate(core2(atoms%jmtd, 0:maxval(lmaxc), maxval(nindxcr), atoms%ntype))
225 272 : allocate(eig_c(0:maxval(lmaxc), maxval(nindxcr), atoms%ntype))
226 :
227 202304 : core1 = 0; core2 = 0
228 184 : nindxcr = 0
229 40 : DO itype = 1, atoms%ntype
230 24 : z = atoms%zatom(itype)
231 24 : dxx = atoms%dx(itype)
232 24 : bmu = 0.0
233 24 : call atoms%econf(itype)%get_core(nst,nprnc,kappa,occ_h)
234 : IF ((bmu > 99.)) THEN
235 : occ(1:nst) = input%jspins*occ_h(1:nst, jspin)
236 : ELSE
237 144 : occ(1:nst) = occ_h(1:nst, 1)
238 : END IF
239 24 : rnot = atoms%rmsh(1, itype)
240 24 : d = exp(atoms%dx(itype))
241 24 : ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1)
242 24 : ncmsh = min(ncmsh, atoms%msh)
243 24 : rn = rnot*(d**(ncmsh - 1))
244 24 : IF (fmpi%irank == 0) THEN
245 12 : WRITE (oUnit, FMT=8000) z, rnot, dxx, atoms%jri(itype)
246 : END IF
247 18992 : DO j = 1, atoms%jri(itype)
248 18992 : vrd(j) = vr0(j, itype, jspin)
249 : END DO
250 :
251 24 : if (input%l_core_confpot) THEN
252 :
253 : ! linear extension of the potential with slope t1 / a.u.
254 24 : t1 = 0.125
255 24 : t2 = vrd(atoms%jri(itype))/atoms%rmt(itype) - atoms%rmt(itype)*t1
256 24 : rr = atoms%rmt(itype)
257 : else
258 0 : t2 = vrd(atoms%jri(itype))/(atoms%jri(itype) - ncmsh)
259 : endif
260 24 : IF (atoms%jri(itype) < ncmsh) THEN
261 2596 : DO i = atoms%jri(itype) + 1, ncmsh
262 2596 : if (input%l_core_confpot) THEN
263 2572 : rr = d*rr
264 2572 : vrd(i) = rr*(t2 + rr*t1)
265 : else
266 0 : vrd(i) = vrd(atoms%jri(itype)) + t2*(i - atoms%jri(itype))
267 : endif
268 : !
269 : END DO
270 : END IF
271 :
272 24 : nst = atoms%econf(itype)%num_core_states
273 :
274 160 : DO korb = 1, nst
275 144 : IF (occ(korb) > 0) THEN
276 120 : fn = nprnc(korb)
277 120 : fj = iabs(kappa(korb)) - 0.5
278 120 : weight = 2*fj + 1.
279 : IF (bmu > 99.) weight = occ(korb)
280 120 : fl = fj + (0.5)*isign(1, kappa(korb))
281 120 : e = -2*(z/(fn + fl))**2
282 120 : CALL differ(fn, fl, fj, c, z, dxx, rnot, rn, d, ncmsh, vrd, e, a, b, ierr)
283 :
284 120 : nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1
285 :
286 : core1(:atoms%jri(itype), NINT(fl), nindxcr(NINT(fl), itype), itype)&
287 93776 : = a(:atoms%jri(itype))
288 : core2(:atoms%jri(itype), NINT(fl), nindxcr(NINT(fl), itype), itype)&
289 93776 : = b(:atoms%jri(itype))
290 :
291 120 : eig_c(NINT(fl), nindxcr(NINT(fl), itype), itype) = e
292 :
293 120 : IF (fmpi%irank == 0) THEN
294 60 : WRITE (oUnit, FMT=8010) fn, fl, fj, e, weight
295 : END IF
296 240 : IF (ierr /= 0) call judft_error('error in core-level routine')
297 : ENDIF
298 : END DO
299 :
300 : END DO
301 :
302 : 8000 FORMAT(/, /, 10x, 'z=', f4.0, 5x, 'r(1)=', e14.6, 5x, 'dx=', f8.6, 5x,&
303 : 'm.t.index=', i4, /, 15x, 'n', 4x, 'l', 5x, 'j', 4x, 'energy', 7x,&
304 : 'weight')
305 : 8010 FORMAT(12x, 2f5.0, f6.1, f10.4, f10.0)
306 :
307 16 : call timestop("calcorewf")
308 16 : END SUBROUTINE calcorewf
309 :
310 12 : SUBROUTINE core_init(input, atoms, lmaxcd, maxindxc)
311 :
312 : USE m_types
313 : USE m_constants
314 : USE m_intgr, ONLY: intgr3, intgr0, intgr1
315 : USE m_differ
316 : IMPLICIT NONE
317 :
318 : TYPE(t_input), INTENT(IN) :: input
319 : TYPE(t_atoms), INTENT(IN) :: atoms
320 : INTEGER, INTENT(INOUT) :: maxindxc, lmaxcd
321 :
322 : ! - local scalars -
323 : INTEGER :: itype, korb, ncmsh, nst
324 : REAL :: e, fj, fl, fn, bmu
325 : REAL :: d, dxx, rn, rnot, z
326 :
327 : ! - local arrays -
328 : INTEGER :: kappa(29), nprnc(29)
329 12 : INTEGER :: nindxcr(0:29, atoms%ntype)
330 12 : REAL :: occ(29), occ_h(29, input%jspins)
331 12 : INTEGER :: lmaxc(atoms%ntype)
332 :
333 : ! - intrinsic functions -
334 : INTRINSIC exp, iabs, isign
335 :
336 : ! this loop determines the dimensions
337 :
338 652 : lmaxc = 0; nindxcr = 0
339 32 : DO itype = 1, atoms%ntype
340 20 : z = atoms%zatom(itype)
341 20 : dxx = atoms%dx(itype)
342 20 : bmu = 0.0
343 :
344 20 : call atoms%econf(itype)%get_core(nst,nprnc,kappa,occ_h)
345 124 : occ(1:nst) = occ_h(1:nst, 1)
346 :
347 20 : rnot = atoms%rmsh(1, itype)
348 : d = exp(atoms%dx(itype))
349 : ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1)
350 20 : ncmsh = min(ncmsh, atoms%msh)
351 : rn = rnot*(d**(ncmsh - 1))
352 :
353 20 : nst = atoms%econf(itype)%num_core_states
354 :
355 136 : DO korb = 1, nst
356 124 : IF (occ(korb) > 0) THEN
357 104 : fn = nprnc(korb)
358 104 : fj = iabs(kappa(korb)) - 0.5
359 :
360 104 : fl = fj + (0.5)*isign(1, kappa(korb))
361 104 : e = -2*(z/(fn + fl))**2
362 :
363 104 : nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1
364 104 : lmaxc(itype) = max(lmaxc(itype), NINT(fl))
365 : ENDIF
366 : END DO
367 :
368 : END DO
369 :
370 32 : lmaxcd = maxval(lmaxc)
371 : ! The following commented line seems to be wrong but doesn't produce any problems.
372 : ! I replace it by the two following lines. If there is a system for which the
373 : ! error in corewf is thrown this has to be undone.
374 : ! maxindxc = maxval(nindxcr(0, :), nint((maxval(nindxcr)/2.0)))
375 32 : maxindxc = maxval(nindxcr(0, :))
376 632 : maxindxc = MAX(maxindxc, nint((maxval(nindxcr)/2.0)))
377 :
378 12 : END SUBROUTINE core_init
379 :
380 : END MODULE m_hybrid_core
|
