rdmft.F90 35.9 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,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars,vacuum,&
Matthias Redies's avatar
Matthias Redies committed
12
                 sphhar,sym,field,vTot,vCoul,oneD,noco,xcpot,mpinp,mpdata,hybinp, hybdat,&
13
                 gfinp,hub1inp,results,coreSpecInput,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
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
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
38
   USE m_symmetrizeh
Gregor Michalicek's avatar
Gregor Michalicek committed
39 40
   USE m_bfgs_b2

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

   IMPLICIT NONE

   TYPE(t_mpi),           INTENT(IN)    :: mpi
   TYPE(t_input),         INTENT(IN)    :: input
   TYPE(t_kpts),          INTENT(IN)    :: kpts
   TYPE(t_banddos),       INTENT(IN)    :: banddos
Gregor Michalicek's avatar
Gregor Michalicek committed
52
   TYPE(t_sliceplot),     INTENT(IN)    :: sliceplot
53 54
   TYPE(t_cell),          INTENT(IN)    :: cell
   TYPE(t_atoms),         INTENT(IN)    :: atoms
Gregor Michalicek's avatar
Gregor Michalicek committed
55
   TYPE(t_enpara),        INTENT(INOUT) :: enpara
56 57 58 59
   TYPE(t_stars),         INTENT(IN)    :: stars
   TYPE(t_vacuum),        INTENT(IN)    :: vacuum
   TYPE(t_sphhar),        INTENT(IN)    :: sphhar
   TYPE(t_sym),           INTENT(IN)    :: sym
60
   TYPE(t_field),         INTENT(INOUT) :: field
61
   TYPE(t_potden),        INTENT(INOUT) :: vTot
62
   TYPE(t_potden),        INTENT(INOUT) :: vCoul
63
   TYPE(t_oneD),          INTENT(IN)    :: oneD
Gregor Michalicek's avatar
Gregor Michalicek committed
64
   TYPE(t_noco),          INTENT(INOUT) :: noco
Matthias Redies's avatar
Matthias Redies committed
65
   TYPE(t_xcpot_inbuild), INTENT(INOUT) :: xcpot
Matthias Redies's avatar
Matthias Redies committed
66
   TYPE(t_mpinp),         INTENT(IN)    :: mpinp
67 68 69
   TYPE(t_mpdata),        intent(inout) :: mpdata
   TYPE(t_hybinp),        INTENT(IN)    :: hybinp
   TYPE(t_hybdat),        INTENT(INOUT) :: hybdat
70 71
   TYPE(t_gfinp),         INTENT(IN)    :: gfinp
   TYPE(t_hub1inp),       INTENT(IN)    :: hub1inp
72
   TYPE(t_results),       INTENT(INOUT) :: results
Gregor Michalicek's avatar
Gregor Michalicek committed
73 74
   TYPE(t_coreSpecInput), INTENT(IN)    :: coreSpecInput
   TYPE(t_potden),        INTENT(INOUT) :: outDen
75 76

   INTEGER,               INTENT(IN)    :: eig_id
Gregor Michalicek's avatar
Gregor Michalicek committed
77
   INTEGER,               INTENT(IN)    :: archiveType
Matthias Redies's avatar
Matthias Redies committed
78
   TYPE(t_potden)                       :: EnergyDen
79

80
#ifndef CPP_OLDINTEL
81
   TYPE(t_cdnvalJob)                    :: cdnvalJob
82
   TYPE(t_potden)                       :: singleStateDen, overallDen, overallVCoul, vTotTemp
83 84 85
   TYPE(t_regionCharges)                :: regCharges
   TYPE(t_dos)                          :: dos
   TYPE(t_moments)                      :: moments
86
   TYPE(t_mat)                          :: exMat, zMat, olap, trafo, invtrafo, tmpMat, exMatLAPW
Gregor Michalicek's avatar
Gregor Michalicek committed
87
   TYPE(t_lapw)                         :: lapw
Gregor Michalicek's avatar
Gregor Michalicek committed
88
   INTEGER                              :: ikpt, ikpt_i, iBand, jkpt, jBand, iAtom, i, na, itype, lh, j
89
   INTEGER                              :: jspin, jspmax, jsp, isp, ispin, nbasfcn, nbands
