rdmft.F90 42.2 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 12 13
SUBROUTINE rdmft(eig_id,mpi,fi,enpara,stars,&
                 sphhar,vTot,vCoul,nococonv,xcpot,mpdata,hybdat,&
                 results,archiveType,outDen)
14 15 16

   USE m_types
   USE m_juDFT
17
   USE m_constants
Gregor Michalicek's avatar
Gregor Michalicek committed
18
   USE m_intgr, ONLY : intgr3
Gregor Michalicek's avatar
Gregor Michalicek committed
19
   USE m_eig66_io
20
#ifndef CPP_OLDINTEL
21
   USE m_cdnval
Gregor Michalicek's avatar
Gregor Michalicek committed
22
   USE m_cdngen
23
   USE m_cdn_io
24 25
   USE m_cdncore
   USE m_qfix
26 27 28
   USE m_vgen_coulomb
   USE m_convol
   USE m_intnv
Gregor Michalicek's avatar
Gregor Michalicek committed
29 30 31 32 33

   USE m_mixedbasis
   USE m_coulombmatrix
   USE m_hf_init
   USE m_hf_setup
Matthias Redies's avatar
Matthias Redies committed
34
   USE m_io_hybinp
Gregor Michalicek's avatar
Gregor Michalicek committed
35 36 37
   USE m_symm_hf
   USE m_exchange_valence_hf
   USE m_exchange_core
Gregor Michalicek's avatar
Gregor Michalicek committed
38
   USE m_symmetrizeh
Gregor Michalicek's avatar
Gregor Michalicek committed
39
   USE m_bfgs_b2
40
   USE m_xmlOutput
Gregor Michalicek's avatar
Gregor Michalicek committed
41

42 43
#ifdef CPP_MPI
   USE m_mpi_bc_potden
44
#endif
45
#endif
46 47 48 49

   IMPLICIT NONE

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

   INTEGER,               INTENT(IN)    :: eig_id
Gregor Michalicek's avatar
Gregor Michalicek committed
64
   INTEGER,               INTENT(IN)    :: archiveType
Matthias Redies's avatar
merge  
Matthias Redies committed
65
   TYPE(t_potden)                       :: EnergyDen
66

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

95 96 97 98 99 100
   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
101

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

104
   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
105

