LCOV - code coverage report
Current view: top level - main - cdngen.F90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 107 118 90.7 %
Date: 2024-04-24 04:44:14 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             : MODULE m_cdngen
       7             : #ifdef CPP_MPI
       8             :    USE mpi
       9             : #endif
      10             : CONTAINS
      11             : 
      12         682 : SUBROUTINE cdngen(eig_id,fmpi,input,banddos,sliceplot,vacuum,&
      13             :                   kpts,atoms,sphhar,stars,sym,juphon,gfinp,hub1inp,&
      14             :                   enpara,cell,noco,nococonv,vTot,results ,coreSpecInput,&
      15         682 :                   archiveType, xcpot,outDen,EnergyDen,greensFunction,hub1data,vxc,exc)
      16             : 
      17             :    !*****************************************************
      18             :    !    Charge density generator
      19             :    !    calls cdnval to generate the valence charge and the
      20             :    !    core routines for the core contribution
      21             :    !*****************************************************
      22             :    use m_types_vacdos
      23             :    use m_types_mcd
      24             :    use m_types_slab
      25             :    use m_types_orbcomp
      26             :    use m_types_jdos
      27             :    USE m_types
      28             :    USE m_constants
      29             :    USE m_juDFT
      30             :    !USE m_prpqfftmap
      31             :    USE m_cdnval
      32             :    USE m_plot
      33             :    USE m_cdn_io
      34             :    USE m_wrtdop
      35             :    USE m_cdntot
      36             :    USE m_qfix
      37             :    USE m_xmlOutput
      38             :    USE m_magMultipoles
      39             :    USE m_magmoments
      40             :    USE m_resMoms
      41             :    USE m_cdncore
      42             :    USE m_make_dos
      43             :    !USE m_Ekwritesl
      44             :    !USE m_banddos_io
      45             :    USE m_metagga
      46             :    !USE m_unfold_band_kpts
      47             :    USE m_denMultipoleExp
      48             :    use m_slater
      49             :    USE m_greensfPostProcess
      50             :    USE m_types_greensfContourData
      51             :    USE m_types_eigdos
      52             :    USE m_types_dos
      53             : 
      54             :    USE m_force_sf ! Klueppelberg (force level 3)
      55             : 
      56             :    IMPLICIT NONE
      57             : 
      58             :    ! Type instance arguments
      59             :    TYPE(t_results),INTENT(INOUT)    :: results
      60             :    TYPE(t_mpi),INTENT(IN)           :: fmpi
      61             : 
      62             :     
      63             :    TYPE(t_enpara),INTENT(INOUT)     :: enpara
      64             :    TYPE(t_banddos),INTENT(IN)       :: banddos
      65             :    TYPE(t_sliceplot),INTENT(IN)     :: sliceplot
      66             :    TYPE(t_input),INTENT(IN)         :: input
      67             :    TYPE(t_vacuum),INTENT(IN)        :: vacuum
      68             :    TYPE(t_noco),INTENT(IN)          :: noco
      69             :    TYPE(t_nococonv),INTENT(INOUT)   :: nococonv
      70             :    TYPE(t_sym),INTENT(IN)           :: sym
      71             :    TYPE(t_juphon),INTENT(IN)        :: juphon
      72             :    TYPE(t_stars),INTENT(IN)         :: stars
      73             :    TYPE(t_cell),INTENT(IN)          :: cell
      74             :    TYPE(t_kpts),INTENT(IN)          :: kpts
      75             :    TYPE(t_sphhar),INTENT(IN)        :: sphhar
      76             :    TYPE(t_atoms),INTENT(IN)         :: atoms
      77             :    TYPE(t_coreSpecInput),INTENT(IN) :: coreSpecInput
      78             :    TYPE(t_potden),INTENT(IN)        :: vTot
      79             :    TYPE(t_gfinp),INTENT(IN)         :: gfinp
      80             :    TYPE(t_hub1inp),INTENT(IN)       :: hub1inp
      81             :    TYPE(t_greensf),OPTIONAL,INTENT(INOUT)    :: greensFunction(:)
      82             :    TYPE(t_hub1data),OPTIONAL,INTENT(INOUT)    :: hub1data
      83             :    CLASS(t_xcpot),INTENT(IN)     :: xcpot
      84             :    TYPE(t_potden),INTENT(INOUT)     :: outDen, EnergyDen
      85             :    TYPE(t_potden),INTENT(INOUT),OPTIONAL:: vxc, exc
      86             : 
      87             :    !Scalar Arguments
      88             :    INTEGER, INTENT (IN)             :: eig_id, archiveType
      89             : 
      90             :    ! Local type instances
      91         682 :    TYPE(t_regionCharges)          :: regCharges
      92         682 :    TYPE(t_dos),TARGET             :: dos
      93         682 :    TYPE(t_vacdos),TARGET          :: vacdos
      94         682 :    TYPE(t_moments)                :: moments
      95         682 :    TYPE(t_mcd),TARGET             :: mcd
      96         682 :    TYPE(t_slab),TARGET            :: slab
      97         682 :    TYPE(t_orbcomp),TARGET         :: orbcomp
      98         682 :    TYPE(t_jDOS),TARGET            :: jDOS
      99         682 :    TYPE(t_cdnvalJob)       :: cdnvalJob
     100         682 :    TYPE(t_greensfImagPart) :: greensfImagPart
     101         682 :    TYPE(t_potden)          :: val_den, core_den
     102        2264 :    TYPE(t_greensfContourData) :: contour(gfinp%numberContours)
     103             : 
     104             : 
     105             :    !Local Scalars
     106             :    REAL                  :: fix, qtot, dummy,eFermiPrev
     107             :    INTEGER               :: jspin, ierr
     108             :    INTEGER               :: dim_idx
     109             :    INTEGER               :: i_gf,iContour,n
     110             : 
     111         682 :    TYPE(t_eigdos_list),allocatable :: eigdos(:)
     112             : 
     113             : #ifdef CPP_HDF
     114             :    INTEGER(HID_T)        :: banddosFile_id
     115             : #endif
     116             :    LOGICAL               :: l_error,Perform_metagga
     117             : 
     118             :    ! Initialization section
     119         682 :    CALL regCharges%init(input,atoms)
     120         682 :    CALL moments%init(fmpi,input,sphhar,atoms)
     121             :    !initalize data for DOS
     122       27210 :    if (noco%l_noco) results%eig(:,:,2)=results%eig(:,:,1)
     123         682 :    CALL dos%init(input,atoms,kpts,banddos,results%eig)
     124         682 :    CALL vacdos%init(input,atoms,kpts,banddos,results%eig)
     125         682 :    CALL mcd%init(banddos,input,atoms,kpts,results%eig)
     126         682 :    CALL slab%init(banddos,atoms,cell,input,kpts)
     127         682 :    CALL orbcomp%init(input,banddos,atoms,kpts,results%eig)
     128         682 :    CALL jDOS%init(input,banddos,atoms,kpts,results%eig)
     129             : 
     130         682 :    if (banddos%dos.or.banddos%band.or.input%cdinf) then
     131         100 :      allocate(eigdos(count((/banddos%dos.or.banddos%band.or.input%cdinf,banddos%vacdos,banddos%l_mcd,banddos%l_slab,banddos%l_orb,banddos%l_jDOS/))))
     132          10 :      n=2
     133          10 :      eigdos(1)%p=>dos
     134          10 :      if (banddos%vacdos) THEN
     135           0 :        eigdos(n)%p=>vacdos; n=n+1;
     136             :      endif
     137          10 :      if (banddos%l_mcd) THEN
     138           0 :        eigdos(n)%p=>mcd; n=n+1
     139             :      endif
     140          10 :      if (banddos%l_slab) THEN
     141           0 :        eigdos(n)%p=>slab; n=n+1
     142             :      endif
     143          10 :      if (banddos%l_orb) THEN
     144           2 :        eigdos(n)%p=>orbcomp; n=n+1
     145             :      endif
     146          10 :      if (banddos%l_jdos) eigdos(n)%p=>jDOS
     147             :    endif
     148             : 
     149             : 
     150             : 
     151         682 :    CALL outDen%init(stars,    atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
     152         682 :    CALL EnergyDen%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_EnergyDen)
     153             : 
     154             : 
     155         682 :    IF(PRESENT(greensFunction).AND.gfinp%n.GT.0) THEN
     156             :       !Only calculate the greens function when needed
     157          42 :       IF(ANY(greensFunction(:)%l_calc)) THEN
     158             :          !calculate all contours only once
     159          88 :          DO iContour = 1, gfinp%numberContours
     160          46 :             CALL contour(iContour)%init(gfinp%contour(iContour))
     161          88 :             CALL contour(iContour)%eContour(gfinp%contour(iContour),results%ef,fmpi%irank)
     162             :          ENDDO
     163             :       ENDIF
     164        1060 :       DO i_gf = 1, gfinp%n
     165        1018 :          IF(.NOT.greensFunction(i_gf)%l_calc) CYCLE
     166        1018 :          iContour = gfinp%elem(i_gf)%iContour
     167        1018 :          greensFunction(i_gf)%contour = contour(iContour)
     168        1060 :          CALL greensFunction(i_gf)%reset()
     169             :       ENDDO
     170          42 :       CALL greensfImagPart%init(gfinp,atoms,input,noco,ANY(greensFunction(:)%l_calc),SIZE(fmpi%k_list))
     171             :    ENDIF
     172             : 
     173         682 :    IF(fmpi%irank==0 .AND.PRESENT(hub1data)) THEN
     174        1705 :       hub1data%mag_mom = 0.0
     175     2675575 :       hub1data%cdn_atomic = 0.0
     176             :    ENDIF
     177             : 
     178             : 
     179         682 :    IF (fmpi%irank == 0) CALL openXMLElementNoAttributes('valenceDensity')
     180             : 
     181             :    !In a non-collinear calcuation where the off-diagonal part of the
     182             :    !density matrix in the muffin-tins is calculated, the a- and
     183             :    !b-coef. for both spins are needed at once. Thus, cdnval is only
     184             :    !called once and both spin directions are calculated in a single run.
     185         682 :    CALL timestart("cdngen: cdnval")
     186        1700 :    DO jspin = 1,merge(1,input%jspins,noco%l_mperp.OR.banddos%l_jDOS)
     187        1018 :       CALL cdnvalJob%init(fmpi,input,kpts,noco,results,jspin)
     188        1018 :       IF (sliceplot%slice) CALL cdnvalJob%select_slice(sliceplot,results,input,kpts,noco,jspin)
     189             :       CALL cdnval(eig_id,fmpi,kpts,jspin,noco,nococonv,input,banddos,cell,atoms,enpara,stars,vacuum,&
     190             :                   sphhar,sym,vTot ,cdnvalJob,outDen,regCharges,dos,vacdos,results,moments,gfinp,&
     191        1700 :                   hub1inp,hub1data,coreSpecInput,mcd,slab,orbcomp,jDOS,greensfImagPart)
     192             :    END DO
     193         682 :    CALL timestop("cdngen: cdnval")
     194             : 
     195         682 :    call val_den%copyPotDen(outDen)
     196             :    ! calculate kinetic energy density for MetaGGAs
     197         682 :    if(xcpot%exc_is_metagga()) then
     198             :       CALL calc_EnergyDen(eig_id, fmpi, kpts, noco, nococonv,input, banddos, cell, atoms, enpara, stars,&
     199           0 :                              vacuum,  sphhar, sym, gfinp, hub1inp, vTot,   results, EnergyDen)
     200             :    endif
     201             : 
     202         682 :    IF (banddos%dos.or.banddos%band.or.input%cdinf) THEN
     203          10 :       IF (fmpi%irank == 0) THEN
     204           5 :          CALL timestart("cdngen: dos")
     205             :          CALL make_dos(kpts,atoms,vacuum,input,banddos,&
     206           5 :                       sliceplot,noco,sym,cell,results,eigdos )
     207           5 :          CALL timestop("cdngen: dos")
     208             :       END IF
     209             :    END IF
     210             : 
     211         682 :    CALL cdntot(stars,nococonv,atoms,sym,vacuum,input,cell ,outDen,.TRUE.,qtot,dummy,fmpi,.TRUE.)
     212         682 :    IF (fmpi%irank.EQ.0) THEN
     213         341 :       CALL closeXMLElement('valenceDensity')
     214             :    END IF ! fmpi%irank = 0
     215             : 
     216         682 :    IF(PRESENT(greensFunction) .AND.gfinp%n.GT.0) THEN
     217          42 :       IF(greensfImagPart%l_calc) THEN
     218             :          CALL greensfPostProcess(greensFunction,greensfImagPart,atoms,kpts,cell,gfinp,input,sym,noco,fmpi,&
     219          42 :                                  nococonv,vTot,enpara,hub1inp,sphhar,hub1data,results)
     220             :       ELSE
     221           0 :          IF(fmpi%irank.EQ.0) THEN
     222           0 :             WRITE(oUnit,'(/,A)') "Green's Functions are not calculated: "
     223           0 :             WRITE(oUnit,'(A,f12.7,TR5,A,f12.7/)') "lastDistance: ", results%last_distance,&
     224           0 :                                                   "minCalcDistance: ", gfinp%minCalcDistance
     225             :          ENDIF
     226             :       ENDIF
     227             :    ENDIF
     228             : 
     229         682 :    IF (banddos%vacdos.or.banddos%dos.or.banddos%band.or.input%cdinf) THEN
     230          10 :       CALL juDFT_end("Charge density postprocessing done.",fmpi%irank)
     231             :    END IF
     232             : 
     233         672 :    IF (sliceplot%slice) THEN
     234           2 :       IF (fmpi%irank == 0) THEN
     235           4 :          IF(any(noco%l_alignMT)) CALL juDFT_error("Relaxation of SQA and sliceplot not implemented. To perfom a sliceplot of the correct cdn deactivate realaxation.", calledby = "cdngen" )
     236             :          CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
     237           1 :                            0,-1.0,0.0,-1.0,-1.0,.FALSE.,outDen,inFilename='cdn_slice')
     238             :       END IF
     239           2 :       call outDen%distribute(fmpi%mpi_comm)
     240           2 :       CALL juDFT_end("slice OK",fmpi%irank)
     241             :    END IF
     242             : 
     243             :    !IF (sliceplot%iplot.NE.0) THEN
     244             :    !   CALL makeplots(stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, noco,nococonv, outDen, PLOT_OUTDEN_Y_CORE, sliceplot)
     245             :    !END IF
     246             : 
     247         670 :    CALL timestart("cdngen: cdncore")
     248         670 :    if(xcpot%exc_is_MetaGGA()) then
     249             :       CALL cdncore(fmpi ,input,vacuum,noco,nococonv,sym,&
     250           0 :                    stars,cell,sphhar,atoms,vTot,outDen,moments,results, EnergyDen)
     251             :    else
     252             :       CALL cdncore(fmpi ,input,vacuum,noco,nococonv,sym,&
     253         670 :                    stars,cell,sphhar,atoms,vTot,outDen,moments,results)
     254             :    endif
     255         670 :    call core_den%subPotDen(outDen, val_den)
     256         670 :    CALL timestop("cdngen: cdncore")
     257             : 
     258             :    IF(.FALSE.) CALL denMultipoleExp(input, fmpi, atoms, sphhar, stars, sym, juphon, cell,   outDen) ! There should be a switch in the inp file for this
     259         670 :    IF(fmpi%irank.EQ.0) THEN
     260         335 :       IF(input%lResMax>0) CALL resMoms(sym,input,atoms,sphhar,noco,nococonv,outDen,moments%rhoLRes) ! There should be a switch in the inp file for this
     261             :    END IF
     262             : 
     263         670 :    IF(atoms%n_opc>0) THEN
     264          72 :       do jspin=1, input%jspins
     265          72 :          call slater(input,jspin,atoms,vTot%mt(:,0,:,jspin),l_write=fmpi%irank==0)
     266             :       enddo
     267             :    ENDIF
     268             : 
     269         670 :    CALL enpara%calcOutParams(input,atoms,vacuum,regCharges)
     270             : 
     271         670 :    IF (fmpi%irank == 0) CALL openXMLElementNoAttributes('allElectronCharges')
     272         670 :    CALL qfix(fmpi,stars,nococonv,atoms,sym,vacuum,sphhar,input,cell ,outDen,noco%l_noco,.TRUE.,l_par=.TRUE.,force_fix=.TRUE.,fix=fix)
     273         670 :    IF (fmpi%irank == 0) THEN
     274         335 :       CALL closeXMLElement('allElectronCharges')
     275             : 
     276         335 :       IF (input%jspins == 2) THEN
     277             :          !Calculate and write out spin densities at the nucleus and magnetic moments in the spheres
     278         196 :          CALL spinMoments(input,atoms,noco,nococonv,den=outDen)
     279         196 :          CALL orbMoments(input,atoms,noco,nococonv,moments)
     280             : 
     281         196 :          if (sym%nop==1.and..not.input%film) call magMultipoles(sym,juphon,stars, atoms,cell, sphhar, vacuum, input, noco,nococonv,outden)
     282             :          !Generate and save the new nocoinp file if the directions of the local
     283             :          !moments are relaxed or a constraint B-field is calculated.
     284             :       END IF
     285             :    END IF ! fmpi%irank == 0
     286             :    Perform_metagga = Allocated(Energyden%Mt) &
     287         670 :                    .And. (Xcpot%Exc_is_metagga() .Or. Xcpot%Vx_is_metagga())
     288             :    If(Perform_metagga) Then
     289           0 :      IF(any(noco%l_alignMT)) CALL juDFT_error("Relaxation of SQA and metagga not implemented.", calledby = "cdngen" )
     290             :      CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
     291           0 :                            0,-1.0,0.0,-1.0,-1.0,.FALSE.,core_den,inFilename='cdnc')
     292             :    endif
     293             : 
     294             : #ifdef CPP_MPI
     295         670 :    CALL MPI_BCAST(nococonv%alph,atoms%ntype,MPI_DOUBLE_PRECISION,0,fmpi%mpi_comm,ierr)
     296         670 :    CALL MPI_BCAST(nococonv%beta,atoms%ntype,MPI_DOUBLE_PRECISION,0,fmpi%mpi_comm,ierr)
     297         670 :    CALL MPI_BCAST(nococonv%b_con,atoms%ntype*2,MPI_DOUBLE_PRECISION,0,fmpi%mpi_comm,ierr)
     298         670 :    CALL MPI_BCAST(nococonv%qss,3,MPI_DOUBLE_PRECISION,0,fmpi%mpi_comm,ierr)
     299             : #endif
     300         670 :    CALL outDen%distribute(fmpi%mpi_comm)
     301             : 
     302             :    ! Klueppelberg (force level 3)
     303         670 :    IF (input%l_f.AND.(input%f_level.GE.3).AND.(fmpi%irank.EQ.0)) THEN
     304           2 :       DO jspin = 1,input%jspins ! jsp_start, jsp_end
     305           2 :          CALL force_sf_mt(atoms,sphhar,jspin,jspin,fmpi,vtot%mt(:,0:,:,jspin),exc%mt(:,0:,:,1),vxc%mt(:,0:,:,:),outDen%mt(:,0:,:,:),sym,cell)
     306             :       END DO
     307             :    END IF
     308             : 
     309        1340 : END SUBROUTINE cdngen
     310             : 
     311             : END MODULE m_cdngen

Generated by: LCOV version 1.14