rdmft.F90 42.8 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
!--------------------------------------------------------------------------------
! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_rdmft

CONTAINS

11
SUBROUTINE rdmft(eig_id,fmpi,fi,enpara,stars,&
12 13
                 sphhar,vTot,vCoul,nococonv,xcpot,mpdata,hybdat,&
                 results,archiveType,outDen)
14
   use m_types_vacdos
Matthias Redies's avatar
Matthias Redies committed
15
   use m_work_package
16 17
   USE m_types
   USE m_juDFT
18
   USE m_constants
Gregor Michalicek's avatar
Gregor Michalicek committed
19
   USE m_intgr, ONLY : intgr3
Gregor Michalicek's avatar
Gregor Michalicek committed
20
   USE m_eig66_io
21
#ifndef CPP_OLDINTEL
22
   USE m_cdnval
Gregor Michalicek's avatar
Gregor Michalicek committed
23
   USE m_cdngen
24
   USE m_cdn_io
25 26
   USE m_cdncore
   USE m_qfix
27 28 29
   USE m_vgen_coulomb
   USE m_convol
   USE m_intnv
Gregor Michalicek's avatar
Gregor Michalicek committed
30 31 32 33 34

   USE m_mixedbasis
   USE m_coulombmatrix
   USE m_hf_init
   USE m_hf_setup
Matthias Redies's avatar
Matthias Redies committed
35
   USE m_io_hybinp
Gregor Michalicek's avatar
Gregor Michalicek committed
36 37 38
   USE m_symm_hf
   USE m_exchange_valence_hf
   USE m_exchange_core
Gregor Michalicek's avatar
Gregor Michalicek committed
39
   USE m_symmetrizeh
Gregor Michalicek's avatar
Gregor Michalicek committed
40
   USE m_bfgs_b2
41
   USE m_xmlOutput
42
   USE m_types_dos
Matthias Redies's avatar
Matthias Redies committed
43
   use m_calc_cmt
Gregor Michalicek's avatar
Gregor Michalicek committed
44

45 46
#ifdef CPP_MPI
   USE m_mpi_bc_potden
47
#endif
48
#endif
49 50 51

   IMPLICIT NONE

52
   TYPE(t_mpi),           INTENT(IN)    :: fmpi
53
   type(t_fleurinput), intent(in)    :: fi
Gregor Michalicek's avatar
Gregor Michalicek committed
54
   TYPE(t_enpara),        INTENT(INOUT) :: enpara
55 56
   TYPE(t_stars),         INTENT(IN)    :: stars
   TYPE(t_sphhar),        INTENT(IN)    :: sphhar
57
   TYPE(t_potden),        INTENT(INOUT) :: vTot
Gregor Michalicek's avatar
Gregor Michalicek committed
58
   TYPE(t_potden),        INTENT(INOUT) :: vCoul
59 60
   TYPE(t_nococonv),      INTENT(INOUT) :: nococonv
   TYPE(t_xcpot_inbuild), INTENT(IN)    :: xcpot
61 62
   TYPE(t_mpdata),        intent(inout) :: mpdata
   TYPE(t_hybdat),        INTENT(INOUT) :: hybdat
63
   TYPE(t_results),       INTENT(INOUT) :: results
Gregor Michalicek's avatar
Gregor Michalicek committed
64
   TYPE(t_potden),        INTENT(INOUT) :: outDen
65 66

   INTEGER,               INTENT(IN)    :: eig_id
Gregor Michalicek's avatar
Gregor Michalicek committed
67
   INTEGER,               INTENT(IN)    :: archiveType
Matthias Redies's avatar
merge  
Matthias Redies committed
68
   TYPE(t_potden)                       :: EnergyDen
69

70
#ifndef CPP_OLDINTEL
71
   TYPE(t_cdnvalJob)                    :: cdnvalJob
Gregor Michalicek's avatar
Gregor Michalicek committed
72
   TYPE(t_potden)                       :: singleStateDen, overallDen, overallVCoul, vTotTemp
73 74
   TYPE(t_regionCharges)                :: regCharges
   TYPE(t_dos)                          :: dos
75
   TYPE(t_vacdos)                       :: vacdos
76
   TYPE(t_moments)                      :: moments
Gregor Michalicek's avatar
Gregor Michalicek committed
77
   TYPE(t_mat)                          :: exMat, zMat, olap, trafo, invtrafo, tmpMat, exMatLAPW
Gregor Michalicek's avatar
Gregor Michalicek committed
78
   TYPE(t_lapw)                         :: lapw
Matthias Redies's avatar
Matthias Redies committed
79
   type(t_work_package)                 :: work_pack
80
   INTEGER                              :: ikpt, ikpt_i, iBand, jkpt, jBand, iAtom, na, itype, lh, iGrid
Gregor Michalicek's avatar
Gregor Michalicek committed
81
   INTEGER                              :: jspin, jspmax, jsp, isp, ispin, nbasfcn, nbands
Matthias Redies's avatar
Matthias Redies committed
82
   INTEGER                              :: nsymop, ikptf, iterHF, ierr
83
   INTEGER                              :: iState, jState, iStep, numStates, numRelevantStates, convIter
Gregor Michalicek's avatar
Gregor Michalicek committed
84
   INTEGER                              :: maxHistoryLength
85
   INTEGER                              :: lastGroupEnd, currentGroupEnd
Gregor Michalicek's avatar
Gregor Michalicek committed
86 87
   REAL                                 :: fix, potDenInt, fermiEnergyTemp, spinDegenFac
   REAL                                 :: rdmftFunctionalValue, occStateI, gradSum
Gregor Michalicek's avatar
Gregor Michalicek committed
88
   REAL                                 :: exchangeTerm, lagrangeMultiplier, equalityCriterion
Gregor Michalicek's avatar
Gregor Michalicek committed
89 90
   REAL                                 :: mixParam, rdmftEnergy, occSum
   REAL                                 :: sumOcc, addCharge, subCharge, addChargeWeight, subChargeWeight
91 92
   REAL                                 :: rhs, totz, theta, temp
   REAL                                 :: averageParam, averageGrad
Gregor Michalicek's avatar
Gregor Michalicek committed
93
   REAL, PARAMETER                      :: degenEps = 0.00001
94 95
   REAL, PARAMETER                      :: convCrit = 5.0e-6, occMixParam = 0.5
   REAL, PARAMETER                      :: minOcc = 1.0e-13, minOccB = 1.0e-5
Gregor Michalicek's avatar
Gregor Michalicek committed
96
   LOGICAL                              :: converged, l_qfix, l_restart, l_zref
97
   CHARACTER(LEN=20)                    :: filename
