cdnval.F90 13.4 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 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.
!--------------------------------------------------------------------------------

7
MODULE m_cdnval
8 9 10

USE m_juDFT

11
CONTAINS
12

13
SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,stars,&
14
                  vacuum,dimension,sphhar,sym,vTot,oneD,cdnvalJob,den,regCharges,dos,results,&
15
                  moments,coreSpecInput,mcd,slab,orbcomp)
16 17 18 19 20 21 22 23 24 25 26

   !************************************************************************************
   !     This is the FLEUR valence density generator
   !******** ABBREVIATIONS *************************************************************
   !     noccbd   : number of occupied bands
   !     pallst   : if set to .true. bands above the Fermi-Energy are taken into account
   !     ener     : band energy averaged over all bands and k-points,
   !                wheighted with the l-like charge of each atom type
   !     sqal     : l-like charge of each atom type. sum over all k-points and bands
   !************************************************************************************

27
   USE m_types
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
   USE m_eig66_io
   USE m_genMTBasis
   USE m_calcDenCoeffs
   USE m_mcdinit
   USE m_sympsi
   USE m_eparas      ! energy parameters and partial charges
   USE m_qal21       ! off-diagonal part of partial charges
   USE m_abcof
   USE m_nmat        ! calculate density matrix for LDA + U
   USE m_vacden
   USE m_pwden
   USE m_forcea8
   USE m_checkdopall
   USE m_cdnmt       ! calculate the density and orbital moments etc.
   USE m_orbmom      ! coeffd for orbital moments
   USE m_qmtsl       ! These subroutines divide the input%film into vacuum%layers
   USE m_qintsl      ! (slabs) and intergate the DOS in these vacuum%layers
   USE m_orbcomp     ! calculate orbital composition (like p_x,p_y,p_z)
   USE m_abcrot2
   USE m_corespec, only : l_cs    ! calculation of core spectra (EELS)
   USE m_corespec_io, only : corespec_init
   USE m_corespec_eval, only : corespec_gaunt,corespec_rme,corespec_dos,corespec_ddscs
   USE m_xmlOutput
51
#ifdef CPP_MPI
52
   USE m_mpi_col_den ! collect density data from parallel nodes
53
#endif
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72

   IMPLICIT NONE

   TYPE(t_results),       INTENT(INOUT) :: results
   TYPE(t_mpi),           INTENT(IN)    :: mpi
   TYPE(t_dimension),     INTENT(IN)    :: dimension
   TYPE(t_oneD),          INTENT(IN)    :: oneD
   TYPE(t_enpara),        INTENT(IN)    :: enpara
   TYPE(t_banddos),       INTENT(IN)    :: banddos
   TYPE(t_input),         INTENT(IN)    :: input
   TYPE(t_vacuum),        INTENT(IN)    :: vacuum
   TYPE(t_noco),          INTENT(IN)    :: noco
   TYPE(t_sym),           INTENT(IN)    :: sym
   TYPE(t_stars),         INTENT(IN)    :: stars
   TYPE(t_cell),          INTENT(IN)    :: cell
   TYPE(t_kpts),          INTENT(IN)    :: kpts
   TYPE(t_sphhar),        INTENT(IN)    :: sphhar
   TYPE(t_atoms),         INTENT(IN)    :: atoms
   TYPE(t_potden),        INTENT(IN)    :: vTot
73
   TYPE(t_cdnvalJob),     INTENT(IN)    :: cdnvalJob
74 75
   TYPE(t_potden),        INTENT(INOUT) :: den
   TYPE(t_regionCharges), INTENT(INOUT) :: regCharges
76
   TYPE(t_dos),           INTENT(INOUT) :: dos
77
   TYPE(t_moments),       INTENT(INOUT) :: moments
78 79 80
   TYPE(t_coreSpecInput), OPTIONAL, INTENT(IN)    :: coreSpecInput
   TYPE(t_mcd),           OPTIONAL, INTENT(INOUT) :: mcd
   TYPE(t_slab),          OPTIONAL, INTENT(INOUT) :: slab
81
   TYPE(t_orbcomp),       OPTIONAL, INTENT(INOUT) :: orbcomp
