TABLE OF CONTENTS
- ABINIT/m_ppmodel
- m_ppmodel/calc_sig_ppm
- m_ppmodel/cppm1par
- m_ppmodel/cppm2par
- m_ppmodel/cppm3par
- m_ppmodel/cppm4par
- m_ppmodel/cqratio
- m_ppmodel/get_ppm_eigenvalues
- m_ppmodel/getem1_from_ppm
- m_ppmodel/getem1_from_ppm_one_ggp
- m_ppmodel/new_setup_ppmodel
- m_ppmodel/ppm_free
- m_ppmodel/ppm_get_qbz
- m_ppmodel/ppm_init
- m_ppmodel/ppm_mallocq
- m_ppmodel/ppm_nullify
- m_ppmodel/ppm_symmetrizer
- m_ppmodel/ppm_table_free
- m_ppmodel/ppm_times_ket
- m_ppmodel/ppmodel_t
- m_ppmodel/setup_ppmodel
ABINIT/m_ppmodel [ Modules ]
NAME
m_ppmodel
FUNCTION
Module containing the definition of the ppmodel_t used to deal with the plasmonpole technique. Methods to operate on the object are also provided.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG, GMR, VO, LR, RWG, RS) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_ppmodel 24 25 use defs_basis 26 use m_errors 27 use m_abicore 28 use m_array 29 use m_linalg_interfaces 30 use m_distribfft 31 32 use defs_abitypes, only : MPI_type 33 use m_fstrings, only : sjoin, itoa 34 use m_hide_lapack, only : xhegv 35 use m_gwdefs, only : GW_Q0_DEFAULT, czero_gw 36 use m_crystal, only : crystal_t 37 use m_bz_mesh, only : kmesh_t 38 use m_gsphere, only : gsphere_t 39 use m_vcoul, only : vcoul_t 40 use m_qplusg, only : cmod_qpg 41 use m_fft_mesh, only : g2ifft 42 use m_fft, only : fourdp 43 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 44 45 implicit none 46 47 private
m_ppmodel/calc_sig_ppm [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
calc_sig_ppm
FUNCTION
Calculate the contribution to self-energy operator using a plasmon-pole model.
INPUTS
nomega=Number of frequencies. nspinor=Number of spinorial components. npwc=Number of G vectors in the plasmon pole. npwx=number of G vectors in rhotwgp, i.e. no. of G-vectors for the exchange part. theta_mu_minus_e0i= $\theta(\mu-\epsilon_{k-q,b1,s}), defines if the state is occupied or not zcut=Small imaginary part to avoid the divergence. (see related input variable) omegame0i(nomega)=Frequencies used to evaluate \Sigma_c ($\omega$ - $\epsilon_i)$ otq(npwc,dm2_otq)=Plasmon pole parameters for this q-point. PPm<ppmodel_t>=structure gathering info on the Plasmon-pole technique. %model=type plasmon pole model %dm2_botsq= 1 if model==3, =npwc if model== 4, 1 for all the other cases %dm2_otq= 1 if model==3, =1 if model== 4, 1 for all the other cases %dm_eig=npwc if model=3, 0 otherwise botsq(npwc,dm2_botsq)=Plasmon pole parameters for this q-point. eig(dm_eig,dm_eig)=The eigvectors of the symmetrized inverse dielectric matrix for this q point (first index for G, second index for bands) rhotwgp(npwx)=oscillator matrix elements divided by |q+G| i.e. $\frac{\langle b1 k-q s | e^{-i(q+G)r | b2 k s \rangle}{|q+G|}$
OUTPUT
sigcme(nomega) (to be described), only relevant if ppm3 or ppm4
ket(npwc,nomega):
=== model==1,2 ====
ket(G,omega) = Sum_G2 conjg(rhotw(G)) * Omega(G,G2) * rhotw(G2)
---------------------------------------------------
2 omegatw(G,G2) (omega-E_i + omegatw(G,G2)(2f-1))
SOURCE
2106 subroutine calc_sig_ppm(PPm,nspinor,npwc,nomega,rhotwgp,botsq,otq,& 2107 & omegame0i,zcut,theta_mu_minus_e0i,eig,npwx,ket,sigcme) 2108 2109 !Arguments ------------------------------------ 2110 !scalars 2111 integer,intent(in) :: nomega,npwc,npwx,nspinor 2112 real(dp),intent(in) :: theta_mu_minus_e0i,zcut 2113 type(ppmodel_t),intent(in) :: PPm 2114 !arrays 2115 real(dp),intent(in) :: omegame0i(nomega) 2116 complex(gwpc),intent(in) :: botsq(npwc,PPm%dm2_botsq),eig(PPm%dm_eig,PPm%dm_eig),otq(npwc,PPm%dm2_otq) 2117 complex(gwpc),intent(in) :: rhotwgp(npwx*nspinor) 2118 complex(gwpc),intent(inout) :: ket(npwc*nspinor,nomega) 2119 complex(gwpc),intent(out) :: sigcme(nomega) 2120 2121 !Local variables------------------------------- 2122 !scalars 2123 integer :: ig,igp,ii,ios,ispinor,spadc,spadx 2124 real(dp) :: den,ff,inv_den,omegame0i_io,otw,twofm1,twofm1_zcut 2125 complex(gwpc) :: ct,num,numf,rhotwgdp_igp 2126 logical :: fully_occupied,totally_empty 2127 !character(len=500) :: msg 2128 !arrays 2129 complex(gwpc),allocatable :: rhotwgdpcc(:) 2130 2131 !************************************************************************* 2132 2133 SELECT CASE (PPm%model) 2134 CASE (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE) 2135 fully_occupied =(ABS(theta_mu_minus_e0i-one)<0.001) 2136 totally_empty =(ABS(theta_mu_minus_e0i )<0.001) 2137 2138 do ispinor=1,nspinor 2139 spadx=(ispinor-1)*npwx; spadc=(ispinor-1)*npwc 2140 2141 if (.not.totally_empty) then 2142 ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(\mu-e_s) / (\omega+\omegat_{G1G2}-e_s-i\delta) 2143 twofm1_zcut=zcut 2144 !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den) 2145 do ios=1,nomega 2146 omegame0i_io=omegame0i(ios) 2147 do igp=1,npwc 2148 rhotwgdp_igp=rhotwgp(spadx+igp) 2149 do ig=1,npwc 2150 otw=DBLE(otq(ig,igp)) !in principle otw -> otw - ieta 2151 num = botsq(ig,igp)*rhotwgdp_igp 2152 den = omegame0i_io+otw 2153 if (den**2>zcut**2) then 2154 ket(spadc+ig,ios)=ket(spadc+ig,ios) + num/(den*otw)*theta_mu_minus_e0i 2155 else 2156 ket(spadc+ig,ios)=ket(spadc+ig,ios) + & 2157 & num*CMPLX(den,twofm1_zcut)/((den**2+twofm1_zcut**2)*otw)*theta_mu_minus_e0i 2158 end if 2159 end do !ig 2160 end do !igp 2161 end do !ios 2162 end if !not totally empty 2163 2164 if (.not.(fully_occupied)) then 2165 ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(e_s-\mu) / (\omega-\omegat_{G1G2}-e_s+i\delta) 2166 twofm1_zcut=-zcut 2167 !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den) 2168 do ios=1,nomega 2169 omegame0i_io=omegame0i(ios) 2170 do igp=1,npwc 2171 rhotwgdp_igp=rhotwgp(spadx+igp) 2172 do ig=1,npwc 2173 otw=DBLE(otq(ig,igp)) !in principle otw -> otw - ieta 2174 num = botsq(ig,igp)*rhotwgdp_igp 2175 den=omegame0i_io-otw 2176 if (den**2>zcut**2) then 2177 ket(spadc+ig,ios)=ket(spadc+ig,ios) + num/(den*otw)*(one-theta_mu_minus_e0i) 2178 else 2179 ket(spadc+ig,ios)=ket(spadc+ig,ios) + & 2180 & num*CMPLX(den,twofm1_zcut)/((den**2+twofm1_zcut**2)*otw)*(one-theta_mu_minus_e0i) 2181 end if 2182 end do !ig 2183 end do !igp 2184 ! 2185 end do !ios 2186 end if !not fully occupied 2187 2188 end do !ispinor 2189 2190 ket=ket*half 2191 2192 CASE (PPM_LINDEN_HORSH,PPM_ENGEL_FARID) 2193 ABI_CHECK(nspinor == 1, "nspinor/=1 not allowed") 2194 2195 ! rho-twiddle(G) is formed, introduce rhotwgdpcc, for speed reason 2196 ABI_MALLOC(rhotwgdpcc,(npwx)) 2197 2198 ff=theta_mu_minus_e0i ! occupation number f (include poles if ...) 2199 twofm1=two*ff-one ! 2f-1 2200 twofm1_zcut=twofm1*zcut 2201 rhotwgdpcc(:)=CONJG(rhotwgp(:)) 2202 2203 do ios=1,nomega 2204 omegame0i_io=omegame0i(ios) 2205 ct=czero_gw 2206 do ii=1,npwc ! Loop over the DM bands 2207 num=czero_gw 2208 2209 SELECT CASE (PPm%model) 2210 CASE (PPM_LINDEN_HORSH) 2211 ! Calculate \beta (eq. 106 pag 47) 2212 do ig=1,npwc 2213 num=num+rhotwgdpcc(ig)*eig(ig,ii) 2214 end do 2215 numf=num*CONJG(num) !MG this means that we cannot do SCGW 2216 numf=numf*botsq(ii,1) 2217 2218 CASE (PPM_ENGEL_FARID) 2219 do ig=1,npwc 2220 num=num+rhotwgdpcc(ig)*botsq(ig,ii) 2221 end do 2222 numf=num*CONJG(num) !MG this means that we cannot do SCGW 2223 2224 CASE DEFAULT 2225 ABI_ERROR("Wrong PPm%model") 2226 END SELECT 2227 2228 !numf=num*CONJG(num) !MG this means that we cannot do SCGW 2229 !if (PPm%model==3) numf=numf*botsq(ii,1) 2230 2231 otw=DBLE(otq(ii,1)) ! in principle otw -> otw - ieta 2232 den=omegame0i_io+otw*twofm1 2233 2234 if (den**2>zcut**2) then 2235 inv_den=one/den 2236 ct=ct+numf*inv_den 2237 else 2238 inv_den=one/((den**2+twofm1_zcut**2)) 2239 ct=ct+numf*CMPLX(den,twofm1_zcut)*inv_den 2240 end if 2241 2242 end do !ii DM bands 2243 sigcme(ios)=ct*half 2244 end do !ios 2245 ABI_FREE(rhotwgdpcc) 2246 2247 CASE DEFAULT 2248 ABI_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model))) 2249 END SELECT 2250 2251 end subroutine calc_sig_ppm
m_ppmodel/cppm1par [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
cppm1par
FUNCTION
Calculate the plasmon-pole parameters big-omega-twiddle-squared and omega-twiddle from epsilon-twiddle^-1 calculated for nomega (usually 2) frequencies omega=0 and omega=iE0.
INPUTS
epsm1(npwc,npwc,nomega)=dielectric matrix at nomega frequencies. npwc=number of plane waves nomega=number of frequencies (usually 2) omega(nomega)=frequencies omegaplasma=input variable or Drude plasma frequency
OUTPUT
bigomegatwsq(npwc,npwc)=parameter of the plasmon-pole model (see gwa.pdf file) omegatw(npwc,npwc)=parameter of the plasmon-pole model (see gwa.pdf file)
TODO
Calculation can be done in place.
SOURCE
1171 subroutine cppm1par(npwc,nomega,omega,omegaplasma,epsm1,omegatw,bigomegatwsq) 1172 1173 !Arguments ------------------------------------ 1174 !scalars 1175 integer,intent(in) :: nomega,npwc 1176 real(dp),intent(in) :: omegaplasma 1177 !arrays 1178 complex(gwpc),intent(in) :: epsm1(npwc,npwc,nomega) 1179 complex(dpc),intent(in) :: omega(nomega) 1180 complex(gwpc),intent(out) :: bigomegatwsq(npwc,npwc) 1181 complex(gwpc),intent(out) :: omegatw(npwc,npwc) 1182 1183 !Local variables------------------------------- 1184 !scalars 1185 integer :: ig,igp,io,io0,ioe0 1186 real(dp) :: e0,minomega 1187 character(len=500) :: msg 1188 complex(gwpc) :: AA,omegatwsq,diff,ratio 1189 complex(gwpc) :: epsm1_io0,epsm1_ioe0 1190 1191 ! ************************************************************************* 1192 1193 DBG_ENTER("COLL") 1194 ! 1195 ! === Find omega=0 and omega=imag (closest to omegaplasma) to fit the ppm parameters === 1196 minomega=1.0d-3; io0=0 1197 do io=1,nomega 1198 if (ABS(omega(io))<minomega) then 1199 io0=io; minomega=ABS(omega(io)) 1200 end if 1201 end do 1202 ABI_CHECK(io0/=0,"omega=0 not found") 1203 1204 minomega=1.0d-3; e0=200.0; ioe0=0 1205 do io=1,nomega 1206 if (REAL(omega(io))<minomega.and.AIMAG(omega(io))>minomega) then 1207 if (ABS(AIMAG(omega(io))-omegaplasma)<ABS(e0-omegaplasma)) then 1208 ioe0=io; e0=AIMAG(omega(io)) 1209 end if 1210 end if 1211 end do 1212 1213 write(msg,'(a,f9.4,a)')' Imaginary frequency for fit located at: ',e0*Ha_eV,' [eV] ' 1214 call wrtout(std_out,msg,'COLL') 1215 ABI_CHECK(ioe0/=0,"Imaginary omega not found") 1216 ! 1217 ! ================================================================ 1218 ! === Calculate plasmon-pole A parameter A=epsilon^-1(0)-delta === 1219 ! ================================================================ 1220 do ig=1,npwc 1221 do igp=1,npwc 1222 epsm1_io0 = epsm1(ig,igp,io0) 1223 epsm1_ioe0 = epsm1(ig,igp,ioe0) 1224 1225 AA=epsm1_io0 1226 if (ig==igp) AA=AA-one 1227 1228 ! === Calculate plasmon-pole omega-twiddle-square parameter === 1229 ! XG201009 Strangely, the next formula does not work with gcc43-debug 1230 ! omegatwsq=(AA/(epsm1_io0-epsm1_ioe0)-one)*e0**2 1231 ! This seems to be due to precision issue at the level of division by a complex whose norm squared 1232 ! is below the smallest representable number. 1233 ! After many trials, I have decided to shift the difference by a small number ... well, not so small ... 1234 ! for numerical issues 1235 diff=epsm1_io0-epsm1_ioe0 1236 diff=diff+cmplx(tol10,tol10) 1237 ratio=AA/diff 1238 omegatwsq=(ratio-cone)*e0**2 1239 ! 1240 ! If omega-twiddle-squared is negative,set omega-twiddle-squared to 1.0 (a reasonable way of treating 1241 ! such terms, in which epsilon**-1 was originally increasing along this part of the imaginary axis) 1242 ! (note: originally these terms were ignored in Sigma; this was changed on 6 March 1990.) 1243 1244 if (REAL(omegatwsq)<=zero) omegatwsq=one 1245 ! 1246 ! Get omega-twiddle 1247 ! * Neglect the imag part (if any) in omega-twiddle-squared 1248 omegatw(ig,igp)=SQRT(REAL(omegatwsq)) 1249 ! 1250 ! Get big-omega-twiddle-squared=-omega-twiddle-squared AA 1251 bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2 1252 1253 end do !igp 1254 end do !ig 1255 1256 write(msg,'(2a,f15.12,2a,2i5,a)')ch10,& 1257 & ' cppm1par : omega twiddle minval [eV] = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,& 1258 & ' omega twiddle min location = ',MINLOC(ABS(omegatw)),ch10 1259 call wrtout(std_out,msg,'COLL') 1260 1261 DBG_EXIT("COLL") 1262 1263 end subroutine cppm1par
m_ppmodel/cppm2par [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
cppm2par
FUNCTION
Calculate plasmon-pole parameters of the Hybertsen and Louie model (PRB 34, 5390 (1986) [[cite:Hybertsen1986]])
INPUTS
qpt(3)=The coordinates of the q-point in the IBZ. epsm1(npwc,npwc)=symmetrized inverse dielectric (static limit is used) gmet(3,3)=metric in reciprocal space ngfftf(18)=contain all needed information about the 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft npwc=number of plane waves in epsm1 rhor(nfftf)=charge density on the real space FFT grid nfftf= total number of points in the fine FFT mesh (for this processor) invalid_freq: what to do when PPM omega is found negative or imaginary 0) drop it (default as specified in Hybersen-Louie original GW paper) 1) set to 1 hartree 2) set to infinity
OUTPUT
bigomegatwsq(npwc,npwc)= squared bare plasma frequencies
\Omega^2_{G1 G2}(q) = 4\pi \frac {(q+G1).(q+G2)}/{|q+G1|^2} n(G1-G2)
omegatw(npwc,npwc)= plasmon frequencies \tilde\omega_{G1 G2}(q) where:
\tilde\omega^2_{G1 G2}(q) =
\frac {\Omega^2_{G1 G2}(q)} {\delta_{G1 G2}-\tilde\epsilon^{-1}_{G1 G2} (q, \omega=0)}
SOURCE
1297 subroutine cppm2par(qpt,npwc,epsm1,ngfftf,gvec,gprimd,rhor,nfftf,gmet,bigomegatwsq,omegatw,invalid_freq) 1298 1299 !Arguments ------------------------------------ 1300 !scalars 1301 integer,intent(in) :: npwc,nfftf,invalid_freq 1302 !arrays 1303 integer,intent(in) :: gvec(3,npwc) 1304 integer,intent(in) :: ngfftf(18) 1305 real(dp),intent(in) :: qpt(3),gmet(3,3),gprimd(3,3) 1306 real(dp),intent(in) :: rhor(nfftf) 1307 complex(gwpc),intent(in) :: epsm1(npwc,npwc) 1308 complex(gwpc),intent(out) :: bigomegatwsq(npwc,npwc) 1309 complex(gwpc),intent(out) :: omegatw(npwc,npwc) 1310 1311 !Local variables------------------------------- 1312 !scalars 1313 integer,parameter :: paral_kgb0=0 1314 integer :: ig,igp,nimwp,ngfft1,ngfft2,ngfft3,gmgp_idx,ierr 1315 real(dp) :: lambda,phi,AA 1316 logical,parameter :: use_symmetrized=.TRUE.,check_imppf=.FALSE. 1317 character(len=500) :: msg 1318 type(MPI_type) :: MPI_enreg_seq 1319 !arrays 1320 real(dp) :: qlist(3,1) 1321 real(dp),allocatable :: tmp_rhor(:),qratio(:,:),qplusg(:),rhog_dp(:,:) 1322 complex(gwpc),allocatable :: omegatwsq(:,:) 1323 complex(gwpc),allocatable :: rhog(:),rhogg(:,:),temp(:,:) !MG these should be double precision TODO 1324 !************************************************************************* 1325 1326 call initmpi_seq(MPI_enreg_seq) 1327 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftf(2),ngfftf(3),'all') 1328 ! 1329 ! === Calculate qratio(npwec,npvec) = (q+G).(q+Gp)/|q+G|^2 === 1330 ABI_MALLOC_OR_DIE(qratio,(npwc,npwc), ierr) 1331 1332 call cqratio(npwc,gvec,qpt,gmet,gprimd,qratio) 1333 ! 1334 ! === Compute the density in G space rhor(R)--> rhog(G) === 1335 ABI_MALLOC(rhog_dp,(2,nfftf)) 1336 ABI_MALLOC(rhog,(nfftf)) 1337 ngfft1=ngfftf(1) 1338 ngfft2=ngfftf(2) 1339 ngfft3=ngfftf(3) 1340 1341 ABI_MALLOC(tmp_rhor,(nfftf)) 1342 tmp_rhor=rhor ! To avoid having to use intent(inout). 1343 call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,1,ngfftf,0) 1344 ABI_FREE(tmp_rhor) 1345 1346 rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf)) 1347 ! 1348 ! Calculate the FFT index of each (G-Gp) vector and assign 1349 ! the value of the correspondent density simultaneously 1350 ABI_MALLOC_OR_DIE(rhogg,(npwc,npwc), ierr) 1351 1352 ierr=0 1353 do ig=1,npwc 1354 do igp=1,npwc 1355 gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf) 1356 if (gmgp_idx/=0) then 1357 rhogg(ig,igp)=rhog(gmgp_idx) 1358 else 1359 ierr=ierr+1 1360 rhogg(ig,igp)=czero 1361 end if 1362 end do 1363 end do 1364 1365 if (ierr/=0) then 1366 write(msg,'(a,i4,3a)')& 1367 & 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,& 1368 & 'Enlarge the FFT mesh to get rid of this problem. ' 1369 ABI_WARNING(msg) 1370 end if 1371 1372 rhogg=four_pi*rhogg 1373 ABI_FREE(rhog_dp) 1374 ABI_FREE(rhog) 1375 ! 1376 ! Calculate GPP parameters 1377 ! unsymmetrized epsm1 -> epsm1=|q+Gp|/|q+G|*epsm1 1378 ABI_MALLOC(qplusg,(npwc)) 1379 ABI_MALLOC(temp,(npwc,npwc)) 1380 ABI_MALLOC_OR_DIE(omegatwsq,(npwc,npwc), ierr) 1381 1382 temp = -epsm1(:,:) 1383 ! 1384 ! RS still not obvious for me whether one shall use the symmetrized inverse DM or the unsymmetrized one 1385 ! the default here is to use the symmetrized one, I must discuss this with XG 1386 ! 1387 ! MG it turns out that using the symmetrized inverse DM in the plasmon-pole 1388 ! equations give the same results for the squared plasmon frequencies omegatwsq while the 1389 ! squared bare plasma frequencies bigomegatwsq related to the symmetrized dielectric matrix 1390 ! are obtained multipling by |q+G1|/|q+G2| 1391 ! 1392 if (.not.use_symmetrized) then 1393 qlist(:,1) = qpt 1394 call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q 1395 do ig=1,npwc 1396 do igp=1,npwc 1397 temp(ig,igp)=qplusg(igp)/qplusg(ig)*temp(ig,igp) 1398 end do 1399 end do 1400 end if 1401 1402 nimwp=0 1403 do ig=1,npwc 1404 temp(ig,ig)=temp(ig,ig)+one 1405 do igp=1,npwc 1406 bigomegatwsq(ig,igp) = rhogg(ig,igp)*qratio(ig,igp) 1407 omegatwsq(ig,igp)=bigomegatwsq(ig,igp)/temp(ig,igp) 1408 ! 1409 ! Set to an arbitrary value the omegawsq which become negative or imaginary 1410 ! in principle these correspond to cases where the imaginary part of epsm1 does not have 1411 ! a well defined peak. The imaginary part of epsm1 in these cases oscillates with a small amplitude 1412 ! since the amplitude A_GGpr=-pi/2*bigomegatwsq/omegatw, 1413 ! it follows that bigomegatwsq shall be set to zero for these cases 1414 if ( REAL(omegatwsq(ig,igp))<= tol12 .or. AIMAG(omegatwsq(ig,igp))**2*tol12> REAL(omegatwsq(ig,igp))**2) then 1415 nimwp=nimwp+1 1416 1417 if ( invalid_freq == 1 ) then 1418 ! set omegatwsq to 1 hartree 1419 omegatwsq(ig,igp)=cone 1420 AA = epsm1(ig,igp) 1421 if ( ig == igp ) AA = AA - one 1422 omegatw(ig,igp)=SQRT(REAL(omegatwsq(ig,igp))) 1423 bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2 1424 elseif ( invalid_freq == 2 ) then 1425 ! set omegatwsq to infinity 1426 omegatwsq(ig,igp)=cone/tol6 1427 AA = epsm1(ig,igp) 1428 if ( ig == igp ) AA = AA - one 1429 omegatw(ig,igp)=SQRT(REAL(omegatwsq(ig,igp))) 1430 bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2 1431 else 1432 ! simply ignore all cases of omegatw with imaginary values 1433 bigomegatwsq(ig,igp)=(0.,0.) 1434 omegatw(ig,igp)=(ten,0.) 1435 end if 1436 if (check_imppf) then 1437 write(msg,'(a,2i8)')' Imaginary plasmon frequency at : ',ig,igp 1438 call wrtout(std_out,msg,'COLL') 1439 end if 1440 else 1441 ! this part has been added to deal with systems without inversion symmetry 1442 ! this new implementation gives the same results as the previous one if 1443 ! omegatwsq is a pure real number and has the advantage of being an improved 1444 ! approach for systems without an inversion center. 1445 lambda=ABS(omegatwsq(ig,igp)) 1446 phi=ATAN(AIMAG(omegatwsq(ig,igp))/REAL(omegatwsq(ig,igp))) 1447 omegatw(ig,igp)=SQRT(lambda/COS(phi)) 1448 bigomegatwsq(ig,igp)=bigomegatwsq(ig,igp)*(1.-(0.,1.)*TAN(phi)) 1449 ! Uncomment the following line and comment the previous to restore the old version. 1450 !omegatw(ig,igp)=sqrt(real(omegatwsq(ig,igp))) 1451 end if 1452 end do 1453 end do 1454 1455 write(msg,'(a,3f12.6,a,i8,a,i8)')& 1456 & ' at q-point : ',qpt,& 1457 & ' number of imaginary plasmonpole frequencies = ',nimwp,' / ',npwc**2 1458 call wrtout(std_out,msg,'COLL') 1459 1460 write(msg,'(2a,f12.8,2a,3i5)')ch10,& 1461 & ' cppm2par : omega twiddle minval [eV] = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,& 1462 & ' omega twiddle min location = ',MINLOC(ABS(omegatw)) 1463 call wrtout(std_out,msg,'COLL') 1464 1465 call destroy_mpi_enreg(MPI_enreg_seq) 1466 1467 ABI_FREE(omegatwsq) 1468 ABI_FREE(rhogg) 1469 ABI_FREE(temp) 1470 ABI_FREE(qplusg) 1471 ABI_FREE(qratio) 1472 1473 end subroutine cppm2par
m_ppmodel/cppm3par [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
cppm3par
FUNCTION
Calculate the plasmon-pole parameters using the von Linden-Horsh model (PRB 37, 8351, 1988) [[cite:vonderLinden1988]] (see also Pag 22 of Quasiparticle Calculations in Solids [[cite:Aulbur2001]].
INPUTS
epsm1(npwc,npwc))= symmetrized inverse dielectric ngfftf(18)=contain all needed information about 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft npwc=number of plane waves in epsm1 qratio=(q+G1).(q+G2)/(|q+G1|.|q+G2|) rhor(nfftf)=charge density on the real space FFT grid nfftf=number of points in the FFT grid (for this processor) gvec(3,npwc)= G vectors in reduced coordinates
OUTPUT
omegatw(npwc,npwc)= plasmon pole positions
bigomegatwsq(npwc,npwc)=(E_{q,ii}^{-1}-1)*omegatw
where E^{-1} is the eigenvalue of the inverse dielectric matrix
eigtot(npwc,npwc)=the eigvectors of the symmetrized inverse dielectric matrix
(first index for G, second index for bands)
SOURCE
1504 subroutine cppm3par(qpt,npwc,epsm1,ngfftf,gvec,gprimd,rhor,nfftf,bigomegatwsq,omegatw,eigtot) 1505 1506 !Arguments ------------------------------------ 1507 !scalars 1508 integer,intent(in) :: nfftf,npwc 1509 !arrays 1510 integer,intent(in) :: gvec(3,npwc),ngfftf(18) 1511 real(dp),intent(in) :: qpt(3),gprimd(3,3),rhor(nfftf) 1512 complex(gwpc),intent(in) :: epsm1(npwc,npwc) 1513 complex(gwpc),intent(out) :: bigomegatwsq(npwc,1),eigtot(npwc,npwc) 1514 complex(gwpc),intent(out) :: omegatw(npwc) 1515 1516 !Local variables------------------------------- 1517 !TODO these should be dp 1518 !scalars 1519 integer,parameter :: paral_kgb0=0 1520 integer :: idx,ierr,ig,igp,ii,jj,ngfft1,ngfft2,ngfft3,gmgp_idx 1521 real(dp) :: num,qpg_dot_qpgp 1522 complex(dpc) :: conjg_eig 1523 logical :: qiszero 1524 character(len=500) :: msg 1525 type(MPI_type) :: MPI_enreg_seq 1526 !arrays 1527 real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),qlist(3,1) 1528 real(dp),allocatable :: eigval(:),qplusg(:),rhog_dp(:,:),zhpev2(:),tmp_rhor(:) 1529 complex(dpc),allocatable :: eigvec(:,:),matr(:),mm(:,:),rhog(:),rhogg(:,:) 1530 complex(dpc),allocatable :: zhpev1(:),zz(:) 1531 1532 !************************************************************************* 1533 1534 ! Fake MPI_type for the sequential part. 1535 call initmpi_seq(MPI_enreg_seq) 1536 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftf(2),ngfftf(3),'all') 1537 1538 qiszero = (ALL(ABS(qpt)<1.0e-3)) 1539 1540 b1=two_pi*gprimd(:,1) 1541 b2=two_pi*gprimd(:,2) 1542 b3=two_pi*gprimd(:,3) 1543 1544 ngfft1=ngfftf(1) 1545 ngfft2=ngfftf(2) 1546 ngfft3=ngfftf(3) 1547 1548 ABI_MALLOC(rhog_dp,(2,nfftf)) 1549 ABI_MALLOC(rhog,(nfftf)) 1550 1551 ABI_MALLOC_OR_DIE(rhogg,(npwc,npwc), ierr) 1552 ! 1553 ! === Compute the density in G space rhog(r)--> rho(G) === 1554 ! FIXME this has to be fixed, rho(G) should be passed instead of doing FFT for each q 1555 1556 ABI_MALLOC(tmp_rhor,(nfftf)) 1557 tmp_rhor=rhor ! To avoid having to use intent(inout). 1558 call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,1,ngfftf,0) 1559 ABI_FREE(tmp_rhor) 1560 1561 rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf)) 1562 ! 1563 ! Calculate the FFT index of each (G-Gp) vector and assign the value 1564 ! of the correspondent density simultaneously 1565 ierr=0 1566 do ig=1,npwc 1567 do igp=1,npwc 1568 gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf) 1569 if (gmgp_idx/=0) then 1570 rhogg(ig,igp)=rhog(gmgp_idx) 1571 else 1572 ierr=ierr+1 1573 rhogg(ig,igp)=czero 1574 end if 1575 end do 1576 end do 1577 1578 if (ierr/=0) then 1579 write(msg,'(a,i4,3a)')& 1580 & 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,& 1581 & 'Enlarge the FFT mesh to get rid of this problem. ' 1582 ABI_WARNING(msg) 1583 end if 1584 ! 1585 ! mm(G,Gp) = (q+G) \cdot (q+Gp) n(G-Gp) 1586 ABI_MALLOC_OR_DIE(mm,(npwc,npwc), ierr) 1587 1588 do ig=1,npwc 1589 if (qiszero) then 1590 ! To be discussed with Riad, here we should use the small q 1591 ! to be consistent and consider the limit q-->0 1592 gpq(:)=gvec(:,ig) 1593 else 1594 gpq(:)=gvec(:,ig)+qpt 1595 end if 1596 do igp=1,npwc 1597 if (qiszero) then 1598 gppq(:)=gvec(:,igp) 1599 else 1600 gppq(:)=gvec(:,igp)+qpt 1601 end if 1602 qpg_dot_qpgp=zero 1603 do ii=1,3 1604 qpg_dot_qpgp=qpg_dot_qpgp+& 1605 & ( gpq(1)*b1(ii) +gpq(2)*b2(ii) +gpq(3)*b3(ii))*& 1606 & (gppq(1)*b1(ii)+gppq(2)*b2(ii)+gppq(3)*b3(ii)) 1607 end do 1608 mm(ig,igp)=rhogg(ig,igp)*qpg_dot_qpgp 1609 end do !igp 1610 end do !ig 1611 1612 ABI_FREE(rhog_dp) 1613 ABI_FREE(rhog) 1614 ! === Now we have rhogg,rho0 === 1615 ! 1616 ! Calculate the dielectric matrix eigenvalues and vectors 1617 ! Use only the static epsm1 i.e., only the w=0 part (eps(:,:,1,:)) 1618 ABI_MALLOC(eigval,(npwc)) 1619 1620 ABI_MALLOC_OR_DIE(eigvec,(npwc,npwc), ierr) 1621 1622 ABI_MALLOC(zz,(npwc)) 1623 zz=czero 1624 1625 ABI_MALLOC(qplusg,(npwc)) 1626 1627 ! Store the susceptibility matrix in upper mode before calling zhpev. 1628 ABI_MALLOC_OR_DIE(matr,(npwc*(npwc+1)/2), ierr) 1629 1630 idx=1 1631 do ii=1,npwc 1632 do jj=1,ii 1633 matr(idx)=epsm1(jj,ii); idx=idx+1 1634 end do 1635 end do 1636 1637 ABI_MALLOC(zhpev2,(3*npwc-2)) 1638 ABI_MALLOC(zhpev1,(2*npwc-1)) 1639 1640 call ZHPEV('V','U',npwc,matr,eigval,eigvec,npwc,zhpev1,zhpev2,ierr) 1641 ABI_FREE(matr) 1642 ABI_FREE(zhpev2) 1643 ABI_FREE(zhpev1) 1644 1645 if (ierr<0) then 1646 write (msg,'(2a,i4,a)')& 1647 & ' Failed to calculate the eigenvalues and eigenvectors of the dielectric matrix ',ch10,& 1648 & ierr*(-1),'-th argument in the matrix has an illegal value. ' 1649 ABI_ERROR(msg) 1650 end if 1651 1652 if (ierr>0) then 1653 write(msg,'(3a,i4,2a)')& 1654 & ' Failed to calculate the eigenvalues and eigenvectors of the dielectric matrix ',ch10,& 1655 & ' the algorithm failed to converge; ierr = ', ierr,ch10,& 1656 & ' off-diagonal elements of an intermediate tridiagonal form did not converge to zero. ' 1657 ABI_ERROR(msg) 1658 end if 1659 ! 1660 ! Calculate the PPM parameters and the eigenpotentials needed for 1661 ! the calculation of the generalized overlap matrix 1662 ! Note: the eigenpotentials has to be calculated on the FFT (G-Gp) index 1663 ! 1664 ! Save eigenvectors of \tilde\epsilon^{-1} 1665 ! MG well it is better to save \Theta otherwise 1666 ! we have to calculare \Theta for each band, spin, k-point but oh well 1667 eigtot=eigvec 1668 1669 qlist(:,1) = qpt 1670 call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q 1671 ! 1672 ! Basic Equation: 1673 ! 1674 ! \Theta_{q,ii}(G)=\Psi_{q,ii}(G)/|q+G| 1675 ! where \Psi_{q,ii}(G) is the eigenvector of \tilde\epsilon^{-1} 1676 1677 ! \tilde\omega_{ii,q}^2= 4\pi (1-eigenval(ii,q))) 1678 ! \sum_{G,Gp} \Theta^*_{q,ii}(G) (q+G)\cdot(q+Gp) n(G-Gp) \Theta_{q,ii}(Gp) 1679 1680 do ii=1,npwc !DM band 1681 ! Calculate \Theta_{q,ii}(G) 1682 ! why the first element is not modified? if the problem is the small value of qplusg(1) 1683 ! we could multiply by sqrt(mod((q+G)(q+G'))) and then add the sing at the end 1684 if (qiszero)then 1685 eigvec(2:,ii)=eigvec(2:,ii)/qplusg(2:) 1686 else 1687 eigvec(:,ii)=eigvec(:,ii)/qplusg(:) 1688 end if 1689 do ig=1,npwc 1690 conjg_eig=CONJG(eigvec(ig,ii)) 1691 do igp=1,npwc 1692 if(qiszero .and. ig==1 .and. igp==1)then 1693 zz(ii)=zz(ii)+conjg_eig*rhogg(ig,igp)*eigvec(igp,ii) 1694 else 1695 zz(ii)=zz(ii)+conjg_eig*mm(ig,igp)*eigvec(igp,ii) 1696 end if 1697 end do 1698 end do 1699 1700 num=one-eigval(ii) 1701 if (num<=zero) then 1702 ! here I think we should set bigomegatwsq=0 and omegatw to an arbitrary value 1703 ! maybe we can output a warning TO BE discussed with Riad 1704 if (ABS(num)<1.0d-4) then 1705 num=1.0d-5 1706 else 1707 ABI_ERROR("One or more imaginary plasmon pole energies") 1708 end if 1709 end if 1710 1711 omegatw(ii)=SQRT(4*pi*REAL(zz(ii))/num) 1712 ! this should be \alpha = 2\pi omegatw * (1-eigenval) 1713 ! MG check this, in the review I found a factor 2\pi, maybe it is reintroduced later 1714 bigomegatwsq(ii,1)=num*omegatw(ii) 1715 end do 1716 1717 ABI_FREE(rhogg) 1718 ABI_FREE(mm) 1719 ABI_FREE(eigval) 1720 ABI_FREE(zz) 1721 ABI_FREE(eigvec) 1722 ABI_FREE(qplusg) 1723 1724 call destroy_mpi_enreg(MPI_enreg_seq) 1725 1726 write(msg,'(2a,f12.8,2a,3i5)')ch10,& 1727 & ' cppm3par : omega twiddle minval [eV] = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,& 1728 & ' omega twiddle min location = ',MINLOC(ABS(omegatw)) 1729 call wrtout(std_out,msg,'COLL') 1730 1731 end subroutine cppm3par
m_ppmodel/cppm4par [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
cppm4par
FUNCTION
Calculate the plasmon-pole parameters using Engel and Farid model (PRB47,15931,1993) [[cite:Engel1993]]. See also Quasiparticle Calculations in Solids [[cite:Aulbur2001]] p. 23
INPUTS
qpt(3)=Reduced coordinates of the q-point. epsm1(npwc,npwc)=symmetrized inverse dielectric matrix. gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) ngfftf(18)=contain all needed information about 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft npwc=number of plane waves in epsm1 rhor(nfftf)=charge density on the real space FFT grid gvec(3,npwc)=G vectors in reduced coordinated
OUTPUT
bigomegatwsq(npwc,npwc)=plasmon-pole strength omegatw(npwc)=plasmon-pole frequencies
SOURCE
1759 subroutine cppm4par(qpt,npwc,epsm1,ngfftf,gvec,gprimd,rhor,nfftf,bigomegatwsq,omegatw) 1760 1761 !Arguments ------------------------------------ 1762 !scalars 1763 integer,intent(in) :: nfftf,npwc 1764 !arrays 1765 integer,intent(in) :: gvec(3,npwc),ngfftf(18) 1766 real(dp),intent(in) :: gprimd(3,3),qpt(3) 1767 real(dp),intent(in) :: rhor(nfftf) 1768 complex(gwpc),intent(in) :: epsm1(npwc,npwc) 1769 complex(gwpc),intent(out) :: bigomegatwsq(npwc,npwc),omegatw(npwc) 1770 1771 !Local variables------------------------------- 1772 !scalars 1773 integer,parameter :: paral_kgb0=0 1774 integer :: ierr,ig,igp,ii,ngfft1,ngfft2,ngfft3,gmgp_idx 1775 real(dp) :: qpg_dot_qpgp 1776 character(len=500) :: msg 1777 character(len=80) :: bar 1778 type(MPI_type) :: MPI_enreg_seq 1779 !arrays 1780 real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),qlist(3,1) 1781 real(dp),allocatable :: eigval(:),qplusg(:),rhog_dp(:,:) 1782 real(dp),allocatable :: tmp_rhor(:) 1783 complex(dpc),allocatable :: chi(:,:),chitmp(:,:),chitmps(:,:),eigvec(:,:) 1784 complex(dpc),allocatable :: mm(:,:),mtemp(:,:),rhog(:) 1785 complex(dpc),allocatable :: rhogg(:,:),tmp1(:),zz2(:,:) 1786 1787 !************************************************************************* 1788 1789 DBG_ENTER("COLL") 1790 1791 call initmpi_seq(MPI_enreg_seq) 1792 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftf(2),ngfftf(3),'all') 1793 1794 b1=two_pi*gprimd(:,1) 1795 b2=two_pi*gprimd(:,2) 1796 b3=two_pi*gprimd(:,3) 1797 ! 1798 ! === Calculate density in G space rhog(G) === 1799 ABI_MALLOC(rhog_dp,(2,nfftf)) 1800 ABI_MALLOC(rhog,(nfftf)) 1801 1802 ABI_MALLOC_OR_DIE(rhogg,(npwc,npwc), ierr) 1803 ! 1804 ! Conduct FFT tho(r)-->rhog(G) 1805 ! FIXME this has to be fixed, rho(G) should be passed instead of doing FFT for each q 1806 ABI_MALLOC(tmp_rhor,(nfftf)) 1807 tmp_rhor=rhor ! To avoid having to use intent(inout). 1808 call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,1,ngfftf,0) 1809 ABI_FREE(tmp_rhor) 1810 1811 rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf)) 1812 ABI_FREE(rhog_dp) 1813 ! 1814 ! Calculate the FFT index of each (G-Gp) vector and assign the value 1815 ! of the correspondent density simultaneously 1816 ngfft1=ngfftf(1) 1817 ngfft2=ngfftf(2) 1818 ngfft3=ngfftf(3) 1819 1820 ierr=0 1821 do ig=1,npwc 1822 do igp=1,npwc 1823 gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf) 1824 if (gmgp_idx/=0) then 1825 rhogg(ig,igp)=rhog(gmgp_idx) 1826 else 1827 ierr=ierr+1 1828 rhogg(ig,igp)=czero 1829 end if 1830 end do 1831 end do 1832 1833 if (ierr/=0) then 1834 write(msg,'(a,i0,3a)')& 1835 & 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,& 1836 & 'Enlarge the FFT mesh to get rid of this problem. ' 1837 ABI_WARNING(msg) 1838 end if 1839 1840 ABI_FREE(rhog) 1841 ! 1842 ! Now we have rhogg, calculate the M matrix (q+G1).(q+G2) n(G1-G2) 1843 ABI_MALLOC_OR_DIE(mm,(npwc,npwc), ierr) 1844 1845 do ig=1,npwc 1846 gpq(:)=gvec(:,ig)+qpt 1847 do igp=1,npwc 1848 gppq(:)=gvec(:,igp)+qpt 1849 qpg_dot_qpgp=zero 1850 do ii=1,3 1851 qpg_dot_qpgp=qpg_dot_qpgp+& 1852 & ( gpq(1)*b1(ii) +gpq(2)*b2(ii) +gpq(3)*b3(ii))*& 1853 & (gppq(1)*b1(ii)+gppq(2)*b2(ii)+gppq(3)*b3(ii)) 1854 end do 1855 mm(ig,igp)=rhogg(ig,igp)*qpg_dot_qpgp 1856 end do !igp 1857 end do !ig 1858 1859 !MG TODO too much memory in chi, we can do all this stuff inside a loop 1860 ABI_MALLOC_OR_DIE(chitmp,(npwc,npwc), ierr) 1861 1862 ABI_MALLOC_OR_DIE(chi,(npwc,npwc), ierr) 1863 1864 ABI_MALLOC(qplusg,(npwc)) 1865 ! 1866 ! Extract the full polarizability from \tilde \epsilon^{-1} 1867 ! \tilde\epsilon^{-1}_{G1 G2} = \delta_{G1 G2} + 4\pi \frac{\chi_{G1 G2}}{|q+G1| |q+G2|} 1868 chitmp(:,:)=epsm1(:,:) 1869 1870 qlist(:,1) = qpt 1871 call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q 1872 1873 do ig=1,npwc 1874 chitmp(ig,ig)=chitmp(ig,ig)-one 1875 end do 1876 do ig=1,npwc 1877 do igp=1,npwc 1878 chi(ig,igp)=chitmp(ig,igp)*qplusg(ig)*qplusg(igp)/four_pi 1879 end do 1880 end do 1881 ABI_FREE(chitmp) 1882 ! 1883 ! === Solve chi*X = Lambda M*X where Lambda=-1/em(q)**2 === 1884 ABI_MALLOC(eigval,(npwc)) 1885 1886 ABI_MALLOC_OR_DIE(eigvec,(npwc,npwc), ierr) 1887 ABI_MALLOC_OR_DIE(mtemp,(npwc,npwc), ierr) 1888 ABI_MALLOC_OR_DIE(chitmps,(npwc,npwc), ierr) 1889 ! 1890 ! Copy chi and mm into working arrays 1891 chitmps(:,:)=chi(:,:) 1892 mtemp(:,:)=mm(:,:) 1893 1894 call xhegv(1,"Vectors","Upper",npwc,chitmps,mtemp,eigval) 1895 ! 1896 ! Eigenvectors are normalized as : X_i^* M X_j = \delta_{ij} 1897 eigvec(:,:)=chitmps(:,:) 1898 ABI_FREE(mtemp) 1899 ABI_FREE(chitmps) 1900 ! 1901 ! === Calculate the plasmon pole parameters === 1902 ABI_MALLOC(tmp1,(npwc)) 1903 1904 ABI_MALLOC_OR_DIE(zz2,(npwc,npwc), ierr) 1905 ! 1906 ! good check: 1907 ! the lowest plasmon energy on gamma should be 1908 ! close to experimental plasma energy within an error of 10 % 1909 ! this error can be reduced further if one includes the non local 1910 ! commutators in the calculation of the polarizability at q==0 1911 zz2(:,:)=(0.0,0.0) 1912 1913 qlist(:,1) = qpt 1914 call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q 1915 1916 do ii=1,npwc 1917 ! 1918 ! keeping in mind that the above matrix is negative definite 1919 ! we might have a small problem with the eigval that correspond to large G vectors 1920 ! i.e. DM band index, where the eigevalues become very small with 1921 ! possibility of being small positive numbers (due to numerical problems) 1922 ! thus as a caution one can use the following condition 1923 ! this will not affect the result since such a huge plasmon energy give almost zero 1924 ! contribution to the self correlation energy 1925 1926 if (eigval(ii)>=zero) then 1927 eigval(ii) = -1.0d-4 1928 if (eigval(ii)>1.0d-3) then 1929 eigval(ii) = -1.0d-22 1930 write(msg,'(a,i6,a,es16.6)')& 1931 & ' Imaginary plasmon pole eigenenergy, eigenvector number ',ii,' with eigval',eigval(ii),ch10 1932 ABI_ERROR(msg) 1933 end if 1934 end if 1935 ! 1936 ! === Save plasmon energies === 1937 omegatw(ii)=SQRT(-1/eigval(ii)) 1938 ! 1939 ! Calculate and save scaled plasmon-pole eigenvectors 1940 ! defined as \sqrt{4\pi} \frac{Mx}{\sqrt{\tilde\omega} |q+G|} 1941 tmp1(:)=eigvec(:,ii) 1942 1943 do ig=1,npwc 1944 do igp=1,npwc 1945 zz2(ig,ii)=zz2(ig,ii)+mm(ig,igp)*tmp1(igp) ! z--->y 1946 end do 1947 bigomegatwsq(ig,ii)=SQRT(four_pi)*zz2(ig,ii)/SQRT(omegatw(ii)) 1948 bigomegatwsq(ig,ii)=bigomegatwsq(ig,ii)/qplusg(ig) 1949 end do 1950 1951 end do 1952 1953 ABI_FREE(tmp1) 1954 ABI_FREE(eigvec) 1955 ABI_FREE(eigval) 1956 ABI_FREE(zz2) 1957 1958 call destroy_mpi_enreg(MPI_enreg_seq) 1959 1960 ABI_FREE(qplusg) 1961 ABI_FREE(chi) 1962 ABI_FREE(rhogg) 1963 ABI_FREE(mm) 1964 1965 bar=REPEAT('-',80) 1966 write(msg,'(3a)')bar,ch10,' plasmon energies vs q vector shown for lowest 10 bands ' 1967 call wrtout(std_out,msg,'COLL') 1968 write(msg,'(2x,5x,10f7.3)')(REAL(omegatw(ig))*Ha_eV, ig=1,10) 1969 call wrtout(std_out,msg,'COLL') 1970 write(msg,'(a)')bar 1971 call wrtout(std_out,msg,'COLL') 1972 1973 write(msg,'(2a,f12.8,2a,3i5)')ch10,& 1974 & ' cppm4par : omega twiddle minval [eV] = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,& 1975 & ' omega twiddle min location = ',MINLOC(ABS(omegatw)) 1976 call wrtout(std_out,msg,'COLL') 1977 1978 DBG_EXIT("COLL") 1979 1980 end subroutine cppm4par
m_ppmodel/cqratio [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
cqratio
FUNCTION
Calculate qratio(G,Gp,q)= (q+G)\cdot(q+Gp) / |q+G|^2 needed for Hybertsen-Louie and Plasmonpole model
INPUTS
npwc=number of planewaves considered (used for the correlation part) gvec(3,npwc)=reduced coordinates of the plane waves q(3)=coordinates of q points gmet(3,3)=metric in reciprocal space gprimd(3,3)=reciprocal lattice vectors
OUTPUT
qratio(npwc,npwc)=(q+G).(q+Gp)
SOURCE
2004 subroutine cqratio(npwc,gvec,q,gmet,gprimd,qratio) 2005 2006 !Arguments ------------------------------------ 2007 !scalars 2008 integer,intent(in) :: npwc 2009 !arrays 2010 integer,intent(in) :: gvec(3,npwc) 2011 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),q(3) 2012 real(dp),intent(out) :: qratio(npwc,npwc) 2013 2014 !Local variables ------------------------------ 2015 !scalars 2016 integer :: ig,igp,ii 2017 real(dp) :: qpg_dot_qpgp 2018 !arrays 2019 real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),norm(npwc) 2020 2021 !************************************************************************ 2022 2023 b1=two_pi*gprimd(:,1); b2=two_pi*gprimd(:,2); b3=two_pi*gprimd(:,3) 2024 2025 norm(:)=zero; qratio=zero 2026 2027 !FIXME this loops have to be rewritten!!!! 2028 do ig=1,npwc 2029 gpq(:)=gvec(:,ig)+q 2030 norm(ig)=two_pi*SQRT(DOT_PRODUCT(gpq,MATMUL(gmet,gpq))) 2031 ! norm(ig)=normv(gpq,gmet,'g') 2032 end do 2033 do ig=1,npwc 2034 gpq(:)=gvec(:,ig)+q 2035 do igp=1,npwc 2036 gppq(:)=gvec(:,igp)+q 2037 qpg_dot_qpgp=zero 2038 ! qpg_dot_qpgp=vdotw(gpq,gppq,gmet,'g') 2039 do ii=1,3 2040 qpg_dot_qpgp=qpg_dot_qpgp+& 2041 & ( gpq(1)*b1(ii) + gpq(2)*b2(ii) + gpq(3)*b3(ii))*& 2042 & (gppq(1)*b1(ii) + gppq(2)*b2(ii) +gppq(3)*b3(ii)) 2043 end do 2044 2045 ! Now calculate qratio = (q+G).(q+Gp)/|q+G|^2 2046 ! when |q+G|^2 and (q+G).(q+Gp) are both zero set (q+G).(q+Gp)/|q+G|^2 = 1 2047 ! when |q+G|^2 is zero and |q+Gp| is not zero set (q+G).(q+Gp)/|q+G|^2 = 0 2048 if (norm(ig)<0.001) then 2049 if (norm(igp)<0.001) then ! Case q=0 and G=Gp=0 2050 qratio(ig,igp)=one 2051 else ! Case q=0 and G=0 and Gp !=0 2052 qratio(ig,igp)=zero 2053 end if 2054 else if (norm(igp)<0.001) then ! Case q=0 and G= !0 and Gp=0 2055 qratio(ig,igp)=zero 2056 else 2057 qratio(ig,igp)=qpg_dot_qpgp/norm(ig)**2 2058 end if 2059 2060 end do 2061 end do 2062 2063 end subroutine cqratio
m_ppmodel/get_ppm_eigenvalues [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
get_ppm_eigenvalues
FUNCTION
Constructs the inverse dielectri matrix starting from the plasmon-pole parameters and calculates the frequency-dependent eigenvalues for each of the nomega frequencies specifies in the array omega.
INPUTS
OUTPUT
SIDE EFFECTS
SOURCE
1034 subroutine get_ppm_eigenvalues(PPm,iqibz,zcut,nomega,omega,Vcp,eigenvalues) 1035 1036 !Arguments ------------------------------------ 1037 !scalars 1038 integer,intent(in) :: iqibz,nomega 1039 type(ppmodel_t),intent(in) :: PPm 1040 type(vcoul_t),intent(in) :: Vcp 1041 real(dp),intent(in) :: zcut 1042 !arrays 1043 complex(dpc),intent(in) :: omega(nomega) 1044 complex(dpc),intent(out) :: eigenvalues(PPm%npwc,nomega) 1045 1046 !Local variables------------------------------- 1047 !scalars 1048 integer :: info,lwork,negw,ig1,ig2,idx,sdim,iomega,ierr 1049 character(len=500) :: msg 1050 !arrays 1051 real(dp),allocatable :: ww(:),rwork(:) 1052 complex(dpc),allocatable :: work(:),Adpp(:),eigvec(:,:),wwc(:),vs(:,:),Afull(:,:) 1053 complex(dpc),allocatable :: em1q(:,:,:) 1054 logical,allocatable :: bwork(:) 1055 logical :: sortcplx !BUG in abilint 1056 1057 ! ************************************************************************* 1058 1059 ABI_CHECK(PPm%mqmem/=0,'mqmem==0 not implemented') 1060 1061 ABI_MALLOC(em1q,(PPm%npwc,PPm%npwc,nomega)) 1062 1063 call getem1_from_PPm(PPm,PPm%npwc,iqibz,zcut,nomega,omega,Vcp,em1q) 1064 1065 do iomega=1,nomega 1066 ! 1067 if (ABS(REAL(omega(iomega)))>0.00001) then 1068 !if (.TRUE.) then 1069 ! === Eigenvalues for a generic complex matrix === 1070 1071 lwork=4*2*PPm%npwc 1072 ABI_MALLOC(wwc,(PPm%npwc)) 1073 ABI_MALLOC(work,(lwork)) 1074 ABI_MALLOC(rwork,(PPm%npwc)) 1075 ABI_MALLOC(bwork,(PPm%npwc)) 1076 ABI_MALLOC(vs,(PPm%npwc,PPm%npwc)) 1077 ABI_MALLOC(Afull,(PPm%npwc,PPm%npwc)) 1078 Afull=em1q(:,:,iomega) 1079 1080 !for the moment no sort, maybe here I should sort using the real part? 1081 call ZGEES('V','N',sortcplx,PPm%npwc,Afull,PPm%npwc,sdim,wwc,vs,PPm%npwc,work,lwork,rwork,bwork,info) 1082 if (info/=0) then 1083 write(msg,'(2a,i10)')' get_ppm_eigenvalues : Error in ZGEES, diagonalizing complex matrix, info = ',info 1084 call wrtout(std_out,msg,'COLL') 1085 end if 1086 1087 eigenvalues(:,iomega)=wwc(:) 1088 1089 ABI_FREE(wwc) 1090 ABI_FREE(work) 1091 ABI_FREE(rwork) 1092 ABI_FREE(bwork) 1093 ABI_FREE(vs) 1094 ABI_FREE(Afull) 1095 1096 else 1097 ! === Hermitian Case === 1098 lwork=2*PPm%npwc-1 1099 ABI_MALLOC(ww,(PPm%npwc)) 1100 ABI_MALLOC(work,(lwork)) 1101 ABI_MALLOC(rwork,(3*PPm%npwc-2)) 1102 ABI_MALLOC(eigvec,(PPm%npwc,PPm%npwc)) 1103 1104 ABI_MALLOC_OR_DIE(Adpp,(PPm%npwc*(PPm%npwc+1)/2), ierr) 1105 write(std_out,*) 'in hermitian' 1106 1107 idx=0 1108 do ig2=1,PPm%npwc 1109 do ig1=1,ig2 1110 idx=idx+1 1111 Adpp(idx)=em1q(ig1,ig2,iomega) 1112 end do 1113 end do 1114 1115 ! For the moment we require also the eigenvectors. 1116 call ZHPEV('V','U',PPm%npwc,Adpp,ww,eigvec,PPm%npwc,work,rwork,info) 1117 1118 if (info/=0) then 1119 write(msg,'(2a,i10)')' get_ppm_eigenvalues : Error diagonalizing matrix, info = ',info 1120 call wrtout(std_out,msg,'COLL') 1121 end if 1122 negw = (COUNT((REAL(ww)<tol6))) 1123 if (negw/=0) then 1124 write(msg,'(a,i0,a,i0,a,f8.4)')& 1125 & 'Found negative eigenvalues. No. ',negw,' at iqibz= ',iqibz,' minval= ',MINVAL(REAL(ww)) 1126 ABI_WARNING(msg) 1127 end if 1128 1129 eigenvalues(:,iomega)=ww(:) 1130 1131 ABI_FREE(ww) 1132 ABI_FREE(work) 1133 ABI_FREE(rwork) 1134 ABI_FREE(eigvec) 1135 ABI_FREE(Adpp) 1136 end if 1137 ! 1138 end do !iomega 1139 1140 ABI_FREE(em1q) 1141 1142 end subroutine get_ppm_eigenvalues
m_ppmodel/getem1_from_ppm [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
getem1_from_ppm
FUNCTION
Calculate the symmetrized inverse dielectric matrix from the parameters of the plasmon-pole model.
INPUTS
OUTPUT
SIDE EFFECTS
SOURCE
786 subroutine getem1_from_ppm(PPm,mpwc,iqibz,zcut,nomega,omega,Vcp,em1q,only_ig1,only_ig2) 787 788 !Arguments ------------------------------------ 789 !scalars 790 integer,intent(in) :: mpwc,iqibz,nomega 791 type(ppmodel_t),intent(in) :: PPm 792 type(vcoul_t),intent(in) :: Vcp 793 real(dp),intent(in) :: zcut 794 integer,optional,intent(in) :: only_ig1,only_ig2 795 !arrays 796 complex(dpc),intent(in) :: omega(nomega) 797 complex(dpc),intent(out) :: em1q(mpwc,mpwc,nomega) 798 799 !Local variables------------------------------- 800 !scalars 801 integer :: ig1,ig2,io,idm,ig1_min,ig2_min,ig1_max,ig2_max 802 real(dp) :: den 803 complex(dpc) :: qpg1,qpg2,ug1,ug2 804 complex(dpc) :: delta,em1ggp,otw,zzpq,yg1,yg2,bot1,bot2,chig1g2 805 !character(len=500) :: msg 806 807 ! ************************************************************************* 808 809 ABI_CHECK(PPm%mqmem/=0,'mqmem==0 not implemented') 810 811 !TODO zcut should be an entry in PPm 812 delta=CMPLX(zero,zcut) 813 814 ! To save memory, a particular combination of 815 ! ig1 and ig2 can be selected 816 ig1_min = 1 817 ig2_min = 1 818 ig1_max = PPm%npwc 819 ig2_max = PPm%npwc 820 if (present(only_ig1)) then 821 ig1_min = only_ig1 822 ig1_max = only_ig1 823 end if 824 if (present(only_ig2)) then 825 ig2_min = only_ig2 826 ig2_max = only_ig2 827 end if 828 829 select case (PPm%model) 830 case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE) 831 do io=1,nomega 832 ! 833 do ig2=ig2_min,ig2_max 834 do ig1=ig1_min,ig1_max 835 !den = omega(io)**2-REAL(PPm%omegatw(iqibz)%vals(ig1,ig2)**2) 836 !if (den**2<zcut**2) den = omega(io)**2-REAL( (PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 ) 837 den = omega(io)**2 - REAL( (PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 ) 838 em1ggp = PPm%bigomegatwsq(iqibz)%vals(ig1,ig2)/den 839 if (ig1==ig2) em1ggp=em1ggp+one 840 em1q(ig1,ig2,io)=em1ggp 841 !em1q(ig1,ig2,io)=em1ggp*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz) 842 end do 843 end do 844 ! 845 end do !io 846 847 case (PPM_LINDEN_HORSH) 848 !TODO Check coefficients 849 do io=1,nomega 850 ! 851 do ig2=ig2_min,ig2_max 852 do ig1=ig1_min,ig1_max 853 ! 854 em1ggp=czero 855 do idm=1,PPm%npwc 856 !den=omega(io)**2-(PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 857 !em1w(io)=em1w(io)+eigvec(ig1,idm,iqibz)*conjg(eigvec(ig2,idm,iqibz))*bigomegatwsq(ig1,ig2,iqibz)/den 858 ug1 = PPm%eigpot(iqibz)%vals(ig1,idm) 859 ug2 = PPm%eigpot(iqibz)%vals(ig2,idm) 860 otw = PPm%bigomegatwsq(iqibz)%vals(idm,1)*PPm%omegatw(iqibz)%vals(idm,1) 861 zzpq=PPm%bigomegatwsq(iqibz)%vals(idm,1) 862 den=half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) )) 863 em1ggp=em1ggp+ug1*CONJG(ug2)*den 864 !eigenvalues(idm,io)=one + half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) )) 865 end do 866 if (ig2==ig1) em1ggp=em1ggp+one 867 em1q(ig1,ig2,io)=em1ggp 868 end do !ig1 869 end do !ig2 870 ! 871 end do !iomega 872 873 case (PPM_ENGEL_FARID) 874 ! Make e^-1 875 do io=1,nomega 876 ! 877 do ig2=ig2_min,ig2_max 878 qpg2=one/Vcp%vc_sqrt(ig2,iqibz) 879 do ig1=ig1_min,ig1_max 880 qpg1=one/Vcp%vc_sqrt(ig1,iqibz) 881 882 chig1g2=czero 883 do idm=1,PPm%npwc 884 otw =PPm%omegatw(iqibz)%vals(idm,1) 885 bot1=PPm%bigomegatwsq(iqibz)%vals(ig1,idm) 886 bot2=PPm%bigomegatwsq(iqibz)%vals(ig2,idm) 887 yg1=SQRT(otw/four_pi)*qpg1*bot1 888 yg2=SQRT(otw/four_pi)*qpg2*bot2 889 chig1g2=chig1g2 + yg1*CONJG(yg2)/(omega(io)**2-(otw-delta)**2) 890 end do 891 892 em1ggp=four_pi*chig1g2/(qpg1*qpg2) 893 if (ig1==ig2) em1ggp=em1ggp+one 894 em1q(ig1,ig2,io)=em1ggp !*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz) 895 end do !ig1 896 end do !ig2 897 ! 898 end do !iomega 899 900 case default 901 ABI_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model))) 902 end select 903 904 end subroutine getem1_from_ppm
m_ppmodel/getem1_from_ppm_one_ggp [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
getem1_from_ppm
FUNCTION
Same as getem1_from_ppm, but does it for a single set of G,G' vectors
INPUTS
OUTPUT
SIDE EFFECTS
SOURCE
924 subroutine getem1_from_ppm_one_ggp(PPm,iqibz,zcut,nomega,omega,Vcp,em1q,ig1,ig2) 925 926 !Arguments ------------------------------------ 927 !scalars 928 integer,intent(in) :: iqibz,nomega 929 type(ppmodel_t),intent(in) :: PPm 930 type(vcoul_t),intent(in) :: Vcp 931 real(dp),intent(in) :: zcut 932 integer, intent(in) :: ig1,ig2 933 !arrays 934 complex(dpc),intent(in) :: omega(nomega) 935 complex(dpc),intent(out) :: em1q(nomega) 936 937 !Local variables------------------------------- 938 !scalars 939 integer :: io,idm !,ig1_min,ig2_min,ig2_max 940 real(dp) :: den 941 complex(dpc) :: qpg1,qpg2,ug1,ug2 942 complex(dpc) :: delta,em1ggp,otw,zzpq,yg1,yg2,bot1,bot2,chig1g2 943 !character(len=500) :: msg 944 945 ! ************************************************************************* 946 947 ABI_CHECK(PPm%mqmem/=0,'mqmem==0 not implemented') 948 949 !TODO zcut should be an entry in PPm 950 delta=CMPLX(zero,zcut) 951 952 select case (PPm%model) 953 954 case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE) 955 do io=1,nomega 956 !den = omega(io)**2-REAL(PPm%omegatw(iqibz)%vals(ig1,ig2)**2) 957 !if (den**2<zcut**2) den = omega(io)**2-REAL( (PPm%omegatw(iqibz)%value(ig1,ig2)-delta)**2 ) 958 den = omega(io)**2-REAL( (PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 ) 959 em1ggp = PPm%bigomegatwsq(iqibz)%vals(ig1,ig2)/den 960 if (ig1==ig2) em1ggp=em1ggp+one 961 em1q(io)=em1ggp 962 !em1q(io)=em1ggp*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz) 963 end do !io 964 965 case (PPM_LINDEN_HORSH) 966 !TODO Check coefficients 967 do io=1,nomega 968 em1ggp=czero 969 do idm=1,PPm%npwc 970 !den=omega(io)**2-(PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 971 !em1w(io)=em1w(io)+eigvec(ig1,idm,iqibz)*conjg(eigvec(ig2,idm,iqibz))*bigomegatwsq(ig1,ig2,iqibz)/den 972 ug1 =PPm%eigpot(iqibz)%vals(ig1,idm) 973 ug2 =PPm%eigpot(iqibz)%vals(ig2,idm) 974 otw =PPm%bigomegatwsq(iqibz)%vals(idm,1)*PPm%omegatw(iqibz)%vals(idm,1) 975 zzpq=PPm%bigomegatwsq(iqibz)%vals(idm,1) 976 den=half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) )) 977 em1ggp=em1ggp+ug1*CONJG(ug2)*den 978 !eigenvalues(idm,io)=one + half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) )) 979 end do 980 981 if (ig2==ig1) em1ggp=em1ggp+one 982 em1q(io)=em1ggp 983 end do !iomega 984 985 case (PPM_ENGEL_FARID) 986 ! Make e^-1 987 do io=1,nomega 988 ! 989 qpg2=one/Vcp%vc_sqrt(ig2,iqibz) 990 qpg1=one/Vcp%vc_sqrt(ig1,iqibz) 991 992 chig1g2=czero 993 do idm=1,PPm%npwc 994 otw =PPm%omegatw(iqibz)%vals(idm,1) 995 bot1=PPm%bigomegatwsq(iqibz)%vals(ig1,idm) 996 bot2=PPm%bigomegatwsq(iqibz)%vals(ig2,idm) 997 yg1=SQRT(otw/four_pi)*qpg1*bot1 998 yg2=SQRT(otw/four_pi)*qpg2*bot2 999 chig1g2=chig1g2 + yg1*CONJG(yg2)/(omega(io)**2-(otw-delta)**2) 1000 end do 1001 1002 em1ggp=four_pi*chig1g2/(qpg1*qpg2) 1003 if (ig1==ig2) em1ggp=em1ggp+one 1004 em1q(io)=em1ggp !*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz) 1005 1006 end do !iomega 1007 1008 case default 1009 ABI_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model))) 1010 end select 1011 1012 end subroutine getem1_from_ppm_one_ggp
m_ppmodel/new_setup_ppmodel [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
new_setup_ppmodel
FUNCTION
Initialize some values of several arrays of the PPm datastructure that are used in case of plasmonpole calculations Just a wrapper around different plasmonpole routines.
INPUTS
Cryst<crystal_t>=Info on the unit cell and crystal symmetries. Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix iq_ibz=Index of the q-point in the BZ. npwe=number of G vectors for the correlation part nomega=number of frequencies in $\epsilon^{-1}$ omega=frequencies in epsm1_ggw epsm1_ggw(npwe,npwe,nomega)=the inverse dielctric matrix ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) nfftf=the number of points in the FFT mesh (for this processor) rhor_tot(nfftf)=the total charge in real space PPm<ppmodel_t>:
SIDE EFFECTS
PPm<ppmodel_t>: == if ppmodel 1 or 2 == %omegatw and %bigomegatwsq == if ppmodel 3 == %omegatw, %bigomegatwsq and %eigpot == if ppmodel 4 == %omegatw and %bigomegatwsq
NOTES
* FFT parallelism not implemented. * TODO: rhor_tot should be replaced by rhog_tot
SOURCE
2433 subroutine new_setup_ppmodel(PPm,iq_ibz,Cryst,Qmesh,npwe,nomega,omega,epsm1_ggw,nfftf,gvec,ngfftf,rhor_tot) 2434 2435 !Arguments ------------------------------------ 2436 !scalars 2437 integer,intent(in) :: nfftf,npwe,nomega,iq_ibz 2438 type(kmesh_t),intent(in) :: Qmesh 2439 type(crystal_t),intent(in) :: Cryst 2440 type(ppmodel_t),intent(inout) :: PPm 2441 !arrays 2442 integer,intent(in) :: gvec(3,npwe),ngfftf(18) 2443 real(dp),intent(in) :: rhor_tot(nfftf) 2444 complex(dpc),intent(in) :: omega(nomega) 2445 complex(gwpc),intent(in) :: epsm1_ggw(npwe,npwe,nomega) 2446 2447 !Local variables------------------------------- 2448 !scalars 2449 real(dp) :: n_at_G_zero 2450 character(len=500) :: msg 2451 !scalars 2452 real(dp) :: qpt(3) 2453 ! ************************************************************************* 2454 2455 DBG_ENTER("COLL") 2456 2457 !@ppmodel_t 2458 ! 2459 if (PPm%has_q(iq_ibz) /= PPM_TAB_ALLOCATED) then 2460 ABI_ERROR("ppmodel tables are not allocated") 2461 end if 2462 2463 qpt = Qmesh%ibz(:,iq_ibz) 2464 PPm%has_q(iq_ibz) = PPM_TAB_STORED 2465 ! 2466 ! Calculate plasmonpole parameters 2467 ! TODO ppmodel==1 by default, should be set to 0 if AC and CD 2468 SELECT CASE (PPm%model) 2469 2470 CASE (PPM_NONE) 2471 ABI_COMMENT('Skipping Plasmonpole model calculation') 2472 2473 CASE (PPM_GODBY_NEEDS) 2474 ! Note: the q-dependency enters only through epsilon^-1. 2475 call cppm1par(npwe,nomega,omega,PPm%drude_plsmf,epsm1_ggw,& 2476 & PPm%omegatw(iq_ibz)%vals,PPm%bigomegatwsq(iq_ibz)%vals) 2477 2478 CASE (PPM_HYBERTSEN_LOUIE) 2479 call cppm2par(qpt,npwe,epsm1_ggw(:,:,1),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,Cryst%gmet,& 2480 & PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals,PPm%invalid_freq) 2481 ! 2482 ! Quick-and-dirty change of the plasma frequency. Never executed in standard runs. 2483 if (PPm%force_plsmf>tol6) then ! Integrate the real-space density 2484 n_at_G_zero = SUM(rhor_tot(:))/nfftf 2485 ! Change the prefactor 2486 write(msg,'(2(a,es16.8))') 'Forced ppmfreq:',PPm%force_plsmf*Ha_eV,' nelec/ucvol:',n_at_G_zero 2487 ABI_WARNING(msg) 2488 PPm%force_plsmf = (PPm%force_plsmf**2)/(four_pi*n_at_G_zero) 2489 PPm%bigomegatwsq(iq_ibz)%vals = PPm%force_plsmf * PPm%bigomegatwsq(iq_ibz)%vals 2490 PPm%omegatw(iq_ibz)%vals = PPm%force_plsmf * PPm%omegatw(iq_ibz)%vals 2491 write(msg,'(a,es16.8)') 'Plasma frequency forced in HL ppmodel, new prefactor is:',PPm%force_plsmf 2492 ABI_WARNING(msg) 2493 end if 2494 2495 CASE (PPM_LINDEN_HORSH) 2496 ! TODO Check better double precision, this routine is in a messy state 2497 call cppm3par(qpt,npwe,epsm1_ggw(:,:,1),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,& 2498 & PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals(:,1),PPm%eigpot(iq_ibz)%vals) 2499 2500 CASE (PPM_ENGEL_FARID) ! TODO Check better double precision, this routine is in a messy state 2501 if ((ALL(ABS(qpt)<1.0e-3))) qpt = GW_Q0_DEFAULT ! FIXME 2502 2503 call cppm4par(qpt,npwe,epsm1_ggw(:,:,1),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,& 2504 & PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals(:,1)) 2505 2506 CASE DEFAULT 2507 ABI_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model))) 2508 END SELECT 2509 2510 DBG_EXIT("COLL") 2511 2512 end subroutine new_setup_ppmodel
m_ppmodel/ppm_free [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
ppm_free
FUNCTION
Deallocate all associated pointers defined in a variable of type ppmodel_t.
SIDE EFFECTS
PPm<ppmodel_t>=All dynamic memory is released.
SOURCE
360 subroutine ppm_free(PPm) 361 362 !Arguments ------------------------------------ 363 type(ppmodel_t),intent(inout) :: PPm 364 365 !Local variables------------------------------- 366 !scalars 367 integer :: dim_q,iq_ibz 368 369 ! ********************************************************************* 370 371 !@ppmodel_t 372 373 ! Be careful here to avoid dangling pointers. 374 if (associated(PPm%bigomegatwsq_qbz)) then 375 select case (PPm%bigomegatwsq_qbz_stat) 376 case (PPM_ISALLOCATED) 377 call array_free(PPm%bigomegatwsq_qbz) 378 case (PPM_ISPOINTER) 379 nullify(PPm%bigomegatwsq_qbz) 380 end select 381 end if 382 383 if (associated(PPm%omegatw_qbz)) then 384 select case (PPm%omegatw_qbz_stat) 385 case (PPM_ISALLOCATED) 386 call array_free(PPm%omegatw_qbz) 387 case (PPM_ISPOINTER) 388 nullify(PPm%omegatw_qbz) 389 end select 390 end if 391 392 if (associated(PPm%eigpot_qbz)) then 393 select case (PPm%eigpot_qbz_stat) 394 case (PPM_ISALLOCATED) 395 call array_free(PPm%eigpot_qbz) 396 case (PPM_ISPOINTER) 397 nullify(PPm%eigpot_qbz) 398 end select 399 end if 400 401 dim_q=PPm%nqibz; if (PPm%mqmem==0) dim_q=1 402 403 #if 0 404 do iq_ibz=1,dim_q 405 call ppm_table_free(PPm,iq_ibz) 406 end do 407 ABI_FREE(PPm%bigomegatwsq) 408 ABI_FREE(PPm%omegatw) 409 ABI_FREE(PPm%eigpot) 410 411 #else 412 if (allocated(PPm%bigomegatwsq)) then 413 do iq_ibz=1,dim_q 414 call array_free(PPm%bigomegatwsq(iq_ibz)) 415 end do 416 ABI_FREE(PPm%bigomegatwsq) 417 end if 418 ! 419 if (allocated(PPm%omegatw)) then 420 do iq_ibz=1,dim_q 421 call array_free(PPm%omegatw(iq_ibz)) 422 end do 423 ABI_FREE(PPm%omegatw) 424 end if 425 ! 426 if (allocated(PPm%eigpot)) then 427 do iq_ibz=1,dim_q 428 call array_free(PPm%eigpot(iq_ibz)) 429 end do 430 ABI_FREE(PPm%eigpot) 431 end if 432 #endif 433 434 ! logical flags must be deallocated here. 435 ABI_SFREE(PPm%keep_q) 436 ABI_SFREE(PPm%has_q) 437 438 end subroutine ppm_free
m_ppmodel/ppm_get_qbz [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
ppm_get_qbz
FUNCTION
Calculates the plasmonpole matrix elements in the full BZ zone.
INPUTS
PPm<ppmodel_t>=data type containing information on the plasmonpole technique Gsph<gsphere_t>=data related to the G-sphere %grottb %phmSGt Qmesh<kmesh_t>=Info on the q-mesh %nbz=number if q-points in the BZ %tab(nbz)=index of the symmeric q-point in the IBZ, for each point in the BZ %tabo(nbz)=the operation that rotates q_ibz onto \pm q_bz (depending on tabi) %tabi(nbz)=-1 if time-reversal has to be considered, 1 otherwise iq_bz=Index of the q-point in the BZ where PPmodel parameters have to be symmetrized
OUTPUT
botsq otq eig (only if PPm%ppmodel==3)
NOTES
In the present implementation we are not considering a possible umklapp vector G0. In this case,indeed, the equation is different since we have to consider G-G0. There is however a check in sigma * Remember the symmetry properties of \tilde\espilon^{-1} If q_bz=Sq_ibz+G0: $\epsilon^{-1}_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau}\epsilon^{-1}_{G1,G2)}(q) If time-reversal symmetry can be used then : $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau}\epsilon^{-1}_{-S^{-1}(G1+Go),-S^{-1}(G2+G0)}^*(q) * Notice that eig is used only if PPm%model==3
SOURCE
219 subroutine ppm_get_qbz(PPm,Gsph,Qmesh,iq_bz,botsq,otq,eig) 220 221 !Arguments ------------------------------------ 222 !scalars 223 integer,intent(in) :: iq_bz 224 type(ppmodel_t),target,intent(inout) :: PPm 225 type(gsphere_t),target,intent(in) :: Gsph 226 type(kmesh_t),intent(in) :: Qmesh 227 complex(gwpc),intent(out) :: botsq(:,:),otq(:,:),eig(:,:) 228 229 !Local variables------------------------------- 230 !scalars 231 integer :: ii,jj,iq_ibz,itim_q,isym_q,iq_curr,isg1,isg2 232 !character(len=500) :: msg 233 !arrays 234 integer, ABI_CONTIGUOUS pointer :: grottb(:) 235 complex(gwpc),ABI_CONTIGUOUS pointer :: phsgt(:),bigomegatwsq(:,:),omegatw(:,:) 236 ! ********************************************************************* 237 238 !@ppmodel_t 239 ! Save the index of the q-point for checking purpose. 240 PPm%iq_bz = iq_bz 241 242 ! Here there is a problem with the small q, still cannot use BZ methods 243 iq_ibz=Qmesh%tab(iq_bz) 244 isym_q=Qmesh%tabo(iq_bz) 245 itim_q=(3-Qmesh%tabi(iq_bz))/2 246 247 !call Qmesh%get_bz_item(iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred) 248 249 iq_curr=iq_ibz; if (PPm%mqmem==0) iq_curr=1 250 251 grottb => Gsph%rottb (1:PPm%npwc,itim_q,isym_q) 252 phsgt => Gsph%phmSGt(1:PPm%npwc,isym_q) 253 254 bigomegatwsq => PPm%bigomegatwsq(iq_curr)%vals 255 omegatw => PPm%omegatw(iq_curr)%vals 256 257 ! Symmetrize the PPM parameters 258 SELECT CASE (PPm%model) 259 CASE (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE) 260 ! Plasmon pole frequencies otq are invariant under symmetry 261 !$omp parallel do private(isg1, isg2) 262 do jj=1,PPm%npwc 263 isg2 = grottb(jj) 264 do ii=1,PPm%npwc 265 isg1 = grottb(ii) 266 botsq(isg1,isg2) = bigomegatwsq(ii,jj)*phsgt(ii)*CONJG(phsgt(jj)) 267 otq (isg1,isg2) = omegatw(ii,jj) 268 end do 269 end do 270 271 CASE (PPM_LINDEN_HORSH) 272 ! * For notations see pag 22 of Quasiparticle Calculations in solid (Aulbur et al) 273 ! If q_bz=Sq_ibz+G0 then: 274 ! 275 ! $\omega^2_{ii}(q_bz) = \omega^2_{ii}(q)$ (otq array) 276 ! $\alpha_{ii}(q_bz) = \alpha_{ii}(q)$ (botq array 277 ! $\Phi_{SG-G0}(q_bz) = \Phi_{G}(q) e^{-iSG.t}$ (eigenvectors of e^{-1}, eig array) 278 ! 279 do ii=1,PPm%npwc ! DM bands index 280 botsq(ii,1) = bigomegatwsq(ii,1) 281 otq (ii,1) = omegatw (ii,1) 282 do jj=1,PPm%npwc 283 eig(grottb(jj),ii) = PPm%eigpot(iq_curr)%vals(jj,ii) * phsgt(jj) 284 end do 285 end do 286 if (itim_q==2) eig=CONJG(eig) ! Time-reversal 287 288 CASE (PPM_ENGEL_FARID) 289 ! * For notations see pag 23 of Quasiparticle Calculations in solid (Aulbur et al) 290 ! If q_bz=Sq_ibz+G0 then: 291 ! 292 ! $\omega^2_{ii}(q_bz) = \omega^2_{ii}(q)$ (otq array) 293 ! $y_{SG-G0}(q_bz) = y_{G}(q) e^{-iSG.t}$ (y=Lx) 294 ! 295 do ii=1,PPm%npwc ! DM bands index 296 otq(ii,1) = omegatw(ii,1) 297 do jj=1,PPm%npwc 298 botsq(grottb(jj),ii) = bigomegatwsq(jj,ii)*phsgt(jj) 299 end do 300 end do 301 302 CASE DEFAULT 303 ABI_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model))) 304 END SELECT 305 ! 306 ! * Take into account time-reversal symmetry. 307 if (itim_q==2) then 308 !$omp parallel workshare 309 botsq=CONJG(botsq) 310 !$omp end parallel workshare 311 end if 312 313 end subroutine ppm_get_qbz
m_ppmodel/ppm_init [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
ppm_init
FUNCTION
Initialize dimensions and other useful variables related to the PPmodel
INPUTS
ppmodel= drude_plsmf=
SIDE EFFECTS
PPm<ppmodel_t>=Arrays needed to store the ppmodel parameters are allocated with proper dimensions according to the plasmon-pole model.
SOURCE
532 subroutine ppm_init(PPm,mqmem,nqibz,npwe,ppmodel,drude_plsmf,invalid_freq) 533 534 !Arguments ------------------------------------ 535 integer,intent(in) :: mqmem,nqibz,npwe,ppmodel,invalid_freq 536 real(dp),intent(in) :: drude_plsmf 537 type(ppmodel_t),intent(out) :: PPm 538 539 !Local variables------------------------------- 540 !scalars 541 integer :: dim_q,iq_ibz,ierr 542 logical :: ltest 543 !character(len=500) :: msg 544 ! ********************************************************************* 545 546 DBG_ENTER("COLL") 547 548 !@ppmodel_t 549 call ppm_nullify(PPm) 550 551 PPm%invalid_freq=invalid_freq 552 PPm%nqibz = nqibz 553 PPm%mqmem = mqmem 554 ltest = (PPm%mqmem==0.or.PPm%mqmem==PPm%nqibz) 555 ! write(std_out,'(a,I0)') ' PPm%mqmem = ',PPm%mqmem 556 ! write(std_out,'(a,I0)') ' PPm%nqibz = ',PPm%nqibz 557 ! ABI_CHECK(ltest,'Wrong value for mqmem') 558 559 PPm%npwc = npwe 560 PPm%model = ppmodel 561 PPm%drude_plsmf = drude_plsmf 562 PPm%userho = 0 563 if (ANY(ppmodel==(/PPM_HYBERTSEN_LOUIE, PPM_LINDEN_HORSH, PPM_ENGEL_FARID/))) PPm%userho=1 564 565 ABI_MALLOC(PPm%keep_q,(nqibz)) 566 PPm%keep_q=.FALSE.; if (PPm%mqmem>0) PPm%keep_q=.TRUE. 567 568 ABI_MALLOC(PPm%has_q,(nqibz)) 569 PPm%has_q=PPM_NOTAB 570 ! 571 ! Full q-mesh is stored or out-of-memory solution. 572 dim_q=PPm%nqibz; if (PPm%mqmem==0) dim_q=1 573 574 ABI_MALLOC(PPm%bigomegatwsq, (dim_q)) 575 ABI_MALLOC(PPm%omegatw,(dim_q)) 576 ABI_MALLOC(PPm%eigpot,(dim_q)) 577 578 SELECT CASE (PPm%model) 579 580 CASE (PPM_NONE) 581 ABI_WARNING("Called with ppmodel==0") 582 PPm%dm2_botsq = 0 583 PPm%dm2_otq = 0 584 PPm%dm_eig = 0 585 RETURN 586 587 CASE (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE) 588 PPm%dm2_botsq = PPm%npwc 589 PPm%dm2_otq = PPm%npwc 590 PPm%dm_eig = 1 ! Should be set to 0, but g95 doesnt like zero-sized arrays 591 592 CASE (PPM_LINDEN_HORSH) 593 PPm%dm2_botsq = 1 594 PPm%dm2_otq = 1 595 PPm%dm_eig = PPm%npwc 596 597 CASE (PPM_ENGEL_FARID) 598 PPm%dm2_botsq = PPm%npwc 599 PPm%dm2_otq = 1 600 PPm%dm_eig = 1 ! Should be set to 0, but g95 doesnt like zero-sized arrays 601 602 CASE DEFAULT 603 ABI_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model))) 604 END SELECT 605 ! 606 ! Allocate tables depending on the value of keep_q. 607 do iq_ibz=1,dim_q 608 !%if (keep_q(iq_ibz)) then 609 #if 1 610 ABI_MALLOC_OR_DIE(PPm%bigomegatwsq(iq_ibz)%vals, (PPm%npwc,PPm%dm2_botsq), ierr) 611 ABI_MALLOC_OR_DIE(PPm%omegatw(iq_ibz)%vals, (PPm%npwc,PPm%dm2_otq), ierr) 612 ABI_MALLOC_OR_DIE(PPm%eigpot(iq_ibz)%vals, (PPm%dm_eig,PPm%dm_eig), ierr) 613 !%endif 614 #else 615 call ppm_mallocq(PPm,iq_ibz) 616 #endif 617 !%end if 618 end do 619 620 DBG_EXIT("COLL") 621 622 end subroutine ppm_init
m_ppmodel/ppm_mallocq [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
ppm_mallocq
FUNCTION
Allocate the ppmodel tables for the selected q-point in the IBZ. INPUT iq_ibz=Index of the q-point in the IBZ.
SIDE EFFECTS
PPm<ppmodel_t>=PPm tables are allocated.
SOURCE
458 subroutine ppm_mallocq(PPm,iq_ibz) 459 460 !Arguments ------------------------------------ 461 integer,intent(in) :: iq_ibz 462 type(ppmodel_t),intent(inout) :: PPm 463 464 ! ********************************************************************* 465 466 !@ppmodel_t 467 ABI_MALLOC(PPm%bigomegatwsq(iq_ibz)%vals, (PPm%npwc,PPm%dm2_botsq)) 468 469 ABI_MALLOC(PPm%omegatw(iq_ibz)%vals, (PPm%npwc,PPm%dm2_otq)) 470 471 ABI_MALLOC(PPm%eigpot(iq_ibz)%vals, (PPm%dm_eig,PPm%dm_eig)) 472 473 PPm%has_q(iq_ibz) = PPM_TAB_ALLOCATED 474 475 end subroutine ppm_mallocq
m_ppmodel/ppm_nullify [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
ppm_nullify
FUNCTION
Nullify dynamic entities in a ppmodel_t object.
INPUTS
OUTPUT
SOURCE
331 subroutine ppm_nullify(PPm) 332 333 !Arguments ------------------------------------ 334 type(ppmodel_t),intent(inout) :: PPm 335 ! ********************************************************************* 336 337 !@ppmodel_t 338 ! types 339 nullify(Ppm%bigomegatwsq_qbz) 340 nullify(Ppm%omegatw_qbz) 341 nullify(Ppm%eigpot_qbz) 342 343 end subroutine ppm_nullify
m_ppmodel/ppm_symmetrizer [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
ppm_symmetrizer
FUNCTION
Symmetrize the plasmonpole matrix elements in the full BZ zone.
INPUTS
iq_bz=Index of the q-point in the BZ where the PPmodel parameters are wanted. Gsph<gsphere_t>=data related to the G-sphere. Cryst<crystal_t>=Info on the unit cell and crystal symmetries. Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix iq_ibz=Index of the q-point in the BZ. npwe=number of G vectors for the correlation part nomega=number of frequencies in $\epsilon^{-1}$ omega=frequencies in epsm1_ggw epsm1_ggw(npwe,npwe,nomega)=the inverse dielctric matrix ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) nfftf=the number of points in the FFT mesh (for this processor) rhor_tot(nfftf)=the total charge in real space
SIDE EFFECTS
PPm<ppmodel_t>=data type containing information on the plasmonpole technique. Internal tables are modified so that they (point|store) the plasmon-pole parameters for the specified q-point in the BZ.
SOURCE
2285 subroutine ppm_symmetrizer(PPm,iq_bz,Cryst,Qmesh,Gsph,npwe,nomega,omega,epsm1_ggw,nfftf,ngfftf,rhor_tot) 2286 2287 !Arguments ------------------------------------ 2288 !scalars 2289 integer,intent(in) :: nfftf,npwe,nomega,iq_bz 2290 type(crystal_t),intent(in) :: Cryst 2291 type(gsphere_t),intent(in) :: Gsph 2292 type(kmesh_t),intent(in) :: Qmesh 2293 type(ppmodel_t),target,intent(inout) :: PPm 2294 !arrays 2295 integer,intent(in) :: ngfftf(18) 2296 real(dp),intent(in) :: rhor_tot(nfftf) 2297 complex(dpc),intent(in) :: omega(nomega) 2298 complex(gwpc),intent(in) :: epsm1_ggw(npwe,npwe,nomega) 2299 2300 !Local variables------------------------------- 2301 !scalars 2302 integer :: iq_ibz,itim_q,isym_q,iq_curr,ierr 2303 logical :: q_isirred 2304 !character(len=500) :: msg 2305 !arrays 2306 real(dp) :: qbz(3) 2307 ! ********************************************************************* 2308 2309 !@ppmodel_t 2310 ! 2311 ! Save the index of the q-point in the BZ for checking purpose. 2312 PPm%iq_bz = iq_bz 2313 ! 2314 ! Here there is a problem with the small q, still cannot use BZ methods 2315 !iq_ibz=Qmesh%tab(iq_bz) 2316 !isym_q=Qmesh%tabo(iq_bz) 2317 !itim_q=(3-Qmesh%tabi(iq_bz))/2 2318 2319 call qmesh%get_bz_item(iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred) 2320 iq_curr=iq_ibz; if (PPm%mqmem==0) iq_curr=1 2321 ! 2322 ! ======================================================= 2323 ! ==== Branching for in-core or out-of-core solution ==== 2324 ! ======================================================= 2325 ! 2326 if (PPm%has_q(iq_ibz)==PPM_NOTAB) then 2327 ! Allocate the tables for this q_ibz 2328 call ppm_mallocq(PPm,iq_ibz) 2329 end if 2330 2331 if (PPm%has_q(iq_ibz)==PPM_TAB_ALLOCATED) then 2332 ! Calculate the ppmodel tables for this q_ibz 2333 call new_setup_ppmodel(PPm,iq_ibz,Cryst,Qmesh,npwe,nomega,omega,epsm1_ggw,nfftf,Gsph%gvec,ngfftf,rhor_tot) !Optional 2334 end if 2335 2336 if (q_isirred) then 2337 ! Symmetrization is not needed. 2338 if (PPm%bigomegatwsq_qbz_stat==PPM_ISALLOCATED) then 2339 ABI_FREE(PPm%bigomegatwsq_qbz%vals) 2340 end if 2341 if (PPm%omegatw_qbz_stat==PPM_ISALLOCATED) then 2342 ABI_FREE(PPm%omegatw_qbz%vals) 2343 end if 2344 if (PPm%eigpot_qbz_stat==PPM_ISALLOCATED) then 2345 ABI_FREE(PPm%eigpot_qbz%vals) 2346 end if 2347 2348 ! Point the data in memory and change the status. 2349 PPm%bigomegatwsq_qbz => PPm%bigomegatwsq(iq_ibz) 2350 PPm%bigomegatwsq_qbz_stat = PPM_ISPOINTER 2351 2352 PPm%omegatw_qbz => PPm%omegatw(iq_ibz) 2353 PPm%omegatw_qbz_stat = PPM_ISPOINTER 2354 2355 PPm%eigpot_qbz => PPm%eigpot(iq_ibz) 2356 PPm%eigpot_qbz_stat = PPM_ISPOINTER 2357 else 2358 ! Allocate memory if not done yet. 2359 if (PPm%bigomegatwsq_qbz_stat==PPM_ISPOINTER) then 2360 nullify(PPm%bigomegatwsq_qbz) 2361 ABI_MALLOC_OR_DIE(PPm%bigomegatwsq_qbz%vals, (PPm%npwc,PPm%dm2_botsq), ierr) 2362 PPm%bigomegatwsq_qbz_stat=PPM_ISALLOCATED 2363 end if 2364 2365 if (PPm%omegatw_qbz_stat==PPM_ISPOINTER) then 2366 nullify(PPm%omegatw_qbz) 2367 ABI_MALLOC_OR_DIE(PPm%omegatw_qbz%vals,(PPm%npwc,PPm%dm2_otq), ierr) 2368 PPm%omegatw_qbz_stat=PPM_ISALLOCATED 2369 end if 2370 2371 if (PPm%eigpot_qbz_stat==PPM_ISPOINTER) then 2372 nullify(PPm%eigpot_qbz) 2373 ABI_MALLOC_OR_DIE(PPm%eigpot_qbz%vals,(PPm%dm_eig,PPm%dm_eig), ierr) 2374 PPm%eigpot_qbz_stat=PPM_ISALLOCATED 2375 end if 2376 ! 2377 ! Calculate new table for this q-point in the BZ. 2378 ! Beware: Dimensions should not change. 2379 !botsq => PPm%bigomegatwsq_qbz%vals 2380 !otq => PPm%omegatw_qbz%vals 2381 !eig => PPm%eigpot_qbz%vals 2382 2383 call ppm_get_qbz(PPm,Gsph,Qmesh,iq_bz,PPm%bigomegatwsq_qbz%vals,& 2384 & PPm%omegatw_qbz%vals,PPm%eigpot_qbz%vals) 2385 2386 ! Release the table in the IBZ if required. 2387 if (.not.PPm%keep_q(iq_ibz)) call ppm_table_free(PPm,iq_ibz) 2388 end if 2389 2390 end subroutine ppm_symmetrizer
m_ppmodel/ppm_table_free [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
ppm_table_free
FUNCTION
Free the ppmodel tables for the selected q-point in the IBZ. INPUT iq_ibz=Index of the q-point in the IBZ.
SIDE EFFECTS
PPm<ppmodel_t>=PPm tables are deallocated.
SOURCE
495 subroutine ppm_table_free(PPm,iq_ibz) 496 497 !Arguments ------------------------------------ 498 integer,intent(in) :: iq_ibz 499 type(ppmodel_t),intent(inout) :: PPm 500 501 ! ********************************************************************* 502 503 !@ppmodel_t 504 if (allocated(PPm%bigomegatwsq)) call array_free(PPm%bigomegatwsq(iq_ibz)) 505 if (allocated(PPm%omegatw)) call array_free(PPm%omegatw(iq_ibz)) 506 if (allocated(PPm%eigpot)) call array_free(PPm%eigpot(iq_ibz)) 507 508 PPm%has_q(iq_ibz) = PPM_NOTAB 509 510 end subroutine ppm_table_free
m_ppmodel/ppm_times_ket [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
ppm_times_ket
FUNCTION
Calculate the contribution to self-energy operator using a plasmon-pole model.
INPUTS
nomega=Number of frequencies. nspinor=Number of spinorial components. npwc=Number of G vectors in the plasmon pole. npwx=number of G vectors in rhotwgp, i.e. no. of G-vectors for the exchange part. theta_mu_minus_e0i= $\theta(\mu-\epsilon_{k-q,b1,s}), defines if the state is occupied or not zcut=Small imaginary part to avoid the divergence. (see related input variable) omegame0i(nomega)=Frequencies used to evaluate \Sigma_c ($\omega$ - $\epsilon_i)$ PPm<ppmodel_t>=structure gathering info on the Plasmon-pole technique. %model=type plasmon pole model %dm2_botsq= 1 if model==3, =npwc if model== 4, 1 for all the other cases %dm2_otq= 1 if model==3, =1 if model== 4, 1 for all the other cases %dm_eig=npwc if model=3, 0 otherwise %botsq(npwc,dm2_botsq)=Plasmon pole parameters for this q-point. %eig(dm_eig,dm_eig)=The eigvectors of the symmetrized inverse dielectric matrix for this q point (first index for G, second index for bands) %otq(npwc,dm2_otq)=Plasmon pole parameters for this q-point. rhotwgp(npwx)=oscillator matrix elements divided by |q+G| i.e. $\frac{\langle b1 k-q s | e^{-i(q+G)r | b2 k s \rangle}{|q+G|}$
OUTPUT
sigcme(nomega) (to be described), only relevant if ppm3 or ppm4
ket(npwc,nomega):
=== model==1,2 ====
ket(G,omega) = Sum_G2 conjg(rhotw(G)) * Omega(G,G2) * rhotw(G2)
---------------------------------------------------
2 omegatw(G,G2) (omega-E_i + omegatw(G,G2)(2f-1))
SOURCE
2556 subroutine ppm_times_ket(PPm,nspinor,npwc,nomega,rhotwgp,omegame0i,zcut,theta_mu_minus_e0i,npwx,ket,sigcme) 2557 2558 !Arguments ------------------------------------ 2559 !scalars 2560 integer,intent(in) :: nomega,npwc,npwx,nspinor 2561 real(dp),intent(in) :: theta_mu_minus_e0i,zcut 2562 type(ppmodel_t),intent(in) :: PPm 2563 !arrays 2564 real(dp),intent(in) :: omegame0i(nomega) 2565 complex(gwpc),intent(in) :: rhotwgp(npwx*nspinor) 2566 complex(gwpc),intent(inout) :: ket(npwc*nspinor,nomega) 2567 complex(gwpc),intent(out) :: sigcme(nomega) 2568 2569 !Local variables------------------------------- 2570 !scalars 2571 integer :: ig,igp,ii,ios,ispinor,spadc,spadx 2572 real(dp) :: den,ff,inv_den,omegame0i_io,otw,twofm1,twofm1_zcut 2573 complex(gwpc) :: ct,num,numf,rhotwgdp_igp 2574 logical :: fully_occupied,totally_empty 2575 !arrays 2576 complex(gwpc),allocatable :: rhotwgdpcc(:) 2577 complex(gwpc), ABI_CONTIGUOUS pointer :: botsq(:,:),eig(:,:),otq(:,:) 2578 2579 !************************************************************************* 2580 2581 SELECT CASE (PPm%model) 2582 CASE (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE) 2583 botsq => PPm%bigomegatwsq_qbz%vals !(1:npwc,1:PPm%dm2_botsq) 2584 otq => PPm%omegatw_qbz%vals !(1:npwc,1:PPm%dm2_otq) 2585 fully_occupied =(ABS(theta_mu_minus_e0i-one)<0.001) 2586 totally_empty =(ABS(theta_mu_minus_e0i )<0.001) 2587 2588 do ispinor=1,nspinor 2589 spadx=(ispinor-1)*npwx; spadc=(ispinor-1)*npwc 2590 2591 if (.not.totally_empty) then 2592 ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(\mu-e_s) / (\omega+\omegat_{G1G2}-e_s-i\delta) 2593 twofm1_zcut=zcut 2594 !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den) 2595 do ios=1,nomega 2596 omegame0i_io=omegame0i(ios) 2597 do igp=1,npwc 2598 rhotwgdp_igp=rhotwgp(spadx+igp) 2599 do ig=1,npwc 2600 otw=DBLE(otq(ig,igp)) !in principle otw -> otw - ieta 2601 num = botsq(ig,igp)*rhotwgdp_igp 2602 den = omegame0i_io+otw 2603 if (den**2>zcut**2) then 2604 ket(spadc+ig,ios)=ket(spadc+ig,ios) + num/(den*otw)*theta_mu_minus_e0i 2605 else 2606 ket(spadc+ig,ios)=ket(spadc+ig,ios) + & 2607 & num*CMPLX(den,twofm1_zcut)/((den**2+twofm1_zcut**2)*otw)*theta_mu_minus_e0i 2608 end if 2609 end do !ig 2610 end do !igp 2611 end do !ios 2612 end if !not totally empty 2613 2614 if (.not.(fully_occupied)) then 2615 ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(e_s-\mu) / (\omega-\omegat_{G1G2}-e_s+i\delta) 2616 twofm1_zcut=-zcut 2617 !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den) 2618 do ios=1,nomega 2619 omegame0i_io=omegame0i(ios) 2620 do igp=1,npwc 2621 rhotwgdp_igp=rhotwgp(spadx+igp) 2622 do ig=1,npwc 2623 otw=DBLE(otq(ig,igp)) !in principle otw -> otw - ieta 2624 num = botsq(ig,igp)*rhotwgdp_igp 2625 den=omegame0i_io-otw 2626 if (den**2>zcut**2) then 2627 ket(spadc+ig,ios)=ket(spadc+ig,ios) + num/(den*otw)*(one-theta_mu_minus_e0i) 2628 else 2629 ket(spadc+ig,ios)=ket(spadc+ig,ios) + & 2630 & num*CMPLX(den,twofm1_zcut)/((den**2+twofm1_zcut**2)*otw)*(one-theta_mu_minus_e0i) 2631 end if 2632 end do !ig 2633 end do !igp 2634 end do !ios 2635 end if !not fully occupied 2636 2637 end do !ispinor 2638 2639 !$omp parallel workshare 2640 ket=ket*half 2641 !$omp end parallel workshare 2642 2643 CASE (PPM_LINDEN_HORSH,PPM_ENGEL_FARID) 2644 ABI_CHECK(nspinor == 1, "nspinor/=1 not allowed") 2645 2646 botsq => PPm%bigomegatwsq_qbz%vals !(1:npwc,1:PPm%dm2_botsq) 2647 otq => PPm%omegatw_qbz%vals !(1:npwc,1:PPm%dm2_otq) 2648 eig => PPm%eigpot_qbz%vals !(1:PPm%dm_eig,1:PPm%dm_eig) 2649 2650 ! rho-twiddle(G) is formed, introduce rhotwgdpcc, for speed reason 2651 ABI_MALLOC(rhotwgdpcc,(npwx)) 2652 2653 ff=theta_mu_minus_e0i ! occupation number f (include poles if ...) 2654 twofm1=two*ff-one ! 2f-1 2655 twofm1_zcut=twofm1*zcut 2656 rhotwgdpcc(:)=CONJG(rhotwgp(:)) 2657 2658 do ios=1,nomega 2659 omegame0i_io=omegame0i(ios) 2660 ct=czero_gw 2661 do ii=1,npwc ! Loop over the DM bands 2662 num=czero_gw 2663 2664 SELECT CASE (PPm%model) 2665 CASE (PPM_LINDEN_HORSH) 2666 ! Calculate \beta (eq. 106 pag 47) 2667 do ig=1,npwc 2668 num=num+rhotwgdpcc(ig)*eig(ig,ii) 2669 end do 2670 numf=num*CONJG(num) !MG this means that we cannot do SCGW 2671 numf=numf*botsq(ii,1) 2672 2673 CASE (PPM_ENGEL_FARID) 2674 do ig=1,npwc 2675 num=num+rhotwgdpcc(ig)*botsq(ig,ii) 2676 end do 2677 numf=num*CONJG(num) !MG this means that we cannot do SCGW 2678 2679 CASE DEFAULT 2680 ABI_ERROR("Wrong PPm%model") 2681 END SELECT 2682 2683 !numf=num*CONJG(num) !MG this means that we cannot do SCGW 2684 !if (PPm%model==3) numf=numf*botsq(ii,1) 2685 2686 otw=DBLE(otq(ii,1)) ! in principle otw -> otw - ieta 2687 den=omegame0i_io+otw*twofm1 2688 2689 if (den**2>zcut**2) then 2690 inv_den=one/den 2691 ct=ct+numf*inv_den 2692 else 2693 inv_den=one/((den**2+twofm1_zcut**2)) 2694 ct=ct+numf*CMPLX(den,twofm1_zcut)*inv_den 2695 end if 2696 2697 end do !ii DM bands 2698 sigcme(ios)=ct*half 2699 end do !ios 2700 ABI_FREE(rhotwgdpcc) 2701 2702 CASE DEFAULT 2703 ABI_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model))) 2704 END SELECT 2705 2706 end subroutine ppm_times_ket
m_ppmodel/ppmodel_t [ Types ]
[ Top ] [ m_ppmodel ] [ Types ]
NAME
ppmodel_t
FUNCTION
For the GW part of ABINIT, this structure datatype gathers all the information on the Plasmonpole technique used in the calculations
NOTES
If you modify this datatype, please check there there is no creation/destruction/copy routine, declared in another part of ABINIT, that might need to take into account your modification. Procedures that should be modified to reflect the changes done in the datatype have the tag @ppmodel_t.
SOURCE
82 type,public :: ppmodel_t 83 84 !integers 85 integer :: dm2_botsq 86 ! =npwc if ppmodel=1,2; =1 if ppmodel=3,4 87 88 integer :: dm_eig 89 ! =npwc if ppmodel=3; =0 if ppmodel=1,2,4 90 91 integer :: dm2_otq 92 ! =npwc if ppmodel=1,2; =1 if ppmodel=3,4 93 94 integer :: invalid_freq 95 ! what to do when PPM frequencies are invalid 96 97 integer :: model 98 ! The type of Plasmonpole model 99 100 integer :: mqmem 101 ! =nqibz if in-core solutions, =0 for out-of-core for which the last dimension in the PPm arrays has size 1. 102 103 integer :: nqibz 104 ! Number of q-points in the IBZ 105 106 integer :: npwc 107 ! Number of G vectors in $\tilde \epsilon $ 108 109 integer :: userho 110 ! 1 if the ppmodel requires rho(G). 111 112 integer :: iq_bz=0 113 ! The index of the q-point in the BZ that is referenced by the internal pointer. 114 115 integer :: bigomegatwsq_qbz_stat = PPM_ISPOINTER 116 ! Status of bigomegatwsq. 117 118 integer :: omegatw_qbz_stat = PPM_ISPOINTER 119 ! Status of omegatw_qbz. 120 121 integer :: eigpot_qbz_stat = PPM_ISPOINTER 122 ! Status of new_eigpot_qbz_stat. 123 124 real(dp) :: drude_plsmf 125 ! Drude plasma frequency 126 127 real(dp) :: force_plsmf=zero 128 ! Force plasma frequency to that set by ppmfreq (not used if set zero). 129 130 ! arrays 131 logical,allocatable :: keep_q(:) 132 ! .TRUE. if the ppmodel tables for this q in the IBZ are kept in memory. 133 134 integer,allocatable :: has_q(:) 135 ! Flag defining the status of the tables for the different q. See the PPM_TAB flags. 136 137 type(array2_gwpc_t),pointer :: bigomegatwsq_qbz => null() 138 ! (Points|Stores) the symmetrized plasmon pole parameters $\tilde\Omega^2_{G Gp}(q_bz)$. 139 140 type(array2_gwpc_t),pointer :: omegatw_qbz => null() 141 ! (Points|Stores) the symmetrized plasmon pole parameters $\tilde\omega_{G Gp}(q_bz)$. 142 143 type(array2_gwpc_t),pointer :: eigpot_qbz => null() 144 ! (Points|Stores) the eigvectors of the symmetrized inverse dielectric matrix 145 146 type(array2_gwpc_t),allocatable :: bigomegatwsq(:) 147 ! bigomegatwsq(nqibz)%value(npwc,dm2_botsq) 148 ! Plasmon pole parameters $\tilde\Omega^2_{G Gp}(q)$. 149 150 type(array2_gwpc_t),allocatable :: omegatw(:) 151 ! omegatw(nqibz)%value(npwc,dm2_otq) 152 ! Plasmon pole parameters $\tilde\omega_{G Gp}(q)$. 153 154 type(array2_gwpc_t),allocatable :: eigpot(:) 155 ! eigpot(nqibz)%value(dm_eig,dm_eig) 156 ! Eigvectors of the symmetrized inverse dielectric matrix 157 158 end type ppmodel_t 159 160 public :: ppm_get_qbz ! Symmetrize the PPm parameters in the BZ. 161 public :: ppm_nullify ! Nullify all pointers 162 public :: ppm_init ! Initialize dimensions and pointers 163 public :: ppm_free ! Destruction method. 164 public :: setup_ppmodel ! Main Driver 165 public :: new_setup_ppmodel ! Main Driver 166 public :: getem1_from_PPm ! Reconstruct e^{-1}(w) from PPm. 167 public :: getem1_from_PPm_one_ggp ! Reconstruct e^{-1}(w) from PPm for one G,G' pair 168 public :: get_ppm_eigenvalues 169 public :: calc_sig_ppm ! Matrix elements of the self-energy with ppmodel. 170 public :: ppm_times_ket ! Matrix elements of the self-energy with ppmodel. 171 public :: ppm_symmetrizer 172 public :: cqratio
m_ppmodel/setup_ppmodel [ Functions ]
[ Top ] [ m_ppmodel ] [ Functions ]
NAME
setup_ppmodel
FUNCTION
Initialize some values of several arrays of the PPm datastructure that are used in case of plasmonpole calculations This is a wrapper around different plasmonpole routines.
INPUTS
Cryst<crystal_t>=Info on the unit cell and crystal symmetries. Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix %nibz=number of irreducible q-points %ibz(3,%nibz)=the irred q-point npwe=number of G vectors for the correlation part nomega=number of frequencies in $\epsilon^{-1}$ omega=frequencies in epsm1 epsm1=the inverse dielctric matrix ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) nfftf=the number of points in the FFT mesh (for this processor) rhor_tot(nfftf)=the total charge in real space PPm<ppmodel_t>: %ppmodel=the type of plasmonpole model
OUTPUT
SIDE EFFECTS
PPm<ppmodel_t>: == if ppmodel 1 or 2 == %omegatw and %bigomegatwsq == if ppmodel 3 == %omegatw, %bigomegatwsq and %eigpot == if ppmodel 4 == %omegatw and %bigomegatwsq
NOTES
* FFT parallelism not implemented. * TODO: rhor_tot should be replaced by rhog_tot
SOURCE
669 subroutine setup_ppmodel(PPm,Cryst,Qmesh,npwe,nomega,omega,epsm1,nfftf,gvec,ngfftf,rhor_tot,& 670 & iqiA) !Optional 671 672 !Arguments ------------------------------------ 673 !scalars 674 integer,intent(in) :: nfftf,npwe,nomega 675 integer,intent(in),optional :: iqiA 676 type(kmesh_t),intent(in) :: Qmesh 677 type(crystal_t),intent(in) :: Cryst 678 type(ppmodel_t),intent(inout) :: PPm 679 !arrays 680 integer,intent(in) :: gvec(3,npwe),ngfftf(18) 681 real(dp),intent(in) :: rhor_tot(nfftf) 682 complex(dpc),intent(in) :: omega(nomega) 683 complex(gwpc),intent(in) :: epsm1(:,:,:,:) 684 685 !Local variables------------------------------- 686 !scalars 687 integer :: nqiA,iq_ibz 688 real(dp) :: n_at_G_zero 689 logical :: single_q 690 character(len=500) :: msg 691 !scalars 692 real(dp) :: qpt(3) 693 ! ************************************************************************* 694 695 DBG_ENTER("COLL") 696 697 !@ppmodel_t 698 ! 699 ! === if iqiA is present, then consider only one qpoint to save memory === 700 ! * This means the object has been already initialized 701 nqiA=Qmesh%nibz; single_q=.FALSE. 702 if (PRESENT(iqiA)) then 703 nqiA=1; single_q=.TRUE. 704 end if 705 ! 706 ! Allocate plasmonpole parameters 707 ! TODO ppmodel==1 by default, should be set to 0 if AC and CD 708 SELECT CASE (PPm%model) 709 710 CASE (PPM_NONE) 711 ABI_COMMENT(' Skipping Plasmompole model calculation') 712 713 CASE (PPM_GODBY_NEEDS) 714 ! * Note: the q-dependency enters only through epsilon^-1. 715 do iq_ibz=1,nqiA 716 call cppm1par(npwe,nomega,omega,PPm%drude_plsmf,& 717 & epsm1(:,:,:,iq_ibz),PPm%omegatw(iq_ibz)%vals,PPm%bigomegatwsq(iq_ibz)%vals) 718 end do 719 720 CASE (PPM_HYBERTSEN_LOUIE) 721 do iq_ibz=1,nqiA 722 qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA) 723 724 call cppm2par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,Cryst%gmet,& 725 & PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals,PPm%invalid_freq) 726 end do 727 ! 728 ! Quick-and-dirty change of the plasma frequency. Never executed in standard runs. 729 if (PPm%force_plsmf>tol6) then ! Integrate the real-space density 730 n_at_G_zero = SUM(rhor_tot(:))/nfftf 731 ! Change the prefactor 732 write(msg,'(2(a,es16.8))') 'Forced ppmfreq:',PPm%force_plsmf*Ha_eV,' nelec/ucvol:',n_at_G_zero 733 ABI_WARNING(msg) 734 PPm%force_plsmf = (PPm%force_plsmf**2)/(four_pi*n_at_G_zero) 735 do iq_ibz=1,PPm%nqibz 736 PPm%bigomegatwsq(iq_ibz)%vals = PPm%force_plsmf * PPm%bigomegatwsq(iq_ibz)%vals 737 PPm%omegatw(iq_ibz)%vals = PPm%force_plsmf * PPm%omegatw(iq_ibz)%vals 738 end do 739 write(msg,'(a,es16.8)') 'Plasma frequency forced in HL ppmodel, new prefactor is:',PPm%force_plsmf 740 ABI_WARNING(msg) 741 end if 742 743 CASE (PPM_LINDEN_HORSH) ! TODO Check better double precision, this routine is in a messy state 744 do iq_ibz=1,nqiA 745 qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA) 746 747 call cppm3par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,& 748 & PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals(:,1),PPm%eigpot(iq_ibz)%vals) 749 end do 750 751 CASE (PPM_ENGEL_FARID) ! TODO Check better double precision, this routine is in a messy state 752 do iq_ibz=1,nqiA 753 qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA) 754 if ((ALL(ABS(qpt)<1.0e-3))) qpt = GW_Q0_DEFAULT ! FIXME 755 756 call cppm4par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,& 757 & PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals(:,1)) 758 end do 759 760 CASE DEFAULT 761 ABI_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model))) 762 END SELECT 763 764 DBG_EXIT("COLL") 765 766 end subroutine setup_ppmodel