TABLE OF CONTENTS
ABINIT/dfptnl_exc3 [ Functions ]
NAME
dfptnl_exc3
FUNCTION
Compute the third-order xc energy. Take into account the contribution of the term \frac{d}{d \lambda} \frac{\delta^2 E_{Hxc}}{\delta n(r) \delta n(r\prim)} (seventh term of Eq. (110) of X. Gonze, PRA 52, 1096 (1995) [[cite:Gonze1995]]). The following is essentially the 4th and the 3rd terms of PRB 71,125107 [[cite:Veithen2005]]. However, we treat the nonlinear xc core correction in a slightly different way, in order to keep the symmetry between pert1,pert2 and pert3. It helps for debugging. Namely, here we consider Exc as a functional depending on the TOTAL density (core+valence), so it does not depend explicitely on the perturbation, and the term given above is always zero. The "lost" terms are recovered adding the derivative of the core densities for EVERY perturbations, and not only to 'pert2', as in the first case.
COPYRIGHT
Copyright (C) 2020-2024 ABINIT group (LB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cplex= if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX k3xc(nfftf,nk3xc)=third-order exchange-correlation kernel mpi_enreg=MPI-parallelisation information nk3xc=second dimension of the array k3xc nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90) nfftotf=total number of real space fine grid points nspden = number of spin-density components rho1r1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i1pert) rho1r2(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i2pert) rho1r3(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i3pert) ucvol=volume of the unit cell xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc (i1pert) xccc3d2(cplex*n3xccc)=3D change in core charge density, see n3xccc (i2pert) xccc3d3(cplex*n3xccc)=3D change in core charge density, see n3xccc (i3pert)
OUTPUT
exc3(2)=real and imaginray part of the exchange correlation energy term
SIDE EFFECTS
SOURCE
1177 subroutine dfptnl_exc3(cplex,exc3,k3xc,mpi_enreg,nk3xc,nfftf,nfftotf,nspden,rho1r1,rho2r1,rho3r1,ucvol,xccc3d1,xccc3d2,xccc3d3) 1178 1179 !Arguments ------------------------------------ 1180 !scalars 1181 integer,intent(in) :: cplex,nk3xc,nfftf,nfftotf,nspden 1182 real(dp),intent(in) :: ucvol 1183 type(MPI_type),intent(inout) :: mpi_enreg 1184 1185 !arrays 1186 real(dp),intent(in) :: k3xc(nfftf,nk3xc) 1187 real(dp),intent(in) :: rho1r1(cplex*nfftf,nspden),rho2r1(cplex*nfftf,nspden),rho3r1(cplex*nfftf,nspden) 1188 real(dp),intent(in) :: xccc3d1(cplex*nfftf),xccc3d2(cplex*nfftf),xccc3d3(cplex*nfftf) 1189 real(dp),intent(out) :: exc3(2) 1190 1191 !Local variables------------------------------- 1192 !scalars 1193 integer :: ifft,ifft_im,ifft_re 1194 real(dp) :: rho2ur,rho2ui,rho2dr,rho2di,rho3ur,rho3ui,rho3dr,rho3di 1195 ! character(len=1000) :: msg 1196 !arrays 1197 real(dp),allocatable :: rho1r1_tot(:,:),xc_tmp(:,:) 1198 1199 !*********************************************************************** 1200 1201 DBG_ENTER("COLL") 1202 1203 ABI_MALLOC(xc_tmp,(cplex*nfftf,nspden)) 1204 ABI_MALLOC(rho1r1_tot,(cplex*nfftf,nspden)) 1205 if (nspden==1)then 1206 1207 if (cplex==1) then 1208 do ifft=1,nfftf 1209 rho1r1_tot(ifft,1) = rho1r1(ifft,1) + xccc3d1(ifft) 1210 rho2ur = rho2r1(ifft,1)+xccc3d2(ifft) 1211 rho3ur = rho3r1(ifft,1)+xccc3d3(ifft) 1212 xc_tmp(ifft,1)= k3xc(ifft,1)*rho2ur*rho3ur 1213 end do 1214 else 1215 do ifft=1,nfftf 1216 ifft_re = 2*ifft-1 ! Real part 1217 ifft_im = 2*ifft ! Imaginary part 1218 rho1r1_tot(ifft_re,1) = rho1r1(ifft_re,1) + xccc3d1(ifft_re) 1219 rho1r1_tot(ifft_im,1) = rho1r1(ifft_im,1) + xccc3d1(ifft_im) 1220 rho2ur = rho2r1(ifft_re,1)+xccc3d2(ifft_re) 1221 rho2ui = rho2r1(ifft_im,1)+xccc3d2(ifft_im) 1222 rho3ur = rho3r1(ifft_re,1)+xccc3d3(ifft_re) 1223 rho3ui = rho3r1(ifft_im,1)+xccc3d3(ifft_im) 1224 xc_tmp(ifft_re,1)= k3xc(ifft,1)*(rho2ur*rho3ur-rho2ui*rho3ui) 1225 xc_tmp(ifft_im,1)= k3xc(ifft,1)*(rho2ur*rho3ui+rho2ui*rho3ur) 1226 end do 1227 end if 1228 1229 else if (nspden==2) then 1230 1231 ! Remember : rhor(...,1) = total density ( n_up(r) + n_down(r) ) 1232 ! rhor(...,2) = up density ( n_up(r) ) 1233 ! But : pot(...,1) = up potential ( v_up(r) ) 1234 ! pot(...,2) = down potential ( v_down(r) ) 1235 1236 if (cplex==1) then 1237 do ifft=1,nfftf 1238 rho1r1_tot(ifft,1) = rho1r1(ifft,1) + xccc3d1(ifft) ! 1 tot 1239 rho1r1_tot(ifft,2) = rho1r1(ifft,2) + xccc3d1(ifft)*half ! 1 up 1240 rho2ur = rho2r1(ifft,2) + xccc3d2(ifft)*half ! 2 up 1241 rho2dr = rho2r1(ifft,1)-rho2r1(ifft,2) + xccc3d2(ifft)*half ! 2 down 1242 rho3ur = rho3r1(ifft,2) + xccc3d3(ifft)*half ! 3 up 1243 rho3dr = rho3r1(ifft,1)-rho3r1(ifft,2) + xccc3d3(ifft)*half ! 3 down 1244 ! uuu uud 1245 xc_tmp(ifft,1)= k3xc(ifft,1)*rho2ur*rho3ur + k3xc(ifft,2)*rho2ur*rho3dr + & 1246 ! udu udd 1247 & k3xc(ifft,2)*rho2dr*rho3ur + k3xc(ifft,3)*rho2dr*rho3dr 1248 ! duu dud 1249 xc_tmp(ifft,2)= k3xc(ifft,2)*rho2ur*rho3ur + k3xc(ifft,3)*rho2ur*rho3dr + & 1250 ! ddu ddd 1251 & k3xc(ifft,3)*rho2dr*rho3ur + k3xc(ifft,4)*rho2dr*rho3dr 1252 end do 1253 1254 else ! cplex = 2 1255 1256 do ifft=1,nfftf 1257 ifft_re = 2*ifft-1 ! Real part 1258 ifft_im = 2*ifft ! Imaginary part 1259 rho1r1_tot(ifft_re,1) = rho1r1(ifft_re,1) + xccc3d1(ifft_re) ! 1 tot re 1260 rho1r1_tot(ifft_im,1) = rho1r1(ifft_im,1) + xccc3d1(ifft_im) ! 1 tot im 1261 rho1r1_tot(ifft_re,2) = rho1r1(ifft_re,2) + xccc3d1(ifft_re)*half ! 1 up re 1262 rho1r1_tot(ifft_im,2) = rho1r1(ifft_im,2) + xccc3d1(ifft_im)*half ! 1 up im 1263 rho2ur = rho2r1(ifft_re,2) + xccc3d2(ifft_re)*half ! 2 up re 1264 rho2ur = rho2r1(ifft_im,2) + xccc3d2(ifft_im)*half ! 2 up im 1265 rho2dr = rho2r1(ifft_re,1)-rho2r1(ifft_re,2) + xccc3d2(ifft_re)*half ! 2 down re 1266 rho2dr = rho2r1(ifft_im,1)-rho2r1(ifft_im,2) + xccc3d2(ifft_im)*half ! 2 down im 1267 rho3ur = rho3r1(ifft_re,2) + xccc3d3(ifft_re)*half ! 3 up re 1268 rho3ur = rho3r1(ifft_im,2) + xccc3d3(ifft_im)*half ! 3 up im 1269 rho3dr = rho3r1(ifft_re,1)-rho3r1(ifft_re,2) + xccc3d3(ifft_re)*half ! 3 down re 1270 rho3dr = rho3r1(ifft_im,1)-rho3r1(ifft_im,2) + xccc3d3(ifft_im)*half ! 3 down im 1271 ! Real part: 1272 ! uuu uud 1273 xc_tmp(ifft_re,1)= k3xc(ifft,1)*(rho2ur*rho3ur-rho2ui*rho3ui) + k3xc(ifft,2)*(rho2ur*rho3dr-rho2ui*rho3di) + & 1274 ! udu udd 1275 & k3xc(ifft,2)*(rho2dr*rho3ur-rho2di*rho3ui) + k3xc(ifft,3)*(rho2dr*rho3dr-rho2di*rho3di) 1276 ! duu dud 1277 xc_tmp(ifft_re,2)= k3xc(ifft,2)*(rho2ur*rho3ur-rho2ui*rho3ui) + k3xc(ifft,3)*(rho2ur*rho3dr-rho2ui*rho3di) + & 1278 ! ddu ddd 1279 & k3xc(ifft,3)*(rho2dr*rho3ur-rho2di*rho3ui) + k3xc(ifft,4)*(rho2dr*rho3dr-rho2di*rho3di) 1280 ! Imaginary part: 1281 ! uuu uud 1282 xc_tmp(ifft_im,1)= k3xc(ifft,1)*(rho2ur*rho3ui+rho2ui*rho3ur) + k3xc(ifft,2)*(rho2ur*rho3di+rho2ui*rho3dr) + & 1283 ! udu udd 1284 & k3xc(ifft,2)*(rho2dr*rho3ui+rho2di*rho3ur) + k3xc(ifft,3)*(rho2dr*rho3di+rho2di*rho3dr) 1285 ! duu dud 1286 xc_tmp(ifft_im,2)= k3xc(ifft,2)*(rho2ur*rho3ui+rho2ui*rho3ur) + k3xc(ifft,3)*(rho2ur*rho3di+rho2ui*rho3dr) + & 1287 ! ddu ddd 1288 & k3xc(ifft,3)*(rho2dr*rho3ui+rho2di*rho3ur) + k3xc(ifft,4)*(rho2dr*rho3di+rho2di*rho3dr) 1289 end do 1290 1291 end if 1292 1293 else 1294 ABI_BUG('DFPTNL_PERT is implemented only for nspden=1 or 2') 1295 end if 1296 1297 call dotprod_vn(cplex,rho1r1_tot,exc3(1),exc3(2),nfftf,nfftotf,nspden,2,xc_tmp,ucvol,mpi_comm_sphgrid=mpi_enreg%comm_fft) 1298 ABI_FREE(xc_tmp) 1299 ABI_FREE(rho1r1_tot) 1300 1301 DBG_EXIT("COLL") 1302 1303 end subroutine dfptnl_exc3
ABINIT/dfptnl_pert [ Functions ]
NAME
dfptnl_pert
FUNCTION
Compute the linear response part to the 3dte. The main inputs are : - GS WFs and Hamiltonian (cg,gs_hamkq) - 1st-order WFs for three perturbations i1pert/i1dir,i2pert/i2dir,i3pert/i3dir (cg1,cg2,cg3) - 1st-order potentials for i2pert (vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert) - 1st-order WFs DDK,DDE and 2nd-order WF DKDE (ddk_f)
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (LB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = array for planewave coefficients of wavefunctions cg1 = first derivative of cg with respect the perturbation i1pert cg2 = first derivative of cg with respect the perturbation i2pert cg3 = first derivative of cg with respect the perturbation i3pert cplex= if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q k3xc(nfftf,nk3xc)=third-order exchange-correlation kernel indsy1(4,nsym1,natom)=indirect indexing array for atom labels i1dir,i2dir,i3dir=directions of the corresponding perturbations i1pert,i2pert,i3pert = type of perturbation that has to be computed kg(3,mpw*mkmem_rbz)=reduced planewave coordinates mband = maximum number of bands mgfft=maximum size of 1D FFTs mkmem_rbz = maximum number of k points which can fit in core memory mk1mem = maximum number of k points for first-order WF which can fit in core memory mpert =maximum number of ipert mpi_enreg=MPI-parallelisation information mpsang= 1+maximum angular momentum for nonlocal pseudopotentials mpw = maximum number of planewaves in basis sphere (large number) natom = number of atoms in unit cell nattyp(ntypat)= # atoms of each type. nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90) nfftotf=total number of real space fine grid points ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90) nkpt = number of k points nk3xc=second dimension of the array k3xc nspden = number of spin-density components nspinor = number of spinorial components of the wavefunctions nsppol = number of channels for spin-polarization (1 or 2) nsym1=number of symmetry elements in space group consistent with the perturbation npwarr(nkpt) = array holding npw for each k point occ(mband*nkpt*nsppol) = occupation number for each band and k pawang <type(pawang_type)>=paw angular mesh and related data pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data pawrhoij0(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS pawrhoij1_i1pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i1pert) pawrhoij1_i2pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i2pert) pawrhoij1_i3pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i3pert) paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh paw_an1_i2pert(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh (i2pert) paw_ij1_i2pert(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels (i2pert) ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information psps <type(pseudopotential_type)> = variables related to pseudopotentials rho1r1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i1pert) rho1r2(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i2pert) rho1r3(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i3pert) rprimd(3,3) = dimensional primitive translations (bohr) symaf1(nsym)=(anti)ferromagnetic part of symmetry operations symrc1(3,3,nsym)=symmetries of group in terms of operations on reciprocal space primitive translations ucvol=volume of the unit cell vtrial(nfftf,nspden)=GS Vtrial(r). vhartr1_i2pert(cplex*nfftf,nspden)=firs-order hartree potential vtrial1_i2pert(cplex*nfft,nspden)=firs-order local potential vxc1_i2pert(cplex*nfft,nspden)=firs-order exchange-correlation potential ddk_f = wf files xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc (i1pert) xccc3d2(cplex*n3xccc)=3D change in core charge density, see n3xccc (i2pert) xccc3d3(cplex*n3xccc)=3D change in core charge density, see n3xccc (i3pert) xred(3,natom) = reduced atomic coordinates
OUTPUT
d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor = \sum_{i=1}^9 d3etot_i d3etot_1(2,3,mpert,3,mpert,3,mpert) = 1st term of d3etot d3etot_2(2,3,mpert,3,mpert,3,mpert) = 2nd term of d3etot d3etot_3(2,3,mpert,3,mpert,3,mpert) = 3rd term of d3etot d3etot_4(2,3,mpert,3,mpert,3,mpert) = 4th term of d3etot d3etot_5(2,3,mpert,3,mpert,3,mpert) = 5th term of d3etot d3etot_6(2,3,mpert,3,mpert,3,mpert) = 6th term of d3etot d3etot_7(2,3,mpert,3,mpert,3,mpert) = 7th term of d3etot d3etot_8(2,3,mpert,3,mpert,3,mpert) = 8th term of d3etot d3etot_9(2,3,mpert,3,mpert,3,mpert) = 9th term of d3etot
SIDE EFFECTS
TO DO!
SOURCE
183 subroutine dfptnl_pert(atindx,cg,cg1,cg2,cg3,cplex,dtfil,dtset,d3etot,eigen0,gs_hamkq,k3xc,indsy1,i1dir,i2dir,i3dir,& 184 & i1pert,i2pert,i3pert,kg,mband,mgfft,mkmem_rbz,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,nattyp,nfftf,nfftotf,ngfftf,nkpt,nk3xc,& 185 & nspden,nspinor,nsppol,nsym1,npwarr,occ,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawtab,& 186 & pawrhoij0,pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,& 187 & paw_an0,paw_an1_i2pert,paw_ij1_i2pert,ph1d,psps,rho1r1,rho2r1,rho3r1,rprimd,symaf1,symrc1,& 188 & ucvol,vtrial,vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert,ddk_f,xccc3d1,xccc3d2,xccc3d3,xred,& 189 & d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9) 190 191 !Arguments ------------------------------------ 192 !scalars 193 integer,intent(in) :: cplex,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband,mgfft 194 integer,intent(in) :: mk1mem,mkmem_rbz,mpert,mpsang,mpw,natom,nfftf,nfftotf,nkpt,nspden,nsym1 195 integer,intent(in) :: nk3xc,nspinor,nsppol 196 real(dp),intent(in) :: ucvol 197 type(MPI_type),intent(inout) :: mpi_enreg 198 type(datafiles_type),intent(in) :: dtfil 199 type(dataset_type),intent(in) :: dtset 200 type(pseudopotential_type),intent(in) :: psps 201 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 202 type(pawang_type),intent(inout) :: pawang,pawang1 203 type(pawfgr_type),intent(in) :: pawfgr 204 type(wfk_t),intent(inout) :: ddk_f(5) 205 206 !arrays 207 integer,intent(in) :: atindx(natom),kg(3,mpw*mkmem_rbz),nattyp(psps%ntypat),ngfftf(18),npwarr(nkpt) 208 integer,intent(in) :: indsy1(4,nsym1,dtset%natom),symaf1(nsym1),symrc1(3,3,nsym1) 209 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) 210 real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol) 211 real(dp),intent(in) :: cg2(2,mpw*nspinor*mband*mk1mem*nsppol) 212 real(dp),intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol) 213 real(dp),intent(in) :: eigen0(dtset%mband*dtset%nkpt*dtset%nsppol) 214 real(dp),intent(in) :: k3xc(nfftf,nk3xc) 215 real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom) 216 real(dp),intent(in) :: rho1r1(cplex*nfftf,dtset%nspden),rho2r1(cplex*nfftf,dtset%nspden) 217 real(dp),intent(in) :: rho3r1(cplex*nfftf,dtset%nspden),rprimd(3,3) 218 real(dp),intent(in) :: vtrial(cplex*nfftf,nspden) 219 real(dp),intent(in) :: xccc3d1(cplex*nfftf),xccc3d2(cplex*nfftf),xccc3d3(cplex*nfftf),xred(3,natom) 220 real(dp),intent(in) :: vxc1_i2pert(cplex*nfftf,nspden),vhartr1_i2pert(cplex*nfftf) 221 real(dp),intent(inout) :: vtrial1_i2pert(cplex*nfftf,nspden),d3etot(2,3,mpert,3,mpert,3,mpert) 222 real(dp),intent(inout) :: d3etot_1(2,3,mpert,3,mpert,3,mpert) 223 real(dp),intent(inout) :: d3etot_2(2,3,mpert,3,mpert,3,mpert) 224 real(dp),intent(inout) :: d3etot_3(2,3,mpert,3,mpert,3,mpert) 225 real(dp),intent(inout) :: d3etot_4(2,3,mpert,3,mpert,3,mpert) 226 real(dp),intent(inout) :: d3etot_5(2,3,mpert,3,mpert,3,mpert) 227 real(dp),intent(inout) :: d3etot_6(2,3,mpert,3,mpert,3,mpert) 228 real(dp),intent(inout) :: d3etot_7(2,3,mpert,3,mpert,3,mpert) 229 real(dp),intent(inout) :: d3etot_8(2,3,mpert,3,mpert,3,mpert) 230 real(dp),intent(inout) :: d3etot_9(2,3,mpert,3,mpert,3,mpert) 231 type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom*psps%usepaw) 232 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 233 type(pawrhoij_type),intent(in) :: pawrhoij0(natom*psps%usepaw) 234 type(pawrhoij_type),intent(in),target :: pawrhoij1_i1pert(natom*psps%usepaw) 235 type(pawrhoij_type),intent(in) :: pawrhoij1_i2pert(natom*psps%usepaw) 236 type(pawrhoij_type),intent(in),target :: pawrhoij1_i3pert(natom*psps%usepaw) 237 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 238 type(paw_an_type),intent(in) :: paw_an0(natom*psps%usepaw) 239 type(paw_an_type),intent(inout) :: paw_an1_i2pert(natom*psps%usepaw) 240 type(paw_ij_type),intent(inout) :: paw_ij1_i2pert(natom*psps%usepaw) 241 242 !Local variables------------------------------- 243 !scalars 244 logical :: has_cprj_jband,compute_conjugate,compute_rho21 245 integer,parameter :: level=52,tim_nonlop=0 246 integer :: bandtot,choice,counter,cplex_cprj,cplex_loc,cplex_rhoij,cpopt,dimffnl1,iband,icg0,ider,ierr 247 integer :: idir0,idir_getgh2c,idir_phon,idir_elfd,ipert_phon,ipert_elfd 248 integer :: ia,iatm,ibg,ii,ikg,ikg1,ikpt,ilm,isppol,istwf_k,jband 249 integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nkpg1,nnlout,nsp,nspden_rhoij,npert_phon,npw_k,npw1_k,nzlmopt 250 integer :: offset_cgi,offset_cgj,offset_eig0,option,paw_opt,qphase_rhoij,debug_mode 251 integer :: signs,size_wf,size_cprj,spaceComm,typat_ipert_phon,usepaw,useylmgr1 252 real(dp) :: arg,dot1i,dot1r,dot2i,dot2r,doti,dotr,e3tot,lagi,lagi_paw,lagr,lagr_paw 253 real(dp) :: sumi,sum_psi1H1psi1,sum_psi1H1psi1_i 254 real(dp) :: sum_lambda1psi1psi1,sum_lambda1psi1psi1_i 255 real(dp) :: sum_psi0H2psi1a,sum_psi0H2psi1a_i,sum_psi0H2psi1b,sum_psi0H2psi1b_i 256 real(dp) :: sum_lambda1psi0S1psi1,sum_lambda1psi0S1psi1_i 257 character(len=1000) :: msg 258 !arrays 259 integer,allocatable :: kg_k(:,:),kg1_k(:,:) 260 real(dp) :: buffer(10),eHxc21_paw(2),eHxc21_nhat(2),exc3(2),exc3_paw(2),kpt(3),eig0_k(mband) 261 real(dp) :: enlout1(2),enlout2(2) 262 real(dp) :: rmet(3,3),wtk_k 263 real(dp),allocatable :: cgi(:,:),cgj(:,:),cg_jband(:,:,:),cwavef1(:,:),cwavef2(:,:),cwavef3(:,:),dkinpw(:) 264 real(dp),allocatable :: eig1_k_i2pert(:),eig1_k_stored(:) 265 real(dp),allocatable :: chi_ij(:,:,:,:),cwave_right(:,:),cwave_left(:,:),dudk(:,:),dudkde(:,:),dummy_array(:),dummy_array2(:,:) 266 real(dp),allocatable :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:) 267 real(dp),allocatable :: h_cwave(:,:),iddk(:,:),kinpw1(:),kpg_k(:,:),kpg1_k(:,:),nhat21(:,:),occ_k(:) 268 real(dp),allocatable :: phkxred(:,:),ph3d(:,:,:),s_cwave(:,:) 269 real(dp),allocatable :: vlocal(:,:,:,:),vlocal1_i2pert(:,:,:,:),v_i2pert(:,:),wfraug(:,:,:,:) 270 real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:) 271 real(dp),allocatable :: ylm_k(:,:),ylm1_k(:,:),ylmgr1_k(:,:,:) 272 type(pawcprj_type),allocatable :: cwaveprj0(:,:),cwaveprj1(:,:) 273 type(pawcprj_type),target :: cprj_empty(0,0) 274 type(pawcprj_type),allocatable,target :: cprj_jband(:,:) 275 type(pawrhoij_type),allocatable,target :: pawrhoij21(:) 276 type(pawrhoij_type),pointer :: pawrhoij21_unsym(:),pawrhoij11(:) 277 type(paw_ij_type),allocatable :: paw_ij_tmp(:) 278 type(rf_hamiltonian_type) :: rf_hamkq_i2pert 279 280 !*********************************************************************** 281 282 DBG_ENTER("COLL") 283 284 ABI_UNUSED(dtfil%ireadwf) 285 286 287 me = mpi_enreg%me 288 spaceComm=mpi_enreg%comm_cell 289 290 npert_phon = 0 291 if(i1pert<=natom) npert_phon = npert_phon + 1 292 if(i2pert<=natom) npert_phon = npert_phon + 1 293 if(i3pert<=natom) npert_phon = npert_phon + 1 294 if (npert_phon>1) then 295 ABI_ERROR("dfptnl_pert is available with at most one phonon perturbation. Change your input!") 296 end if 297 298 usepaw = psps%usepaw 299 size_cprj = nspinor 300 301 call init_rf_hamiltonian(cplex,gs_hamkq,i2pert,rf_hamkq_i2pert,paw_ij1=paw_ij1_i2pert,has_e1kbsc=.true.) 302 303 ABI_MALLOC(dummy_array,(0)) 304 ABI_MALLOC(dummy_array2,(0,0)) 305 306 !Acivate computation of rho^(2:1) and related energy derivatives if needed 307 compute_rho21 = .false. 308 if (usepaw==1.and.npert_phon==1.and.(i1pert<=natom.or.i3pert<=natom)) then ! so i2pert==natom+2 309 compute_rho21 = .true. 310 if (i1pert<=natom) then 311 ipert_phon = i1pert 312 idir_phon = i1dir 313 ipert_elfd = i3pert 314 idir_elfd = i3dir 315 pawrhoij11 => pawrhoij1_i3pert 316 else if (i3pert<=natom) then 317 ipert_phon = i3pert 318 idir_phon = i3dir 319 ipert_elfd = i1pert 320 idir_elfd = i1dir 321 pawrhoij11 => pawrhoij1_i1pert 322 end if 323 ABI_MALLOC(pawrhoij21,(natom)) 324 call pawrhoij_nullify(pawrhoij21) 325 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 326 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 327 call pawrhoij_alloc(pawrhoij21,cplex_rhoij,nspden_rhoij,nspinor,dtset%nsppol,dtset%typat,& 328 & qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 329 ABI_MALLOC(cwaveprj0,(natom,size_cprj)) 330 ABI_MALLOC(cwaveprj1,(natom,size_cprj)) 331 call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj) 332 call pawcprj_alloc(cwaveprj1,1,gs_hamkq%dimcprj) 333 ! if (paral_atom) then 334 ! ABI_MALLOC(pawrhoij1_unsym,(natom)) 335 ! call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 336 !& nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 337 ! call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,& 338 !& dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1) 339 ! else 340 pawrhoij21_unsym => pawrhoij21 341 call pawrhoij_init_unpacked(pawrhoij21_unsym) 342 ! end if 343 ! Compute phkxred : 344 ABI_MALLOC(phkxred,(2,natom)) 345 do ia=1,natom 346 iatm=min(atindx(ia),natom) 347 arg=two_pi*(kpt(1)*xred(1,ia)+kpt(2)*xred(2,ia)+kpt(3)*xred(3,ia)) 348 phkxred(1,iatm)=cos(arg);phkxred(2,iatm)=sin(arg) 349 end do 350 ABI_MALLOC(paw_ij_tmp,(natom)) 351 call paw_ij_nullify(paw_ij_tmp) 352 cplex_loc=1;nsp=1 ! Force nsppol/nspden to 1 because Dij^(1) due to electric field is spin-independent 353 call paw_ij_init(paw_ij_tmp,cplex_loc,dtset%nspinor,nsp,nsp,dtset%pawspnorb,natom,psps%ntypat,& 354 & dtset%typat,pawtab,has_dijfr=1,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 355 call paw_ij_reset_flags(paw_ij_tmp,all=.True.) 356 call pawdijfr(gs_hamkq%gprimd,idir_elfd,ipert_elfd,natom,natom,nfftf,ngfftf,nsp,nsp,psps%ntypat,& 357 & 1,paw_ij_tmp,pawang,pawfgrtab,pawrad,pawtab,cplex_loc,& 358 & (/zero,zero,zero/),rprimd,ucvol,dummy_array2,dummy_array2,dummy_array2,xred,& 359 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 360 ABI_MALLOC(chi_ij,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,dtset%nspinor**2,cplex_loc)) 361 call pawdij2e1kb(paw_ij_tmp,1,mpi_enreg%comm_atom,mpi_enreg%my_atmtab,e1kbfr=chi_ij) 362 call paw_ij_free(paw_ij_tmp) 363 ABI_FREE(paw_ij_tmp) 364 else 365 ABI_MALLOC(chi_ij,(0,0,0,0)) 366 ABI_MALLOC(phkxred,(0,0)) 367 ABI_MALLOC(pawrhoij21,(0)) 368 pawrhoij21_unsym => pawrhoij21 369 ABI_MALLOC(cwaveprj0,(0,0)) 370 ABI_MALLOC(cwaveprj1,(0,0)) 371 end if 372 373 nnlout = 0 374 bandtot = 0 375 icg0 = 0 376 option = 2 377 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 378 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 379 380 ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc)) 381 ABI_MALLOC(vlocal1_i2pert,(cplex*n4,n5,n6,gs_hamkq%nvloc)) 382 383 ABI_MALLOC(wfraug,(2,n4,n5,n6)) 384 385 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 386 387 sumi = zero 388 389 !Set up the Ylm for each k point 390 ABI_MALLOC(ylm,(dtset%mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)) 391 ABI_MALLOC(ylmgr,(dtset%mpw*mkmem_rbz,9,psps%mpsang*psps%mpsang*psps%useylm)) 392 if (psps%useylm==1) then 393 option=2 394 call initylmg(gs_hamkq%gprimd,kg,dtset%kptns,mkmem_rbz,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,& 395 dtset%nkpt,npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr) 396 end if 397 398 !Set up the spherical harmonics (Ylm) at k+q 399 useylmgr1=0; option=0 400 if (psps%useylm==1.and. & 401 & (i2pert==natom+1.or.i2pert==natom+3.or.i2pert==natom+4.or.(usepaw==1.and.i2pert==natom+2))) then 402 useylmgr1=1; option=1 403 end if 404 ABI_MALLOC(ylm1,(dtset%mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)) 405 ABI_MALLOC(ylmgr1,(dtset%mpw*mkmem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 406 !To change the following when q/=0 407 if (psps%useylm==1) then 408 call initylmg(gs_hamkq%gprimd,kg,dtset%kptns,mkmem_rbz,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,& 409 dtset%nkpt,npwarr,dtset%nsppol,option,rprimd,ylm1,ylmgr1) 410 end if 411 412 debug_mode = 0 413 if (dtset%nonlinear_info>3.or.dtset%nonlinear_info==2) debug_mode = 1 414 415 !Real parts 416 sum_psi1H1psi1 = zero 417 sum_lambda1psi1psi1 = zero 418 sum_lambda1psi0S1psi1 = zero 419 sum_psi0H2psi1a = zero 420 sum_psi0H2psi1b = zero 421 !Imaginary parts 422 sum_psi1H1psi1_i = zero 423 sum_lambda1psi1psi1_i = zero 424 sum_lambda1psi0S1psi1_i = zero 425 sum_psi0H2psi1a_i = zero 426 sum_psi0H2psi1b_i = zero 427 428 compute_conjugate = .false. 429 !We have to compute < u^(ip1) | H^(ip2) | u^(ip3) > 430 !For some cases, we want to apply H^(ip2) on < u^(ip1) |, not on | u^(ip3) > (see below) 431 if (i3pert<=natom) then ! As npert_phon<=1, this implies that i1pert=natom+2 and i2pert=natom+2 432 compute_conjugate = .true. 433 end if 434 435 !Loop over spins 436 do isppol = 1, nsppol 437 438 ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 439 ! Same thing for vlocal from vtrial Also take into account the spin. 440 call rf_transgrid_and_pack(isppol,nspden,usepaw,cplex,nfftf,dtset%nfft,dtset%ngfft,& 441 & gs_hamkq%nvloc,pawfgr,mpi_enreg,vtrial,vtrial1_i2pert,vlocal,vlocal1_i2pert) 442 443 ! Continue to initialize the Hamiltonian 444 call gs_hamkq%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.) 445 call rf_hamkq_i2pert%load_spin(isppol,vlocal1=vlocal1_i2pert,with_nonlocal=.true.) 446 447 ! Loop over k-points 448 449 ikg = 0 450 ikg1 = 0 451 452 do ikpt = 1, nkpt 453 454 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,mpi_enreg%me)) then 455 cycle ! Skip the rest of the k-point loop 456 end if 457 458 counter = 100*ikpt 459 460 nband_k = dtset%nband(ikpt+(isppol-1)*nkpt) 461 npw_k = npwarr(ikpt) 462 npw1_k = npw_k ! To change for q/=0 463 istwf_k = dtset%istwfk(ikpt) 464 ABI_MALLOC(occ_k,(nband_k)) 465 occ_k(:) = occ(1+bandtot:nband_k+bandtot) 466 wtk_k = dtset%wtk(ikpt) 467 468 size_wf = nspinor*npw_k 469 470 kpt(:) = dtset%kptns(:,ikpt) 471 472 ABI_MALLOC(cwavef1,(2,npw_k*nspinor)) 473 ABI_MALLOC(cwavef3,(2,npw_k*nspinor)) 474 if (compute_rho21) then 475 ABI_MALLOC(cwavef2,(2,npw_k*nspinor)) 476 end if 477 478 ABI_MALLOC(kg_k,(3,npw_k)) 479 ABI_MALLOC(kg1_k,(3,npw1_k)) 480 ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm)) 481 ABI_MALLOC(ylm1_k,(npw1_k,mpsang*mpsang*psps%useylm)) 482 ABI_MALLOC(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 483 484 ! Get (k+G) wave vectors and associated spherical harmonics 485 kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg) 486 if (psps%useylm==1) then 487 do ilm=1,mpsang*mpsang 488 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 489 end do 490 end if 491 492 ! Get (k+q+G) wave vectors and associated spherical harmonics 493 kg1_k(:,1:npw1_k)=kg(:,1+ikg1:npw1_k+ikg1) ! To change for q/=0 494 if (psps%useylm==1) then 495 do ilm=1,psps%mpsang*psps%mpsang 496 ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm) 497 end do 498 if (useylmgr1==1) then 499 do ilm=1,psps%mpsang*psps%mpsang 500 do ii=1,3 501 ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm) 502 end do 503 end do 504 end if 505 end if 506 507 ! Compute (k+G) vectors 508 nkpg=0;if(i2pert>=1.and.i2pert<=natom) nkpg=3*dtset%nloalg(3) 509 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 510 if (nkpg>0) then 511 call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k) 512 end if 513 514 ! Compute (k+q+G) vectors 515 nkpg1=0;if(i2pert>=1.and.i2pert<=natom) nkpg1=3*dtset%nloalg(3) 516 ABI_MALLOC(kpg1_k,(npw1_k,nkpg1)) 517 if (nkpg1>0) then 518 call mkkpg(kg1_k,kpg1_k,kpt,nkpg1,npw1_k) 519 end if 520 521 ! ===== Preparation of the non-local contributions 522 523 ! Compute nonlocal form factors ffnl1 at (k+q+G) 524 !-- Atomic displacement perturbation 525 if (i2pert<=natom) then 526 ider=0;idir0=0 527 !-- Electric field perturbation 528 else if (i2pert==natom+2) then 529 if (psps%usepaw==1) then 530 ider=1;idir0=i2dir 531 else 532 ider=0;idir0=0 533 end if 534 end if 535 if (compute_rho21) then ! compute_rho21 implies i2pert==natom+2 536 ider=1; idir0=4 537 end if 538 539 ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms 540 dimffnl1=1+ider 541 if (ider==1.and.(idir0==0.or.idir0==4)) dimffnl1=2+2*psps%useylm 542 ! if (ider==2.and.idir0==4) dimffnl1=3+7*psps%useylm 543 ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat)) 544 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,& 545 & psps%indlmn,kg1_k,kpg1_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,& 546 & npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k) 547 548 ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms 549 if (compute_rho21.and.debug_mode/=0) then 550 ABI_MALLOC(ffnl1_test,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat)) 551 idir0 = 0 ! for nonlop with signs = 1 552 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1_test,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,& 553 & psps%indlmn,kg1_k,kpg1_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,& 554 & npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k) 555 end if 556 557 ! ===== Preparation of the kinetic contributions 558 559 ! Note that not all these arrays should be allocated in the general case when wtk_k vanishes 560 561 ! Compute (1/2) (2 Pi)**2 (k+q+G)**2: 562 ABI_MALLOC(kinpw1,(npw1_k)) 563 kinpw1(:)=zero 564 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gs_hamkq%gmet,kg1_k,kinpw1,kpt,npw1_k,0,0) 565 566 ABI_MALLOC(dkinpw,(npw_k)) ! 1st derivative (1st direction) 567 dkinpw(:)=zero 568 569 !===== Load the k/k+q dependent parts of the Hamiltonian 570 571 ! Load k-dependent part in the Hamiltonian datastructure 572 ABI_MALLOC(ph3d,(2,npw_k,gs_hamkq%matblk)) 573 call gs_hamkq%load_k(kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,& 574 & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.) 575 call gs_hamkq%load_k(ffnl_k=ffnl1,kpt_k=kpt,npw_k=npw1_k,istwf_k=istwf_k,& 576 & kinpw_k=kinpw1,kg_k=kg1_k,kpg_k=kpg1_k,compute_gbound=.true.) 577 ! end if 578 579 ! Load k-dependent part in the 1st-order Hamiltonian datastructure 580 call rf_hamkq_i2pert%load_k(npw_k=npw_k,dkinpw_k=dkinpw) 581 582 ABI_MALLOC_OR_DIE(dudk, (2,nband_k*size_wf), ierr) 583 ABI_MALLOC_OR_DIE(dudkde,(2,nband_k*size_wf), ierr) 584 ABI_MALLOC_OR_DIE(eig1_k_i2pert,(2*nband_k), ierr) 585 ABI_MALLOC(eig1_k_stored,(2*nband_k**2)) 586 ABI_MALLOC(cgi,(2,size_wf)) 587 ABI_MALLOC(cwave_right,(2,size_wf)) 588 ABI_MALLOC(cwave_left,(2,size_wf)) 589 590 ! ************************************************************************************************** 591 ! Read dudk and dudkde 592 ! ************************************************************************************************** 593 594 do iband = 1,nband_k 595 if (occ_k(iband)>tol10) then 596 597 ! Read dude file 598 call ddk_f(1)%read_bks(iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert) 599 ! Copy eig1_k_i2pert in "eig1_k_stored" 600 eig1_k_stored(1+(iband-1)*2*nband_k:2*nband_k+(iband-1)*2*nband_k)=eig1_k_i2pert(:) 601 602 if (i2pert==natom+2) then 603 ! Read dudk file 604 call ddk_f(2)%read_bks(iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert) 605 offset_cgi = (iband-1)*size_wf+icg0 606 cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi) 607 ! Copy cwave_right in "dudk" 608 dudk(:,1+(iband-1)*size_wf:iband*size_wf)=cwave_right(:,:) 609 610 ! Read dudkde file 611 call ddk_f(3)%read_bks(iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert) 612 offset_cgi = (iband-1)*size_wf+icg0 613 cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi) 614 ! Copy cwave_right in "dudkde" 615 dudkde(:,1+(iband-1)*size_wf:iband*size_wf)=cwave_right(:,:) 616 end if 617 618 end if 619 end do 620 621 ABI_MALLOC(cgj,(2,size_wf)) 622 ABI_MALLOC(iddk,(2,size_wf)) 623 624 offset_eig0 = mband*(ikpt-1+nkpt*(isppol-1)) 625 eig0_k(:) = eigen0(1+offset_eig0:mband+offset_eig0) 626 627 ABI_MALLOC_OR_DIE(h_cwave,(2,size_wf), ierr) 628 ABI_MALLOC_OR_DIE(s_cwave,(2,size_wf), ierr) 629 630 ! Allocate work spaces when debug_mode is activated 631 has_cprj_jband=.false. 632 if (debug_mode/=0) then ! Only for test purposes 633 ABI_MALLOC(cg_jband,(2,size_wf*nband_k,2)) 634 cg_jband(:,:,1) = cg(:,1+icg0:size_wf*nband_k+icg0) 635 if (i2pert==natom+2) then ! Note the multiplication by "i" 636 cg_jband(1,:,2) = -dudk(2,1:size_wf*nband_k) 637 cg_jband(2,:,2) = dudk(1,1:size_wf*nband_k) 638 end if 639 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 640 ABI_MALLOC(cprj_jband,(natom,size_cprj*nband_k)) 641 has_cprj_jband=.true. 642 else 643 ABI_MALLOC(cprj_jband,(natom,0)) 644 end if 645 else 646 ABI_MALLOC(cg_jband,(2,0,2)) 647 ABI_MALLOC(cprj_jband,(natom,0)) 648 end if 649 650 ! Loop over bands 651 do jband = 1,nband_k 652 ! Skip bands not treated by current proc 653 if((mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me)) cycle 654 655 if (occ_k(jband)>tol10) then 656 657 ! tol_test = tol8 658 offset_cgj = (jband-1)*size_wf+icg0 659 cgj(:,:) = cg(:,1+offset_cgj:size_wf+offset_cgj) 660 661 ! ************************************************************************************************** 662 ! Compute < u^(1) | ( H^(1) - eps^(0) S^(1) ) | u^(1) > 663 ! ************************************************************************************************** 664 665 eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k) 666 cwavef1(:,:) = cg1(:,1+offset_cgj:size_wf+offset_cgj) 667 cwavef3(:,:) = cg3(:,1+offset_cgj:size_wf+offset_cgj) 668 if (i2pert==natom+2) then ! Note the multiplication by i 669 iddk(1,:) = -dudkde(2,1+(jband-1)*size_wf:jband*size_wf) 670 iddk(2,:) = dudkde(1,1+(jband-1)*size_wf:jband*size_wf) 671 else 672 iddk(:,:) = zero 673 end if 674 cwave_right(:,:) = cwavef3(:,:) 675 cwave_left(:,:) = cwavef1(:,:) 676 if (compute_conjugate) then 677 cwave_right(:,:) = cwavef1(:,:) 678 cwave_left(:,:) = cwavef3(:,:) 679 end if 680 681 ! Compute : < u^(ip1) | ( H^(ip2) - eps^(0) S^(ip2) ) | u^(ip3) > 682 ! or : < u^(ip3) | ( H^(ip2) - eps^(0) S^(ip2) ) | u^(ip1) > 683 call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_right,cprj_empty,h_cwave,s_cwave,eig0_k,eig1_k_i2pert,& 684 & jband,gs_hamkq,iddk,i2dir,i2pert,ikpt,isppol,mkmem_rbz,mpi_enreg,nband_k,nsppol,& 685 & debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf) 686 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_left,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 687 688 if (usepaw==1.and.i2pert/=natom+2) then ! S^(1) is zero for ipert=natom+2 689 call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwave_left,s_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 690 dotr = dotr - eig0_k(jband)*dot2r 691 doti = doti - eig0_k(jband)*dot2i 692 end if 693 if (compute_conjugate) doti = -doti 694 695 ! ************************************************************************************************** 696 ! Compute sum_i Lambda_ij^(1) < u_j^(1) | u_i^(1)> 697 ! ************************************************************************************************** 698 699 ! eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k) 700 lagr = zero ; lagi = zero 701 lagr_paw = zero ; lagi_paw = zero 702 703 do iband = 1, nband_k 704 if(occ_k(iband)>tol10) then 705 706 offset_cgi = (iband-1)*size_wf+icg0 707 cwavef3(:,:) = cg3(:,1+offset_cgi:size_wf+offset_cgi) 708 709 if(debug_mode/=0) then 710 eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k) 711 end if 712 713 ! Get Lambda_ij^(1) = < u_i^(0) | H^(1) - eps^(0) S^(1) | u_j^(0) > (see dfpt_cgwf.F90) 714 dot1r = eig1_k_i2pert(2*iband-1) 715 dot1i = eig1_k_i2pert(2*iband ) 716 717 ! Compute < u_j^(1) | S^(0) | u_i^(1) > 718 if (usepaw==1) then 719 ibg = 0 720 call getgsc(cwavef3,cprj_empty,gs_hamkq,s_cwave,ibg,0,0,ikpt,isppol,& 721 & size_wf,size_cprj,size_wf,mpi_enreg,natom,-1,npw_k,nspinor,select_k=KPRIME_H_KPRIME) 722 else 723 s_cwave(:,:) = cwavef3(:,:) 724 end if 725 726 call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwavef1,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft) 727 lagr = lagr + dot1r*dot2r - dot1i*dot2i 728 lagi = lagi + dot1r*dot2i + dot1i*dot2r 729 730 ! Compute < u_j^(0) | S^(1) | u_i^(1) > 731 if (usepaw==1.and.i1pert<=natom) then ! S^(1) is zero for ipert=natom+2 732 733 cpopt=-1+5*gs_hamkq%usecprj ; choice=2 ; signs=2 ; paw_opt=3 734 call nonlop(choice,cpopt,cprj_empty,dummy_array,gs_hamkq,i1dir,(/zero/),mpi_enreg,1,nnlout,& 735 & paw_opt,signs,s_cwave,tim_nonlop,cwavef3,dummy_array2,iatom_only=i1pert) 736 call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cgj,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft) 737 lagr_paw = lagr_paw + dot1r*dot2r - dot1i*dot2i 738 lagi_paw = lagi_paw + dot1r*dot2i + dot1i*dot2r 739 740 end if 741 742 ! Compute < u_j^(1) | S^(1) | u_i^(0) > 743 if (usepaw==1.and.i3pert<=natom) then ! S^(1) is zero for ipert=natom+2 744 745 cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi) 746 cpopt=-1+5*gs_hamkq%usecprj ; choice=2 ; signs=2 ; paw_opt=3 747 call nonlop(choice,cpopt,cprj_empty,dummy_array,gs_hamkq,i3dir,(/zero/),mpi_enreg,1,nnlout,& 748 & paw_opt,signs,s_cwave,tim_nonlop,cgi,dummy_array2,iatom_only=i3pert) 749 call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwavef1,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft) 750 lagr_paw = lagr_paw + dot1r*dot2r - dot1i*dot2i 751 lagi_paw = lagi_paw + dot1r*dot2i + dot1i*dot2r 752 753 end if 754 755 end if 756 end do ! iband 757 758 ! ************************************************************************************************** 759 ! Sum all band_by_band contributions 760 ! ************************************************************************************************** 761 762 ! Real part 763 sum_psi1H1psi1 = sum_psi1H1psi1 + dtset%wtk(ikpt)*occ_k(jband)*dotr 764 sum_lambda1psi1psi1 = sum_lambda1psi1psi1 - dtset%wtk(ikpt)*occ_k(jband)*lagr 765 sum_lambda1psi0S1psi1 = sum_lambda1psi0S1psi1 - dtset%wtk(ikpt)*occ_k(jband)*lagr_paw 766 767 ! Imaginary part 768 sum_psi1H1psi1_i = sum_psi1H1psi1_i + dtset%wtk(ikpt)*occ_k(jband)*doti 769 sum_lambda1psi1psi1_i = sum_lambda1psi1psi1_i - dtset%wtk(ikpt)*occ_k(jband)*lagi 770 sum_lambda1psi0S1psi1_i = sum_lambda1psi0S1psi1_i - dtset%wtk(ikpt)*occ_k(jband)*lagi_paw 771 772 ! ************************************************************************************************** 773 ! If compute_rho21 : accumulate rhoij and compute term with H_KV^(2) 774 ! ************************************************************************************************** 775 776 if (compute_rho21) then 777 778 if (i1pert<=natom) then ! If true, i3pert==natom+2 779 cwave_right = cg3(:,1+offset_cgj:size_wf+offset_cgj) 780 else if (i3pert<=natom) then ! If true, i1pert==natom+2 781 cwave_right = cg1(:,1+offset_cgj:size_wf+offset_cgj) 782 end if 783 choice = 2 784 cpopt = 0 785 call getcprj(choice,cpopt,cgj,cwaveprj0,& 786 & ffnl1,idir_phon,psps%indlmn,gs_hamkq%istwf_k,kg_k,kpg_k,kpt,psps%lmnmax,& 787 & mgfft,mpi_enreg,natom,nattyp,dtset%ngfft,dtset%nloalg,& 788 & npw_k,nspinor,psps%ntypat,phkxred,ph1d,ph3d,ucvol,psps%useylm) 789 call getcprj(choice,cpopt,cwave_right,cwaveprj1,& 790 & ffnl1,idir_phon,psps%indlmn,gs_hamkq%istwf_k,kg_k,kpg_k,kpt,psps%lmnmax,& 791 & mgfft,mpi_enreg,natom,nattyp,dtset%ngfft,dtset%nloalg,& 792 & npw_k,nspinor,psps%ntypat,phkxred,ph1d,ph3d,ucvol,psps%useylm) 793 794 cplex_cprj=2;if (gs_hamkq%istwf_k>1) cplex_cprj=1 795 call paw_dfptnl_accrhoij(atindx,cplex_cprj,cwaveprj0,cwaveprj0,cwaveprj1,cwaveprj1,i1pert,i3pert,isppol,natom,natom,& 796 & nspinor,occ_k(jband),pawrhoij21_unsym,wtk_k) 797 !& comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 798 799 ! Compute < psi^(0) | H_KV^(pert1pert3) | psi^(pert2) > + < psi^(pert2) | H_KV^(pert1pert3) | psi^(0) > 800 cwavef2(:,:) = cg2(:,1+offset_cgj:size_wf+offset_cgj) 801 802 ! Read dkk file (for tests only) 803 if (debug_mode/=0) then 804 if(idir_elfd==i2dir) then 805 call ddk_f(2)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk) 806 else 807 call ddk_f(4)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk) 808 end if 809 cg_jband(1,1+size_wf*(jband-1):size_wf*jband,2) = -iddk(2,1:size_wf) 810 cg_jband(2,1+size_wf*(jband-1):size_wf*jband,2) = iddk(1,1:size_wf) 811 end if 812 ! Read dkde file 813 if(idir_elfd==i2dir) then 814 call ddk_f(3)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk) 815 else 816 call ddk_f(5)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk) 817 end if 818 s_cwave = iddk 819 iddk(1,:) = -s_cwave(2,:) 820 iddk(2,:) = s_cwave(1,:) 821 call rf2_getidir(idir_phon,idir_elfd,idir_getgh2c) 822 call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwavef2,cprj_empty,s_cwave,dummy_array2,eig0_k,eig1_k_i2pert,& 823 & jband,gs_hamkq,iddk,idir_getgh2c,ipert_phon+natom+11,ikpt,isppol,mkmem_rbz,mpi_enreg,nband_k,nsppol,& 824 & debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf,enl=chi_ij,ffnl1=ffnl1,ffnl1_test=ffnl1_test) 825 call dotprod_g(enlout1(1),enlout1(2),gs_hamkq%istwf_k,npw_k*nspinor,2,cgj,s_cwave,& 826 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 827 sum_psi0H2psi1a = sum_psi0H2psi1a + half*dtset%wtk(ikpt)*occ_k(jband)*enlout1(1) ! be careful : factor 0.5 828 sum_psi0H2psi1a_i = sum_psi0H2psi1a_i + half*dtset%wtk(ikpt)*occ_k(jband)*enlout1(2) ! be careful : factor 0.5 829 830 ! Read ddk file 831 if(idir_elfd==i2dir) then 832 call ddk_f(2)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk) 833 else 834 call ddk_f(4)%read_bks(jband, ikpt, isppol, xmpio_single, cg_bks=iddk) 835 end if 836 s_cwave = iddk 837 iddk(1,:) = -s_cwave(2,:) 838 iddk(2,:) = s_cwave(1,:) 839 call rf2_apply_hamiltonian(cg_jband,cprj_jband,cgj,cprj_empty,s_cwave,dummy_array2,eig0_k,eig1_k_i2pert,& 840 & jband,gs_hamkq,iddk,idir_getgh2c,ipert_phon+natom+11,ikpt,isppol,mkmem_rbz,mpi_enreg,nband_k,nsppol,& 841 & debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf,enl=chi_ij,ffnl1=ffnl1,ffnl1_test=ffnl1_test) 842 call dotprod_g(enlout2(1),enlout2(2),gs_hamkq%istwf_k,npw_k*nspinor,2,cwavef2,s_cwave,& 843 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 844 sum_psi0H2psi1b = sum_psi0H2psi1b + half*dtset%wtk(ikpt)*occ_k(jband)*enlout2(1) ! be careful : factor 0.5 845 sum_psi0H2psi1b_i = sum_psi0H2psi1b_i + half*dtset%wtk(ikpt)*occ_k(jband)*enlout2(2) ! be careful : factor 0.5 846 847 end if ! end if compute_rho21 848 849 end if 850 end do ! end loop over jband 851 852 ! ************************************************************************************************** 853 ! END OF BAND LOOP 854 ! ************************************************************************************************** 855 856 ABI_FREE(cgi) 857 ABI_FREE(cgj) 858 ABI_FREE(iddk) 859 ABI_FREE(h_cwave) 860 ABI_FREE(s_cwave) 861 862 ABI_FREE(cwave_right) 863 ABI_FREE(cwave_left) 864 ABI_FREE(eig1_k_i2pert) 865 ABI_FREE(eig1_k_stored) 866 867 bandtot = bandtot + nband_k 868 icg0 = icg0 + npw_k*nspinor*nband_k 869 ikg = ikg + npw_k 870 ikg1 = ikg1 + npw1_k 871 872 ABI_FREE(cwavef1) 873 if (compute_rho21) then 874 ABI_FREE(cwavef2) 875 end if 876 ABI_FREE(cwavef3) 877 ABI_FREE(dkinpw) 878 ABI_FREE(kg_k) 879 ABI_FREE(kg1_k) 880 ABI_FREE(kinpw1) 881 ABI_FREE(dudk) 882 ABI_FREE(dudkde) 883 ABI_FREE(ylm_k) 884 ABI_FREE(ylm1_k) 885 ABI_FREE(ylmgr1_k) 886 ABI_FREE(ffnl1) 887 if (compute_rho21.and.debug_mode/=0) then 888 ABI_FREE(ffnl1_test) 889 end if 890 ABI_FREE(kpg_k) 891 ABI_FREE(kpg1_k) 892 ABI_FREE(cg_jband) 893 ABI_FREE(occ_k) 894 ABI_FREE(ph3d) 895 if (has_cprj_jband) call pawcprj_free(cprj_jband) 896 ABI_FREE(cprj_jband) 897 898 end do ! end loop over k-points 899 900 end do ! end loop over spins 901 902 call rf_hamkq_i2pert%free() 903 904 ! ************************************************************************************************** 905 ! GATHER BAND-BY-BAND AND XC CONTRIBUTIONS 906 ! ************************************************************************************************** 907 908 if (xmpi_paral == 1) then 909 910 ! Real parts 911 buffer(1) = sum_psi1H1psi1 912 buffer(2) = sum_lambda1psi1psi1 913 buffer(3) = sum_lambda1psi0S1psi1 914 buffer(4) = sum_psi0H2psi1a 915 buffer(5) = sum_psi0H2psi1b 916 917 ! Imaginary parts 918 buffer(6) = sum_psi1H1psi1_i 919 buffer(7) = sum_lambda1psi1psi1_i 920 buffer(8) = sum_lambda1psi0S1psi1_i 921 buffer(9) = sum_psi0H2psi1a_i 922 buffer(10) = sum_psi0H2psi1b_i 923 924 call xmpi_sum(buffer,spaceComm,ierr) 925 926 ! Real parts 927 sum_psi1H1psi1 = buffer(1) 928 sum_lambda1psi1psi1 = buffer(2) 929 sum_lambda1psi0S1psi1 = buffer(3) 930 sum_psi0H2psi1a = buffer(4) 931 sum_psi0H2psi1b = buffer(5) 932 933 ! Imaginary parts 934 sum_psi1H1psi1_i = buffer(6) 935 sum_lambda1psi1psi1_i = buffer(7) 936 sum_lambda1psi0S1psi1_i = buffer(8) 937 sum_psi0H2psi1a_i = buffer(9) 938 sum_psi0H2psi1b_i = buffer(10) 939 940 ! Accumulate PAW occupancies 941 if (compute_rho21) then 942 call pawrhoij_mpisum_unpacked(pawrhoij21_unsym,spaceComm) 943 end if 944 945 end if 946 947 ! ************************************************************************************************** 948 ! Compute E_xc^(3) (NOTE : E_H^(3) = 0) 949 ! ************************************************************************************************** 950 951 call dfptnl_exc3(cplex,exc3,k3xc,mpi_enreg,nk3xc,nfftf,nfftotf,nspden,rho1r1,rho2r1,rho3r1,ucvol,xccc3d1,xccc3d2,xccc3d3) 952 953 exc3_paw = zero 954 if (usepaw==1) then 955 956 call paw_dfptnl_energy(exc3_paw,dtset%ixc,natom,natom,psps%ntypat,paw_an0,pawang,dtset%pawprtvol,pawrad,& 957 & pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,pawtab,dtset%pawxcdev,mpi_enreg%my_atmtab,mpi_enreg%comm_atom) 958 959 end if 960 961 ! ************************************************************************************************** 962 ! Compute E_Hxc^(2:1) 963 ! ************************************************************************************************** 964 965 eHxc21_paw = zero 966 eHxc21_nhat = zero 967 if (compute_rho21) then 968 969 if (pawfgr%nfft/=nfftf) then 970 write(msg,'(2(a,i10))') 'pawfgr%nfft/=nfftf : pawfgr%nfft=',pawfgr%nfft,' nfftf = ',nfftf 971 ABI_ERROR(msg) 972 end if 973 974 call pawnhatfr(0,idir_phon,ipert_phon,natom,dtset%natom,nspden,psps%ntypat,& 975 & pawang,pawfgrtab,pawrhoij11,pawtab,rprimd) 976 977 ABI_MALLOC(nhat21,(cplex*nfftf,nspden)) 978 call pawmkrho(0,arg,cplex,gs_hamkq%gprimd,idir_phon,indsy1,ipert_phon,mpi_enreg,& 979 & natom,natom,nspden,nsym1,psps%ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,& 980 & dtset%pawprtvol,pawrhoij21,pawrhoij21_unsym,pawtab,dtset%qptn,dummy_array2,dummy_array2,& 981 & dummy_array2,rprimd,symaf1,symrc1,dtset%typat,ucvol,dtset%usewvl,xred,& 982 & pawang_sym=pawang1,pawnhat=nhat21,pawrhoij0=pawrhoij0) 983 984 ! if (paral_atom) then 985 ! call pawrhoij_free(pawrhoij21_unsym) 986 ! ABI_FREE(pawrhoij21_unsym) 987 ! end if 988 nzlmopt = 0 989 call pawdfptenergy(eHxc21_paw,i2pert,ipert_phon,dtset%ixc,natom,dtset%natom,dtset%ntypat,& 990 & nzlmopt,nzlmopt,paw_an0,paw_an1_i2pert,paw_ij1_i2pert,pawang,dtset%pawprtvol,& 991 & pawrad,pawrhoij1_i2pert,pawrhoij21,pawtab,dtset%pawxcdev,dtset%xclevel) 992 !& mpi_atmtab=my_atmtab,comm_atom=my_comm_atom 993 994 ABI_MALLOC(v_i2pert,(cplex*nfftf,nspden)) 995 v_i2pert(:,1) = vhartr1_i2pert(:) 996 if(nspden>1) then 997 v_i2pert(:,2) = vhartr1_i2pert(:) 998 end if 999 1000 typat_ipert_phon = dtset%typat(ipert_phon) 1001 if (debug_mode/=0) then 1002 write(msg,'(2(a,i6))') ' DFPTNL_PERT : pawtab(',typat_ipert_phon,')%usexcnhat = ',pawtab(ipert_phon)%usexcnhat 1003 call wrtout(std_out,msg,'COLL') 1004 end if 1005 if (pawtab(typat_ipert_phon)%usexcnhat>0) then 1006 v_i2pert(:,:) = v_i2pert(:,:) + vxc1_i2pert(:,:) 1007 end if 1008 1009 call dotprod_vn(cplex,nhat21,eHxc21_nhat(1),eHxc21_nhat(2),nfftf,nfftotf,nspden,2,v_i2pert,& 1010 & ucvol,mpi_comm_sphgrid=mpi_enreg%comm_fft) 1011 1012 ABI_FREE(v_i2pert) 1013 ABI_FREE(nhat21) 1014 1015 end if 1016 1017 ! ************************************************************************************************** 1018 ! ALL TERMS HAVE BEEN COMPUTED 1019 ! ************************************************************************************************** 1020 1021 !Real part 1022 e3tot = sum_psi1H1psi1 + sum_lambda1psi1psi1 + sum_lambda1psi0S1psi1 + sum_psi0H2psi1a + sum_psi0H2psi1b 1023 e3tot = e3tot + half * (eHxc21_paw(1)+eHxc21_nhat(1)) + sixth * (exc3(1) + exc3_paw(1)) 1024 1025 !Imaginary part 1026 sumi = sum_psi1H1psi1_i + sum_lambda1psi1psi1_i + sum_lambda1psi0S1psi1_i + sum_psi0H2psi1a_i + sum_psi0H2psi1b_i 1027 sumi = sumi + half * (eHxc21_paw(2)+eHxc21_nhat(2)) + sixth * (exc3(2) + exc3_paw(2)) 1028 1029 !Real parts 1030 d3etot(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = e3tot 1031 d3etot_1(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi1H1psi1 1032 d3etot_2(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi1psi1 1033 d3etot_3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi0S1psi1 1034 d3etot_4(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1a 1035 d3etot_5(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1b 1036 d3etot_6(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_paw(1) 1037 d3etot_7(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_nhat(1) 1038 d3etot_8(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3(1) 1039 d3etot_9(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3_paw(1) 1040 !Imaginary parts 1041 d3etot(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi 1042 d3etot_1(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi1H1psi1_i 1043 d3etot_2(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi1psi1_i 1044 d3etot_3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi0S1psi1_i 1045 d3etot_4(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1a_i 1046 d3etot_5(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1b_i 1047 d3etot_6(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_paw(2) 1048 d3etot_7(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_nhat(2) 1049 d3etot_8(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3(2) 1050 d3etot_9(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3_paw(2) 1051 1052 !Before printing, set small contributions to zero 1053 !Real parts 1054 if (abs(sum_psi1H1psi1) <tol8) sum_psi1H1psi1 = zero 1055 if (abs(sum_lambda1psi1psi1) <tol8) sum_lambda1psi1psi1 = zero 1056 if (abs(sum_lambda1psi0S1psi1)<tol8) sum_lambda1psi0S1psi1 = zero 1057 if (abs(sum_psi0H2psi1a) <tol8) sum_psi0H2psi1a = zero 1058 if (abs(sum_psi0H2psi1b) <tol8) sum_psi0H2psi1b = zero 1059 if (abs(eHxc21_paw(1)) <tol8) eHxc21_paw(1) = zero 1060 if (abs(eHxc21_nhat(1)) <tol8) eHxc21_nhat(1) = zero 1061 if (abs(exc3(1)) <tol8) exc3(1) = zero 1062 if (abs(exc3_paw(1)) <tol8) exc3_paw(1) = zero 1063 if (abs(e3tot) <tol8) e3tot = zero 1064 1065 !Imaginary parts 1066 if (abs(sum_psi1H1psi1_i) <tol8) sum_psi1H1psi1_i = zero 1067 if (abs(sum_lambda1psi1psi1_i) <tol8) sum_lambda1psi1psi1_i = zero 1068 if (abs(sum_lambda1psi0S1psi1_i)<tol8) sum_lambda1psi0S1psi1_i = zero 1069 if (abs(sum_psi0H2psi1a_i) <tol8) sum_psi0H2psi1a_i = zero 1070 if (abs(sum_psi0H2psi1b_i) <tol8) sum_psi0H2psi1b_i = zero 1071 if (abs(eHxc21_paw(2)) <tol8) eHxc21_paw(2) = zero 1072 if (abs(eHxc21_nhat(2)) <tol8) eHxc21_nhat(2) = zero 1073 if (abs(exc3(2)) <tol8) exc3(2) = zero 1074 if (abs(exc3_paw(2)) <tol8) exc3_paw(2) = zero 1075 if (abs(sumi) <tol8) sumi = zero 1076 1077 write(msg,'(2a,3(a,i2,a,i1))') ch10,'NONLINEAR : ',& 1078 ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir 1079 call wrtout(std_out,msg,'COLL') 1080 call wrtout(ab_out,msg,'COLL') 1081 if (dtset%nonlinear_info>0) then 1082 write(msg,'(10(a,2(a,f18.8)),a)') & 1083 ch10,' sum_psi1H1psi1 = ',sum_psi1H1psi1, ',',sum_psi1H1psi1_i,& 1084 ch10,' sum_lambda1psi1psi1 = ',sum_lambda1psi1psi1, ',',sum_lambda1psi1psi1_i,& 1085 ch10,' sum_lambda1psi0S1psi1 = ',sum_lambda1psi0S1psi1, ',',sum_lambda1psi0S1psi1_i,& 1086 ch10,' sum_psi0H2psi1a = ',sum_psi0H2psi1a, ',',sum_psi0H2psi1a_i,& 1087 ch10,' sum_psi0H2psi1b = ',sum_psi0H2psi1b, ',',sum_psi0H2psi1b_i,& 1088 ch10,' eHxc21_paw/2 = ',half*eHxc21_paw(1), ',',half*eHxc21_paw(2),& 1089 ch10,' eHxc21_nhat/2 = ',half*eHxc21_nhat(1), ',',half*eHxc21_nhat(2),& 1090 ch10,' exc3/6 = ',sixth*exc3(1), ',',sixth*exc3(2),& 1091 ch10,' exc3_paw/6 = ',sixth*exc3_paw(1), ',',sixth*exc3_paw(2),& 1092 ch10,' >>>>>>>>>>>>>>> e3tot = ',e3tot, ',',sumi,ch10 1093 call wrtout(std_out,msg,'COLL') 1094 call wrtout(ab_out,msg,'COLL') 1095 end if 1096 1097 if (compute_rho21) then 1098 call pawcprj_free(cwaveprj0) 1099 call pawcprj_free(cwaveprj1) 1100 call pawrhoij_free(pawrhoij21) 1101 end if 1102 ABI_FREE(chi_ij) 1103 ABI_FREE(dummy_array) 1104 ABI_FREE(dummy_array2) 1105 ABI_FREE(phkxred) 1106 ABI_FREE(cwaveprj0) 1107 ABI_FREE(cwaveprj1) 1108 ABI_FREE(pawrhoij21) 1109 ABI_FREE(ylm) 1110 ABI_FREE(ylm1) 1111 ABI_FREE(ylmgr) 1112 ABI_FREE(ylmgr1) 1113 ABI_FREE(vlocal) 1114 ABI_FREE(vlocal1_i2pert) 1115 ABI_FREE(wfraug) 1116 1117 DBG_EXIT("COLL") 1118 1119 end subroutine dfptnl_pert
ABINIT/m_dfptnl_pert [ Modules ]
NAME
m_dfptnl_pert
FUNCTION
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group () 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_dfptnl_pert 23 24 use defs_basis 25 use m_abicore 26 use m_wffile 27 use m_wfk 28 use m_xmpi 29 use m_hamiltonian 30 use m_errors 31 use m_rf2 32 use m_kg 33 use m_dtset 34 use m_dtfil 35 36 use defs_datatypes, only : pseudopotential_type 37 use defs_abitypes, only : MPI_type 38 use m_cgtools, only : dotprod_g,sqnorm_g,dotprod_vn 39 use m_pawang, only : pawang_type 40 use m_pawfgrtab, only : pawfgrtab_type 41 use m_pawrad, only : pawrad_type 42 use m_pawtab, only : pawtab_type 43 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 44 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify,paw_ij_reset_flags 45 use m_pawdij, only : pawdijfr 46 use m_pawfgr, only : pawfgr_type 47 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc , pawrhoij_nullify, pawrhoij_free,& 48 & pawrhoij_init_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_inquire_dim 49 use m_paw_an, only : paw_an_type 50 use m_paw_mkrho, only : pawmkrho 51 use m_paw_nhat, only : pawnhatfr 52 use m_paw_dfpt, only : pawdfptenergy 53 use m_paw_dfptnl, only : paw_dfptnl_accrhoij,paw_dfptnl_energy 54 use m_initylmg, only : initylmg 55 use m_mkffnl, only : mkffnl 56 use m_getghc, only : getgsc 57 use m_getgh1c, only : rf_transgrid_and_pack 58 use m_mpinfo, only : proc_distrb_cycle 59 use m_nonlop, only : nonlop 60 use m_fourier_interpol, only : transgrid 61 use m_cgprj, only : getcprj 62 63 implicit none 64 65 private