rdmft.F90 29.7 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

Gregor Michalicek's avatar
Gregor Michalicek committed
11
SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars,vacuum,dimension,&
Gregor Michalicek's avatar
Gregor Michalicek committed
12
                 sphhar,sym,field,vTot,vCoul,oneD,noco,xcpot,hybrid,results,coreSpecInput,archiveType,outDen)
13 14 15

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

   USE m_mixedbasis
   USE m_coulombmatrix
   USE m_hf_init
   USE m_hf_setup
   USE m_io_hybrid
   USE m_symm_hf
   USE m_exchange_valence_hf
   USE m_exchange_core
Gregor Michalicek's avatar
Gregor Michalicek committed
36
   USE m_symmetrizeh
Gregor Michalicek's avatar
Gregor Michalicek committed
37 38
   USE m_bfgs_b2

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

   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
50
   TYPE(t_sliceplot),     INTENT(IN)    :: sliceplot
51 52
   TYPE(t_cell),          INTENT(IN)    :: cell
   TYPE(t_atoms),         INTENT(IN)    :: atoms
Gregor Michalicek's avatar
Gregor Michalicek committed
53
   TYPE(t_enpara),        INTENT(INOUT) :: enpara
54 55 56 57 58
   TYPE(t_stars),         INTENT(IN)    :: stars
   TYPE(t_vacuum),        INTENT(IN)    :: vacuum
   TYPE(t_dimension),     INTENT(IN)    :: dimension
   TYPE(t_sphhar),        INTENT(IN)    :: sphhar
   TYPE(t_sym),           INTENT(IN)    :: sym
59
   TYPE(t_field),         INTENT(INOUT) :: field
60
   TYPE(t_potden),        INTENT(INOUT) :: vTot
Gregor Michalicek's avatar
Gregor Michalicek committed
61
   TYPE(t_potden),        INTENT(INOUT) :: vCoul
62
   TYPE(t_oneD),          INTENT(IN)    :: oneD
Gregor Michalicek's avatar
Gregor Michalicek committed
63
   TYPE(t_noco),          INTENT(INOUT) :: noco
Matthias Redies's avatar
merge  
Matthias Redies committed
64
   TYPE(t_xcpot_inbuild), INTENT(INOUT) :: xcpot
Gregor Michalicek's avatar
Gregor Michalicek committed
65
   TYPE(t_hybrid),        INTENT(INOUT) :: hybrid
66
   TYPE(t_results),       INTENT(INOUT) :: results
Gregor Michalicek's avatar
Gregor Michalicek committed
67 68
   TYPE(t_coreSpecInput), INTENT(IN)    :: coreSpecInput
   TYPE(t_potden),        INTENT(INOUT) :: outDen
69 70

   INTEGER,               INTENT(IN)    :: eig_id
Gregor Michalicek's avatar
Gregor Michalicek committed
71
   INTEGER,               INTENT(IN)    :: archiveType
Matthias Redies's avatar
merge  
Matthias Redies committed
72
   TYPE(t_potden)                       :: EnergyDen
73

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

Gregor Michalicek's avatar
Gregor Michalicek committed
98 99 100 101 102 103 104 105 106 107
   INTEGER                              :: nsest(dimension%neigd) ! probably too large
   INTEGER                              :: indx_sest(dimension%neigd,dimension%neigd) ! probably too large
   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)

   REAL                                 :: wl_iks(dimension%neigd,kpts%nkptf)

