Line data Source code
1 : module m_psqpw
2 : ! ***********************************************************
3 : ! generates the fourier coefficients of pseudo charge density
4 : !
5 : ! For yukawa_residual = .true. the subroutines calculate the
6 : ! pseudo charge density for the generation of the Yukawa
7 : ! potential instead of the Coulomb potential. This is used in
8 : ! the preconditioning of the SCF iteration for metallic systems.
9 : !
10 : ! references:
11 : ! for both the Coulomb and Yukawa cases:
12 : ! F. Tran, P. Blaha: Phys. Rev. B 83, 235118 (2011)
13 : ! or see the original paper for the normal Coulomb case only:
14 : ! M. Weinert: J. Math. Phys. 22(11) (1981) p.2434 eq. (10)-(15)
15 : ! ***********************************************************
16 :
17 : contains
18 :
19 682 : subroutine psqpw( fmpi, atoms, sphhar, stars, vacuum, cell, input, sym, juphon, &
20 682 : & den, ispin, l_xyav, potdenType, psq, sigma_disc, rhoimag, stars2, iDtype, iDir, rho0, qpw0, iDir2, mat2ord )
21 :
22 : #ifdef CPP_MPI
23 : use mpi
24 : #endif
25 : use m_constants
26 : use m_phasy1
27 : use m_mpmom
28 : use m_sphbes
29 : use m_qsf
30 : USE m_mpi_reduce_tool
31 :
32 :
33 : use m_types
34 : use m_DoubleFactorial
35 : use m_SphBessel
36 : implicit none
37 :
38 : type(t_mpi), intent(in) :: fmpi
39 : type(t_atoms), intent(in) :: atoms
40 : type(t_sphhar), intent(in) :: sphhar
41 : type(t_stars), intent(in) :: stars
42 : type(t_vacuum), intent(in) :: vacuum
43 : type(t_cell), intent(in) :: cell
44 : type(t_input), intent(in) :: input
45 : type(t_sym), intent(in) :: sym
46 : type(t_juphon), intent(in) :: juphon
47 : type(t_potden), intent(in) :: den
48 :
49 : logical, intent(in) :: l_xyav
50 : !complex, intent(in) :: qpw(stars%ng3)
51 : !real, intent(in) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
52 : !real, intent(in) :: rht(vacuum%nmzd,2)
53 : integer, intent(in) :: potdenType, ispin
54 : complex, intent(out) :: psq(stars%ng3)
55 : complex, intent(out) :: sigma_disc(2)
56 :
57 : REAL, OPTIONAL, INTENT(IN) :: rhoimag(atoms%jmtd,0:sphhar%nlhd,atoms%ntype), rho0(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
58 :
59 : TYPE(t_stars), OPTIONAL, INTENT(IN) :: stars2
60 : COMPLEX, OPTIONAL, INTENT(IN) :: qpw0(:)
61 :
62 : INTEGER, OPTIONAL, INTENT(IN) :: iDtype, iDir ! DFPT: Type and direction of displaced atom
63 : INTEGER, OPTIONAL, INTENT(IN) :: iDir2
64 : COMPLEX, OPTIONAL, INTENT(IN) :: mat2ord(5,3,3)
65 :
66 : complex :: psint, sa, sl, sm, qvac, fact, fc
67 : real :: f, fpo, gz, p, rmtl, s, fJ, gr, g
68 : integer :: ivac, k, l, n, n1, nc, ncvn, lm, ll1, nd, m, nz, kStart, kEnd
69 682 : complex :: psq_local(stars%ng3)
70 682 : complex :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
71 682 : complex :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
72 682 : complex :: q2(vacuum%nmzd)
73 682 : real :: q2r(vacuum%nmzd),q2i(vacuum%nmzd)
74 682 : real :: pn(0:atoms%lmaxd,atoms%ntype)
75 1908 : real :: aj(0:atoms%lmaxd+maxval(atoms%ncv)+1)
76 682 : real, allocatable, dimension(:) :: il, kl
77 682 : real :: g0(atoms%ntype)
78 682 : complex :: qpw(stars%ng3)
79 682 : real :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
80 682 : complex :: rht(vacuum%nmzd,2)
81 : LOGICAL :: l_dfptvgen ! If this is true, we handle things differently!
82 :
83 : #ifdef CPP_MPI
84 : integer :: ierr
85 : #endif
86 :
87 682 : l_dfptvgen = PRESENT(stars2)
88 1502016 : qpw = den%pw(:,ispin)
89 23987348 : rho = den%mt(:,:,:,ispin)
90 37830 : IF (input%film) rht = den%vac(:,1,:,ispin)
91 682 : sigma_disc = cmplx(0.0,0.0)
92 :
93 : ! Calculate multipole moments
94 682 : call timestart("mpmom")
95 682 : IF (.NOT.l_dfptvgen) THEN
96 682 : call mpmom( input, fmpi, atoms, sphhar, stars, sym, juphon, cell, qpw, rho, potdenType, qlm )
97 0 : ELSE IF (PRESENT(iDir2)) THEN
98 : ! DFPT case:
99 : ! Additional contributions to qlm due to surface corrections.
100 : call mpmom( input, fmpi, atoms, sphhar, stars, sym, juphon, cell, qpw, rho, potdenType, qlm, &
101 0 : & rhoimag=rhoimag, stars2=stars2, iDtype=iDtype, iDir=iDir, rho0=rho0, qpw0=qpw0, iDir2=iDir2, mat2ord=mat2ord )
102 : ELSE
103 : ! DFPT case:
104 : ! Additional contributions to qlm due to surface corrections.
105 : call mpmom( input, fmpi, atoms, sphhar, stars, sym, juphon, cell, qpw, rho, potdenType, qlm, &
106 0 : & rhoimag=rhoimag, stars2=stars2, iDtype=iDtype, iDir=iDir, rho0=rho0, qpw0=qpw0 )
107 : END IF
108 682 : call timestop("mpmom")
109 :
110 1502016 : psq(:) = cmplx( 0.0, 0.0 )
111 1502016 : psq_local(:) = cmplx( 0.0, 0.0 )
112 : #ifdef CPP_MPI
113 682 : call MPI_BCAST( qpw, size(qpw), MPI_DOUBLE_COMPLEX, 0, fmpi%mpi_comm, ierr )
114 682 : nd = ( 2 * atoms%lmaxd + 1 ) * ( atoms%lmaxd + 1 ) * atoms%ntype
115 682 : call MPI_BCAST( qlm, nd, MPI_DOUBLE_COMPLEX, 0, fmpi%MPI_COMM, ierr )
116 : #endif
117 :
118 : ! prefactor in (A10) (Coulomb case) or (A11) (Yukawa case)
119 : ! nc(n) is the integer p in the paper; ncv(n) is l + p
120 : ! Coulomb case: pn(l,n) = (2 * l + 2 * p + 3)!! / ( (2 * l + 1)!! * R ** (ncv(n) + 1) ),
121 : ! Yukawa case: pn(l,n) = lambda ** (l + p + 1) / ( i_{l+p+1}(lambda * R) * (2 * l + 1)!! )
122 : ! g0 is the prefactor for the q=0 component in (A13)
123 13558 : pn = 0.
124 1908 : do n = 1, atoms%ntype
125 1908 : if ( potdenType /= POTDEN_TYPE_POTYUK ) then
126 12088 : do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
127 12088 : pn(l, n) = DoubleFactorial( atoms%ncv(n) + 1, l ) / ( atoms%rmt(n) ** ( atoms%ncv(n) + 1 ) )
128 : end do
129 : else
130 30 : allocate( il(0:atoms%ncv(n)+1), kl(0:atoms%ncv(n)+1) )
131 30 : call ModSphBessel( il(0:), kl(0:), input%preconditioning_param * atoms%rmt(n), atoms%ncv(n) + 1 )
132 30 : g0(n) = ( input%preconditioning_param * atoms%rmt(n) ) ** ( atoms%ncv(n) + 1 ) / DoubleFactorial( atoms%ncv(n) + 1 ) / il( atoms%ncv(n) + 1 ) !p=ncv(n)
133 300 : do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
134 300 : pn(l, n) = input%preconditioning_param ** ( atoms%ncv(n) + 1 ) / il( atoms%ncv(n) + 1 ) / DoubleFactorial( l )
135 : end do
136 30 : deallocate( il, kl )
137 : end if
138 : end do
139 :
140 : ! q=0 term: see (A12) (Coulomb case) or (A13) (Yukawa case)
141 2728 : if( fmpi%irank == 0 .AND. norm2(stars%center)<=1e-8) then
142 : s = 0.
143 954 : do n = 1, atoms%ntype
144 954 : if ( potdenType /= POTDEN_TYPE_POTYUK ) then
145 598 : s = s + atoms%neq(n) * real( qlm(0,0,n) )
146 : else
147 15 : s = s + atoms%neq(n) * real( qlm(0,0,n) ) * g0(n)
148 : end if
149 : end do
150 : !if( fmpi%irank == 0 ) then
151 341 : psq_local(1) = qpw(1) + ( sfp_const / cell%omtil ) * s
152 : end if
153 :
154 : ! q/=0 term: see (A10) (Coulomb case) or (A11) (Yukawa case)
155 682 : fpo = 1. / cell%omtil
156 :
157 682 : call timestart("loop in psqpw")
158 :
159 2728 : CALL calcIndexBounds(fmpi, MERGE(2,1,norm2(stars%center)<=1e-8), stars%ng3, kStart, kEnd)
160 : !$OMP parallel do default( NONE ) &
161 : !$OMP SHARED(atoms,stars,sym,cell,kStart,kEnd,psq_local,qpw,qlm,pn,fpo,iDir,iDir2,iDtype) &
162 682 : !$OMP private( pylm, sa, n, ncvn, aj, sl, l, n1, ll1, sm, m, lm )
163 : do k = kStart, kEnd
164 : call phasy1( atoms, stars, sym, cell, k, pylm )
165 : sa = 0.0
166 : do n = 1, atoms%ntype
167 : ncvn = atoms%ncv(n)
168 : call sphbes( ncvn + 1 , stars%sk3(k) * atoms%rmt(n), aj )
169 : sl = 0.
170 : do l = 0, atoms%lmax(n)
171 : IF(l.LT.ncvn) THEN
172 : n1 = ncvn - l + 1
173 : ll1 = l * ( l + 1 ) + 1
174 : sm = 0.
175 : do m = -l, l
176 : lm = ll1 + m
177 : sm = sm + qlm(m,l,n) * conjg( pylm(lm,n) )
178 : end do
179 : END IF
180 : sl = sl + pn(l,n) / ( stars%sk3(k) ** n1 ) * aj( ncvn + 1 ) * sm
181 : end do
182 : sa = sa + atoms%neq(n) * sl
183 : IF (PRESENT(iDir2)) THEN
184 : IF (iDir2==iDir) THEN
185 : IF (n==iDtype.OR.0==iDtype) THEN
186 : sa = sa + atoms%zatom(n) * 5 * pn(2,n) * conjg(pylm(1,n)) * aj( ncvn + 1 ) / sfp_const / ( stars%sk3(k) ** (ncvn-1) )
187 : END IF
188 : END IF
189 : END IF
190 : end do
191 : psq_local(k) = qpw(k) + fpo * sa
192 : end do
193 : !$omp end parallel do
194 :
195 682 : CALL mpi_sum_reduce(psq_local,psq,fmpi%MPI_COMM)
196 :
197 682 : call timestop("loop in psqpw")
198 :
199 2728 : IF (.NOT.norm2(stars%center)<1e-8) RETURN ! TODO: Change this!
200 : !IF (l_dfptvgen) RETURN
201 : ! TODO: As of yet unclear if we need to do somecorrection here for
202 : ! DFPT, to make up for possible charge shifts from/to the vacuum.
203 :
204 682 : if ( fmpi%irank == 0 ) then
205 341 : if ( potdenType == POTDEN_TYPE_POTYUK ) return
206 : ! Check: integral of the pseudo charge density within the slab
207 338 : if ( input%film ) then
208 37 : psint = psq(1) * stars%nstr(1) * vacuum%dvac
209 138057 : do k = 2, stars%ng3
210 138057 : if ( stars%ig2(k) == 1 ) then
211 730 : gz = stars%kv3(3,k) * cell%bmat(3,3)
212 730 : f = 2. * sin( gz * cell%z1 ) / gz
213 730 : psint = psint + stars%nstr(k) * psq(k) * f
214 : end if
215 : end do
216 37 : psint = cell%area * psint
217 : else if ( .not. input%film ) then
218 301 : psint = psq(1) * stars%nstr(1) * cell%omtil
219 : end if
220 338 : write(oUnit, fmt=8000 ) psint
221 : 8000 format (/,10x,'integral of pseudo charge density inside the slab='&
222 : & ,5x,2f11.6)
223 338 : if ( .not. input%film .or. potdenType == POTDEN_TYPE_POTYUK ) return
224 :
225 : ! Normalized pseudo density
226 37 : qvac = cmplx(0.0,0.0)
227 84 : do ivac = 1, vacuum%nvac
228 47 : if (.not.l_dfptvgen) then
229 11844 : call qsf( vacuum%delz, real(rht(:,ivac)), q2r, vacuum%nmz, 0 )
230 11797 : q2 = q2r
231 : else
232 0 : call qsf( vacuum%delz, real(rht(:,ivac)), q2r, vacuum%nmz, 0 )
233 0 : call qsf( vacuum%delz,aimag(rht(:,ivac)), q2i, vacuum%nmz, 0 )
234 0 : q2 = q2r + ImagUnit*q2i
235 : end if
236 47 : q2(1) = q2(1) * cell%area
237 84 : qvac = qvac + q2(1) * 2. / real( vacuum%nvac )
238 : end do
239 : !TODO: reactivate electric fields
240 : !qvac = qvac - 2 * input%sigma
241 37 : if ( l_xyav ) return
242 37 : fact = ( qvac + psint ) / ( stars%nstr(1) * cell%vol )
243 37 : psq(1) = psq(1) - fact
244 : ! Find discontinuity at the boundaries
245 138094 : do k = 1, stars%ng3
246 138094 : if ( stars%ig2(k) == 1 ) then
247 767 : gz = stars%kv3(3,k) * cell%bmat(3,3)
248 767 : fc = cmplx(cos(gz * cell%z1),sin(gz * cell%z1))
249 767 : sigma_disc(1) = sigma_disc(1) - stars%nstr(k) * psq(k) * fc
250 767 : sigma_disc(2) = sigma_disc(2) + stars%nstr(k) * psq(k) * conjg(fc)
251 : end if
252 : end do
253 : sigma_disc(1) = sigma_disc(1) + rht(1,1)
254 : sigma_disc(2) = sigma_disc(2) - rht(1,2)
255 37 : sigma_disc = 0.0
256 : 8010 format (/,10x,' 1000 * normalization const. ='&
257 : & ,5x,2f11.6)
258 : end if ! fmpi%irank == 0
259 :
260 682 : end subroutine psqpw
261 :
262 : end module m_psqpw
|
