rdmft.F90 7.94 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 44 45 46 47 48 49
   TYPE(t_potden),        INTENT(IN)    :: vTot
   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 66 67
   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))

68 69 70
   converged = .FALSE.

   ! Calculate all single state densities
71 72 73
   CALL regCharges%init(input,atoms)
   CALL dos%init(input,atoms,dimension,kpts,vacuum)
   CALL moments%init(input,atoms)
74 75
   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)
76 77 78 79 80 81
   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))
82
   DO jspin = 1, input%jspins
83
      jsp = MERGE(1,jspin,noco%l_noco)
84
      DO ikpt = 1, kpts%nkpt
85
         DO iBand = 1, results%neig(ikpt,jsp)
86
            ! Construct cdnvalJob object for this state
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
            ! (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
102 103

            ! Call cdnval to construct density
104 105 106 107 108
            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)
109 110

            ! Store the density on disc (These are probably way too many densities to keep them in memory)
111 112
            filename = ''
            WRITE(filename,'(a,i1.1,a,i4.4,a,i5.5)') 'cdn-', jsp, '-', ikpt, '-', iBand
113 114 115 116
            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
117 118 119 120 121 122 123
         END DO
      END DO
   END DO

   DO WHILE (.NOT.converged)

      ! Calculate overall density with current occupation numbers (don't forget core electron density)
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
      jspmax = input%jspins
      IF (noco%l_mperp) jspmax = 1
      DO jspin = 1,jspmax
         CALL cdnvalJob%init(mpi,input,kpts,banddos,noco,results,jspin)
         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
      CALL cdncore(results,mpi,dimension,oneD,input,vacuum,noco,sym,&
                   stars,cell,sphhar,atoms,vTot,overallDen,moments)
      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
139 140

      ! Calculate Coulomb potential for overall density (+including external potential)
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
      CALL overallDen%sum_both_spin()!workden)
      CALL overallVCoul%resetPotDen()
      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?!

      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
161

162 163 164 165 166
               ! For all states calculate integral over Coulomb potential times single state density
               potDenInt = 0.0
               CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,&              ! Is there a problem with a second spin?!
                           overallVCoul,singleStateDen,potDenInt)
               overallVCoulSSDen(iBand,ikpt,jsp) = potDenInt
167

168 169 170 171 172
               ! For all states calculate Integral over other potential contributions times single state density

            END DO
         END DO
      END DO
173 174 175 176 177 178 179 180 181 182 183 184 185

      ! Construct exchange matrix in the basis of eigenstates

      ! Optimize occupation numbers

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

   END DO ! WHILE (.NOT.converged)

   ! Calculate final overall density

   ! Calculate total energy

186 187 188
END SUBROUTINE rdmft

END MODULE m_rdmft