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 : MODULE m_mix
7 :
8 : !------------------------------------------------------------------------
9 : ! mixing of charge densities or potentials:
10 : ! IMIX = 0 : linear mixing
11 : ! IMIX = 3 : Broyden's First method
12 : ! IMIX = 5 : Broyden's Second method
13 : ! IMIX = 7 : Generalized Anderson method
14 : !------------------------------------------------------------------------
15 :
16 : contains
17 :
18 654 : SUBROUTINE mix_charge( field, fmpi, l_writehistory,&
19 : stars, atoms, sphhar, vacuum, input, sym, juphon, cell, noco, nococonv,&
20 : archiveType, xcpot, iteration, inDen, outDen, results, l_runhia, sliceplot,&
21 : inDenIm, outDenIm, dfpt_tag)
22 :
23 : use m_juDFT
24 : use m_constants
25 : use m_cdn_io
26 : use m_stmix
27 : use m_broyden
28 : use m_qfix
29 : use m_types
30 : use m_umix
31 : use m_vmix
32 : use m_checkMMPmat
33 : USE m_kerker
34 : use m_pulay
35 : use m_a_pulay
36 : use m_types_mixvector
37 : USE m_distance
38 : use m_mixing_history
39 : use m_RelaxSpinAxisMagn
40 : USE m_plot
41 : implicit none
42 :
43 :
44 : type(t_input), intent(in) :: input
45 : type(t_vacuum), intent(in) :: vacuum
46 : type(t_noco), intent(in) :: noco
47 : type(t_nococonv), intent(in) :: nococonv
48 : TYPE(t_sym),TARGET,INTENT(in) :: sym
49 : TYPE(t_juphon), INTENT(in) :: juphon
50 : TYPE(t_stars),TARGET,INTENT(in) :: stars
51 : TYPE(t_cell),TARGET,INTENT(in) :: cell
52 : TYPE(t_sphhar),TARGET,INTENT(in) :: sphhar
53 : type(t_field), intent(inout) :: field
54 : TYPE(t_sliceplot), INTENT(IN) :: sliceplot
55 :
56 : type(t_mpi), intent(in) :: fmpi
57 : TYPE(t_atoms),TARGET,INTENT(in) :: atoms
58 : class(t_xcpot), intent(in) :: xcpot
59 : type(t_potden), intent(inout) :: outDen
60 : type(t_results), intent(inout) :: results
61 : type(t_potden), intent(inout) :: inDen
62 : integer, intent(in) :: archiveType
63 : integer, intent(inout) :: iteration
64 : LOGICAL, INTENT(IN) :: l_writehistory
65 : LOGICAL, INTENT(IN) :: l_runhia
66 :
67 : type(t_potden), OPTIONAL, INTENT(INOUT) :: inDenIm, outDenIm
68 :
69 : CHARACTER(len=20), OPTIONAL, INTENT(IN) :: dfpt_tag
70 :
71 : real :: fix
72 : type(t_potden) :: resDen, vYukawa
73 654 : TYPE(t_mixvector),ALLOCATABLE :: sm(:), fsm(:)
74 654 : TYPE(t_mixvector) :: fsm_mag
75 : LOGICAL :: l_densitymatrix, l_firstItU, l_firstItV, l_densitymatrixV, ldavLinMix, l_dfpt
76 : INTEGER :: it,maxiter
77 : INTEGER :: indStart_noDenmatmixing, indEnd_noDenmatmixing
78 :
79 :
80 654 : CALL timestart("Charge Density Mixing")
81 654 : l_densitymatrix = .FALSE.
82 654 : l_densitymatrixV = .FALSE.
83 654 : l_firstItU = .FALSE.
84 654 : l_firstItV = .FALSE.
85 654 : ldavLinMix = .FALSE.
86 654 : l_dfpt = PRESENT(dfpt_tag)
87 : !The density/potential matrices for DFT+U are split into two parts
88 : ! 1:atoms%n_u Are the elements for normal DFT+U
89 : ! atoms%n_u+1:atoms%n_u+atoms%n_hia are the elements for DFT+Hubbard 1
90 : ! atoms%n_u+atoms%n_hia+1:atoms%n_u+atoms%n_hia+atoms%n_opc are the elements for DFT+OPC
91 : !The latter are never mixed and held constant
92 654 : indStart_noDenmatmixing = atoms%n_u + 1
93 654 : indEnd_noDenmatmixing = atoms%n_u + atoms%n_hia + atoms%n_opc
94 :
95 654 : IF (atoms%n_u>0) THEN
96 1610 : l_firstItU = ALL(inDen%mmpMat(:,:,1:atoms%n_u,:)==0.0)
97 62 : l_densitymatrix=.NOT.input%ldaulinmix.AND..NOT.l_firstItU
98 62 : IF (fmpi%irank==0) CALL u_mix(input,atoms,noco,inDen%mmpMat,outDen%mmpMat)
99 : ENDIF
100 :
101 654 : IF (atoms%n_v.GT.0) THEN
102 0 : l_firstItV = ALL(inDen%nIJ_llp_mmp(:,:,:,:).EQ.0.0)
103 0 : l_densitymatrixV = .NOT.ldavlinmix.AND..NOT.l_firstItV
104 0 : IF (fmpi%irank==0) CALL v_mix(input,atoms,noco,inDen%nIJ_llp_mmp,outDen%nIJ_llp_mmp)
105 : ENDIF
106 :
107 654 : CALL timestart("Reading of distances")
108 654 : IF (iteration==1) CALL mixvector_reset(.TRUE.)
109 654 : CALL mixvector_init(fmpi%mpi_comm,l_densitymatrix,l_densitymatrixV,input,vacuum,noco,stars,cell,sphhar,atoms,sym,l_dfpt)
110 654 : CALL timestart("read history")
111 654 : IF (.NOT.l_dfpt) THEN
112 654 : CALL mixing_history_open(fmpi,input%maxiter)
113 : ELSE
114 0 : CALL mixing_history_open(fmpi,input%maxiter,dfpt_tag)
115 : END IF
116 :
117 654 : CALL timestop("read history")
118 654 : maxiter=MERGE(1,input%maxiter,input%imix==0)
119 654 : IF (.NOT.l_dfpt) THEN
120 654 : CALL mixing_history(input%imix,maxiter,inden,outden,sm,fsm,it,vacuum%nmzxyd)
121 : ELSE
122 0 : IF (iteration==1) CALL dfpt_mixing_history_reset()
123 0 : CALL mixing_history(input%imix,maxiter,inden,outden,sm,fsm,it,vacuum%nmzxyd,inDenIm,outDenIm)
124 : END IF
125 :
126 654 : IF (.NOT.l_dfpt) THEN
127 654 : CALL distance(fmpi%irank,cell%vol,input%jspins,vacuum%nmzxyd,fsm(it),inDen,outDen,results,fsm_Mag)
128 : ELSE
129 : ! TODO: For now dfpt_distance handles Re/Im separately. Maybe that is not all to necessary.
130 0 : CALL dfpt_distance(fmpi%irank,cell%vol,input%jspins,vacuum%nmzxyd,fsm(it),inDen,outDen,inDenIm,outDenIm,results,fsm_Mag)
131 : END IF
132 654 : CALL timestop("Reading of distances")
133 :
134 654 : IF (.NOT.l_dfpt) THEN
135 : ! Preconditioner for relaxation of Magnetic moments
136 654 : call precond_noco(it,vacuum,sphhar,stars,sym ,cell,noco,nococonv,input,atoms,inden,outden,fsm(it))
137 : END IF
138 :
139 : ! KERKER PRECONDITIONER
140 654 : IF( input%preconditioning_param /= 0 .AND. .NOT.l_dfpt) THEN
141 6 : CALL timestart("Preconditioner")
142 : CALL kerker( field, fmpi, &
143 : stars, atoms, sphhar, vacuum, input, sym, juphon, cell, noco, nococonv,&
144 6 : inDen, outDen, fsm(it) )
145 : !Store modified density in history
146 6 : CALL mixing_history_store(fsm(it))
147 6 : CALL timestop("Preconditioner")
148 : END IF
149 :
150 654 : if (atoms%n_u>0.and.fmpi%irank.ne.0.and.input%ldaulinmix) inden%mmpMat(:,:,:atoms%n_u,:)=0.0
151 : if (atoms%n_v>0.and.fmpi%irank.ne.0.and.ldavlinmix) inden%nIJ_llp_mmp(:,:,:,:) = 0.0
152 :
153 : !mixing of the densities
154 654 : CALL timestart("Mixing")
155 0 : SELECT CASE(input%imix)
156 : CASE(0)
157 0 : IF (fmpi%irank==0) WRITE(oUnit, fmt='(a,f10.5,a,f10.5)' ) &
158 0 : 'STRAIGHT MIXING: alpha=',input%alpha," spin-mixing=",MERGE(input%alpha*input%spinf,0.,input%jspins>1)
159 0 : CALL stmix(atoms,input,noco,fsm(it),fsm_mag,sm(it))
160 : CASE(3,5)
161 0 : CALL judft_error("Broyden 1/2 method not implemented")
162 : CASE(7)
163 654 : IF (fmpi%irank==0) WRITE(oUnit, fmt='(a,f10.5,a,i0)' ) &
164 327 : 'GENERALIZED ANDERSON MIXING: alpha=',input%alpha," History-length=",it-1
165 654 : Call broyden(input%alpha,fsm,sm,l_dfpt)
166 : CASE(9)
167 0 : IF (fmpi%irank==0) WRITE(oUnit, fmt='(a,f10.5,a,i0,a,i0)' ) &
168 0 : 'PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
169 0 : CALL pulay(input%alpha,fsm,sm,0,l_dfpt)
170 : CASE(11)
171 0 : IF (fmpi%irank==0) WRITE(oUnit, fmt='(a,f10.5,a,i0,a,i0)' ) &
172 0 : 'PERIODIC PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
173 0 : CALL pulay(input%alpha,fsm,sm,input%maxiter,l_dfpt)
174 : CASE(13)
175 0 : IF (fmpi%irank==0) WRITE(oUnit, fmt='(a,f10.5,a,i0,a,i0)' ) &
176 0 : 'RESTARTED PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
177 0 : CALL pulay(input%alpha,fsm,sm,0,l_dfpt)
178 0 : IF (it==input%maxiter) CALL mixing_history_limit(0) !Restarting Pulay
179 : CASE(15)
180 0 : IF (fmpi%irank==0) WRITE(oUnit, fmt='(a,f10.5,a,i0,a,i0)' ) &
181 0 : 'ADAPTED PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
182 0 : CALL a_pulay(input%alpha,fsm,sm,l_dfpt)
183 : CASE DEFAULT
184 654 : CALL judft_error("Unknown Mixing schema")
185 : END SELECT
186 654 : CALL timestop("Mixing")
187 :
188 654 : CALL timestart("Postprocessing")
189 : !extracte mixed density
190 47833352 : inDen%pw=0.0;inDen%mt=0.0
191 654 : IF (l_dfpt) inDenIm%mt=0.0
192 16449876 : IF (ALLOCATED(inDen%vac)) inden%vac=0.0
193 12602 : IF (ALLOCATED(inDen%mmpMat).AND.l_densitymatrix) inden%mmpMat(:,:,:atoms%n_u,:)=0.0
194 654 : IF (.NOT.l_dfpt) THEN
195 654 : CALL sm(it)%to_density(inDen,vacuum%nmzxyd)
196 : ELSE
197 0 : CALL sm(it)%to_density(inDen,vacuum%nmzxyd,inDenIm)
198 : END IF
199 :
200 654 : IF (atoms%n_u>0.AND.l_firstItU) THEN
201 : !No density matrix was present
202 : !but is now created...
203 922 : inden%mmpMAT(:,:,:atoms%n_u,:)=outden%mmpMat(:,:,:atoms%n_u,:)
204 : !Delete the history without U
205 4 : CALL mixing_history_reset(fmpi)
206 4 : CALL mixvector_reset()
207 : ENDIF
208 :
209 654 : IF (atoms%n_v>0.AND.l_firstItV) THEN
210 : !No density matrix was present
211 : !but is now created...
212 0 : inden%nIJ_llp_mmp(:,:,:,:) = outden%nIJ_llp_mmp(:,:,:,:)
213 : !Delete the history without U
214 0 : CALL mixing_history_reset(fmpi)
215 0 : CALL mixvector_reset()
216 : ENDIF
217 :
218 :
219 654 : IF(atoms%n_u>0) THEN
220 : !When the mixing of the density matrix is done together
221 : !with the charge density depending on the mixing scheme
222 : !it can become unstable
223 : !Check whether the mixed density matrix makes sense
224 : !And correct invalid elements
225 62 : CALL checkMMPmat(1,atoms%n_u,l_densitymatrix,fmpi,atoms,input,outden,inden)
226 : ENDIF
227 :
228 654 : IF(atoms%n_hia+atoms%n_opc>0) THEN
229 : !For LDA+HIA we don't use any mixing of the density matrices we just pass it on
230 : inDen%mmpMat(:,:,indStart_noDenmatmixing:indEnd_noDenmatmixing,:) = &
231 5544 : outDen%mmpMat(:,:,indStart_noDenmatmixing:indEnd_noDenmatmixing,:)
232 24 : IF(l_runhia) THEN
233 0 : iteration = 1
234 0 : CALL mixing_history_reset(fmpi)
235 0 : CALL mixvector_reset()
236 : ENDIF
237 : ENDIF
238 :
239 654 : if(iteration == 1 .and. xcpot%vx_is_MetaGGA()) then
240 0 : CALL mixing_history_reset(fmpi)
241 0 : CALL mixvector_reset()
242 : endif
243 :
244 654 : call timestart("qfix")
245 : !fix charge of the new density
246 654 : IF (fmpi%irank==0.AND..NOT.l_dfpt) CALL qfix(fmpi,stars,nococonv,atoms,sym,vacuum, sphhar,input,cell ,inDen,noco%l_noco,.FALSE.,.FALSE.,.FALSE., fix)
247 654 : call timestop("qfix")
248 :
249 654 : IF(vacuum%nvac.EQ.1) THEN
250 54 : IF (sym%invs) THEN
251 923716 : inDen%vac(:,:,2,:) = CONJG(inDen%vac(:,:,1,:))
252 : ELSE
253 5632520 : inDen%vac(:,:,2,:) = inDen%vac(:,:,1,:)
254 : END IF
255 : END IF
256 :
257 654 : call timestart("Density output")
258 : !write out mixed density (but not for a plotting run)
259 654 : IF ((fmpi%irank==0).AND.(sliceplot%iplot==0).AND..NOT.l_dfpt) THEN
260 : CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,CDN_INPUT_DEN_const,&
261 327 : 1,results%last_distance,results%ef,results%last_mmpmatDistance,results%last_occDistance,.TRUE.,inDen)
262 327 : ELSE IF ((fmpi%irank==0).AND.(sliceplot%iplot==0).AND.l_dfpt) THEN
263 : CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,CDN_INPUT_DEN_const,&
264 0 : 1,results%last_distance,results%ef,results%last_mmpmatDistance,results%last_occDistance,.TRUE.,inDen,inFilename=dfpt_tag,denIm=inDenIm)
265 : END IF
266 :
267 : #ifdef CPP_HDF
268 : ! TODO: Could be a neat option for DFPT as well.
269 654 : IF (fmpi%irank==0.and.judft_was_argument("-last_extra").AND..NOT.l_dfpt) THEN
270 0 : CALL system("rm cdn_last.hdf")
271 : CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,CDN_INPUT_DEN_const,&
272 : 1,results%last_distance,results%ef,results%last_mmpmatDistance,results%last_occDistance,.TRUE.,&
273 0 : inDen,inFilename='cdn_last')
274 0 : CALL writeCoreDensity(input,atoms,inDen%mtCore,inDen%tec,inDen%qint,'cdn_last')
275 : END IF
276 : #endif
277 654 : call timestop("Density output")
278 654 : inDen%iter = inDen%iter + 1
279 :
280 654 : IF (.NOT.l_dfpt) THEN
281 654 : IF (l_writehistory.AND.input%imix.NE.0) CALL mixing_history_close(fmpi)
282 : ELSE
283 0 : IF (l_writehistory.AND.input%imix.NE.0) CALL mixing_history_close(fmpi,dfpt_tag)
284 : END IF
285 :
286 654 : CALL timestop("Postprocessing")
287 654 : CALL timestop("Charge Density Mixing")
288 9150 : END SUBROUTINE mix_charge
289 :
290 : END MODULE m_mix
|
