Line data Source code
1 : module m_vmts
2 : #ifdef CPP_MPI
3 : use mpi
4 : #endif
5 : contains
6 :
7 696 : subroutine vmts( input, fmpi, stars, sphhar, atoms, sym, cell, 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 :
54 : LOGICAL, INTENT(IN) :: dosf
55 : complex, intent(in) :: vpw(:)!(stars%ng3,input%jspins)
56 : real, intent(in) :: rho(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
57 : integer, intent(in) :: potdenType
58 : real, intent(out) :: vr(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
59 : REAL, OPTIONAL, INTENT(IN) :: rhoIm(:,0:,:)
60 : REAL, OPTIONAL, INTENT(OUT) :: vrIm(:,0:,:)
61 : INTEGER, OPTIONAL, INTENT(IN) :: iDtype, iDir
62 : INTEGER, OPTIONAL, INTENT(IN) :: iDir2
63 : COMPLEX, OPTIONAL, INTENT(IN) :: mat2ord(5,3,3)
64 :
65 : complex :: cp, sm
66 : integer :: i, jm, k, l, lh, n, nd, lm, m, imax, lmax, iMem, ptsym
67 : integer :: maxBunchSize, numBunches, iBunch, firstStar, iTempArray
68 : complex :: temp
69 696 : complex :: vtl(0:sphhar%nlhd, atoms%ntype)
70 696 : complex :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
71 : real :: green_factor, termsR, pref
72 696 : real :: green_1 (1:atoms%jmtd), green_2 (1:atoms%jmtd)
73 696 : real :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd)
74 696 : real :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)!, integral_3 (1:atoms%jmtd)
75 696 : real :: sbf(0:atoms%lmaxd)
76 696 : real, allocatable, dimension(:,:) :: il, kl
77 : LOGICAL :: l_dfptvgen
78 696 : COMPLEX, ALLOCATABLE :: vtlStars(:,:,:), vtlLocal(:,:)
79 : TYPE(t_parallelLoop) :: mpiLoop, ompLoop
80 :
81 696 : l_dfptvgen = PRESENT(rhoIm)
82 :
83 : ! SPHERE BOUNDARY CONTRIBUTION to the coefficients calculated from the values
84 : ! of the interstitial Coulomb / Yukawa potential on the sphere boundary
85 :
86 2784 : ALLOCATE (vtlLocal(0:sphhar%nlhd,atoms%ntype))
87 35144 : vtlLocal(:,:) = cmplx(0.0,0.0)
88 :
89 2784 : firstStar = MERGE(2,1,norm2(stars%center)<=1e-8)
90 696 : maxBunchSize = 2*getNumberOfThreads() ! This bunch size is kind of a magic number determined from some
91 : ! naive performance tests for a 64 atom unit cell
92 696 : CALL calcNumberComputationBunches(firstStar, stars%ng3, maxBunchSize, numBunches)
93 :
94 696 : CALL mpiLoop%init(fmpi%irank,fmpi%isize,0,numBunches-1)
95 :
96 3480 : ALLOCATE(vtlStars(0:sphhar%nlhd,atoms%ntype,maxBunchSize))
97 141272 : vtlStars = CMPLX(0.0,0.0)
98 :
99 : ! q/=0 components
100 195887 : DO iBunch = mpiLoop%bunchMinIndex, mpiLoop%bunchMaxIndex
101 195191 : CALL ompLoop%init(iBunch,numBunches,firstStar,stars%ng3)
102 : !$OMP parallel do default( NONE ) &
103 : !$OMP SHARED(ompLoop,atoms,stars,sym,cell,sphhar,vpw,vtlStars,potdenType) &
104 195887 : !$OMP private(iTempArray,cp,pylm,n,sbf,nd,lh,l,sm,jm,m,lm)
105 : do k = ompLoop%bunchMinIndex, ompLoop%bunchMaxIndex
106 : iTempArray = k - ompLoop%bunchMinIndex + 1
107 : cp = vpw(k) * stars%nstr(k)
108 : call phasy1( atoms, stars, sym, cell, k, pylm )
109 : do n = 1, atoms%ntype
110 : call sphbes( atoms%lmax(n), stars%sk3(k) * atoms%rmt(n), sbf )
111 : nd = sym%ntypsy(atoms%firstAtom(n))
112 : do lh = 0, sphhar%nlh(nd)
113 : l = sphhar%llh(lh,nd)
114 : sm = (0.0,0.0)
115 : do jm = 1, sphhar%nmem(lh,nd)
116 : m = sphhar%mlh(jm,lh,nd)
117 : lm = l * ( l + 1 ) + m + 1
118 : sm = sm + conjg( sphhar%clnu(jm,lh,nd) ) * pylm(lm,n)
119 : end do
120 : vtlStars(lh,n,iTempArray) = vtlStars(lh,n,iTempArray) + cp * sbf(l) * sm
121 : end do
122 : end do
123 : end do
124 : !$OMP END PARALLEL DO
125 : END DO
126 :
127 : !$OMP parallel do default( NONE ) &
128 : !$OMP SHARED(atoms,sym,sphhar,maxBunchSize,vtlStars,vtlLocal) &
129 696 : !$OMP private(nd,lh,temp,iTempArray)
130 : DO n = 1, atoms%ntype
131 : nd = sym%ntypsy(atoms%firstAtom(n))
132 : do lh = 0, sphhar%nlh(nd)
133 : temp = CMPLX(0.0,0.0)
134 : DO iTempArray = 1, maxBunchSize
135 : temp = temp + vtlStars(lh,n,iTempArray)
136 : END DO
137 : vtlLocal(lh,n) = temp
138 : END DO
139 : END DO
140 : !$OMP END PARALLEL DO
141 :
142 : ! q=0 component
143 2784 : if ( fmpi%irank == 0 .AND. norm2(stars%center)<=1e-8 ) then
144 975 : DO n = 1, atoms%ntype
145 975 : vtlLocal(0,n) = vtlLocal(0,n) + sfp_const * vpw(1)
146 : END DO
147 : end if
148 :
149 696 : CALL mpi_sum_reduce(vtlLocal, vtl,fmpi%mpi_comm)
150 :
151 : ! SPHERE INTERIOR CONTRIBUTION to the coefficients calculated from the
152 : ! values of the sphere Coulomb/Yukawa potential on the sphere boundary
153 :
154 696 : if( fmpi%irank == 0 ) then
155 348 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
156 18 : allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
157 : end if
158 :
159 975 : do n = 1, atoms%ntype
160 627 : nd = sym%ntypsy(atoms%firstAtom(n))
161 627 : imax = atoms%jri(n)
162 627 : lmax = sphhar%llh(sphhar%nlh(nd), nd)
163 627 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
164 : !do concurrent (i = 1:imax)
165 11370 : do i = 1,imax
166 11370 : call ModSphBessel( il(0:,i), kl(0:,i), input%preconditioning_param * atoms%rmsh(i,n), lmax )
167 : end do
168 : end if
169 16806 : do lh = 0, sphhar%nlh(nd)
170 15831 : l = sphhar%llh(lh,nd)
171 15831 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
172 156906 : green_1(1:imax) = il(l,1:imax)
173 156906 : green_2(1:imax) = kl(l,1:imax)
174 207 : green_factor = fpi_const * input%preconditioning_param
175 : else
176 11166770 : green_1(1:imax) = atoms%rmsh(1:imax,n) ** l
177 11166770 : green_2(1:imax) = 1.0 / ( green_1(1:imax) * atoms%rmsh(1:imax,n) )
178 15624 : green_factor = fpi_const / ( 2 * l + 1 )
179 : end if
180 :
181 11323676 : integrand_1(1:imax) = green_1(1:imax) * rho(1:imax,lh,n)
182 11323676 : integrand_2(1:imax) = green_2(1:imax) * rho(1:imax,lh,n)
183 15831 : if (.not.dosf) THEN
184 15669 : call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
185 15669 : call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
186 : else
187 162 : call sfint(integrand_1(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_1(1:imax))
188 162 : call sfint(integrand_2(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_2(1:imax))
189 : end if
190 15831 : termsR = integral_2(imax) + ( vtl(lh,n) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
191 : vr(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
192 11323676 : + green_2(1:imax) * integral_1(1:imax) )
193 16458 : IF (l_dfptvgen) THEN
194 : ! Integrate the imaginary part of the density perturbation as well.
195 0 : integrand_1(1:imax) = green_1(1:imax) * rhoIm(1:imax,lh,n)
196 0 : integrand_2(1:imax) = green_2(1:imax) * rhoIm(1:imax,lh,n)
197 0 : call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
198 0 : call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
199 0 : termsR = integral_2(imax) + ( AIMAG(vtl(lh,n)) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
200 : vrIm(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
201 0 : + green_2(1:imax) * integral_1(1:imax) )
202 : END IF
203 : end do
204 : end do
205 348 : if ( potdenType == POTDEN_TYPE_POTYUK ) then
206 3 : deallocate( il, kl )
207 : end if
208 :
209 348 : if ( potdenType /= POTDEN_TYPE_POTYUK .AND. potdenType /= POTDEN_TYPE_CRYSTALFIELD) then
210 344 : IF (.NOT.l_dfptvgen) THEN
211 951 : do n = 1, atoms%ntype
212 417810 : 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)
213 : end do
214 0 : ELSE IF (.NOT.PRESENT(iDir2)) THEN
215 : ! DFPT case:
216 : ! l=1 contributions from the Coulomb singularity instead of l=0 (1/r -> 1/r^2)
217 0 : DO n = MERGE(1,iDtype,iDtype==0), MERGE(atoms%ntype,iDtype,iDtype==0)
218 0 : ptsym = sym%ntypsy(atoms%firstAtom(n))
219 0 : pref = MERGE(atoms%zatom(n),-atoms%zatom(n),iDtype==0)
220 0 : DO lh = 1, 3
221 0 : l = sphhar%llh(lh, ptsym)
222 0 : DO iMem = 1, sphhar%nmem(lh, ptsym)
223 0 : m = sphhar%mlh(iMem, lh, ptsym)
224 0 : lm = l*(l+1) + m + 1
225 : vr(1:atoms%jri(n),lh,n) = vr(1:atoms%jri(n),lh,n) + &
226 : conjg(sphhar%clnu(iMem, lh, ptsym)) * c_im(iDir, lm - 1) * pref * &
227 0 : ( 1 - (atoms%rmsh(1:atoms%jri(n), n) / atoms%rmt(n))**3) / atoms%rmsh(1:atoms%jri(n),n)**2
228 : END DO
229 : END DO
230 : END DO
231 : ELSE
232 : ! DFPT 2nd order case:
233 : ! l=2 contributions from the Coulomb singularity instead of l=0 (1/r -> 1/r^3)
234 0 : DO n = 1, atoms%ntype!MERGE(1,iDtype,iDtype==0), MERGE(atoms%ntype,iDtype,iDtype==0)
235 0 : ptsym = sym%ntypsy(atoms%firstAtom(n))
236 0 : pref = -atoms%zatom(n)!MERGE(atoms%zatom(n),-atoms%zatom(n),iDtype==0)
237 0 : DO lh = 4, 8
238 0 : l = sphhar%llh(lh, ptsym)
239 0 : DO iMem = 1, sphhar%nmem(lh, ptsym)
240 0 : m = sphhar%mlh(iMem, lh, ptsym)
241 0 : lm = l*(l+1) + m + 1
242 0 : IF ((n.EQ.iDtype).OR.(0.EQ.iDtype)) vr(1:atoms%jri(n),lh,n) = vr(1:atoms%jri(n),lh,n) + &
243 : conjg(sphhar%clnu(iMem, lh, ptsym)) * mat2ord(lm-4,iDir2,iDir) * pref * &
244 0 : ( 1 - (atoms%rmsh(1:atoms%jri(n), n) / atoms%rmt(n))**5) / atoms%rmsh(1:atoms%jri(n),n)**3
245 : END DO
246 : END DO
247 : END DO
248 : END IF
249 : end if
250 : end if
251 :
252 1392 : end subroutine vmts
253 :
254 : end module m_vmts
|
