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 15
                  vacuum,dimension,sphhar,sym,obsolete,vTot,oneD,cdnvalJob,den,regCharges,dos,results,&
                  moments,orbcomp,coreSpecInput,mcd,slab)
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

   !************************************************************************************
   !     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
   !************************************************************************************

   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_types
   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 73

   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_obsolete),      INTENT(IN)    :: obsolete
   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
74
   TYPE(t_cdnvalJob),     INTENT(IN)    :: cdnvalJob
75 76
   TYPE(t_potden),        INTENT(INOUT) :: den
   TYPE(t_regionCharges), INTENT(INOUT) :: regCharges
77
   TYPE(t_dos),           INTENT(INOUT) :: dos
78
   TYPE(t_moments),       INTENT(INOUT) :: moments
79
   TYPE(t_orbcomp),       INTENT(INOUT) :: orbcomp
80 81 82
   TYPE(t_coreSpecInput), OPTIONAL, INTENT(IN)    :: coreSpecInput
   TYPE(t_mcd),           OPTIONAL, INTENT(INOUT) :: mcd
   TYPE(t_slab),          OPTIONAL, INTENT(INOUT) :: slab
83 84 85

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

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

91 92
   ! Local Scalars
   INTEGER :: ikpt,jsp_start,jsp_end,ispin,jsp
93
   INTEGER :: iErr,nbands,noccbd,iType
94
   INTEGER :: skip_t,skip_tt,nStart,nEnd,nbasfcn
95
   LOGICAL :: l_orbcomprot, l_real, l_write, l_dosNdir, l_corespec
96 97 98 99 100 101 102 103 104 105 106 107 108 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 143

   ! 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
   TYPE (t_zMat)             :: zMat
   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")
144
   IF (banddos%l_mcd.AND..NOT.PRESENT(mcd)) CALL juDFT_error("mcd is missing",calledby ="cdnval")
145 146

   ! calculation of core spectra (EELS) initializations -start-
147 148 149 150 151 152 153
   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
154 155 156 157 158 159 160
   ! 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
161 162
8000 FORMAT (/,/,10x,'valence density: spin=',i2)

163 164 165 166 167 168 169
   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)
170 171
      IF (l_coreSpec) CALL corespec_rme(atoms,input,iType,dimension%nstd,input%jspins,jspin,results%ef,&
                                        dimension%msh,vTot%mt(:,0,:,:),f,g)
172 173
   END DO
   DEALLOCATE (f,g,flo)
174

175 176
   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
177 178
   ALLOCATE (we(MAXVAL(cdnvalJob%noccbd(:))))
   ALLOCATE (eig(MAXVAL(cdnvalJob%noccbd(:))))
179
   jsp = MERGE(1,jspin,noco%l_noco)
180

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

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

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

198
      IF (cdnvalJob%l_evp) THEN
199 200 201 202 203 204 205 206
         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)
207
#ifdef CPP_MPI
208
      CALL MPI_BARRIER(mpi%mpi_comm,iErr) ! Synchronizes the RMA operations
209
#endif
210

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

      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,&
219
                    jspin,lapw,noccbd,we,eig,den,results,force%f_b8,zMat,dos)
220
         ! charge of each valence state in this k-point of the SBZ in the layer interstitial region of the film
221
         IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat)
222 223 224
         ! valence density in the vacuum region
         IF (input%film) THEN
            CALL vacden(vacuum,dimension,stars,oneD, kpts,input,sym,cell,atoms,noco,banddos,&
225
                        gVacMap,we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac0,eig,den,zMat,dos)
226 227
         END IF
      END IF
228
      IF (input%film) CALL regCharges%sumBandsVac(vacuum,dos,noccbd,ikpt,jsp_start,jsp_end,eig,we)
229 230 231 232 233 234 235 236 237 238 239 240 241

      ! 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,&
242
                     skip_t,cdnvalJob%l_evp,eigVecCoeffs,usdus,regCharges,dos,banddos%l_mcd,mcd)
243

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

         ! 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
248
            IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
249 250 251 252

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

253
            CALL orb_comp(ispin,ikpt,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
254 255 256 257 258 259 260
         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)
261 262
         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)
263 264 265
      END DO ! end loop over ispin
      IF (noco%l_mperp) CALL denCoeffsOffdiag%calcCoefficients(atoms,sphhar,sym,eigVecCoeffs,we,noccbd)

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

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

279 280 281
   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)
282
      IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins)
283 284 285 286 287 288 289 290 291 292
      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
293

294 295 296
   CALL timestop("cdnval")

END SUBROUTINE cdnval
297

298
END MODULE m_cdnval