Line data Source code
1 : MODULE m_greensfTorque
2 :
3 : USE m_types
4 : USE m_juDFT
5 : USE m_constants
6 : USE m_intgr
7 : USE m_gaunt
8 : USE m_xmlOutput
9 : USE m_lattHarmsSphHarmsConv
10 : #ifdef CPP_MPI
11 : USE mpi
12 : #endif
13 :
14 : IMPLICIT NONE
15 :
16 : CONTAINS
17 :
18 0 : SUBROUTINE greensfTorque(greensFunction,gfinp,fmpi,sphhar,atoms,sym,noco,nococonv,input,f,g,flo,vTot)
19 :
20 : !--------------------------------------------------------------------------
21 : ! This Subroutine implements the formula:
22 : ! alpha 1 -> Ef i -> alpha ->->
23 : ! J = - - Im Tr int dx int dE B (x) sigma G (x,x,E)
24 : ! i pi sigma -infinity xc ii
25 : !
26 : ! For the evaluation of the torque at site i
27 : !--------------------------------------------------------------------------
28 :
29 : TYPE(t_greensf), INTENT(IN) :: greensFunction(:)
30 : TYPE(t_gfinp), INTENT(IN) :: gfinp
31 : TYPE(t_mpi), INTENT(IN) :: fmpi
32 : TYPE(t_sphhar), INTENT(IN) :: sphhar
33 : TYPE(t_atoms), INTENT(IN) :: atoms
34 : TYPE(t_sym), INTENT(IN) :: sym
35 : TYPE(t_noco), INTENT(IN) :: noco
36 : TYPE(t_nococonv), INTENT(IN) :: nococonv
37 : TYPE(t_input), INTENT(IN) :: input
38 : REAL, INTENT(IN) :: f(:,:,0:,:,:)
39 : REAL, INTENT(IN) :: g(:,:,0:,:,:)
40 : REAL, INTENT(IN) :: flo(:,:,:,:,:)
41 : TYPE(t_potden), INTENT(IN) :: vTot
42 :
43 : INTEGER :: l,lp,iContour,iGrid,ispin,iTorque,atomType,index_task,extra,ierr
44 : INTEGER :: lh,mu,m,mp,iz,ipm,jr,alpha,lhmu,index,index_start,index_end,n,i_gf
45 : COMPLEX :: phaseFactor, weight
46 : REAL :: realIntegral
47 : CHARACTER(LEN=20) :: attributes(5)
48 :
49 0 : TYPE(t_greensf) :: gf_rot
50 :
51 0 : REAL, ALLOCATABLE :: torque(:,:),rtmp(:)
52 0 : COMPLEX, ALLOCATABLE :: bxc(:,:,:)
53 0 : COMPLEX, ALLOCATABLE :: integrand(:,:)
54 0 : COMPLEX, ALLOCATABLE :: g_ii(:,:), mag_ii(:,:)
55 0 : COMPLEX, ALLOCATABLE :: vlm(:,:,:)
56 0 : INTEGER, ALLOCATABLE :: gf_indices(:)
57 :
58 0 : ALLOCATE(gf_indices(SUM(gfinp%numTorqueElems(:))), source=-1)
59 0 : ALLOCATE(bxc(atoms%jmtd,atoms%lmaxd*(atoms%lmaxd+2)+1,atoms%ntype), source=cmplx_0)
60 0 : CALL timestart("Green's Function Torque: init")
61 :
62 0 : DO atomType = 1, atoms%ntype
63 0 : IF(gfinp%numTorqueElems(atomType)==0) CYCLE
64 :
65 : iContour = -1
66 0 : DO iTorque = 1, gfinp%numTorqueElems(atomType)
67 0 : gf_indices(SUM(gfinp%numTorqueElems(:atomType-1))+iTorque) = gfinp%torqueElem(atomType,iTorque)
68 :
69 : !Check that its actually right
70 0 : IF(greensFunction(gfinp%torqueElem(atomType,iTorque))%elem%atomType.NE.atomType.OR.&
71 : greensFunction(gfinp%torqueElem(atomType,iTorque))%elem%atomTypep.NE.atomType) THEN
72 0 : CALL juDFT_error("Provided greensFunction for wrong atomType", calledby="greensFunctionTorque")
73 : ENDIF
74 :
75 0 : IF(iContour == -1) THEN
76 0 : iContour = greensFunction(gfinp%torqueElem(atomType,iTorque))%elem%iContour
77 0 : ELSE IF(greensFunction(gfinp%torqueElem(atomType,iTorque))%elem%iContour/=iContour) THEN
78 0 : CALL juDFT_error("Provided different energy contours", calledby="greensFunctionTorque")
79 : ENDIF
80 : ENDDO
81 : !Get Bxc from the total potential (local frame)
82 : !TODO: FFN components
83 0 : ALLOCATE(vlm(atoms%jmtd,atoms%lmaxd*(atoms%lmaxd+2)+1,input%jspins),source=cmplx_0)
84 0 : vlm = cmplx_0
85 0 : DO ispin = 1, input%jspins
86 0 : CALL lattHarmsRepToSphHarms(sym, atoms, sphhar, atomType, vTot%mt(:,0:,atomType,ispin), vlm(:,:,ispin))
87 : ENDDO
88 : !Get the Bxc part of the potential
89 0 : bxc(:,:,atomType) = (vlm(:,:,1) - vlm(:,:,2))/2.0
90 0 : DEALLOCATE(vlm)
91 :
92 : !L=0 of potential has an additional rescaling of r/sqrt(4pi)
93 : bxc(:atoms%jri(atomType),1,atomType) = bxc(:atoms%jri(atomType),1,atomType) &
94 0 : * sfp_const/atoms%rmsh(:atoms%jri(atomType),atomType)
95 : ENDDO
96 :
97 0 : IF(ANY(gf_indices<1) .OR. ANY(gf_indices>SIZE(greensFunction))) THEN
98 0 : CALL juDFT_error("Invalid index in greensFunction mapping array", calledby="greensFunctionTorque")
99 : ENDIF
100 :
101 : #ifdef CPP_MPI
102 0 : IF(fmpi%isize > 1) THEN
103 : !Just distribute the individual gf elements over the ranks
104 0 : index_task = FLOOR(REAL(SIZE(gf_indices))/(fmpi%isize))
105 0 : extra = SIZE(gf_indices) - index_task*fmpi%isize
106 0 : index_start = fmpi%irank*index_task + 1 + extra
107 0 : index_end = (fmpi%irank+1)*index_task + extra
108 0 : IF(fmpi%irank < extra) THEN
109 0 : index_start = index_start - (extra - fmpi%irank)
110 0 : index_end = index_end - (extra - fmpi%irank - 1)
111 : ENDIF
112 : ELSE
113 0 : index_start = 1
114 0 : index_end = SIZE(gf_indices)
115 : ENDIF
116 : #else
117 : index_start = 1
118 : index_end = SIZE(gf_indices)
119 : #endif
120 :
121 :
122 0 : CALL timestop("Green's Function Torque: init")
123 0 : CALL timestart("Green's Function Torque: Integration")
124 :
125 0 : ALLOCATE(torque(3,atoms%ntype), source=0.0)
126 0 : !$ phaseFactor = gaunt1(0,0,0,0,0,0,atoms%lmaxd)
127 0 : DO index = index_start, index_end
128 0 : IF(index.LT.1 .OR. index.GT.SIZE(gf_indices)) CYCLE
129 0 : i_gf = gf_indices(index)
130 :
131 0 : gf_rot = greensFunction(i_gf)
132 0 : l = gf_rot%elem%l
133 0 : lp = gf_rot%elem%lp
134 0 : atomType = gf_rot%elem%atomType
135 :
136 : !Rotate the greens function into the global real space frame
137 0 : IF(noco%l_noco) THEN
138 0 : CALL gf_rot%rotate_euler_angles(atoms,nococonv%alph(atomType),nococonv%beta(atomType),0.0)
139 0 : ELSE IF(noco%l_soc) THEN
140 0 : CALL gf_rot%rotate_euler_angles(atoms,nococonv%phi,nococonv%theta,0.0)
141 : ENDIF
142 :
143 : #ifndef CPP_NOTYPEPROCINOMP
144 : !$OMP parallel default(none) &
145 : !$OMP shared(sphhar,atoms,input,gf_rot,f,g,flo,bxc) &
146 : !$OMP shared(l,lp,atomType,torque) &
147 : !$OMP private(lh,m,mu,mp,lhmu,phaseFactor,weight,ispin,ipm,iz,alpha,jr) &
148 0 : !$OMP private(realIntegral,integrand,g_ii,mag_ii)
149 : #endif
150 : ALLOCATE(integrand(atoms%jmtd,3),source=cmplx_0)
151 : ALLOCATE(g_ii(atoms%jmtd,gf_rot%contour%nz),source=cmplx_0)
152 : ALLOCATE(mag_ii(atoms%jmtd,gf_rot%contour%nz),source=cmplx_0)
153 : #ifndef CPP_NOTYPEPROCINOMP
154 : !$OMP do collapse(2)
155 : #endif
156 : DO lh = 0, atoms%lmaxd
157 : DO m = -l, l
158 : IF(MOD(lh+l+lp,2) .NE. 0) CYCLE
159 : IF(lh.GT.l+lp) CYCLE
160 : IF(lh.LT.abs(l-lp)) CYCLE
161 : DO mu = -lh, lh
162 : lhmu = lh * (lh+1) + mu + 1
163 : mp = m + mu
164 : IF(ABS(mp).GT.lp) CYCLE
165 : phaseFactor = gaunt1(lp,lh,l,mp,mu,m,atoms%lmaxd)
166 : IF(ABS(phaseFactor).LT.1e-12) CYCLE
167 : integrand = cmplx_0
168 : DO ipm = 1, 2
169 : DO alpha = 1, 3 !(x,y,z)
170 : IF (alpha.EQ.1) THEN
171 : !magnetization in x-direction
172 : CALL gf_rot%getRadial(atoms,m,mp,ipm==2,3,f,g,flo,mag_ii)
173 : CALL gf_rot%getRadial(atoms,mp,m,ipm==2,3,f,g,flo,g_ii)
174 : mag_ii = mag_ii + conjg(g_ii)
175 : ELSE IF (alpha.EQ.2) THEN
176 : !magnetization in y-direction
177 : CALL gf_rot%getRadial(atoms,m,mp,ipm==2,3,f,g,flo,mag_ii)
178 : CALL gf_rot%getRadial(atoms,mp,m,ipm==2,3,f,g,flo,g_ii)
179 : mag_ii = ImagUnit * (mag_ii - conjg(g_ii))
180 : ELSE
181 : !magnetization in z-direction
182 : mag_ii = cmplx_0
183 : DO ispin = 1, input%jspins
184 : CALL gf_rot%getRadial(atoms,m,mp,ipm==2,ispin,f,g,flo,g_ii)
185 : mag_ii = mag_ii + (-1)**(ispin-1) * g_ii
186 : ENDDO
187 : ENDIF
188 :
189 : DO iz = 1, SIZE(mag_ii,2)
190 : weight = gf_rot%contour%de(iz) * phaseFactor
191 : weight = MERGE(weight, conjg(weight), ipm==1)
192 :
193 : DO jr = 1, atoms%jri(atomType)
194 : integrand(jr,alpha) = integrand(jr,alpha) + ImagUnit/tpi_const * (-1)**(ipm-1) * mag_ii(jr,iz) &
195 : * MERGE(bxc(jr,lhmu,atomType),conjg(bxc(jr,lhmu,atomType)), ipm==1) * weight
196 : ENDDO
197 : ENDDO
198 : ENDDO
199 : ENDDO
200 :
201 : DO alpha = 1, 3 !(x,y,z)
202 : CALL intgr3(REAL(integrand(:,alpha)),atoms%rmsh(:,atomType),atoms%dx(atomType),atoms%jri(atomType),realIntegral)
203 : #ifndef CPP_NOTYPEPROCINOMP
204 : !$OMP critical
205 : torque(alpha,atomType) = torque(alpha,atomType) + realIntegral
206 : !$OMP end critical
207 : #else
208 : torque(alpha,atomType) = torque(alpha,atomType) + realIntegral
209 : #endif
210 : ENDDO
211 : ENDDO
212 : ENDDO
213 : ENDDO
214 : #ifndef CPP_NOTYPEPROCINOMP
215 : !$OMP end do
216 : DEALLOCATE(integrand,g_ii,mag_ii)
217 : !$OMP end parallel
218 : #else
219 : DEALLOCATE(integrand,g_ii,mag_ii)
220 : #endif
221 :
222 : ENDDO
223 0 : CALL timestop("Green's Function Torque: Integration")
224 :
225 : #ifdef CPP_MPI
226 : !Collect the torque to rank 0
227 0 : n = SIZE(torque)
228 0 : ALLOCATE(rtmp(n))
229 0 : CALL MPI_REDUCE(torque,rtmp,n,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
230 0 : IF(fmpi%irank.EQ.0) CALL dcopy(n,rtmp,1,torque,1)
231 0 : DEALLOCATE(rtmp)
232 : #endif
233 :
234 0 : IF(fmpi%irank.EQ.0) THEN
235 0 : CALL openXMLElementNoAttributes('noncollinearTorque')
236 0 : WRITE(oUnit,'(/,A)') 'Torque Calculation (noco):'
237 0 : WRITE(oUnit,'(/,A)') '---------------------------'
238 0 : DO atomType = 1, atoms%ntype
239 0 : IF(gfinp%numTorqueElems(atomType)==0) CYCLE
240 0 : WRITE(oUnit,'(A,I4,A,3f16.8,A)') ' atom: ', atomType, ' torque: ', torque(:,atomType) * hartree_to_ev_const * 1000, ' meV'
241 :
242 0 : attributes = ''
243 0 : WRITE(attributes(1),'(i0)') atomType
244 0 : WRITE(attributes(2),'(f14.8)') torque(1,atomType) * hartree_to_ev_const * 1000
245 0 : WRITE(attributes(3),'(f14.8)') torque(2,atomType) * hartree_to_ev_const * 1000
246 0 : WRITE(attributes(4),'(f14.8)') torque(3,atomType) * hartree_to_ev_const * 1000
247 0 : WRITE(attributes(5),'(a3)') 'meV'
248 : CALL writeXMLElementForm('torque',['atomType','sigma_x ','sigma_y ','sigma_z ','units '],&
249 0 : attributes,reshape([8,7,7,7,4,6,14,14,14,3],[5,2]))
250 : ENDDO
251 0 : CALL closeXMLElement('noncollinearTorque')
252 : ENDIF
253 :
254 0 : END SUBROUTINE greensfTorque
255 :
256 : END MODULE m_greensfTorque
|
