cdncore.F90 5.46 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 12
SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
                   stars,cell,sphhar,atoms,vTot,outDen,moments,results)
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
   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
46 47 48 49 50

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

51 52 53 54
   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)
55 56 57 58 59

   results%seigc = 0.0
   IF (mpi%irank.EQ.0) THEN
      DO jspin = 1,input%jspins
         DO n = 1,atoms%ntype
60
            moments%svdn(n,jspin) = outDen%mt(1,0,n,jspin) / (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
         END DO
      END DO
   END IF

   IF (input%kcrel.EQ.0) THEN
      ! Generate input file ecore for subsequent GW calculation
      ! 11.2.2004 Arno Schindlmayr
      IF ((input%gw.eq.1 .or. input%gw.eq.3).AND.(mpi%irank.EQ.0)) THEN
         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
         IF (mpi%irank.EQ.0) THEN
            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

85 86 87 88 89 90
   !add in core density
   IF (mpi%irank.EQ.0) THEN
      IF (input%kcrel.EQ.0) THEN
         DO jspin = 1,input%jspins
            CALL cored(input,jspin,atoms,outDen%mt,dimension,sphhar,vTot%mt(:,0,:,jspin), qint,rh,tec,seig)
            rhTemp(:,:,jspin) = rh(:,:,jspin)
91
            results%seigc = results%seigc + seig
92 93 94 95
         END DO
      ELSE
         CALL coredr(input,atoms,seig, outDen%mt,dimension,sphhar,vTot%mt(:,0,:,:),qint,rh)
         results%seigc = results%seigc + seig
96
      END IF
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
   END IF
   DO jspin = 1,input%jspins
      IF (mpi%irank.EQ.0) THEN
         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
      IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
         IF (jspin.EQ.2) THEN
            !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
               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
               !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
               outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(iType))*sin(noco%beta(iType)),0.0)
               !imaginary part rho_21
               outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(iType))*sin(noco%beta(iType)))
121
            END DO
122
            !pk non-collinear (end)
123
         END IF
124 125 126 127 128 129 130 131
      ELSE
         IF (input%ctail) THEN
            !+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)
         ELSE IF (mpi%irank.EQ.0) THEN
            DO iType = 1,atoms%ntype
132
               outDen%pw(1,jspin) = outDen%pw(1,jspin) + qint(iType,jspin) / (input%jspins * cell%volint)
133
            END DO
134
         END IF
135 136
      END IF
   END DO
137 138 139 140 141 142 143 144 145 146 147

   IF (input%kcrel.EQ.0) THEN
      IF (mpi%irank.EQ.0) THEN
         CALL writeCoreDensity(input,atoms,dimension,rhTemp,tec,qint)
      END IF
      IF ((input%gw.eq.1 .or. input%gw.eq.3).AND.(mpi%irank.EQ.0)) CLOSE(15)
   END IF

END SUBROUTINE cdncore

END MODULE m_cdncore