Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2018 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_rdmft
8 :
9 : CONTAINS
10 :
11 0 : SUBROUTINE rdmft(eig_id,fmpi,fi,enpara,stars,&
12 : sphhar,vTot,vCoul,nococonv,xcpot,mpdata,hybdat,&
13 : results,archiveType,outDen)
14 : use m_types_vacdos
15 : use m_work_package
16 : USE m_types
17 : USE m_juDFT
18 : USE m_constants
19 : USE m_intgr, ONLY : intgr3
20 : USE m_eig66_io
21 : #ifndef CPP_OLDINTEL
22 : USE m_cdnval
23 : USE m_cdngen
24 : USE m_cdn_io
25 : USE m_cdncore
26 : USE m_qfix
27 : USE m_vgen_coulomb
28 : USE m_convol
29 : USE m_intnv
30 :
31 : USE m_mixedbasis
32 : USE m_coulombmatrix
33 : USE m_hf_init
34 : USE m_hf_setup
35 : USE m_io_hybrid
36 : USE m_symm_hf
37 : USE m_exchange_valence_hf
38 : USE m_exchange_core
39 : USE m_symmetrizeh
40 : USE m_bfgs_b2
41 : USE m_xmlOutput
42 : USE m_types_dos
43 : use m_calc_cmt
44 :
45 : #endif
46 :
47 : IMPLICIT NONE
48 :
49 : TYPE(t_mpi), INTENT(IN) :: fmpi
50 : type(t_fleurinput), intent(in) :: fi
51 : TYPE(t_enpara), INTENT(INOUT) :: enpara
52 : TYPE(t_stars), INTENT(IN) :: stars
53 : TYPE(t_sphhar), INTENT(IN) :: sphhar
54 : TYPE(t_potden), INTENT(INOUT) :: vTot
55 : TYPE(t_potden), INTENT(INOUT) :: vCoul
56 : TYPE(t_nococonv), INTENT(INOUT) :: nococonv
57 : TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
58 : TYPE(t_mpdata), intent(inout) :: mpdata
59 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
60 : TYPE(t_results), INTENT(INOUT) :: results
61 : TYPE(t_potden), INTENT(INOUT) :: outDen
62 :
63 : INTEGER, INTENT(IN) :: eig_id
64 : INTEGER, INTENT(IN) :: archiveType
65 0 : TYPE(t_potden) :: EnergyDen
66 :
67 : #ifndef CPP_OLDINTEL
68 0 : TYPE(t_cdnvalJob) :: cdnvalJob
69 0 : TYPE(t_potden) :: singleStateDen, overallDen, overallVCoul, vTotTemp
70 0 : TYPE(t_regionCharges) :: regCharges
71 0 : TYPE(t_dos) :: dos
72 0 : TYPE(t_vacdos) :: vacdos
73 0 : TYPE(t_moments) :: moments
74 0 : TYPE(t_mat) :: exMat, zMat, olap, trafo, invtrafo, tmpMat, exMatLAPW
75 0 : TYPE(t_lapw) :: lapw
76 0 : type(t_work_package) :: work_pack
77 : INTEGER :: ikpt, ikpt_i, iBand, jkpt, jBand, iAtom, na, itype, lh, iGrid
78 : INTEGER :: jspin, jspmax, jsp, isp, ispin, nbasfcn, nbands
79 : INTEGER :: nsymop, ikptf, iterHF, ierr
80 : INTEGER :: iState, jState, iStep, numStates, numRelevantStates, convIter
81 : INTEGER :: maxHistoryLength
82 : INTEGER :: lastGroupEnd, currentGroupEnd
83 : REAL :: fix, potDenInt, fermiEnergyTemp, tempDistance, spinDegenFac
84 : REAL :: rdmftFunctionalValue, occStateI, gradSum
85 : REAL :: exchangeTerm, lagrangeMultiplier, equalityCriterion
86 : REAL :: mixParam, rdmftEnergy, occSum
87 : REAL :: sumOcc, addCharge, subCharge, addChargeWeight, subChargeWeight
88 : REAL :: rhs, totz, theta, temp
89 : REAL :: averageParam, averageGrad
90 : REAL, PARAMETER :: degenEps = 0.00001
91 : REAL, PARAMETER :: convCrit = 5.0e-6, occMixParam = 0.5
92 : REAL, PARAMETER :: minOcc = 1.0e-13, minOccB = 1.0e-5
93 : LOGICAL :: converged, l_qfix, l_restart
94 : CHARACTER(LEN=20) :: filename
95 : CHARACTER(LEN=20) :: attributes(3)
96 :
97 0 : INTEGER :: nsest(fi%input%neig) ! probably too large
98 0 : INTEGER :: rrot(3,3,fi%sym%nsym)
99 0 : INTEGER :: psym(fi%sym%nsym) ! Note: psym is only filled up to index nsymop
100 0 : INTEGER :: lowestState(fi%kpts%nkpt,fi%input%jspins)
101 0 : INTEGER :: highestState(fi%kpts%nkpt,fi%input%jspins)
102 0 : INTEGER :: neigTemp(fi%kpts%nkpt,fi%input%jspins)
103 :
104 0 : REAL :: wl_iks(fi%input%neig,fi%kpts%nkptf)
105 :
106 0 : REAL :: vmd(fi%atoms%ntype), zintn_r(fi%atoms%ntype), dpj(fi%atoms%jmtd), mt(fi%atoms%jmtd,fi%atoms%ntype)
107 :
108 0 : REAL, ALLOCATABLE :: overallVCoulSSDen(:,:,:)
109 0 : REAL, ALLOCATABLE :: vTotSSDen(:,:,:)
110 0 : REAL, ALLOCATABLE :: dEdOcc(:,:,:)
111 :
112 0 : REAL, ALLOCATABLE :: zintn_rSSDen(:,:,:)
113 0 : REAL, ALLOCATABLE :: vmdSSDen(:,:,:)
114 :
115 :
116 0 : REAL, ALLOCATABLE :: exDiag(:,:,:)
117 0 : REAL, ALLOCATABLE :: eig_irr(:,:)
118 :
119 0 : REAL, ALLOCATABLE :: gradient(:), lastGradient(:)
120 0 : REAL, ALLOCATABLE :: parameters(:), lastParameters(:)
121 0 : REAL, ALLOCATABLE :: minConstraints(:)
122 0 : REAL, ALLOCATABLE :: maxConstraints(:)
123 0 : REAL, ALLOCATABLE :: gradientCorrections(:,:)
124 0 : REAL, ALLOCATABLE :: paramCorrections(:,:)
125 0 : REAL, ALLOCATABLE :: equalityLinCombi(:)
126 :
127 0 : REAL, ALLOCATABLE :: occupationVec(:)
128 :
129 0 : INTEGER, ALLOCATABLE :: paramGroup(:)
130 0 : INTEGER, ALLOCATABLE :: indx_sest(:,:)
131 0 : INTEGER, ALLOCATABLE :: parent(:)
132 : INTEGER, ALLOCATABLE :: pointer_EIBZ(:)
133 0 : INTEGER, ALLOCATABLE :: n_q(:)
134 0 : complex, allocatable :: cmt_nk(:,:,:)
135 :
136 0 : LOGICAL, ALLOCATABLE :: enabledConstraints(:)
137 : type(t_hybmpi) :: glob_mpi
138 :
139 0 : complex :: c_phase(fi%input%neig)
140 : complex :: sigma_loc(2)
141 :
142 : #endif
143 :
144 : #ifndef CPP_OLDINTEL
145 :
146 0 : WRITE(*,*) 'entered RDMFT subroutine'
147 :
148 : ! General initializations
149 0 : mixParam = 0.001
150 0 : lagrangeMultiplier = 0.1 !results%ef
151 0 : spinDegenFac = 2.0 / fi%input%jspins ! This factor is used to compensate the missing second spin in non-spinpolarized calculations
152 :
153 0 : neigTemp(:,:) = results%neig(:,:)
154 :
155 : ! Determine which states have to be considered
156 0 : lowestState(:,:) = -1
157 0 : highestState(:,:) = -1
158 0 : DO jspin = 1, fi%input%jspins
159 0 : jsp = MERGE(1,jspin,fi%noco%l_noco)
160 0 : DO ikpt = 1, fi%kpts%nkpt
161 : ! determine lowest state
162 0 : DO iBand = 1, results%neig(ikpt,jsp)
163 0 : IF((results%w_iks(iBand,ikpt,jspin) / fi%kpts%wtkpt(ikpt)).LT.(1.0-fi%input%rdmftOccEps)) THEN
164 0 : lowestState(ikpt,jspin) = iBand
165 0 : EXIT
166 : END IF
167 : END DO
168 0 : lowestState(ikpt,jspin) = lowestState(ikpt,jspin) - fi%input%rdmftStatesBelow
169 0 : DO iBand = lowestState(ikpt,jspin)-1, 1, -1
170 0 : IF((results%eig(iBand+1,ikpt,jsp) - results%eig(iBand,ikpt,jsp)).GT.degenEps) THEN
171 : EXIT
172 : END IF
173 0 : lowestState(ikpt,jspin) = iBand
174 : END DO
175 0 : lowestState(ikpt,jspin) = MAX(lowestState(ikpt,jspin),1)
176 :
177 : ! determine highest state
178 0 : DO iBand = results%neig(ikpt,jsp)-1, 1, -1
179 0 : IF((results%w_iks(iBand,ikpt,jspin) / fi%kpts%wtkpt(ikpt)).GT.(0.0+fi%input%rdmftOccEps)) THEN
180 0 : highestState(ikpt,jspin) = iBand
181 0 : EXIT
182 : END IF
183 : END DO
184 0 : highestState(ikpt,jspin) = highestState(ikpt,jspin) + fi%input%rdmftStatesAbove
185 0 : IF((results%neig(ikpt,jsp)-1).LT.highestState(ikpt,jspin)) THEN
186 0 : WRITE(oUnit,*) 'Error: Not enough states calculated:'
187 0 : WRITE(oUnit,*) 'ikpt, jsp: ', ikpt, jsp
188 0 : WRITE(oUnit,*) 'highestState(ikpt,jspin): ', highestState(ikpt,jspin)
189 0 : WRITE(oUnit,*) 'results%neig(ikpt,jsp): ', results%neig(ikpt,jsp)
190 0 : CALL juDFT_error('Not enough states calculated', calledby = 'rdmft')
191 : END IF
192 0 : DO iBand = highestState(ikpt,jspin)+1, results%neig(ikpt,jsp)
193 0 : IF((results%eig(iBand,ikpt,jsp) - results%eig(iBand-1,ikpt,jsp)).GT.degenEps) THEN
194 : EXIT
195 : END IF
196 0 : highestState(ikpt,jspin) = iBand
197 : END DO
198 0 : IF(highestState(ikpt,jspin).EQ.results%neig(ikpt,jsp)) THEN
199 0 : WRITE(oUnit,*) 'Error: Highest state is degenerate:'
200 0 : WRITE(oUnit,*) 'ikpt, jsp: ', ikpt, jsp
201 0 : WRITE(oUnit,*) 'highestState(ikpt,jspin): ', highestState(ikpt,jspin)
202 0 : WRITE(oUnit,*) 'results%eig(highestState(ikpt,jspin),ikpt,jsp): ', results%eig(highestState(ikpt,jspin),ikpt,jsp)
203 0 : WRITE(oUnit,*) 'results%eig(highestState(ikpt,jspin)-1,ikpt,jsp): ', results%eig(highestState(ikpt,jspin)-1,ikpt,jsp)
204 0 : CALL juDFT_error('Highest state is degenerate', calledby = 'rdmft')
205 : END IF
206 : END DO
207 : END DO
208 :
209 0 : IF (ANY(results%w_iksRDMFT(:,:,:).NE.0.0)) THEN
210 0 : results%w_iks(:,:,:) = results%w_iksRDMFT(:,:,:)
211 : END IF
212 : ! results%w_iksRDMFT(:,:,:) = results%w_iks(:,:,:)
213 :
214 : ! Move occupations of relevant states well into allowed region
215 0 : numRelevantStates = SUM(highestState(:,:)) - SUM(lowestState(:,:)) + fi%input%jspins*fi%kpts%nkpt
216 0 : ALLOCATE(occupationVec(numRelevantStates))
217 0 : occupationVec(:) = 0.0
218 0 : sumOcc = 0.0
219 0 : addCharge = 0.0
220 0 : subCharge = 0.0
221 0 : addChargeWeight = 0.0
222 0 : subChargeWeight = 0.0
223 0 : iState = 0
224 0 : DO jspin = 1, fi%input%jspins
225 0 : jsp = MERGE(1,jspin,fi%noco%l_noco)
226 0 : DO ikpt = 1, fi%kpts%nkpt
227 0 : DO iBand = lowestState(ikpt,jspin), highestState(ikpt,jspin)
228 0 : iState = iState + 1
229 0 : occupationVec(iState) = results%w_iks(iBand,ikpt,jsp) / (fi%kpts%wtkpt(ikpt))
230 0 : sumOcc = sumOcc + results%w_iks(iBand,ikpt,jsp)
231 0 : IF(occupationVec(iState).LT.minOccB) THEN
232 0 : addCharge = addCharge + (minOccB-occupationVec(iState))*fi%kpts%wtkpt(ikpt)
233 0 : addChargeWeight = addChargeWeight + fi%kpts%wtkpt(ikpt)
234 : END IF
235 0 : IF(occupationVec(iState).GT.(1.0-minOccB)) THEN
236 0 : subCharge = subCharge + (occupationVec(iState)-(1.0-minOccB))*fi%kpts%wtkpt(ikpt)
237 0 : subChargeWeight = subChargeWeight + fi%kpts%wtkpt(ikpt)
238 : END IF
239 : END DO
240 : END DO
241 : END DO
242 0 : iState = 0
243 0 : DO jspin = 1, fi%input%jspins
244 0 : jsp = MERGE(1,jspin,fi%noco%l_noco)
245 0 : DO ikpt = 1, fi%kpts%nkpt
246 0 : DO iBand = lowestState(ikpt,jspin), highestState(ikpt,jspin)
247 0 : iState = iState + 1
248 0 : IF(occupationVec(iState).LT.minOccB) THEN
249 0 : occupationVec(iState) = occupationVec(iState) + 0.5*(subCharge+addCharge)*(fi%kpts%wtkpt(ikpt)/addChargeWeight)
250 : END IF
251 0 : IF(occupationVec(iState).GT.(1.0-minOccB)) THEN
252 0 : occupationVec(iState) = occupationVec(iState) - 0.5*(subCharge+addCharge)*(fi%kpts%wtkpt(ikpt)/subChargeWeight)
253 : END IF
254 : ! results%w_iks(iBand,ikpt,jsp) = occupationVec(iState) * fi%kpts%wtkpt(ikpt)
255 : END DO
256 : END DO
257 : END DO
258 0 : DEALLOCATE(occupationVec)
259 :
260 : ! Some more initializations
261 :
262 0 : results%neig(:,:) = highestState(:,:) + 1
263 :
264 0 : ALLOCATE(overallVCoulSSDen(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
265 0 : ALLOCATE(vTotSSDen(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
266 0 : ALLOCATE(dEdOcc(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
267 :
268 0 : ALLOCATE(zintn_rSSDen(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
269 0 : ALLOCATE(vmdSSDen(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
270 :
271 0 : zintn_rSSDen(:,:,:) = 0.0
272 0 : vmdSSDen(:,:,:) = 0.0
273 :
274 0 : CALL regCharges%init(fi%input,fi%atoms)
275 0 : CALL dos%init(fi%input,fi%atoms,fi%kpts,fi%banddos,results%eig)
276 0 : CALL vacdos%init(fi%input,fi%atoms,fi%kpts,fi%banddos,results%eig)
277 0 : CALL moments%init(fmpi,fi%input,sphhar,fi%atoms)
278 0 : CALL overallDen%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_DEN)
279 0 : CALL overallVCoul%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_POTCOUL)
280 0 : IF (ALLOCATED(vTot%pw_w)) DEALLOCATE (vTot%pw_w)
281 0 : ALLOCATE(vTot%pw_w(SIZE(overallDen%pw,1),fi%input%jspins))
282 0 : DO jspin = 1, fi%input%jspins
283 0 : CALL convol(stars,vTot%pw_w(:,jspin),vTot%pw(:,jspin))
284 : END DO
285 :
286 0 : CALL vTotTemp%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_POTTOT)
287 0 : vTotTemp = vTot
288 :
289 0 : DO jsp = 1,SIZE(vTot%mt,4)
290 0 : DO iAtom = 1,fi%atoms%ntype
291 0 : vTotTemp%mt(:fi%atoms%jri(iAtom),0,iAtom,jsp) = sfp_const * vTot%mt(:fi%atoms%jri(iAtom),0,iAtom,jsp) / fi%atoms%rmsh(:fi%atoms%jri(iAtom),iAtom)
292 : END DO
293 : END DO
294 :
295 0 : vCoul%pw_w = CMPLX(0.0,0.0)
296 0 : DO jspin = 1, fi%input%jspins
297 0 : CALL convol(stars,vCoul%pw_w(:,jspin),vCoul%pw(:,jspin))
298 : END DO
299 :
300 0 : vTotSSDen = 0.0
301 :
302 : ! Calculate all single state densities
303 :
304 0 : numStates = 0
305 0 : DO jspin = 1, fi%input%jspins
306 0 : jsp = MERGE(1,jspin,fi%noco%l_noco)
307 :
308 0 : CALL cdnvalJob%init(fmpi,fi%input,fi%kpts,fi%noco,results,jsp)
309 :
310 0 : DO ikpt_i = 1, SIZE(fmpi%k_list)
311 0 : ikpt= fmpi%k_list(ikpt_i)
312 0 : DO iBand = 1, highestState(ikpt,jsp)
313 0 : numStates = numStates + 1
314 : ! Construct cdnvalJob object for this state
315 : ! (Reasonable parallelization is not yet done - should be placed over the loops enclosing this section)
316 :
317 0 : cdnvalJob%k_list=[ikpt]
318 : ! cdnvalJob%ev_list=[iBand]
319 0 : cdnvalJob%weights(:,:) = 0.0
320 0 : cdnvalJob%weights(iBand,ikpt) = spinDegenFac
321 :
322 : ! Call cdnval to construct density
323 0 : WRITE(*,*) 'Note: some optional flags may have to be reset in rdmft before the cdnval call'
324 0 : WRITE(*,*) 'This is not yet implemented!'
325 0 : CALL singleStateDen%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_DEN)
326 : CALL cdnval(eig_id,fmpi,fi%kpts,jsp,fi%noco,nococonv,fi%input,fi%banddos,fi%cell,fi%atoms,enpara,stars,fi%vacuum,&
327 : sphhar,fi%sym,vTot, cdnvalJob,singleStateDen,regCharges,dos,vacdos,results,moments,&
328 0 : fi%gfinp,fi%hub1inp)
329 :
330 : ! Store the density on disc (These are probably way too many densities to keep them in memory)
331 0 : filename = ''
332 0 : WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
333 0 : IF (fmpi%irank.EQ.0) THEN
334 : CALL writeDensity(stars,fi%noco,fi%vacuum,fi%atoms,fi%cell,sphhar,fi%input,fi%sym, CDN_ARCHIVE_TYPE_CDN_const,CDN_input_DEN_const,&
335 0 : 0,-1.0,0.0,-1.0,-1.0,.FALSE.,singleStateDen,inFilename=TRIM(ADJUSTL(filename)))
336 : END IF
337 0 : CALL singleStateDen%distribute(fmpi%mpi_comm)
338 :
339 : ! For each state calculate Integral over KS effective potential times single state density
340 0 : potDenInt = 0.0
341 0 : CALL int_nv(jsp,stars,fi%vacuum,fi%atoms,sphhar,fi%cell,fi%sym,fi%input,vTotTemp,singleStateDen,potDenInt)
342 0 : vTotSSDen(iBand,ikpt,jsp) = potDenInt
343 :
344 0 : mt(:,:) = 0.0
345 0 : DO iType = 1, fi%atoms%ntype
346 0 : DO iGrid = 1, fi%atoms%jri(iType)
347 0 : mt(iGrid,iType) = singleStateDen%mt(iGrid,0,iType,jsp)
348 : END DO
349 :
350 0 : DO iGrid = 1, fi%atoms%jri(iType)
351 0 : dpj(iGrid) = mt(iGrid,iType)/fi%atoms%rmsh(iGrid,iType)
352 : END DO
353 :
354 0 : CALL intgr3(dpj,fi%atoms%rmsh(1,iType),fi%atoms%dx(iType),fi%atoms%jri(iType),rhs)
355 :
356 0 : zintn_r(iType) = fi%atoms%neq(iType)*fi%atoms%zatom(iType)*sfp_const*rhs/2.0
357 0 : zintn_rSSDen(iBand,ikpt,jsp) = zintn_rSSDen(iBand,ikpt,jsp) + zintn_r(iType)
358 :
359 0 : CALL intgr3(mt(1,iType),fi%atoms%rmsh(1,iType),fi%atoms%dx(iType),fi%atoms%jri(iType),totz)
360 :
361 : ! vmd(iType) = fi%atoms%rmt(iType)*vCoul%mt(fi%atoms%jri(iType),0,iType,1)/sfp_const + fi%atoms%zatom(iType) - totz*sfp_const
362 0 : vmd(iType) = -totz*sfp_const
363 0 : vmd(iType) = -fi%atoms%neq(iType)*fi%atoms%zatom(iType)*vmd(iType)/ (2.0*fi%atoms%rmt(iType))
364 :
365 0 : vmdSSDen(iBand,ikpt,jsp) = vmdSSDen(iBand,ikpt,jsp) + vmd(iType)
366 : END DO
367 :
368 : END DO
369 : END DO
370 : END DO
371 :
372 : ! Initializations for exchange contributions
373 :
374 0 : WRITE(*,*) 'RDMFT: HF initializations start'
375 :
376 0 : IF(ALLOCATED(hybdat%nbands)) DEALLOCATE(hybdat%nbands)
377 0 : IF(ALLOCATED(hybdat%nobd)) DEALLOCATE(hybdat%nobd)
378 0 : IF(ALLOCATED(hybdat%nbasm)) DEALLOCATE(hybdat%nbasm)
379 0 : IF(ALLOCATED(hybdat%div_vv)) DEALLOCATE(hybdat%div_vv)
380 0 : ALLOCATE(hybdat%nbasm(fi%kpts%nkptf))
381 0 : ALLOCATE(hybdat%div_vv(fi%input%neig,fi%kpts%nkpt,fi%input%jspins))
382 :
383 0 : iterHF = 0
384 0 : hybdat%l_calhf = .TRUE.
385 :
386 0 : CALL mixedbasis(fi%atoms,fi%kpts,fi%input,fi%cell,xcpot,fi%mpinp,mpdata,fi%hybinp, hybdat,enpara,fmpi,vTot, iterHF)
387 :
388 : !allocate coulomb matrix
389 0 : IF (.NOT.ALLOCATED(hybdat%coul)) ALLOCATE(hybdat%coul(fi%kpts%nkpt))
390 0 : DO ikpt = 1, fi%kpts%nkpt
391 0 : CALL hybdat%coul(ikpt)%alloc(fi, mpdata%num_radbasfn, mpdata%n_g, ikpt)
392 : END DO
393 :
394 0 : CALL glob_mpi%copy_mpi(fmpi)
395 0 : call work_pack%init(fi, hybdat, mpdata, glob_mpi, jsp, glob_mpi%rank, glob_mpi%size)
396 :
397 0 : CALL coulombmatrix(fmpi, fi, mpdata, hybdat, xcpot)
398 :
399 0 : DO ikpt = 1, fi%kpts%nkpt
400 0 : CALL hybdat%coul(ikpt)%mpi_bcast(fi, fmpi%mpi_comm, 0)
401 : END DO
402 :
403 0 : CALL hf_init(mpdata,fi,hybdat)
404 :
405 0 : WRITE(*,*) 'RDMFT: HF initializations end'
406 :
407 : maxHistoryLength = 7*numStates
408 : maxHistoryLength = 5*numStates
409 0 : maxHistoryLength = 23
410 :
411 0 : ALLOCATE(parent(fi%kpts%nkptf))
412 0 : ALLOCATE(exDiag(fi%input%neig,ikpt,fi%input%jspins))
413 0 : ALLOCATE(lastGradient(numStates+1))
414 0 : ALLOCATE(lastParameters(numStates+1))
415 0 : lastGradient = 0.0
416 0 : lastParameters = 0.0
417 0 : ALLOCATE(gradientCorrections(numStates+1,maxHistoryLength))
418 0 : ALLOCATE(paramCorrections(numStates+1,maxHistoryLength))
419 0 : gradientCorrections = 0.0
420 0 : paramCorrections = 0.0
421 0 : istep = 0
422 :
423 : ! Occupation number optimization loop
424 :
425 0 : convIter = 0
426 :
427 0 : converged = .FALSE.
428 0 : DO WHILE (.NOT.converged)
429 :
430 0 : WRITE(*,*) 'RDMFT: convergence loop start'
431 0 : convIter = convIter + 1
432 0 : WRITE(*,'(a,i7)') 'convIter: ', convIter
433 :
434 0 : DO jspin = 1, fi%input%jspins
435 0 : DO ikpt = 1,fi%kpts%nkpt
436 0 : WRITE(*,*) 'jspin, ikpt: ', jspin, ikpt
437 0 : WRITE(*,'(10f12.7)') results%w_iks(1:10,ikpt,jspin)
438 : END DO
439 : END DO
440 :
441 : ! Calculate overall density with current occupation numbers (don't forget core electron density)
442 0 : CALL overallDen%resetPotDen()
443 0 : jspmax = fi%input%jspins
444 0 : IF (fi%noco%l_mperp) jspmax = 1
445 :
446 0 : DO jspin = 1,jspmax
447 0 : CALL cdnvalJob%init(fmpi,fi%input,fi%kpts,fi%noco,results,jspin)
448 : CALL cdnval(eig_id,fmpi,fi%kpts,jspin,fi%noco,nococonv,fi%input,fi%banddos,fi%cell,fi%atoms,enpara,stars,fi%vacuum,&
449 : sphhar,fi%sym,vTot, cdnvalJob,overallDen,regCharges,dos,vacdos,results,moments,&
450 0 : fi%gfinp,fi%hub1inp)
451 : END DO
452 :
453 : CALL cdncore(fmpi, fi%input,fi%vacuum,fi%noco,nococonv,fi%sym,&
454 0 : stars,fi%cell,sphhar,fi%atoms,vTot,overallDen,moments,results)
455 0 : IF (fmpi%irank.EQ.0) THEN
456 : CALL qfix(fmpi,stars,nococonv,fi%atoms,fi%sym,fi%vacuum,sphhar,fi%input,fi%cell, overallDen,&
457 0 : fi%noco%l_noco,.TRUE.,l_par=.FALSE.,force_fix=.TRUE.,fix=fix)
458 : END IF
459 0 : CALL overallDen%distribute(fmpi%mpi_comm)
460 :
461 : ! Calculate Coulomb potential for overall density (+including external potential)
462 0 : CALL overallDen%sum_both_spin()!workden)
463 0 : CALL overallVCoul%resetPotDen()
464 0 : ALLOCATE(overallVCoul%pw_w(size(overallVCoul%pw,1),size(overallVCoul%pw,2)))
465 0 : overallVCoul%pw_w(:,:) = 0.0
466 0 : sigma_loc = cmplx(0.0,0.0)
467 0 : CALL vgen_coulomb(1,fmpi, fi%input,fi%field,fi%vacuum,fi%sym,fi%juphon,stars,fi%cell,sphhar,fi%atoms,.FALSE.,overallDen,overallVCoul,sigma_loc)
468 0 : CALL convol(stars,overallVCoul%pw_w(:,1),overallVCoul%pw(:,1)) ! Is there a problem with a second spin?!
469 0 : CALL overallVCoul%distribute(fmpi%mpi_comm)
470 :
471 0 : overallVCoulSSDen = 0.0
472 0 : DO jspin = 1, fi%input%jspins
473 0 : jsp = MERGE(1,jspin,fi%noco%l_noco)
474 0 : DO ikpt = 1, fi%kpts%nkpt
475 0 : DO iBand = 1, highestState(ikpt,jsp)
476 : ! Read the single-state density from disc
477 0 : filename = ''
478 0 : WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
479 0 : IF (fmpi%irank.EQ.0) THEN
480 : CALL readDensity(stars,fi%noco,fi%vacuum,fi%atoms,fi%cell,sphhar,fi%input,fi%sym, CDN_ARCHIVE_TYPE_CDN_const,&
481 0 : CDN_input_DEN_const,0,fermiEnergyTemp,tempDistance,l_qfix,singleStateDen,inFilename=TRIM(ADJUSTL(filename)))
482 0 : CALL singleStateDen%sum_both_spin()!workden)
483 : END IF
484 0 : CALL singleStateDen%distribute(fmpi%mpi_comm)
485 :
486 : ! For each state calculate integral over Coulomb potential times single state density
487 0 : potDenInt = 0.0
488 0 : CALL int_nv(1,stars,fi%vacuum,fi%atoms,sphhar,fi%cell,fi%sym,fi%input, vCoul,singleStateDen,potDenInt) ! Is there a problem with a second spin?!
489 0 : overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
490 : END DO
491 : END DO
492 : END DO
493 :
494 : ! Construct exchange matrix diagonal
495 :
496 0 : exDiag = 0.0
497 0 : DO jspin = 1, fi%input%jspins
498 0 : jsp = MERGE(1,jspin,fi%noco%l_noco)
499 : ! remove weights(wtkpt) in w_iks
500 0 : wl_iks = 0.0
501 0 : DO ikptf = 1, fi%kpts%nkptf
502 0 : ikpt = fi%kpts%bkp(ikptf)
503 0 : DO iBand=1, results%neig(ikpt,jsp)
504 0 : wl_iks(iBand,ikptf) = results%w_iks(iBand,ikpt,jspin) / (fi%kpts%wtkpt(ikpt))!*fi%kpts%nkptf) Last term to be included after applying functional
505 0 : wl_iks(iBand,ikptf) = sqrt(wl_iks(iBand,ikptf)) ! For Müller functional
506 0 : wl_iks(iBand,ikptf) = wl_iks(iBand,ikptf) / fi%kpts%nkptf ! This is the last part of two lines above
507 : END DO
508 : END DO
509 :
510 0 : IF(ALLOCATED(eig_irr)) DEALLOCATE (eig_irr)
511 0 : IF(ALLOCATED(hybdat%pntgptd)) DEALLOCATE (hybdat%pntgptd)
512 0 : IF(ALLOCATED(hybdat%pntgpt)) DEALLOCATE (hybdat%pntgpt)
513 0 : IF(ALLOCATED(hybdat%prodm)) DEALLOCATE (hybdat%prodm)
514 0 : IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE (hybdat%nindxp1)
515 :
516 0 : results%neig(:,:) = neigTemp(:,:)
517 :
518 : CALL HF_setup(mpdata,fi,fmpi,nococonv,results,jspin,enpara,&
519 0 : hybdat,vTot%mt(:,0,:,:),eig_irr)
520 :
521 0 : results%neig(:,:) = highestState(:,:) + 1
522 :
523 0 : DO ikpt = 1, fi%kpts%nkpt
524 :
525 0 : CALL lapw%init(fi%input,fi%noco,nococonv,fi%kpts,fi%atoms,fi%sym,ikpt,fi%cell)
526 :
527 : nbasfcn = 0
528 0 : IF(fi%noco%l_noco) then
529 0 : nbasfcn = lapw%nv(1) + lapw%nv(2) + 2*fi%atoms%nlotot
530 : ELSE
531 0 : nbasfcn = lapw%nv(1) + fi%atoms%nlotot
532 : END IF
533 :
534 0 : parent = 0
535 0 : CALL zMat%init(fi%sym%invs,nbasfcn,fi%input%neig)
536 :
537 0 : if(ikpt /= fi%kpts%bkp(ikpt)) call juDFT_error("We should be reading the parent z-mat here!")
538 0 : call read_z(fi%atoms, fi%cell, hybdat, fi%kpts, fi%sym, fi%noco, nococonv, fi%input, ikpt, jsp, zMat, c_phase=c_phase)
539 0 : allocate(cmt_nk(hybdat%nbands(ikpt,jsp), hybdat%maxlmindx, fi%atoms%nat), stat=ierr)
540 0 : if(ierr /= 0) call judft_error("can't allocate cmt_nk")
541 : call calc_cmt(fi%atoms, fi%cell, fi%input, fi%noco, nococonv, fi%hybinp, hybdat, mpdata, fi%kpts, &
542 0 : fi%sym, zMat, jsp, ikpt, c_phase, cmt_nk)
543 :
544 0 : ALLOCATE (indx_sest(hybdat%nbands(ikpt,jsp), hybdat%nbands(ikpt,jsp)))
545 0 : indx_sest = 0
546 :
547 0 : call symm_hf_init(fi,ikpt,nsymop,rrot,psym)
548 : call symm_hf(fi,ikpt,hybdat,results,work_pack%k_packs(ikpt)%submpi, eig_irr,mpdata, c_phase,&
549 0 : rrot,nsymop,psym,n_q,parent,nsest,indx_sest, jsp)
550 :
551 0 : exMat%l_real=fi%sym%invs
552 : CALL exchange_valence_hf(work_pack%k_packs(ikpt),fi,fmpi,zMat, mpdata,jspin,hybdat,lapw,&
553 : eig_irr,results,n_q,wl_iks,xcpot,nococonv,stars,nsest,indx_sest,&
554 0 : cmt_nk, exMat)
555 0 : deallocate(cmt_nk)
556 : CALL exchange_vccv1(ikpt,fi, mpdata,hybdat,jspin,lapw,glob_mpi,nsymop,nsest,indx_sest,&
557 0 : 1.0,results,cmt_nk,exMat)
558 :
559 0 : DEALLOCATE(indx_sest)
560 :
561 : !Start of workaround for increased functionality of fi%symmetrizeh (call it))
562 :
563 0 : nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*fi%atoms%nlotot,lapw%nv(1)+fi%atoms%nlotot,fi%noco%l_noco)
564 :
565 0 : CALL read_eig(hybdat%eig_id,ikpt,jspin, smat=olap)
566 :
567 0 : zMat%matsize2 = hybdat%nbands(ikpt,jsp) ! reduce "visible matsize" for the following computations
568 :
569 0 : CALL olap%multiply(zMat,trafo)
570 :
571 0 : CALL invtrafo%alloc(olap%l_real,hybdat%nbands(ikpt,jsp),nbasfcn)
572 0 : CALL trafo%TRANSPOSE(invtrafo)
573 :
574 0 : DO iBand = 1, hybdat%nbands(ikpt,jsp)
575 0 : DO jBand = 1, iBand-1
576 0 : IF (exMat%l_real) THEN
577 0 : exMat%data_r(iBand,jBand)=exMat%data_r(jBand,iBand)
578 : ELSE
579 0 : exMat%data_c(iBand,jBand)=conjg(exMat%data_c(jBand,iBand))
580 : END IF
581 : END DO
582 : END DO
583 :
584 0 : CALL exMat%multiply(invtrafo,tmpMat)
585 0 : CALL trafo%multiply(tmpMat,exMatLAPW)
586 :
587 0 : call symmetrizeh(fi%atoms,fi%kpts%bkf(:,ikpt),jspin,lapw,fi%sym,fi%cell,nsymop,psym,exMatLAPW)
588 :
589 0 : IF (.NOT.exMatLAPW%l_real) exMatLAPW%data_c=conjg(exMatLAPW%data_c)
590 0 : zMat%matsize1=MIN(zMat%matsize1,exMatLAPW%matsize2)
591 :
592 0 : CALL exMatLAPW%multiply(zMat,tmpMat)
593 :
594 0 : DO iBand = 1, highestState(ikpt,jsp)
595 0 : IF (zMat%l_real) THEN
596 0 : exDiag(iBand,ikpt,jspin) = dot_product(zMat%data_r(:zMat%matsize1,iband),tmpMat%data_r(:,iband))
597 : ELSE
598 0 : exDiag(iBand,ikpt,jspin) = REAL(dot_product(zMat%data_c(:zMat%matsize1,iband),tmpMat%data_c(:,iband)))
599 : END IF
600 : END DO
601 :
602 : !End of workaround for increased functionality of fi%symmetrizeh (call it))
603 :
604 : ! DO iBand = 1, highestState(ikpt,jsp)
605 : ! IF (exMat%l_real) THEN
606 : ! exDiag(iBand,ikpt,jspin) = exMat%data_r(iBand,iBand)
607 : ! ELSE
608 : ! exDiag(iBand,ikpt,jspin) = REAL(exMat%data_c(iBand,iBand))
609 : ! END IF
610 : ! END DO
611 : END DO
612 : END DO
613 :
614 : ! Calculate total energy derivative with respect to occupations (dEdOcc)
615 :
616 0 : occSum = 0.0
617 0 : gradSum = 0.0
618 0 : DO ispin = 1, fi%input%jspins
619 0 : isp = MERGE(1,ispin,fi%noco%l_noco)
620 : ! CALL cdnvalJob%init(fmpi,fi%input,fi%kpts,fi%noco,results,isp,fi%banddos=fi%banddos)
621 0 : DO ikpt = 1, fi%kpts%nkpt
622 0 : DO iBand = 1, highestState(ikpt,isp)
623 0 : occStateI = results%w_iks(iBand,ikpt,isp) / (fi%kpts%wtkpt(ikpt))!*fi%kpts%nkptf)
624 0 : occStateI = MAX(occStateI,minOcc)
625 :
626 0 : occStateI = MIN(occStateI,1.0-minOcc)
627 :
628 : ! IF(occStateI.LT.1.0e-7) occStateI = 5.0e-4 ! This is preliminary. I have to discuss what do do here.
629 : ! occStateI = cdnvalJob%weights(iBand,ikpt)
630 0 : rdmftFunctionalValue = 0.5*0.5*SQRT(1.0/occStateI) ! for Müller functional derivative
631 :
632 :
633 : !!! Test start
634 : ! occStateI = MIN(occStateI,1.0-minOcc)
635 : ! rdmftFunctionalValue = 0.5 * 0.5*pi_const*COS(0.5*pi_const*occStateI)
636 :
637 : ! rdmftFunctionalValue = ASIN(SQRT(occStateI)) * 2.0 / pi_const
638 : ! rdmftFunctionalValue = (SIN(rdmftFunctionalValue*pi_const/2.0)) * (COS(rdmftFunctionalValue*pi_const/2.0))**2.0 * pi_const
639 : !!! Test 2:
640 : ! occStateI = MIN(occStateI,1.0-minOcc)
641 : ! rdmftFunctionalValue = ASIN(SQRT(occStateI)) * 2.0 / pi_const
642 : ! rdmftFunctionalValue = COS(rdmftFunctionalValue*pi_const/2.0) * pi_const / 4.0
643 : !!! Test end
644 :
645 0 : occSum = occSum + results%w_iks(iBand,ikpt,isp)
646 :
647 0 : exchangeTerm = - rdmftFunctionalValue * exDiag(iBand,ikpt,isp) * fi%kpts%wtkpt(ikpt) * spinDegenFac !*fi%kpts%nkptf
648 :
649 : dEdOcc(iBand,ikpt,isp) = +((spinDegenFac * results%eig(iBand,ikpt,isp)) - vTotSSDen(iBand,ikpt,isp) + &
650 : overallVCoulSSDen(iBand,ikpt,isp) - &
651 0 : zintn_rSSDen(iBand,ikpt,isp) + vmdSSDen(iBand,ikpt,isp) + exchangeTerm )
652 :
653 0 : theta = ASIN(SQRT(occStateI))! * 2.0 / pi_const
654 0 : dEdOcc(iBand,ikpt,isp) = 2.0 * sin(theta) * cos(theta) * dEdOcc(iBand,ikpt,isp)
655 : ! dEdOcc(iBand,ikpt,isp) = dEdOcc(iBand,ikpt,isp) + 2.0 * COS(theta) * exDiag(iBand,ikpt,isp) * fi%kpts%wtkpt(ikpt) * spinDegenFac
656 :
657 0 : WRITE(*,*) 'ENERGY GRADIENT CONTRIBUTIONS'
658 0 : WRITE(*,*) 'ispin, ikpt, iBand, weight', ispin, ikpt, iBand, results%w_iks(iBand,ikpt,isp)
659 0 : WRITE(*,*) 'results%eig(iBand,ikpt,isp)', results%eig(iBand,ikpt,isp)
660 0 : WRITE(*,*) 'vTotSSDen(iBand,ikpt,isp)', vTotSSDen(iBand,ikpt,isp)
661 0 : WRITE(*,*) 'overallVCoulSSDen(iBand,ikpt,isp)', overallVCoulSSDen(iBand,ikpt,isp)
662 0 : WRITE(*,*) 'zintn_rSSDen(iBand,ikpt,isp)', zintn_rSSDen(iBand,ikpt,isp)
663 0 : WRITE(*,*) 'vmdSSDen(iBand,ikpt,isp)', vmdSSDen(iBand,ikpt,isp)
664 0 : WRITE(*,*) 'exchangeTerm', exchangeTerm
665 0 : WRITE(*,*) 'exDiag(iBand,ikpt,isp)', exDiag(iBand,ikpt,isp)
666 0 : WRITE(*,*) 'rdmftFunctionalValue', rdmftFunctionalValue
667 :
668 0 : gradSum = gradSum + dEdOcc(iBand,ikpt,isp) !* results%w_iks(iBand,ikpt,isp))
669 : END DO
670 : END DO
671 : END DO
672 0 : lagrangeMultiplier = -gradSum / numStates ! occSum !(fi%input%zelec/(2.0/REAL(fi%input%jspins)))
673 :
674 0 : WRITE(*,*) 'lagrangeMultiplier: ', lagrangeMultiplier
675 :
676 : ! Optimize occupation numbers
677 :
678 0 : ALLOCATE (gradient(numStates+1))
679 0 : ALLOCATE (parameters(numStates+1))
680 0 : ALLOCATE (paramGroup(numStates+1))
681 0 : ALLOCATE (equalityLinCombi(numStates+1))
682 0 : ALLOCATE (minConstraints(numStates+1))
683 0 : ALLOCATE (maxConstraints(numStates+1))
684 0 : ALLOCATE (enabledConstraints(numStates+1))
685 0 : enabledConstraints = .FALSE.
686 0 : enabledConstraints(1:numStates) = .TRUE.
687 0 : minConstraints(:) = 0.0
688 0 : maxConstraints(:) = 1.0
689 0 : equalityLinCombi = 0.0
690 0 : gradient = 0.0
691 0 : parameters = 0.0
692 0 : iState = 0
693 0 : DO ispin = 1, fi%input%jspins
694 0 : isp = MERGE(1,ispin,fi%noco%l_noco)
695 0 : DO ikpt = 1, fi%kpts%nkpt
696 0 : DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
697 0 : iState = iState + 1
698 0 : occStateI = results%w_iks(iBand,ikpt,isp) / fi%kpts%wtkpt(ikpt)
699 :
700 0 : occStateI = MAX(occStateI,minOcc)
701 0 : occStateI = MIN(occStateI,1.0-minOcc)
702 :
703 0 : theta = ASIN(SQRT(occStateI))! * 2.0 / pi_const
704 :
705 0 : WRITE(7865,'(i7,4f15.10)') iState, occStateI, theta, sin(theta), cos(theta)
706 :
707 : ! occStateI = MAX(occStateI,minOcc)
708 0 : equalityLinCombi(iState) = fi%kpts%wtkpt(ikpt)
709 :
710 : ! dEdOcc(iBand,ikpt,isp) = dEdOcc(iBand,ikpt,isp) + lagrangeMultiplier
711 :
712 : ! dEdOcc(iBand,ikpt,isp) = 2.0 * sin(theta) * cos(theta) * dEdOcc(iBand,ikpt,isp)
713 :
714 0 : gradient(iState) = dEdOcc(iBand,ikpt,isp) + lagrangeMultiplier
715 :
716 0 : gradient(numStates+1) = gradient(numStates+1) + occStateI * fi%kpts%wtkpt(ikpt)
717 : ! parameters(iState) = theta !occStateI
718 0 : parameters(iState) = occStateI
719 0 : IF(iBand.EQ.1) THEN
720 0 : IF(iState.EQ.1) THEN
721 0 : paramGroup(iState) = 1
722 : ELSE
723 0 : paramGroup(iState) = paramGroup(iState-1) + 1
724 : END IF
725 : ELSE
726 0 : IF ((results%eig(iBand,ikpt,jsp) - results%eig(iBand-1,ikpt,jsp)).GT.degenEps) THEN
727 0 : paramGroup(iState) = paramGroup(iState-1) + 1
728 : ELSE
729 0 : paramGroup(iState) = paramGroup(iState-1)
730 : END IF
731 : END IF
732 : END DO
733 : END DO
734 : END DO
735 0 : equalityCriterion = fi%input%zelec/(2.0/REAL(fi%input%jspins))
736 0 : gradient(numStates+1) = gradient(numStates+1) - equalityCriterion ! This should actually always be 0.0
737 0 : parameters(numStates+1) = lagrangeMultiplier
738 0 : paramGroup(numStates+1) = paramGroup(numStates) + 1
739 :
740 0 : currentGroupEnd = 0
741 0 : lastGroupEnd = 0
742 0 : DO iState = 2, numStates + 1
743 0 : IF (paramGroup(iState).NE.paramGroup(iState-1)) THEN
744 0 : currentGroupEnd = iState - 1
745 0 : averageParam = 0.0
746 0 : averageGrad = 0.0
747 0 : DO jState = lastGroupEnd + 1, currentGroupEnd
748 0 : averageParam = averageParam + parameters(jState)
749 0 : averageGrad = averageGrad + gradient(jState)
750 : END DO
751 0 : averageParam = averageParam / (currentGroupEnd - lastGroupEnd)
752 0 : averageGrad = averageGrad / (currentGroupEnd - lastGroupEnd)
753 0 : DO jState = lastGroupEnd + 1, currentGroupEnd
754 0 : temp = ABS(parameters(jState) - averageParam) / (ABS(averageParam) + degenEps)
755 0 : IF (temp.GT.degenEps) THEN
756 0 : WRITE(*,'(a,i0,a,f15.10,a,f15.10)') 'iState: ', jState, 'average: ', averageParam, 'parameter: ', parameters(jState)
757 0 : CALL juDFT_error('parameter difference in paramGroup is very large (rdmft-1)')
758 : END IF
759 0 : temp = ABS(gradient(jState) - averageGrad) / (ABS(averageGrad) + degenEps)
760 0 : IF (temp.GT.degenEps) THEN
761 0 : WRITE(*,'(a,i0,a,f15.10,a,f15.10)') 'iState: ', jState, 'average: ', averageGrad, 'gradient: ', gradient(jState)
762 0 : CALL juDFT_error('Gradient difference in paramGroup is very large (rdmft-1)')
763 : END IF
764 0 : parameters(jState) = averageParam
765 0 : gradient(jState) = averageGrad
766 : END DO
767 : lastGroupEnd = currentGroupEnd
768 : END IF
769 : END DO
770 :
771 0 : WRITE(*,*) 'gradient(numStates+1): ', gradient(numStates+1)
772 :
773 : ! mixParam = 0.01 / MAXVAL(ABS(gradient(:numStates)))
774 : ! mixParam = MIN(0.0002,mixParam)
775 0 : mixParam = 0.001
776 0 : WRITE(*,*) 'mixParam: ', mixParam
777 :
778 : CALL bfgs_b2(numStates+1,gradient,lastGradient,minConstraints,maxConstraints,enabledConstraints,parameters,&
779 : lastParameters,equalityLinCombi,equalityCriterion,maxHistoryLength,paramCorrections,&
780 0 : gradientCorrections,iStep,mixParam,converged,convCrit)
781 :
782 :
783 0 : currentGroupEnd = 0
784 0 : lastGroupEnd = 0
785 0 : DO iState = 2, numStates + 1
786 0 : IF (paramGroup(iState).NE.paramGroup(iState-1)) THEN
787 0 : currentGroupEnd = iState - 1
788 0 : averageParam = 0.0
789 0 : averageGrad = 0.0
790 0 : DO jState = lastGroupEnd + 1, currentGroupEnd
791 0 : averageParam = averageParam + parameters(jState)
792 0 : averageGrad = averageGrad + gradient(jState)
793 : END DO
794 0 : averageParam = averageParam / (currentGroupEnd - lastGroupEnd)
795 0 : averageGrad = averageGrad / (currentGroupEnd - lastGroupEnd)
796 0 : DO jState = lastGroupEnd + 1, currentGroupEnd
797 0 : temp = ABS(parameters(jState) - averageParam) / (ABS(averageParam) + degenEps)
798 0 : IF (temp.GT.degenEps) THEN
799 0 : WRITE(*,'(a,i0,a,f15.10,a,f15.10)') 'iState: ', jState, 'average: ', averageParam, 'parameter: ', parameters(jState)
800 0 : CALL juDFT_error('parameter difference in paramGroup is very large (rdmft-2)')
801 : END IF
802 0 : temp = ABS(gradient(jState) - averageGrad) / (ABS(averageGrad) + degenEps)
803 0 : IF (temp.GT.degenEps) THEN
804 0 : WRITE(*,'(a,i0,a,f15.10,a,f15.10)') 'iState: ', jState, 'average: ', averageGrad, 'gradient: ', gradient(jState)
805 0 : CALL juDFT_error ('Gradient difference in paramGroup is very large (rdmft-2)')
806 : END IF
807 0 : parameters(jState) = averageParam
808 0 : gradient(jState) = averageGrad
809 : END DO
810 : lastGroupEnd = currentGroupEnd
811 : END IF
812 : END DO
813 :
814 :
815 0 : WRITE(3555,*) 'Occupation numbers:'
816 0 : iState = 0
817 0 : DO ispin = 1, fi%input%jspins
818 0 : isp = MERGE(1,ispin,fi%noco%l_noco)
819 0 : DO ikpt = 1, fi%kpts%nkpt
820 0 : DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
821 0 : iState = iState + 1
822 : ! parameters(iState) = (SIN(parameters(iState)*0.5*pi_const))**2.0
823 0 : WRITE(3555,'(3i7,f15.10)') iBand, ikpt, isp, parameters(iState)
824 0 : results%w_iks(iBand,ikpt,isp) = parameters(iState) * fi%kpts%wtkpt(ikpt)
825 : ! results%w_iks(iBand,ikpt,isp) = MERGE(parameters(iState) * fi%kpts%wtkpt(ikpt),0.0,parameters(iState).GT.minOcc)
826 : END DO
827 : END DO
828 : END DO
829 :
830 0 : WRITE(3555,'(a,f15.10)') 'total occupation: ', SUM(parameters(:numStates))
831 :
832 0 : DEALLOCATE (enabledConstraints,maxConstraints,minConstraints)
833 0 : DEALLOCATE (parameters,gradient,paramGroup,equalityLinCombi)
834 :
835 :
836 : END DO ! WHILE (.NOT.converged)
837 :
838 : ! Only mix a part of the newly determined occupation numbers into the occupation numbers from the
839 : ! previous iteration.
840 : ! WRITE(*,*) 'Test: mixing in only a part of newly determined occupation numbers!'
841 : ! DO ispin = 1, fi%input%jspins
842 : ! isp = MERGE(1,ispin,fi%noco%l_noco)
843 : ! DO ikpt = 1, fi%kpts%nkpt
844 : ! DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
845 : ! results%w_iks(iBand,ikpt,isp) = (1.0-occMixParam) * results%w_iksRDMFT(iBand,ikpt,isp) + &
846 : ! occMixParam * results%w_iks(iBand,ikpt,isp)
847 : ! END DO
848 : ! END DO
849 : ! END DO
850 :
851 :
852 0 : WRITE(2503,*) 'convIter: ', convIter
853 :
854 0 : WRITE(*,*) 'RDMFT: convergence loop end'
855 :
856 0 : hybdat%l_calhf = .FALSE.
857 :
858 : ! Calculate final overall density
859 :
860 : !I think we need most of cdngen at this place so I just use cdngen
861 0 : CALL outDen%resetPotDen()
862 : CALL cdngen(eig_id,fmpi,fi%input,fi%banddos,fi%sliceplot,fi%vacuum,fi%kpts,fi%atoms,sphhar,stars,fi%sym,fi%juphon,fi%gfinp,fi%hub1inp,&
863 0 : enpara,fi%cell,fi%noco,nococonv,vTot,results, fi%corespecinput,archiveType,xcpot,outDen, EnergyDen)
864 :
865 : ! Calculate RDMFT energy
866 0 : rdmftEnergy = 0.0
867 0 : DO ispin = 1, fi%input%jspins
868 0 : isp = MERGE(1,ispin,fi%noco%l_noco)
869 : ! CALL cdnvalJob%init(fmpi,fi%input,fi%kpts,fi%noco,results,isp,fi%banddos=fi%banddos)
870 0 : DO ikpt = 1, fi%kpts%nkpt
871 0 : DO iBand = 1, highestState(ikpt,isp)
872 0 : occStateI = results%w_iks(iBand,ikpt,isp) / (fi%kpts%wtkpt(ikpt))!*fi%kpts%nkptf)
873 0 : rdmftFunctionalValue = 0.5*SQRT(occStateI) ! for Müller functional
874 :
875 0 : exchangeTerm = -rdmftFunctionalValue * exDiag(iBand,ikpt,isp) * fi%kpts%wtkpt(ikpt) * spinDegenFac !*fi%kpts%nkptf
876 :
877 : rdmftEnergy = rdmftEnergy + exchangeTerm + &
878 : occStateI * ((spinDegenFac*results%eig(iBand,ikpt,isp)) - vTotSSDen(iBand,ikpt,isp) + &
879 0 : 0.5*overallVCoulSSDen(iBand,ikpt,isp))
880 :
881 0 : WRITE(2505,*) 'ENERGY CONTRIBUTIONS'
882 0 : WRITE(2505,*) 'ispin, ikpt, iBand, weight', ispin, ikpt, iBand, results%w_iks(iBand,ikpt,isp)
883 0 : WRITE(2505,*) 'results%eig(iBand,ikpt,isp)', occStateI * spinDegenFac * results%eig(iBand,ikpt,isp)
884 0 : WRITE(2505,*) 'vTotSSDen(iBand,ikpt,isp)', - occStateI * vTotSSDen(iBand,ikpt,isp)
885 0 : WRITE(2505,*) 'overallVCoulSSDen(iBand,ikpt,isp)', occStateI * overallVCoulSSDen(iBand,ikpt,isp)
886 0 : WRITE(2505,*) 'exchangeTerm', exchangeTerm
887 0 : WRITE(2505,*) 'exDiag(iBand,ikpt,isp)', exDiag(iBand,ikpt,isp)
888 0 : WRITE(2505,*) 'rdmftFunctionalValue', rdmftFunctionalValue
889 :
890 : END DO
891 : END DO
892 : END DO
893 :
894 0 : results%w_iksRDMFT(:,:,:) = results%w_iks(:,:,:)
895 0 : results%neig(:,:) = neigTemp(:,:)
896 :
897 :
898 : ! Madelung term (taken from totale):
899 :
900 0 : mt=0.0
901 0 : DO iType = 1, fi%atoms%ntype
902 0 : DO iGrid = 1, fi%atoms%jri(iType)
903 0 : mt(iGrid,iType) = outDen%mt(iGrid,0,iType,1) + outDen%mt(iGrid,0,iType,fi%input%jspins)
904 : END DO
905 : END DO
906 0 : IF (fi%input%jspins.EQ.1) mt=mt/2.0 !we just added the same value twice
907 :
908 0 : DO iType = 1, fi%atoms%ntype
909 0 : DO iGrid = 1,fi%atoms%jri(iType)
910 0 : dpj(iGrid) = mt(iGrid,iType)/fi%atoms%rmsh(iGrid,iType)
911 : END DO
912 0 : CALL intgr3(dpj,fi%atoms%rmsh(1,iType),fi%atoms%dx(iType),fi%atoms%jri(iType),rhs)
913 :
914 0 : zintn_r(iType) = fi%atoms%neq(iType)*fi%atoms%zatom(iType)*sfp_const*rhs/2.
915 0 : rdmftEnergy = rdmftEnergy - zintn_r(iType)
916 :
917 0 : CALL intgr3(mt(1,iType),fi%atoms%rmsh(1,iType),fi%atoms%dx(iType),fi%atoms%jri(iType),totz)
918 :
919 0 : vmd(iType) = fi%atoms%rmt(iType)*vCoul%mt(fi%atoms%jri(iType),0,iType,1)/sfp_const + fi%atoms%zatom(iType) - totz*sfp_const
920 0 : vmd(iType) = -fi%atoms%neq(iType)*fi%atoms%zatom(iType)*vmd(iType)/ (2.0*fi%atoms%rmt(iType))
921 :
922 0 : rdmftEnergy = rdmftEnergy + vmd(iType)
923 :
924 0 : WRITE(2505,*) '======================================='
925 0 : WRITE(2505,*) 'iType: ', iType
926 0 : WRITE(2505,*) 'zintn_r(iType): ', zintn_r(iType)
927 0 : WRITE(2505,*) 'vmd(iType): ', vmd(iType)
928 0 : WRITE(2505,*) '======================================='
929 :
930 : END DO
931 :
932 : ! Output
933 0 : WRITE(oUnit,'(a,f20.10,a)') 'RDMFT energy: ', rdmftEnergy, ' Htr'
934 0 : CALL openXMLElementPoly('rdmft',(/'energy'/),(/rdmftEnergy/))
935 0 : DO ispin = 1, fi%input%jspins
936 0 : isp = MERGE(1,ispin,fi%noco%l_noco)
937 0 : DO ikpt = 1, fi%kpts%nkpt
938 0 : CALL openXMLElementPoly('occupations',(/'spin ','kpoint'/),(/ispin,ikpt/))
939 0 : DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
940 0 : occStateI = results%w_iks(iBand,ikpt,isp) / (fi%kpts%wtkpt(ikpt))!*fi%kpts%nkptf)
941 0 : WRITE(attributes(1),'(i0)') iBand
942 0 : WRITE(attributes(2),'(f18.10)') results%eig(iBand,ikpt,isp)
943 0 : WRITE(attributes(3),'(f18.10)') occStateI
944 0 : CALL writeXMLElement('state',(/'index ','energy ','occupation'/),attributes)
945 : END DO
946 0 : CALL closeXMLElement('occupations')
947 : END DO
948 : END DO
949 0 : CALL closeXMLElement('rdmft')
950 :
951 0 : WRITE(2505,*) '======================================='
952 0 : WRITE(2505,*) 'convIter: ', convIter
953 0 : WRITE(2505,'(a,f20.10,a)') 'RDMFT energy: ', rdmftEnergy, ' Htr'
954 0 : WRITE(2505,*) '======================================='
955 0 : WRITE(2505,*) '======================================='
956 0 : WRITE(2505,*) '======================================='
957 :
958 : #endif
959 0 : END SUBROUTINE rdmft
960 :
961 : END MODULE m_rdmft
|