106 107
   REAL, ALLOCATABLE                    :: overallVCoulSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vTotSSDen(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
108 109
   REAL, ALLOCATABLE                    :: dEdOcc(:,:,:)

Gregor Michalicek's avatar
Gregor Michalicek committed
110 111 112 113
   REAL, ALLOCATABLE                    :: zintn_rSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vmdSSDen(:,:,:)


Gregor Michalicek's avatar
Gregor Michalicek committed
114 115 116 117 118 119 120 121 122 123 124 125
   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(:)
126

127
   INTEGER, ALLOCATABLE                 :: paramGroup(:)
128
   INTEGER, ALLOCATABLE                 :: indx_sest(:,:)
129
   INTEGER, ALLOCATABLE                 :: parent(:)
Gregor Michalicek's avatar
Gregor Michalicek committed
130 131 132 133
   INTEGER, ALLOCATABLE                 :: pointer_EIBZ(:)
   INTEGER, ALLOCATABLE                 :: n_q(:)

   LOGICAL, ALLOCATABLE                 :: enabledConstraints(:)
134
   type(t_hybmpi)    :: hybmpi
Gregor Michalicek's avatar
Gregor Michalicek committed
135

136
   complex :: c_phase(fi%input%neig)
137

Gregor Michalicek's avatar
Gregor Michalicek committed
138
#endif
139

140
#ifndef CPP_OLDINTEL
Gregor Michalicek's avatar
Gregor Michalicek committed
141 142 143

   WRITE(*,*) 'entered RDMFT subroutine'

144
   ! General initializations
145
   mixParam = 0.001
Gregor Michalicek's avatar
Gregor Michalicek committed
146
   lagrangeMultiplier = 0.1 !results%ef
147
   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
148 149 150 151 152 153

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

   ! Determine which states have to be considered
   lowestState(:,:) = -1
   highestState(:,:) = -1
154 155 156
   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
157 158
         ! determine lowest state
         DO iBand = 1, results%neig(ikpt,jsp)
159
            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
160 161 162 163
               lowestState(ikpt,jspin) = iBand
               EXIT
            END IF
         END DO
164
         lowestState(ikpt,jspin) = lowestState(ikpt,jspin) - fi%input%rdmftStatesBelow
Gregor Michalicek's avatar
Gregor Michalicek committed
165 166 167 168 169 170 171 172 173 174
         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
175
            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
176 177 178 179
               highestState(ikpt,jspin) = iBand
               EXIT
            END IF
         END DO
180
         highestState(ikpt,jspin) = highestState(ikpt,jspin) + fi%input%rdmftStatesAbove
Gregor Michalicek's avatar
Gregor Michalicek committed
181
         IF((results%neig(ikpt,jsp)-1).LT.highestState(ikpt,jspin)) THEN
182 183 184 185
            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
186 187 188 189 190 191 192 193 194
            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
195 196 197 198 199
            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
200 201 202 203 204
            CALL juDFT_error('Highest state is degenerate', calledby = 'rdmft')
         END IF
      END DO
   END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
205 206 207
   IF (ANY(results%w_iksRDMFT(:,:,:).NE.0.0)) THEN
      results%w_iks(:,:,:) = results%w_iksRDMFT(:,:,:)
   END IF
208
!   results%w_iksRDMFT(:,:,:) = results%w_iks(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
209

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

   ! Some more initializations

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

260 261 262
   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))
263

264 265
   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
266

Gregor Michalicek's avatar
Gregor Michalicek committed
267 268 269
   zintn_rSSDen(:,:,:) = 0.0
   vmdSSDen(:,:,:) = 0.0

270 271 272 273 274
   CALL regCharges%init(fi%input,fi%atoms)
   CALL dos%init(fi%input,fi%atoms,fi%kpts,fi%vacuum)
   CALL moments%init(mpi,fi%input,sphhar,fi%atoms)
   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)
275
   IF (ALLOCATED(vTot%pw_w)) DEALLOCATE (vTot%pw_w)
276 277
   ALLOCATE(vTot%pw_w(SIZE(overallDen%pw,1),fi%input%jspins))
   DO jspin = 1, fi%input%jspins
278 279 280
      CALL convol(stars,vTot%pw_w(:,jspin),vTot%pw(:,jspin),stars%ufft)
   END DO

281
   CALL vTotTemp%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_POTTOT)
Gregor Michalicek's avatar
Gregor Michalicek committed
282 283 284
   vTotTemp = vTot

   DO jsp = 1,SIZE(vTot%mt,4)
285 286
      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
287 288 289 290
      END DO
   END DO

   vCoul%pw_w = CMPLX(0.0,0.0)
291
   DO jspin = 1, fi%input%jspins
Gregor Michalicek's avatar
Gregor Michalicek committed
292 293 294
      CALL convol(stars,vCoul%pw_w(:,jspin),vCoul%pw(:,jspin),stars%ufft)
   END DO

295 296
   vTotSSDen = 0.0

297
   ! Calculate all single state densities
298

Gregor Michalicek's avatar
Gregor Michalicek committed
299
   numStates = 0
300 301
   DO jspin = 1, fi%input%jspins
      jsp = MERGE(1,jspin,fi%noco%l_noco)
Gregor Michalicek's avatar
Gregor Michalicek committed
302

