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 12 : 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 8 : 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 :
127 : !$OMP parallel PRIVATE(rhodummy,rhodumms,j,rhodummyR,lh,itype,na) DEFAULT(none) &
128 : !$OMP SHARED(noco,den,zeros,atoms,sphhar,input,sym,l_flip,scalespin,toGlobal,nococonv) &
129 14 : !$OMP FIRSTPRIVATE(rotAngleTheta,rotAnglePhi)
130 : !$OMP do
131 : DO itype = 1, atoms%ntype
132 : na = atoms%firstAtom(itype)
133 : IF (l_flip(itype).AND.(.NOT.scaleSpin(itype))) THEN
134 : ! spherical and non-spherical m.t. charge density
135 : DO lh = 0,sphhar%nlh(sym%ntypsy(na))
136 : DO j = 1,atoms%jri(itype)
137 : IF (any(noco%l_unrestrictMT)) THEN
138 : rhodummy=CMPLX(den%mt(j,lh,itype,3),den%mt(j,lh,itype,4))
139 : !CALL rot_den_mat(zeros(itype),rotAngleTheta(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy)
140 : !CALL rot_den_mat(rotAnglePhi(itype),zeros(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy)
141 : call nococonv%rotdenmat(rotAnglePhi(itype),rotAngleTheta(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy, toGlobal)
142 : den%mt(j,lh,itype,3)=REAL(rhodummy)
143 : den%mt(j,lh,itype,4)=AIMAG(rhodummy)
144 : ELSE
145 : IF (rotAngleTheta(itype).EQ.(pimach()).AND.rotAnglePhi(itype).EQ.0) THEN
146 : rhodummyR = den%mt(j,lh,itype,1)
147 : den%mt(j,lh,itype,1) = den%mt(j,lh,itype,input%jspins)
148 : den%mt(j,lh,itype,input%jspins) = rhodummyR
149 : ELSE
150 : !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.
151 : CALL judft_error("l_mtNocoPot=F in combination with spin flips different from flipSpinTheta=Pi and flipSpinPhi=0 is currently not supported.",&
152 : calledby="flipcdn")
153 : END IF
154 : END IF
155 : END DO
156 : END DO
157 : ELSE IF (l_flip(itype).AND.scaleSpin(itype)) THEN
158 : 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")
159 : DO lh = 0,sphhar%nlh(sym%ntypsy(na))
160 : DO j = 1,atoms%jri(itype)
161 : rhodummy = den%mt(j,lh,itype,1) + den%mt(j,lh,itype,input%jspins)
162 : rhodumms = den%mt(j,lh,itype,1) - den%mt(j,lh,itype,input%jspins)
163 : den%mt(j,lh,itype,1) = 0.5 * (rhodummy + atoms%bmu(itype)*rhodumms)
164 : den%mt(j,lh,itype,input%jspins) = 0.5 * (rhodummy - atoms%bmu(itype)*rhodumms )
165 : END DO
166 : END DO
167 : END IF
168 : END DO
169 : !$OMP end do
170 : !$OMP end parallel
171 :
172 14 : IF (input%l_onlyMtStDen) THEN
173 : !!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.
174 0 : den%pw(:,2)=(den%pw(:,1)+den%pw(:,2))*0.5 !mean value
175 0 : den%pw(:,1)=den%pw(:,2)
176 0 : IF (noco%l_noco) THEN
177 0 : den%pw(:,3)=CMPLX(0.0,0.0)
178 : END IF
179 : END IF
180 :
181 :
182 : ! for LDA+U: flip density matrix
183 2450 : IF (input%lflip.AND.ANY(ABS(den%mmpMat) > 1e-12).AND.atoms%n_u+atoms%n_hia+atoms%n_opc>0) THEN
184 0 : DO i_u = 1, atoms%n_u+atoms%n_hia+atoms%n_opc
185 0 : if(i_u> atoms%n_u+atoms%n_hia) then
186 0 : itype = atoms%lda_u(i_u)%atomType
187 0 : l = atoms%lda_u(i_u)%l
188 : else
189 0 : itype = atoms%lda_opc(i_u-atoms%n_u-atoms%n_hia)%atomType
190 0 : l = atoms%lda_opc(i_u-atoms%n_u-atoms%n_hia)%l
191 : endif
192 0 : IF (l_flip(itype).AND.(.NOT.scaleSpin(itype))) THEN
193 0 : IF (any(noco%l_unrestrictMT)) THEN
194 : den%mmpMat(:,:,i_u,:) = rotMMPmat(den%mmpMat(:,:,i_u,:),rotAnglePhi(itype),rotAngleTheta(itype),0.0,&
195 0 : l,inverse=toGlobal,real_space_rotation=.FALSE., spin_rotation=.TRUE.)
196 : ELSE
197 0 : IF (rotAngleTheta(itype).EQ.(pimach()).AND.rotAnglePhi(itype).EQ.0) THEN
198 0 : DO m = -lmaxU_const,lmaxU_const
199 0 : DO mp = -lmaxU_const,lmaxU_const
200 0 : rhodummyR = den%mmpMat(m,mp,i_u,1)
201 0 : den%mmpMat(m,mp,i_u,1) = den%mmpMat(m,mp,i_u,input%jspins)
202 0 : den%mmpMat(m,mp,i_u,input%jspins) = rhodummyR
203 : ENDDO
204 : ENDDO
205 : ELSE
206 : !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.
207 : CALL judft_error("l_mtNocoPot=F in combination with spin flips different from flipSpinTheta=Pi, flipSpinPhi=0 is currently not supported.",&
208 0 : calledby="flipcdn")
209 : END IF
210 : END IF
211 0 : ELSE IF (l_flip(itype).AND.(scaleSpin(itype))) THEN
212 0 : DO m = -lmaxU_const,lmaxU_const
213 0 : DO mp = -lmaxU_const,lmaxU_const
214 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")
215 0 : rhodummy = den%mmpMat(m,mp,i_u,1) + den%mmpMat(m,mp,i_u,input%jspins)
216 0 : rhodumms = den%mmpMat(m,mp,i_u,1) - den%mmpMat(m,mp,i_u,input%jspins)
217 0 : den%mmpMat(m,mp,i_u,1) = 0.5 * (rhodummy + atoms%bmu(itype) * rhodumms)
218 0 : den%mmpMat(m,mp,i_u,input%jspins) = 0.5 * (rhodummy - atoms%bmu(itype) * rhodumms)
219 : END DO
220 : END DO
221 : END IF
222 : END DO
223 : END IF
224 :
225 : ! write the spin-polarized density
226 14 : IF(input%lflip) CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,CDN_INPUT_DEN_const,&
227 2 : 1,-1.0,0.0,-1.0,-1.0,.FALSE.,den)
228 5894386 : IF(opt) optDen%mt=den%mt
229 :
230 : ! read enpara and flip lines
231 14 : INQUIRE(file='enpara',exist=n_exist)
232 14 : IF (n_exist) THEN
233 0 : OPEN(40,file ='enpara',status='old',form='formatted')
234 :
235 0 : j = 2
236 0 : DO itype = 1, atoms%ntype
237 0 : j = j + 1
238 0 : IF (atoms%nlo(itype)>0) j = j + 2
239 : END DO
240 0 : IF (input%film) j = j + 1
241 0 : ALLOCATE (clines(2*j))
242 0 : DO i = 1, 2*j
243 0 : READ (40,'(a)') clines(i)
244 : END DO
245 :
246 0 : REWIND 40
247 0 : i = 0
248 0 : DO ispin = 1,input%jspins
249 0 : i = i + 2
250 0 : WRITE (40,'(a)') TRIM(clines(i-1))
251 0 : WRITE (40,'(a)') TRIM(clines(i))
252 0 : DO itype = 1, atoms%ntype
253 0 : i = i + 1
254 0 : m = i
255 0 : IF (l_flip(itype)) m = MOD(i+j,2*j)
256 0 : IF (m==0) m = 2*j
257 0 : WRITE (40,'(a)') TRIM(clines(m))
258 0 : IF (atoms%nlo(itype)>0) THEN
259 0 : WRITE (40,'(a)') TRIM(clines(m+1))
260 0 : WRITE (40,'(a)') TRIM(clines(m+2))
261 0 : i = i + 2
262 : END IF
263 : END DO
264 0 : IF (input%film) THEN
265 0 : i = i + 1
266 0 : WRITE (40,'(a)') TRIM(clines(i))
267 : END IF
268 : END DO
269 0 : DEALLOCATE (clines)
270 0 : CLOSE(40)
271 : END IF
272 :
273 :
274 14 : IF(PRESENT(optDen)) optDen=den
275 :
276 14 : END SUBROUTINE flipcdn
277 :
278 : END MODULE m_flipcdn
|
