cdncore.F90 6.02 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

7 8 9 10
MODULE m_cdncore

CONTAINS

11
SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
Matthias Redies's avatar
Matthias Redies committed
12
                   stars,cell,sphhar,atoms,vTot,outDen,moments,results, EnergyDen)
13 14

   USE m_constants
Matthias Redies's avatar
Matthias Redies committed
15
   USE m_judft
16 17 18 19 20 21 22 23 24 25 26 27 28 29
   USE m_cdn_io
   USE m_cdnovlp
   USE m_cored
   USE m_coredr
   USE m_types
   USE m_xmlOutput
   USE m_magMoms
   USE m_orbMagMoms
#ifdef CPP_MPI
   USE m_mpi_bc_coreden
#endif

   IMPLICIT NONE

30

Matthias Redies's avatar
Matthias Redies committed
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
   TYPE(t_mpi),        INTENT(IN)              :: mpi
   TYPE(t_dimension),  INTENT(IN)              :: dimension
   TYPE(t_oneD),       INTENT(IN)              :: oneD
   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_sphhar),     INTENT(IN)              :: sphhar
   TYPE(t_atoms),      INTENT(IN)              :: atoms
   TYPE(t_potden),     INTENT(IN)              :: vTot
   TYPE(t_potden),     INTENT(INOUT)           :: outDen
   TYPE(t_moments),    INTENT(INOUT)           :: moments
   TYPE(t_results),    INTENT(INOUT)           :: results
   TYPE(t_potden),     INTENT(INOUT), OPTIONAL :: EnergyDen
47 48 49 50 51

   INTEGER                          :: jspin, n, iType
   REAL                             :: seig, rhoint, momint
   LOGICAL, PARAMETER               :: l_st=.FALSE.

52 53 54 55
   REAL                             :: rh(dimension%msh,atoms%ntype,input%jspins)
   REAL                             :: qint(atoms%ntype,input%jspins)
   REAL                             :: tec(atoms%ntype,input%jspins)
   REAL                             :: rhTemp(dimension%msh,atoms%ntype,input%jspins)
56

Matthias Redies's avatar
Matthias Redies committed
57

58
   results%seigc = 0.0
Matthias Redies's avatar
Matthias Redies committed
59
   IF (mpi%irank==0) THEN
60 61
      DO jspin = 1,input%jspins
         DO n = 1,atoms%ntype
