Line data Source code
1 : module m_vmts
2 : #ifdef CPP_MPI
3 : use mpi
4 : #endif
5 : contains
6 :
7 682 : subroutine vmts( input, fmpi, stars, sphhar, atoms, sym, cell, juphon, dosf, vpw, rho, potdenType, vr, rhoIm, vrIm, iDtype, iDir, iDir2, mat2ord )
8 :
9 : !-------------------------------------------------------------------------
10 : ! This subroutine calculates the lattice harmonics expansion coefficients
11 : ! of the Coulomb / Yukawa potential for all atom types.
12 : !
13 : ! In more detail:
14 : ! the radius-dependent coefficient of the potential is
15 : ! V_{lm}(r) = \int_0^R G_l(r,r') \rho_{lm}(r') r'^2 dr' + V_{lm}(R) P_l(r)
16 : ! where
17 : ! [sphere boundary contribution - first part of code:]
18 : ! V_{lm}(R) is the spherical harmonics expansion of the interstitial
19 : ! potential,
20 : ! P_l(r) is the derivative of the spherical Green's function on the sphere
21 : ! boundary (times a constant),
22 : ! [sphere interior contribution - second part of code:]
23 : ! G_l(r,r') ) = green_factor u_1(r_<) u_2(r_>) is the spherical Green's
24 : ! function with homogeneous solutions u_1 and u_2 and
25 : ! r_< = min(r,r') and r_> = max(r,r')
26 : ! The integral is split in a part where r'=r_< and a part where r'=r_> and
27 : ! the integral from r to R is split in \int_0^R - \int_0^r. Resulting
28 : ! terms depending solely on R (not on r) are summarised in the variable
29 : ! termsR.
30 : !
31 : ! More information and equations can be found in
32 : ! F. Tran, P. Blaha: Phys. Rev. B 83, 235118(2011)
33 : !-------------------------------------------------------------------------
34 :
35 : use m_constants
36 : use m_types
37 : use m_mpi_reduce_tool
38 : use m_intgr, only : intgr2, sfint
39 : use m_phasy1
40 : use m_sphbes
41 :
42 : use m_SphBessel
43 : !$ use omp_lib
44 : implicit none
45 :
46 : type(t_input), intent(in) :: input
47 : type(t_mpi), intent(in) :: fmpi
48 : type(t_stars), intent(in) :: stars
49 : type(t_sphhar), intent(in) :: sphhar
50 : type(t_atoms), intent(in) :: atoms
51 : type(t_sym), intent(in) :: sym
52 : type(t_cell), intent(in) :: cell
53 : type(t_juphon), intent(in) :: juphon
54 :
55 : LOGICAL, INTENT(IN) :: dosf
56 : complex, intent(in) :: vpw(:)!(stars%ng3,input%jspins)
57 : real, intent(in) :: rho(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
58 : integer, intent(in) :: potdenType
59 : real, intent(out) :: vr(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
60 : REAL, OPTIONAL, INTENT(IN) :: rhoIm(:,0:,:)
61 : REAL, OPTIONAL, INTENT(OUT) :: vrIm(:,0:,:)
62 : INTEGER, OPTIONAL, INTENT(IN) :: iDtype, iDir
63 : INTEGER, OPTIONAL, INTENT(IN) :: iDir2
64 : COMPLEX, OPTIONAL, INTENT(IN) :: mat2ord(5,3,3)
65 :
66 : complex :: cp, sm
67 : integer :: i, jm, k, l, lh, n, nd, lm, m, imax, lmax, iMem, ptsym
68 : integer :: maxBunchSize, numBunches, iBunch, firstStar, iTempArray
69 : complex :: temp
70 682 : complex :: vtl(0:sphhar%nlhd, atoms%ntype)
71 682 : complex :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
72 : real :: green_factor, termsR, pref
73 682 : real :: green_1 (1:atoms%jmtd), green_2 (1:atoms%jmtd)
74 682 : real :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd)
75 682 : real :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)!, integral_3 (1:atoms%jmtd)
76 682 : real :: sbf(0:atoms%lmaxd)
77 682 : real, allocatable, dimension(:,:) :: il, kl
78 : LOGICAL :: l_dfptvgen
79 682 : COMPLEX, ALLOCATABLE :: vtlStars(:,:,:), vtlLocal(:,:)
80 : TYPE(t_parallelLoop) :: mpiLoop, ompLoop
81 :
82 682 : l_dfptvgen = PRESENT(rhoIm)
83 :
84 : ! SPHERE BOUNDARY CONTRIBUTION to the coefficients calculated from the values
85 : ! of the interstitial Coulomb / Yukawa potential on the sphere boundary
86 :
87 1364 : ALLOCATE (vtlLocal(0:sphhar%nlhd,atoms%ntype))
88 34318 : vtlLocal(:,:) = cmplx(0.0,0.0)
89 :
90 2728 : firstStar = MERGE(2,1,norm2(stars%center)<=1e-8)
91 682 : maxBunchSize = 2*getNumberOfThreads() ! This bunch size is kind of a magic number determined from some
92 : ! naive performance tests for a 64 atom unit cell
93 682 : CALL calcNumberComputationBunches(firstStar, stars%ng3, maxBunchSize, numBunches)
94 :
95 682 : CALL mpiLoop%init(fmpi%irank,fmpi%isize,0,numBunches-1)
96 :
97 2046 : ALLOCATE(vtlStars(0:sphhar%nlhd,atoms%ntype,maxBunchSize))
98 137954 : vtlStars = CMPLX(0.0,0.0)
99 :
100 : ! q/=0 components
101 188327 : DO iBunch = mpiLoop%bunchMinIndex, mpiLoop%bunchMaxIndex
102 187645 : CALL ompLoop%init(iBunch,numBunches,firstStar,stars%ng3)
103 : !$OMP parallel do default( NONE ) &
104 : !$OMP SHARED(ompLoop,atoms,stars,sym,cell,sphhar,vpw,vtlStars,potdenType) &
105 188327 : !$OMP private(iTempArray,cp,pylm,n,sbf,nd,lh,l,sm,jm,m,lm)
106 : do k = ompLoop%bunchMinIndex, ompLoop%bunchMaxIndex
107 : iTempArray = k - ompLoop%bunchMinIndex + 1
108 : cp = vpw(k) * stars%nstr(k)
109 : call phasy1( atoms, stars, sym, cell, k, pylm )
110 : do n = 1, atoms%ntype
111 : call sphbes( atoms%lmax(n), stars%sk3(k) * atoms%rmt(n), sbf )
112 : nd = sym%ntypsy(atoms%firstAtom(n))
113 : do lh = 0, sphhar%nlh(nd)
114 : l = sphhar%llh(lh,nd)
115 : sm = (0.0,0.0)
116 : do jm = 1, sphhar%nmem(lh,nd)
117 : m = sphhar%mlh(jm,lh,nd)
118 : lm = l * ( l + 1 ) + m + 1
119 : sm = sm + conjg( sphhar%clnu(jm,lh,nd) ) * pylm(lm,n)
120 : end do
121 : vtlStars(lh,n,iTempArray) = vtlStars(lh,n,iTempArray) + cp * sbf(l) * sm
122 : end do
123 : end do
124 : end do
125 : !$OMP END PARALLEL DO
126 : END DO
127 :
128 : !$OMP parallel do default( NONE ) &
129 : !$OMP SHARED(atoms,sym,sphhar,maxBunchSize,vtlStars,vtlLocal) &
130 682 : !$OMP private(nd,lh,temp,iTempArray)
131 : DO n = 1, atoms%ntype
132 : nd = sym%ntypsy(atoms%firstAtom(n))
133 : do lh = 0, sphhar%nlh(nd)
134 : temp = CMPLX(0.0,0.0)
135 : DO iTempArray = 1, maxBunchSize
136 : temp = temp + vtlStars(lh,n,iTempArray)
137 : END DO
138 : vtlLocal(lh,n) = temp
139 : END DO
140 : END DO
141 : !$OMP END PARALLEL DO
142 :
143 : ! q=0 component
144 2728 : if ( fmpi%irank == 0 .AND. norm2(stars%center)<=1e-8 ) then
145 954 : DO n = 1, atoms%ntype
146 954 : vtlLocal(0,n) = vtlLocal(0,n) + sfp_const * vpw(1)
147 : END DO
148 : end if
149 :
150 682 : CALL mpi_sum_reduce(vtlLocal, vtl,fmpi%mpi_comm)
151 :
152 : ! SPHERE INTERIOR CONTRIBUTION to the coefficients calculated from the
153 : ! values of the sphere Coulomb/Yukawa potential on the sphere boundary
154 :
155 682 : if( fmpi%irank == 0 ) then
156 341 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
157 6 : allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
158 : end if
159 :
160 954 : do n = 1, atoms%ntype
161 613 : nd = sym%ntypsy(atoms%firstAtom(n))
162 613 : imax = atoms%jri(n)
163 613 : lmax = sphhar%llh(sphhar%nlh(nd), nd)
164 613 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
165 : !do concurrent (i = 1:imax)
166 11370 : do i = 1,imax
167 11370 : call ModSphBessel( il(0:,i), kl(0:,i), input%preconditioning_param * atoms%rmsh(i,n), lmax )
168 : end do
169 : end if
170 16477 : do lh = 0, sphhar%nlh(nd)
171 15523 : l = sphhar%llh(lh,nd)
172 15523 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
173 156906 : green_1(1:imax) = il(l,1:imax)
174 156906 : green_2(1:imax) = kl(l,1:imax)
175 207 : green_factor = fpi_const * input%preconditioning_param
176 : else
177 11053930 : green_1(1:imax) = atoms%rmsh(1:imax,n) ** l
178 11053930 : green_2(1:imax) = 1.0 / ( green_1(1:imax) * atoms%rmsh(1:imax,n) )
179 15316 : green_factor = fpi_const / ( 2 * l + 1 )
180 : end if
181 :
182 11210836 : integrand_1(1:imax) = green_1(1:imax) * rho(1:imax,lh,n)
183 11210836 : integrand_2(1:imax) = green_2(1:imax) * rho(1:imax,lh,n)
184 15523 : if (.not.dosf) THEN
185 15361 : call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
186 15361 : call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
187 : else
188 162 : call sfint(integrand_1(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_1(1:imax))
189 162 : call sfint(integrand_2(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_2(1:imax))
190 : end if
191 15523 : termsR = integral_2(imax) + ( vtl(lh,n) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
192 : vr(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
193 11210836 : + green_2(1:imax) * integral_1(1:imax) )
194 16136 : IF (l_dfptvgen) THEN
195 : ! Integrate the imaginary part of the density perturbation as well.
196 0 : integrand_1(1:imax) = green_1(1:imax) * rhoIm(1:imax,lh,n)
197 0 : integrand_2(1:imax) = green_2(1:imax) * rhoIm(1:imax,lh,n)
198 0 : call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
199 0 : call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
200 0 : termsR = integral_2(imax) + ( AIMAG(vtl(lh,n)) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
201 : vrIm(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
202 0 : + green_2(1:imax) * integral_1(1:imax) )
203 : END IF
204 : end do
205 : end do
206 341 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
207 3 : deallocate( il, kl )
208 : end if
209 :
210 341 : if ( potdenType /= POTDEN_TYPE_POTYUK .AND. potdenType /= POTDEN_TYPE_CRYSTALFIELD) then
211 337 : IF (.NOT.l_dfptvgen) THEN
212 930 : do n = 1, atoms%ntype
213 412483 : vr(1:atoms%jri(n),0,n) = vr(1:atoms%jri(n),0,n) - sfp_const * ( 1.0 / atoms%rmsh(1:atoms%jri(n),n) - 1.0 / atoms%rmt(n) ) * atoms%zatom(n)
214 : end do
215 0 : ELSE IF (juphon%l_phonon) THEN
216 0 : IF (.NOT.PRESENT(iDir2)) THEN
217 : ! DFPT(-phonon) case:
218 : ! l=1 contributions from the Coulomb singularity instead of l=0 (1/r -> 1/r^2)
219 0 : DO n = MERGE(1,iDtype,iDtype==0), MERGE(atoms%ntype,iDtype,iDtype==0)
220 0 : ptsym = sym%ntypsy(atoms%firstAtom(n))
221 0 : pref = MERGE(atoms%zatom(n),-atoms%zatom(n),iDtype==0)
222 0 : DO lh = 1, 3
223 0 : l = sphhar%llh(lh, ptsym)
224 0 : DO iMem = 1, sphhar%nmem(lh, ptsym)
225 0 : m = sphhar%mlh(iMem, lh, ptsym)
226 0 : lm = l*(l+1) + m + 1
227 : vr(1:atoms%jri(n),lh,n) = vr(1:atoms%jri(n),lh,n) + &
228 : conjg(sphhar%clnu(iMem, lh, ptsym)) * c_im(iDir, lm - 1) * pref * &
229 0 : ( 1 - (atoms%rmsh(1:atoms%jri(n), n) / atoms%rmt(n))**3) / atoms%rmsh(1:atoms%jri(n),n)**2
230 : END DO
231 : END DO
232 : END DO
233 : ELSE
234 : ! DFPT 2nd order case:
235 : ! l=2 contributions from the Coulomb singularity instead of l=0 (1/r -> 1/r^3)
236 0 : DO n = 1, atoms%ntype!MERGE(1,iDtype,iDtype==0), MERGE(atoms%ntype,iDtype,iDtype==0)
237 0 : ptsym = sym%ntypsy(atoms%firstAtom(n))
238 0 : pref = -atoms%zatom(n)!MERGE(atoms%zatom(n),-atoms%zatom(n),iDtype==0)
239 0 : DO lh = 4, 8
240 0 : l = sphhar%llh(lh, ptsym)
241 0 : DO iMem = 1, sphhar%nmem(lh, ptsym)
242 0 : m = sphhar%mlh(iMem, lh, ptsym)
243 0 : lm = l*(l+1) + m + 1
244 0 : IF ((n.EQ.iDtype).OR.(0.EQ.iDtype)) vr(1:atoms%jri(n),lh,n) = vr(1:atoms%jri(n),lh,n) + &
245 : conjg(sphhar%clnu(iMem, lh, ptsym)) * mat2ord(lm-4,iDir2,iDir) * pref * &
246 0 : ( 1 - (atoms%rmsh(1:atoms%jri(n), n) / atoms%rmt(n))**5) / atoms%rmsh(1:atoms%jri(n),n)**3
247 : END DO
248 : END DO
249 : END DO
250 : END IF
251 : END IF
252 : end if
253 : end if
254 :
255 1364 : end subroutine vmts
256 :
257 : end module m_vmts
|