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 111 : 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 111 : TYPE(t_potden) :: den
50 111 : 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 111 : REAL, ALLOCATABLE :: vbar(:,:)
60 222 : REAL, ALLOCATABLE :: rat(:,:),eig(:,:,:),sigm(:)
61 111 : REAL, ALLOCATABLE :: rh(:,:,:),rh1(:,:,:),rhoss(:,:)
62 : REAL :: vacpar(2),occ_state(2)
63 222 : INTEGER lnum(29,atoms%ntype),nst(atoms%ntype)
64 111 : INTEGER, ALLOCATABLE :: state_indices(:)
65 111 : INTEGER jrc(atoms%ntype)
66 111 : LOGICAL l_found(0:3),llo_found(atoms%nlod),l_st
67 111 : REAL,ALLOCATABLE :: occ(:,:)
68 111 : COMPLEX,ALLOCATABLE :: pw_tmp(:,:)
69 111 : COMPLEX,ALLOCATABLE :: mmpmat_tmp(:,:,:,:)
70 :
71 444 : 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 111 : vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
80 :
81 444 : ALLOCATE ( rat(atoms%msh,atoms%ntype),eig(29,input%jspins,atoms%ntype) )
82 444 : ALLOCATE ( rh(atoms%msh,atoms%ntype,input%jspins),rh1(atoms%msh,atoms%ntype,input%jspins) )
83 111 : ALLOCATE ( vbar(2,atoms%ntype),sigm(vacuum%nmzd) )
84 222 : ALLOCATE ( rhoss(atoms%msh,input%jspins) )
85 :
86 287718 : rh = 0.0
87 153324 : rhoss = 0.0
88 :
89 111 : CALL timestart("stden - init")
90 :
91 111 : 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 57 : 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 57 : WRITE (oUnit,FMT=8000)
106 : 8000 FORMAT (/,/,/,' superposition of atomic densities',/,/,' original atomic densities:',/)
107 165 : DO n = 1,atoms%ntype
108 108 : r = atoms%rmsh(1,n)
109 108 : d = EXP(atoms%dx(n))
110 108 : jrc(n) = 0
111 90803 : DO WHILE (r < atoms%rmt(n) + 20.0)
112 90638 : IF (jrc(n) > atoms%msh) CALL juDFT_error("increase msh in fl7para!",calledby ="stden")
113 90638 : jrc(n) = jrc(n) + 1
114 90638 : rat(jrc(n),n) = r
115 90638 : r = r*d
116 : END DO
117 : END DO
118 :
119 : ! Generate the atomic charge densities
120 165 : DO n = 1,atoms%ntype
121 108 : z = atoms%zatom(n)
122 108 : r = atoms%rmt(n)
123 108 : h = atoms%dx(n)
124 108 : jr = atoms%jri(n)
125 108 : IF (input%jspins.EQ.2) THEN
126 56 : bm = atoms%bmu(n)
127 : ELSE
128 : bm = 0.
129 : END IF
130 2514 : occ=atoms%econf(n)%Occupation(:,:)
131 : ! check whether this atom has been done already
132 136 : DO n1 = 1, n - 1
133 64 : IF (ABS(z-atoms%zatom(n1)).GT.del) CYCLE
134 37 : IF (ABS(r-atoms%rmt(n1)).GT.del) CYCLE
135 37 : IF (ABS(h-atoms%dx(n1)).GT.del) CYCLE
136 37 : IF (ABS(bm-atoms%bmu(n1)).GT.del) CYCLE
137 746 : IF (ANY(ABS(occ(:,:)-atoms%econf(n1)%Occupation(:,:))>del)) CYCLE
138 36 : IF (jr.NE.atoms%jri(n1)) CYCLE
139 91 : DO ispin = 1, input%jspins
140 45481 : DO i = 1,jrc(n) ! atoms%msh
141 45445 : rh(i,n,ispin) = rh(i,n1,ispin)
142 : END DO
143 : END DO
144 36 : nst(n) = nst(n1)
145 36 : vbar(1,n) = vbar(1,n1)
146 36 : vbar(input%jspins,n) = vbar(input%jspins,n1)
147 350 : DO i = 1, nst(n1)
148 314 : lnum(i,n) = lnum(i,n1)
149 314 : eig(i,1,n) = eig(i,1,n1)
150 350 : eig(i,input%jspins,n) = eig(i,input%jspins,n1)
151 : END DO
152 100 : GO TO 70
153 : END DO
154 : !---> new atom
155 72 : rnot = atoms%rmsh(1,n)
156 72 : 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 72 : rhoss,nst(n),lnum(1,n),eig(1,1,n),vbar(1,n),.true.)
166 181 : DO ispin = 1, input%jspins
167 92868 : DO i = 1, jrc(n) ! atoms%msh
168 92796 : rh(i,n,ispin) = rhoss(i,ispin)
169 : END DO
170 : END DO
171 : END IF
172 : ! list atomic density
173 72 : iza = atoms%zatom(n) + 0.0001
174 72 : WRITE (oUnit,FMT=8030) namat_const(iza)
175 : 8030 FORMAT (/,/,' atom: ',a2,/)
176 : 8040 FORMAT (4 (3x,i5,f8.5,f12.6))
177 57 : 70 CONTINUE
178 : END DO
179 :
180 : !roa+
181 : ! use cdnovlp to generate total density out of atom densities
182 146 : DO ispin = 1, input%jspins
183 310 : DO n = 1,atoms%ntype
184 164 : nat = atoms%firstAtom(n)
185 138241 : DO i = 1, jrc(n)
186 138241 : rh1(i,n,ispin) = rh(i,n,ispin)*fpi_const*rat(i,n)**2
187 : END DO
188 8249 : rh1(jrc(n):atoms%msh,n,ispin) = 0.0
189 : ! prepare spherical mt charge
190 115848 : DO i = 1,atoms%jri(n)
191 115848 : 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 6778 : DO k = 1,sphhar%nlh(sym%ntypsy(nat))
195 4698816 : DO j = 1,atoms%jri(n)
196 4698652 : 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 111 : CALL timestop("stden - init")
204 :
205 111 : if (noco%l_noco) THEN
206 46400 : rh1(:,:,1)=0.5*(rh1(:,:,1)+rh1(:,:,2))
207 46400 : rh1(:,:,2)=rh1(:,:,1)
208 : ENDIF
209 286 : DO ispin = 1, input%jspins
210 : CALL cdnovlp(fmpi,sphhar,stars,atoms,sym,vacuum,&
211 : cell,input ,l_st,ispin,rh1(:,:,ispin),&
212 286 : den%pw,den%mt,den%vac,.TRUE.)
213 : !roa-
214 : END DO
215 :
216 111 : if (noco%l_noco) THEN
217 267028 : den%pw(:,1)=(den%pw(:,1)+den%pw(:,2))*0.5
218 267028 : den%pw(:,2)=den%pw(:,1)
219 40 : if (input%film) THEN
220 0 : den%vac(:,:,:,1)=(den%vac(:,:,:,1)+den%vac(:,:,:,2))*0.5
221 0 : den%vac(:,:,:,2)=den%vac(:,:,:,1)
222 : endif
223 : endif
224 :
225 : ! Check the normalization of total density
226 111 : CALL timestart("stden - qfix")
227 111 : CALL qfix(fmpi,stars,nococonv,atoms,sym,vacuum,sphhar,input,cell ,den,.FALSE.,.FALSE.,l_par=.FALSE.,force_fix=.TRUE.,fix=fix)
228 111 : CALL timestop("stden - qfix")
229 111 : CALL timestart("stden - finalize")
230 : !Rotate density into global frame if l_alignSQA
231 318 : IF (any(noco%l_alignMT)) then
232 2 : allocate(nococonv%beta(atoms%ntype),nococonv%alph(atoms%ntype))
233 8 : nococonv%beta=noco%beta_inp
234 8 : nococonv%alph=noco%alph_inp
235 2 : CALL toGlobalSpinFrame(noco, nococonv, vacuum, sphhar, stars, sym, cell, input, atoms, Den)
236 2 : Allocate(pw_tmp(size(den%pw,1),3))
237 18434 : pw_tmp=0.0
238 12290 : pw_tmp(:,:size(den%pw,2))=den%pw
239 2 : call move_alloc(pw_tmp,den%pw)
240 : ENDIF
241 111 : IF(input%ldauSpinoffd) THEN
242 0 : allocate(mmpmat_tmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,SIZE(den%mmpmat,3),3))
243 0 : mmpmat_tmp = cmplx_0
244 0 : mmpmat_tmp(:,:,:,:size(den%mmpmat,4)) = den%mmpmat
245 0 : call move_alloc(mmpmat_tmp,den%mmpmat)
246 : ENDIF
247 :
248 :
249 111 : IF(atoms%n_u>0.and.fmpi%irank==0.and.input%ldauInitialGuess) THEN
250 : !Create initial guess for the density matrix based on the occupations in the inp.xml file
251 0 : WRITE(*,*) "Creating initial guess for LDA+U density matrix"
252 0 : den%mmpMat = cmplx_0
253 0 : do i_u = 1, atoms%n_u
254 0 : l = atoms%lda_u(i_u)%l
255 0 : atomType = atoms%lda_u(i_u)%atomType
256 :
257 0 : state_indices = atoms%econf(atomType)%get_states_for_orbital(l)
258 0 : do istate = 1, atoms%econf(atomType)%num_states
259 0 : if(state_indices(istate)<0) exit
260 0 : occ_state=atoms%econf(atomType)%occupation(state_indices(istate),:)
261 0 : kappa = atoms%econf(atomType)%kappa(state_indices(istate))
262 :
263 0 : DO ispin=1,2
264 :
265 0 : j_state = abs(kappa)-0.5
266 0 : mj_state = -j_state
267 0 : DO WHILE(mj_state <= MERGE(l,l+1,kappa>0))
268 0 : ms = MERGE(0.5,-0.5,ispin==1)
269 0 : mj = mj_state * SIGN(1.,ms)
270 0 : m = mj - ms
271 0 : IF(ABS(m)<=l) then
272 0 : cl = clebsch(REAL(l),0.5,REAL(m),ms,j_state,mj)**2
273 0 : 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))
274 0 : occ_state(ispin) = MAX(occ_state(ispin)-cl,0.0)
275 : endif
276 0 : mj_state = mj_state + 1
277 : ENDDO
278 : ENDDO
279 : enddo
280 0 : IF(.not.noco%l_soc) THEN
281 : !Average -m and m
282 0 : DO ispin = 1, input%jspins
283 0 : DO m = 1,l
284 0 : den%mmpMat(m,m,i_u,ispin) = (den%mmpMat(m,m,i_u,ispin) + den%mmpMat(-m,-m,i_u,ispin))/2.0
285 0 : den%mmpMat(-m,-m,i_u,ispin) = den%mmpMat(m,m,i_u,ispin)
286 : enddo
287 : enddo
288 : endif
289 : !Rotate into the global real frame
290 0 : IF(noco%l_noco) THEN
291 0 : do ispin =1, input%jspins
292 0 : den%mmpMat(:,:,i_u,ispin) = rotMMPmat(den%mmpMat(:,:,i_u,ispin),noco%alph_inp(atomType),noco%beta_inp(atomType),0.0,l)
293 : enddo
294 0 : ELSE IF(noco%l_soc) THEN
295 0 : do ispin =1, input%jspins
296 0 : den%mmpMat(:,:,i_u,ispin) = rotMMPmat(den%mmpMat(:,:,i_u,ispin),noco%phi_inp,noco%theta_inp,0.0,l)
297 : enddo
298 : ENDIF
299 : enddo
300 : ENDIF
301 :
302 111 : IF (fmpi%irank == 0) THEN
303 165 : z=SUM(atoms%neq(:)*atoms%zatom(:))
304 57 : IF (ABS(fix*z-z)>0.5) THEN
305 : CALL judft_warn("Starting density not charge neutral",hint= &
306 0 : "Your electronic configuration might be broken",calledby="stden.f90")
307 : END IF
308 :
309 : ! Write superposed density onto density file
310 57 : den%iter = 0
311 :
312 57 : if (noco%l_noco) THEN
313 133514 : den%pw(:,1)=(den%pw(:,1)+den%pw(:,2))*0.5
314 133514 : den%pw(:,2)=den%pw(:,1)
315 : endif
316 57 : CALL timestart("stden - write density")
317 : CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,&
318 : merge(CDN_ARCHIVE_TYPE_FFN_const,CDN_ARCHIVE_TYPE_CDN1_const,any(noco%l_alignMT)),&
319 219 : CDN_INPUT_DEN_const,1,-1.0,0.0,-1.0,-1.0,.TRUE.,den)
320 57 : CALL timestop("stden - write density")
321 : ! Check continuity
322 57 : IF (input%vchk) THEN
323 10 : DO ispin = 1, input%jspins
324 5 : WRITE (oUnit,'(a8,i2)') 'spin No.',ispin
325 : CALL checkDOPAll(input,sphhar,stars,atoms,sym,vacuum ,&
326 10 : cell,den,ispin)
327 : END DO ! ispin = 1, input%jspins
328 : END IF ! input%vchk
329 :
330 : ! set up parameters for enpara-file
331 57 : IF ((juDFT_was_argument("-genEnpara"))) THEN
332 0 : CALL enpara%init_enpara(atoms,input%jspins,input%film)
333 :
334 0 : enpara%lchange = .TRUE.
335 0 : enpara%llochg = .TRUE.
336 :
337 0 : DO ispin = 1, input%jspins
338 : ! vacpar is taken as highest occupied level from atomic eigenvalues
339 : ! vacpar (+0.3) serves as fermi level for film or bulk
340 0 : vacpar(1) = -999.9
341 0 : DO n = 1,atoms%ntype
342 0 : vacpar(1) = MAX( vacpar(1),eig(nst(n),ispin,n) )
343 : END DO
344 0 : IF (.NOT.input%film) vacpar(1) = vacpar(1) + 0.4
345 0 : vacpar(2) = vacpar(1)
346 0 : DO n = 1,atoms%ntype
347 0 : enpara%skiplo(n,ispin) = 0
348 0 : DO i = 0, 3
349 0 : l_found(i) = .FALSE.
350 0 : enpara%el0(i,n,ispin) = vacpar(1)
351 : END DO
352 0 : DO i = 1, atoms%nlo(n)
353 0 : llo_found(i) = .FALSE.
354 0 : enpara%ello0(i,n,ispin) = vacpar(1) - 1.5
355 : END DO
356 :
357 : ! take energy-parameters from atomic calculation
358 0 : DO i = nst(n), 1, -1
359 0 : IF (.NOT.input%film) eig(i,ispin,n) = eig(i,ispin,n) + 0.4
360 0 : IF (.NOT.l_found(lnum(i,n)).AND.(lnum(i,n).LE.3)) THEN
361 0 : enpara%el0(lnum(i,n),n,ispin) = eig(i,ispin,n)
362 0 : IF (enpara%el0(lnum(i,n),n,ispin).LT.input_ellow) THEN
363 0 : enpara%el0(lnum(i,n),n,ispin) = vacpar(1)
364 0 : l_found(lnum(i,n)) = .TRUE.
365 : END IF
366 0 : IF (enpara%el0(lnum(i,n),n,ispin).LT.input_elup) THEN
367 0 : l_found(lnum(i,n)) = .TRUE.
368 : END IF
369 : ELSE
370 0 : IF (l_found(lnum(i,n)).AND.(atoms%nlo(n).GT.0)) THEN
371 0 : DO ilo = 1, atoms%nlo(n)
372 0 : IF (atoms%llo(ilo,n).EQ.lnum(i,n)) THEN
373 0 : IF (.NOT.llo_found(ilo)) THEN
374 0 : enpara%ello0(ilo,n,ispin) = eig(i,ispin,n)
375 0 : IF ((enpara%ello0(ilo,n,ispin).GT.input_elup).OR.&
376 : (enpara%ello0(ilo,n,ispin).LT.input_ellow)) THEN
377 0 : enpara%ello0(ilo,n,ispin)= vacpar(1)
378 0 : ELSE IF (atoms%l_dulo(ilo,n)) THEN
379 0 : enpara%ello0(ilo,n,ispin)= enpara%el0(atoms%llo(ilo,n),n,ispin)
380 : ELSE
381 0 : enpara%skiplo(n,ispin) = enpara%skiplo(n,ispin) + 2*atoms%llo(ilo,n)+1
382 : END IF
383 0 : llo_found(ilo) = .TRUE.
384 : IF ((enpara%el0(atoms%llo(ilo,n),n,ispin)-enpara%ello0(ilo,n,ispin).LT.0.5)&
385 0 : .AND.(atoms%llo(ilo,n).GE.0)) THEN
386 0 : enpara%el0(atoms%llo(ilo,n),n,ispin) = vacpar(1)
387 0 : IF (atoms%l_dulo(ilo,n)) enpara%ello0(ilo,n,ispin)= enpara%el0(atoms%llo(ilo,n),n,ispin)
388 : END IF
389 : END IF
390 : END IF
391 : END DO ! ilo = 1, atoms%nlo(n)
392 : END IF
393 : END IF ! .NOT.l_found(lnum(i,n)).AND.(lnum(i,n).LE.3)
394 : END DO ! i = nst(n), 1, -1
395 : END DO ! atom types
396 :
397 0 : IF (input%film) THEN
398 : ! get guess for vacuum parameters
399 : ! YM : in 1D case should be modified, but also not that bad like this
400 : !
401 : ! generate coulomb potential by integrating inward to z1
402 :
403 0 : DO ivac = 1, vacuum%nvac
404 0 : DO i=1,vacuum%nmz
405 0 : sigm(i) = (i-1)*vacuum%delz*REAL(den%vac(i,1,ivac,ispin))
406 : END DO
407 0 : CALL qsf(vacuum%delz,sigm,vacpar(ivac),vacuum%nmz,0)
408 0 : denz1 = REAL(den%vac(1,1,ivac,ispin)) ! get estimate for potential at vacuum boundary
409 0 : CALL xcpot%get_vxc(1,denz1,vacpot,vacxpot)
410 : ! seems to be the best choice for 1D not to substract vacpar
411 0 : vacpot = vacpot - fpi_const*vacpar(ivac)
412 0 : vacpar(ivac) = vacpot(1,1)
413 : END DO
414 0 : IF (vacuum%nvac.EQ.1) vacpar(2) = vacpar(1)
415 : END IF
416 :
417 0 : enpara%enmix = 1.0
418 :
419 :
420 0 : enpara%evac0(:,ispin)=vacpar(:SIZE(enpara%evac0,1))
421 :
422 : END DO ! ispin
423 0 : CALL enpara%WRITE(atoms,input%jspins,input%film)
424 : END IF
425 : END IF ! fmpi%irank == 0
426 :
427 111 : DEALLOCATE ( rat,eig )
428 111 : DEALLOCATE ( rh,rh1)
429 111 : DEALLOCATE ( vbar,sigm )
430 111 : DEALLOCATE ( rhoss )
431 111 : deallocate(den%vac)
432 :
433 111 : CALL timestop("stden - finalize")
434 :
435 111 : END SUBROUTINE stden
436 :
437 : END MODULE m_stden
|
