flipcdn.f90 11 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_flipcdn
8
!     *******************************************************
9 10
!     this subroutine reads the charge density and flips the
!     magnetic moment within the m.t.sphere for each atom
11 12
!     according to the variable nflip. This variable is read in
!     the main program
Robin Hilgers's avatar
Robin Hilgers committed
13
!  TODO; (Test)           nflip = -1 : flip spin in sphere
Robin Hilgers's avatar
Robin Hilgers committed
14
!  TODO:           nflip = -2 : scale spin by bmu(n)
15 16
!             nflip = any: no spin flip
!                            r.pentcheva,kfa,Feb'96
Daniel Wortmann's avatar
Daniel Wortmann committed
17 18
!
!     Extension to multiple U per atom type by G.M. 2017
19 20
!
!     Removed integer nflip switch and added angles phi/theta
Robin Hilgers's avatar
Robin Hilgers committed
21
!     (and an additional spin scale switch)
22 23
!     which defines spin flip for each atom individually.
!     => Magnetisation axis can now be chosen independet
Robin Hilgers's avatar
Robin Hilgers committed
24
!     of spin quantization axis.
25
!     R. Hilgers, Okt. 2019
26
!     *******************************************************
27 28
CONTAINS

29
SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,phi,theta,optDen)
30
   USE m_rotdenmat
31 32 33 34 35 36 37 38
   USE m_constants
   USE m_cdn_io
   USE m_types

   IMPLICIT NONE

   TYPE(t_stars),INTENT(IN)    :: stars
   TYPE(t_vacuum),INTENT(IN)   :: vacuum
39
   TYPE(t_atoms),INTENT(IN)    :: atoms
40
   TYPE(t_sphhar),INTENT(IN)   :: sphhar
41
   TYPE(t_input),INTENT(IN)    :: input
42
   TYPE(t_sym),INTENT(IN)      :: sym
43
   TYPE(t_noco),INTENT(IN)     :: noco
44 45
   TYPE(t_oneD),INTENT(IN)     :: oneD
   TYPE(t_cell),INTENT(IN)     :: cell
46 47 48
   REAL, OPTIONAL, INTENT(IN)  :: phi(atoms%ntype)
   REAL, OPTIONAL, INTENT(IN)  :: theta(atoms%ntype)
   TYPE(t_potden), OPTIONAL,INTENT(INOUT) :: optDen
49 50 51 52 53

   ! Local type instance
   TYPE(t_potden)            :: den

   ! Local Scalars
Robin Hilgers's avatar
Robin Hilgers committed
54
   COMPLEX                   :: rhodummy, imPart12, realPart12
Robin Hilgers's avatar
Robin Hilgers committed
55
   REAL                      :: rhodumms,fermiEnergyTemp, realPart1, realPart2, imPart1,imPart2, rhodummyR, rotAnglePhi(atoms%ntype),rotAngleTheta(atoms%ntype),zeros(atoms%ntype)
56
   INTEGER                   :: i,nt,j,lh,na,mp,ispin,urec,itype,m,i_u,k
57
   INTEGER                   :: archiveType
58
   LOGICAL                   :: n_exist,l_qfix,l_error, l_flip(atoms%ntype), scaleSpin(atoms%ntype)
59 60 61
   ! Local Arrays
   CHARACTER(len=80), ALLOCATABLE :: clines(:)

62
   zeros=0.0
63

64
!Flipcdn by optional given angle if lflip is false but routine is called.
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
   DO k=1, atoms%ntype
      IF(.NOT.input%lflip.AND.PRESENT(phi).AND.present(theta)) THEN
         rotAnglePhi(k)=phi(k)
         rotAngleTheta(k)=theta(k)
         scaleSpin(k)=.FALSE.
      ELSE IF (input%lflip) THEN
!Rotation triggerd by lflip.
         rotAnglePhi(k)=atoms%flipSpinPhi(k)
         rotAngleTheta(k)=atoms%flipSpinTheta(k)
         scaleSpin(k)=atoms%flipSpinScale(k)
      ELSE
         CALL judft_error("You shouldn't be here. There went something wrong.",calledby="flipcdn")
      END IF
   END DO



82
   DO itype=1, atoms%ntype
83
      l_flip(itype)=MERGE(.TRUE.,.FALSE.,(rotAnglePhi(itype).NE.0.0) .OR.(rotAngleTheta(itype).NE.0.0))
84
   END DO
85
   !rot_den_mat(alph,beta,rho11,rho22,rho21)
86
   archiveType = MERGE(CDN_ARCHIVE_TYPE_NOCO_const,CDN_ARCHIVE_TYPE_CDN1_const,noco%l_noco)
87
   IF(.NOT.PRESENT(optDen)) THEN
88
      CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
89
      ! read the charge density
90 91
      CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
                       CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den)