98
   CHARACTER(LEN=20)                    :: attributes(3)
99

100 101 102 103 104 105
   INTEGER                              :: nsest(fi%input%neig) ! probably too large
   INTEGER                              :: rrot(3,3,fi%sym%nsym)
   INTEGER                              :: psym(fi%sym%nsym) ! Note: psym is only filled up to index nsymop
   INTEGER                              :: lowestState(fi%kpts%nkpt,fi%input%jspins)
   INTEGER                              :: highestState(fi%kpts%nkpt,fi%input%jspins)
   INTEGER                              :: neigTemp(fi%kpts%nkpt,fi%input%jspins)
Gregor Michalicek's avatar
Gregor Michalicek committed
106

107
   REAL                                 :: wl_iks(fi%input%neig,fi%kpts%nkptf)
Gregor Michalicek's avatar
Gregor Michalicek committed
108

109
   REAL                                 :: vmd(fi%atoms%ntype), zintn_r(fi%atoms%ntype), dpj(fi%atoms%jmtd), mt(fi%atoms%jmtd,fi%atoms%ntype)
Gregor Michalicek's avatar
Gregor Michalicek committed
110

111 112
   REAL, ALLOCATABLE                    :: overallVCoulSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vTotSSDen(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
113 114
   REAL, ALLOCATABLE                    :: dEdOcc(:,:,:)

Gregor Michalicek's avatar
Gregor Michalicek committed
115 116 117 118
   REAL, ALLOCATABLE                    :: zintn_rSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vmdSSDen(:,:,:)


Gregor Michalicek's avatar
Gregor Michalicek committed
119 120 121 122 123 124 125 126 127 128 129 130
   REAL, ALLOCATABLE                    :: exDiag(:,:,:)
   REAL, ALLOCATABLE                    :: eig_irr(:,:)

   REAL, ALLOCATABLE                    :: gradient(:), lastGradient(:)
   REAL, ALLOCATABLE                    :: parameters(:), lastParameters(:)
   REAL, ALLOCATABLE                    :: minConstraints(:)
   REAL, ALLOCATABLE                    :: maxConstraints(:)
   REAL, ALLOCATABLE                    :: gradientCorrections(:,:)
   REAL, ALLOCATABLE                    :: paramCorrections(:,:)
   REAL, ALLOCATABLE                    :: equalityLinCombi(:)

   REAL, ALLOCATABLE                    :: occupationVec(:)
131

132
   INTEGER, ALLOCATABLE                 :: paramGroup(:)
133
   INTEGER, ALLOCATABLE                 :: indx_sest(:,:)
134
   INTEGER, ALLOCATABLE                 :: parent(:)
Gregor Michalicek's avatar
Gregor Michalicek committed
135 136
   INTEGER, ALLOCATABLE                 :: pointer_EIBZ(:)
   INTEGER, ALLOCATABLE                 :: n_q(:)
Matthias Redies's avatar
Matthias Redies committed
137
   complex, allocatable                 :: cmt_nk(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
138 139

   LOGICAL, ALLOCATABLE                 :: enabledConstraints(:)
Matthias Redies's avatar
Matthias Redies committed
140
   type(t_hybmpi)    :: glob_mpi
Gregor Michalicek's avatar
Gregor Michalicek committed
141

142
   complex :: c_phase(fi%input%neig)
143

Gregor Michalicek's avatar
Gregor Michalicek committed
144
#endif
145

146
#ifndef CPP_OLDINTEL
Gregor Michalicek's avatar
Gregor Michalicek committed
147 148 149

   WRITE(*,*) 'entered RDMFT subroutine'

150
   ! General initializations
151
   mixParam = 0.001
Gregor Michalicek's avatar
Gregor Michalicek committed
152
   lagrangeMultiplier = 0.1 !results%ef
153
   spinDegenFac = 2.0 / fi%input%jspins ! This factor is used to compensate the missing second spin in non-spinpolarized calculations
Gregor Michalicek's avatar
Gregor Michalicek committed
154 155 156 157 158 159

   neigTemp(:,:) = results%neig(:,:)

   ! Determine which states have to be considered
   lowestState(:,:) = -1
   highestState(:,:) = -1
160 161 162
   DO jspin = 1, fi%input%jspins
      jsp = MERGE(1,jspin,fi%noco%l_noco)
      DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
163 164
         ! determine lowest state
         DO iBand = 1, results%neig(ikpt,jsp)
165
            IF((results%w_iks(iBand,ikpt,jspin) / fi%kpts%wtkpt(ikpt)).LT.(1.0-fi%input%rdmftOccEps)) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
166 167 168 169
               lowestState(ikpt,jspin) = iBand
               EXIT
            END IF
         END DO
170
         lowestState(ikpt,jspin) = lowestState(ikpt,jspin) - fi%input%rdmftStatesBelow
Gregor Michalicek's avatar
Gregor Michalicek committed
171 172 173 174 175 176 177 178 179 180
         DO iBand = lowestState(ikpt,jspin)-1, 1, -1
            IF((results%eig(iBand+1,ikpt,jsp) - results%eig(iBand,ikpt,jsp)).GT.degenEps) THEN
               EXIT
            END IF
            lowestState(ikpt,jspin) = iBand
         END DO
         lowestState(ikpt,jspin) = MAX(lowestState(ikpt,jspin),1)

         ! determine highest state
         DO iBand = results%neig(ikpt,jsp)-1, 1, -1
181
            IF((results%w_iks(iBand,ikpt,jspin) / fi%kpts%wtkpt(ikpt)).GT.(0.0+fi%input%rdmftOccEps)) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
182 183 184 185
               highestState(ikpt,jspin) = iBand
               EXIT
            END IF
         END DO
186
         highestState(ikpt,jspin) = highestState(ikpt,jspin) + fi%input%rdmftStatesAbove
Gregor Michalicek's avatar
Gregor Michalicek committed
187
         IF((results%neig(ikpt,jsp)-1).LT.highestState(ikpt,jspin)) THEN
188 189 190 191
            WRITE(oUnit,*) 'Error: Not enough states calculated:'
            WRITE(oUnit,*) 'ikpt, jsp: ', ikpt, jsp
            WRITE(oUnit,*) 'highestState(ikpt,jspin): ', highestState(ikpt,jspin)
            WRITE(oUnit,*) 'results%neig(ikpt,jsp): ', results%neig(ikpt,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
192 193 194 195 196 197 198 199 200
            CALL juDFT_error('Not enough states calculated', calledby = 'rdmft')
         END IF
         DO iBand = highestState(ikpt,jspin)+1, results%neig(ikpt,jsp)
            IF((results%eig(iBand,ikpt,jsp) - results%eig(iBand-1,ikpt,jsp)).GT.degenEps) THEN
               EXIT
            END IF
            highestState(ikpt,jspin) = iBand
         END DO
         IF(highestState(ikpt,jspin).EQ.results%neig(ikpt,jsp)) THEN
201 202 203 204 205
            WRITE(oUnit,*) 'Error: Highest state is degenerate:'
            WRITE(oUnit,*) 'ikpt, jsp: ', ikpt, jsp
            WRITE(oUnit,*) 'highestState(ikpt,jspin): ', highestState(ikpt,jspin)
            WRITE(oUnit,*) 'results%eig(highestState(ikpt,jspin),ikpt,jsp): ', results%eig(highestState(ikpt,jspin),ikpt,jsp)
            WRITE(oUnit,*) 'results%eig(highestState(ikpt,jspin)-1,ikpt,jsp): ', results%eig(highestState(ikpt,jspin)-1,ikpt,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
206 207 208 209 210
            CALL juDFT_error('Highest state is degenerate', calledby = 'rdmft')
         END IF
      END DO
   END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
211 212 213
   IF (ANY(results%w_iksRDMFT(:,:,:).NE.0.0)) THEN
      results%w_iks(:,:,:) = results%w_iksRDMFT(:,:,:)
   END IF
214
!   results%w_iksRDMFT(:,:,:) = results%w_iks(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
215

Gregor Michalicek's avatar
Gregor Michalicek committed
216
   ! Move occupations of relevant states well into allowed region
217
   numRelevantStates = SUM(highestState(:,:)) - SUM(lowestState(:,:)) + fi%input%jspins*fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
218 219 220 221 222 223 224 225
   ALLOCATE(occupationVec(numRelevantStates))
   occupationVec(:) = 0.0
   sumOcc = 0.0
   addCharge = 0.0
   subCharge = 0.0
   addChargeWeight = 0.0
   subChargeWeight = 0.0
   iState = 0
226 227 228
   DO jspin = 1, fi%input%jspins
      jsp = MERGE(1,jspin,fi%noco%l_noco)
      DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
229 230
         DO iBand = lowestState(ikpt,jspin), highestState(ikpt,jspin)
            iState = iState + 1
231
            occupationVec(iState) = results%w_iks(iBand,ikpt,jsp) / (fi%kpts%wtkpt(ikpt))
Gregor Michalicek's avatar
Gregor Michalicek committed
232
            sumOcc = sumOcc + results%w_iks(iBand,ikpt,jsp)
233 234
            IF(occupationVec(iState).LT.minOccB) THEN
               addCharge = addCharge + (minOccB-occupationVec(iState))*fi%kpts%wtkpt(ikpt)
235
               addChargeWeight = addChargeWeight + fi%kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
236
            END IF
237 238
            IF(occupationVec(iState).GT.(1.0-minOccB)) THEN
               subCharge = subCharge + (occupationVec(iState)-(1.0-minOccB))*fi%kpts%wtkpt(ikpt)
239
               subChargeWeight = subChargeWeight + fi%kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
240 241 242 243 244
            END IF
         END DO
      END DO
   END DO
   iState = 0
245 246 247
   DO jspin = 1, fi%input%jspins
      jsp = MERGE(1,jspin,fi%noco%l_noco)
      DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
248 249
         DO iBand = lowestState(ikpt,jspin), highestState(ikpt,jspin)
            iState = iState + 1
250
            IF(occupationVec(iState).LT.minOccB) THEN
251
               occupationVec(iState) = occupationVec(iState) + 0.5*(subCharge+addCharge)*(fi%kpts%wtkpt(ikpt)/addChargeWeight)
Gregor Michalicek's avatar
Gregor Michalicek committed
252
            END IF
253
            IF(occupationVec(iState).GT.(1.0-minOccB)) THEN
254
               occupationVec(iState) = occupationVec(iState) - 0.5*(subCharge+addCharge)*(fi%kpts%wtkpt(ikpt)/subChargeWeight)
Gregor Michalicek's avatar
Gregor Michalicek committed
255
            END IF
256
!            results%w_iks(iBand,ikpt,jsp) = occupationVec(iState) * fi%kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
257 258 259 260 261 262 263 264 265
         END DO
      END DO
   END DO
   DEALLOCATE(occupationVec)

   ! Some more initializations

   results%neig(:,:) = highestState(:,:) + 1

266 267 268
   ALLOCATE(overallVCoulSSDen(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
   ALLOCATE(vTotSSDen(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
   ALLOCATE(dEdOcc(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
269

270 271
   ALLOCATE(zintn_rSSDen(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
   ALLOCATE(vmdSSDen(MAXVAL(results%neig(1:fi%kpts%nkpt,1:fi%input%jspins)),fi%kpts%nkpt,fi%input%jspins))
Matthias Redies's avatar
Matthias Redies committed
272

Gregor Michalicek's avatar
Gregor Michalicek committed
273 274 275
   zintn_rSSDen(:,:,:) = 0.0
   vmdSSDen(:,:,:) = 0.0

276
   CALL regCharges%init(fi%input,fi%atoms)
277 278
   CALL dos%init(fi%input,fi%atoms,fi%kpts,fi%banddos,results%eig)
   CALL vacdos%init(fi%input,fi%atoms,fi%kpts,fi%banddos,results%eig)
279
   CALL moments%init(fmpi,fi%input,sphhar,fi%atoms)
280 281
   CALL overallDen%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_DEN)
   CALL overallVCoul%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_POTCOUL)
282
   IF (ALLOCATED(vTot%pw_w)) DEALLOCATE (vTot%pw_w)
283 284
   ALLOCATE(vTot%pw_w(SIZE(overallDen%pw,1),fi%input%jspins))
   DO jspin = 1, fi%input%jspins
285 286 287
      CALL convol(stars,vTot%pw_w(:,jspin),vTot%pw(:,jspin),stars%ufft)
   END DO

288
   CALL vTotTemp%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_POTTOT)
Gregor Michalicek's avatar
Gregor Michalicek committed
289 290 291
   vTotTemp = vTot

   DO jsp = 1,SIZE(vTot%mt,4)
292 293
      DO iAtom = 1,fi%atoms%ntype
         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)
Gregor Michalicek's avatar
Gregor Michalicek committed
294 295 296 297
      END DO
   END DO

   vCoul%pw_w = CMPLX(0.0,0.0)
298
   DO jspin = 1, fi%input%jspins
Gregor Michalicek's avatar
Gregor Michalicek committed
299 300 301
      CALL convol(stars,vCoul%pw_w(:,jspin),vCoul%pw(:,jspin),stars%ufft)
   END DO

302 303
   vTotSSDen = 0.0

304
   ! Calculate all single state densities
305

Gregor Michalicek's avatar
Gregor Michalicek committed
306
   numStates = 0
307 308
   DO jspin = 1, fi%input%jspins
      jsp = MERGE(1,jspin,fi%noco%l_noco)
Gregor Michalicek's avatar
Gregor Michalicek committed
309

310
      CALL cdnvalJob%init(fmpi,fi%input,fi%kpts,fi%noco,results,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
311

312 313
      DO ikpt_i = 1, SIZE(fmpi%k_list)
         ikpt= fmpi%k_list(ikpt_i)
Gregor Michalicek's avatar
Gregor Michalicek committed
314
         DO iBand = 1, highestState(ikpt,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
315
            numStates = numStates + 1
316
            ! Construct cdnvalJob object for this state
317
            ! (Reasonable parallelization is not yet done - should be placed over the loops enclosing this section)
Gregor Michalicek's avatar
Gregor Michalicek committed
318

319
            cdnvalJob%k_list=[ikpt]
Gregor Michalicek's avatar
Gregor Michalicek committed
320 321
!            cdnvalJob%ev_list=[iBand]
            cdnvalJob%weights(:,:) = 0.0
322
            cdnvalJob%weights(iBand,ikpt) = spinDegenFac
323

324
            ! Call cdnval to construct density
325 326
            WRITE(*,*) 'Note: some optional flags may have to be reset in rdmft before the cdnval call'
            WRITE(*,*) 'This is not yet implemented!'
327
            CALL singleStateDen%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_DEN)
328
            CALL cdnval(eig_id,fmpi,fi%kpts,jsp,fi%noco,nococonv,fi%input,fi%banddos,fi%cell,fi%atoms,enpara,stars,fi%vacuum,&
329
                        sphhar,fi%sym,vTot,fi%oned,cdnvalJob,singleStateDen,regCharges,dos,vacdos,results,moments,&
330
                        fi%gfinp,fi%hub1inp)
331 332

            ! Store the density on disc (These are probably way too many densities to keep them in memory)
333 334
            filename = ''
            WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
335
            IF (fmpi%irank.EQ.0) THEN
336
               CALL writeDensity(stars,fi%noco,fi%vacuum,fi%atoms,fi%cell,sphhar,fi%input,fi%sym,fi%oned,CDN_ARCHIVE_TYPE_CDN_const,CDN_input_DEN_const,&
337
                                 0,-1.0,0.0,-1.0,-1.0,.FALSE.,singleStateDen,TRIM(ADJUSTL(filename)))
338
            END IF
339
#ifdef CPP_MPI
340
            CALL mpi_bc_potden(fmpi,stars,sphhar,fi%atoms,fi%input,fi%vacuum,fi%oned,fi%noco,singleStateDen)
341 342
#endif
            ! For each state calculate Integral over KS effective potential times single state density
343
            potDenInt = 0.0
344
            CALL int_nv(jsp,stars,fi%vacuum,fi%atoms,sphhar,fi%cell,fi%sym,fi%input,fi%oned,vTotTemp,singleStateDen,potDenInt)
345
            vTotSSDen(iBand,ikpt,jsp) = potDenInt
Gregor Michalicek's avatar
Gregor Michalicek committed
346 347

            mt(:,:) = 0.0
348
            DO iType = 1, fi%atoms%ntype
349 350
               DO iGrid = 1, fi%atoms%jri(iType)
                  mt(iGrid,iType) = singleStateDen%mt(iGrid,0,iType,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
351
               END DO
Matthias Redies's avatar
Matthias Redies committed
352

353 354
               DO iGrid = 1, fi%atoms%jri(iType)
                  dpj(iGrid) = mt(iGrid,iType)/fi%atoms%rmsh(iGrid,iType)
Gregor Michalicek's avatar
Gregor Michalicek committed
355
               END DO
356

357
               CALL intgr3(dpj,fi%atoms%rmsh(1,iType),fi%atoms%dx(iType),fi%atoms%jri(iType),rhs)
Gregor Michalicek's avatar
Gregor Michalicek committed
358

359
               zintn_r(iType) = fi%atoms%neq(iType)*fi%atoms%zatom(iType)*sfp_const*rhs/2.0
Gregor Michalicek's avatar
Gregor Michalicek committed
360 361
               zintn_rSSDen(iBand,ikpt,jsp) = zintn_rSSDen(iBand,ikpt,jsp) + zintn_r(iType)

362
               CALL intgr3(mt(1,iType),fi%atoms%rmsh(1,iType),fi%atoms%dx(iType),fi%atoms%jri(iType),totz)
Gregor Michalicek's avatar
Gregor Michalicek committed
363

364
!               vmd(iType) = fi%atoms%rmt(iType)*vCoul%mt(fi%atoms%jri(iType),0,iType,1)/sfp_const + fi%atoms%zatom(iType) - totz*sfp_const
Gregor Michalicek's avatar
Gregor Michalicek committed
365
               vmd(iType) = -totz*sfp_const
366
               vmd(iType) = -fi%atoms%neq(iType)*fi%atoms%zatom(iType)*vmd(iType)/ (2.0*fi%atoms%rmt(iType))
Gregor Michalicek's avatar
Gregor Michalicek committed
367 368 369 370

               vmdSSDen(iBand,ikpt,jsp) = vmdSSDen(iBand,ikpt,jsp) + vmd(iType)
            END DO

371 372 373 374
         END DO
      END DO
   END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
375 376
   ! Initializations for exchange contributions

Gregor Michalicek's avatar
Gregor Michalicek committed
377 378
   WRITE(*,*) 'RDMFT: HF initializations start'

379 380 381
   IF(ALLOCATED(hybdat%ne_eig)) DEALLOCATE(hybdat%ne_eig)
   IF(ALLOCATED(hybdat%nbands)) DEALLOCATE(hybdat%nbands)
   IF(ALLOCATED(hybdat%nobd)) DEALLOCATE(hybdat%nobd)
382
   IF(ALLOCATED(hybdat%nbasm)) DEALLOCATE(hybdat%nbasm)
383
   IF(ALLOCATED(hybdat%div_vv)) DEALLOCATE(hybdat%div_vv)
384 385 386
   ALLOCATE(hybdat%ne_eig(fi%kpts%nkpt),hybdat%nbands(fi%kpts%nkpt),hybdat%nobd(fi%kpts%nkptf,fi%input%jspins))
   ALLOCATE(hybdat%nbasm(fi%kpts%nkptf))
   ALLOCATE(hybdat%div_vv(fi%input%neig,fi%kpts%nkpt,fi%input%jspins))
387

388
   l_zref = (fi%sym%zrfs.AND.(SUM(ABS(fi%kpts%bk(3,:fi%kpts%nkpt))).LT.1e-9).AND..NOT.fi%noco%l_noco)
Gregor Michalicek's avatar
Gregor Michalicek committed
389
   iterHF = 0
390
   hybdat%l_calhf = .TRUE.
Gregor Michalicek's avatar
Gregor Michalicek committed
391

392
   CALL mixedbasis(fi%atoms,fi%kpts,fi%input,fi%cell,xcpot,fi%mpinp,mpdata,fi%hybinp, hybdat,enpara,fmpi,vTot, iterHF)
Gregor Michalicek's avatar
Gregor Michalicek committed
393

Matthias Redies's avatar
Matthias Redies committed
394
   !allocate coulomb matrix
395 396 397 398
   IF (.NOT.ALLOCATED(hybdat%coul)) ALLOCATE(hybdat%coul(fi%kpts%nkpt))
   DO ikpt = 1, fi%kpts%nkpt
      CALL hybdat%coul(ikpt)%alloc(fi, mpdata%num_radbasfn, mpdata%n_g, ikpt)
   END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
399

400
   CALL glob_mpi%copy_mpi(fmpi)
Matthias Redies's avatar
Matthias Redies committed
401
   call work_pack%init(fi, hybdat, glob_mpi, jsp, glob_mpi%rank, glob_mpi%size)
Matthias Redies's avatar
Matthias Redies committed
402

Matthias Redies's avatar
Matthias Redies committed
403
   CALL coulombmatrix(fmpi, fi, mpdata, hybdat, xcpot)
404 405


406
   DO ikpt = 1, fi%kpts%nkpt
407
      CALL hybdat%coul(ikpt)%mpi_bc(fi, fmpi%mpi_comm, 0)
408
   END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
409

410
   CALL hf_init(eig_id,mpdata,fi,hybdat)
Gregor Michalicek's avatar
Gregor Michalicek committed
411

Gregor Michalicek's avatar
Gregor Michalicek committed
412
   WRITE(*,*) 'RDMFT: HF initializations end'
Gregor Michalicek's avatar
Gregor Michalicek committed
413

Gregor Michalicek's avatar
Gregor Michalicek committed
414 415
   maxHistoryLength = 7*numStates
   maxHistoryLength = 5*numStates
416
   maxHistoryLength = 23
Gregor Michalicek's avatar
Gregor Michalicek committed
417

418 419
   ALLOCATE(parent(fi%kpts%nkptf))
   ALLOCATE(exDiag(fi%input%neig,ikpt,fi%input%jspins))
Gregor Michalicek's avatar
Gregor Michalicek committed
420 421 422 423 424 425 426 427 428 429 430
   ALLOCATE(lastGradient(numStates+1))
   ALLOCATE(lastParameters(numStates+1))
   lastGradient = 0.0
   lastParameters = 0.0
   ALLOCATE(gradientCorrections(numStates+1,maxHistoryLength))
   ALLOCATE(paramCorrections(numStates+1,maxHistoryLength))
   gradientCorrections = 0.0
   paramCorrections = 0.0
   istep = 0

   ! Occupation number optimization loop
431

Gregor Michalicek's avatar
Gregor Michalicek committed
432 433
   convIter = 0

434
   converged = .FALSE.
435 436
   DO WHILE (.NOT.converged)

Gregor Michalicek's avatar
Gregor Michalicek committed
437
      WRITE(*,*) 'RDMFT: convergence loop start'
Gregor Michalicek's avatar
Gregor Michalicek committed
438 439
      convIter = convIter + 1
      WRITE(*,'(a,i7)') 'convIter: ', convIter
Gregor Michalicek's avatar
Gregor Michalicek committed
440

441 442
      DO jspin = 1, fi%input%jspins
         DO ikpt = 1,fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
443
            WRITE(*,*) 'jspin, ikpt: ', jspin, ikpt
444
            WRITE(*,'(10f12.7)') results%w_iks(1:10,ikpt,jspin)
Gregor Michalicek's avatar
Gregor Michalicek committed
445 446 447
         END DO
      END DO

448
      ! Calculate overall density with current occupation numbers (don't forget core electron density)
Gregor Michalicek's avatar
Gregor Michalicek committed
449
      CALL overallDen%resetPotDen()
450 451
      jspmax = fi%input%jspins
      IF (fi%noco%l_mperp) jspmax = 1
Gregor Michalicek's avatar
Gregor Michalicek committed
452

453
      DO jspin = 1,jspmax
454 455
         CALL cdnvalJob%init(fmpi,fi%input,fi%kpts,fi%noco,results,jspin)
         CALL cdnval(eig_id,fmpi,fi%kpts,jspin,fi%noco,nococonv,fi%input,fi%banddos,fi%cell,fi%atoms,enpara,stars,fi%vacuum,&
456
                     sphhar,fi%sym,vTot,fi%oned,cdnvalJob,overallDen,regCharges,dos,vacdos,results,moments,&
457
                     fi%gfinp,fi%hub1inp)
458
      END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
459

460
      CALL cdncore(fmpi,fi%oned,fi%input,fi%vacuum,fi%noco,nococonv,fi%sym,&
461
                   stars,fi%cell,sphhar,fi%atoms,vTot,overallDen,moments,results)
462 463
      IF (fmpi%irank.EQ.0) THEN
         CALL qfix(fmpi,stars,fi%atoms,fi%sym,fi%vacuum,sphhar,fi%input,fi%cell,fi%oned,overallDen,&
Uliana Alekseeva's avatar
Uliana Alekseeva committed
464
                   fi%noco%l_noco,.TRUE.,l_par=.FALSE.,force_fix=.TRUE.,fix=fix)
465 466
      END IF
#ifdef CPP_MPI
467
      CALL mpi_bc_potden(fmpi,stars,sphhar,fi%atoms,fi%input,fi%vacuum,fi%oned,fi%noco,overallDen)
468
#endif
469 470

      ! Calculate Coulomb potential for overall density (+including external potential)
471 472
      CALL overallDen%sum_both_spin()!workden)
      CALL overallVCoul%resetPotDen()
Gregor Michalicek's avatar
Gregor Michalicek committed
473 474
      ALLOCATE(overallVCoul%pw_w(size(overallVCoul%pw,1),size(overallVCoul%pw,2)))
      overallVCoul%pw_w(:,:) = 0.0
475
      CALL vgen_coulomb(1,fmpi,fi%oned,fi%input,fi%field,fi%vacuum,fi%sym,stars,fi%cell,sphhar,fi%atoms,.FALSE.,overallDen,overallVCoul)
Gregor Michalicek's avatar
Gregor Michalicek committed
476
      CALL convol(stars,overallVCoul%pw_w(:,1),overallVCoul%pw(:,1),stars%ufft)   ! Is there a problem with a second spin?!
477
#ifdef CPP_MPI
478
      CALL mpi_bc_potden(fmpi,stars,sphhar,fi%atoms,fi%input,fi%vacuum,fi%oned,fi%noco,overallVCoul)
479
#endif
480

481
      overallVCoulSSDen = 0.0
482 483 484
      DO jspin = 1, fi%input%jspins
         jsp = MERGE(1,jspin,fi%noco%l_noco)
         DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
485
            DO iBand = 1, highestState(ikpt,jsp)
486 487 488
               ! Read the single-state density from disc
               filename = ''
               WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
489
               IF (fmpi%irank.EQ.0) THEN
490 491
                  CALL readDensity(stars,fi%noco,fi%vacuum,fi%atoms,fi%cell,sphhar,fi%input,fi%sym,fi%oned,CDN_ARCHIVE_TYPE_CDN_const,&
                                   CDN_input_DEN_const,0,fermiEnergyTemp,l_qfix,singleStateDen,TRIM(ADJUSTL(filename)))
492 493 494
                  CALL singleStateDen%sum_both_spin()!workden)
               END IF
#ifdef CPP_MPI
495
               CALL mpi_bc_potden(fmpi,stars,sphhar,fi%atoms,fi%input,fi%vacuum,fi%oned,fi%noco,singleStateDen)
496
#endif
497

498
               ! For each state calculate integral over Coulomb potential times single state density
499
               potDenInt = 0.0
500
               CALL int_nv(1,stars,fi%vacuum,fi%atoms,sphhar,fi%cell,fi%sym,fi%input,fi%oned,vCoul,singleStateDen,potDenInt) ! Is there a problem with a second spin?!
501 502 503 504
               overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
            END DO
         END DO
      END DO
505

Gregor Michalicek's avatar
Gregor Michalicek committed
506 507 508
      ! Construct exchange matrix diagonal

      exDiag = 0.0
509 510
      DO jspin = 1, fi%input%jspins
         jsp = MERGE(1,jspin,fi%noco%l_noco)
Gregor Michalicek's avatar
Gregor Michalicek committed
511 512
         ! remove weights(wtkpt) in w_iks
         wl_iks = 0.0
513
         DO ikptf = 1, fi%kpts%nkptf
514
            ikpt = fi%kpts%bkp(ikptf)
Gregor Michalicek's avatar
Gregor Michalicek committed
515
            DO iBand=1, results%neig(ikpt,jsp)
516
               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
Gregor Michalicek's avatar
Gregor Michalicek committed
517
               wl_iks(iBand,ikptf) = sqrt(wl_iks(iBand,ikptf)) ! For Müller functional
518
               wl_iks(iBand,ikptf) = wl_iks(iBand,ikptf) / fi%kpts%nkptf ! This is the last part of two lines above
Gregor Michalicek's avatar
Gregor Michalicek committed
519 520 521 522 523 524 525 526 527 528
            END DO
         END DO

         IF(ALLOCATED(eig_irr)) DEALLOCATE (eig_irr)
         IF(ALLOCATED(hybdat%kveclo_eig)) DEALLOCATE (hybdat%kveclo_eig)
         IF(ALLOCATED(hybdat%pntgptd)) DEALLOCATE (hybdat%pntgptd)
         IF(ALLOCATED(hybdat%pntgpt)) DEALLOCATE (hybdat%pntgpt)
         IF(ALLOCATED(hybdat%prodm)) DEALLOCATE (hybdat%prodm)
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE (hybdat%nindxp1)

Gregor Michalicek's avatar
Gregor Michalicek committed
529 530
         results%neig(:,:) = neigTemp(:,:)

531
         CALL HF_setup(mpdata,fi,fmpi,nococonv,results,jspin,enpara,&
Matthias Redies's avatar
Matthias Redies committed
532
                       hybdat,vTot%mt(:,0,:,:),eig_irr)
Gregor Michalicek's avatar
Gregor Michalicek committed
533

Gregor Michalicek's avatar
Gregor Michalicek committed
534 535
         results%neig(:,:) = highestState(:,:) + 1

536
         DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
537

538
            CALL lapw%init(fi%input,fi%noco,nococonv,fi%kpts,fi%atoms,fi%sym,ikpt,fi%cell,l_zref)
Gregor Michalicek's avatar
Gregor Michalicek committed
539

540 541 542 543 544 545 546
            nbasfcn = 0
            IF(fi%noco%l_noco) then
               nbasfcn = lapw%nv(1) + lapw%nv(2) + 2*fi%atoms%nlotot
            ELSE
               nbasfcn = lapw%nv(1) + fi%atoms%nlotot
            END IF

547
            parent = 0
548
            CALL zMat%init(fi%sym%invs,nbasfcn,fi%input%neig)
Matthias Redies's avatar
Matthias Redies committed
549

550 551
            if(ikpt /= fi%kpts%bkp(ikpt)) call juDFT_error("We should be reading the parent z-mat here!")
            call read_z(fi%atoms, fi%cell, hybdat, fi%kpts, fi%sym, fi%noco, nococonv,  fi%input, ikpt, jsp, zMat, c_phase=c_phase)
Matthias Redies's avatar
Matthias Redies committed
552 553 554 555
            allocate(cmt_nk(hybdat%nbands(ikpt), hybdat%maxlmindx, fi%atoms%nat), stat=ierr)
            if(ierr  /= 0) call judft_error("can't allocate cmt_nk")
            call calc_cmt(fi%atoms, fi%cell, fi%input, fi%noco, nococonv, fi%hybinp, hybdat, mpdata, fi%kpts, &
                        fi%sym, fi%oneD, zMat, jsp, ikpt, c_phase, cmt_nk)
556

557 558 559
            ALLOCATE (indx_sest(hybdat%nbands(ikpt), hybdat%nbands(ikpt)))
            indx_sest = 0

Matthias Redies's avatar
Matthias Redies committed
560
            call symm_hf_init(fi,ikpt,nsymop,rrot,psym)
Matthias Redies's avatar
Matthias Redies committed
561
            call symm_hf(fi,ikpt,hybdat,work_pack%k_packs(ikpt)%submpi, eig_irr,mpdata,lapw,nococonv, zMat, c_phase,jspin,&
Matthias Redies's avatar
Matthias Redies committed
562
                         rrot,nsymop,psym,n_q,parent,nsest,indx_sest)
Gregor Michalicek's avatar
Gregor Michalicek committed
563

564
            exMat%l_real=fi%sym%invs
Matthias Redies's avatar
Matthias Redies committed
565
            CALL exchange_valence_hf(work_pack%k_packs(ikpt),fi,zMat, mpdata,jspin,hybdat,lapw,&
Matthias Redies's avatar
Matthias Redies committed
566
                                     eig_irr,results,n_q,wl_iks,xcpot,nococonv,stars,nsest,indx_sest,&
Matthias Redies's avatar
Matthias Redies committed
567 568
                                     fmpi, cmt_nk, exMat)
            deallocate(cmt_nk)
Matthias Redies's avatar
Matthias Redies committed
569
            CALL exchange_vccv1(ikpt,fi, nococonv,mpdata,hybdat,jspin,lapw,glob_mpi,nsymop,nsest,indx_sest,&
Matthias Redies's avatar
Matthias Redies committed
570
                                1.0,results,cmt_nk,exMat)
Gregor Michalicek's avatar
Gregor Michalicek committed
571

572 573
            DEALLOCATE(indx_sest)

574
            !Start of workaround for increased functionality of fi%symmetrizeh (call it))
Gregor Michalicek's avatar
Gregor Michalicek committed
575

576
            nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*fi%atoms%nlotot,lapw%nv(1)+fi%atoms%nlotot,fi%noco%l_noco)
Gregor Michalicek's avatar
Gregor Michalicek committed
577

Matthias Redies's avatar
Matthias Redies committed
578
            CALL read_eig(hybdat%eig_id,ikpt,jspin, smat=olap)
Gregor Michalicek's avatar
Gregor Michalicek committed
579

580
            zMat%matsize2 = hybdat%nbands(ikpt) ! reduce "visible matsize" for the following computations
Gregor Michalicek's avatar
Gregor Michalicek committed
581 582 583

            CALL olap%multiply(zMat,trafo)

584
            CALL invtrafo%alloc(olap%l_real,hybdat%nbands(ikpt),nbasfcn)
Gregor Michalicek's avatar
Gregor Michalicek committed
585
            CALL trafo%TRANSPOSE(invtrafo)
Matthias Redies's avatar
Matthias Redies committed
586

587 588
            DO iBand = 1, hybdat%nbands(ikpt)
               DO jBand = 1, iBand-1
Gregor Michalicek's avatar
Gregor Michalicek committed
589
                  IF (exMat%l_real) THEN
590
                     exMat%data_r(iBand,jBand)=exMat%data_r(jBand,iBand)
Gregor Michalicek's avatar
Gregor Michalicek committed
591
                  ELSE
592
                     exMat%data_c(iBand,jBand)=conjg(exMat%data_c(jBand,iBand))
Gregor Michalicek's avatar
Gregor Michalicek committed
593 594 595 596 597 598 599
                  END IF
               END DO
            END DO

            CALL exMat%multiply(invtrafo,tmpMat)
            CALL trafo%multiply(tmpMat,exMatLAPW)

600
            call symmetrizeh(fi%atoms,fi%kpts%bkf(:,ikpt),jspin,lapw,fi%sym,hybdat%kveclo_eig,fi%cell,nsymop,psym,exMatLAPW)
Gregor Michalicek's avatar
Gregor Michalicek committed
601

602
            IF (.NOT.exMatLAPW%l_real) exMatLAPW%data_c=conjg(exMatLAPW%data_c)
Gregor Michalicek's avatar
Gregor Michalicek committed
603 604 605 606
            zMat%matsize1=MIN(zMat%matsize1,exMatLAPW%matsize2)

            CALL exMatLAPW%multiply(zMat,tmpMat)

Gregor Michalicek's avatar
Gregor Michalicek committed
607
            DO iBand = 1, highestState(ikpt,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
608 609
               IF (zMat%l_real) THEN
                  exDiag(iBand,ikpt,jspin) = dot_product(zMat%data_r(:zMat%matsize1,iband),tmpMat%data_r(:,iband))
Gregor Michalicek's avatar
Gregor Michalicek committed
610
               ELSE
Gregor Michalicek's avatar
Gregor Michalicek committed
611
                  exDiag(iBand,ikpt,jspin) = REAL(dot_product(zMat%data_c(:zMat%matsize1,iband),tmpMat%data_c(:,iband)))
Gregor Michalicek's avatar
Gregor Michalicek committed
612 613
               END IF
            END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
614

615
            !End of workaround for increased functionality of fi%symmetrizeh (call it))
Gregor Michalicek's avatar
Gregor Michalicek committed
616 617 618 619 620 621 622 623

!            DO iBand = 1, highestState(ikpt,jsp)
!               IF (exMat%l_real) THEN
!                  exDiag(iBand,ikpt,jspin) = exMat%data_r(iBand,iBand)
!               ELSE
!                  exDiag(iBand,ikpt,jspin) = REAL(exMat%data_c(iBand,iBand))
!               END IF
!            END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
624 625 626 627 628
         END DO
      END DO

      ! Calculate total energy derivative with respect to occupations (dEdOcc)

Gregor Michalicek's avatar
Gregor Michalicek committed
629
      occSum = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
630
      gradSum = 0.0
631 632
      DO ispin = 1, fi%input%jspins
         isp = MERGE(1,ispin,fi%noco%l_noco)
633
!         CALL cdnvalJob%init(fmpi,fi%input,fi%kpts,fi%noco,results,isp,fi%banddos=fi%banddos)
634
         DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
635
            DO iBand = 1, highestState(ikpt,isp)
636
               occStateI = results%w_iks(iBand,ikpt,isp) / (fi%kpts%wtkpt(ikpt))!*fi%kpts%nkptf)
637
               occStateI = MAX(occStateI,minOcc)
Gregor Michalicek's avatar
Gregor Michalicek committed
638 639 640

               occStateI = MIN(occStateI,1.0-minOcc)

Gregor Michalicek's avatar
Gregor Michalicek committed
641
!               IF(occStateI.LT.1.0e-7) occStateI = 5.0e-4 ! This is preliminary. I have to discuss what do do here.
Gregor Michalicek's avatar
Gregor Michalicek committed
642
!               occStateI = cdnvalJob%weights(iBand,ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
643 644
               rdmftFunctionalValue = 0.5*0.5*SQRT(1.0/occStateI) ! for Müller functional derivative

Gregor Michalicek's avatar
Gregor Michalicek committed
645 646 647 648 649 650 651 652 653 654 655 656 657 658 659

               !!! Test start
!                  occStateI = MIN(occStateI,1.0-minOcc)
!                  rdmftFunctionalValue = 0.5 * 0.5*pi_const*COS(0.5*pi_const*occStateI)

!                  rdmftFunctionalValue = ASIN(SQRT(occStateI)) * 2.0 / pi_const
!                  rdmftFunctionalValue = (SIN(rdmftFunctionalValue*pi_const/2.0)) * (COS(rdmftFunctionalValue*pi_const/2.0))**2.0 * pi_const
               !!! Test 2:
!                  occStateI = MIN(occStateI,1.0-minOcc)
!                  rdmftFunctionalValue = ASIN(SQRT(occStateI)) * 2.0 / pi_const
!                  rdmftFunctionalValue = COS(rdmftFunctionalValue*pi_const/2.0) * pi_const / 4.0
               !!! Test end

               occSum = occSum + results%w_iks(iBand,ikpt,isp)

660
               exchangeTerm = - rdmftFunctionalValue * exDiag(iBand,ikpt,isp) * fi%kpts%wtkpt(ikpt) * spinDegenFac !*fi%kpts%nkptf
Gregor Michalicek's avatar
Gregor Michalicek committed
661 662

               dEdOcc(iBand,ikpt,isp) = +((spinDegenFac * results%eig(iBand,ikpt,isp)) - vTotSSDen(iBand,ikpt,isp) + &
Gregor Michalicek's avatar
Gregor Michalicek committed
663 664 665 666 667
                                          overallVCoulSSDen(iBand,ikpt,isp) - &
                                          zintn_rSSDen(iBand,ikpt,isp) + vmdSSDen(iBand,ikpt,isp) + exchangeTerm )

               theta = ASIN(SQRT(occStateI))! * 2.0 /  pi_const
               dEdOcc(iBand,ikpt,isp) = 2.0 * sin(theta) * cos(theta) * dEdOcc(iBand,ikpt,isp)
668
!               dEdOcc(iBand,ikpt,isp) = dEdOcc(iBand,ikpt,isp) + 2.0 * COS(theta) * exDiag(iBand,ikpt,isp) * fi%kpts%wtkpt(ikpt) * spinDegenFac
Gregor Michalicek's avatar
Gregor Michalicek committed
669

Gregor Michalicek's avatar
Gregor Michalicek committed
670
               WRITE(*,*) 'ENERGY GRADIENT CONTRIBUTIONS'
Gregor Michalicek's avatar
Gregor Michalicek committed
671
               WRITE(*,*) 'ispin, ikpt, iBand, weight', ispin, ikpt, iBand, results%w_iks(iBand,ikpt,isp)
Gregor Michalicek's avatar
Gregor Michalicek committed
672 673 674
               WRITE(*,*) 'results%eig(iBand,ikpt,isp)', results%eig(iBand,ikpt,isp)
               WRITE(*,*) 'vTotSSDen(iBand,ikpt,isp)', vTotSSDen(iBand,ikpt,isp)
               WRITE(*,*) 'overallVCoulSSDen(iBand,ikpt,isp)', overallVCoulSSDen(iBand,ikpt,isp)
Gregor Michalicek's avatar
Gregor Michalicek committed
675 676
               WRITE(*,*) 'zintn_rSSDen(iBand,ikpt,isp)', zintn_rSSDen(iBand,ikpt,isp)
               WRITE(*,*) 'vmdSSDen(iBand,ikpt,isp)', vmdSSDen(iBand,ikpt,isp)
Gregor Michalicek's avatar
Gregor Michalicek committed
677 678 679
               WRITE(*,*) 'exchangeTerm', exchangeTerm
               WRITE(*,*) 'exDiag(iBand,ikpt,isp)', exDiag(iBand,ikpt,isp)
               WRITE(*,*) 'rdmftFunctionalValue', rdmftFunctionalValue
Gregor Michalicek's avatar
Gregor Michalicek committed
680

Gregor Michalicek's avatar
Gregor Michalicek committed
681
               gradSum = gradSum + dEdOcc(iBand,ikpt,isp) !* results%w_iks(iBand,ikpt,isp))
Gregor Michalicek's avatar
Gregor Michalicek committed
682 683 684
            END DO
         END DO
      END DO
685
      lagrangeMultiplier = -gradSum / numStates ! occSum  !(fi%input%zelec/(2.0/REAL(fi%input%jspins)))
Gregor Michalicek's avatar
Gregor Michalicek committed
686

Gregor Michalicek's avatar
Gregor Michalicek committed
687
   WRITE(*,*) 'lagrangeMultiplier: ', lagrangeMultiplier
Gregor Michalicek's avatar
Gregor Michalicek committed
688

689 690
      ! Optimize occupation numbers

Gregor Michalicek's avatar
Gregor Michalicek committed
691 692
      ALLOCATE (gradient(numStates+1))
      ALLOCATE (parameters(numStates+1))
693
      ALLOCATE (paramGroup(numStates+1))
Gregor Michalicek's avatar
Gregor Michalicek committed
694 695 696 697 698 699 700 701 702 703 704 705
      ALLOCATE (equalityLinCombi(numStates+1))
      ALLOCATE (minConstraints(numStates+1))
      ALLOCATE (maxConstraints(numStates+1))
      ALLOCATE (enabledConstraints(numStates+1))
      enabledConstraints = .FALSE.
      enabledConstraints(1:numStates) = .TRUE.
      minConstraints(:) = 0.0
      maxConstraints(:) = 1.0
      equalityLinCombi = 0.0
      gradient = 0.0
      parameters = 0.0
      iState = 0
706 707 708
      DO ispin = 1, fi%input%jspins
         isp = MERGE(1,ispin,fi%noco%l_noco)
         DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
709 710
            DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
               iState = iState + 1
711
               occStateI = results%w_iks(iBand,ikpt,isp) / fi%kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
712

713
               occStateI = MAX(occStateI,minOcc)
Gregor Michalicek's avatar
Gregor Michalicek committed
714 715 716 717
               occStateI = MIN(occStateI,1.0-minOcc)

               theta = ASIN(SQRT(occStateI))! * 2.0 /  pi_const

718
               WRITE(7865,'(i7,4f15.10)') iState, occStateI, theta, sin(theta), cos(theta)
Gregor Michalicek's avatar
Gregor Michalicek committed
719 720

!               occStateI = MAX(occStateI,minOcc)
721
               equalityLinCombi(iState) = fi%kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
722 723 724 725 726

!               dEdOcc(iBand,ikpt,isp) = dEdOcc(iBand,ikpt,isp) + lagrangeMultiplier

!               dEdOcc(iBand,ikpt,isp) = 2.0 * sin(theta) * cos(theta) * dEdOcc(iBand,ikpt,isp)

Gregor Michalicek's avatar
Gregor Michalicek committed