303
      CALL cdnvalJob%init(mpi,fi%input,fi%kpts,fi%noco,results,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
304

305 306
      DO ikpt_i = 1, SIZE(mpi%k_list)
         ikpt= mpi%k_list(ikpt_i)
Gregor Michalicek's avatar
Gregor Michalicek committed
307
         DO iBand = 1, highestState(ikpt,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
308
            numStates = numStates + 1
309
            ! Construct cdnvalJob object for this state
310
            ! (Reasonable parallelization is not yet done - should be placed over the loops enclosing this section)
Gregor Michalicek's avatar
Gregor Michalicek committed
311

312
            cdnvalJob%k_list=[ikpt]
Gregor Michalicek's avatar
Gregor Michalicek committed
313 314
!            cdnvalJob%ev_list=[iBand]
            cdnvalJob%weights(:,:) = 0.0
315
            cdnvalJob%weights(iBand,ikpt) = spinDegenFac
316

317
            ! Call cdnval to construct density
318 319
            WRITE(*,*) 'Note: some optional flags may have to be reset in rdmft before the cdnval call'
            WRITE(*,*) 'This is not yet implemented!'
320 321 322 323
            CALL singleStateDen%init(stars,fi%atoms,sphhar,fi%vacuum,fi%noco,fi%input%jspins,POTDEN_TYPE_DEN)
            CALL cdnval(eig_id,mpi,fi%kpts,jsp,fi%noco,nococonv,fi%input,fi%banddos,fi%cell,fi%atoms,enpara,stars,fi%vacuum,&
                        sphhar,fi%sym,vTot,fi%oned,cdnvalJob,singleStateDen,regCharges,dos,results,moments,&
                        fi%gfinp,fi%hub1inp)
324 325

            ! Store the density on disc (These are probably way too many densities to keep them in memory)
326 327
            filename = ''
            WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
328
            IF (mpi%irank.EQ.0) THEN
329
               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,&
330 331
                                 0,-1.0,0.0,.FALSE.,singleStateDen,TRIM(ADJUSTL(filename)))
            END IF
332
#ifdef CPP_MPI
333
            CALL mpi_bc_potden(mpi,stars,sphhar,fi%atoms,fi%input,fi%vacuum,fi%oned,fi%noco,singleStateDen)
334 335
#endif
            ! For each state calculate Integral over KS effective potential times single state density
336
            potDenInt = 0.0
337
            CALL int_nv(jsp,stars,fi%vacuum,fi%atoms,sphhar,fi%cell,fi%sym,fi%input,fi%oned,vTotTemp,singleStateDen,potDenInt)
338
            vTotSSDen(iBand,ikpt,jsp) = potDenInt
Gregor Michalicek's avatar
Gregor Michalicek committed
339 340

            mt(:,:) = 0.0
341
            DO iType = 1, fi%atoms%ntype
342 343
               DO iGrid = 1, fi%atoms%jri(iType)
                  mt(iGrid,iType) = singleStateDen%mt(iGrid,0,iType,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
344
               END DO
Matthias Redies's avatar
Matthias Redies committed
345

346 347
               DO iGrid = 1, fi%atoms%jri(iType)
                  dpj(iGrid) = mt(iGrid,iType)/fi%atoms%rmsh(iGrid,iType)
Gregor Michalicek's avatar
Gregor Michalicek committed
348
               END DO
349

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

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

355
               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
356

357
!               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
358
               vmd(iType) = -totz*sfp_const
359
               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
360 361 362 363

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

364 365 366 367
         END DO
      END DO
   END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
368 369
   ! Initializations for exchange contributions

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

372 373 374
   IF(ALLOCATED(hybdat%ne_eig)) DEALLOCATE(hybdat%ne_eig)
   IF(ALLOCATED(hybdat%nbands)) DEALLOCATE(hybdat%nbands)
   IF(ALLOCATED(hybdat%nobd)) DEALLOCATE(hybdat%nobd)
375
   IF(ALLOCATED(hybdat%nbasm)) DEALLOCATE(hybdat%nbasm)
376
   IF(ALLOCATED(hybdat%div_vv)) DEALLOCATE(hybdat%div_vv)
377 378 379
   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))
380

381
   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
382
   iterHF = 0
383
   hybdat%l_calhf = .TRUE.