Gregor Michalicek's avatar
Gregor Michalicek committed
90
   INTEGER                              :: nsymop, nkpt_EIBZ, ikptf, iterHF, mnobd
Gregor Michalicek's avatar
Gregor Michalicek committed
91 92
   INTEGER                              :: iState, iStep, numStates, numRelevantStates, convIter
   INTEGER                              :: maxHistoryLength
93 94
   REAL                                 :: fix, potDenInt, fermiEnergyTemp, spinDegenFac
   REAL                                 :: rdmftFunctionalValue, occStateI, gradSum
95
   REAL                                 :: exchangeTerm, lagrangeMultiplier, equalityCriterion
Gregor Michalicek's avatar
Gregor Michalicek committed
96 97 98
   REAL                                 :: mixParam, rdmftEnergy, occSum
   REAL                                 :: sumOcc, addCharge, subCharge, addChargeWeight, subChargeWeight
   REAL                                 :: rhs, totz, theta
Gregor Michalicek's avatar
Gregor Michalicek committed
99
   REAL, PARAMETER                      :: degenEps = 0.00001
Gregor Michalicek's avatar
Gregor Michalicek committed
100 101
   REAL, PARAMETER                      :: convCrit = 5.0e-6
   REAL, PARAMETER                      :: minOcc = 1.0e-13
Gregor Michalicek's avatar
Gregor Michalicek committed
102
   LOGICAL                              :: converged, l_qfix, l_restart, l_zref
103
   CHARACTER(LEN=20)                    :: filename
104

105 106
   INTEGER                              :: nsest(input%neig) ! probably too large
   INTEGER                              :: indx_sest(input%neig,input%neig) ! probably too large
Gregor Michalicek's avatar
Gregor Michalicek committed
107 108 109 110 111 112
   INTEGER                              :: rrot(3,3,sym%nsym)
   INTEGER                              :: psym(sym%nsym) ! Note: psym is only filled up to index nsymop
   INTEGER                              :: lowestState(kpts%nkpt,input%jspins)
   INTEGER                              :: highestState(kpts%nkpt,input%jspins)
   INTEGER                              :: neigTemp(kpts%nkpt,input%jspins)

113
   REAL                                 :: wl_iks(input%neig,kpts%nkptf)
Gregor Michalicek's avatar
Gregor Michalicek committed
114

Gregor Michalicek's avatar
Gregor Michalicek committed
115 116
   REAL                                 :: vmd(atoms%ntype), zintn_r(atoms%ntype), dpj(atoms%jmtd), mt(atoms%jmtd,atoms%ntype)