82 83 84

   ! Scalar Arguments
   INTEGER,               INTENT(IN)    :: eig_id, jspin
85 86

#ifdef CPP_MPI
87
   INCLUDE 'mpif.h'
88
#endif
89

90 91
   ! Local Scalars
   INTEGER :: ikpt,jsp_start,jsp_end,ispin,jsp
92
   INTEGER :: iErr,nbands,noccbd,iType
93
   INTEGER :: skip_t,skip_tt,nStart,nEnd,nbasfcn
94
   LOGICAL :: l_orbcomprot, l_real, l_write, l_dosNdir, l_corespec
95 96 97 98 99 100 101 102 103 104 105 106 107

   ! Local Arrays
   REAL,    ALLOCATABLE :: we(:)
   REAL,    ALLOCATABLE :: eig(:)
   REAL,    ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions

   TYPE (t_lapw)             :: lapw
   TYPE (t_orb)              :: orb
   TYPE (t_denCoeffs)        :: denCoeffs
   TYPE (t_denCoeffsOffdiag) :: denCoeffsOffdiag
   TYPE (t_force)            :: force
   TYPE (t_eigVecCoeffs)     :: eigVecCoeffs
   TYPE (t_usdus)            :: usdus
108
   TYPE (t_mat)              :: zMat
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
   TYPE (t_gVacMap)          :: gVacMap

   CALL timestart("cdnval")

   l_real = sym%invs.AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco)
   l_dosNdir = banddos%dos.AND.(banddos%ndir.EQ.-3)

   IF (noco%l_mperp) THEN
      ! when the off-diag. part of the desinsity matrix, i.e. m_x and
      ! m_y, is calculated inside the muffin-tins (l_mperp = T), cdnval
      ! is called only once. therefore, several spin loops have been
      ! added. if l_mperp = F, these loops run only from jspin - jspin.
      jsp_start = 1
      jsp_end   = 2
   ELSE
      jsp_start = jspin
      jsp_end   = jspin
   END IF

   ALLOCATE (f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end)) ! Deallocation before mpi_col_den
   ALLOCATE (g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end))
   ALLOCATE (flo(atoms%jmtd,2,atoms%nlod,dimension%jspd))

   ! Initializations
   CALL usdus%init(atoms,input%jspins)
   CALL denCoeffs%init(atoms,sphhar,jsp_start,jsp_end)
   ! The last entry in denCoeffsOffdiag%init is l_fmpl. It is meant as a switch to a plot of the full magnet.
   ! density without the atomic sphere approximation for the magnet. density. It is not completely implemented (lo's missing).
   CALL denCoeffsOffdiag%init(atoms,noco,sphhar,.FALSE.)
   CALL force%init1(input,atoms)
   CALL orb%init(atoms,noco,jsp_start,jsp_end)

   IF (denCoeffsOffdiag%l_fmpl.AND.(.NOT.noco%l_mperp)) CALL juDFT_error("for fmpl set noco%l_mperp = T!" ,calledby ="cdnval")
   IF (l_dosNdir.AND.oneD%odi%d1) CALL juDFT_error("layer-resolved feature does not work with 1D",calledby ="cdnval")
143
   IF (banddos%l_mcd.AND..NOT.PRESENT(mcd)) CALL juDFT_error("mcd is missing",calledby ="cdnval")
144 145

   ! calculation of core spectra (EELS) initializations -start-
146 147 148 149 150 151 152
   l_coreSpec = .FALSE.
   IF (PRESENT(coreSpecInput)) THEN
      CALL corespec_init(input,atoms,coreSpecInput)
      IF(l_cs.AND.(mpi%isize.NE.1)) CALL juDFT_error('EELS + MPI not implemented', calledby = 'cdnval')
      IF(l_cs.AND.jspin.EQ.1) CALL corespec_gaunt()
      l_coreSpec = l_cs
   END IF
153 154 155 156 157 158 159
   ! calculation of core spectra (EELS) initializations -end-

   IF (mpi%irank==0) THEN
      WRITE (6,FMT=8000) jspin
      WRITE (16,FMT=8000) jspin
      CALL openXMLElementPoly('mtCharges',(/'spin'/),(/jspin/))
   END IF
160 161
8000 FORMAT (/,/,10x,'valence density: spin=',i2)

162 163 164 165 166 167 168
   l_write = input%cdinf.AND.mpi%irank==0
   DO iType = 1, atoms%ntype
      DO ispin = jsp_start, jsp_end
         CALL genMTBasis(atoms,enpara,vTot,mpi,iType,ispin,l_write,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin))
      END DO
      IF (noco%l_mperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
      IF (banddos%l_mcd) CALL mcd_init(atoms,input,dimension,vTot%mt(:,0,:,:),g,f,mcd,iType,jspin)
169 170
      IF (l_coreSpec) CALL corespec_rme(atoms,input,iType,dimension%nstd,input%jspins,jspin,results%ef,&
                                        dimension%msh,vTot%mt(:,0,:,:),f,g)
171 172
   END DO
   DEALLOCATE (f,g,flo)
173

174 175
   skip_tt = dot_product(enpara%skiplo(:atoms%ntype,jspin),atoms%neq(:atoms%ntype))
   IF (noco%l_soc.OR.noco%l_noco) skip_tt = 2 * skip_tt
176 177
   ALLOCATE (we(MAXVAL(cdnvalJob%noccbd(:))))
   ALLOCATE (eig(MAXVAL(cdnvalJob%noccbd(:))))
178
   jsp = MERGE(1,jspin,noco%l_noco)
179

180
   DO ikpt = cdnvalJob%ikptStart, cdnvalJob%nkptExtended, cdnvalJob%ikptIncrement
181

182
      IF (ikpt.GT.kpts%nkpt) THEN
183
#ifdef CPP_MPI
184
         CALL MPI_BARRIER(mpi%mpi_comm,iErr) ! Synchronizes the RMA operations
185
#endif
186 187 188 189 190
         EXIT
      END IF

      CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi)
      skip_t = skip_tt