108 109
   REAL, ALLOCATABLE                    :: overallVCoulSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vTotSSDen(:,:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
110 111 112 113 114 115 116 117 118 119 120 121 122 123
   REAL, ALLOCATABLE                    :: dEdOcc(:,:,:)

   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(:)
124

Gregor Michalicek's avatar
Gregor Michalicek committed
125 126 127 128 129 130 131
   INTEGER, ALLOCATABLE                 :: parent(:)
   INTEGER, ALLOCATABLE                 :: pointer_EIBZ(:)
   INTEGER, ALLOCATABLE                 :: n_q(:)

   LOGICAL, ALLOCATABLE                 :: enabledConstraints(:)

#endif
132

133
#ifndef CPP_OLDINTEL
Gregor Michalicek's avatar
Gregor Michalicek committed
134 135 136

   WRITE(*,*) 'entered RDMFT subroutine'

137
   ! General initializations
Gregor Michalicek's avatar
Gregor Michalicek committed
138
   mixParam = 0.0001
Gregor Michalicek's avatar
Gregor Michalicek committed
139
   lagrangeMultiplier = 0.1 !results%ef
Gregor Michalicek's avatar
Gregor Michalicek committed
140
   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
141 142 143 144 145 146 147 148 149 150 151 152 153 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 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249

   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

   ! 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
   tempOcc = 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)
            IF(occupationVec(iState).LT.0.01) THEN
               addCharge = addCharge + (0.01-occupationVec(iState))*kpts%wtkpt(ikpt)
               addChargeWeight = addChargeWeight + kpts%wtkpt(ikpt)
            END IF
            IF(occupationVec(iState).GT.0.99) THEN
               subCharge = subCharge + (occupationVec(iState)-0.99)*kpts%wtkpt(ikpt)
               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
            IF(occupationVec(iState).LT.0.01) THEN
               occupationVec(iState) = occupationVec(iState) + 0.5*(subCharge+addCharge)*(kpts%wtkpt(ikpt)/addChargeWeight)
            END IF
            IF(occupationVec(iState).GT.0.99) THEN
               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)
            tempOcc = tempOcc + occupationVec(iState) * kpts%wtkpt(ikpt)
         END DO
      END DO
   END DO
   DEALLOCATE(occupationVec)

   ! Some more initializations

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

250 251
   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
252
   ALLOCATE(dEdOcc(MAXVAL(results%neig(1:kpts%nkpt,1:input%jspins)),kpts%nkpt,input%jspins))
253

254 255 256
   CALL regCharges%init(input,atoms)
   CALL dos%init(input,atoms,dimension,kpts,vacuum)
   CALL moments%init(input,atoms)
257 258
   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)
259
   IF (ALLOCATED(vTot%pw_w)) DEALLOCATE (vTot%pw_w)
Gregor Michalicek's avatar
Gregor Michalicek committed
260
   ALLOCATE(vTot%pw_w(SIZE(overallDen%pw,1),input%jspins))
261 262 263 264
   DO jspin = 1, input%jspins
      CALL convol(stars,vTot%pw_w(:,jspin),vTot%pw(:,jspin),stars%ufft)
   END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
265 266 267 268 269 270 271 272 273 274 275 276 277 278
   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

279 280
   vTotSSDen = 0.0

281
   ! Calculate all single state densities
282 283 284 285 286 287
   cdnvalJob%l_evp = .FALSE.
   cdnvalJob%nkptExtended = kpts%nkpt
   ALLOCATE(cdnvalJob%noccbd(kpts%nkpt))
   ALLOCATE(cdnvalJob%nStart(kpts%nkpt))
   ALLOCATE(cdnvalJob%nEnd(kpts%nkpt))
   ALLOCATE(cdnvalJob%weights(1,kpts%nkpt))
Gregor Michalicek's avatar
Gregor Michalicek committed
288
   numStates = 0
289
   DO jspin = 1, input%jspins
290
      jsp = MERGE(1,jspin,noco%l_noco)
291
      DO ikpt = 1, kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
292 293
         DO iBand = 1, highestState(ikpt,jsp)
            numStates = numStates + 1
294
            ! Construct cdnvalJob object for this state
295 296 297 298 299 300 301 302
            ! (Reasonable parallelization is not yet done - should be placed over the loops enclosing this section)
            cdnvalJob%ikptStart = ikpt
            cdnvalJob%ikptIncrement = kpts%nkpt
            IF(mpi%irank.EQ.0) THEN
               cdnvalJob%noccbd = 1
               cdnvalJob%nStart = iBand
               cdnvalJob%nEnd = iBand
               cdnvalJob%weights = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
303
               cdnvalJob%weights(1,ikpt) = spinDegenFac ! Note the doubling of the weight for non-spinpolarized calculations.