Gregor Michalicek's avatar
Gregor Michalicek committed
384

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

Matthias Redies's avatar
Matthias Redies committed
387
   !allocate coulomb matrix
388 389 390 391
   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
392

393
   CALL coulombmatrix(mpi, fi, mpdata, hybdat, xcpot, [(ikpt,ikpt=1,fi%kpts%nkpt)])
Matthias Redies's avatar
Matthias Redies committed
394
   
395 396 397 398
   CALL hybmpi%copy_mpi(mpi)
   DO ikpt = 1, fi%kpts%nkpt
      CALL hybdat%coul(ikpt)%mpi_ibc(fi, hybmpi, 0)
   END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
399

400
   CALL hf_init(eig_id,mpdata,fi,hybdat)
Gregor Michalicek's avatar
Gregor Michalicek committed
401

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

Gregor Michalicek's avatar
Gregor Michalicek committed
404 405
   maxHistoryLength = 7*numStates
   maxHistoryLength = 5*numStates
406
   maxHistoryLength = 23
Gregor Michalicek's avatar
Gregor Michalicek committed
407

408 409
   ALLOCATE(parent(fi%kpts%nkptf))
   ALLOCATE(exDiag(fi%input%neig,ikpt,fi%input%jspins))
Gregor Michalicek's avatar
Gregor Michalicek committed
410 411 412 413 414 415 416 417 418 419 420
   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
421

Gregor Michalicek's avatar
Gregor Michalicek committed
422 423
   convIter = 0

424
   converged = .FALSE.
425 426
   DO WHILE (.NOT.converged)

Gregor Michalicek's avatar
Gregor Michalicek committed
427
      WRITE(*,*) 'RDMFT: convergence loop start'
Gregor Michalicek's avatar
Gregor Michalicek committed
428 429
      convIter = convIter + 1
      WRITE(*,'(a,i7)') 'convIter: ', convIter
Gregor Michalicek's avatar
Gregor Michalicek committed
430

431 432
      DO jspin = 1, fi%input%jspins
         DO ikpt = 1,fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
433
            WRITE(*,*) 'jspin, ikpt: ', jspin, ikpt
434
            WRITE(*,'(10f12.7)') results%w_iks(1:10,ikpt,jspin)
Gregor Michalicek's avatar
Gregor Michalicek committed
435 436 437
         END DO
      END DO

438
      ! Calculate overall density with current occupation numbers (don't forget core electron density)
Gregor Michalicek's avatar
Gregor Michalicek committed
439
      CALL overallDen%resetPotDen()
440 441
      jspmax = fi%input%jspins
      IF (fi%noco%l_mperp) jspmax = 1
Gregor Michalicek's avatar
Gregor Michalicek committed
442

443
      DO jspin = 1,jspmax
444 445 446 447
         CALL cdnvalJob%init(mpi,fi%input,fi%kpts,fi%noco,results,jspin)
         CALL cdnval(eig_id,mpi,fi%kpts,jspin,fi%noco,nococonv,fi%input,fi%banddos,fi%cell,fi%atoms,enpara,stars,fi%vacuum,&
                     sphhar,fi%sym,vTot,fi%oned,cdnvalJob,overallDen,regCharges,dos,results,moments,&
                     fi%gfinp,fi%hub1inp)
448
      END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
449

450 451
      CALL cdncore(mpi,fi%oned,fi%input,fi%vacuum,fi%noco,nococonv,fi%sym,&
                   stars,fi%cell,sphhar,fi%atoms,vTot,overallDen,moments,results)
452
      IF (mpi%irank.EQ.0) THEN
Uliana Alekseeva's avatar
Uliana Alekseeva committed
453 454
         CALL qfix(mpi,stars,fi%atoms,fi%sym,fi%vacuum,sphhar,fi%input,fi%cell,fi%oned,overallDen,&
                   fi%noco%l_noco,.TRUE.,l_par=.FALSE.,force_fix=.TRUE.,fix=fix)
455 456
      END IF
#ifdef CPP_MPI
457
      CALL mpi_bc_potden(mpi,stars,sphhar,fi%atoms,fi%input,fi%vacuum,fi%oned,fi%noco,overallDen)
458
#endif
459 460

      ! Calculate Coulomb potential for overall density (+including external potential)
461 462
      CALL overallDen%sum_both_spin()!workden)
      CALL overallVCoul%resetPotDen()