62
            moments%svdn(n,jspin) = outDen%mt(1,0,n,jspin) / (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
63 64 65 66
         END DO
      END DO
   END IF

Matthias Redies's avatar
Matthias Redies committed
67
   IF (input%kcrel==0) THEN
68 69
      ! Generate input file ecore for subsequent GW calculation
      ! 11.2.2004 Arno Schindlmayr
Matthias Redies's avatar
Matthias Redies committed
70
      IF ((input%gw==1 .or. input%gw==3).AND.(mpi%irank==0)) THEN
71 72 73 74 75 76 77
         OPEN (15,file='ecore',status='unknown', action='write',form='unformatted')
      END IF

      rh = 0.0
      tec = 0.0
      qint = 0.0
      IF (input%frcor) THEN
Matthias Redies's avatar
Matthias Redies committed
78
         IF (mpi%irank==0) THEN
79 80 81 82 83 84 85 86
            CALL readCoreDensity(input,atoms,dimension,rh,tec,qint)
         END IF
#ifdef CPP_MPI
         CALL mpi_bc_coreDen(mpi,atoms,input,dimension,rh,tec,qint)
#endif
      END IF
   END IF

87
   !add in core density
Matthias Redies's avatar
Matthias Redies committed
88 89
   IF (mpi%irank==0) THEN
      IF (input%kcrel==0) THEN
90
         DO jspin = 1,input%jspins
91 92 93 94 95 96
            IF(PRESENT(EnergyDen)) THEN
               CALL cored(input,jspin,atoms,outDen%mt,dimension,sphhar,vTot%mt(:,0,:,jspin), qint,rh ,tec,seig, EnergyDen%mt)
            ELSE
               CALL cored(input,jspin,atoms,outDen%mt,dimension,sphhar,vTot%mt(:,0,:,jspin), qint,rh ,tec,seig)
            ENDIF

97
            rhTemp(:,:,jspin) = rh(:,:,jspin)
98
            results%seigc = results%seigc + seig
99 100
         END DO
      ELSE
101
         IF(PRESENT(EnergyDen)) call juDFT_error("Energyden not implemented for relativistic core calculations")
102 103
         CALL coredr(input,atoms,seig, outDen%mt,dimension,sphhar,vTot%mt(:,0,:,:),qint,rh)
         results%seigc = results%seigc + seig
104
      END IF
105 106
   END IF
   DO jspin = 1,input%jspins
Matthias Redies's avatar
Matthias Redies committed
107
      IF (mpi%irank==0) THEN
108 109 110 111
         DO n = 1,atoms%ntype
            moments%stdn(n,jspin) = outDen%mt(1,0,n,jspin) / (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
         END DO
      END IF
Matthias Redies's avatar
Matthias Redies committed
112 113 114 115
      IF ((noco%l_noco).AND.(mpi%irank==0)) THEN
         IF (jspin==2) THEN

            IF(PRESENT(EnergyDen)) call juDFT_error("Energyden not implemented for noco")
116 117 118 119 120
            !pk non-collinear (start)
            !add the coretail-charge to the constant interstitial
            !charge (star 0), taking into account the direction of
            !magnetisation of this atom
            DO iType = 1,atoms%ntype
Matthias Redies's avatar
Matthias Redies committed
121 122
               rhoint = (qint(iType,1) + qint(iType,2)) /(cell%volint * input%jspins * 2.0)
               momint = (qint(iType,1) - qint(iType,2)) /(cell%volint * input%jspins * 2.0)
123 124 125 126 127
               !rho_11
               outDen%pw(1,1) = outDen%pw(1,1) + rhoint + momint*cos(noco%beta(iType))
               !rho_22
               outDen%pw(1,2) = outDen%pw(1,2) + rhoint - momint*cos(noco%beta(iType))
               !real part rho_21
Matthias Redies's avatar
Matthias Redies committed
128
               outDen%pw(1,3) = outDen%pw(1,3) + cmplx( 0.5*momint *cos(noco%alph(iType))*sin(noco%beta(iType)),&
129
               !imaginary part rho_21
Matthias Redies's avatar
Matthias Redies committed
130
                                                       -0.5*momint *sin(noco%alph(iType))*sin(noco%beta(iType)))
131
            END DO
132
            !pk non-collinear (end)
133
         END IF
134 135
      ELSE
         IF (input%ctail) THEN
Matthias Redies's avatar
Matthias Redies committed
136
            IF(PRESENT(EnergyDen)) call juDFT_error("Energyden not implemented for ctail")
137 138 139 140
            !+gu hope this works as well
            CALL cdnovlp(mpi,sphhar,stars,atoms,sym,dimension,vacuum,&
                         cell,input,oneD,l_st,jspin,rh(:,:,jspin),&
                         outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz)
Matthias Redies's avatar
Matthias Redies committed
141
         ELSE IF (mpi%irank==0) THEN
142
            DO iType = 1,atoms%ntype
143
               outDen%pw(1,jspin) = outDen%pw(1,jspin) + qint(iType,jspin) / (input%jspins * cell%volint)
144
            END DO
145
         END IF
146 147
      END IF
   END DO
148

Matthias Redies's avatar
Matthias Redies committed
149 150
   IF (input%kcrel==0) THEN
      IF (mpi%irank==0) THEN
151 152
         CALL writeCoreDensity(input,atoms,dimension,rhTemp,tec,qint)
      END IF
Matthias Redies's avatar
Matthias Redies committed
153
      IF ((input%gw==1 .or. input%gw==3).AND.(mpi%irank==0)) CLOSE(15)
154 155 156 157 158
   END IF

END SUBROUTINE cdncore

END MODULE m_cdncore