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 16 : 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 26 : 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 18 : 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,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_potden), DIMENSION(3), INTENT(INOUT) :: aVec
137 : TYPE(t_potden), INTENT(IN) :: vScal
138 : TYPE(t_potden), INTENT(OUT) :: vCorr
139 :
140 2 : TYPE(t_potden) :: div, phi, divloc
141 212 : TYPE(t_potden), DIMENSION(3) :: cvec, corrB
142 2 : TYPE(t_atoms) :: atloc
143 : INTEGER :: n, jr, lh, lhmax, jcut, nat ,i
144 : REAL :: xp(3,(atoms%lmaxd+1+mod(atoms%lmaxd+1,2))*(2*atoms%lmaxd+1))
145 2 : REAL, ALLOCATABLE :: intden(:,:)
146 2 : TYPE(t_gradients) :: tmp_grad
147 : complex :: sigma_loc(2)
148 :
149 : CALL div%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,atoms%ntype, &
150 : atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN, &
151 2 : vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
152 6 : ALLOCATE(div%pw_w,mold=div%pw)
153 6146 : div%pw_w = CMPLX(0.0,0.0)
154 :
155 2 : CALL timestart("Building divergence")
156 2 : CALL divergence(fmpi,input,stars,atoms,sphhar,vacuum,sym,cell,aVec,div)
157 2 : CALL timestop("Building divergence")
158 :
159 : ! Local atoms variable with no charges;
160 : ! needed for the potential generation from the divergence.
161 6 : atloc=atoms
162 6 : atloc%zatom=0.0
163 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)
164 4 : ALLOCATE(phi%pw_w(SIZE(phi%pw,1),size(phi%pw,2)))
165 6146 : phi%pw_w = CMPLX(0.0,0.0)
166 :
167 2 : CALL timestart("Building potential")
168 2 : sigma_loc = cmplx(0.0,0.0)
169 2 : CALL vgen_coulomb(1,fmpi ,input,field,vacuum,sym,juphon,stars,cell,sphhar,atloc,.TRUE.,div,phi,sigma_loc)
170 2 : CALL timestop("Building potential")
171 :
172 8 : DO i=1,3
173 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)
174 18 : ALLOCATE(cvec(i)%pw_w,mold=cvec(i)%pw)
175 18440 : cvec(i)%pw_w=CMPLX(0.0,0.0)
176 : ENDDO
177 :
178 2 : CALL timestart("Building correction field")
179 2 : CALL divpotgrad(input,stars,atloc,sphhar,vacuum,sym,cell,phi,cvec)
180 2 : CALL timestop("Building correction field")
181 :
182 2 : CALL init_pw_grid(stars,sym,cell)
183 8 : DO i=1,3
184 6 : CALL pw_to_grid(.FALSE.,1,.FALSE.,stars,cell,cvec(i)%pw,tmp_grad,rho=intden)
185 18438 : cvec(i)%pw=CMPLX(0.0,0.0)
186 18438 : cvec(i)%pw_w=CMPLX(0.0,0.0)
187 6 : CALL pw_from_grid(stars,intden,cvec(i)%pw,cvec(i)%pw_w)
188 8 : DEALLOCATE(intden)
189 : END DO !i
190 2 : CALL finish_pw_grid()
191 :
192 8 : DO i=1,3
193 : 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.,&
194 6 : POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
195 18 : ALLOCATE(corrB(i)%pw_w,mold=corrB(i)%pw)
196 18438 : corrB(i)%pw_w=CMPLX(0.0,0.0)
197 8 : CALL corrB(i)%addPotDen(aVec(i),cvec(i))
198 : ENDDO
199 :
200 2 : CALL BfieldtoVmat(sym, stars, atoms, sphhar, vacuum, vScal, corrB(1), corrB(2), corrB(3), vCorr)
201 :
202 8 : END SUBROUTINE sourcefree
203 :
204 0 : SUBROUTINE correctPot(vTot,c)
205 : USE m_types
206 :
207 : ! Takes a vectorial quantity c and saves its components into the appro
208 : ! priate components of the potential matrix V.
209 : !
210 : ! An initial V_mat = V*Id_(2x2) + sigma_vec*B_vec will become
211 : ! V_mat = V*Id_(2x2) + sigma_vec*(B_vec+C_vec)
212 : !
213 : ! TODO: Both quantities are assumed to be in the global frame of refe-
214 : ! rence. Make sure this is true. Also: consider SOC, LDA+U etc.
215 :
216 : TYPE(t_potden), INTENT(INOUT) :: vTot
217 : TYPE(t_potden), DIMENSION(3), INTENT(IN) :: c
218 :
219 0 : REAL :: pwr(SIZE(vTot%pw(1:,3))), pwi(SIZE(vTot%pw(1:,3)))
220 :
221 : !vtot%mt(:,0:,:,1)=(vtot%mt(:,0:,:,1)+vtot%mt(:,0:,:,2))/2
222 : !vtot%mt(:,0:,:,2)=vtot%mt(:,0:,:,1)
223 : !vtot%mt(:,0:,:,3:4)=0.0
224 :
225 : !vTot%pw(1:,1)=(vTot%pw(1:,1)+vTot%pw(1:,2))/2
226 : !vTot%pw(1:,2)=vTot%pw(1:,1)
227 : !vTot%pw(1:,3)=CMPLX(0.0,0.0)
228 :
229 : !vTot%pw_w(1:,1)=(vTot%pw_w(1:,1)+vTot%pw_w(1:,2))/2
230 : !vTot%pw_w(1:,2)=vTot%pw_w(1:,1)
231 : !vTot%pw_w(1:,3)=CMPLX(0.0,0.0)
232 :
233 : !vTot%mt(:,0:,:,1)=vTot%mt(:,0:,:,1)+b(3)%mt(:,0:,:,1)/2
234 : !vTot%mt(:,0:,:,2)=vTot%mt(:,0:,:,2)-b(3)%mt(:,0:,:,1)/2
235 : !vTot%mt(:,0:,:,3)=b(1)%mt(:,0:,:,1)/2
236 : !vTot%mt(:,0:,:,4)=b(2)%mt(:,0:,:,1)/2
237 :
238 : !vTot%pw(1:,1)=vTot%pw(1:,1)+b(3)%pw(1:,1)/2
239 : !vTot%pw(1:,2)=vTot%pw(1:,2)-b(3)%pw(1:,1)/2
240 : !vTot%pw(1:,3)=(b(1)%pw(1:,1)+ImagUnit*b(2)%pw(1:,1))/2
241 :
242 : !vTot%pw_w(1:,1)=vTot%pw_w(1:,1)+b(3)%pw_w(1:,1)/2
243 : !vTot%pw_w(1:,2)=vTot%pw_w(1:,2)-b(3)%pw_w(1:,1)/2
244 : !vTot%pw_w(1:,3)=(b(1)%pw_w(1:,1)+ImagUnit*b(2)%pw_w(1:,1))/2
245 :
246 0 : vTot%mt(:,0:,:,1)=c(3)%mt(:,0:,:,1)+vTot%mt(:,0:,:,1)
247 0 : vTot%mt(:,0:,:,2)=-c(3)%mt(:,0:,:,1)+vTot%mt(:,0:,:,2)
248 0 : vTot%mt(:,0:,:,3)=c(1)%mt(:,0:,:,1)+vTot%mt(:,0:,:,3)
249 0 : vTot%mt(:,0:,:,4)=c(2)%mt(:,0:,:,1)+vTot%mt(:,0:,:,4)
250 :
251 0 : vTot%pw(1:,1)=vTot%pw(1:,1)+c(3)%pw(1:,1)
252 0 : vTot%pw(1:,2)=vTot%pw(1:,2)-c(3)%pw(1:,1)
253 0 : pwr=REAL(vTot%pw(1:,3)) + REAL(c(1)%pw(1:,1)) - AIMAG(c(2)%pw(1:,1))
254 0 : pwi=AIMAG(vTot%pw(1:,3)) + AIMAG(c(1)%pw(1:,1)) + REAL(c(2)%pw(1:,1))
255 0 : vTot%pw(1:,3)=CMPLX(pwr,pwi)
256 :
257 : !vTot%pw_w(1:,1)=vTot%pw_w(1:,1)+c(3)%pw_w(1:,1)
258 : !vTot%pw_w(1:,2)=vTot%pw_w(1:,2)-c(3)%pw_w(1:,1)
259 : !pwr=REAL(vTot%pw_w(1:,3)) + REAL(c(1)%pw_w(1:,1)) - AIMAG(c(2)%pw_w(1:,1))
260 : !pwi=AIMAG(vTot%pw_w(1:,3)) + AIMAG(c(1)%pw_w(1:,1)) + REAL(c(2)%pw_w(1:,1))
261 : !vTot%pw_w(1:,3)=CMPLX(pwr,pwi)
262 :
263 : !vTot%vacz(1:,1:,1)=vTot%vacz(1:,1:,1)+c(3)%vacz(1:,1:,1)
264 : !vTot%vacz(1:,1:,2)=vTot%vacz(1:,1:,2)-c(3)%vacz(1:,1:,1)
265 : !vTot%vacz(1:,1:,3)=vTot%vacz(1:,1:,3)+c(1)%vacz(1:,1:,1)
266 : !vTot%vacz(1:,1:,4)=vTot%vacz(1:,1:,4)+c(2)%vacz(1:,1:,1)
267 :
268 : !vTot%vacxy(1:,1:,1:,1)=vTot%vacxy(1:,1:,1:,1)+c(3)%vacxy(1:,1:,1:,1)
269 : !vTot%vacxy(1:,1:,1:,2)=vTot%vacxy(1:,1:,1:,2)-c(3)%vacxy(1:,1:,1:,1)
270 : !vTot%vacxy(1:,1:,1:,3)=vTot%vacxy(1:,1:,1:,3)+c(1)%vacxy(1:,1:,1:,1)
271 : !vTot%vacxy(1:,1:,1:,4)=vTot%vacxy(1:,1:,1:,4)+c(2)%vacxy(1:,1:,1:,1)
272 :
273 0 : END SUBROUTINE correctPot
274 :
275 : END MODULE m_xcBfield
|