Gregor Michalicek's avatar
Gregor Michalicek committed
463 464
      ALLOCATE(overallVCoul%pw_w(size(overallVCoul%pw,1),size(overallVCoul%pw,2)))
      overallVCoul%pw_w(:,:) = 0.0
465
      CALL vgen_coulomb(1,mpi,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
466
      CALL convol(stars,overallVCoul%pw_w(:,1),overallVCoul%pw(:,1),stars%ufft)   ! Is there a problem with a second spin?!
467
#ifdef CPP_MPI
468
      CALL mpi_bc_potden(mpi,stars,sphhar,fi%atoms,fi%input,fi%vacuum,fi%oned,fi%noco,overallVCoul)
469
#endif
470

471
      overallVCoulSSDen = 0.0
472 473 474
      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
475
            DO iBand = 1, highestState(ikpt,jsp)
476 477 478 479
               ! Read the single-state density from disc
               filename = ''
               WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
               IF (mpi%irank.EQ.0) THEN
480 481
                  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)))
482 483 484
                  CALL singleStateDen%sum_both_spin()!workden)
               END IF
#ifdef CPP_MPI
485
               CALL mpi_bc_potden(mpi,stars,sphhar,fi%atoms,fi%input,fi%vacuum,fi%oned,fi%noco,singleStateDen)
486
#endif
487

488
               ! For each state calculate integral over Coulomb potential times single state density
489
               potDenInt = 0.0
490
               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?!
491 492 493 494
               overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
            END DO
         END DO
      END DO
495

Gregor Michalicek's avatar
Gregor Michalicek committed
496 497 498
      ! Construct exchange matrix diagonal

      exDiag = 0.0
499 500
      DO jspin = 1, fi%input%jspins
         jsp = MERGE(1,jspin,fi%noco%l_noco)
Gregor Michalicek's avatar
Gregor Michalicek committed
501 502
         ! remove weights(wtkpt) in w_iks
         wl_iks = 0.0
503
         DO ikptf = 1, fi%kpts%nkptf
504
            ikpt = fi%kpts%bkp(ikptf)
Gregor Michalicek's avatar
Gregor Michalicek committed
505
            DO iBand=1, results%neig(ikpt,jsp)
506
               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
507
               wl_iks(iBand,ikptf) = sqrt(wl_iks(iBand,ikptf)) ! For Müller functional
508
               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
509 510 511 512 513 514 515 516 517 518
            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
519 520
         results%neig(:,:) = neigTemp(:,:)

521 522 523
         CALL HF_setup(mpdata,fi%hybinp,fi%input,fi%sym,fi%kpts,fi%atoms,mpi,fi%noco,nococonv,&
                       fi%cell,fi%oned,results,jspin,enpara,&
                       hybdat,fi%sym%invs,vTot%mt(:,0,:,:),eig_irr)
Gregor Michalicek's avatar
Gregor Michalicek committed
524

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

527
         DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
528

529
            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
530

531 532 533 534 535 536 537
            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

538
            parent = 0
539
            CALL zMat%init(fi%sym%invs,nbasfcn,fi%input%neig)
Matthias Redies's avatar
Matthias Redies committed
540

541 542
            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)
543

544 545 546
            ALLOCATE (indx_sest(hybdat%nbands(ikpt), hybdat%nbands(ikpt)))
            indx_sest = 0

Matthias Redies's avatar
Matthias Redies committed
547
            call symm_hf_init(fi,ikpt,nsymop,rrot,psym)
