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_flipcdn
8 : ! *******************************************************
9 : ! this subroutine reads the charge density and flips the
10 : ! magnetic moment within the m.t.sphere for each atom
11 : ! according to the variable nflip. This variable is read in
12 : ! the main program
13 : ! TODO; (Test) nflip = -1 : flip spin in sphere
14 : ! TODO: nflip = -2 : scale spin by bmu(n)
15 : ! nflip = any: no spin flip
16 : ! r.pentcheva,kfa,Feb'96
17 : !
18 : ! Extension to multiple U per atom type by G.M. 2017
19 : !
20 : ! Removed integer nflip switch and added angles phi/theta
21 : ! (and an additional spin scale switch)
22 : ! which defines spin flip for each atom individually.
23 : ! => Magnetisation axis can now be chosen independet
24 : ! of spin quantization axis.
25 : ! R. Hilgers, Okt. 2019
26 : ! *******************************************************
27 : CONTAINS
28 :
29 14 : SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco ,cell,phi,theta,optDen,toGlobal)
30 : !USE m_rotdenmat
31 : USE m_rotMMPmat
32 : USE m_constants
33 : USE m_cdn_io
34 : USE m_types
35 :
36 : IMPLICIT NONE
37 :
38 : TYPE(t_stars),INTENT(IN) :: stars
39 : TYPE(t_vacuum),INTENT(IN) :: vacuum
40 : TYPE(t_atoms),INTENT(IN) :: atoms
41 : TYPE(t_sphhar),INTENT(IN) :: sphhar
42 : TYPE(t_input),INTENT(IN) :: input
43 : TYPE(t_sym),INTENT(IN) :: sym
44 : TYPE(t_noco),INTENT(IN) :: noco
45 :
46 : TYPE(t_cell),INTENT(IN) :: cell
47 : REAL, OPTIONAL, INTENT(IN) :: phi(atoms%ntype)
48 : REAL, OPTIONAL, INTENT(IN) :: theta(atoms%ntype)
49 : TYPE(t_potden), OPTIONAL,INTENT(INOUT) :: optDen
50 : LOGICAL,OPTIONAL,INTENT(IN) :: toGlobal
51 :
52 : ! Local type instance
53 14 : TYPE(t_potden) :: den
54 56 : TYPE(t_nococonv) :: nococonv
55 :
56 : ! Local Scalars
57 : COMPLEX :: rhodummy, imPart12, realPart12
58 14 : REAL :: rhodumms,fermiEnergyTemp, realPart1, realPart2, imPart1,imPart2, rhodummyR, rotAnglePhi(atoms%ntype),rotAngleTheta(atoms%ntype),zeros(atoms%ntype)
59 : REAL :: tempDistance
60 : INTEGER :: i,nt,j,lh,na,mp,ispin,urec,itype,m,i_u,k,l
61 : INTEGER :: archiveType
62 14 : LOGICAL :: n_exist,l_qfix,l_error, l_flip(atoms%ntype), scaleSpin(atoms%ntype),opt
63 : ! Local Arrays
64 14 : CHARACTER(len=80), ALLOCATABLE :: clines(:)
65 14 : REAL,ALLOCATABLE :: mt_tmp(:,:,:,:)
66 14 : COMPLEX,ALLOCATABLE :: mmpMat_tmp(:,:,:,:)
67 40 : zeros=0.0
68 :
69 : !Flipcdn by optional given angle if lflip is false but routine is called.
70 40 : DO k=1, atoms%ntype
71 40 : IF(.NOT.input%lflip.AND.PRESENT(phi).AND.present(theta)) THEN
72 24 : rotAnglePhi(k)=phi(k)
73 24 : rotAngleTheta(k)=theta(k)
74 24 : scaleSpin(k)=.FALSE.
75 2 : ELSE IF (input%lflip) THEN
76 : !Rotation triggerd by lflip.
77 2 : rotAnglePhi(k)=atoms%flipSpinPhi(k)
78 2 : rotAngleTheta(k)=atoms%flipSpinTheta(k)
79 2 : scaleSpin(k)=atoms%flipSpinScale(k)
80 : ELSE
81 0 : CALL judft_error("You shouldn't be here. There went something wrong.",calledby="flipcdn")
82 : END IF
83 : END DO
84 :
85 :
86 40 : DO itype=1, atoms%ntype
87 61 : l_flip(itype)=MERGE(.TRUE.,.FALSE.,(rotAnglePhi(itype).NE.0.0) .OR.(rotAngleTheta(itype).NE.0.0))
88 : END DO
89 :
90 : !rot_den_mat(alph,beta,rho11,rho22,rho21)
91 14 : IF (any(noco%l_unrestrictMT)) THEN
92 14 : archiveType=CDN_ARCHIVE_TYPE_FFN_const
93 0 : ELSE IF (noco%l_noco) THEN
94 0 : archiveType=CDN_ARCHIVE_TYPE_NOCO_const
95 : ELSE
96 0 : archiveType=CDN_ARCHIVE_TYPE_CDN1_const
97 : END IF
98 :
99 :
100 14 : IF(.NOT.PRESENT(optDen)) THEN
101 2 : opt=.FALSE.
102 2 : CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
103 : ! read the charge density
104 : CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,&
105 2 : CDN_INPUT_DEN_const,0,fermiEnergyTemp,tempDistance,l_qfix,den)
106 : ELSE
107 12 : den=optDen
108 12 : opt=.TRUE.
109 : END IF
110 :
111 : ! flip cdn for each atom with rotation angles given
112 14 : if (any(noco%l_unrestrictMT).and.size(den%mt,4)<4) then
113 : !So far the density was collinear in spheres, now we make it FFN ready
114 2 : CALL move_alloc(den%mt,mt_tmp)
115 8 : allocate(den%mt(size(mt_tmp,1),0:size(mt_tmp,2)-1,size(mt_tmp,3),4))
116 491198 : den%mt(:,:,:,1:2)=mt_tmp
117 491198 : den%mt(:,:,:,3:)=0.0
118 :
119 2 : if (allocated(den%mmpMat)) then
120 2 : CALL move_alloc(den%mmpMat,mmpMat_tmp)
121 2 : allocate(den%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,SIZE(mmpMat_tmp,3),3))
122 234 : den%mmpMat(:,:,:,1:2)=mmpMat_tmp
123 116 : den%mmpMat(:,:,:,3)=0.0
124 : endif
125 : endif
126 : !$OMP parallel PRIVATE(rhodummy,rhodumms,j,rhodummyR,lh,itype,na) DEFAULT(none) &
127 : !$OMP SHARED(noco,den,zeros,atoms,sphhar,input,sym,l_flip,scalespin,toGlobal,nococonv) &
128 14 : !$OMP FIRSTPRIVATE(rotAngleTheta,rotAnglePhi)
129 : !$OMP do
130 : DO itype = 1, atoms%ntype
131 : na = atoms%firstAtom(itype)
132 : IF (l_flip(itype).AND.(.NOT.scaleSpin(itype))) THEN
133 : ! spherical and non-spherical m.t. charge density
134 : DO lh = 0,sphhar%nlh(sym%ntypsy(na))
135 : DO j = 1,atoms%jri(itype)
136 : IF (any(noco%l_unrestrictMT)) THEN
137 : rhodummy=CMPLX(den%mt(j,lh,itype,3),den%mt(j,lh,itype,4))
138 : !CALL rot_den_mat(zeros(itype),rotAngleTheta(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy)
139 : !CALL rot_den_mat(rotAnglePhi(itype),zeros(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy)
140 : call nococonv%rotdenmat(rotAnglePhi(itype),rotAngleTheta(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy, toGlobal)
141 : den%mt(j,lh,itype,3)=REAL(rhodummy)
142 : den%mt(j,lh,itype,4)=AIMAG(rhodummy)
143 : ELSE
144 : IF (rotAngleTheta(itype).EQ.(pimach()).AND.rotAnglePhi(itype).EQ.0) THEN
145 : rhodummyR = den%mt(j,lh,itype,1)
146 : den%mt(j,lh,itype,1) = den%mt(j,lh,itype,input%jspins)
147 : den%mt(j,lh,itype,input%jspins) = rhodummyR
148 : ELSE
149 : !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.
150 : CALL judft_error("l_mtNocoPot=F in combination with spin flips different from flipSpinTheta=Pi and flipSpinPhi=0 is currently not supported.",&
151 : calledby="flipcdn")
152 : END IF
153 : END IF
154 : END DO
155 : END DO
156 : ELSE IF (l_flip(itype).AND.scaleSpin(itype)) THEN
157 : 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")
158 : DO lh = 0,sphhar%nlh(sym%ntypsy(na))
159 : DO j = 1,atoms%jri(itype)
160 : rhodummy = den%mt(j,lh,itype,1) + den%mt(j,lh,itype,input%jspins)
161 : rhodumms = den%mt(j,lh,itype,1) - den%mt(j,lh,itype,input%jspins)
162 : den%mt(j,lh,itype,1) = 0.5 * (rhodummy + atoms%bmu(itype)*rhodumms)
163 : den%mt(j,lh,itype,input%jspins) = 0.5 * (rhodummy - atoms%bmu(itype)*rhodumms )
164 : END DO
165 : END DO
166 : END IF
167 : END DO
168 : !$OMP end do
169 : !$OMP end parallel
170 14 : IF (input%l_onlyMtStDen) THEN
171 :
172 : !!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.
173 0 : den%pw(:,2)=(den%pw(:,1)+den%pw(:,2))*0.5 !mean value
174 0 : den%pw(:,1)=den%pw(:,2)
175 0 : IF (noco%l_noco) THEN
176 0 : den%pw(:,3)=CMPLX(0.0,0.0)
177 : END IF
178 : END IF
179 : ! for LDA+U: flip density matrix
180 14 : if (allocated(den%mmpMat)) THEN
181 2450 : IF (input%lflip.AND.ANY(ABS(den%mmpMat) > 1e-12).AND.atoms%n_u+atoms%n_hia+atoms%n_opc>0) THEN
182 0 : DO i_u = 1, atoms%n_u+atoms%n_hia+atoms%n_opc
183 0 : if(i_u> atoms%n_u+atoms%n_hia) then
184 0 : itype = atoms%lda_u(i_u)%atomType
185 0 : l = atoms%lda_u(i_u)%l
186 : else
187 0 : itype = atoms%lda_opc(i_u-atoms%n_u-atoms%n_hia)%atomType
188 0 : l = atoms%lda_opc(i_u-atoms%n_u-atoms%n_hia)%l
189 : endif
190 0 : IF (l_flip(itype).AND.(.NOT.scaleSpin(itype))) THEN
191 0 : IF (any(noco%l_unrestrictMT)) THEN
192 : den%mmpMat(:,:,i_u,:) = rotMMPmat(den%mmpMat(:,:,i_u,:),rotAnglePhi(itype),rotAngleTheta(itype),0.0,&
193 0 : l,inverse=toGlobal,real_space_rotation=.FALSE., spin_rotation=.TRUE.)
194 : ELSE
195 0 : IF (rotAngleTheta(itype).EQ.(pimach()).AND.rotAnglePhi(itype).EQ.0) THEN
196 0 : DO m = -lmaxU_const,lmaxU_const
197 0 : DO mp = -lmaxU_const,lmaxU_const
198 0 : rhodummyR = den%mmpMat(m,mp,i_u,1)
199 0 : den%mmpMat(m,mp,i_u,1) = den%mmpMat(m,mp,i_u,input%jspins)
200 0 : den%mmpMat(m,mp,i_u,input%jspins) = rhodummyR
201 : ENDDO
202 : ENDDO
203 : ELSE
204 : !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.
205 : CALL judft_error("l_mtNocoPot=F in combination with spin flips different from flipSpinTheta=Pi, flipSpinPhi=0 is currently not supported.",&
206 0 : calledby="flipcdn")
207 : END IF
208 : END IF
209 0 : ELSE IF (l_flip(itype).AND.(scaleSpin(itype))) THEN
210 0 : DO m = -lmaxU_const,lmaxU_const
211 0 : DO mp = -lmaxU_const,lmaxU_const
212 0 : 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")
213 0 : rhodummy = den%mmpMat(m,mp,i_u,1) + den%mmpMat(m,mp,i_u,input%jspins)
214 0 : rhodumms = den%mmpMat(m,mp,i_u,1) - den%mmpMat(m,mp,i_u,input%jspins)
215 0 : den%mmpMat(m,mp,i_u,1) = 0.5 * (rhodummy + atoms%bmu(itype) * rhodumms)
216 0 : den%mmpMat(m,mp,i_u,input%jspins) = 0.5 * (rhodummy - atoms%bmu(itype) * rhodumms)
217 : END DO
218 : END DO
219 : END IF
220 : END DO
221 : END IF
222 : ENDIF
223 : ! write the spin-polarized density
224 14 : IF(input%lflip) CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,CDN_INPUT_DEN_const,&
225 2 : 1,-1.0,0.0,-1.0,-1.0,.FALSE.,den)
226 5894386 : IF(opt) optDen%mt=den%mt
227 : ! read enpara and flip lines
228 14 : INQUIRE(file='enpara',exist=n_exist)
229 14 : IF (n_exist) THEN
230 0 : OPEN(40,file ='enpara',status='old',form='formatted')
231 :
232 0 : j = 2
233 0 : DO itype = 1, atoms%ntype
234 0 : j = j + 1
235 0 : IF (atoms%nlo(itype)>0) j = j + 2
236 : END DO
237 0 : IF (input%film) j = j + 1
238 0 : ALLOCATE (clines(2*j))
239 0 : DO i = 1, 2*j
240 0 : READ (40,'(a)') clines(i)
241 : END DO
242 :
243 0 : REWIND 40
244 0 : i = 0
245 0 : DO ispin = 1,input%jspins
246 0 : i = i + 2
247 0 : WRITE (40,'(a)') TRIM(clines(i-1))
248 0 : WRITE (40,'(a)') TRIM(clines(i))
249 0 : DO itype = 1, atoms%ntype
250 0 : i = i + 1
251 0 : m = i
252 0 : IF (l_flip(itype)) m = MOD(i+j,2*j)
253 0 : IF (m==0) m = 2*j
254 0 : WRITE (40,'(a)') TRIM(clines(m))
255 0 : IF (atoms%nlo(itype)>0) THEN
256 0 : WRITE (40,'(a)') TRIM(clines(m+1))
257 0 : WRITE (40,'(a)') TRIM(clines(m+2))
258 0 : i = i + 2
259 : END IF
260 : END DO
261 0 : IF (input%film) THEN
262 0 : i = i + 1
263 0 : WRITE (40,'(a)') TRIM(clines(i))
264 : END IF
265 : END DO
266 0 : DEALLOCATE (clines)
267 0 : CLOSE(40)
268 : END IF
269 :
270 :
271 14 : IF(PRESENT(optDen)) optDen=den
272 :
273 14 : END SUBROUTINE flipcdn
274 :
275 : END MODULE m_flipcdn
|