191 192 193 194
      noccbd = cdnvalJob%noccbd(ikpt)
      nStart = cdnvalJob%nStart(ikpt)
      nEnd = cdnvalJob%nEnd(ikpt)
      we(1:noccbd) = cdnvalJob%weights(1:noccbd,ikpt)
195 196
      eig(1:noccbd) = results%eig(nStart:nEnd,ikpt,jsp)

197
      IF (cdnvalJob%l_evp) THEN
198 199 200 201 202 203 204 205
         IF (nStart > skip_tt) skip_t = 0
         IF (nEnd <= skip_tt) skip_t = noccbd
         IF ((nStart <= skip_tt).AND.(nEnd > skip_tt)) skip_t = mod(skip_tt,noccbd)
      END IF

      nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
      CALL zMat%init(l_real,nbasfcn,noccbd)
      CALL read_eig(eig_id,ikpt,jsp,n_start=nStart,n_end=nEnd,neig=nbands,zmat=zMat)
206
#ifdef CPP_MPI
207
      CALL MPI_BARRIER(mpi%mpi_comm,iErr) ! Synchronizes the RMA operations
208
#endif
209

210
      IF (noccbd.LE.0) CYCLE ! Note: This jump has to be after the MPI_BARRIER is called
211 212 213 214 215 216 217

      CALL gVacMap%init(dimension,sym,atoms,vacuum,stars,lapw,input,cell,kpts,enpara,vTot,ikpt,jspin)

      ! valence density in the interstitial and vacuum region has to be called only once (if jspin=1) in the non-collinear case
      IF (.NOT.((jspin.EQ.2).AND.noco%l_noco)) THEN
         ! valence density in the interstitial region
         CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,&
218
                    jspin,lapw,noccbd,we,eig,den,results,force%f_b8,zMat,dos)
219
         ! charge of each valence state in this k-point of the SBZ in the layer interstitial region of the film
220
         IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat)
221 222 223
         ! valence density in the vacuum region
         IF (input%film) THEN
            CALL vacden(vacuum,dimension,stars,oneD, kpts,input,sym,cell,atoms,noco,banddos,&
224
                        gVacMap,we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac,eig,den,zMat,dos)
225 226
         END IF
      END IF