117 118
   REAL, ALLOCATABLE                    :: overallVCoulSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vTotSSDen(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
119 120
   REAL, ALLOCATABLE                    :: dEdOcc(:,:,:)

Gregor Michalicek's avatar
Gregor Michalicek committed
121 122 123 124
   REAL, ALLOCATABLE                    :: zintn_rSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vmdSSDen(:,:,:)


Gregor Michalicek's avatar
Gregor Michalicek committed
125 126 127 128 129 130 131 132 133 134 135 136
   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(:)
137

Gregor Michalicek's avatar
Gregor Michalicek committed
138 139 140 141 142 143 144
   INTEGER, ALLOCATABLE                 :: parent(:)
   INTEGER, ALLOCATABLE                 :: pointer_EIBZ(:)
   INTEGER, ALLOCATABLE                 :: n_q(:)

   LOGICAL, ALLOCATABLE                 :: enabledConstraints(:)

#endif
145

146
#ifndef CPP_OLDINTEL
147 148 149

   WRITE(*,*) 'entered RDMFT subroutine'

150
   ! General initializations
Gregor Michalicek's avatar
Gregor Michalicek committed
151
   mixParam = 0.0001
152
   lagrangeMultiplier = 0.1 !results%ef
153
   spinDegenFac = 2.0 / 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 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210

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

   ! Determine which states have to be considered
   lowestState(:,:) = -1
   highestState(:,:) = -1
   DO jspin = 1, input%jspins
      jsp = MERGE(1,jspin,noco%l_noco)
      DO ikpt = 1, kpts%nkpt
         ! determine lowest state
         DO iBand = 1, results%neig(ikpt,jsp)
            IF((results%w_iks(iBand,ikpt,jspin) / kpts%wtkpt(ikpt)).LT.(1.0-input%rdmftOccEps)) THEN
               lowestState(ikpt,jspin) = iBand
               EXIT
            END IF
         END DO
         lowestState(ikpt,jspin) = lowestState(ikpt,jspin) - input%rdmftStatesBelow
         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
            IF((results%w_iks(iBand,ikpt,jspin) / kpts%wtkpt(ikpt)).GT.(0.0+input%rdmftOccEps)) THEN
               highestState(ikpt,jspin) = iBand
               EXIT
            END IF
         END DO
         highestState(ikpt,jspin) = highestState(ikpt,jspin) + input%rdmftStatesAbove
         IF((results%neig(ikpt,jsp)-1).LT.highestState(ikpt,jspin)) THEN
            WRITE(6,*) 'Error: Not enough states calculated:'
            WRITE(6,*) 'ikpt, jsp: ', ikpt, jsp
            WRITE(6,*) 'highestState(ikpt,jspin): ', highestState(ikpt,jspin)
            WRITE(6,*) 'results%neig(ikpt,jsp): ', results%neig(ikpt,jsp)
            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
            WRITE(6,*) 'Error: Highest state is degenerate:'
            WRITE(6,*) 'ikpt, jsp: ', ikpt, jsp
            WRITE(6,*) 'highestState(ikpt,jspin): ', highestState(ikpt,jspin)
            WRITE(6,*) 'results%eig(highestState(ikpt,jspin),ikpt,jsp): ', results%eig(highestState(ikpt,jspin),ikpt,jsp)
            WRITE(6,*) 'results%eig(highestState(ikpt,jspin)-1,ikpt,jsp): ', results%eig(highestState(ikpt,jspin)-1,ikpt,jsp)
            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 214
   IF (ANY(results%w_iksRDMFT(:,:,:).NE.0.0)) THEN
      results%w_iks(:,:,:) = results%w_iksRDMFT(:,:,:)
   END IF

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

   ! Some more initializations

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

265 266
   ALLOCATE(overallVCoulSSDen(MAXVAL(results%neig(1:kpts%nkpt,1:input%jspins)),kpts%nkpt,input%jspins))
   ALLOCATE(vTotSSDen(MAXVAL(results%neig(1:kpts%nkpt,1:input%jspins)),kpts%nkpt,input%jspins))
Gregor Michalicek's avatar
Gregor Michalicek committed
267
   ALLOCATE(dEdOcc(MAXVAL(results%neig(1:kpts%nkpt,1:input%jspins)),kpts%nkpt,input%jspins))
268

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

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

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

286 287 288 289 290 291 292 293 294 295 296 297 298 299
   CALL vTotTemp%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTTOT)
   vTotTemp = vTot

   DO jsp = 1,SIZE(vTot%mt,4)
      DO iAtom = 1,atoms%ntype
         vTotTemp%mt(:atoms%jri(iAtom),0,iAtom,jsp)  = sfp_const * vTot%mt(:atoms%jri(iAtom),0,iAtom,jsp) / atoms%rmsh(:atoms%jri(iAtom),iAtom)
      END DO
   END DO

   vCoul%pw_w = CMPLX(0.0,0.0)
   DO jspin = 1, input%jspins
      CALL convol(stars,vCoul%pw_w(:,jspin),vCoul%pw(:,jspin),stars%ufft)
   END DO

300 301
   vTotSSDen = 0.0

302
   ! Calculate all single state densities
303

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

      CALL cdnvalJob%init(mpi,input,kpts,noco,results,jsp)

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

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

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

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

            mt(:,:) = 0.0
            DO iType = 1, atoms%ntype
               DO i = 1, atoms%jri(iType)
                  mt(i,iType) = singleStateDen%mt(i,0,iType,jsp)
               END DO
Matthias Redies's avatar
Matthias Redies committed
350

Gregor Michalicek's avatar
Gregor Michalicek committed
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367
               DO j = 1,atoms%jri(iType)
                  dpj(j) = mt(j,iType)/atoms%rmsh(j,iType)
               END DO
               CALL intgr3(dpj,atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),rhs)

               zintn_r(iType) = atoms%neq(iType)*atoms%zatom(iType)*sfp_const*rhs/2.0
               zintn_rSSDen(iBand,ikpt,jsp) = zintn_rSSDen(iBand,ikpt,jsp) + zintn_r(iType)

               CALL intgr3(mt(1,iType),atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),totz)