304 305 306 307 308 309
            ELSE
               cdnvalJob%noccbd = 0
               cdnvalJob%nStart = 1
               cdnvalJob%nEnd = 0
               cdnvalJob%weights = 0.0
            END IF
310 311

            ! Call cdnval to construct density
312 313
            WRITE(*,*) 'Note: some optional flags may have to be reset in rdmft before the cdnval call'
            WRITE(*,*) 'This is not yet implemented!'
314
            CALL singleStateDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
315 316
            CALL cdnval(eig_id,mpi,kpts,jsp,noco,input,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
                        sphhar,sym,vTot,oneD,cdnvalJob,singleStateDen,regCharges,dos,results,moments)
317 318

            ! Store the density on disc (These are probably way too many densities to keep them in memory)
319 320
            filename = ''
            WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
321 322 323 324
            IF (mpi%irank.EQ.0) THEN
               CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
                                 0,-1.0,0.0,.FALSE.,singleStateDen,TRIM(ADJUSTL(filename)))
            END IF
325 326 327 328
#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
329
            potDenInt = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
330
            CALL int_nv(jsp,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vTotTemp,singleStateDen,potDenInt)
331
            vTotSSDen(iBand,ikpt,jsp) = potDenInt
332 333 334 335
         END DO
      END DO
   END DO

Gregor Michalicek's avatar
Gregor Michalicek committed
336 337
   ! Initializations for exchange contributions

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

340 341 342 343 344
   IF(ALLOCATED(hybrid%ne_eig)) DEALLOCATE(hybrid%ne_eig)
   IF(ALLOCATED(hybrid%nbands)) DEALLOCATE(hybrid%nbands)
   IF(ALLOCATED(hybrid%nobd)) DEALLOCATE(hybrid%nobd)
   IF(ALLOCATED(hybrid%nbasm)) DEALLOCATE(hybrid%nbasm)
   IF(ALLOCATED(hybrid%div_vv)) DEALLOCATE(hybrid%div_vv)
Gregor Michalicek's avatar
Gregor Michalicek committed
345 346 347 348 349 350 351 352
   ALLOCATE(hybrid%ne_eig(kpts%nkpt),hybrid%nbands(kpts%nkpt),hybrid%nobd(kpts%nkptf))
   ALLOCATE(hybrid%nbasm(kpts%nkptf))
   ALLOCATE(hybrid%div_vv(DIMENSION%neigd,kpts%nkpt,input%jspins))
   l_restart = .FALSE.
   l_zref = (sym%zrfs.AND.(SUM(ABS(kpts%bk(3,:kpts%nkpt))).LT.1e-9).AND..NOT.noco%l_noco)
   iterHF = 0
   hybrid%l_calhf = .TRUE.

Gregor Michalicek's avatar
Gregor Michalicek committed
353
!   CALL open_hybrid_io1(DIMENSION,sym%invs)
Gregor Michalicek's avatar
Gregor Michalicek committed
354 355 356 357 358 359 360 361 362

   CALL mixedbasis(atoms,kpts,dimension,input,cell,sym,xcpot,hybrid,enpara,mpi,vTot,l_restart)

   CALL open_hybrid_io2(hybrid,DIMENSION,atoms,sym%invs)

   CALL coulombmatrix(mpi,atoms,kpts,cell,sym,hybrid,xcpot,l_restart)

   CALL hf_init(hybrid,kpts,atoms,input,DIMENSION,hybdat,sym%invs)

Gregor Michalicek's avatar
Gregor Michalicek committed
363
   WRITE(*,*) 'RDMFT: HF initializations end'
Gregor Michalicek's avatar
Gregor Michalicek committed
364 365 366 367 368 369 370

   ALLOCATE(parent(kpts%nkptf))
   ALLOCATE(exDiag(dimension%neigd,ikpt,input%jspins))
   ALLOCATE(lastGradient(numStates+1))
   ALLOCATE(lastParameters(numStates+1))
   lastGradient = 0.0
   lastParameters = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
371
   maxHistoryLength = 17!7
