LCOV - code coverage report
Current view: top level - cdn - cdnval.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 101 110 91.8 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_cdnval
       8             : 
       9             : USE m_juDFT
      10             : 
      11             : CONTAINS
      12             : 
      13         608 : 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             :    !     This is the FLEUR valence density generator
      19             :    !******** ABBREVIATIONS *************************************************************
      20             :    !     noccbd   : number of occupied bands
      21             :    !     pallst   : if set to .true. bands above the Fermi-Energy are taken into account
      22             :    !     ener     : band energy averaged over all bands and k-points,
      23             :    !                wheighted with the l-like charge of each atom type
      24             :    !     sqal     : l-like charge of each atom type. sum over all k-points and bands
      25             :    !************************************************************************************
      26             : 
      27             :    USE m_types
      28             :    USE m_eig66_io
      29             :    USE m_genMTBasis
      30             :    USE m_calcDenCoeffs
      31             :    USE m_mcdinit
      32             :    USE m_sympsi
      33             :    USE m_eparas      ! energy parameters and partial charges
      34             :    USE m_qal21       ! off-diagonal part of partial charges
      35             :    USE m_abcof
      36             :    USE m_nmat        ! calculate density matrix for LDA + U
      37             :    USE m_vacden
      38             :    USE m_pwden
      39             :    USE m_forcea8
      40             :    USE m_checkdopall
      41             :    USE m_cdnmt       ! calculate the density and orbital moments etc.
      42             :    USE m_orbmom      ! coeffd for orbital moments
      43             :    USE m_qmtsl       ! These subroutines divide the input%film into vacuum%layers
      44             :    USE m_qintsl      ! (slabs) and intergate the DOS in these vacuum%layers
      45             :    USE m_orbcomp     ! calculate orbital composition (like p_x,p_y,p_z)
      46             :    USE m_abcrot2
      47             :    USE m_corespec, only : l_cs    ! calculation of core spectra (EELS)
      48             :    USE m_corespec_io, only : corespec_init
      49             :    USE m_corespec_eval, only : corespec_gaunt,corespec_rme,corespec_dos,corespec_ddscs
      50             :    USE m_xmlOutput
      51             : #ifdef CPP_MPI
      52             :    USE m_mpi_col_den ! collect density data from parallel nodes
      53             : #endif
      54             : 
      55             :    IMPLICIT NONE
      56             : 
      57             :    TYPE(t_results),       INTENT(INOUT) :: results
      58             :    TYPE(t_mpi),           INTENT(IN)    :: mpi
      59             :    TYPE(t_dimension),     INTENT(IN)    :: dimension
      60             :    TYPE(t_oneD),          INTENT(IN)    :: oneD
      61             :    TYPE(t_enpara),        INTENT(IN)    :: enpara
      62             :    TYPE(t_banddos),       INTENT(IN)    :: banddos
      63             :    TYPE(t_input),         INTENT(IN)    :: input
      64             :    TYPE(t_vacuum),        INTENT(IN)    :: vacuum
      65             :    TYPE(t_noco),          INTENT(IN)    :: noco
      66             :    TYPE(t_sym),           INTENT(IN)    :: sym
      67             :    TYPE(t_stars),         INTENT(IN)    :: stars
      68             :    TYPE(t_cell),          INTENT(IN)    :: cell
      69             :    TYPE(t_kpts),          INTENT(IN)    :: kpts
      70             :    TYPE(t_sphhar),        INTENT(IN)    :: sphhar
      71             :    TYPE(t_atoms),         INTENT(IN)    :: atoms
      72             :    TYPE(t_potden),        INTENT(IN)    :: vTot
      73             :    TYPE(t_cdnvalJob),     INTENT(IN)    :: cdnvalJob
      74             :    TYPE(t_potden),        INTENT(INOUT) :: den
      75             :    TYPE(t_regionCharges), INTENT(INOUT) :: regCharges
      76             :    TYPE(t_dos),           INTENT(INOUT) :: dos
      77             :    TYPE(t_moments),       INTENT(INOUT) :: moments
      78             :    TYPE(t_coreSpecInput), OPTIONAL, INTENT(IN)    :: coreSpecInput
      79             :    TYPE(t_mcd),           OPTIONAL, INTENT(INOUT) :: mcd
      80             :    TYPE(t_slab),          OPTIONAL, INTENT(INOUT) :: slab
      81             :    TYPE(t_orbcomp),       OPTIONAL, INTENT(INOUT) :: orbcomp
      82             : 
      83             :    ! Scalar Arguments
      84             :    INTEGER,               INTENT(IN)    :: eig_id, jspin
      85             : 
      86             : #ifdef CPP_MPI
      87             :    INCLUDE 'mpif.h'
      88             : #endif
      89             : 
      90             :    ! Local Scalars
      91             :    INTEGER :: ikpt,ikpt_i,jsp_start,jsp_end,ispin,jsp
      92             :    INTEGER :: iErr,nbands,noccbd,iType
      93             :    INTEGER :: skip_t,skip_tt,nbasfcn
      94             :    LOGICAL :: l_orbcomprot, l_real, l_dosNdir, l_corespec
      95             : 
      96             :    ! Local Arrays
      97        1216 :    REAL,ALLOCATABLE :: we(:),eig(:)
      98         608 :    INTEGER,ALLOCATABLE :: ev_list(:)
      99        1216 :    REAL,    ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions
     100             : 
     101         608 :    TYPE (t_lapw)             :: lapw
     102         608 :    TYPE (t_orb)              :: orb
     103         608 :    TYPE (t_denCoeffs)        :: denCoeffs
     104         608 :    TYPE (t_denCoeffsOffdiag) :: denCoeffsOffdiag
     105         608 :    TYPE (t_force)            :: force
     106         608 :    TYPE (t_eigVecCoeffs)     :: eigVecCoeffs
     107         608 :    TYPE (t_usdus)            :: usdus
     108         608 :    TYPE (t_mat)              :: zMat
     109         608 :    TYPE (t_gVacMap)          :: gVacMap
     110             : 
     111         608 :    CALL timestart("cdnval")
     112             : 
     113         608 :    l_real = sym%invs.AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco)
     114         608 :    l_dosNdir = banddos%dos.AND.(banddos%ndir.EQ.-3)
     115             : 
     116         608 :    IF (noco%l_mperp) THEN
     117             :       ! when the off-diag. part of the desinsity matrix, i.e. m_x and
     118             :       ! m_y, is calculated inside the muffin-tins (l_mperp = T), cdnval
     119             :       ! is called only once. therefore, several spin loops have been
     120             :       ! added. if l_mperp = F, these loops run only from jspin - jspin.
     121           0 :       jsp_start = 1
     122           0 :       jsp_end   = 2
     123             :    ELSE
     124         608 :       jsp_start = jspin
     125         608 :       jsp_end   = jspin
     126             :    END IF
     127             : 
     128         608 :    ALLOCATE (f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end)) ! Deallocation before mpi_col_den
     129         608 :    ALLOCATE (g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end))
     130         608 :    ALLOCATE (flo(atoms%jmtd,2,atoms%nlod,input%jspins))
     131             : 
     132             :    ! Initializations
     133         608 :    CALL usdus%init(atoms,input%jspins)
     134         608 :    CALL denCoeffs%init(atoms,sphhar,jsp_start,jsp_end)
     135             :    ! The last entry in denCoeffsOffdiag%init is l_fmpl. It is meant as a switch to a plot of the full magnet.
     136             :    ! density without the atomic sphere approximation for the magnet. density. It is not completely implemented (lo's missing).
     137         608 :    CALL denCoeffsOffdiag%init(atoms,noco,sphhar,noco%l_mtnocopot)
     138         608 :    CALL force%init1(input,atoms)
     139         608 :    CALL orb%init(atoms,noco,jsp_start,jsp_end)
     140             : 
     141         608 :    IF (denCoeffsOffdiag%l_fmpl.AND.(.NOT.noco%l_mperp)) CALL juDFT_error("for fmpl set noco%l_mperp = T!" ,calledby ="cdnval")
     142         608 :    IF (l_dosNdir.AND.oneD%odi%d1) CALL juDFT_error("layer-resolved feature does not work with 1D",calledby ="cdnval")
     143         608 :    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         608 :    l_coreSpec = .FALSE.
     147         608 :    IF (PRESENT(coreSpecInput)) THEN
     148         608 :       CALL corespec_init(input,atoms,coreSpecInput)
     149         608 :       IF(l_cs.AND.(mpi%isize.NE.1)) CALL juDFT_error('EELS + MPI not implemented', calledby = 'cdnval')
     150         608 :       IF(l_cs.AND.jspin.EQ.1) CALL corespec_gaunt()
     151         608 :       l_coreSpec = l_cs
     152             :    END IF
     153             :    ! calculation of core spectra (EELS) initializations -end-
     154             : 
     155         608 :    IF (mpi%irank==0) THEN
     156         304 :       WRITE (6,FMT=8000) jspin
     157         304 :       CALL openXMLElementPoly('mtCharges',(/'spin'/),(/jspin/))
     158             :    END IF
     159             : 8000 FORMAT (/,/,10x,'valence density: spin=',i2)
     160             : 
     161        1772 :    DO iType = 1, atoms%ntype
     162        2328 :       DO ispin = jsp_start, jsp_end
     163        2328 :          CALL genMTBasis(atoms,enpara,vTot,mpi,iType,ispin,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin))
     164             :       END DO
     165        1164 :       IF (noco%l_mperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
     166        1164 :       IF (banddos%l_mcd) CALL mcd_init(atoms,input,dimension,vTot%mt(:,0,:,:),g,f,mcd,iType,jspin)
     167        1164 :       IF (l_coreSpec) CALL corespec_rme(atoms,input,iType,dimension%nstd,input%jspins,jspin,results%ef,&
     168         608 :                                         dimension%msh,vTot%mt(:,0,:,:),f,g)
     169             :    END DO
     170         608 :    DEALLOCATE (f,g,flo)
     171             : 
     172         608 :    skip_tt = dot_product(enpara%skiplo(:atoms%ntype,jspin),atoms%neq(:atoms%ntype))
     173         608 :    IF (noco%l_soc.OR.noco%l_noco) skip_tt = 2 * skip_tt
     174             : 
     175         608 :    jsp = MERGE(1,jspin,noco%l_noco)
     176             : 
     177        2460 :    DO ikpt_i = 1,size(cdnvalJob%k_list)
     178        1852 :       ikpt=cdnvalJob%k_list(ikpt_i)
     179             : 
     180        1852 :       CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi)
     181        1852 :       skip_t = skip_tt
     182        1852 :       ev_list=cdnvaljob%compact_ev_list(ikpt_i,banddos%dos)
     183        1852 :       noccbd = SIZE(ev_list)
     184        1852 :       we  = cdnvalJob%weights(ev_list,ikpt)
     185        1852 :       eig = results%eig(ev_list,ikpt,jsp)
     186             : 
     187        1852 :       IF (cdnvalJob%l_evp) THEN
     188        1548 :          IF (minval(ev_list) > skip_tt) skip_t = 0
     189       19426 :          IF (maxval(ev_list) <= skip_tt) skip_t = noccbd
     190       19426 :          IF ((minval(ev_list) <= skip_tt).AND.(maxval(ev_list) > skip_tt)) skip_t = mod(skip_tt,noccbd)
     191             :       END IF
     192             : 
     193        1852 :       nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
     194        1852 :       CALL zMat%init(l_real,nbasfcn,noccbd)
     195        1852 :       CALL read_eig(eig_id,ikpt,jsp,list=ev_list,neig=nbands,zmat=zMat)
     196             : #ifdef CPP_MPI
     197        1852 :       CALL MPI_BARRIER(mpi%mpi_comm,iErr) ! Synchronizes the RMA operations
     198             : #endif
     199             : 
     200        1852 :       IF (noccbd.LE.0) CYCLE ! Note: This jump has to be after the MPI_BARRIER is called
     201             : 
     202        1848 :       CALL gVacMap%init(dimension,sym,atoms,vacuum,stars,lapw,input,cell,kpts,enpara,vTot,ikpt,jspin)
     203             : 
     204             :       ! valence density in the interstitial and vacuum region has to be called only once (if jspin=1) in the non-collinear case
     205        1848 :       IF (.NOT.((jspin.EQ.2).AND.noco%l_noco)) THEN
     206             :          ! valence density in the interstitial region
     207             :          CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,&
     208        1184 :                     jspin,lapw,noccbd,we,eig,den,results,force%f_b8,zMat,dos)
     209             :          ! charge of each valence state in this k-point of the SBZ in the layer interstitial region of the film
     210        1184 :          IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat)
     211             :          ! valence density in the vacuum region
     212        1184 :          IF (input%film) THEN
     213             :             CALL vacden(vacuum,dimension,stars,oneD, kpts,input,sym,cell,atoms,noco,banddos,&
     214          28 :                         gVacMap,we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac,eig,den,zMat,dos)
     215             :          END IF
     216             :       END IF
     217        1848 :       IF (input%film) CALL regCharges%sumBandsVac(vacuum,dos,noccbd,ikpt,jsp_start,jsp_end,eig,we)
     218             : 
     219             :       ! valence density in the atomic spheres
     220        1848 :       CALL eigVecCoeffs%init(input,DIMENSION,atoms,noco,jspin,noccbd)
     221        3696 :       DO ispin = jsp_start, jsp_end
     222        1848 :          IF (input%l_f) CALL force%init2(noccbd,input,atoms)
     223             :          CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,noco,ispin,oneD,&
     224             :                     eigVecCoeffs%acof(:,0:,:,ispin),eigVecCoeffs%bcof(:,0:,:,ispin),&
     225        1848 :                     eigVecCoeffs%ccof(-atoms%llod:,:,:,:,ispin),zMat,eig,force)
     226        1848 :          IF (atoms%n_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,den%mmpMat(:,:,:,jspin))
     227             : 
     228             :          ! perform Brillouin zone integration and summation over the
     229             :          ! bands in order to determine the energy parameters for each atom and angular momentum
     230             :          CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,&
     231        1848 :                      skip_t,cdnvalJob%l_evp,eigVecCoeffs,usdus,regCharges,dos,banddos%l_mcd,mcd)
     232             : 
     233        1848 :          IF (noco%l_mperp.AND.(ispin==jsp_end)) CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos)
     234             : 
     235             :          ! layer charge of each valence state in this k-point of the SBZ from the mt-sphere region of the film
     236        1848 :          IF (l_dosNdir) THEN
     237           0 :             IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
     238             : 
     239           0 :             IF (banddos%l_orb.AND.ANY((/banddos%alpha,banddos%beta,banddos%gamma/).NE.0.0)) THEN
     240           0 :                CALL abcrot2(atoms,banddos,noccbd,eigVecCoeffs,ispin) ! rotate ab-coeffs
     241           0 :                IF (PRESENT(orbcomp)) CALL orb_comp(ispin,ikpt,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
     242             :             END IF
     243             :          ENDIF
     244        1848 :          CALL calcDenCoeffs(atoms,sphhar,sym,we,noccbd,eigVecCoeffs,ispin,denCoeffs)
     245             : 
     246        1848 :          IF (noco%l_soc) CALL orbmom(atoms,noccbd,we,ispin,eigVecCoeffs,orb)
     247        1848 :          IF (input%l_f) CALL force%addContribsA21A12(input,atoms,dimension,sym,cell,oneD,enpara,&
     248          12 :                                                      usdus,eigVecCoeffs,noccbd,ispin,eig,we,results)
     249        1848 :          IF(l_coreSpec) CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,dimension%neigd,&
     250        1848 :                                           noccbd,results%ef,banddos%sig_dos,eig,we,eigVecCoeffs)
     251             :       END DO ! end loop over ispin
     252        1848 :       IF (noco%l_mperp) CALL denCoeffsOffdiag%calcCoefficients(atoms,sphhar,sym,eigVecCoeffs,we,noccbd)
     253             : 
     254        2456 :       IF ((banddos%dos.OR.banddos%vacdos.OR.input%cdinf).AND.(banddos%ndir.GT.0)) THEN
     255             :          ! since z is no longer an argument of cdninf sympsi has to be called here!
     256           0 :          CALL sympsi(lapw,jspin,sym,dimension,nbands,cell,eig,noco,dos%ksym(:,ikpt,jspin),dos%jsym(:,ikpt,jspin),zMat)
     257             :       END IF
     258             :    END DO ! end of k-point loop
     259             : 
     260             : #ifdef CPP_MPI
     261        1216 :    DO ispin = jsp_start,jsp_end
     262             :       CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,input,noco,ispin,regCharges,dos,&
     263        1216 :                        results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin),mcd,slab,orbcomp)
     264             :    END DO
     265             : #endif
     266             : 
     267             :    CALL cdnmt(mpi,input%jspins,atoms,sphhar,noco,jsp_start,jsp_end,&
     268         608 :                  enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt)
     269         608 :    IF (mpi%irank==0) THEN
     270         304 :       IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins)
     271         608 :       DO ispin = jsp_start,jsp_end
     272         304 :          IF (input%cdinf) THEN
     273           0 :             WRITE (6,FMT=8210) ispin
     274             : 8210        FORMAT (/,5x,'check continuity of cdn for spin=',i2)
     275           0 :             CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,den,ispin)
     276             :          END IF
     277         608 :          IF (input%l_f) CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,force,results)
     278             :       END DO
     279         304 :       CALL closeXMLElement('mtCharges')
     280             :    END IF
     281             : 
     282         608 :    CALL timestop("cdnval")
     283             : 
     284         608 : END SUBROUTINE cdnval
     285             : 
     286             : END MODULE m_cdnval

Generated by: LCOV version 1.13