227
      IF (input%film) CALL regCharges%sumBandsVac(vacuum,dos,noccbd,ikpt,jsp_start,jsp_end,eig,we)
228 229 230 231 232 233 234 235 236 237 238 239 240

      ! valence density in the atomic spheres
      CALL eigVecCoeffs%init(dimension,atoms,noco,jspin,noccbd)
      DO ispin = jsp_start, jsp_end
         IF (input%l_f) CALL force%init2(noccbd,input,atoms)
         CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,noco,ispin,oneD,&
                    eigVecCoeffs%acof(:,0:,:,ispin),eigVecCoeffs%bcof(:,0:,:,ispin),&
                    eigVecCoeffs%ccof(-atoms%llod:,:,:,:,ispin),zMat,eig,force)
         IF (atoms%n_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,den%mmpMat(:,:,:,jspin))

         ! perform Brillouin zone integration and summation over the
         ! bands in order to determine the energy parameters for each atom and angular momentum
         CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,&
241
                     skip_t,cdnvalJob%l_evp,eigVecCoeffs,usdus,regCharges,dos,banddos%l_mcd,mcd)
242

243
         IF (noco%l_mperp.AND.(ispin==jsp_end)) CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos)
244 245 246

         ! layer charge of each valence state in this k-point of the SBZ from the mt-sphere region of the film
         IF (l_dosNdir) THEN
247
            IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
248 249 250 251

            INQUIRE (file='orbcomprot',exist=l_orbcomprot)
            IF (l_orbcomprot) CALL abcrot2(atoms,noccbd,eigVecCoeffs,ispin) ! rotate ab-coeffs

252
            IF (PRESENT(orbcomp)) CALL orb_comp(ispin,ikpt,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
253 254 255 256 257 258 259
         END IF

         CALL calcDenCoeffs(atoms,sphhar,sym,we,noccbd,eigVecCoeffs,ispin,denCoeffs)

         IF (noco%l_soc) CALL orbmom(atoms,noccbd,we,ispin,eigVecCoeffs,orb)
         IF (input%l_f) CALL force%addContribsA21A12(input,atoms,dimension,sym,cell,oneD,enpara,&
                                                     usdus,eigVecCoeffs,noccbd,ispin,eig,we,results)
260 261
         IF(l_coreSpec) CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,dimension%neigd,&
                                          noccbd,results%ef,banddos%sig_dos,eig,we,eigVecCoeffs)
262 263 264
      END DO ! end loop over ispin
      IF (noco%l_mperp) CALL denCoeffsOffdiag%calcCoefficients(atoms,sphhar,sym,eigVecCoeffs,we,noccbd)

265
      IF ((banddos%dos.OR.banddos%vacdos.OR.input%cdinf).AND.(banddos%ndir.GT.0)) THEN
266
         ! since z is no longer an argument of cdninf sympsi has to be called here!
267
         CALL sympsi(lapw,jspin,sym,dimension,nbands,cell,eig,noco,dos%ksym(:,ikpt,jspin),dos%jsym(:,ikpt,jspin),zMat)
268 269
      END IF
   END DO ! end of k-point loop
270

271
#ifdef CPP_MPI
272
   DO ispin = jsp_start,jsp_end
273 274
      CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,input,noco,ispin,regCharges,dos,&
                       results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin),mcd,slab,orbcomp)
275
   END DO
276
#endif
277

278 279 280
   IF (mpi%irank==0) THEN
      CALL cdnmt(dimension%jspd,atoms,sphhar,noco,jsp_start,jsp_end,&
                 enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt)
281
      IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins)
282 283 284 285 286 287 288 289 290 291
      DO ispin = jsp_start,jsp_end
         IF (input%cdinf) THEN
            WRITE (6,FMT=8210) ispin
8210        FORMAT (/,5x,'check continuity of cdn for spin=',i2)
            CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,den,ispin)
         END IF
         IF (input%l_f) CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,force,results)
      END DO
      CALL closeXMLElement('mtCharges')
   END IF
292

293 294 295
   CALL timestop("cdnval")

END SUBROUTINE cdnval
296

297
END MODULE m_cdnval