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 1768 : subroutine psqpw( fmpi, atoms, sphhar, stars, vacuum, cell, input, sym, juphon, &
20 1768 : & den, ispin, l_xyav, potdenType, psq, sigma_disc, rhoimag, stars2, iDtype, iDir, rho0, iDir2 )
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 : type(t_potden),optional,intent(in) :: rhoimag, rho0
58 :
59 : TYPE(t_stars), OPTIONAL, INTENT(IN) :: stars2
60 :
61 : INTEGER, OPTIONAL, INTENT(IN) :: iDtype, iDir ! DFPT: Type and direction of displaced atom
62 : INTEGER, OPTIONAL, INTENT(IN) :: iDir2
63 :
64 : complex :: psint, sa, sl, sm, qvac, fact, fc
65 : real :: f, fpo, gz, p, rmtl, s, fJ, gr, g
66 : integer :: ivac, k, l, n, n1, nc, ncvn, lm, ll1, nd, m, nz, kStart, kEnd
67 1768 : complex :: psq_local(stars%ng3)
68 1768 : complex :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
69 1768 : complex :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
70 1768 : complex :: q2(vacuum%nmzd)
71 1768 : real :: q2r(vacuum%nmzd),q2i(vacuum%nmzd)
72 1768 : real :: pn(0:atoms%lmaxd,atoms%ntype)
73 4618 : real :: aj(0:atoms%lmaxd+maxval(atoms%ncv)+1)
74 1768 : real, allocatable, dimension(:) :: il, kl
75 1768 : real :: g0(atoms%ntype)
76 1768 : complex :: qpw(stars%ng3)
77 1768 : real :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
78 1768 : complex :: rht(vacuum%nmzd,2)
79 : LOGICAL :: l_dfptvgen ! If this is true, we handle things differently!
80 :
81 : #ifdef CPP_MPI
82 : integer :: ierr
83 : #endif
84 :
85 1768 : l_dfptvgen = PRESENT(stars2)
86 4329175 : qpw = den%pw(:,ispin)
87 100062550 : rho = den%mt(:,:,:,ispin)
88 260298 : IF (input%film) rht = den%vac(:,1,:,ispin)
89 1768 : sigma_disc = cmplx(0.0,0.0)
90 :
91 : ! Calculate multipole moments
92 1768 : call timestart("mpmom")
93 : ! DFPT case:
94 : ! Additional contributions to qlm due to surface corrections.
95 : call mpmom( input, fmpi, atoms, sphhar, stars, sym, juphon, cell, qpw, rho, potdenType, qlm, ispin, &
96 1768 : & rhoimag=rhoimag, stars2=stars2, iDtype=iDtype, iDir=iDir, rho0=rho0, iDir2=iDir2 )
97 1768 : call timestop("mpmom")
98 :
99 4329175 : psq(:) = cmplx( 0.0, 0.0 )
100 4329175 : psq_local(:) = cmplx( 0.0, 0.0 )
101 : #ifdef CPP_MPI
102 1768 : call MPI_BCAST( qpw, size(qpw), MPI_DOUBLE_COMPLEX, 0, fmpi%mpi_comm, ierr )
103 1768 : nd = ( 2 * atoms%lmaxd + 1 ) * ( atoms%lmaxd + 1 ) * atoms%ntype
104 1768 : call MPI_BCAST( qlm, nd, MPI_DOUBLE_COMPLEX, 0, fmpi%MPI_COMM, ierr )
105 : #endif
106 :
107 : ! prefactor in (A10) (Coulomb case) or (A11) (Yukawa case)
108 : ! nc(n) is the integer p in the paper; ncv(n) is l + p
109 : ! Coulomb case: pn(l,n) = (2 * l + 2 * p + 3)!! / ( (2 * l + 1)!! * R ** (ncv(n) + 1) ),
110 : ! Yukawa case: pn(l,n) = lambda ** (l + p + 1) / ( i_{l+p+1}(lambda * R) * (2 * l + 1)!! )
111 : ! g0 is the prefactor for the q=0 component in (A13)
112 28320 : pn = 0.
113 4618 : do n = 1, atoms%ntype
114 4618 : if ( potdenType /= POTDEN_TYPE_POTYUK ) then
115 23924 : do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
116 23924 : pn(l, n) = DoubleFactorial( atoms%ncv(n) + 1, l ) / ( atoms%rmt(n) ** ( atoms%ncv(n) + 1 ) )
117 : end do
118 : else
119 40 : allocate( il(0:atoms%ncv(n)+1), kl(0:atoms%ncv(n)+1) )
120 40 : call ModSphBessel( il(0:), kl(0:), input%preconditioning_param * atoms%rmt(n), atoms%ncv(n) + 1 )
121 40 : 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)
122 360 : do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
123 360 : pn(l, n) = input%preconditioning_param ** ( atoms%ncv(n) + 1 ) / il( atoms%ncv(n) + 1 ) / DoubleFactorial( l )
124 : end do
125 40 : deallocate( il, kl )
126 : end if
127 : end do
128 :
129 : ! q=0 term: see (A12) (Coulomb case) or (A13) (Yukawa case)
130 7072 : if( fmpi%irank == 0 .AND. norm2(stars%center)<=1e-8) then
131 : s = 0.
132 2282 : do n = 1, atoms%ntype
133 2282 : if ( potdenType /= POTDEN_TYPE_POTYUK ) then
134 1385 : s = s + atoms%neq(n) * real( qlm(0,0,n) )
135 : else
136 25 : s = s + atoms%neq(n) * real( qlm(0,0,n) ) * g0(n)
137 : end if
138 : end do
139 : !if( fmpi%irank == 0 ) then
140 872 : psq_local(1) = qpw(1) + ( sfp_const / cell%omtil ) * s
141 : end if
142 :
143 : ! q/=0 term: see (A10) (Coulomb case) or (A11) (Yukawa case)
144 1768 : fpo = 1. / cell%omtil
145 :
146 1768 : call timestart("loop in psqpw")
147 :
148 7482 : CALL calcIndexBounds(fmpi, MERGE(2,1,norm2(stars%center)<=1e-8), stars%ng3, kStart, kEnd)
149 : !$OMP parallel do default( NONE ) &
150 : !$OMP SHARED(atoms,stars,sym,cell,kStart,kEnd,psq_local,qpw,qlm,pn,fpo,iDir,iDir2,iDtype) &
151 1768 : !$OMP private( pylm, sa, n, ncvn, aj, sl, l, n1, ll1, sm, m, lm )
152 : do k = kStart, kEnd
153 : call phasy1( atoms, stars, sym, cell, k, pylm )
154 : sa = 0.0
155 : do n = 1, atoms%ntype
156 : ncvn = atoms%ncv(n)
157 : call sphbes( ncvn + 1 , stars%sk3(k) * atoms%rmt(n), aj )
158 : sl = 0.
159 : do l = 0, atoms%lmax(n)
160 : IF(l.LT.ncvn) THEN
161 : n1 = ncvn - l + 1
162 : ll1 = l * ( l + 1 ) + 1
163 : sm = 0.
164 : do m = -l, l
165 : lm = ll1 + m
166 : sm = sm + qlm(m,l,n) * conjg( pylm(lm,n) )
167 : end do
168 : END IF
169 : sl = sl + pn(l,n) / ( stars%sk3(k) ** n1 ) * aj( ncvn + 1 ) * sm
170 : end do
171 : sa = sa + atoms%neq(n) * sl
172 : IF (PRESENT(iDir2)) THEN
173 : IF (iDir2==iDir) THEN
174 : IF (n==iDtype.OR.0==iDtype) THEN
175 : sa = sa + atoms%zatom(n) * 5 * pn(2,n) * conjg(pylm(1,n)) * aj( ncvn + 1 ) / sfp_const / ( stars%sk3(k) ** (ncvn-1) )
176 : END IF
177 : END IF
178 : END IF
179 : end do
180 : psq_local(k) = qpw(k) + fpo * sa
181 : end do
182 : !$omp end parallel do
183 :
184 1768 : CALL mpi_sum_reduce(psq_local,psq,fmpi%MPI_COMM)
185 :
186 1768 : call timestop("loop in psqpw")
187 :
188 7072 : IF (.NOT.norm2(stars%center)<1e-8) RETURN ! TODO: Change this!
189 : !IF (l_dfptvgen) RETURN
190 : ! TODO: As of yet unclear if we need to do somecorrection here for
191 : ! DFPT, to make up for possible charge shifts from/to the vacuum.
192 :
193 1358 : if ( fmpi%irank == 0 ) then
194 872 : if ( potdenType == POTDEN_TYPE_POTYUK ) return
195 : ! Check: integral of the pseudo charge density within the slab
196 864 : if ( input%film ) then
197 259 : psint = psq(1) * stars%nstr(1) * vacuum%dvac
198 1325603 : do k = 2, stars%ng3
199 1325603 : if ( stars%ig2(k) == 1 ) then
200 10644 : gz = stars%kv3(3,k) * cell%bmat(3,3)
201 10644 : f = 2. * sin( gz * cell%z1 ) / gz
202 10644 : psint = psint + stars%nstr(k) * psq(k) * f
203 : end if
204 : end do
205 259 : psint = cell%area * psint
206 : else if ( .not. input%film ) then
207 605 : psint = psq(1) * stars%nstr(1) * cell%omtil
208 : end if
209 864 : write(oUnit, fmt=8000 ) psint
210 : 8000 format (/,10x,'integral of pseudo charge density inside the slab='&
211 : & ,5x,2f11.6)
212 864 : if ( .not. input%film .or. potdenType == POTDEN_TYPE_POTYUK ) return
213 :
214 : ! Normalized pseudo density
215 259 : qvac = cmplx(0.0,0.0)
216 752 : do ivac = 1, vacuum%nvac
217 493 : if (.not.l_dfptvgen) then
218 18396 : call qsf( vacuum%delz, real(rht(:,ivac)), q2r, vacuum%nmz, 0 )
219 18323 : q2 = q2r
220 : else
221 105840 : call qsf( vacuum%delz, real(rht(:,ivac)), q2r, vacuum%nmz, 0 )
222 105420 : call qsf( vacuum%delz,aimag(rht(:,ivac)), q2i, vacuum%nmz, 0 )
223 105420 : q2 = q2r + ImagUnit*q2i
224 : end if
225 493 : q2(1) = q2(1) * cell%area
226 752 : qvac = qvac + q2(1) * 2. / real( vacuum%nvac )
227 : end do
228 : !TODO: reactivate electric fields
229 : !qvac = qvac - 2 * input%sigma
230 259 : if ( l_xyav ) return
231 259 : fact = ( qvac + psint ) / ( stars%nstr(1) * cell%vol )
232 259 : psq(1) = psq(1) - fact
233 : ! Find discontinuity at the boundaries
234 1325862 : do k = 1, stars%ng3
235 1325862 : if ( stars%ig2(k) == 1 ) then
236 10903 : gz = stars%kv3(3,k) * cell%bmat(3,3)
237 10903 : fc = cmplx(cos(gz * cell%z1),sin(gz * cell%z1))
238 10903 : sigma_disc(1) = sigma_disc(1) - stars%nstr(k) * psq(k) * fc
239 10903 : sigma_disc(2) = sigma_disc(2) + stars%nstr(k) * psq(k) * conjg(fc)
240 : end if
241 : end do
242 : sigma_disc(1) = sigma_disc(1) + rht(1,1)
243 : sigma_disc(2) = sigma_disc(2) - rht(1,2)
244 259 : sigma_disc = 0.0
245 : 8010 format (/,10x,' 1000 * normalization const. ='&
246 : & ,5x,2f11.6)
247 : end if ! fmpi%irank == 0
248 :
249 : end subroutine psqpw
250 :
251 : end module m_psqpw
|