Gregor Michalicek's avatar
Gregor Michalicek committed
372 373 374 375 376 377 378
   ALLOCATE(gradientCorrections(numStates+1,maxHistoryLength))
   ALLOCATE(paramCorrections(numStates+1,maxHistoryLength))
   gradientCorrections = 0.0
   paramCorrections = 0.0
   istep = 0

   ! Occupation number optimization loop
379

380
   converged = .FALSE.
381 382
   DO WHILE (.NOT.converged)

Gregor Michalicek's avatar
Gregor Michalicek committed
383
      WRITE(*,*) 'RDMFT: convergence loop start'
Gregor Michalicek's avatar
Gregor Michalicek committed
384 385 386 387 388 389 390 391

      DO jspin = 1, input%jspins
         DO ikpt = 1,kpts%nkpt
            WRITE(*,*) 'jspin, ikpt: ', jspin, ikpt
            WRITE(*,'(8f10.5)') results%w_iks(1:10,ikpt,jspin)
         END DO
      END DO

392
      ! Calculate overall density with current occupation numbers (don't forget core electron density)
Gregor Michalicek's avatar
Gregor Michalicek committed
393
      CALL overallDen%resetPotDen()
394 395 396
      jspmax = input%jspins
      IF (noco%l_mperp) jspmax = 1
      DO jspin = 1,jspmax
397
         CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin,banddos=banddos)
398 399 400
         CALL cdnval(eig_id,mpi,kpts,jsp,noco,input,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
                     sphhar,sym,vTot,oneD,cdnvalJob,overallDen,regCharges,dos,results,moments)
      END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
401

402 403
      CALL cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
                   stars,cell,sphhar,atoms,vTot,overallDen,moments,results)
404
      IF (mpi%irank.EQ.0) THEN
405
         CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,overallDen,noco%l_noco,.TRUE.,.true.,fix)
406 407 408 409
      END IF
#ifdef CPP_MPI
      CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallDen)
#endif
410 411

      ! Calculate Coulomb potential for overall density (+including external potential)
412 413
      CALL overallDen%sum_both_spin()!workden)
      CALL overallVCoul%resetPotDen()
Gregor Michalicek's avatar
Gregor Michalicek committed
414 415 416 417
      ALLOCATE(overallVCoul%pw_w(size(overallVCoul%pw,1),size(overallVCoul%pw,2)))
      overallVCoul%pw_w(:,:) = 0.0
      CALL vgen_coulomb(1,mpi,DIMENSION,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,overallDen,overallVCoul)
      CALL convol(stars,overallVCoul%pw_w(:,1),overallVCoul%pw(:,1),stars%ufft)   ! Is there a problem with a second spin?!
418 419 420
#ifdef CPP_MPI
      CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallVCoul)
#endif
421

422
      overallVCoulSSDen = 0.0
423 424 425
      DO jspin = 1, input%jspins
         jsp = MERGE(1,jspin,noco%l_noco)
         DO ikpt = 1, kpts%nkpt
Gregor Michalicek's avatar
Gregor Michalicek committed
426
            DO iBand = 1, highestState(ikpt,jsp)
427 428 429 430 431 432 433 434 435 436 437
               ! 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
                  CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN_const,&
                                   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
438

439
               ! For each state calculate integral over Coulomb potential times single state density
440
               potDenInt = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
441
               CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vCoul,singleStateDen,potDenInt) ! Is there a problem with a second spin?!
442 443 444 445
               overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
            END DO
         END DO
      END DO
446

Gregor Michalicek's avatar
Gregor Michalicek committed
447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470
      ! 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)
         IF(ALLOCATED(hybdat%prod)) DEALLOCATE (hybdat%prod)
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE (hybdat%nindxp1)

Gregor Michalicek's avatar
Gregor Michalicek committed
471 472
         results%neig(:,:) = neigTemp(:,:)

Gregor Michalicek's avatar
Gregor Michalicek committed
473 474 475
         CALL HF_setup(hybrid,input,sym,kpts,dimension,atoms,mpi,noco,cell,oneD,results,jspin,enpara,eig_id,&
                       hybdat,iterHF,sym%invs,vTot%mt(:,0,:,:),eig_irr)

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