Matthias Redies's avatar
Matthias Redies committed
548
            call symm_hf(fi,ikpt,hybdat,eig_irr,mpdata,lapw,nococonv, zMat, c_phase,jspin,&
Matthias Redies's avatar
Matthias Redies committed
549
                         rrot,nsymop,psym,n_q,parent,pointer_EIBZ,nsest,indx_sest)
Gregor Michalicek's avatar
Gregor Michalicek committed
550

551
            exMat%l_real=fi%sym%invs
Matthias Redies's avatar
Matthias Redies committed
552
            CALL exchange_valence_hf(ikpt,fi,zMat, c_phase,mpdata,jspin,hybdat,lapw,&
553
                                     eig_irr,results,pointer_EIBZ,n_q,wl_iks,xcpot,nococonv,stars,nsest,indx_sest,&
Gregor Michalicek's avatar
Gregor Michalicek committed
554
                                     mpi,exMat)
555 556
            CALL exchange_vccv1(ikpt,fi%input,fi%atoms,fi%cell, fi%kpts, fi%sym, fi%noco,nococonv, fi%oned,&
                                mpdata,fi%hybinp,hybdat,jspin,lapw,nsymop,nsest,indx_sest,mpi,&
Matthias Redies's avatar
Matthias Redies committed
557
                                1.0,results,exMat)
Gregor Michalicek's avatar
Gregor Michalicek committed
558

559 560
            DEALLOCATE(indx_sest)

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

563
            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
564

565 566
            CALL olap%alloc(fi%sym%invs,nbasfcn)
            CALL read_olap(olap, fi%kpts%nkpt*(jspin-1)+ikpt, nbasfcn)
Gregor Michalicek's avatar
Gregor Michalicek committed
567

568
            zMat%matsize2 = hybdat%nbands(ikpt) ! reduce "visible matsize" for the following computations
Gregor Michalicek's avatar
Gregor Michalicek committed
569 570 571

            CALL olap%multiply(zMat,trafo)

572
            CALL invtrafo%alloc(olap%l_real,hybdat%nbands(ikpt),nbasfcn)
Gregor Michalicek's avatar
Gregor Michalicek committed
573
            CALL trafo%TRANSPOSE(invtrafo)
Matthias Redies's avatar
Matthias Redies committed
574

575 576
            DO iBand = 1, hybdat%nbands(ikpt)
               DO jBand = 1, iBand-1
Gregor Michalicek's avatar
Gregor Michalicek committed
577
                  IF (exMat%l_real) THEN
578
                     exMat%data_r(iBand,jBand)=exMat%data_r(jBand,iBand)
Gregor Michalicek's avatar
Gregor Michalicek committed
579
                  ELSE
580
                     exMat%data_c(iBand,jBand)=conjg(exMat%data_c(jBand,iBand))
Gregor Michalicek's avatar
Gregor Michalicek committed
581 582 583 584 585 586 587
                  END IF
               END DO
            END DO

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

588
            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
589

590
            IF (.NOT.exMatLAPW%l_real) exMatLAPW%data_c=conjg(exMatLAPW%data_c)
Gregor Michalicek's avatar
Gregor Michalicek committed
591 592 593 594
            zMat%matsize1=MIN(zMat%matsize1,exMatLAPW%matsize2)

            CALL exMatLAPW%multiply(zMat,tmpMat)

Gregor Michalicek's avatar
Gregor Michalicek committed
595
            DO iBand = 1, highestState(ikpt,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
596 597
               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
598
               ELSE
Gregor Michalicek's avatar
Gregor Michalicek committed
599
                  exDiag(iBand,ikpt,jspin) = REAL(dot_product(zMat%data_c(:zMat%matsize1,iband),tmpMat%data_c(:,iband)))
Gregor Michalicek's avatar
Gregor Michalicek committed
600 601
               END IF
            END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
602

603
            !End of workaround for increased functionality of fi%symmetrizeh (call it))
