Line data Source code
1 : MODULE m_stden
2 : USE m_juDFT
3 : REAL,PARAMETER :: input_ellow=-2.0
4 : REAL,PARAMETER :: input_elup=1.0
5 : ! ************************************************************
6 : ! generate flapw starting density by superposition of
7 : ! atomic densities. the non-spherical terms inside
8 : ! the spheres are obtained by a least squares fit
9 : ! and the interstitial and vacuum warping terms are calculated
10 : ! by a fast fourier transform.
11 : ! e. wimmer nov. 1984 c.l.fu diagonized 1987
12 : ! *************************************************************
13 :
14 : CONTAINS
15 :
16 :
17 136 : SUBROUTINE stden(fmpi,sphhar,stars,atoms,sym,vacuum,&
18 : input,cell,xcpot,noco )
19 :
20 : USE m_juDFT_init
21 : USE m_types
22 : USE m_types_xcpot_inbuild
23 : USE m_constants
24 : USE m_qsf
25 : USE m_checkdopall
26 : USE m_cdnovlp
27 : USE m_cdn_io
28 : USE m_qfix
29 : USE m_atom2
30 : USE m_clebsch
31 : USE m_rotMMPmat
32 : USE m_RelaxSpinAxisMagn
33 : IMPLICIT NONE
34 :
35 : TYPE(t_mpi),INTENT(IN) :: fmpi
36 : TYPE(t_atoms),INTENT(IN) :: atoms
37 :
38 : TYPE(t_sphhar),INTENT(IN) :: sphhar
39 : TYPE(t_sym),INTENT(IN) :: sym
40 : TYPE(t_stars),INTENT(IN) :: stars
41 : TYPE(t_noco),INTENT(IN) :: noco
42 :
43 : TYPE(t_input),INTENT(IN) :: input
44 : TYPE(t_vacuum),INTENT(IN) :: vacuum
45 : TYPE(t_cell),INTENT(IN) :: cell
46 : CLASS(t_xcpot),INTENT(IN) :: xcpot
47 :
48 : ! Local type instances
49 136 : TYPE(t_potden) :: den
50 136 : TYPE(t_enpara) :: enpara
51 : ! Local Scalars
52 : REAL d,del,fix,h,r,rnot,z,bm,qdel,va,cl,j_state,mj,mj_state,ms
53 : REAL denz1(1,1),vacxpot(1,1),vacpot(1,1)
54 : INTEGER i,ivac,iza,j,jr,k,n,n1,ispin
55 : INTEGER nw,ilo,natot,nat,l,atomType,istate,i_u,kappa,m
56 :
57 :
58 : ! Local Arrays
59 136 : REAL, ALLOCATABLE :: vbar(:,:)
60 272 : REAL, ALLOCATABLE :: rat(:,:),eig(:,:,:),sigm(:)
61 136 : REAL, ALLOCATABLE :: rh(:,:,:),rh1(:,:,:),rhoss(:,:)
62 : REAL :: vacpar(2),occ_state(2)
63 272 : INTEGER lnum(29,atoms%ntype),nst(atoms%ntype)
64 136 : INTEGER, ALLOCATABLE :: state_indices(:)
65 136 : INTEGER jrc(atoms%ntype)
66 136 : LOGICAL l_found(0:3),llo_found(atoms%nlod),l_st
67 136 : REAL,ALLOCATABLE :: occ(:,:)
68 136 : COMPLEX,ALLOCATABLE :: pw_tmp(:,:)
69 136 : COMPLEX,ALLOCATABLE :: mmpmat_tmp(:,:,:,:)
70 :
71 544 : TYPE(t_nococonv) :: nococonv
72 : ! Data statements
73 : DATA del/1.e-6/
74 : PARAMETER (l_st=.true.)
75 :
76 : !use the init_potden_simple routine to prevent extra dimensions from noco calculations
77 : CALL den%init(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,atoms%ntype,&
78 : atoms%n_denmat,atoms%n_vPairs,input%jspins,.FALSE.,.FALSE.,POTDEN_TYPE_DEN,&
79 136 : vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
80 :
81 952 : ALLOCATE ( rat(atoms%msh,atoms%ntype),eig(29,input%jspins,atoms%ntype) )
82 1088 : ALLOCATE ( rh(atoms%msh,atoms%ntype,input%jspins),rh1(atoms%msh,atoms%ntype,input%jspins) )
83 680 : ALLOCATE ( vbar(2,atoms%ntype),sigm(vacuum%nmzd) )
84 544 : ALLOCATE ( rhoss(atoms%msh,input%jspins) )
85 :
86 357494 : rh = 0.0
87 208450 : rhoss = 0.0
88 :
89 136 : CALL timestart("stden - init")
90 :
91 136 : IF (fmpi%irank == 0) THEN
92 : ! if sigma is not 0.0, then divide this charge among all atoms
93 : !TODO: reactivate efields
94 : !IF ( ABS(input%sigma).LT. 1.e-6) THEN
95 : IF (1==1) THEN
96 68 : qdel = 0.0
97 : ELSE
98 : natot = 0
99 : DO n = 1, atoms%ntype
100 : IF (atoms%zatom(n).GE.1.0) natot = natot + atoms%neq(n)
101 : END DO
102 : !qdel = 2.*input%sigma/natot
103 : END IF
104 :
105 68 : WRITE (oUnit,FMT=8000)
106 : 8000 FORMAT (/,/,/,' superposition of atomic densities',/,/,' original atomic densities:',/)
107 189 : DO n = 1,atoms%ntype
108 121 : r = atoms%rmsh(1,n)
109 121 : d = EXP(atoms%dx(n))
110 121 : jrc(n) = 0
111 103036 : DO WHILE (r < atoms%rmt(n) + 20.0)
112 102847 : IF (jrc(n) > atoms%msh) CALL juDFT_error("increase msh in fl7para!",calledby ="stden")
113 102847 : jrc(n) = jrc(n) + 1
114 102847 : rat(jrc(n),n) = r
115 102847 : r = r*d
116 : END DO
117 : END DO
118 :
119 : ! Generate the atomic charge densities
120 189 : DO n = 1,atoms%ntype
121 121 : z = atoms%zatom(n)
122 121 : r = atoms%rmt(n)
123 121 : h = atoms%dx(n)
124 121 : jr = atoms%jri(n)
125 121 : IF (input%jspins.EQ.2) THEN
126 77 : bm = atoms%bmu(n)
127 : ELSE
128 : bm = 0.
129 : END IF
130 2958 : occ=atoms%econf(n)%Occupation(:,:)
131 : ! check whether this atom has been done already
132 150 : DO n1 = 1, n - 1
133 66 : IF (ABS(z-atoms%zatom(n1)).GT.del) CYCLE
134 38 : IF (ABS(r-atoms%rmt(n1)).GT.del) CYCLE
135 38 : IF (ABS(h-atoms%dx(n1)).GT.del) CYCLE
136 38 : IF (ABS(bm-atoms%bmu(n1)).GT.del) CYCLE
137 793 : IF (ANY(ABS(occ(:,:)-atoms%econf(n1)%Occupation(:,:))>del)) CYCLE
138 37 : IF (jr.NE.atoms%jri(n1)) CYCLE
139 97 : DO ispin = 1, input%jspins
140 49664 : DO i = 1,jrc(n) ! atoms%msh
141 49627 : rh(i,n,ispin) = rh(i,n1,ispin)
142 : END DO
143 : END DO
144 37 : nst(n) = nst(n1)
145 37 : vbar(1,n) = vbar(1,n1)
146 37 : vbar(input%jspins,n) = vbar(input%jspins,n1)
147 373 : DO i = 1, nst(n1)
148 336 : lnum(i,n) = lnum(i,n1)
149 336 : eig(i,1,n) = eig(i,1,n1)
150 373 : eig(i,input%jspins,n) = eig(i,input%jspins,n1)
151 : END DO
152 113 : GO TO 70
153 : END DO
154 : !---> new atom
155 84 : rnot = atoms%rmsh(1,n)
156 84 : IF (z.LT.1.0) THEN
157 0 : va = max(z,1.e-8)/(input%jspins*sfp_const*atoms%volmts(n))
158 0 : DO ispin = 1, input%jspins
159 0 : DO i = 1,jrc(n) ! atoms%msh
160 0 : rh(i,n,ispin) = va/rat(i,n)**2
161 : END DO
162 : END DO
163 : ELSE
164 : CALL atom2(atoms,xcpot,input,n,jrc(n),rnot,qdel,&
165 84 : rhoss,nst(n),lnum(1,n),eig(1,1,n),vbar(1,n),.true.)
166 222 : DO ispin = 1, input%jspins
167 121041 : DO i = 1, jrc(n) ! atoms%msh
168 120957 : rh(i,n,ispin) = rhoss(i,ispin)
169 : END DO
170 : END DO
171 : END IF
172 : ! list atomic density
173 84 : iza = atoms%zatom(n) + 0.0001
174 84 : WRITE (oUnit,FMT=8030) namat_const(iza)
175 : 8030 FORMAT (/,/,' atom: ',a2,/)
176 : 8040 FORMAT (4 (3x,i5,f8.5,f12.6))
177 68 : 70 CONTINUE
178 : END DO
179 :
180 : !roa+
181 : ! use cdnovlp to generate total density out of atom densities
182 184 : DO ispin = 1, input%jspins
183 382 : DO n = 1,atoms%ntype
184 198 : nat = atoms%firstAtom(n)
185 170584 : DO i = 1, jrc(n)
186 170584 : rh1(i,n,ispin) = rh(i,n,ispin)*fpi_const*rat(i,n)**2
187 : END DO
188 8375 : rh1(jrc(n):atoms%msh,n,ispin) = 0.0
189 : ! prepare spherical mt charge
190 143532 : DO i = 1,atoms%jri(n)
191 143532 : den%mt(i,0,n,ispin) = rh(i,n,ispin)*sfp_const*atoms%rmsh(i,n)**2
192 : END DO
193 : ! reset nonspherical mt charge
194 7967 : DO k = 1,sphhar%nlh(sym%ntypsy(nat))
195 5623494 : DO j = 1,atoms%jri(n)
196 5623296 : den%mt(j,k,n,ispin) = 0.e0
197 : END DO
198 : END DO
199 : END DO
200 : END DO ! ispin
201 : END IF ! fmpi%irank == 0
202 :
203 136 : CALL timestop("stden - init")
204 :
205 368 : DO ispin = 1, input%jspins
206 : CALL cdnovlp(fmpi,sphhar,stars,atoms,sym,vacuum,&
207 : cell,input ,l_st,ispin,rh1(:,:,ispin),&
208 368 : den%pw,den%mt,den%vac)
209 : !roa-
210 : END DO
211 :
212 136 : if (noco%l_noco) THEN
213 289022 : den%pw(:,1)=(den%pw(:,1)+den%pw(:,2))*0.5
214 289022 : den%pw(:,2)=den%pw(:,1)
215 46 : if (input%film) THEN
216 0 : den%vac(:,:,:,1)=(den%vac(:,:,:,1)+den%vac(:,:,:,2))*0.5
217 0 : den%vac(:,:,:,2)=den%vac(:,:,:,1)
218 : endif
219 : endif
220 :
221 : ! Check the normalization of total density
222 136 : CALL timestart("stden - qfix")
223 136 : CALL qfix(fmpi,stars,nococonv,atoms,sym,vacuum,sphhar,input,cell ,den,.FALSE.,.FALSE.,l_par=.FALSE.,force_fix=.TRUE.,fix=fix)
224 136 : CALL timestop("stden - qfix")
225 136 : CALL timestart("stden - finalize")
226 : !Rotate density into global frame if l_alignSQA
227 374 : IF (any(noco%l_alignMT)) then
228 8 : allocate(nococonv%beta(atoms%ntype),nococonv%alph(atoms%ntype))
229 8 : nococonv%beta=noco%beta_inp
230 8 : nococonv%alph=noco%alph_inp
231 2 : CALL toGlobalSpinFrame(noco, nococonv, vacuum, sphhar, stars, sym, cell, input, atoms, Den)
232 6 : Allocate(pw_tmp(size(den%pw,1),3))
233 18434 : pw_tmp=0.0
234 12290 : pw_tmp(:,:size(den%pw,2))=den%pw
235 2 : call move_alloc(pw_tmp,den%pw)
236 : ENDIF
237 136 : IF(input%ldauSpinoffd) THEN
238 0 : allocate(mmpmat_tmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,SIZE(den%mmpmat,3),3))
239 0 : mmpmat_tmp = cmplx_0
240 0 : mmpmat_tmp(:,:,:,:size(den%mmpmat,4)) = den%mmpmat
241 0 : call move_alloc(mmpmat_tmp,den%mmpmat)
242 : ENDIF
243 :
244 :
245 136 : IF(atoms%n_u>0.and.fmpi%irank==0.and.input%ldauInitialGuess) THEN
246 : !Create initial guess for the density matrix based on the occupations in the inp.xml file
247 3 : WRITE(*,*) "Creating initial guess for LDA+U density matrix"
248 351 : den%mmpMat = cmplx_0
249 6 : do i_u = 1, atoms%n_u
250 3 : l = atoms%lda_u(i_u)%l
251 3 : atomType = atoms%lda_u(i_u)%atomType
252 :
253 66 : state_indices = atoms%econf(atomType)%get_states_for_orbital(l)
254 9 : do istate = 1, atoms%econf(atomType)%num_states
255 9 : if(state_indices(istate)<0) exit
256 18 : occ_state=atoms%econf(atomType)%occupation(state_indices(istate),:)
257 6 : kappa = atoms%econf(atomType)%kappa(state_indices(istate))
258 :
259 21 : DO ispin=1,2
260 :
261 12 : j_state = abs(kappa)-0.5
262 12 : mj_state = -j_state
263 102 : DO WHILE(mj_state <= MERGE(l,l+1,kappa>0))
264 84 : ms = MERGE(0.5,-0.5,ispin==1)
265 84 : mj = mj_state * SIGN(1.,ms)
266 84 : m = mj - ms
267 84 : IF(ABS(m)<=l) then
268 78 : cl = clebsch(REAL(l),0.5,REAL(m),ms,j_state,mj)**2
269 78 : den%mmpMat(m,m,i_u,MIN(ispin,input%jspins)) = den%mmpMat(m,m,i_u,MIN(ispin,input%jspins)) + MIN(cl,occ_state(ispin))
270 78 : occ_state(ispin) = MAX(occ_state(ispin)-cl,0.0)
271 : endif
272 84 : mj_state = mj_state + 1
273 : ENDDO
274 : ENDDO
275 : enddo
276 3 : IF(.not.noco%l_soc) THEN
277 : !Average -m and m
278 0 : DO ispin = 1, input%jspins
279 0 : DO m = 1,l
280 0 : den%mmpMat(m,m,i_u,ispin) = (den%mmpMat(m,m,i_u,ispin) + den%mmpMat(-m,-m,i_u,ispin))/2.0
281 0 : den%mmpMat(-m,-m,i_u,ispin) = den%mmpMat(m,m,i_u,ispin)
282 : enddo
283 : enddo
284 : endif
285 : !Rotate into the global real frame
286 6 : IF(noco%l_noco) THEN
287 3 : do ispin =1, input%jspins
288 115 : den%mmpMat(:,:,i_u,ispin) = rotMMPmat(den%mmpMat(:,:,i_u,ispin),noco%alph_inp(atomType),noco%beta_inp(atomType),0.0,l)
289 : enddo
290 2 : ELSE IF(noco%l_soc) THEN
291 6 : do ispin =1, input%jspins
292 230 : den%mmpMat(:,:,i_u,ispin) = rotMMPmat(den%mmpMat(:,:,i_u,ispin),noco%phi_inp,noco%theta_inp,0.0,l)
293 : enddo
294 : ENDIF
295 : enddo
296 : ENDIF
297 :
298 136 : IF (fmpi%irank == 0) THEN
299 189 : z=SUM(atoms%neq(:)*atoms%zatom(:))
300 68 : IF (ABS(fix*z-z)>0.5) THEN
301 : CALL judft_warn("Starting density not charge neutral",hint= &
302 0 : "Your electronic configuration might be broken",calledby="stden.f90")
303 : END IF
304 :
305 : ! Write superposed density onto density file
306 68 : den%iter = 0
307 :
308 68 : if (noco%l_noco) THEN
309 144511 : den%pw(:,1)=(den%pw(:,1)+den%pw(:,2))*0.5
310 144511 : den%pw(:,2)=den%pw(:,1)
311 : endif
312 68 : CALL timestart("stden - write density")
313 : CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,&
314 : merge(CDN_ARCHIVE_TYPE_FFN_const,CDN_ARCHIVE_TYPE_CDN1_const,any(noco%l_alignMT)),&
315 254 : CDN_INPUT_DEN_const,1,-1.0,0.0,-1.0,-1.0,.TRUE.,den)
316 68 : CALL timestop("stden - write density")
317 : ! Check continuity
318 68 : IF (input%vchk) THEN
319 8 : DO ispin = 1, input%jspins
320 4 : WRITE (oUnit,'(a8,i2)') 'spin No.',ispin
321 : CALL checkDOPAll(input,sphhar,stars,atoms,sym,vacuum ,&
322 8 : cell,den,ispin)
323 : END DO ! ispin = 1, input%jspins
324 : END IF ! input%vchk
325 :
326 : ! set up parameters for enpara-file
327 68 : IF ((juDFT_was_argument("-genEnpara"))) THEN
328 0 : CALL enpara%init_enpara(atoms,input%jspins,input%film)
329 :
330 0 : enpara%lchange = .TRUE.
331 0 : enpara%llochg = .TRUE.
332 :
333 0 : DO ispin = 1, input%jspins
334 : ! vacpar is taken as highest occupied level from atomic eigenvalues
335 : ! vacpar (+0.3) serves as fermi level for film or bulk
336 0 : vacpar(1) = -999.9
337 0 : DO n = 1,atoms%ntype
338 0 : vacpar(1) = MAX( vacpar(1),eig(nst(n),ispin,n) )
339 : END DO
340 0 : IF (.NOT.input%film) vacpar(1) = vacpar(1) + 0.4
341 0 : vacpar(2) = vacpar(1)
342 0 : DO n = 1,atoms%ntype
343 0 : enpara%skiplo(n,ispin) = 0
344 0 : DO i = 0, 3
345 0 : l_found(i) = .FALSE.
346 0 : enpara%el0(i,n,ispin) = vacpar(1)
347 : END DO
348 0 : DO i = 1, atoms%nlo(n)
349 0 : llo_found(i) = .FALSE.
350 0 : enpara%ello0(i,n,ispin) = vacpar(1) - 1.5
351 : END DO
352 :
353 : ! take energy-parameters from atomic calculation
354 0 : DO i = nst(n), 1, -1
355 0 : IF (.NOT.input%film) eig(i,ispin,n) = eig(i,ispin,n) + 0.4
356 0 : IF (.NOT.l_found(lnum(i,n)).AND.(lnum(i,n).LE.3)) THEN
357 0 : enpara%el0(lnum(i,n),n,ispin) = eig(i,ispin,n)
358 0 : IF (enpara%el0(lnum(i,n),n,ispin).LT.input_ellow) THEN
359 0 : enpara%el0(lnum(i,n),n,ispin) = vacpar(1)
360 0 : l_found(lnum(i,n)) = .TRUE.
361 : END IF
362 0 : IF (enpara%el0(lnum(i,n),n,ispin).LT.input_elup) THEN
363 0 : l_found(lnum(i,n)) = .TRUE.
364 : END IF
365 : ELSE
366 0 : IF (l_found(lnum(i,n)).AND.(atoms%nlo(n).GT.0)) THEN
367 0 : DO ilo = 1, atoms%nlo(n)
368 0 : IF (atoms%llo(ilo,n).EQ.lnum(i,n)) THEN
369 0 : IF (.NOT.llo_found(ilo)) THEN
370 0 : enpara%ello0(ilo,n,ispin) = eig(i,ispin,n)
371 0 : IF ((enpara%ello0(ilo,n,ispin).GT.input_elup).OR.&
372 : (enpara%ello0(ilo,n,ispin).LT.input_ellow)) THEN
373 0 : enpara%ello0(ilo,n,ispin)= vacpar(1)
374 0 : ELSE IF (atoms%l_dulo(ilo,n)) THEN
375 0 : enpara%ello0(ilo,n,ispin)= enpara%el0(atoms%llo(ilo,n),n,ispin)
376 : ELSE
377 0 : enpara%skiplo(n,ispin) = enpara%skiplo(n,ispin) + 2*atoms%llo(ilo,n)+1
378 : END IF
379 0 : llo_found(ilo) = .TRUE.
380 : IF ((enpara%el0(atoms%llo(ilo,n),n,ispin)-enpara%ello0(ilo,n,ispin).LT.0.5)&
381 0 : .AND.(atoms%llo(ilo,n).GE.0)) THEN
382 0 : enpara%el0(atoms%llo(ilo,n),n,ispin) = vacpar(1)
383 0 : IF (atoms%l_dulo(ilo,n)) enpara%ello0(ilo,n,ispin)= enpara%el0(atoms%llo(ilo,n),n,ispin)
384 : END IF
385 : END IF
386 : END IF
387 : END DO ! ilo = 1, atoms%nlo(n)
388 : END IF
389 : END IF ! .NOT.l_found(lnum(i,n)).AND.(lnum(i,n).LE.3)
390 : END DO ! i = nst(n), 1, -1
391 : END DO ! atom types
392 :
393 0 : IF (input%film) THEN
394 : ! get guess for vacuum parameters
395 : ! YM : in 1D case should be modified, but also not that bad like this
396 : !
397 : ! generate coulomb potential by integrating inward to z1
398 :
399 0 : DO ivac = 1, vacuum%nvac
400 0 : DO i=1,vacuum%nmz
401 0 : sigm(i) = (i-1)*vacuum%delz*REAL(den%vac(i,1,ivac,ispin))
402 : END DO
403 0 : CALL qsf(vacuum%delz,sigm,vacpar(ivac),vacuum%nmz,0)
404 0 : denz1 = REAL(den%vac(1,1,ivac,ispin)) ! get estimate for potential at vacuum boundary
405 0 : CALL xcpot%get_vxc(1,denz1,vacpot,vacxpot)
406 : ! seems to be the best choice for 1D not to substract vacpar
407 0 : vacpot = vacpot - fpi_const*vacpar(ivac)
408 0 : vacpar(ivac) = vacpot(1,1)
409 : END DO
410 0 : IF (vacuum%nvac.EQ.1) vacpar(2) = vacpar(1)
411 : END IF
412 :
413 0 : enpara%enmix = 1.0
414 :
415 :
416 0 : enpara%evac0(:,ispin)=vacpar(:SIZE(enpara%evac0,1))
417 :
418 : END DO ! ispin
419 0 : CALL enpara%WRITE(atoms,input%jspins,input%film)
420 : END IF
421 : END IF ! fmpi%irank == 0
422 :
423 136 : DEALLOCATE ( rat,eig )
424 136 : DEALLOCATE ( rh,rh1)
425 136 : DEALLOCATE ( vbar,sigm )
426 136 : DEALLOCATE ( rhoss )
427 136 : deallocate(den%vac)
428 :
429 136 : CALL timestop("stden - finalize")
430 :
431 136 : END SUBROUTINE stden
432 :
433 : END MODULE m_stden
|