Gregor Michalicek's avatar
Gregor Michalicek committed
478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493
         mnobd = MAXVAL(hybrid%nobd)

         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)
            CALL symm_hf(kpts,ikpt,sym,dimension,hybdat,eig_irr,atoms,hybrid,cell,lapw,jspin,mpi,&
                         rrot,nsymop,psym,nkpt_EIBZ,n_q,parent,pointer_EIBZ,nsest,indx_sest)

            exMat%l_real=sym%invs
            CALL exchange_valence_hf(ikpt,kpts,nkpt_EIBZ, sym,atoms,hybrid,cell,dimension,input,jspin,hybdat,mnobd,lapw,&
                                     eig_irr,results,parent,pointer_EIBZ,n_q,wl_iks,iterHF,xcpot,noco,nsest,indx_sest,&
                                     mpi,exMat)
            CALL exchange_vccv1(ikpt,atoms,hybrid,hybdat,dimension,jspin,lapw,nsymop,nsest,indx_sest,mpi,1.0,results,exMat)
Gregor Michalicek's avatar
Gregor Michalicek committed
494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548

            !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

            CALL zMat%init(olap%l_real,nbasfcn,dimension%neigd)

            CALL read_eig(eig_id,ikpt,jspin,n_start=1,n_end=hybrid%nbands(ikpt),neig=nbands,zmat=zMat)

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

            CALL olap%multiply(zMat,trafo)

            CALL invtrafo%alloc(olap%l_real,hybrid%nbands(ikpt),nbasfcn)
            CALL trafo%TRANSPOSE(invtrafo)
            IF(.NOT.invtrafo%l_real) invtrafo%data_c = CONJG(invtrafo%data_c)

            DO i = 1, hybrid%nbands(ikpt)
               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)

            CALL symmetrizeh(atoms,kpts%bkf(:,ikpt),dimension,jspin,lapw,sym,hybdat%kveclo_eig,cell,nsymop,psym,exMatLAPW)

            IF (.NOT.exMatLAPW%l_real) exMatLAPW%data_c=conjg(exMatLAPW%data_c) 
            zMat%matsize1=MIN(zMat%matsize1,exMatLAPW%matsize2)

            CALL exMatLAPW%multiply(zMat,tmpMat)

Gregor Michalicek's avatar
Gregor Michalicek committed
549
            DO iBand = 1, highestState(ikpt,jsp)
Gregor Michalicek's avatar
Gregor Michalicek committed
550 551
               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
552
               ELSE
Gregor Michalicek's avatar
Gregor Michalicek committed
553
                  exDiag(iBand,ikpt,jspin) = REAL(dot_product(zMat%data_c(:zMat%matsize1,iband),tmpMat%data_c(:,iband)))
Gregor Michalicek's avatar
Gregor Michalicek committed
554 555
               END IF
            END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
556 557 558 559 560 561 562 563 564 565

            !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
566 567 568 569 570
         END DO
      END DO

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

Gregor Michalicek's avatar
Gregor Michalicek committed
571
      gradSum = 0.0
Gregor Michalicek's avatar
Gregor Michalicek committed
572 573 574 575 576 577
      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)
578
               occStateI = MAX(occStateI,minOcc)
Gregor Michalicek's avatar
Gregor Michalicek committed
579
!               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
580
!               occStateI = cdnvalJob%weights(iBand,ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
581 582 583 584 585 586
               rdmftFunctionalValue = 0.5*0.5*SQRT(1.0/occStateI) ! for Müller functional derivative

               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) + &
                                              overallVCoulSSDen(iBand,ikpt,isp) + exchangeTerm)
Gregor Michalicek's avatar
Gregor Michalicek committed
587

Gregor Michalicek's avatar
Gregor Michalicek committed
588 589 590 591 592 593 594 595
               WRITE(*,*) 'ENERGY GRADIENT CONTRIBUTIONS'
               WRITE(*,*) 'ispin, ikpt, iBand', ispin, ikpt, iBand
               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)
               WRITE(*,*) 'exchangeTerm', exchangeTerm
               WRITE(*,*) 'exDiag(iBand,ikpt,isp)', exDiag(iBand,ikpt,isp)
               WRITE(*,*) 'rdmftFunctionalValue', rdmftFunctionalValue
