rdmft.F90 8.73 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
!--------------------------------------------------------------------------------
! 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

SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
12
                 sphhar,sym,field,vTot,oneD,noco,results)
13 14 15

   USE m_types
   USE m_juDFT
16
   USE m_constants
17
#ifndef CPP_OLDINTEL
18 19
   USE m_cdnval
   USE m_cdn_io
20 21
   USE m_cdncore
   USE m_qfix
22 23 24
   USE m_vgen_coulomb
   USE m_convol
   USE m_intnv
25 26
#ifdef CPP_MPI
   USE m_mpi_bc_potden
27
#endif
28
#endif
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

   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
   TYPE(t_cell),          INTENT(IN)    :: cell
   TYPE(t_atoms),         INTENT(IN)    :: atoms
   TYPE(t_enpara),        INTENT(IN)    :: enpara
   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
44
   TYPE(t_field),         INTENT(INOUT) :: field
45
   TYPE(t_potden),        INTENT(INOUT) :: vTot
46 47 48 49 50 51
   TYPE(t_oneD),          INTENT(IN)    :: oneD
   TYPE(t_noco),          INTENT(IN)    :: noco
   TYPE(t_results),       INTENT(INOUT) :: results

   INTEGER,               INTENT(IN)    :: eig_id

52
#ifndef CPP_OLDINTEL
53
   TYPE(t_cdnvalJob)                    :: cdnvalJob
54
   TYPE(t_potden)                       :: singleStateDen, overallDen, overallVCoul
55 56 57
   TYPE(t_regionCharges)                :: regCharges
   TYPE(t_dos)                          :: dos
   TYPE(t_moments)                      :: moments
58
   INTEGER                              :: jspin, ikpt, iBand, jsp, jspmax
59 60
   REAL                                 :: fix, potDenInt, fermiEnergyTemp
   LOGICAL                              :: converged, l_qfix
61
   CHARACTER(LEN=20)                    :: filename
62

63 64
   REAL, ALLOCATABLE                    :: overallVCoulSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vTotSSDen(:,:,:)
65
#endif
66

67 68
   CALL juDFT_error('rdmft not yet implemented!', calledby = 'rdmft')

69
#ifndef CPP_OLDINTEL
70
   ! General initializations
71 72 73
   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))

74 75 76
   CALL regCharges%init(input,atoms)
   CALL dos%init(input,atoms,dimension,kpts,vacuum)
   CALL moments%init(input,atoms)
77 78
   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)
79 80 81 82 83 84 85
   ALLOCATE(overallVCoul%pw_w(SIZE(overallDen%pw,1),1))
   IF (ALLOCATED(vTot%pw_w)) DEALLOCATE (vTot%pw_w)
   ALLOCATE(vTot%pw_w(SIZE(overallDen%pw,1),1))
   DO jspin = 1, input%jspins
      CALL convol(stars,vTot%pw_w(:,jspin),vTot%pw(:,jspin),stars%ufft)
   END DO

86 87
   vTotSSDen = 0.0

88
   ! Calculate all single state densities
89 90 91 92 93 94
   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))
95
   DO jspin = 1, input%jspins
96
      jsp = MERGE(1,jspin,noco%l_noco)
97
      DO ikpt = 1, kpts%nkpt
98
         DO iBand = 1, results%neig(ikpt,jsp)
99
            ! Construct cdnvalJob object for this state
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
            ! (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
               cdnvalJob%weights(1,ikpt) = 1.0
            ELSE
               cdnvalJob%noccbd = 0
               cdnvalJob%nStart = 1
               cdnvalJob%nEnd = 0
               cdnvalJob%weights = 0.0
            END IF
115 116

            ! Call cdnval to construct density
117 118
            WRITE(*,*) 'Note: some optional flags may have to be reset in rdmft before the cdnval call'
            WRITE(*,*) 'This is not yet implemented!'
119
            CALL singleStateDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
120 121
            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)
122 123

            ! Store the density on disc (These are probably way too many densities to keep them in memory)
124 125
            filename = ''
            WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
126 127 128 129
            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
130 131 132
#ifdef CPP_MPI
            CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,singleStateDen)
#endif
133

134
            ! For each state calculate Integral over KS effective potential times single state density
135 136 137
            potDenInt = 0.0
            CALL int_nv(jsp,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,vTot,singleStateDen,potDenInt)
            vTotSSDen(iBand,ikpt,jsp) = potDenInt
138 139 140 141
         END DO
      END DO
   END DO

142 143 144
   ! Construct exchange matrix in the basis of eigenstates
   ! TODO!!!!!

145
   converged = .FALSE.
146 147 148
   DO WHILE (.NOT.converged)

      ! Calculate overall density with current occupation numbers (don't forget core electron density)
149 150 151
      jspmax = input%jspins
      IF (noco%l_mperp) jspmax = 1
      DO jspin = 1,jspmax
152
         CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin,banddos=banddos)
153 154 155
         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
156 157
      CALL cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
                   stars,cell,sphhar,atoms,vTot,overallDen,moments,results)
158
      IF (mpi%irank.EQ.0) THEN
159
         CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,overallDen,noco%l_noco,.TRUE.,.true.,fix)
160 161 162 163
      END IF
#ifdef CPP_MPI
      CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallDen)
#endif
164 165

      ! Calculate Coulomb potential for overall density (+including external potential)
166 167
      CALL overallDen%sum_both_spin()!workden)
      CALL overallVCoul%resetPotDen()
168
      CALL vgen_coulomb( 1, mpi, DIMENSION, oneD, input, field, vacuum, sym, stars, cell, &
169
           sphhar, atoms, overallDen, overallVCoul )
170
      CALL convol(stars, overallVCoul%pw_w(:,1), overallVCoul%pw(:,1), stars%ufft)   ! Is there a problem with a second spin?!
171 172 173
#ifdef CPP_MPI
      CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,overallVCoul)
#endif
174

175
      overallVCoulSSDen = 0.0
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
      DO jspin = 1, input%jspins
         jsp = MERGE(1,jspin,noco%l_noco)
         DO ikpt = 1, kpts%nkpt
            DO iBand = 1, results%neig(ikpt,jsp)
               ! 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
191

192
               ! For each state calculate integral over Coulomb potential times single state density
193
               potDenInt = 0.0
194
               CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,overallVCoul,singleStateDen,potDenInt) ! Is there a problem with a second spin?!
195 196 197 198
               overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
            END DO
         END DO
      END DO
199 200 201 202 203 204 205 206 207 208

      ! Optimize occupation numbers

      ! Check convergence of occupation numbers and set "converged" flag

   END DO ! WHILE (.NOT.converged)

   ! Calculate final overall density

   ! Calculate total energy
209
#endif
210 211 212
END SUBROUTINE rdmft

END MODULE m_rdmft