Gregor Michalicek's avatar
Gregor Michalicek committed
604 605 606 607 608 609 610 611

!            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
612 613 614 615 616
         END DO
      END DO

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

Gregor Michalicek's avatar
Gregor Michalicek committed
617
      occSum = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
618
      gradSum = 0.0
619 620 621 622
      DO ispin = 1, fi%input%jspins
         isp = MERGE(1,ispin,fi%noco%l_noco)
!         CALL cdnvalJob%init(mpi,fi%input,fi%kpts,fi%noco,results,isp,fi%banddos=fi%banddos)
         DO ikpt = 1, fi%kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
623
            DO iBand = 1, highestState(ikpt,isp)
624
               occStateI = results%w_iks(iBand,ikpt,isp) / (fi%kpts%wtkpt(ikpt))!*fi%kpts%nkptf)
625
               occStateI = MAX(occStateI,minOcc)
Gregor Michalicek's avatar
Gregor Michalicek committed
626 627 628

               occStateI = MIN(occStateI,1.0-minOcc)

Gregor Michalicek's avatar
Gregor Michalicek committed
629
!               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
630
!               occStateI = cdnvalJob%weights(iBand,ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
631 632
               rdmftFunctionalValue = 0.5*0.5*SQRT(1.0/occStateI) ! for Müller functional derivative

Gregor Michalicek's avatar
Gregor Michalicek committed
633 634 635 636 637 638 639 640 641 642 643 644 645 646 647

               !!! 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)

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

               dEdOcc(iBand,ikpt,isp) = +((spinDegenFac * results%eig(iBand,ikpt,isp)) - vTotSSDen(iBand,ikpt,isp) + &
Gregor Michalicek's avatar
Gregor Michalicek committed
651 652 653 654 655
                                          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)
656
!               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
657

Gregor Michalicek's avatar
Gregor Michalicek committed
658
               WRITE(*,*) 'ENERGY GRADIENT CONTRIBUTIONS'
Gregor Michalicek's avatar
Gregor Michalicek committed
659
               WRITE(*,*) 'ispin, ikpt, iBand, weight', ispin, ikpt, iBand, results%w_iks(iBand,ikpt,isp)
Gregor Michalicek's avatar
Gregor Michalicek committed
660 661 662
               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
663 664
               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
665 666 667
               WRITE(*,*) 'exchangeTerm', exchangeTerm
               WRITE(*,*) 'exDiag(iBand,ikpt,isp)', exDiag(iBand,ikpt,isp)
               WRITE(*,*) 'rdmftFunctionalValue', rdmftFunctionalValue
Gregor Michalicek's avatar
Gregor Michalicek committed
668

Gregor Michalicek's avatar
Gregor Michalicek committed
669
               gradSum = gradSum + dEdOcc(iBand,ikpt,isp) !* results%w_iks(iBand,ikpt,isp))
Gregor Michalicek's avatar
Gregor Michalicek committed
670 671 672
            END DO
         END DO
      END DO
673
      lagrangeMultiplier = -gradSum / numStates ! occSum  !(fi%input%zelec/(2.0/REAL(fi%input%jspins)))
Gregor Michalicek's avatar
Gregor Michalicek committed
674

Gregor Michalicek's avatar
Gregor Michalicek committed
675
   WRITE(*,*) 'lagrangeMultiplier: ', lagrangeMultiplier
Gregor Michalicek's avatar
Gregor Michalicek committed
676

677 678
      ! Optimize occupation numbers

Gregor Michalicek's avatar
Gregor Michalicek committed
679 680
      ALLOCATE (gradient(numStates+1))
      ALLOCATE (parameters(numStates+1))
681
      ALLOCATE (paramGroup(numStates+1))
Gregor Michalicek's avatar
Gregor Michalicek committed
682 683 684 685 686 687 688 689 690 691 692 693
      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
694 695 696
      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
697 698
            DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
               iState = iState + 1
699
               occStateI = results%w_iks(iBand,ikpt,isp) / fi%kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