Gregor Michalicek's avatar
Gregor Michalicek committed
596

Gregor Michalicek's avatar
Gregor Michalicek committed
597 598

               gradSum = gradSum + dEdOcc(iBand,ikpt,isp) ! * results%w_iks(iBand,ikpt,isp)
Gregor Michalicek's avatar
Gregor Michalicek committed
599 600 601
            END DO
         END DO
      END DO
Gregor Michalicek's avatar
Gregor Michalicek committed
602
      lagrangeMultiplier = -gradSum / numStates !(input%zelec/(2.0/REAL(input%jspins)))
Gregor Michalicek's avatar
Gregor Michalicek committed
603

Gregor Michalicek's avatar
Gregor Michalicek committed
604
   WRITE(*,*) 'lagrangeMultiplier: ', lagrangeMultiplier
Gregor Michalicek's avatar
Gregor Michalicek committed
605

606 607
      ! Optimize occupation numbers

Gregor Michalicek's avatar
Gregor Michalicek committed
608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627
      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)
628
               occStateI = MAX(occStateI,minOcc)
Gregor Michalicek's avatar
Gregor Michalicek committed
629
               equalityLinCombi(iState) = kpts%wtkpt(ikpt)
Gregor Michalicek's avatar
Gregor Michalicek committed
630
               gradient(iState) = dEdOcc(iBand,ikpt,isp) + lagrangeMultiplier
Gregor Michalicek's avatar
Gregor Michalicek committed
631 632 633 634 635 636 637
               gradient(numStates+1) = gradient(numStates+1) + occStateI * kpts%wtkpt(ikpt)
               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
Gregor Michalicek's avatar
Gregor Michalicek committed
638
      parameters(numStates+1) = lagrangeMultiplier
Gregor Michalicek's avatar
Gregor Michalicek committed
639 640 641 642 643 644 645 646 647 648 649 650 651 652

      mixParam = 0.01 / MAXVAL(ABS(gradient(:numStates)))
      WRITE(*,*) 'mixParam: ', mixParam

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

      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
653
               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
654 655 656 657 658 659 660
            END DO
         END DO
      END DO

      DEALLOCATE (enabledConstraints,maxConstraints,minConstraints)
      DEALLOCATE (parameters,gradient,equalityLinCombi)

661 662 663

   END DO ! WHILE (.NOT.converged)

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

Gregor Michalicek's avatar
Gregor Michalicek committed
666
   hybrid%l_calhf = .FALSE.
Gregor Michalicek's avatar
Gregor Michalicek committed
667

668 669
   ! Calculate final overall density

Gregor Michalicek's avatar
Gregor Michalicek committed
670 671 672
   !I think we need most of cdngen at this place so I just use cdngen
   CALL outDen%resetPotDen()
   CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,DIMENSION,kpts,atoms,sphhar,stars,sym,&
Matthias Redies's avatar
merge  
Matthias Redies committed
673
               enpara,cell,noco,vTot,results,oneD,coreSpecInput,archiveType,xcpot,outDen, EnergyDen)
Gregor Michalicek's avatar
Gregor Michalicek committed
674 675 676 677 678 679 680 681 682

   ! 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)
Gregor Michalicek's avatar
Gregor Michalicek committed
683
            rdmftFunctionalValue = 0.5*SQRT(occStateI) ! for Müller functional
Gregor Michalicek's avatar
Gregor Michalicek committed
684

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

Gregor Michalicek's avatar
Gregor Michalicek committed
687 688 689
            rdmftEnergy = rdmftEnergy + exchangeTerm + &
                          occStateI * ((spinDegenFac*results%eig(iBand,ikpt,isp)) - vTotSSDen(iBand,ikpt,isp) + &
                                       overallVCoulSSDen(iBand,ikpt,isp))
Gregor Michalicek's avatar
Gregor Michalicek committed
690 691 692 693 694 695 696 697
         END DO
      END DO
   END DO

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

   WRITE(6,'(a,f20.10,a)') 'RDMFT energy: ', rdmftEnergy, ' Htr'

698
#endif
699 700 701
END SUBROUTINE rdmft

END MODULE m_rdmft