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 346 : subroutine psqpw( mpi, atoms, sphhar, stars, vacuum, cell, input, sym, oneD, &
20 346 : & qpw, rho, rht, l_xyav, potdenType, psq )
21 :
22 : #include"cpp_double.h"
23 : use m_constants
24 : use m_phasy1
25 : use m_mpmom
26 : use m_sphbes
27 : use m_qsf
28 : use m_od_phasy
29 : use m_od_cylbes
30 : use m_types
31 : use m_DoubleFactorial
32 : use m_SphBessel
33 : implicit none
34 :
35 : type(t_mpi), intent(in) :: mpi
36 : type(t_atoms), intent(in) :: atoms
37 : type(t_sphhar), intent(in) :: sphhar
38 : type(t_stars), intent(in) :: stars
39 : type(t_vacuum), intent(in) :: vacuum
40 : type(t_cell), intent(in) :: cell
41 : type(t_input), intent(in) :: input
42 : type(t_sym), intent(in) :: sym
43 : type(t_oneD), intent(in) :: oneD
44 : logical, intent(in) :: l_xyav
45 : complex, intent(in) :: qpw(stars%ng3)
46 : real, intent(in) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
47 : real, intent(in) :: rht(vacuum%nmzd,2)
48 : integer, intent(in) :: potdenType
49 : complex, intent(out) :: psq(stars%ng3)
50 :
51 : complex :: psint, sa, sl, sm
52 : real :: f, fact, fpo, gz, p, qvac, rmtl, s, fJ, gr, g
53 : integer :: ivac, k, l, n, n1, nc, ncvn, lm, ll1, nd, m, nz
54 346 : complex :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
55 692 : complex :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
56 692 : real :: q2(vacuum%nmzd)
57 692 : real :: pn(0:atoms%lmaxd,atoms%ntype)
58 346 : real :: aj(0:atoms%lmaxd+maxval(atoms%ncv)+1)
59 692 : real :: rht1(vacuum%nmz)
60 346 : real, allocatable, dimension(:) :: il, kl
61 692 : real :: g0(atoms%ntype)
62 : #ifdef CPP_MPI
63 : include 'mpif.h'
64 : integer :: ierr(3)
65 346 : complex, allocatable :: c_b(:)
66 : #endif
67 :
68 : ! Calculate multipole moments
69 346 : call mpmom( input, mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, potdenType, qlm )
70 : #ifdef CPP_MPI
71 348950 : psq(:) = cmplx( 0.0, 0.0 )
72 346 : call MPI_BCAST( qpw, size(qpw), CPP_MPI_COMPLEX, 0, mpi%mpi_comm, ierr )
73 346 : nd = ( 2 * atoms%lmaxd + 1 ) * ( atoms%lmaxd + 1 ) * atoms%ntype
74 346 : call MPI_BCAST( qlm, nd, CPP_MPI_COMPLEX, 0, mpi%MPI_COMM, ierr )
75 : #endif
76 :
77 : ! prefactor in (A10) (Coulomb case) or (A11) (Yukawa case)
78 : ! nc(n) is the integer p in the paper; ncv(n) is l + p
79 : ! Coulomb case: pn(l,n) = (2 * l + 2 * p + 3)!! / ( (2 * l + 1)!! * R ** (ncv(n) + 1) ),
80 : ! Yukawa case: pn(l,n) = lambda ** (l + p + 1) / ( i_{l+p+1}(lambda * R) * (2 * l + 1)!! )
81 : ! g0 is the prefactor for the q=0 component in (A13)
82 1034 : pn = 0.
83 1722 : do n = 1, atoms%ntype
84 1034 : if ( potdenType /= POTDEN_TYPE_POTYUK ) then
85 6736 : do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
86 658 : pn(l, n) = DoubleFactorial( atoms%ncv(n) + 1, l ) / ( atoms%rmt(n) ** ( atoms%ncv(n) + 1 ) )
87 : end do
88 : else
89 30 : allocate( il(0:atoms%ncv(n)+1), kl(0:atoms%ncv(n)+1) )
90 30 : call ModSphBessel( il(0:), kl(0:), input%preconditioning_param * atoms%rmt(n), atoms%ncv(n) + 1 )
91 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)
92 300 : do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
93 300 : pn(l, n) = input%preconditioning_param ** ( atoms%ncv(n) + 1 ) / il( atoms%ncv(n) + 1 ) / DoubleFactorial( l )
94 : end do
95 30 : deallocate( il, kl )
96 : end if
97 : end do
98 :
99 : ! q=0 term: see (A12) (Coulomb case) or (A13) (Yukawa case)
100 346 : if( mpi%irank == 0 ) then
101 173 : s = 0.
102 517 : do n = 1, atoms%ntype
103 517 : if ( potdenType /= POTDEN_TYPE_POTYUK ) then
104 329 : s = s + atoms%neq(n) * real( qlm(0,0,n) )
105 : else
106 15 : s = s + atoms%neq(n) * real( qlm(0,0,n) ) * g0(n)
107 : end if
108 : end do
109 : !if( mpi%irank == 0 ) then
110 173 : psq(1) = qpw(1) + ( sfp_const / cell%omtil ) * s
111 : end if
112 :
113 : ! q/=0 term: see (A10) (Coulomb case) or (A11) (Yukawa case)
114 346 : fpo = 1. / cell%omtil
115 174475 : !$omp parallel do default( shared ) private( pylm, sa, n, ncvn, aj, sl, l, n1, ll1, sm, m, lm )
116 : do k = mpi%irank+2, stars%ng3, mpi%isize
117 174129 : if ( .not. oneD%odi%d1 ) then
118 174129 : call phasy1( atoms, stars, sym, cell, k, pylm )
119 : else
120 : call od_phasy( atoms%ntype, stars%ng3, atoms%nat, atoms%lmaxd, atoms%ntype, atoms%neq, &
121 0 : atoms%lmax, atoms%taual, cell%bmat, stars%kv3, k, oneD%odi, oneD%ods, pylm )
122 : end if
123 174129 : sa = 0.
124 463910 : do n = 1, atoms%ntype
125 289781 : ncvn = atoms%ncv(n)
126 289781 : call sphbes( ncvn + 1 , stars%sk3(k) * atoms%rmt(n), aj )
127 289781 : sl = 0.
128 2969192 : do l = 0, atoms%lmax(n)
129 2679411 : if ( l >= ncvn ) go to 60
130 2679411 : n1 = ncvn - l + 1
131 2679411 : ll1 = l * ( l + 1 ) + 1
132 2679411 : sm = 0.
133 27860592 : do m = -l, l
134 25181181 : lm = ll1 + m
135 27860592 : sm = sm + qlm(m,l,n) * conjg( pylm(lm,n) )
136 : end do
137 2969192 : 60 sl = sl + pn(l,n) / ( stars%sk3(k) ** n1 ) * aj( ncvn + 1 ) * sm
138 : end do
139 463910 : sa = sa + atoms%neq(n) * sl
140 : end do
141 174821 : psq(k) = qpw(k) + fpo * sa
142 : end do
143 : !$omp end parallel do
144 : #ifdef CPP_MPI
145 346 : allocate( c_b(stars%ng3) )
146 346 : call MPI_REDUCE( psq, c_b, stars%ng3, CPP_MPI_COMPLEX, MPI_SUM, 0, mpi%MPI_COMM, ierr )
147 346 : if ( mpi%irank == 0 ) then
148 174475 : psq(:stars%ng3) = c_b(:stars%ng3)
149 : end if
150 346 : deallocate( c_b )
151 : #endif
152 :
153 346 : if ( mpi%irank == 0 ) then
154 173 : if ( potdenType == POTDEN_TYPE_POTYUK ) return
155 : ! Check: integral of the pseudo charge density within the slab
156 170 : if ( input%film .and. .not. oneD%odi%d1 ) then
157 5 : psint = psq(1) * stars%nstr(1) * vacuum%dvac
158 32588 : do k = 2, stars%ng3
159 16294 : if ( stars%ig2(k) == 1 ) then
160 125 : gz = stars%kv3(3,k) * cell%bmat(3,3)
161 125 : f = 2. * sin( gz * cell%z1 ) / gz
162 125 : psint = psint + stars%nstr(k) * psq(k) * f
163 : end if
164 : end do
165 5 : psint = cell%area * psint
166 165 : else if ( input%film .and. oneD%odi%d1 ) then
167 0 : psint = (0.0, 0.0)
168 0 : do k = 2, stars%ng3
169 0 : if ( stars%kv3(3,k) == 0 ) then
170 : g = ( stars%kv3(1,k) * cell%bmat(1,1) + stars%kv3(2,k) * cell%bmat(2,1) ) ** 2 + &
171 0 : & ( stars%kv3(1,k) * cell%bmat(1,2) + stars%kv3(2,k) * cell%bmat(2,2) ) ** 2
172 0 : gr = sqrt( g )
173 0 : call od_cylbes( 1, gr * cell%z1, fJ )
174 0 : f = 2 * cell%vol * fJ / ( gr * cell%z1 )
175 0 : psint = psint + stars%nstr(k) * psq(k) * f
176 : end if
177 : end do
178 0 : psint = psint + psq(1) * stars%nstr(1) * cell%vol
179 165 : else if ( .not. input%film ) then
180 165 : psint = psq(1) * stars%nstr(1) * cell%omtil
181 : end if
182 170 : write( 6, fmt=8000 ) psint
183 : 8000 format (/,10x,'integral of pseudo charge density inside the slab='&
184 : & ,5x,2f11.6)
185 170 : if ( .not. input%film .or. potdenType == POTDEN_TYPE_POTYUK ) return
186 :
187 : ! Normalized pseudo density
188 5 : if ( .not. oneD%odi%d1 ) then
189 5 : qvac = 0.0
190 10 : do ivac = 1, vacuum%nvac
191 5 : call qsf( vacuum%delz, rht(1,ivac), q2, vacuum%nmz, 0 )
192 5 : q2(1) = q2(1) * cell%area
193 10 : qvac = qvac + q2(1) * 2. / real( vacuum%nvac )
194 : end do
195 5 : qvac = qvac - 2 * input%sigma
196 : else
197 0 : qvac = 0.0
198 0 : do nz = 1, vacuum%nmz
199 0 : rht1(nz) = ( cell%z1 + ( nz - 1 ) * vacuum%delz ) * rht(nz,vacuum%nvac)
200 : end do
201 0 : call qsf( vacuum%delz, rht1(1), q2, vacuum%nmz, 0 )
202 0 : qvac = cell%area * q2(1)
203 : end if
204 5 : if ( l_xyav ) return
205 5 : fact = ( qvac + psint ) / ( stars%nstr(1) * cell%vol )
206 5 : psq(1) = psq(1) - fact
207 5 : write( 6, fmt=8010 ) fact * 1000
208 : 8010 format (/,10x,' 1000 * normalization const. ='&
209 : & ,5x,2f11.6)
210 : end if ! mpi%irank == 0
211 :
212 : end subroutine psqpw
213 :
214 : end module m_psqpw
|