!               vmd(iType) = atoms%rmt(iType)*vCoul%mt(atoms%jri(iType),0,iType,1)/sfp_const + atoms%zatom(iType) - totz*sfp_const
               vmd(iType) = -totz*sfp_const
               vmd(iType) = -atoms%neq(iType)*atoms%zatom(iType)*vmd(iType)/ (2.0*atoms%rmt(iType))

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

368 369 370 371
         END DO
      END DO
   END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
372 373
   ! Initializations for exchange contributions

374 375
   WRITE(*,*) 'RDMFT: HF initializations start'

376 377 378
   IF(ALLOCATED(hybdat%ne_eig)) DEALLOCATE(hybdat%ne_eig)
   IF(ALLOCATED(hybdat%nbands)) DEALLOCATE(hybdat%nbands)
   IF(ALLOCATED(hybdat%nobd)) DEALLOCATE(hybdat%nobd)
379
   IF(ALLOCATED(hybdat%nbasm)) DEALLOCATE(hybdat%nbasm)
380
   IF(ALLOCATED(hybdat%div_vv)) DEALLOCATE(hybdat%div_vv)
381
   ALLOCATE(hybdat%ne_eig(kpts%nkpt),hybdat%nbands(kpts%nkpt),hybdat%nobd(kpts%nkptf,input%jspins))
382
   ALLOCATE(hybdat%nbasm(kpts%nkptf))
383
   ALLOCATE(hybdat%div_vv(input%neig,kpts%nkpt,input%jspins))
384

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

389
!   CALL open_hybinp_io1(sym%invs)
Gregor Michalicek's avatar
Gregor Michalicek committed
390

Matthias Redies's avatar
Matthias Redies committed
391
   CALL mixedbasis(atoms,kpts,input,cell,xcpot,mpinp,mpdata,hybinp, hybdat,enpara,mpi,vTot, iterHF)
Gregor Michalicek's avatar
Gregor Michalicek committed
392

393
   CALL open_hybinp_io2(mpdata, hybinp,hybdat,input,atoms,sym%invs)
Gregor Michalicek's avatar
Gregor Michalicek committed
394

395
   CALL coulombmatrix(mpi,atoms,kpts,cell,sym,mpdata,hybinp,hybdat,xcpot)
Gregor Michalicek's avatar
Gregor Michalicek committed
396

397
   CALL hf_init(mpdata,hybinp,atoms,input,hybdat)
Gregor Michalicek's avatar
Gregor Michalicek committed
398

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

Gregor Michalicek's avatar
Gregor Michalicek committed
401 402 403
   maxHistoryLength = 7*numStates
   maxHistoryLength = 5*numStates

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

Gregor Michalicek's avatar
Gregor Michalicek committed
418 419
   convIter = 0

420
   converged = .FALSE.
421 422
   DO WHILE (.NOT.converged)

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

      DO jspin = 1, input%jspins
         DO ikpt = 1,kpts%nkpt
            WRITE(*,*) 'jspin, ikpt: ', jspin, ikpt
Gregor Michalicek's avatar
Gregor Michalicek committed
430
            WRITE(*,'(10f11.6)') results%w_iks(1:10,ikpt,jspin)
Gregor Michalicek's avatar
Gregor Michalicek committed
431 432 433
         END DO
      END DO

434
      ! Calculate overall density with current occupation numbers (don't forget core electron density)
Gregor Michalicek's avatar
Gregor Michalicek committed
435
      CALL overallDen%resetPotDen()
436 437
      jspmax = input%jspins
      IF (noco%l_mperp) jspmax = 1
Gregor Michalicek's avatar
Gregor Michalicek committed
438

439
      DO jspin = 1,jspmax
440
         CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
441
         CALL cdnval(eig_id,mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,stars,vacuum,&
442 443
                     sphhar,sym,vTot,oneD,cdnvalJob,overallDen,regCharges,dos,results,moments,&
                     gfinp,hub1inp)
444
      END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
445

446
      CALL cdncore(mpi,oneD,input,vacuum,noco,sym,&
447
                   stars,cell,sphhar,atoms,vTot,overallDen,moments,results)
448
      IF (mpi%irank.EQ.0) THEN
449
         CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,overallDen,noco%l_noco,.TRUE.,.true.,fix)