700

701
               occStateI = MAX(occStateI,minOcc)
Gregor Michalicek's avatar
Gregor Michalicek committed
702 703 704 705
               occStateI = MIN(occStateI,1.0-minOcc)

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

706
               WRITE(7865,'(i7,4f15.10)') iState, occStateI, theta, sin(theta), cos(theta)
Gregor Michalicek's avatar
Gregor Michalicek committed
707 708

!               occStateI = MAX(occStateI,minOcc)
709
               equalityLinCombi(iState) = fi%kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
710 711 712 713 714

!               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
715
               gradient(iState) = dEdOcc(iBand,ikpt,isp) + lagrangeMultiplier
Gregor Michalicek's avatar
Gregor Michalicek committed
716

717
               gradient(numStates+1) = gradient(numStates+1) + occStateI * fi%kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
718
!               parameters(iState) = theta !occStateI
Gregor Michalicek's avatar
Gregor Michalicek committed
719
               parameters(iState) = occStateI
720 721 722 723 724 725 726 727 728 729 730 731 732
               IF(iBand.EQ.1) THEN
                  IF(iState.EQ.1) THEN
                     paramGroup(iState) = 1
                  ELSE
                     paramGroup(iState) = paramGroup(iState-1) + 1
                  END IF
               ELSE
                  IF ((results%eig(iBand,ikpt,jsp) - results%eig(iBand-1,ikpt,jsp)).GT.degenEps) THEN
                     paramGroup(iState) = paramGroup(iState-1) + 1
                  ELSE
                     paramGroup(iState) = paramGroup(iState-1)
                  END IF
               END IF
Gregor Michalicek's avatar
Gregor Michalicek committed
733 734 735
            END DO
         END DO
      END DO
736
      equalityCriterion = fi%input%zelec/(2.0/REAL(fi%input%jspins))
Gregor Michalicek's avatar
Gregor Michalicek committed
737
      gradient(numStates+1) = gradient(numStates+1) - equalityCriterion ! This should actually always be 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
738
      parameters(numStates+1) = lagrangeMultiplier
739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770
      paramGroup(numStates+1) = paramGroup(numStates) + 1

      currentGroupEnd = 0
      lastGroupEnd = 0
      DO iState = 2, numStates + 1
         IF (paramGroup(iState).NE.paramGroup(iState-1)) THEN
            currentGroupEnd = iState - 1
            averageParam = 0.0
            averageGrad = 0.0
            DO jState = lastGroupEnd + 1, currentGroupEnd
               averageParam = averageParam + parameters(jState)
               averageGrad = averageGrad + gradient(jState)
            END DO
            averageParam = averageParam / (currentGroupEnd - lastGroupEnd)
            averageGrad = averageGrad / (currentGroupEnd - lastGroupEnd)
            DO jState = lastGroupEnd + 1, currentGroupEnd
               temp = ABS(parameters(jState) - averageParam) / (ABS(averageParam) + degenEps)
               IF (temp.GT.degenEps) THEN
                  WRITE(*,'(a,i0,a,f15.10,a,f15.10)') 'iState: ', jState, 'average: ', averageParam, 'parameter: ', parameters(jState)
                  CALL juDFT_error('parameter difference in paramGroup is very large (rdmft-1)')
               END IF
               temp = ABS(gradient(jState) - averageGrad) / (ABS(averageGrad) + degenEps)
               IF (temp.GT.degenEps) THEN
                  WRITE(*,'(a,i0,a,f15.10,a,f15.10)') 'iState: ', jState, 'average: ', averageGrad, 'gradient: ', gradient(jState)
                  CALL juDFT_error('Gradient difference in paramGroup is very large (rdmft-1)')
               END IF
               parameters(jState) = averageParam
               gradient(jState) = averageGrad
            END DO
            lastGroupEnd = currentGroupEnd
         END IF
      END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
771

Gregor Michalicek's avatar
Gregor Michalicek committed
772 773