Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_types_enpara
8 : USE m_judft
9 : use m_types_enparaxml
10 : IMPLICIT NONE
11 : PRIVATE
12 : TYPE,extends(t_enparaxml):: t_enpara
13 : REAL, ALLOCATABLE :: el0(:,:,:)
14 : REAL :: evac(2,2)
15 : REAL, ALLOCATABLE :: ello0(:,:,:)
16 : REAL, ALLOCATABLE :: el1(:,:,:)
17 : REAL :: evac1(2,2)
18 : REAL, ALLOCATABLE :: ello1(:,:,:)
19 : REAL, ALLOCATABLE :: enmix(:)
20 : INTEGER, ALLOCATABLE :: skiplo(:,:)
21 : LOGICAL, ALLOCATABLE :: lchange(:,:,:)
22 : LOGICAL, ALLOCATABLE :: lchg_v(:,:)
23 : LOGICAL, ALLOCATABLE :: llochg(:,:,:)
24 : REAL :: epara_min
25 : LOGICAL :: ready ! are the enpara's ok for calculation?
26 : LOGICAL :: floating !floating energy parameters are relative to potential
27 : CONTAINS
28 : PROCEDURE :: init_enpara
29 : PROCEDURE :: update
30 : PROCEDURE :: read
31 : PROCEDURE :: write
32 : PROCEDURE :: mix
33 : #ifndef CPP_INPGEN
34 : PROCEDURE :: calcOutParams
35 : #endif
36 : END TYPE t_enpara
37 :
38 :
39 :
40 : PUBLIC:: t_enpara
41 :
42 : CONTAINS
43 160 : SUBROUTINE init_enpara(this,atoms,jspins,film,enparaXML)
44 : USE m_types_setup
45 : USE m_constants
46 : CLASS(t_enpara),INTENT(inout):: this
47 : TYPE(t_atoms),INTENT(IN) :: atoms
48 : INTEGER,INTENT(IN) :: jspins
49 : LOGICAL,INTENT(IN) :: film
50 : TYPE(t_enparaxml),OPTIONAL :: enparaxml
51 :
52 : LOGICAL :: l_enpara
53 :
54 1280 : ALLOCATE(this%el0(0:atoms%lmaxd,atoms%ntype,jspins),this%el1(0:atoms%lmaxd,atoms%ntype,jspins))
55 1280 : ALLOCATE(this%ello0(atoms%nlod,atoms%ntype,jspins),this%ello1(atoms%nlod,atoms%ntype,jspins))
56 5148 : this%el0=-1E99
57 1458 : this%ello0=-1E99
58 1120 : this%evac0=-1E99
59 :
60 :
61 800 : ALLOCATE(this%llochg(atoms%nlod,atoms%ntype,jspins))
62 480 : ALLOCATE(this%lchg_v(2,jspins))
63 640 : ALLOCATE(this%skiplo(atoms%ntype,jspins))
64 800 : ALLOCATE(this%lchange(0:atoms%lmaxd,atoms%ntype,jspins))
65 7570 : this%llochg=.FALSE.;this%lchg_v=.FALSE.;this%lchange=.FALSE.
66 868 : this%skiplo=0
67 480 : ALLOCATE(this%enmix(jspins))
68 428 : this%enmix=0.0
69 :
70 160 : this%ready=.FALSE.
71 160 : this%floating=.FALSE.
72 :
73 640 : ALLOCATE(this%qn_el(0:3,atoms%ntype,jspins))
74 640 : ALLOCATE(this%qn_ello(atoms%nlod,atoms%ntype,jspins))
75 :
76 1120 : this%evac0=eVac0Default_const
77 :
78 160 : inquire(file="enpara",exist=l_enpara)
79 160 : if (l_enpara) then
80 16 : call this%read(atoms,jspins,film,.true.)
81 : else
82 144 : IF (PRESENT(enparaxml)) THEN
83 2534 : this%qn_el=enparaxml%qn_el
84 1516 : this%qn_ello=enparaxml%qn_ello
85 1008 : this%evac0=enparaxml%evac0
86 : END IF
87 : endif
88 :
89 :
90 160 : END SUBROUTINE init_enpara
91 :
92 : !> This subroutine adjusts the energy parameters to the potential. In particular, it
93 : !! calculated them in case of qn_el>-1,qn_ello>-1
94 : !! Before this was done in lodpot.F
95 686 : SUBROUTINE update(enpara,fmpi,atoms,vacuum,input,v,hub1data)
96 : USE m_types_atoms
97 : USE m_types_vacuum
98 : USE m_types_input
99 : USE m_constants
100 : USE m_xmlOutput
101 : USE m_types_potden
102 : USE m_types_hub1data
103 : USE m_find_enpara
104 : USE m_types_parallelLoop
105 : USE m_types_mpi
106 : USE m_mpi_reduce_tool
107 : USE m_mpi_bc_tool
108 : #ifdef CPP_MPI
109 : USE mpi
110 : #endif
111 :
112 : CLASS(t_enpara),INTENT(inout):: enpara
113 : TYPE(t_mpi),INTENT(IN) :: fmpi
114 : TYPE(t_atoms),INTENT(IN) :: atoms
115 : TYPE(t_input),INTENT(IN) :: input
116 : TYPE(t_vacuum),INTENT(IN) :: vacuum
117 : TYPE(t_potden),INTENT(IN) :: v
118 : TYPE(t_hub1data),INTENT(IN) :: hub1data
119 :
120 :
121 : LOGICAL :: l_enpara
122 686 : LOGICAL :: l_done(0:atoms%lmaxd,atoms%ntype,input%jspins)
123 686 : LOGICAL :: lo_done(atoms%nlod,atoms%ntype,input%jspins)
124 686 : LOGICAL :: l_doneLocal(0:atoms%lmaxd,atoms%ntype,input%jspins)
125 686 : LOGICAL :: lo_doneLocal(atoms%nlod,atoms%ntype,input%jspins)
126 : REAL :: vbar,vz0,rj,tr
127 : INTEGER :: n,jsp,l,ilo,j,ivac,irank,i,ierr
128 : CHARACTER(LEN=20) :: attributes(5)
129 686 : REAL :: e_lo(0:3,atoms%ntype)!Store info on branches to do IO after OMP
130 686 : REAL :: e_up(0:3,atoms%ntype)
131 686 : REAL :: elo_lo(atoms%nlod,atoms%ntype)
132 686 : REAL :: elo_up(atoms%nlod,atoms%ntype)
133 686 : REAL :: e_lo_local(0:3,atoms%ntype)!Store info on branches to do IO after OMP
134 686 : REAL :: e_up_local(0:3,atoms%ntype)
135 686 : REAL :: elo_lo_local(atoms%nlod,atoms%ntype)
136 686 : REAL :: elo_up_local(atoms%nlod,atoms%ntype)
137 686 : REAL :: el0Local(0:atoms%lmaxd,atoms%ntype,input%jspins)
138 686 : REAL :: ello0Local(atoms%nlod,atoms%ntype,input%jspins)
139 686 : REAL,ALLOCATABLE :: v_mt(:)
140 : TYPE(t_parallelLoop) mpiLoop
141 :
142 1372 : IF (fmpi%irank == 0) CALL openXMLElement('energyParameters',(/'units'/),(/'Htr'/))
143 :
144 21396 : el0Local = 0.0
145 5616 : ello0Local = 0.0
146 21396 : l_doneLocal = .FALSE.
147 5616 : lo_doneLocal =.FALSE.
148 1768 : DO jsp = 1,input%jspins
149 10452 : e_lo_local = 0.0
150 10452 : e_up_local = 0.0
151 4930 : elo_lo_local = 0.0
152 4930 : elo_up_local = 0.0
153 :
154 1082 : CALL mpiLoop%init(fmpi%irank,fmpi%isize,1,atoms%ntype)
155 : !$OMP PARALLEL DO DEFAULT(none) &
156 : !$OMP SHARED(mpiLoop,atoms,enpara,jsp,l_doneLocal,v,lo_doneLocal,e_lo_local,e_up_local,elo_lo_local,elo_up_local,el0Local,ello0Local) &
157 1082 : !$OMP PRIVATE(n,l,ilo,v_mt)
158 : !! First calculate energy parameter from quantum numbers if these are given...
159 : !! l_done stores the index of those energy parameter updated
160 : DO n = mpiLoop%bunchMinIndex, mpiLoop%bunchMaxIndex
161 : v_mt=v%mt(:,0,n,jsp)
162 : if (atoms%l_nonpolbas(n)) v_mt=(v%mt(:,0,n,1)+v%mt(:,0,n,2))/2
163 : DO l = 0,3
164 : IF( enpara%qn_el(l,n,jsp).ne.0) THEN
165 : l_doneLocal(l,n,jsp) = .TRUE.
166 : el0Local(l,n,jsp)=find_enpara(.FALSE.,l,n,jsp,enpara%qn_el(l,n,jsp),atoms,v_mt,e_lo_local(l,n),e_up_local(l,n))
167 : IF( l .EQ. 3 ) THEN
168 : el0Local(4:,n,jsp) = el0Local(3,n,jsp)
169 : l_doneLocal(4:,n,jsp) = .TRUE.
170 : END IF
171 : ELSE
172 : l_doneLocal(l,n,jsp) = .FALSE.
173 : el0Local(l,n,jsp) = enpara%el0(l,n,jsp)
174 : IF(l .EQ. 3) THEN
175 : el0Local(4:,n,jsp) = enpara%el0(4:,n,jsp)
176 : END IF
177 : END IF
178 : ENDDO ! l
179 : ! Now for the lo's
180 : DO ilo = 1, atoms%nlo(n)
181 : l = atoms%llo(ilo,n)
182 : IF( enpara%qn_ello(ilo,n,jsp).NE.0) THEN
183 : lo_doneLocal(ilo,n,jsp) = .TRUE.
184 : ello0Local(ilo,n,jsp)=find_enpara(.TRUE.,l,n,jsp,enpara%qn_ello(ilo,n,jsp),atoms,v_mt,elo_lo_local(ilo,n),elo_up_local(ilo,n))
185 : ELSE
186 : ello0Local(ilo,n,jsp) = enpara%ello0(ilo,n,jsp)
187 : lo_doneLocal(ilo,n,jsp) = .FALSE.
188 : ENDIF
189 : ENDDO
190 : ENDDO ! n
191 : !$OMP END PARALLEL DO
192 1082 : CALL mpi_sum_reduce(el0Local(0:,:,jsp),enpara%el0(0:,:,jsp),fmpi%mpi_comm)
193 1082 : CALL mpi_sum_reduce(ello0Local(:,:,jsp),enpara%ello0(:,:,jsp),fmpi%mpi_comm)
194 1082 : CALL mpi_sum_reduce(e_lo_local,e_lo,fmpi%mpi_comm)
195 1082 : CALL mpi_sum_reduce(e_up_local,e_up,fmpi%mpi_comm)
196 1082 : CALL mpi_sum_reduce(elo_lo_local,elo_lo,fmpi%mpi_comm)
197 1082 : CALL mpi_sum_reduce(elo_up_local,elo_up,fmpi%mpi_comm)
198 1082 : CALL mpi_lor_reduce(l_doneLocal(:,:,jsp),l_done(:,:,jsp),fmpi%mpi_comm)
199 1082 : CALL mpi_lor_reduce(lo_doneLocal(:,:,jsp),lo_done(:,:,jsp),fmpi%mpi_comm)
200 1082 : CALL mpi_bc(enpara%el0,0,fmpi%mpi_comm)
201 1082 : CALL mpi_bc(enpara%ello0,0,fmpi%mpi_comm)
202 : #ifdef CPP_MPI
203 3246 : CALL MPI_BCAST(e_lo,SIZE(e_lo),MPI_DOUBLE_PRECISION,0,fmpi%mpi_comm,ierr)
204 3246 : CALL MPI_BCAST(e_up,SIZE(e_up),MPI_DOUBLE_PRECISION,0,fmpi%mpi_comm,ierr)
205 3246 : CALL MPI_BCAST(elo_lo,SIZE(elo_lo),MPI_DOUBLE_PRECISION,0,fmpi%mpi_comm,ierr)
206 3246 : CALL MPI_BCAST(elo_up,SIZE(elo_up),MPI_DOUBLE_PRECISION,0,fmpi%mpi_comm,ierr)
207 4328 : CALL MPI_BCAST(l_done,SIZE(l_done),MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
208 4328 : CALL MPI_BCAST(lo_done,SIZE(lo_done),MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
209 : #endif
210 :
211 1082 : IF (fmpi%irank==0) THEN
212 541 : WRITE(oUnit,*)
213 541 : WRITE(oUnit,*) "Updated energy parameters for spin:",jsp
214 : !Same loop for IO
215 1478 : DO n = 1, atoms%ntype
216 4685 : DO l = 0,3
217 4685 : IF( l_done(l,n,jsp)) CALL priv_write(.FALSE.,l,n,jsp,enpara%qn_el(l,n,jsp),e_lo(l,n),e_up(l,n),enpara%el0(l,n,jsp))
218 : ENDDO ! l
219 : ! Now for the lo's
220 2327 : DO ilo = 1, atoms%nlo(n)
221 849 : l = atoms%llo(ilo,n)
222 1786 : IF( lo_done(ilo,n,jsp)) CALL priv_write(.TRUE.,l,n,jsp,enpara%qn_ello(ilo,n,jsp),elo_lo(ilo,n),elo_up(ilo,n),enpara%ello0(ilo,n,jsp))
223 : END DO
224 : END DO
225 : ENDIF
226 :
227 : !! Now check for floating energy parameters (not for those with l_done=T)
228 1082 : IF (enpara%floating) THEN
229 0 : types_loop: DO n = 1,atoms%ntype
230 : !
231 : !---> determine energy parameters if lepr=1. the reference energy
232 : !---> is the value of the l=0 potential at approximately rmt/4.
233 : !
234 0 : j = atoms%jri(n) - (LOG(4.0)/atoms%dx(n)+1.51)
235 0 : rj = atoms%rmt(n)*EXP(atoms%dx(n)* (j-atoms%jri(n)))
236 0 : vbar = v%mt(j,0,n,jsp)/rj
237 0 : IF (fmpi%irank.EQ.0) THEN
238 0 : WRITE (oUnit,'('' spin'',i2,'', atom type'',i3,'' ='',f12.6,'' r='',f8.5)') jsp,n,vbar,rj
239 : ENDIF
240 0 : DO l = 0,atoms%lmax(n)
241 0 : IF ( .NOT.l_done(l,n,jsp) ) THEN
242 0 : enpara%el0(l,n,jsp) = vbar + enpara%el0(l,n,jsp)
243 : END IF
244 : ENDDO
245 0 : IF (atoms%nlo(n).GE.1) THEN
246 0 : DO ilo = 1,atoms%nlo(n)
247 0 : IF ( .NOT. lo_done(ilo,n,jsp) ) THEN
248 0 : enpara%ello0(ilo,n,jsp) = vbar + enpara%ello0(ilo,n,jsp)
249 : !+apw+lo
250 0 : IF (atoms%l_dulo(ilo,n)) THEN
251 0 : enpara%ello0(ilo,n,jsp) = enpara%el0(atoms%llo(ilo,n),n,jsp)
252 : ENDIF
253 : !-apw+lo
254 : END IF
255 : END DO
256 : ENDIF
257 : END DO types_loop
258 : ENDIF
259 1768 : IF (input%film) THEN
260 :
261 116 : INQUIRE (file ='enpara',exist= l_enpara)
262 132 : IF(l_enpara) enpara%evac0(:,jsp) = enpara%evac(:,jsp)
263 :
264 : !---> vacuum energy parameters: for floating: relative to potential
265 : !---> at vacuum-interstitial interface (better for electric field)
266 :
267 272 : DO ivac = 1,vacuum%nvac
268 156 : vz0 = 0.0
269 156 : IF (enpara%floating) THEN
270 0 : vz0 = REAL(v%vac(1,1,ivac,jsp))
271 0 : IF (fmpi%irank.EQ.0) THEN
272 0 : WRITE (oUnit,'('' spin'',i2,'', vacuum '',i3,'' ='',f12.6)') jsp,ivac,vz0
273 : ENDIF
274 : ENDIF
275 156 : enpara%evac(ivac,jsp) = enpara%evac0(ivac,jsp) + vz0
276 156 : IF (.NOT.l_enpara) THEN
277 148 : enpara%evac(ivac,jsp) = REAL(v%vac(vacuum%nmz,1,ivac,jsp)) + enpara%evac0(ivac,jsp)
278 : END IF
279 272 : IF (fmpi%irank.EQ.0) THEN
280 468 : attributes = ''
281 78 : WRITE(attributes(1),'(i0)') ivac
282 78 : WRITE(attributes(2),'(i0)') jsp
283 78 : WRITE(attributes(3),'(f16.10)') REAL(v%vac(1,1,ivac,jsp))
284 78 : WRITE(attributes(4),'(f16.10)') REAL(v%vac(vacuum%nmz,1,ivac,jsp))
285 78 : WRITE(attributes(5),'(f16.10)') enpara%evac(ivac,jsp)
286 : CALL writeXMLElementForm('vacuumEP',(/'vacuum','spin ','vzIR ','vzInf ','value '/),&
287 468 : attributes(1:5),RESHAPE((/6+4,4,4,5,5+13,8,1,16,16,16/),(/5,2/)))
288 : END IF
289 : ENDDO
290 116 : IF (vacuum%nvac.EQ.1) THEN
291 76 : enpara%evac(2,jsp) = enpara%evac(1,jsp)
292 : END IF
293 : END IF
294 : END DO
295 :
296 686 : IF(atoms%n_hia.GT.0 .AND. input%jspins.EQ.2 .AND. hub1data%l_performSpinavg) THEN
297 : !Set the energy parameters to the same value
298 : !We want the shell where Hubbard 1 is applied to
299 : !be non spin-polarized
300 0 : IF(fmpi%irank.EQ.0) THEN
301 0 : WRITE(oUnit,FMT=*)
302 0 : WRITE(oUnit,"(A)") "For Hubbard 1 we treat the correlated shell to be non spin-polarized"
303 : ENDIF
304 0 : DO j = atoms%n_u+1, atoms%n_u+atoms%n_hia
305 0 : l = atoms%lda_u(j)%l
306 0 : n = atoms%lda_u(j)%atomType
307 0 : enpara%el0(l,n,1) = (enpara%el0(l,n,1)+enpara%el0(l,n,2))/2.0
308 0 : enpara%el0(l,n,2) = enpara%el0(l,n,1)
309 : !-------------------------------------------
310 : ! Modify the higher energyParameters for f-orbitals
311 : !-------------------------------------------
312 0 : IF(l.EQ.3) THEN
313 0 : enpara%el0(4:,n,1) = enpara%el0(3,n,1)
314 0 : enpara%el0(4:,n,2) = enpara%el0(3,n,1)
315 : ENDIF
316 0 : IF(fmpi%irank.EQ.0) WRITE(oUnit,"(A,I3,A,I1,A,f16.10)") "New energy parameter atom ", n, " l ", l, "--> ", enpara%el0(l,n,1)
317 : ENDDO
318 : ENDIF
319 :
320 686 : IF(input%ldauAdjEnpara.AND.atoms%n_u+atoms%n_hia>0) THEN
321 : !If requested we can adjust the energy parameters to the LDA+U potential correction
322 0 : IF(irank.EQ.0) THEN
323 0 : WRITE(oUnit,FMT=*)
324 0 : WRITE(oUnit,"(A)") "LDA+U corrections for the energy parameters"
325 : ENDIF
326 0 : DO j = 1, atoms%n_u+atoms%n_hia
327 0 : l = atoms%lda_u(j)%l
328 0 : n = atoms%lda_u(j)%atomType
329 : !Calculate the trace of the LDA+U potential
330 0 : DO jsp = 1, input%jspins
331 0 : tr = 0.0
332 0 : DO i = -l,l
333 0 : tr = tr + REAL(v%mmpMat(i,i,j,jsp))
334 : ENDDO
335 0 : enpara%el0(l,n,jsp) = enpara%el0(l,n,jsp) + tr/REAL(2*l+1)
336 0 : IF(fmpi%irank.EQ.0) WRITE(oUnit,"(A,I3,A,I1,A,I1,A,f16.10)")&
337 0 : "New energy parameter atom ", n, " l ", l, " spin ", jsp,"--> ", enpara%el0(l,n,jsp)
338 : ENDDO
339 : ENDDO
340 : ENDIF
341 :
342 : ! enpara%ready=(ALL(enpara%el0>-1E99).AND.ALL(enpara%ello0>-1E99))
343 26326 : enpara%epara_min=MIN(MINVAL(enpara%el0),MINVAL(enpara%ello0))
344 :
345 686 : IF (fmpi%irank == 0) CALL closeXMLElement('energyParameters')
346 686 : END SUBROUTINE update
347 :
348 16 : SUBROUTINE READ(enpara,atoms,jspins,film,l_required)
349 : USE m_types_setup
350 : USE m_constants
351 : IMPLICIT NONE
352 : CLASS(t_enpara),INTENT(INOUT):: enpara
353 : INTEGER, INTENT (IN) :: jspins
354 : TYPE(t_atoms),INTENT(IN) :: atoms
355 : LOGICAL,INTENT(IN) :: film,l_required
356 :
357 : INTEGER :: n,l,lo,skip_t,io_err,jsp
358 : logical :: l_exist
359 :
360 16 : INQUIRE(file="enpara",exist=l_exist)
361 16 : IF (.NOT.l_exist.AND.l_required) CALL juDFT_error("No enpara file found")
362 32 : IF (.NOT.l_exist) RETURN
363 :
364 16 : OPEN(40,file="enpara",form="formatted",status="old")
365 :
366 38 : DO jsp=1,jspins
367 :
368 : !--> first line contains mixing parameter!
369 :
370 22 : READ (40,FMT ='(48x,f10.6)',END=200) enpara%enmix(jsp)
371 22 : READ (40,*) ! skip next line
372 22 : IF (enpara%enmix(jsp).EQ.0.0) enpara%enmix(jsp) = 1.0
373 22 : WRITE (oUnit,FMT=8001) jsp
374 22 : WRITE (oUnit,FMT=8000)
375 22 : skip_t = 0
376 62 : DO n = 1,atoms%ntype
377 200 : READ (40,FMT=8040,END=200) (enpara%el0(l,n,jsp),l=0,3),&
378 400 : (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp)
379 240 : WRITE (oUnit,FMT=8140) n,(enpara%el0(l,n,jsp),l=0,3),&
380 440 : (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp)
381 : !
382 : !---> energy parameters for the local orbitals
383 : !
384 40 : IF (atoms%nlo(n).GE.1) THEN
385 4 : skip_t = skip_t + enpara%skiplo(n,jsp) * atoms%neq(n)
386 8 : READ (40,FMT=8039,END=200) (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n))
387 8 : READ (40,FMT=8038,END=200) (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n))
388 8 : WRITE (oUnit,FMT=8139) (enpara%ello0(lo,n,jsp),lo=1,atoms%nlo(n))
389 8 : WRITE (oUnit,FMT=8138) (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n))
390 36 : ELSEIF (enpara%skiplo(n,jsp).GT.0) THEN
391 0 : WRITE (oUnit,*) "for atom",n," no LO's were specified"
392 0 : WRITE (oUnit,*) 'but skiplo was set to',enpara%skiplo
393 : CALL juDFT_error("No LO's but skiplo",calledby ="enpara",&
394 0 : hint="If no LO's are set skiplo must be 0 in enpara")
395 : END IF
396 : !Integer values mean we have to set the qn-arrays
397 200 : enpara%qn_el(:,n,jsp)=0
398 200 : DO l=0,3
399 200 : IF (enpara%el0(l,n,jsp)==NINT(enpara%el0(l,n,jsp))) enpara%qn_el(l,n,jsp)=NINT(enpara%el0(l,n,jsp))
400 : ENDDO
401 48 : enpara%qn_ello(:,n,jsp)=0
402 44 : DO l=1,atoms%nlo(n)
403 44 : IF (enpara%ello0(l,n,jsp)==NINT(enpara%ello0(l,n,jsp))) enpara%qn_ello(l,n,jsp)=NINT(enpara%ello0(l,n,jsp))
404 : ENDDO
405 : !
406 : !---> set the energy parameters with l>3 to the value of l=3
407 : !
408 374 : enpara%el0(4:,n,jsp) = enpara%el0(3,n,jsp)
409 : ENDDO ! atoms%ntype
410 :
411 38 : IF (film) THEN
412 56 : enpara%lchg_v = .TRUE.
413 8 : READ (40,FMT=8050,END=200) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp)
414 8 : WRITE (oUnit,FMT=8150) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp)
415 : ENDIF
416 : ! IF (atoms%nlod.GE.1) THEN
417 : ! WRITE (oUnit,FMT=8090) jsp,skip_t
418 : ! WRITE (oUnit,FMT=8091)
419 : ! END IF
420 : END DO
421 :
422 112 : enpara%evac(:,:) = enpara%evac0(:,:)
423 :
424 16 : CLOSE(40)
425 : ! input formats
426 :
427 : 8038 FORMAT (14x,60(l1,8x))
428 : 8039 FORMAT (8x,60f9.5)
429 : 8040 FORMAT (8x,4f9.5,9x,4l1,9x,i3)
430 : 8050 FORMAT (19x,f9.5,9x,l1,15x,f9.5)
431 :
432 : ! output formats
433 :
434 : 8138 FORMAT (' --> change ',60(l1,8x))
435 : 8139 FORMAT (' --> lo ',60f9.5)
436 : 8140 FORMAT (' -->',i3,1x,4f9.5,' change: ',4l1,' skiplo: ',i3)
437 : 8150 FORMAT (' vacuum parameter=',f9.5,' change: ',l1, ' second vacuum=',f9.5)
438 : 8001 FORMAT ('READING enpara for spin: ',i1)
439 : 8000 FORMAT (/,' energy parameters:',/,t10,'s',t20, 'p',t30,'d',t37,'higher l - - -')
440 : 8090 FORMAT ('Spin: ',i1,' -- ',i3,'eigenvalues')
441 : 8091 FORMAT ('will be skipped for energyparameter computation')
442 :
443 16 : RETURN
444 :
445 0 : 200 WRITE (oUnit,*) 'the end of the file enpara has been reached while'
446 0 : WRITE (oUnit,*) 'reading the energy-parameters.'
447 0 : WRITE (oUnit,*) 'possible reason: energy parameters have not been'
448 0 : WRITE (oUnit,*) 'specified for all atom types.'
449 0 : WRITE (oUnit,FMT='(a,i4)') 'the actual number of atom-types is: ntype=',atoms%ntype
450 0 : CALL juDFT_error ("unexpected end of file enpara reached while reading")
451 : END SUBROUTINE read
452 :
453 :
454 24 : SUBROUTINE WRITE(enpara, atoms,jspins,film)
455 :
456 : ! write enpara-file
457 : !
458 : USE m_types_setup
459 : USE m_constants
460 : IMPLICIT NONE
461 : CLASS(t_enpara),INTENT(IN) :: enpara
462 : INTEGER, INTENT (IN) :: jspins
463 : LOGICAL,INTENT(IN) :: film
464 : TYPE(t_atoms),INTENT(IN) :: atoms
465 :
466 : INTEGER n,l,lo,jspin
467 24 : REAL el0Temp(0:3), ello0Temp(1:atoms%nlod)
468 :
469 24 : OPEN(unit=40,file="enpara",form="formatted",status="replace")
470 :
471 70 : DO jspin=1,jspins
472 46 : WRITE (40,FMT=8035) jspin,enpara%enmix(jspin)
473 46 : WRITE (40,FMT=8036)
474 : 8035 FORMAT (5x,'energy parameters for spin ',i1,' mix=',f10.6)
475 : 8036 FORMAT (t6,'atom',t15,'s',t24,'p',t33,'d',t42,'f')
476 136 : DO n = 1,atoms%ntype
477 450 : el0Temp(0:3) = enpara%el0(0:3,n,jspin)
478 450 : DO l = 0, 3
479 450 : IF (enpara%qn_el(l,n,jspin).NE.0) el0Temp(l) = REAL(enpara%qn_el(l,n,jspin))
480 : END DO
481 90 : WRITE (oUnit,FMT=8040) n, (el0Temp(l),l=0,3),&
482 540 : & (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin)
483 90 : WRITE (40,FMT=8040) n, (el0Temp(l),l=0,3),&
484 540 : & (enpara%lchange(l,n,jspin),l=0,3),enpara%skiplo(n,jspin)
485 : !---> energy parameters for the local orbitals
486 136 : IF (atoms%nlo(n).GE.1) THEN
487 4 : ello0Temp(1:atoms%nlo(n)) = enpara%ello0(1:atoms%nlo(n),n,jspin)
488 4 : DO lo = 1, atoms%nlo(n)
489 4 : IF (enpara%qn_ello(lo,n,jspin).NE.0) ello0Temp(lo) = enpara%qn_ello(lo,n,jspin)
490 : END DO
491 2 : WRITE (oUnit,FMT=8039) (ello0Temp(lo),lo=1,atoms%nlo(n))
492 4 : WRITE (oUnit,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n))
493 2 : WRITE (40,FMT=8039) (ello0Temp(lo),lo=1,atoms%nlo(n))
494 4 : WRITE (40,FMT=8038) (enpara%llochg(lo,n,jspin),lo=1,atoms%nlo(n))
495 : END IF
496 :
497 : ENDDO
498 : 8038 FORMAT (' --> change ',60(l1,8x))
499 : 8039 FORMAT (' --> lo ',60f9.5)
500 : 8040 FORMAT (' -->',i3,1x,4f9.5,' change: ',4l1,' skiplo: ',i3)
501 :
502 70 : IF (film) THEN
503 4 : WRITE (40,FMT=8050) enpara%evac(1,jspin),enpara%lchg_v(1,jspin),enpara%evac(2,jspin)
504 4 : WRITE (oUnit,FMT=8050) enpara%evac(1,jspin),enpara%lchg_v(1,jspin),enpara%evac(2,jspin)
505 : 8050 FORMAT (' vacuum parameter=',f9.5,' change: ',l1,&
506 : & ' second vacuum=',f9.5)
507 : ENDIF
508 : ENDDO
509 24 : CLOSE(40)
510 24 : RETURN
511 : END SUBROUTINE WRITE
512 :
513 :
514 :
515 :
516 :
517 664 : SUBROUTINE mix(enpara,fmpi_comm,atoms,vacuum,input,pot)
518 : !------------------------------------------------------------------
519 : USE m_types_atoms
520 : USE m_types_input
521 : USE m_types_vacuum
522 : USE m_types_potden
523 : USE m_constants
524 : #ifdef CPP_MPI
525 : USE mpi
526 : #endif
527 : IMPLICIT NONE
528 : CLASS(t_enpara),INTENT(INOUT) :: enpara
529 : INTEGER,INTENT(IN) :: fmpi_comm
530 : TYPE(t_atoms),INTENT(IN) :: atoms
531 : TYPE(t_vacuum),INTENT(IN) :: vacuum
532 : TYPE(t_input),INTENT(IN) :: input
533 : TYPE(t_potden),INTENT(IN) :: pot
534 :
535 : !REAL, INTENT(IN) :: vr(:,:,:)
536 : !REAL, INTENT(IN) :: vz(vacuum%nmzd,2)
537 :
538 : INTEGER ityp,j,l,lo,jsp,n
539 : REAL vbar,maxdist,maxdist2
540 664 : INTEGER same(atoms%nlod),irank
541 : LOGICAL l_enpara
542 : #ifdef CPP_MPI
543 : INTEGER :: ierr
544 664 : CALL MPI_COMM_RANK(fmpi_comm,irank,ierr)
545 : #else
546 : irank=0
547 : #endif
548 :
549 664 : IF (irank==0) THEN
550 332 : maxdist2=0.0
551 859 : DO jsp=1,SIZE(enpara%el0,3)
552 527 : maxdist=0.0
553 1431 : DO ityp = 1,atoms%ntype
554 : ! look for LO's energy parameters equal to the LAPW (and previous LO) ones
555 1863 : same = 0
556 1733 : DO lo = 1,atoms%nlo(ityp)
557 829 : IF(enpara%el0(atoms%llo(lo,ityp),ityp,jsp).EQ.enpara%ello0(lo,ityp,jsp)) same(lo)=-1
558 2067 : DO l = 1,lo-1
559 334 : IF(atoms%llo(l,ityp).NE.atoms%llo(lo,ityp)) CYCLE
560 831 : IF(enpara%ello0(l,ityp,jsp).EQ.enpara%ello0(lo,ityp,jsp).AND.same(lo).EQ.0) same(lo)=l
561 : ENDDO
562 : ENDDO
563 : !
564 : !---> change energy parameters
565 : !
566 4520 : DO l = 0,3
567 3616 : WRITE(oUnit,*) 'Type:',ityp,' l:',l
568 3616 : WRITE(oUnit,FMT=777) enpara%el0(l,ityp,jsp),enpara%el1(l,ityp,jsp),&
569 7232 : ABS(enpara%el0(l,ityp,jsp)-enpara%el1(l,ityp,jsp))
570 3616 : maxdist=MAX(maxdist,ABS(enpara%el0(l,ityp,jsp)-enpara%el1(l,ityp,jsp)))
571 4520 : IF ( enpara%lchange(l,ityp,jsp) ) THEN
572 344 : maxdist2=MAX(maxdist2,ABS(enpara%el0(l,ityp,jsp)-enpara%el1(l,ityp,jsp)))
573 : enpara%el0(l,ityp,jsp) =(1.0-enpara%enmix(jsp))*enpara%el0(l,ityp,jsp) + &
574 344 : enpara%enmix(jsp)*enpara%el1(l,ityp,jsp)
575 : ENDIF
576 : ENDDO
577 1350 : IF ( enpara%lchange(3,ityp,jsp) ) enpara%el0(4:,ityp,jsp) = enpara%el0(3,ityp,jsp)
578 : !
579 : !---> determine and change local orbital energy parameters
580 : !
581 1733 : DO lo = 1,atoms%nlo(ityp)
582 1733 : IF (atoms%l_dulo(lo,ityp)) THEN
583 0 : enpara%ello0(lo,ityp,jsp) =enpara%el0(atoms%llo(lo,ityp),ityp,jsp)
584 : ELSE
585 829 : IF (enpara%llochg(lo,ityp,jsp) ) THEN
586 0 : IF(same(lo).EQ.-1) THEN
587 0 : enpara%ello0(lo,ityp,jsp) = enpara%el0(atoms%llo(lo,ityp),ityp,jsp)
588 0 : CYCLE
589 0 : ELSE IF(same(lo).GT.0) THEN
590 0 : enpara%ello0(lo,ityp,jsp) = enpara%ello0(same(lo),ityp,jsp)
591 0 : CYCLE
592 : ENDIF
593 : ENDIF
594 829 : WRITE(oUnit,*) 'Type:',ityp,' lo:',lo
595 829 : WRITE(oUnit,FMT=777) enpara%ello0(lo,ityp,jsp),enpara%ello1(lo,ityp,jsp),&
596 1658 : ABS(enpara%ello0(lo,ityp,jsp)-enpara%ello1(lo,ityp,jsp))
597 829 : maxdist=MAX(maxdist,ABS(enpara%ello0(lo,ityp,jsp)-enpara%ello1(lo,ityp,jsp)))
598 829 : IF (enpara%llochg(lo,ityp,jsp) ) THEN
599 0 : maxdist2=MAX(maxdist2,ABS(enpara%ello0(lo,ityp,jsp)-enpara%ello1(lo,ityp,jsp)))
600 : enpara%ello0(lo,ityp,jsp) =(1.0-enpara%enmix(jsp))*enpara%ello0(lo,ityp,jsp)+&
601 0 : enpara%enmix(jsp)*enpara%ello1(lo,ityp,jsp)
602 : ENDIF
603 : END IF
604 : END DO
605 : !Shift if floating energy parameters are used
606 1431 : IF (enpara%floating) THEN
607 0 : j = atoms%jri(ityp) - (LOG(4.0)/atoms%dx(ityp)+1.51)
608 0 : vbar = pot%mt(j,0,ityp,jsp)/( atoms%rmt(ityp)*EXP(atoms%dx(ityp)*(j-atoms%jri(ityp))) )
609 0 : enpara%el0(:,n,jsp)=enpara%el0(:,n,jsp)-vbar
610 : ENDIF
611 : END DO
612 :
613 :
614 527 : IF (input%film) THEN
615 56 : WRITE(oUnit,*) 'Vacuum:'
616 132 : DO n=1,vacuum%nvac
617 76 : WRITE(oUnit,FMT=777) enpara%evac(n,jsp),enpara%evac1(n,jsp),ABS(enpara%evac(n,jsp)-enpara%evac1(n,jsp))
618 76 : maxdist=MAX(maxdist,ABS(enpara%evac(n,jsp)-enpara%evac1(n,jsp)))
619 132 : IF (enpara%lchg_v(n,jsp) ) THEN
620 4 : maxdist2=MAX(maxdist2,ABS(enpara%evac(n,jsp)-enpara%evac1(n,jsp)))
621 : enpara%evac(n,jsp) =(1.0-enpara%enmix(jsp))*enpara%evac(n,jsp)+&
622 4 : enpara%enmix(jsp)*enpara%evac1(n,jsp)
623 : END IF
624 : END DO
625 56 : IF (vacuum%nvac==1) enpara%evac(2,jsp) = enpara%evac(1,jsp)
626 56 : IF (enpara%floating) enpara%evac(:,jsp)=enpara%evac(:,jsp)-REAL(pot%vac(:,1,1,jsp)) ! This is broken!!!
627 : ENDIF
628 859 : WRITE(oUnit,'(a36,f12.6)') 'Max. mismatch of energy parameters:', maxdist
629 : END DO
630 332 : IF (maxdist2>1.0) CALL juDFT_warn&
631 : ("Energy parameter mismatch too large",hint&
632 : ="If any energy parameters calculated from the output "//&
633 : "differ from the input by more than 1Htr, chances are "//&
634 0 : "high that your initial setup was broken.")
635 : ENDIF
636 664 : INQUIRE(file='enpara',exist=l_enpara)
637 664 : IF (irank==0.AND.l_enpara) CALL enpara%WRITE(atoms,input%jspins,input%film)
638 : #ifdef CPP_MPI
639 2656 : CALL MPI_BCAST(enpara%el0,SIZE(enpara%el0),MPI_DOUBLE_PRECISION,0,fmpi_comm,ierr)
640 2656 : CALL MPI_BCAST(enpara%ello0,SIZE(enpara%ello0),MPI_DOUBLE_PRECISION,0,fmpi_comm,ierr)
641 664 : CALL MPI_BCAST(enpara%evac,SIZE(enpara%evac),MPI_DOUBLE_PRECISION,0,fmpi_comm,ierr)
642 : #endif
643 664 : RETURN
644 : 777 FORMAT('Old:',f8.5,' new:',f8.5,' diff:',f8.5)
645 :
646 : END SUBROUTINE mix
647 :
648 : #ifndef CPP_INPGEN
649 670 : SUBROUTINE calcOutParams(enpara,input,atoms,vacuum,regCharges)
650 : USE m_types_setup
651 : USE m_types_regionCharges
652 : IMPLICIT NONE
653 : CLASS(t_enpara),INTENT(INOUT) :: enpara
654 : TYPE(t_input),INTENT(IN) :: input
655 : TYPE(t_atoms),INTENT(IN) :: atoms
656 : TYPE(t_vacuum),INTENT(IN) :: vacuum
657 : TYPE(t_regionCharges),INTENT(IN) :: regCharges
658 :
659 : INTEGER :: ispin, n, ilo, iVac, l
660 :
661 20932 : enpara%el1(:,:,:) = enpara%el0(:,:,:)
662 5526 : enpara%ello1(:,:,:) = enpara%ello0(:,:,:)
663 4690 : enpara%evac1(:,:) = enpara%evac(:,:)
664 :
665 1732 : DO ispin = 1,input%jspins
666 2898 : DO n=1,atoms%ntype
667 9180 : DO l = 0, 3
668 9180 : IF (regCharges%sqal(l,n,ispin).NE.0.0) enpara%el1(l,n,ispin)=regCharges%ener(l,n,ispin)/regCharges%sqal(l,n,ispin)
669 : END DO
670 2898 : IF (atoms%nlo(n)>0) THEN
671 2708 : DO ilo = 1, atoms%nlo(n)
672 2708 : IF (regCharges%sqlo(ilo,n,ispin).NE.0.0) THEN
673 1682 : enpara%ello1(ilo,n,ispin)=regCharges%enerlo(ilo,n,ispin)/regCharges%sqlo(ilo,n,ispin)
674 : END IF
675 : END DO
676 : END IF
677 : END DO
678 1732 : IF (input%film) THEN
679 268 : DO iVac = 1, vacuum%nvac
680 268 : IF (regCharges%svac(iVac,ispin).NE.0.0) THEN
681 154 : enpara%evac1(iVac,ispin)=regCharges%pvac(iVac,ispin)/regCharges%svac(iVac,ispin)
682 : END IF
683 : END DO
684 : END IF
685 : END DO
686 670 : END SUBROUTINE calcOutParams
687 : #endif
688 4211 : SUBROUTINE priv_write(lo,l,n,jsp,nqn,e_lo,e_up,e)
689 : !subroutine to write energy parameters to output
690 : USE m_constants
691 : USE m_xmlOutput
692 : IMPLICIT NONE
693 : LOGICAL,INTENT(IN):: lo
694 : INTEGER,INTENT(IN):: l,n,jsp,nqn
695 : REAL,INTENT(IN) :: e_lo,e_up,e
696 :
697 : CHARACTER(LEN=20) :: attributes(6)
698 4211 : CHARACTER(len=:),ALLOCATABLE:: label
699 : CHARACTER(len=1),PARAMETER,DIMENSION(0:9):: ch=(/'s','p','d','f','g','h','i','j','k','l'/)
700 :
701 29477 : attributes = ''
702 4211 : WRITE(attributes(1),'(i0)') n
703 4211 : WRITE(attributes(2),'(i0)') jsp
704 4211 : WRITE(attributes(3),'(i0,a1)') abs(nqn), ch(l)
705 4211 : WRITE(attributes(4),'(f8.2)') e_lo
706 4211 : WRITE(attributes(5),'(f8.2)') e_up
707 4211 : WRITE(attributes(6),'(f16.10)') e
708 4211 : IF (nqn>0) THEN
709 4207 : IF (lo) THEN
710 843 : label='loAtomicEP'
711 : ELSE
712 3364 : label='atomicEP'
713 : ENDIF
714 : ELSE
715 4 : IF (lo) THEN
716 4 : label='heloAtomicEP'
717 : ELSE
718 0 : label='heAtomicEP'
719 : ENDIF
720 : END IF
721 :
722 : CALL writeXMLElementForm(label,(/'atomType ','spin ','branch ',&
723 : 'branchLowest ','branchHighest','value '/),&
724 29477 : attributes,RESHAPE((/10,4,6,12,13,5,6,1,3,8,8,16/),(/6,2/)))
725 4211 : IF (lo) THEN
726 847 : WRITE(oUnit,'(a6,i5,i2,a1,a12,f6.2,a3,f6.2,a13,f8.4)') ' Atom',n,nqn,ch(l),' branch from',&
727 1694 : e_lo, ' to',e_up,' htr. ; e_l(lo) =',e
728 : ELSE
729 3364 : WRITE(oUnit,'(a6,i5,i2,a1,a12,f6.2,a3,f6.2,a13,f8.4)') ' Atom',n,nqn,ch(l),' branch from',&
730 6728 : e_lo, ' to',e_up,' htr. ; e_l =',e
731 : END IF
732 4211 : END SUBROUTINE priv_write
733 :
734 0 : END MODULE m_types_enpara
|