450 451 452 453
      END IF
#ifdef CPP_MPI
      CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallDen)
#endif
454 455

      ! Calculate Coulomb potential for overall density (+including external potential)
456 457
      CALL overallDen%sum_both_spin()!workden)
      CALL overallVCoul%resetPotDen()
Gregor Michalicek's avatar
Gregor Michalicek committed
458 459
      ALLOCATE(overallVCoul%pw_w(size(overallVCoul%pw,1),size(overallVCoul%pw,2)))
      overallVCoul%pw_w(:,:) = 0.0
460
      CALL vgen_coulomb(1,mpi,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,.FALSE.,overallDen,overallVCoul)
Gregor Michalicek's avatar
Gregor Michalicek committed
461
      CALL convol(stars,overallVCoul%pw_w(:,1),overallVCoul%pw(:,1),stars%ufft)   ! Is there a problem with a second spin?!
462 463 464
#ifdef CPP_MPI
      CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallVCoul)
#endif
465

466
      overallVCoulSSDen = 0.0
467 468 469
      DO jspin = 1, input%jspins
         jsp = MERGE(1,jspin,noco%l_noco)
         DO ikpt = 1, kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
470
            DO iBand = 1, highestState(ikpt,jsp)
471 472 473 474
               ! 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
475
                  CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,&
476 477 478 479 480 481
                                   CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,singleStateDen,TRIM(ADJUSTL(filename)))
                  CALL singleStateDen%sum_both_spin()!workden)
               END IF
#ifdef CPP_MPI
               CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,singleStateDen)
#endif
482

483
               ! For each state calculate integral over Coulomb potential times single state density
484
               potDenInt = 0.0
485
               CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vCoul,singleStateDen,potDenInt) ! Is there a problem with a second spin?!
486 487 488 489
               overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
            END DO
         END DO
      END DO
490

Gregor Michalicek's avatar
Gregor Michalicek committed
491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511
      ! Construct exchange matrix diagonal

      exDiag = 0.0
      DO jspin = 1, input%jspins
         jsp = MERGE(1,jspin,noco%l_noco)
         ! remove weights(wtkpt) in w_iks
         wl_iks = 0.0
         DO ikptf=1,kpts%nkptf
            ikpt = kpts%bkp(ikptf)
            DO iBand=1, results%neig(ikpt,jsp)
               wl_iks(iBand,ikptf) = results%w_iks(iBand,ikpt,jspin) / (kpts%wtkpt(ikpt))!*kpts%nkptf) Last term to be included after applying functional
               wl_iks(iBand,ikptf) = sqrt(wl_iks(iBand,ikptf)) ! For Müller functional
               wl_iks(iBand,ikptf) = wl_iks(iBand,ikptf) / kpts%nkptf ! This is the last part of two lines above
            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)
512

513
         call mpdata%free()
514

Gregor Michalicek's avatar
Gregor Michalicek committed
515 516
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE (hybdat%nindxp1)

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

519
         CALL HF_setup(mpdata,hybinp,input,sym,kpts,atoms,mpi,noco,&
Matthias Redies's avatar
Matthias Redies committed
520
                       cell,oneD,results,jspin,enpara,eig_id,&
521
                       hybdat,sym%invs,vTot%mt(:,0,:,:),eig_irr)
Gregor Michalicek's avatar
Gregor Michalicek committed
522

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

525
         mnobd = MAXVAL(hybdat%nobd(:,jsp))
Gregor Michalicek's avatar
Gregor Michalicek committed
526 527 528 529 530 531 532

         DO ikpt = 1,kpts%nkpt

            CALL lapw%init(input,noco,kpts,atoms,sym,ikpt,cell,l_zref)

            parent = 0
            CALL symm_hf_init(sym,kpts,ikpt,nsymop,rrot,psym)
533
            CALL symm_hf(kpts,ikpt,sym,hybdat,eig_irr,input,atoms,mpdata,hybinp,cell,lapw,jspin,&
Gregor Michalicek's avatar
Gregor Michalicek committed
534 535 536
                         rrot,nsymop,psym,nkpt_EIBZ,n_q,parent,pointer_EIBZ,nsest,indx_sest)

            exMat%l_real=sym%invs
