rdmft.F90 8.67 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
18
   USE m_cdnval
   USE m_cdn_io
19
20
   USE m_cdncore
   USE m_qfix
21
22
23
   USE m_vgen_coulomb
   USE m_convol
   USE m_intnv
24
25
26
#ifdef CPP_MPI
   USE m_mpi_bc_potden
#endif
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

   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
42
   TYPE(t_field),         INTENT(INOUT) :: field
43
   TYPE(t_potden),        INTENT(INOUT) :: vTot
44
45
46
47
48
49
   TYPE(t_oneD),          INTENT(IN)    :: oneD
   TYPE(t_noco),          INTENT(IN)    :: noco
   TYPE(t_results),       INTENT(INOUT) :: results

   INTEGER,               INTENT(IN)    :: eig_id

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

60
61
62
   REAL, ALLOCATABLE                    :: overallVCoulSSDen(:,:,:)
   REAL, ALLOCATABLE                    :: vTotSSDen(:,:,:)

63
64
   CALL juDFT_error('rdmft not yet implemented!', calledby = 'rdmft')

65
   ! General initializations
66
67
68
   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))

69
70
71
   CALL regCharges%init(input,atoms)
   CALL dos%init(input,atoms,dimension,kpts,vacuum)
   CALL moments%init(input,atoms)
72
73
   CALL overallDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
   CALL overallVCoul%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,POTDEN_TYPE_POTCOUL)
74
75
76
77
78
79
80
   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

81
82
   vTotSSDen = 0.0

83
   ! Calculate all single state densities
84
85
86
87
88
89
   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))
90
   DO jspin = 1, input%jspins
91
      jsp = MERGE(1,jspin,noco%l_noco)
92
      DO ikpt = 1, kpts%nkpt
93
         DO iBand = 1, results%neig(ikpt,jsp)
94
            ! Construct cdnvalJob object for this state
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
            ! (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
110
111

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

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

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

137
138
139
   ! Construct exchange matrix in the basis of eigenstates
   ! TODO!!!!!

140
   converged = .FALSE.
141
142
143
   DO WHILE (.NOT.converged)

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

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

170
      overallVCoulSSDen = 0.0
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
      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
186

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

      ! Optimize occupation numbers

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

   END DO ! WHILE (.NOT.converged)

   ! Calculate final overall density

   ! Calculate total energy

205
206
207
END SUBROUTINE rdmft

END MODULE m_rdmft