Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2021 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 : MODULE m_dfpt_sternheimer
7 : USE m_types
8 : USE m_make_stars
9 : USE m_vgen
10 : USE m_dfpt_eigen
11 : USE m_dfpt_cdngen
12 : USE m_dfpt_vgen
13 : USE m_dfpt_fermie
14 : USE m_mix
15 : USE m_constants
16 : USE m_cdn_io
17 : USE m_eig66_io
18 :
19 : IMPLICIT NONE
20 :
21 : CONTAINS
22 0 : SUBROUTINE dfpt_sternheimer(fi, xcpot, sphhar, stars, starsq, nococonv, qpts, fmpi, results, resultsq, enpara, hybdat, &
23 : rho, vTot, grRho, grVtot, grVext, iQ, iDType, iDir, dfpt_tag, eig_id, &
24 : l_real, results1, dfpt_eig_id, dfpt_eig_id2, q_eig_id, &
25 : denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im, sigma_disc, sigma_disc2, &
26 : starsmq, resultsmq, dfpt_eigm_id, dfpt_eigm_id2, qm_eig_id, results1m, vTot1m, vTot1mIm)
27 : TYPE(t_fleurinput), INTENT(IN) :: fi
28 : CLASS(t_xcpot), INTENT(IN) :: xcpot
29 : TYPE(t_sphhar), INTENT(IN) :: sphhar
30 : TYPE(t_stars), INTENT(IN) :: stars
31 : TYPE(t_stars), INTENT(INOUT) :: starsq
32 : TYPE(t_nococonv), INTENT(IN) :: nococonv
33 : TYPE(t_kpts), INTENT(IN) :: qpts
34 : TYPE(t_mpi), INTENT(IN) :: fmpi
35 : TYPE(t_results), INTENT(INOUT) :: results, resultsq, results1
36 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
37 : TYPE(t_enpara), INTENT(INOUT) :: enpara
38 : TYPE(t_potden), INTENT(IN) :: rho, vTot, grRho, grVtot, grVext
39 :
40 : TYPE(t_potden), INTENT(INOUT) :: denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im
41 :
42 : INTEGER, INTENT(IN) :: iQ, iDtype, iDir, eig_id, q_eig_id
43 : LOGICAL, INTENT(IN) :: l_real
44 : CHARACTER(len=20), INTENT(IN) :: dfpt_tag
45 :
46 : INTEGER, INTENT(IN) :: dfpt_eig_id, dfpt_eig_id2
47 :
48 : complex, intent(in) :: sigma_disc(2), sigma_disc2(2)
49 :
50 : TYPE(t_stars), OPTIONAL, INTENT(INOUT) :: starsmq
51 : TYPE(t_results), OPTIONAL, INTENT(INOUT) :: resultsmq, results1m
52 : INTEGER, OPTIONAL, INTENT(IN) :: qm_eig_id
53 : INTEGER, OPTIONAL, INTENT(IN) :: dfpt_eigm_id, dfpt_eigm_id2
54 : TYPE(t_potden), OPTIONAL, INTENT(INOUT) :: vTot1m, vTot1mIm
55 :
56 : #ifdef CPP_MPI
57 : INTEGER :: ierr
58 : #endif
59 :
60 : INTEGER :: archiveType, iter, killcont(6), iterm, realiter
61 : REAL :: bqpt(3), bmqpt(3)
62 : LOGICAL :: l_cont, l_exist, l_lastIter, l_dummy, strho, onedone, final_SH_it, l_minusq, l_existm
63 :
64 : complex :: sigma_loc(2)
65 :
66 0 : TYPE(t_banddos) :: banddosdummy
67 0 : TYPE(t_field) :: field2
68 0 : TYPE(t_potden) :: denOut1, denOut1Im, rho_loc, rho_loc0
69 0 : TYPE(t_potden) :: denIn1m, denIn1mIm, denOut1m, denOut1mIm
70 :
71 0 : l_minusq = PRESENT(starsmq)
72 0 : realiter = 0
73 :
74 : ! killcont can be used to blot out certain contricutions to the
75 : ! perturbed matrices.
76 : ! In this order: V1_pw_pw, T1_pw, S1_pw, V1_MT, ikGH0_MT, ikGS0_MT
77 0 : killcont = [1,1,1,1,1,1]
78 :
79 0 : CALL rho_loc%copyPotDen(rho)
80 0 : CALL rho_loc0%copyPotDen(rho)
81 0 : CALL rho_loc0%resetPotDen()
82 :
83 0 : banddosdummy = fi%banddos
84 :
85 : ! Calculate the q-shifted G-vectors and related quantities like the perturbed step function
86 0 : CALL make_stars(starsq, fi%sym, fi%atoms, fi%vacuum, sphhar, fi%input, fi%cell, fi%noco, fmpi, qpts%bk(:,iQ), iDtype, iDir)!TODO: Theta1 raus fuer efield
87 0 : starsq%ufft = stars%ufft
88 :
89 : ! Initialize the density perturbation; denIn1Im is only for the imaginary MT part
90 0 : CALL denIn1%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN, l_dfpt=.TRUE.)
91 0 : CALL denIn1Im%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN, l_dfpt=.FALSE.)
92 0 : INQUIRE(FILE=TRIM(dfpt_tag)//'.hdf',EXIST=l_exist)
93 :
94 0 : IF (l_minusq) THEN
95 0 : CALL make_stars(starsmq, fi%sym, fi%atoms, fi%vacuum, sphhar, fi%input, fi%cell, fi%noco, fmpi, -qpts%bk(:,iQ), iDtype, iDir)
96 0 : starsmq%ufft = stars%ufft
97 :
98 0 : CALL denIn1m%init(starsmq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN, l_dfpt=.TRUE.)
99 0 : CALL denIn1mIm%init(starsmq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN, l_dfpt=.FALSE.)
100 0 : INQUIRE(FILE=TRIM(dfpt_tag)//'m.hdf',EXIST=l_existm)
101 : END IF
102 :
103 0 : archiveType = CDN_ARCHIVE_TYPE_CDN1_const
104 0 : IF (ANY(fi%noco%l_unrestrictMT)) THEN
105 0 : archiveType = CDN_ARCHIVE_TYPE_FFN_const
106 0 : ELSE IF (fi%noco%l_noco) THEN
107 0 : archiveType = CDN_ARCHIVE_TYPE_NOCO_const
108 : END IF
109 :
110 0 : IF (fmpi%irank == 0) THEN
111 0 : strho = .NOT.l_exist ! There is no density perturbation file yet --> starting density perturbation
112 0 : onedone = .NOT.strho ! Was at least one iteration done yet?
113 0 : final_SH_it = .FALSE. ! Is the density perturbation converged and the last SH run started?
114 : END IF
115 :
116 : #ifdef CPP_MPI
117 0 : CALL MPI_BCAST(strho,1,MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
118 0 : CALL MPI_BCAST(onedone,1,MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
119 0 : CALL MPI_BCAST(final_SH_it,1,MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
120 : #endif
121 :
122 : ! Set the iteration counters to 0 and optionally read in the density perturbation
123 0 : iter = 0
124 0 : iterm = 0
125 0 : l_cont = (iter < fi%input%itmax)
126 :
127 0 : IF (fmpi%irank==0.AND.l_exist) CALL readDensity(starsq, fi%noco, fi%vacuum, fi%atoms, fi%cell, sphhar, &
128 : fi%input, fi%sym, archiveType, CDN_INPUT_DEN_const, 0, &
129 : results%ef, results%last_distance, l_dummy, denIn1, &
130 0 : inFilename=TRIM(dfpt_tag),denIm=denIn1Im)
131 0 : IF (fmpi%irank==0.AND.l_exist.AND.l_minusq) CALL readDensity(starsmq, fi%noco, fi%vacuum, fi%atoms, fi%cell, sphhar, &
132 : fi%input, fi%sym, archiveType, CDN_INPUT_DEN_const, 0, &
133 : results%ef, results%last_distance, l_dummy, denIn1m, &
134 0 : inFilename=TRIM(dfpt_tag)//'m',denIm=denIn1mIm)
135 :
136 0 : IF (fmpi%irank==0.AND..NOT.l_exist) denIn1%iter = 1
137 :
138 : #ifdef CPP_MPI
139 0 : CALL MPI_BCAST(denIn1%iter,1,MPI_INTEGER,0,fmpi%mpi_comm,ierr)
140 : #endif
141 :
142 : ! Initialize the potentials and save the q vector to a local variable
143 0 : CALL vTot1%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.TRUE.)
144 0 : CALL vTot1Im%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.FALSE.)
145 :
146 0 : CALL vC1%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.TRUE.)
147 0 : CALL vC1Im%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.FALSE.)
148 :
149 0 : bqpt = qpts%bk(:, iQ)
150 :
151 0 : IF (l_minusq) THEN
152 0 : CALL vTot1m%init(starsmq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.TRUE.)
153 0 : CALL vTot1mIm%init(starsmq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.FALSE.)
154 :
155 0 : bmqpt = -qpts%bk(:, iQ)
156 : END IF
157 :
158 : l_lastIter = .FALSE.
159 0 : scfloop: DO WHILE (l_cont)
160 0 : IF (onedone) iter = iter + 1
161 : l_lastIter = l_lastIter.OR.(iter.EQ.fi%input%itmax)
162 :
163 0 : CALL timestart("Sternheimer Iteration")
164 0 : IF (fmpi%irank==0.AND.onedone) THEN
165 0 : WRITE (oUnit, FMT=8100) iter
166 : 8100 FORMAT(/, 10x, ' iter= ', i5)
167 : END IF !fmpi%irank==0
168 :
169 0 : CALL denIn1%distribute(fmpi%mpi_comm)
170 0 : CALL denIn1Im%distribute(fmpi%mpi_comm)
171 :
172 0 : IF (l_minusq) THEN
173 0 : CALL denIn1m%distribute(fmpi%mpi_comm)
174 0 : CALL denIn1mIm%distribute(fmpi%mpi_comm)
175 : END IF
176 :
177 : ! Generate the potential perturbation:
178 : ! Vext1 for the starting perturbation
179 : ! Veff1 every other time
180 0 : CALL timestart("Generation of potential perturbation")
181 0 : IF (strho) THEN
182 0 : write(oUnit, *) "vExt1", iDir
183 0 : sigma_loc = cmplx(0.0,0.0)
184 : !IF (iDir==3) sigma_loc = -sigma_disc
185 : CALL dfpt_vgen(hybdat,fi%field,fi%input,xcpot,fi%atoms,sphhar,stars,fi%vacuum,fi%sym,&
186 : fi%juphon,fi%cell,fmpi,fi%noco,nococonv,rho_loc0,vTot,&
187 0 : starsq,denIn1Im,vTot1,.FALSE.,vTot1Im,denIn1,iDtype,iDir,[1,1],sigma_loc)!-?
188 : ! Last variable: [m,n] dictates with [1/0, 1/0], whether or not we take
189 : ! V1*Theta and V*Theta1 into account respectively. For [1,1] all is
190 : ! contained.
191 0 : IF (l_minusq) THEN
192 0 : sigma_loc = cmplx(0.0,0.0)
193 : !IF (iDir==3) sigma_loc = -sigma_disc
194 : CALL dfpt_vgen(hybdat,fi%field,fi%input,xcpot,fi%atoms,sphhar,stars,fi%vacuum,fi%sym,&
195 : fi%juphon,fi%cell,fmpi,fi%noco,nococonv,rho_loc0,vTot,&
196 0 : starsmq,denIn1mIm,vTot1m,.FALSE.,vTot1mIm,denIn1m,iDtype,iDir,[1,1],sigma_loc)!-?
197 : END IF
198 : ELSE
199 0 : write(oUnit, *) "vEff1", iDir
200 0 : sigma_loc = cmplx(0.0,0.0)
201 : !IF (iDir==3) sigma_loc = -sigma_disc2
202 : CALL dfpt_vgen(hybdat,fi%field,fi%input,xcpot,fi%atoms,sphhar,stars,fi%vacuum,fi%sym,&
203 : fi%juphon,fi%cell,fmpi,fi%noco,nococonv,rho_loc,vTot,&
204 0 : starsq,denIn1Im,vTot1,.TRUE.,vTot1Im,denIn1,iDtype,iDir,[1,1],sigma_loc)
205 0 : IF (l_minusq) THEN
206 0 : sigma_loc = cmplx(0.0,0.0)
207 : !IF (iDir==3) sigma_loc = -sigma_disc2
208 : CALL dfpt_vgen(hybdat,fi%field,fi%input,xcpot,fi%atoms,sphhar,stars,fi%vacuum,fi%sym,&
209 : fi%juphon,fi%cell,fmpi,fi%noco,nococonv,rho_loc,vTot,&
210 0 : starsmq,denIn1mIm,vTot1m,.TRUE.,vTot1mIm,denIn1m,iDtype,iDir,[1,1],sigma_loc)
211 : END IF
212 : END IF
213 :
214 : ! For the calculation of the dynamical matrix, we need VC1 additionally
215 0 : IF (final_SH_it) THEN
216 0 : write(oUnit, *) "vC1", iDir
217 0 : sigma_loc = cmplx(0.0,0.0)
218 : !IF (iDir==3) sigma_loc = -sigma_disc2
219 : CALL dfpt_vgen(hybdat,fi%field,fi%input,xcpot,fi%atoms,sphhar,stars,fi%vacuum,fi%sym,&
220 : fi%juphon,fi%cell,fmpi,fi%noco,nococonv,rho_loc,vTot,&
221 0 : starsq,denIn1Im,vC1,.FALSE.,vC1Im,denIn1,iDtype,iDir,[0,0],sigma_loc)
222 : END IF
223 :
224 0 : CALL timestop("Generation of potential perturbation")
225 :
226 : #ifdef CPP_MPI
227 0 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
228 : #endif
229 :
230 : ! Add the potential gradient to the potential perturbation in the displaced MT
231 0 : IF (strho.AND.fi%juphon%l_phonon) THEN
232 0 : vTot1%mt(:,0:,iDtype,1) = vTot1%mt(:,0:,iDtype,1) + grVext%mt(:,0:,iDtype,1)
233 0 : IF (fi%input%jspins==2) vTot1%mt(:,0:,iDtype,2) = vTot1%mt(:,0:,iDtype,2) + grVext%mt(:,0:,iDtype,1)
234 0 : IF (l_minusq) THEN
235 0 : vTot1m%mt(:,0:,iDtype,1) = vTot1m%mt(:,0:,iDtype,1) + grVext%mt(:,0:,iDtype,1)
236 0 : IF (fi%input%jspins==2) vTot1m%mt(:,0:,iDtype,2) = vTot1m%mt(:,0:,iDtype,2) + grVext%mt(:,0:,iDtype,1)
237 : END IF
238 : ELSE
239 0 : IF (fi%juphon%l_phonon) vTot1%mt(:,0:,iDtype,:) = vTot1%mt(:,0:,iDtype,:) + grVtot%mt(:,0:,iDtype,:)
240 0 : realiter = realiter + 1
241 0 : IF (l_minusq) THEN
242 0 : IF (fi%juphon%l_phonon) vTot1m%mt(:,0:,iDtype,:) = vTot1m%mt(:,0:,iDtype,:) + grVtot%mt(:,0:,iDtype,:)
243 : END IF
244 : END IF
245 :
246 0 : CALL vTot1%distribute(fmpi%mpi_comm)
247 :
248 0 : IF (l_minusq) THEN
249 0 : CALL vTot1m%distribute(fmpi%mpi_comm)
250 : END IF
251 :
252 : #ifdef CPP_MPI
253 0 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
254 : #endif
255 :
256 : ! Calculate the perturbed expansion coefficients z1 --> saved to results1
257 0 : CALL timestart("dfpt eigen")
258 0 : IF (.NOT.final_SH_it) THEN
259 : CALL dfpt_eigen(fi, sphhar, results, resultsq, results1, fmpi, enpara, nococonv, starsq, vTot1, vTot1Im, &
260 0 : vTot, rho, bqpt, eig_id, q_eig_id, dfpt_eig_id, iDir, iDtype, killcont, l_real, .TRUE.)
261 : ELSE
262 : CALL dfpt_eigen(fi, sphhar, results, resultsq, results1, fmpi, enpara, nococonv, starsq, vTot1, vTot1Im, &
263 0 : vTot, rho, bqpt, eig_id, q_eig_id, dfpt_eig_id, iDir, iDtype, killcont, l_real, .FALSE., dfpt_eig_id2)
264 : END IF
265 0 : CALL timestop("dfpt eigen")
266 :
267 0 : IF (l_minusq) THEN
268 0 : CALL timestart("dfpt minus eigen")
269 0 : IF (.NOT.final_SH_it) THEN
270 : CALL dfpt_eigen(fi, sphhar, results, resultsmq, results1m, fmpi, enpara, nococonv, starsmq, vTot1m, vTot1mIm, &
271 0 : vTot, rho, bmqpt, eig_id, qm_eig_id, dfpt_eigm_id, iDir, iDtype, killcont, l_real, .TRUE.)
272 : ELSE
273 : CALL dfpt_eigen(fi, sphhar, results, resultsmq, results1m, fmpi, enpara, nococonv, starsmq, vTot1m, vTot1mIm, &
274 0 : vTot, rho, bmqpt, eig_id, qm_eig_id, dfpt_eigm_id, iDir, iDtype, killcont, l_real, .FALSE., dfpt_eigm_id2)
275 : END IF
276 0 : CALL timestop("dfpt minus eigen")
277 : END IF
278 :
279 : #ifdef CPP_MPI
280 0 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
281 : #endif
282 :
283 : ! If q=0, the eigenenergy perturbation is /=0 and if we do not
284 : ! look at a semiconductor or an insulator, there will be a
285 : ! perturbation of the occupation numbers as well.
286 0 : CALL timestart("Fermi energy and occupation derivative")
287 0 : IF (norm2(bqpt)<1e-8.AND.ABS(results%tkb_loc)>1e-12) THEN
288 0 : CALL dfpt_fermie(eig_id,dfpt_eig_id,fmpi,fi%kpts,fi%input,fi%noco,results,results1)
289 : ELSE
290 0 : results1%ef = 0.0
291 0 : results1%w_iks = 0.0
292 : END IF
293 0 : CALL timestop("Fermi energy and occupation derivative")
294 :
295 0 : IF (l_minusq) THEN
296 0 : CALL timestart("Fermi energy and occupation minus derivative")
297 0 : IF (norm2(bqpt)<1e-8) THEN
298 0 : CALL dfpt_fermie(eig_id,dfpt_eigm_id,fmpi,fi%kpts,fi%input,fi%noco,results,results1m)
299 : ELSE
300 0 : results1m%ef = 0.0
301 0 : results1m%w_iks = 0.0
302 : END IF
303 0 : CALL timestop("Fermi energy and occupation minus derivative")
304 : END IF
305 :
306 : #ifdef CPP_MPI
307 0 : CALL MPI_BCAST(results1%ef, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
308 0 : CALL MPI_BCAST(results1%w_iks, SIZE(results1%w_iks), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
309 : #endif
310 :
311 0 : IF (l_minusq) THEN
312 : #ifdef CPP_MPI
313 0 : CALL MPI_BCAST(results1m%ef, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
314 0 : CALL MPI_BCAST(results1m%w_iks, SIZE(results1m%w_iks), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
315 : #endif
316 : END IF
317 :
318 : ! Exit here after the new final eigenvalues have been generated
319 : ! and before a new density perturbation is built
320 0 : IF (final_SH_it) THEN
321 0 : IF (fi%juphon%l_phonon) denIn1%mt(:,0:,iDtype,:) = denIn1%mt(:,0:,iDtype,:) + grRho%mt(:,0:,iDtype,:)
322 0 : CALL denIn1%distribute(fmpi%mpi_comm)
323 :
324 0 : IF (l_minusq) THEN
325 0 : IF (fi%juphon%l_phonon) denIn1m%mt(:,0:,iDtype,:) = denIn1m%mt(:,0:,iDtype,:) + grRho%mt(:,0:,iDtype,:)
326 0 : CALL denIn1m%distribute(fmpi%mpi_comm)
327 : END IF
328 :
329 0 : l_cont = .FALSE.
330 0 : IF (fmpi%irank==0) write(*,*) "Final Sternheimer iteration finished."
331 0 : CALL timestop("Sternheimer Iteration")
332 0 : CYCLE scfloop
333 : END IF
334 :
335 : ! Generate the new density perturbation
336 0 : CALL timestart("generation of new charge density (total)")
337 0 : CALL denOut1%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN, l_dfpt=.TRUE.)
338 0 : CALL denOut1Im%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN, l_dfpt=.FALSE.)
339 0 : denOut1%iter = denIn1%iter
340 0 : IF (.NOT.l_minusq) THEN
341 : CALL dfpt_cdngen(eig_id,dfpt_eig_id,fmpi,fi%input,banddosdummy,fi%vacuum,&
342 : fi%kpts,fi%atoms,sphhar,starsq,fi%sym,fi%juphon,fi%gfinp,fi%hub1inp,&
343 : enpara,fi%cell,fi%noco,nococonv,vTot,results,results1,&
344 0 : archiveType,xcpot,denOut1,denOut1Im,bqpt,iDtype,iDir,l_real)
345 : ELSE
346 : CALL dfpt_cdngen(eig_id,dfpt_eig_id,fmpi,fi%input,banddosdummy,fi%vacuum,&
347 : fi%kpts,fi%atoms,sphhar,starsq,fi%sym,fi%juphon,fi%gfinp,fi%hub1inp,&
348 : enpara,fi%cell,fi%noco,nococonv,vTot,results,results1,&
349 : archiveType,xcpot,denOut1,denOut1Im,bqpt,iDtype,iDir,l_real,&
350 0 : qm_eig_id,dfpt_eigm_id,starsmq,results1m)
351 : END IF
352 0 : CALL timestop("generation of new charge density (total)")
353 :
354 0 : IF (l_minusq) THEN
355 0 : CALL timestart("generation of new charge density (total)")
356 0 : CALL denOut1m%init(starsmq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN, l_dfpt=.TRUE.)
357 0 : CALL denOut1mIm%init(starsmq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN, l_dfpt=.FALSE.)
358 : CALL dfpt_cdngen(eig_id,dfpt_eigm_id,fmpi,fi%input,banddosdummy,fi%vacuum,&
359 : fi%kpts,fi%atoms,sphhar,starsmq,fi%sym,fi%juphon,fi%gfinp,fi%hub1inp,&
360 : enpara,fi%cell,fi%noco,nococonv,vTot,results,results1m,&
361 : archiveType,xcpot,denOut1m,denOut1mIm,-bqpt,iDtype,iDir,l_real,&
362 0 : q_eig_id,dfpt_eig_id,starsq,results1)
363 0 : CALL timestop("generation of new charge density (total)")
364 : END IF
365 :
366 : ! If a starting density perturbation was to be generated, subtract the density
367 : ! gradient and exit here; no mixing yet!
368 0 : IF (strho) THEN
369 0 : strho = .FALSE.
370 0 : denIn1 = denOut1
371 0 : denIn1Im = denOut1Im
372 0 : IF (fi%juphon%l_phonon) denIn1%mt(:,0:,iDtype,:) = denIn1%mt(:,0:,iDtype,:) - grRho%mt(:,0:,iDtype,:)
373 0 : IF (fmpi%irank==0) write(*,*) "Starting perturbation generated."
374 0 : CALL timestop("Sternheimer Iteration")
375 0 : IF (l_minusq) THEN
376 0 : denIn1m = denOut1m
377 0 : denIn1mIm = denOut1mIm
378 0 : IF (fi%juphon%l_phonon) denIn1m%mt(:,0:,iDtype,:) = denIn1m%mt(:,0:,iDtype,:) - grRho%mt(:,0:,iDtype,:)
379 : END IF
380 : CYCLE scfloop
381 : END IF
382 :
383 : ! If a the first full density perturbation was to be generated, subtract the density
384 : ! gradietn and exit here; no mixing yet!
385 0 : IF (.NOT.onedone) THEN
386 0 : onedone = .TRUE.
387 0 : denIn1 = denOut1
388 0 : denIn1Im = denOut1Im
389 0 : IF (fi%juphon%l_phonon) denIn1%mt(:,0:,iDtype,:) = denIn1%mt(:,0:,iDtype,:) - grRho%mt(:,0:,iDtype,:)
390 0 : IF (fmpi%irank==0) write(*,*) "1st 'real' density perturbation generated."
391 0 : CALL timestop("Sternheimer Iteration")
392 0 : IF (l_minusq) THEN
393 0 : denIn1m = denOut1m
394 0 : denIn1mIm = denOut1mIm
395 0 : IF (fi%juphon%l_phonon) denIn1m%mt(:,0:,iDtype,:) = denIn1m%mt(:,0:,iDtype,:) - grRho%mt(:,0:,iDtype,:)
396 : END IF
397 : CYCLE scfloop
398 : END IF
399 :
400 : ! No longer exit here.
401 : !IF (final_SH_it) THEN
402 : ! denIn1 = denOut1
403 : ! denIn1Im = denOut1Im
404 : ! l_cont = .FALSE.
405 : ! IF (fmpi%irank==0) write(*,*) "Final Sternheimer iteration finished."
406 : ! CALL timestop("Sternheimer Iteration")
407 : ! IF (l_minusq) THEN
408 : ! denIn1m = denOut1m
409 : ! denIn1mIm = denOut1mIm
410 : ! END IF
411 : ! CYCLE scfloop
412 : !END IF
413 :
414 0 : CALL denIn1%distribute(fmpi%mpi_comm)
415 0 : CALL denIn1Im%distribute(fmpi%mpi_comm)
416 :
417 0 : IF (l_minusq) THEN
418 0 : CALL denIn1m%distribute(fmpi%mpi_comm)
419 0 : CALL denIn1mIm%distribute(fmpi%mpi_comm)
420 : END IF
421 :
422 : #ifdef CPP_MPI
423 0 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
424 : #endif
425 :
426 0 : field2 = fi%field
427 :
428 : ! First mixing in the 2nd "real" iteration.
429 : ! Remove the gradient from denIn as to not mix identical big numbers.
430 0 : IF (fi%juphon%l_phonon) denIn1%mt(:,0:,iDtype,:) = denIn1%mt(:,0:,iDtype,:) + grRho%mt(:,0:,iDtype,:)
431 0 : CALL denIn1%distribute(fmpi%mpi_comm)
432 :
433 0 : IF (l_minusq) THEN
434 0 : IF (fi%juphon%l_phonon) denIn1m%mt(:,0:,iDtype,:) = denIn1m%mt(:,0:,iDtype,:) + grRho%mt(:,0:,iDtype,:)
435 0 : CALL denIn1m%distribute(fmpi%mpi_comm)
436 : END IF
437 :
438 : #ifdef CPP_MPI
439 0 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
440 : #endif
441 :
442 : ! Mix input and output densities
443 0 : CALL timestart("DFPT mixing")
444 : CALL mix_charge(field2, fmpi, (iter == fi%input%itmax .OR. judft_was_argument("-mix_io")), starsq, &
445 : fi%atoms, sphhar, fi%vacuum, fi%input, fi%sym, fi%juphon, fi%cell, fi%noco, nococonv, &
446 : archiveType, xcpot, iter, denIn1, denOut1, results1, .FALSE., fi%sliceplot,&
447 0 : denIn1Im, denOut1Im, dfpt_tag)
448 0 : CALL timestop("DFPT mixing")
449 :
450 0 : IF (fi%juphon%l_phonon) denIn1%mt(:,0:,iDtype,:) = denIn1%mt(:,0:,iDtype,:) - grRho%mt(:,0:,iDtype,:)
451 :
452 0 : CALL denIn1%distribute(fmpi%mpi_comm)
453 :
454 0 : IF (l_minusq) THEN
455 0 : CALL timestart("DFPT mixing")
456 : CALL mix_charge(field2, fmpi, (iter == fi%input%itmax .OR. judft_was_argument("-mix_io")), starsmq, &
457 : fi%atoms, sphhar, fi%vacuum, fi%input, fi%sym, fi%juphon, fi%cell, fi%noco, nococonv, &
458 : archiveType, xcpot, iterm, denIn1m, denOut1m, results1m, .FALSE., fi%sliceplot,&
459 0 : denIn1mIm, denOut1mIm, dfpt_tag//"m")
460 0 : CALL timestop("DFPT mixing")
461 :
462 0 : IF (fi%juphon%l_phonon) denIn1m%mt(:,0:,iDtype,:) = denIn1m%mt(:,0:,iDtype,:) - grRho%mt(:,0:,iDtype,:)
463 :
464 0 : CALL denIn1m%distribute(fmpi%mpi_comm)
465 : END IF
466 :
467 : #ifdef CPP_MPI
468 0 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
469 : #endif
470 :
471 0 : IF (fmpi%irank==0) THEN
472 0 : WRITE (oUnit, FMT=8130) iter
473 : 8130 FORMAT(/, 5x, '******* it=', i3, ' is completed********', /,/)
474 0 : WRITE (*, *) "Iteration:", iter, " Distance:", results1%last_distance
475 : END IF ! fmpi%irank==0
476 0 : CALL timestop("Sternheimer Iteration")
477 :
478 : #ifdef CPP_MPI
479 0 : CALL MPI_BCAST(results1%last_distance, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
480 0 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
481 : #endif
482 :
483 0 : IF (l_minusq) THEN
484 : #ifdef CPP_MPI
485 0 : CALL MPI_BCAST(results1m%last_distance, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
486 0 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
487 : #endif
488 : END IF
489 :
490 0 : l_cont = l_cont .AND. (iter < fi%input%itmax)
491 0 : l_cont = l_cont .AND. ((fi%input%mindistance <= results1%last_distance))
492 :
493 0 : final_SH_it = fi%input%mindistance > results1%last_distance
494 0 : l_cont = l_cont .OR. final_SH_it ! DO one more iteration so V1, z1 and rho1 match
495 :
496 : END DO scfloop ! DO WHILE (l_cont)
497 :
498 0 : CALL add_usage_data("Iterations", iter)
499 :
500 0 : END SUBROUTINE dfpt_sternheimer
501 : END MODULE m_dfpt_sternheimer
|