537
            CALL exchange_valence_hf(ikpt,kpts,nkpt_EIBZ, sym,atoms,mpdata,hybinp,cell,input,jspin,hybdat,mnobd,lapw,&
538
                                     eig_irr,results,pointer_EIBZ,n_q,wl_iks,xcpot,noco,nsest,indx_sest,&
Gregor Michalicek's avatar
Gregor Michalicek committed
539
                                     mpi,exMat)
540
            CALL exchange_vccv1(ikpt,input,atoms,mpdata,hybinp,hybdat,jspin,lapw,nsymop,nsest,indx_sest,mpi,1.0,results,exMat)
541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562

            !Start of workaround for increased functionality of symmetrizeh (call it))

            nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)

            CALL olap%alloc(sym%invs,nbasfcn)
            CALL read_olap(olap, kpts%nkpt*(jspin-1)+ikpt)
            IF (olap%l_real) THEN
               DO i = 1, nbasfcn
                  DO j = 1, i
                     olap%data_r(i,j) = olap%data_r(j,i)
                  END DO
               END DO
            ELSE
               DO i = 1, nbasfcn
                  DO j = 1, i
                     olap%data_c(i,j) = CONJG(olap%data_c(j,i))
                  END DO
               END DO
               olap%data_c = conjg(olap%data_c)
            END IF

563
            CALL zMat%init(olap%l_real,nbasfcn,input%neig)
564

565
            CALL read_eig(eig_id,ikpt,jspin,list=[(i,i=1,hybdat%nbands(ikpt))],neig=nbands,zmat=zMat)
566 567

!            CALL read_z(zMat,kpts%nkpt*(jspin-1)+ikpt)
568
            zMat%matsize2 = hybdat%nbands(ikpt) ! reduce "visible matsize" for the following computations
569 570 571

            CALL olap%multiply(zMat,trafo)

572
            CALL invtrafo%alloc(olap%l_real,hybdat%nbands(ikpt),nbasfcn)
573 574 575
            CALL trafo%TRANSPOSE(invtrafo)
            IF(.NOT.invtrafo%l_real) invtrafo%data_c = CONJG(invtrafo%data_c)

576
            DO i = 1, hybdat%nbands(ikpt)
577 578 579 580 581 582 583 584 585 586 587 588
               DO j = 1, i-1
                  IF (exMat%l_real) THEN
                     exMat%data_r(i,j)=exMat%data_r(j,i)
                  ELSE
                     exMat%data_c(i,j)=conjg(exMat%data_c(j,i))
                  END IF
               END DO
            END DO

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

Matthias Redies's avatar
Matthias Redies committed
589
            CALL symmetrizeh(atoms,kpts%bkf(:,ikpt),jspin,lapw,sym,hybdat%kveclo_eig,cell,nsymop,psym,exMatLAPW)
590

591
            IF (.NOT.exMatLAPW%l_real) exMatLAPW%data_c=conjg(exMatLAPW%data_c)
592 593 594 595
            zMat%matsize1=MIN(zMat%matsize1,exMatLAPW%matsize2)

            CALL exMatLAPW%multiply(zMat,tmpMat)

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

            !End of workaround for increased functionality of symmetrizeh (call it))

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

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

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

               occStateI = MIN(occStateI,1.0-minOcc)

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

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

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

649 650 651
               exchangeTerm = - rdmftFunctionalValue * exDiag(iBand,ikpt,isp) * kpts%wtkpt(ikpt) * spinDegenFac !*kpts%nkptf

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

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

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

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

678 679
      ! Optimize occupation numbers

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

715
               gradient(iState) = dEdOcc(iBand,ikpt,isp) + lagrangeMultiplier
Gregor Michalicek's avatar
Gregor Michalicek committed
716

Gregor Michalicek's avatar
Gregor Michalicek committed
717
               gradient(numStates+1) = gradient(numStates+1) + occStateI * kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
718
!               parameters(iState) = theta !occStateI
Gregor Michalicek's avatar
Gregor Michalicek committed
719 720 721 722 723 724
               parameters(iState) = occStateI
            END DO
         END DO
      END DO
      equalityCriterion = input%zelec/(2.0/REAL(input%jspins))
      gradient(numStates+1) = gradient(numStates+1) - equalityCriterion ! This should actually always be 0.0
