Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2019 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 : MODULE m_xcBfield
7 : USE m_types
8 : USE m_constants
9 : USE m_plot
10 : USE m_divergence
11 : USE m_juDFT
12 :
13 : IMPLICIT NONE
14 :
15 : !-----------------------------------------------------------------------------
16 : ! This module contains all the operations on exchange-correlation B-fields
17 : ! that are necessary to project out source terms. This way, the whole modifi-
18 : ! cation towards source-free fields can be done by one call, either as a post-
19 : ! process test or in the scf-loop to achieve said fields self-consistently.
20 : !-----------------------------------------------------------------------------
21 :
22 : PUBLIC :: makeVectorField, sourcefree, correctPot
23 :
24 : CONTAINS
25 10 : SUBROUTINE makeVectorField(sym,stars,atoms,sphhar,vacuum,input,noco,nococonv,denmat,factor,vScal,aVec,cell)
26 :
27 : ! Contructs the exchange-correlation magnetic field from the total poten-
28 : ! tial matrix or the magnetic density for the density matrix. Only used for
29 : ! the implementation of source free fields and therefore only applicable in
30 : ! a (fully) non-collinear description of magnetism.
31 : !
32 : ! Assumes the following form for the density/potential matrix:
33 : ! rho_mat = n*Id_(2x2) + conj(sigma_vec)*m_vec
34 : ! V_mat = V*Id_(2x2) + sigma_vec*B_vec
35 : !
36 : ! A_vec is saved as a density type with an additional r^2-factor.
37 : !
38 : ! TODO: How do we constuct B when not only it is saved to the density
39 : ! matrix (SOC, LDA+U etc.)?
40 : TYPE(t_sym), INTENT(IN) :: sym
41 : TYPE(t_stars), INTENT(IN) :: stars
42 : TYPE(t_atoms), INTENT(IN) :: atoms
43 : TYPE(t_sphhar), INTENT(IN) :: sphhar
44 : TYPE(t_vacuum), INTENT(IN) :: vacuum
45 : TYPE(t_input), INTENT(IN) :: input
46 : TYPE(t_noco), INTENT(IN) :: noco
47 : TYPE(t_nococonv), INTENT(IN) :: nococonv
48 : TYPE(t_potden), INTENT(IN) :: denmat
49 : REAL, INTENT(IN) :: factor
50 : TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: aVec
51 : TYPE(t_potden), INTENT(OUT) :: vScal
52 : TYPE(t_cell), INTENT(IN) :: cell
53 :
54 : TYPE(t_gradients) :: tmp_grad
55 140 : TYPE(t_potden), DIMENSION(4) :: dummyDen
56 : INTEGER :: i, itype, ir, lh
57 2 : REAL :: r2(atoms%jmtd), fcut
58 : REAL, ALLOCATABLE :: intden(:,:)
59 :
60 2 : fcut=1.e-12
61 :
62 : ! Initialize and fill a dummy density array, that takes the initial result
63 : ! of matrixsplit.
64 10 : DO i=1, 4
65 : CALL dummyDen(i)%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd, &
66 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE., &
67 8 : POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
68 34 : ALLOCATE(dummyDen(i)%pw_w,mold=dummyDen(i)%pw)
69 : ENDDO
70 :
71 : CALL matrixsplit(sym,stars,atoms,sphhar,vacuum,input,noco,nococonv,factor,denmat, &
72 2 : dummyDen(1),dummyDen(2),dummyDen(3),dummyDen(4))
73 :
74 2 : vScal=dummyDen(1)
75 :
76 1516 : r2=1.0
77 :
78 : ! Initialize and fill the vector field.
79 8 : DO i=1,3
80 : CALL aVec(i)%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd, &
81 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE., &
82 6 : POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
83 24 : ALLOCATE(aVec(i)%pw_w,mold=aVec(i)%pw)
84 736800 : aVec(i)%mt(:,:,:,:) = dummyDen(i+1)%mt(:,:,:,:)
85 18 : DO itype=1,atoms%ntype
86 990 : DO lh=0, sphhar%nlh(sym%ntypsy(atoms%firstAtom(itype)))
87 972 : IF (factor==2.0) THEN
88 736776 : r2=atoms%rmsh(:,itype)**2
89 : END IF
90 972 : IF (i==1) THEN
91 245592 : vScal%mt(:,lh,itype,1) = vScal%mt(:,lh,itype,1)*r2
92 : END IF
93 736788 : aVec(i)%mt(:,lh,itype,1) = aVec(i)%mt(:,lh,itype,1)*r2
94 : END DO !lh
95 : END DO !itype
96 18438 : aVec(i)%pw(1:,:) = dummyDen(i+1)%pw(1:,:)
97 18438 : aVec(i)%pw_w(1:,:) = dummyDen(i+1)%pw_w(1:,:)
98 : !aVec(i)%vacz(1:,1:,:) = dummyDen(i+1)%vacz(1:,1:,:)
99 : !aVec(i)%vacxy(1:,1:,1:,:) = dummyDen(i+1)%vacxy(1:,1:,1:,:)
100 26 : aVec(i)%vac(1:,1:,1:,:) = dummyDen(i+1)%vac(1:,1:,1:,:)
101 : END DO !i
102 :
103 128 : END SUBROUTINE makeVectorField
104 :
105 2 : SUBROUTINE sourcefree(fmpi,field,stars,atoms,sphhar,vacuum,input ,sym,juphon,cell,noco,aVec,vScal,vCorr)
106 : USE m_vgen_coulomb
107 : USE m_gradYlm
108 : USE m_grdchlh
109 : USE m_sphpts
110 : USE m_checkdop
111 : USE m_BfieldtoVmat
112 : USE m_pw_tofrom_grid
113 :
114 : ! Takes a vectorial quantity, i.e. a t_potden variable of dimension 3, and
115 : ! makes it into a source free vector field as follows:
116 : !
117 : ! a) Build the divergence d of the vector field A_vec as d=nabla_vec*A_vec.
118 : ! b) Solve the Poisson equation (nabla_vec*nabla_vec)phi=-4*pi*d for phi.
119 : ! c) Construct an auxiliary vector field C_vec=(nabla_vec phi)/(4*pi).
120 : ! d) Build A_vec_sf=A_vec+C_vec, which is source free by construction.
121 : !
122 : ! Note: A_vec is assumed to be a density with an additional factor of r^2.
123 : ! A field of the same form will also be calculated.
124 :
125 : TYPE(t_mpi), INTENT(IN) :: fmpi
126 : TYPE(t_field), INTENT(IN) :: field
127 : TYPE(t_stars), INTENT(IN) :: stars
128 : TYPE(t_atoms), INTENT(IN) :: atoms
129 : TYPE(t_sphhar), INTENT(IN) :: sphhar
130 : TYPE(t_vacuum), INTENT(IN) :: vacuum
131 : TYPE(t_input), INTENT(IN) :: input
132 :
133 : TYPE(t_sym), INTENT(IN) :: sym
134 : TYPE(t_juphon), INTENT(IN) :: juphon
135 : TYPE(t_cell), INTENT(IN) :: cell
136 : TYPE(t_noco), INTENT(IN) :: noco
137 : TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: aVec
138 : TYPE(t_potden), INTENT(IN) :: vScal
139 : TYPE(t_potden), INTENT(OUT) :: vCorr
140 :
141 2 : TYPE(t_potden) :: div, phi, divloc
142 212 : TYPE(t_potden), DIMENSION(3) :: cvec, corrB
143 2 : TYPE(t_atoms) :: atloc
144 : INTEGER :: n, jr, lh, lhmax, jcut, nat ,i
145 : REAL :: xp(3,(atoms%lmaxd+1+mod(atoms%lmaxd+1,2))*(2*atoms%lmaxd+1))
146 2 : REAL, ALLOCATABLE :: intden(:,:)
147 2 : TYPE(t_gradients) :: tmp_grad
148 : complex :: sigma_loc(2)
149 :
150 : CALL div%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,atoms%ntype, &
151 : atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN, &
152 2 : vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
153 8 : ALLOCATE(div%pw_w,mold=div%pw)
154 6146 : div%pw_w = CMPLX(0.0,0.0)
155 :
156 2 : CALL timestart("Building divergence")
157 2 : CALL divergence(input,stars,atoms,sphhar,vacuum,sym,cell,noco,aVec,div)
158 2 : CALL timestop("Building divergence")
159 :
160 : ! Local atoms variable with no charges;
161 : ! needed for the potential generation from the divergence.
162 6 : atloc=atoms
163 6 : atloc%zatom=0.0
164 2 : CALL phi%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTCOUL,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
165 8 : ALLOCATE(phi%pw_w(SIZE(phi%pw,1),size(phi%pw,2)))
166 6146 : phi%pw_w = CMPLX(0.0,0.0)
167 :
168 2 : CALL timestart("Building potential")
169 2 : sigma_loc = cmplx(0.0,0.0)
170 2 : CALL vgen_coulomb(1,fmpi ,input,field,vacuum,sym,juphon,stars,cell,sphhar,atloc,.TRUE.,div,phi,sigma_loc)
171 2 : CALL timestop("Building potential")
172 :
173 8 : DO i=1,3
174 6 : CALL cvec(i)%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
175 24 : ALLOCATE(cvec(i)%pw_w,mold=cvec(i)%pw)
176 18440 : cvec(i)%pw_w=CMPLX(0.0,0.0)
177 : ENDDO
178 :
179 2 : CALL timestart("Building correction field")
180 2 : CALL divpotgrad(input,stars,atloc,sphhar,vacuum,sym,cell,noco,phi,cvec)
181 2 : CALL timestop("Building correction field")
182 :
183 2 : CALL init_pw_grid(stars,sym,cell)
184 8 : DO i=1,3
185 6 : CALL pw_to_grid(.FALSE.,1,.FALSE.,stars,cell,cvec(i)%pw,tmp_grad,rho=intden)
186 18438 : cvec(i)%pw=CMPLX(0.0,0.0)
187 18438 : cvec(i)%pw_w=CMPLX(0.0,0.0)
188 6 : CALL pw_from_grid(stars,intden,cvec(i)%pw,cvec(i)%pw_w)
189 8 : DEALLOCATE(intden)
190 : END DO !i
191 2 : CALL finish_pw_grid()
192 :
193 8 : DO i=1,3
194 : CALL corrB(i)%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
195 6 : POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
196 24 : ALLOCATE(corrB(i)%pw_w,mold=corrB(i)%pw)
197 18438 : corrB(i)%pw_w=CMPLX(0.0,0.0)
198 8 : CALL corrB(i)%addPotDen(aVec(i),cvec(i))
199 : ENDDO
200 :
201 2 : CALL BfieldtoVmat(sym, stars, atoms, sphhar, vacuum, vScal, corrB(1), corrB(2), corrB(3), vCorr)
202 :
203 8 : END SUBROUTINE sourcefree
204 :
205 0 : SUBROUTINE correctPot(vTot,c)
206 : USE m_types
207 :
208 : ! Takes a vectorial quantity c and saves its components into the appro
209 : ! priate components of the potential matrix V.
210 : !
211 : ! An initial V_mat = V*Id_(2x2) + sigma_vec*B_vec will become
212 : ! V_mat = V*Id_(2x2) + sigma_vec*(B_vec+C_vec)
213 : !
214 : ! TODO: Both quantities are assumed to be in the global frame of refe-
215 : ! rence. Make sure this is true. Also: consider SOC, LDA+U etc.
216 :
217 : TYPE(t_potden), INTENT(INOUT) :: vTot
218 : TYPE(t_potden), DIMENSION(3), INTENT(IN) :: c
219 :
220 0 : REAL :: pwr(SIZE(vTot%pw(1:,3))), pwi(SIZE(vTot%pw(1:,3)))
221 :
222 : !vtot%mt(:,0:,:,1)=(vtot%mt(:,0:,:,1)+vtot%mt(:,0:,:,2))/2
223 : !vtot%mt(:,0:,:,2)=vtot%mt(:,0:,:,1)
224 : !vtot%mt(:,0:,:,3:4)=0.0
225 :
226 : !vTot%pw(1:,1)=(vTot%pw(1:,1)+vTot%pw(1:,2))/2
227 : !vTot%pw(1:,2)=vTot%pw(1:,1)
228 : !vTot%pw(1:,3)=CMPLX(0.0,0.0)
229 :
230 : !vTot%pw_w(1:,1)=(vTot%pw_w(1:,1)+vTot%pw_w(1:,2))/2
231 : !vTot%pw_w(1:,2)=vTot%pw_w(1:,1)
232 : !vTot%pw_w(1:,3)=CMPLX(0.0,0.0)
233 :
234 : !vTot%mt(:,0:,:,1)=vTot%mt(:,0:,:,1)+b(3)%mt(:,0:,:,1)/2
235 : !vTot%mt(:,0:,:,2)=vTot%mt(:,0:,:,2)-b(3)%mt(:,0:,:,1)/2
236 : !vTot%mt(:,0:,:,3)=b(1)%mt(:,0:,:,1)/2
237 : !vTot%mt(:,0:,:,4)=b(2)%mt(:,0:,:,1)/2
238 :
239 : !vTot%pw(1:,1)=vTot%pw(1:,1)+b(3)%pw(1:,1)/2
240 : !vTot%pw(1:,2)=vTot%pw(1:,2)-b(3)%pw(1:,1)/2
241 : !vTot%pw(1:,3)=(b(1)%pw(1:,1)+ImagUnit*b(2)%pw(1:,1))/2
242 :
243 : !vTot%pw_w(1:,1)=vTot%pw_w(1:,1)+b(3)%pw_w(1:,1)/2
244 : !vTot%pw_w(1:,2)=vTot%pw_w(1:,2)-b(3)%pw_w(1:,1)/2
245 : !vTot%pw_w(1:,3)=(b(1)%pw_w(1:,1)+ImagUnit*b(2)%pw_w(1:,1))/2
246 :
247 0 : vTot%mt(:,0:,:,1)=c(3)%mt(:,0:,:,1)+vTot%mt(:,0:,:,1)
248 0 : vTot%mt(:,0:,:,2)=-c(3)%mt(:,0:,:,1)+vTot%mt(:,0:,:,2)
249 0 : vTot%mt(:,0:,:,3)=c(1)%mt(:,0:,:,1)+vTot%mt(:,0:,:,3)
250 0 : vTot%mt(:,0:,:,4)=c(2)%mt(:,0:,:,1)+vTot%mt(:,0:,:,4)
251 :
252 0 : vTot%pw(1:,1)=vTot%pw(1:,1)+c(3)%pw(1:,1)
253 0 : vTot%pw(1:,2)=vTot%pw(1:,2)-c(3)%pw(1:,1)
254 0 : pwr=REAL(vTot%pw(1:,3)) + REAL(c(1)%pw(1:,1)) - AIMAG(c(2)%pw(1:,1))
255 0 : pwi=AIMAG(vTot%pw(1:,3)) + AIMAG(c(1)%pw(1:,1)) + REAL(c(2)%pw(1:,1))
256 0 : vTot%pw(1:,3)=CMPLX(pwr,pwi)
257 :
258 : !vTot%pw_w(1:,1)=vTot%pw_w(1:,1)+c(3)%pw_w(1:,1)
259 : !vTot%pw_w(1:,2)=vTot%pw_w(1:,2)-c(3)%pw_w(1:,1)
260 : !pwr=REAL(vTot%pw_w(1:,3)) + REAL(c(1)%pw_w(1:,1)) - AIMAG(c(2)%pw_w(1:,1))
261 : !pwi=AIMAG(vTot%pw_w(1:,3)) + AIMAG(c(1)%pw_w(1:,1)) + REAL(c(2)%pw_w(1:,1))
262 : !vTot%pw_w(1:,3)=CMPLX(pwr,pwi)
263 :
264 : !vTot%vacz(1:,1:,1)=vTot%vacz(1:,1:,1)+c(3)%vacz(1:,1:,1)
265 : !vTot%vacz(1:,1:,2)=vTot%vacz(1:,1:,2)-c(3)%vacz(1:,1:,1)
266 : !vTot%vacz(1:,1:,3)=vTot%vacz(1:,1:,3)+c(1)%vacz(1:,1:,1)
267 : !vTot%vacz(1:,1:,4)=vTot%vacz(1:,1:,4)+c(2)%vacz(1:,1:,1)
268 :
269 : !vTot%vacxy(1:,1:,1:,1)=vTot%vacxy(1:,1:,1:,1)+c(3)%vacxy(1:,1:,1:,1)
270 : !vTot%vacxy(1:,1:,1:,2)=vTot%vacxy(1:,1:,1:,2)-c(3)%vacxy(1:,1:,1:,1)
271 : !vTot%vacxy(1:,1:,1:,3)=vTot%vacxy(1:,1:,1:,3)+c(1)%vacxy(1:,1:,1:,1)
272 : !vTot%vacxy(1:,1:,1:,4)=vTot%vacxy(1:,1:,1:,4)+c(2)%vacxy(1:,1:,1:,1)
273 :
274 0 : END SUBROUTINE correctPot
275 :
276 : END MODULE m_xcBfield
|