92 93 94
   ELSE
      den=optDen
   END IF
95

96
   ! flip cdn for each atom with rotation angles given
97 98
   na = 1
   DO itype = 1, atoms%ntype
99
      IF (l_flip(itype).AND.(.NOT.scaleSpin(itype))) THEN
100
         ! spherical and non-spherical m.t. charge density
101
!$OMP parallel private(rhodummy,j,rhodummyR)
Robin Hilgers's avatar
Robin Hilgers committed
102
!$OMP do
103
	DO lh = 0,sphhar%nlh(sym%ntypsy(na))
104
            DO j = 1,atoms%jri(itype)
105
                IF (noco%l_mtNocoPot) THEN
106
                   rhodummy=CMPLX(den%mt(j,lh,itype,3),den%mt(j,lh,itype,4))
Robin Hilgers's avatar
Robin Hilgers committed
107 108
                   CALL rot_den_mat(zeros(itype),rotAngleTheta(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy)
                   CALL rot_den_mat(rotAnglePhi(itype),zeros(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy)
109
                   den%mt(j,lh,itype,3)=REAL(rhodummy)
110
		   den%mt(j,lh,itype,4)=AIMAG(rhodummy)
111
                ELSE
112
                   IF (rotAngleTheta(itype).EQ.(pimach()).AND.rotAnglePhi(itype).EQ.0) THEN
113 114 115 116
                      rhodummyR = den%mt(j,lh,itype,1)
                      den%mt(j,lh,itype,1) = den%mt(j,lh,itype,input%jspins)
                      den%mt(j,lh,itype,input%jspins) = rhodummyR
                   ELSE
117
                      !Since in non-noco case the den-matrices are only initialized with two diagonal components we cannot perform flips where off-diagonal elements arise in non-noco case => Only rotations by theta=Pi/2 are allowed.
118
                      CALL judft_error("l_mtNocoPot=F in combination with spin flips different from flipSpinTheta=Pi and flipSpinPhi=0 is currently not supported.",&
119 120 121
                         calledby="flipcdn")
                   END IF
                END IF
122 123
            END DO
         END DO
124
!$OMP end do
Robin Hilgers's avatar
Robin Hilgers committed
125
!$OMP end parallel
126 127
      ELSE IF (l_flip(itype).AND.scaleSpin(itype)) THEN
         IF((rotAngleTheta(itype).NE.(pimach()) .OR.rotAnglePhi(itype).NE.0.0)) CALL judft_error("Spinscaling in combination with flipSpin is currently only implemented using flipSpinTheta=Pi and flipSpinPhi=0.0.",calledby="flipcdn")
128
         DO lh = 0,sphhar%nlh(sym%ntypsy(na))
129 130 131 132 133 134 135 136 137 138 139
            DO j = 1,atoms%jri(itype)
               rhodummy = den%mt(j,lh,itype,1) + den%mt(j,lh,itype,input%jspins)
               rhodumms = den%mt(j,lh,itype,1) - den%mt(j,lh,itype,input%jspins)
               den%mt(j,lh,itype,1) = 0.5 * (rhodummy + atoms%bmu(itype)*rhodumms)
               den%mt(j,lh,itype,input%jspins) = 0.5 * (rhodummy - atoms%bmu(itype)*rhodumms )
            END DO
         END DO
      END IF
         na = na + atoms%neq(itype)
   END DO

140
   IF (input%l_onlyMtStDen) THEN
141 142
!!This Segment takes care that no interstitial magnetization is written in the the density. Meaning: Off diagonal elements of density matrix set to 0 and diagonal elements of density matrix are equal to their mean value.
      den%pw(:,2)=(den%pw(:,1)+den%pw(:,2))*0.5 !mean value
Robin Hilgers's avatar
Robin Hilgers committed
143 144
      den%pw(:,1)=den%pw(:,2)
      IF (noco%l_noco) THEN
145
         den%pw(:,3)=CMPLX(0.0,0.0)
146 147 148 149
      END IF
   END IF


150
   ! for LDA+U: flip density matrix
151 152
   IF (ANY(den%mmpMat(:,:,:,:).NE.0.0).AND.atoms%n_u+atoms%n_hia>0) THEN
      DO i_u = 1, atoms%n_u+atoms%n_hia
153
         itype = atoms%lda_u(i_u)%atomType
154
        IF (l_flip(itype).AND.(.NOT.scaleSpin(itype))) THEN
155 156
            DO m = -3,3
               DO mp = -3,3
157
                  IF (noco%l_mtNocoPot) THEN
158
! Since den%mmpMat is complex but rot_den_mat can only handle real values as diagonals of den_mat a splitting of den%mmpMat in real and imaginary part is performed. Rotations are performed seperatly and added up afterwards.
159 160 161 162 163 164
                    realPart1=REAL(den%mmpMat(m,mp,i_u,1))
                    realPart2=REAL(den%mmpMat(m,mp,i_u,2))
                    realPart12=CMPLX(REAL(den%mmpMat(m,mp,i_u,3)),0)
                    imPart1=AIMAG(den%mmpMat(m,mp,i_u,1))
                    imPart2=AIMAG(den%mmpMat(m,mp,i_u,2))
                    imPart12=CMPLX(0,AIMAG(den%mmpMat(m,mp,i_u,3)))
165
                    CALL rot_den_mat(rotAnglePhi(itype),rotAngleTheta(itype),realPart1,realPart2,&
166
                      realPart12)
167
                    CALL rot_den_mat(rotAnglePhi(itype),rotAngleTheta(itype),imPart1,imPart2,&
168 169 170 171 172
                      imPart12)
                    den%mmpMat(m,mp,i_u,1)=CMPLX(realPart1,imPart1)
                    den%mmpMat(m,mp,i_u,2)=CMPLX(realPart2,imPart2)
                    den%mmpMat(m,mp,i_u,3)=realPart12+imPart12
                 ELSE
173
                    IF (rotAngleTheta(itype).EQ.(pimach()).AND.rotAnglePhi(itype).EQ.0) THEN
174 175 176 177
                       rhodummyR = den%mmpMat(m,mp,i_u,1)
                       den%mmpMat(m,mp,i_u,1) = den%mmpMat(m,mp,i_u,input%jspins)
                       den%mmpMat(m,mp,i_u,input%jspins) = rhodummyR
                    ELSE
178
                       !Since in non-noco case the den-matrices are only initialized with two diagonal components we cannot perform flips where off-diagonal elements arise in non-noco case => Only rotations by Pi degrees are allowed.
179
                       CALL judft_error("l_mtNocoPot=F in combination with spin flips different from flipSpinTheta=Pi, flipSpinPhi=0 is currently not supported.",&
180 181 182 183 184
                         calledby="flipcdn")
                    END IF
                 END IF

END DO
185
            END DO
186
          ELSE IF (l_flip(itype).AND.(scaleSpin(itype))) THEN
187 188
            DO m = -3,3
               DO mp = -3,3
189
                  IF((rotAngleTheta(itype).NE.pimach() .OR.rotAnglePhi(itype).NE.0.0)) CALL judft_error("Spinscaling in combination with flipSpin is currently only implemented using flipSpinTheta=Pi and flipSpinPhi=0.0",calledby="flipcdn")
190 191 192 193 194 195 196 197 198 199 200
                  rhodummy = den%mmpMat(m,mp,i_u,1) + den%mmpMat(m,mp,i_u,input%jspins)
                  rhodumms = den%mmpMat(m,mp,i_u,1) - den%mmpMat(m,mp,i_u,input%jspins)
                  den%mmpMat(m,mp,i_u,1) = 0.5 * (rhodummy + atoms%bmu(itype) * rhodumms)
                  den%mmpMat(m,mp,i_u,input%jspins) = 0.5 * (rhodummy - atoms%bmu(itype) * rhodumms)
               END DO
            END DO
         END IF
      END DO
   END IF

   ! write the spin-polarized density
201
   CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
202
                     1,-1.0,0.0,.FALSE.,den)
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220

   ! read enpara and  flip lines
   INQUIRE(file='enpara',exist=n_exist)
   IF (n_exist) THEN
      OPEN(40,file ='enpara',status='old',form='formatted')

      j = 2
      DO itype = 1, atoms%ntype
         j = j + 1
         IF (atoms%nlo(itype)>0) j = j + 2
      END DO
      IF (input%film) j = j + 1
      ALLOCATE (clines(2*j))
      DO i = 1, 2*j
         READ (40,'(a)') clines(i)
      END DO

      REWIND 40
221
      i = 0
222 223 224 225 226 227 228
      DO ispin = 1,input%jspins
         i = i + 2
         WRITE (40,'(a)') TRIM(clines(i-1))
         WRITE (40,'(a)') TRIM(clines(i))
         DO itype = 1, atoms%ntype
            i = i + 1
            m = i
229
            IF (l_flip(itype)) m = MOD(i+j,2*j)
230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
            IF (m==0) m = 2*j
            WRITE (40,'(a)') TRIM(clines(m))
            IF (atoms%nlo(itype)>0) THEN
               WRITE (40,'(a)') TRIM(clines(m+1))
               WRITE (40,'(a)') TRIM(clines(m+2))
               i = i + 2
            END IF
         END DO
         IF (input%film) THEN
            i = i + 1
            WRITE (40,'(a)') TRIM(clines(i))
         END IF
      END DO
      DEALLOCATE (clines)
      CLOSE(40)
   END IF

247

248
IF(PRESENT(optDen)) optDen=den
249

250 251 252
END SUBROUTINE flipcdn

END MODULE m_flipcdn