725
      parameters(numStates+1) = lagrangeMultiplier
Gregor Michalicek's avatar
Gregor Michalicek committed
726

Gregor Michalicek's avatar
Gregor Michalicek committed
727 728
      WRITE(*,*) 'gradient(numStates+1): ', gradient(numStates+1)

Gregor Michalicek's avatar
Gregor Michalicek committed
729
      mixParam = 0.01 / MAXVAL(ABS(gradient(:numStates)))
Gregor Michalicek's avatar
Gregor Michalicek committed
730
!      mixParam = MIN(0.0002,mixParam)
Gregor Michalicek's avatar
Gregor Michalicek committed
731 732 733 734 735 736
      WRITE(*,*) 'mixParam: ', mixParam

      CALL bfgs_b2(numStates+1,gradient,lastGradient,minConstraints,maxConstraints,enabledConstraints,parameters,&
                   lastParameters,equalityLinCombi,equalityCriterion,maxHistoryLength,paramCorrections,&
                   gradientCorrections,iStep,mixParam,converged,convCrit)

Gregor Michalicek's avatar
Gregor Michalicek committed
737
      WRITE(3555,*) 'Occupation numbers:'
Gregor Michalicek's avatar
Gregor Michalicek committed
738 739 740 741 742 743
      iState = 0
      DO ispin = 1, input%jspins
         isp = MERGE(1,ispin,noco%l_noco)
         DO ikpt = 1, kpts%nkpt
            DO iBand = lowestState(ikpt,isp), highestState(ikpt,isp)
               iState = iState + 1
Gregor Michalicek's avatar
Gregor Michalicek committed
744 745 746 747
!               parameters(iState) = (SIN(parameters(iState)*0.5*pi_const))**2.0
               WRITE(3555,'(3i7,f15.10)') iBand, ikpt, isp, parameters(iState)
               results%w_iks(iBand,ikpt,isp) = parameters(iState) * kpts%wtkpt(ikpt)
!               results%w_iks(iBand,ikpt,isp) = MERGE(parameters(iState) * kpts%wtkpt(ikpt),0.0,parameters(iState).GT.minOcc)
Gregor Michalicek's avatar
Gregor Michalicek committed
748 749 750 751
            END DO
         END DO
      END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
752 753
      WRITE(3555,'(a,f15.10)') 'total occupation: ', SUM(parameters(:numStates))

Gregor Michalicek's avatar
Gregor Michalicek committed
754 755 756
      DEALLOCATE (enabledConstraints,maxConstraints,minConstraints)
      DEALLOCATE (parameters,gradient,equalityLinCombi)

757 758 759

   END DO ! WHILE (.NOT.converged)

Gregor Michalicek's avatar
Gregor Michalicek committed
760 761
   WRITE(2503,*) 'convIter: ', convIter

762
   WRITE(*,*) 'RDMFT: convergence loop end'
Gregor Michalicek's avatar
Gregor Michalicek committed
763

764
   hybdat%l_calhf = .FALSE.
Gregor Michalicek's avatar
Gregor Michalicek committed
765

766 767
   ! Calculate final overall density

Gregor Michalicek's avatar
Gregor Michalicek committed
768 769
   !I think we need most of cdngen at this place so I just use cdngen
   CALL outDen%resetPotDen()
770
   CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,kpts,atoms,sphhar,stars,sym,gfinp,hub1inp,&
Matthias Redies's avatar
Matthias Redies committed
771
               enpara,cell,noco,vTot,results,oneD,coreSpecInput,archiveType,xcpot,outDen, EnergyDen)
Gregor Michalicek's avatar
Gregor Michalicek committed
772 773 774 775 776 777 778 779 780

   ! Calculate RDMFT energy
   rdmftEnergy = 0.0
   DO ispin = 1, input%jspins
      isp = MERGE(1,ispin,noco%l_noco)
!      CALL cdnvalJob%init(mpi,input,kpts,noco,results,isp,banddos=banddos)
      DO ikpt = 1, kpts%nkpt
         DO iBand = 1, highestState(ikpt,isp)
            occStateI = results%w_iks(iBand,ikpt,isp) / (kpts%wtkpt(ikpt))!*kpts%nkptf)
781
            rdmftFunctionalValue = 0.5*SQRT(occStateI) ! for Müller functional
Gregor Michalicek's avatar
Gregor Michalicek committed
782

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

785 786
            rdmftEnergy = rdmftEnergy + exchangeTerm + &
                          occStateI * ((spinDegenFac*results%eig(iBand,ikpt,isp)) - vTotSSDen(iBand,ikpt,isp) + &
Gregor Michalicek's avatar
Gregor Michalicek committed
787 788 789 790 791 792 793 794 795 796 797
                                       0.5*overallVCoulSSDen(iBand,ikpt,isp))

               WRITE(2505,*) 'ENERGY CONTRIBUTIONS'
               WRITE(2505,*) 'ispin, ikpt, iBand, weight', ispin, ikpt, iBand, results%w_iks(iBand,ikpt,isp)
               WRITE(2505,*) 'results%eig(iBand,ikpt,isp)', occStateI * spinDegenFac * results%eig(iBand,ikpt,isp)
               WRITE(2505,*) 'vTotSSDen(iBand,ikpt,isp)', - occStateI * vTotSSDen(iBand,ikpt,isp)
               WRITE(2505,*) 'overallVCoulSSDen(iBand,ikpt,isp)', occStateI * overallVCoulSSDen(iBand,ikpt,isp)
               WRITE(2505,*) 'exchangeTerm', exchangeTerm
               WRITE(2505,*) 'exDiag(iBand,ikpt,isp)', exDiag(iBand,ikpt,isp)
               WRITE(2505,*) 'rdmftFunctionalValue', rdmftFunctionalValue

Gregor Michalicek's avatar
Gregor Michalicek committed
798 799 800 801
         END DO
      END DO
   END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
802
   results%w_iksRDMFT(:,:,:) = results%w_iks(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
803 804
   results%neig(:,:) = neigTemp(:,:)

Gregor Michalicek's avatar
Gregor Michalicek committed
805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840

   ! Madelung term (taken from totale):

   mt=0.0
   DO iType = 1, atoms%ntype
      DO i = 1, atoms%jri(iType)
         mt(i,iType) = outDen%mt(i,0,iType,1) + outDen%mt(i,0,iType,input%jspins)
      END DO
   END DO
   IF (input%jspins.EQ.1) mt=mt/2.0 !we just added the same value twice

   DO iType = 1, atoms%ntype
      DO j = 1,atoms%jri(iType)
         dpj(j) = mt(j,iType)/atoms%rmsh(j,iType)
      END DO
      CALL intgr3(dpj,atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),rhs)

      zintn_r(iType) = atoms%neq(iType)*atoms%zatom(iType)*sfp_const*rhs/2.
      rdmftEnergy = rdmftEnergy - zintn_r(iType)

      CALL intgr3(mt(1,iType),atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),totz)

      vmd(iType) = atoms%rmt(iType)*vCoul%mt(atoms%jri(iType),0,iType,1)/sfp_const + atoms%zatom(iType) - totz*sfp_const
      vmd(iType) = -atoms%neq(iType)*atoms%zatom(iType)*vmd(iType)/ (2.0*atoms%rmt(iType))

      rdmftEnergy = rdmftEnergy + vmd(iType)

      WRITE(2505,*) '======================================='
      WRITE(2505,*) 'iType: ', iType
      WRITE(2505,*) 'zintn_r(iType): ', zintn_r(iType)
      WRITE(2505,*) 'vmd(iType): ', vmd(iType)
      WRITE(2505,*) '======================================='

   END DO


Gregor Michalicek's avatar
Gregor Michalicek committed
841 842
   WRITE(6,'(a,f20.10,a)') 'RDMFT energy: ', rdmftEnergy, ' Htr'

Gregor Michalicek's avatar
Gregor Michalicek committed
843 844 845 846 847 848 849
   WRITE(2505,*) '======================================='
   WRITE(2505,*) 'convIter: ', convIter
   WRITE(2505,'(a,f20.10,a)') 'RDMFT energy: ', rdmftEnergy, ' Htr'
   WRITE(2505,*) '======================================='
   WRITE(2505,*) '======================================='
   WRITE(2505,*) '======================================='

850
#endif
851 852 853
END SUBROUTINE rdmft

END MODULE m_rdmft