TABLE OF CONTENTS
ABINIT/dfpt_eltfrxc [ Functions ]
NAME
dfpt_eltfrxc
FUNCTION
Compute the 2nd derivatives of exchange-correlation energy with respect to all pairs of strain and strain-atomic displacement for the frozen wavefunction contribution to the elastic and internal strain tensors
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DRH, DCA, XG, GMR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
atindx(natom)=index table for atoms ordered by type dtset <type(dataset_type)>=all input variables for this dataset | natom=number of atoms in unit cell | nfft=(effective) number of FFT grid points (for this processor) | nspden=number of spin-density components | ntypat=number of types of atoms in cell. | typat(natom)=integer type for each atom in cell enxc=exchange and correlation energy (hartree) gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere kxc(nfft,nkxc)=exchange and correlation kernel mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid (ngfftf=ngfft for norm-conserving potential runs) nkxc=2nd dimension of kxc n1xccc=dimension of xccc1d ; 0 if no XC core correction is used n3xccc=dimension of xccc3d (0 if no core charge, nfft otherwise) nhat(nfft,nspden*nhatdim)= -PAW only- compensation density pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information psps <type(pseudopotential_type)>=variables related to pseudopotentials rhor(nfft,nspden)=electron density in r space (if spin polarized, array contains total density in first half and spin-up density in second half) (for non-collinear magnetism, first element: total density, 3 next ones: mx,my,mz) rprimd(3,3)=dimensional primitive translation vectors (bohr) usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc vxc(nfft,nspden)=xc potential (spin up in first half and spin down in second half if nspden=2) xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives, for each type of atom, from psp xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xred(3,natom)=reduced coordinates for atoms in unit cell
OUTPUT
eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the elastic tensor
SIDE EFFECTS
NOTES
Much of the code in versions of this routine prior to 4.4.5
has been transfered to its child eltxccore.
PARENTS
respfn
CHILDREN
atm2fft,dfpt_atm2fft,dfpt_mkcore,dfpt_mkvxcstr,dotprod_vn,eltxccore fourdp,metric,paw_spline,pawpsp_cg,pawrad_free,pawrad_init,pawtab_free pawtab_nullify,redgr,timab,xmpi_sum
SOURCE
77 #if defined HAVE_CONFIG_H 78 #include "config.h" 79 #endif 80 81 #include "abi_common.h" 82 83 subroutine dfpt_eltfrxc(atindx,dtset,eltfrxc,enxc,gsqcut,kxc,mpi_enreg,mgfft,& 84 & nattyp,nfft,ngfft,ngfftf,nhat,nkxc,n3xccc,pawtab,ph1d,psps,rhor,rprimd,& 85 & usexcnhat,vxc,xccc3d,xred) 86 87 use defs_basis 88 use defs_datatypes 89 use defs_abitypes 90 use m_profiling_abi 91 use m_xmpi 92 93 use m_geometry, only : metric 94 use m_cgtools, only : dotprod_vn 95 use m_pawtab, only : pawtab_type,pawtab_free,pawtab_nullify 96 use m_pawrad, only : pawrad_type,pawrad_init,pawrad_free 97 use m_pawpsp, only : pawpsp_cg 98 use m_paw_numeric, only : paw_spline 99 100 !This section has been created automatically by the script Abilint (TD). 101 !Do not modify the following lines by hand. 102 #undef ABI_FUNC 103 #define ABI_FUNC 'dfpt_eltfrxc' 104 use interfaces_18_timing 105 use interfaces_53_ffts 106 use interfaces_64_psp 107 use interfaces_72_response, except_this_one => dfpt_eltfrxc 108 !End of the abilint section 109 110 implicit none 111 112 !Arguments ------------------------------------ 113 !type 114 !scalars 115 integer,intent(in) :: mgfft,n3xccc,nfft,nkxc,usexcnhat 116 real(dp),intent(in) :: enxc,gsqcut 117 type(MPI_type),intent(in) :: mpi_enreg 118 type(dataset_type),intent(in) :: dtset 119 type(pseudopotential_type),intent(inout) :: psps 120 !arrays 121 integer,intent(in) :: atindx(dtset%natom),nattyp(dtset%ntypat),ngfft(18) 122 integer,intent(in) :: ngfftf(18) 123 real(dp),intent(in) :: nhat(nfft,dtset%nspden*psps%usepaw) 124 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),rprimd(3,3) 125 real(dp),intent(in) :: vxc(nfft,dtset%nspden),xccc3d(n3xccc) 126 real(dp),intent(in) :: xred(3,dtset%natom) 127 real(dp),intent(in),target :: rhor(nfft,dtset%nspden) 128 real(dp),intent(inout) :: kxc(nfft,nkxc) 129 real(dp),intent(out) :: eltfrxc(6+3*dtset%natom,6) 130 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 131 132 !Local variables------------------------------- 133 !scalars 134 integer,parameter :: mshift=401 135 integer :: cplex,fgga,ia,idir,ielt,ieltx,ierr,ifft,ii,ipert,is1,is2,ispden,ispden_c,jj,ka,kb 136 integer :: kd,kg,n1,n1xccc,n2,n3,n3xccc_loc,optatm,optdyfr,opteltfr,optgr 137 integer :: option,optn,optn2,optstr,optv 138 real(dp) :: d2eacc,d2ecdgs2,d2exdgs2,d2gsds1ds2,d2gstds1ds2,decdgs,dexdgs 139 real(dp) :: dgsds10,dgsds20,dgstds10,dgstds20,rstep,spnorm,tmp0,tmp0t 140 real(dp) :: ucvol,valuei,yp1,ypn 141 type(pawrad_type) :: core_mesh 142 !arrays 143 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 144 real(dp) :: corstr(6),dummy6(0),dummy_in(0,0) 145 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0) 146 real(dp) :: eltfrxc_test1(6+3*dtset%natom,6),eltfrxc_test2(6+3*dtset%natom,6) 147 real(dp) :: gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3),tsec(2) 148 real(dp) :: strn_dummy6(0), strv_dummy6(0) 149 real(dp),allocatable :: d2gm(:,:,:,:),dgm(:,:,:),eltfrxc_tmp(:,:) 150 real(dp),allocatable :: eltfrxc_tmp2(:,:),elt_work(:,:),rho0_redgr(:,:,:) 151 real(dp),allocatable :: vxc10(:,:),vxc10_core(:),vxc10_coreg(:,:) 152 real(dp),allocatable :: vxc1is_core(:),vxc1is_coreg(:,:),vxc_core(:) 153 real(dp),allocatable :: vxc_coreg(:,:),work(:),workgr(:,:),xccc1d(:,:,:) 154 real(dp),allocatable :: xccc3d1(:),xccc3d1_temp(:,:),xcccrc(:) 155 real(dp),pointer :: rhor_(:,:) 156 type(pawtab_type),allocatable :: pawtab_test(:) 157 158 ! ************************************************************************* 159 160 !Initialize variables 161 cplex=1 162 qphon(:)=zero 163 n1=ngfft(1) 164 n2=ngfft(2) 165 n3=ngfft(3) 166 167 n1xccc = psps%n1xccc 168 if(psps%usepaw==0)then 169 ABI_ALLOCATE(xcccrc,(dtset%ntypat)) 170 ABI_ALLOCATE(xccc1d,(n1xccc,6,dtset%ntypat)) 171 xcccrc = psps%xcccrc 172 xccc1d = psps%xccc1d 173 end if 174 175 if (usexcnhat==0.and.dtset%usepaw==1) then 176 ABI_ALLOCATE(rhor_,(nfft,dtset%nspden)) 177 rhor_(:,:) = rhor(:,:)-nhat(:,:) 178 else 179 rhor_ => rhor 180 end if 181 182 !HACK - should be fixed globally 183 if(n1xccc==0) then 184 n3xccc_loc=0 185 else 186 n3xccc_loc=n3xccc 187 end if 188 189 fgga=0 ; if(nkxc==7.or.nkxc==19) fgga=1 190 191 ABI_ALLOCATE(eltfrxc_tmp,(6+3*dtset%natom,6)) 192 ABI_ALLOCATE(eltfrxc_tmp2,(6+3*dtset%natom,6)) 193 ABI_ALLOCATE(vxc10,(nfft,dtset%nspden)) 194 ABI_ALLOCATE(xccc3d1,(cplex*nfft)) 195 196 if(n1xccc/=0) then 197 ABI_ALLOCATE(vxc_core,(nfft)) 198 ABI_ALLOCATE(vxc10_core,(nfft)) 199 ABI_ALLOCATE(vxc1is_core,(nfft)) 200 201 if(dtset%nspden==1) then 202 vxc_core(:)=vxc(:,1) 203 else 204 vxc_core(:)=0.5_dp*(vxc(:,1)+vxc(:,2)) 205 end if 206 end if 207 208 !Compute gmet, gprimd and ucvol from rprimd 209 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 210 211 !For GGA case, prepare quantities needed to evaluate contributions 212 !arising from the strain dependence of the gradient operator itself 213 214 if(fgga==1) then 215 ABI_ALLOCATE(rho0_redgr,(3,nfft,dtset%nspden)) 216 ABI_ALLOCATE(work,(nfft)) 217 ABI_ALLOCATE(workgr,(nfft,3)) 218 219 ! Set up metric tensor derivatives 220 ABI_ALLOCATE(dgm,(3,3,6)) 221 ABI_ALLOCATE(d2gm,(3,3,6,6)) 222 ! Loop over 2nd strain index 223 do is2=1,6 224 kg=idx(2*is2-1);kd=idx(2*is2) 225 do jj = 1,3 226 dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj)) 227 end do 228 229 ! Loop over 1st strain index 230 do is1=1,6 231 ka=idx(2*is1-1);kb=idx(2*is1) 232 d2gm(:,:,is1,is2)=0._dp 233 do jj = 1,3 234 if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 235 & +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj) 236 if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 237 & +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj) 238 if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 239 & +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj) 240 if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 241 & +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj) 242 end do 243 d2gm(:,:,is1,is2)=0.5_dp*d2gm(:,:,is1,is2) 244 end do 245 end do 246 247 ! Compute the reduced gradients of the zero-order charge density. 248 ! Note that in the spin-polarized case, we are computing the reduced 249 ! gradients of 2 X the spin-up or spin-down charge. This simplifies 250 ! subsequent code for the non-spin-polarized case. 251 if(dtset%nspden==1) then 252 work(:)=rhor_(:,1) 253 else 254 work(:)=2.0_dp*rhor_(:,2) 255 end if 256 if(n1xccc/=0) then 257 work(:)=work(:)+xccc3d(:) 258 end if 259 call redgr (work,workgr,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb) 260 do ifft=1,nfft 261 rho0_redgr(:,ifft,1)=workgr(ifft,:) 262 end do 263 if(dtset%nspden==2) then 264 work(:)=2.0_dp*(rhor_(:,1)-rhor_(:,2)) 265 if(n1xccc/=0) then 266 work(:)=work(:)+xccc3d(:) 267 end if 268 call redgr(work,workgr,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb) 269 do ifft=1,nfft 270 rho0_redgr(:,ifft,2)=workgr(ifft,:) 271 end do 272 end if 273 ABI_DEALLOCATE(work) 274 ABI_DEALLOCATE(workgr) 275 end if !GGA 276 277 278 !Null the elastic tensor accumulator 279 eltfrxc(:,:)=zero;eltfrxc_tmp(:,:)=zero;eltfrxc_tmp2(:,:) = zero 280 281 !Normalization factor 282 if(dtset%nspden==1) then 283 spnorm=one 284 else 285 spnorm=half 286 end if 287 288 !Big loop over 2nd strain index 289 do is2=1,6 290 291 ! Translate strain index as needed by dfpt_mkcore below. 292 if(is2<=3) then 293 ipert=dtset%natom+3 294 idir=is2 295 else 296 ipert=dtset%natom+4 297 idir=is2-3 298 end if 299 300 ! Generate first-order core charge for is2 strain if core charges are present. 301 if(n1xccc/=0)then 302 303 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 304 ! Calculation in Reciprocal space for paw or NC with nc_xccc_gspace 305 ABI_ALLOCATE(xccc3d1_temp,(cplex*nfft,1)) 306 xccc3d1_temp = zero 307 308 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,is2,ipert,& 309 & mgfft,psps%mqgrid_vl,dtset%natom,1,nfft,ngfftf,dtset%ntypat,& 310 & ph1d,psps%qgrid_vl,qphon,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 311 & atmrhor1=xccc3d1_temp,optn2_in=1,& 312 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 313 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 314 xccc3d1(:) = xccc3d1_temp(:,1) 315 ABI_DEALLOCATE(xccc3d1_temp) 316 317 else 318 ! Calculation in direct space for norm conserving: 319 call dfpt_mkcore(cplex,idir,ipert,dtset%natom,dtset%ntypat,n1,n1xccc,& 320 & n2,n3,qphon,rprimd,dtset%typat,ucvol,& 321 & xcccrc,xccc1d,xccc3d1,xred) 322 end if 323 else 324 xccc3d1(:)=zero 325 end if 326 327 ! Compute the first-order potentials. 328 ! Standard first-order potential for LDA and GGA with core charge 329 if(fgga==0 .or. (fgga==1 .and. n1xccc/=0)) then 330 option=0 331 call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,& 332 & dummy_in,nkxc,dtset%nspden,n3xccc_loc,option,mpi_enreg%paral_kgb,qphon,rhor,rhor,& 333 & rprimd,dtset%usepaw,usexcnhat,vxc10,xccc3d1) 334 if(n1xccc/=0)then 335 if(dtset%nspden==1) then 336 vxc10_core(:)=vxc10(:,1) 337 vxc1is_core(:)=vxc10(:,1) 338 else 339 vxc10_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2)) 340 vxc1is_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2)) 341 end if 342 end if 343 end if 344 345 ! For GGA, first-order potential with doubled gradient operator strain 346 ! derivative terms needed for elastic tensor but not internal strain. 347 if(fgga==1) then 348 option=2 349 call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,& 350 & dummy_in,nkxc,dtset%nspden,n3xccc_loc,option,mpi_enreg%paral_kgb,qphon,rhor,rhor,& 351 & rprimd,dtset%usepaw,usexcnhat,vxc10,xccc3d1) 352 if(n1xccc/=0)then 353 if(dtset%nspden==1) then 354 vxc10_core(:)=vxc10(:,1) 355 else 356 vxc10_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2)) 357 end if 358 end if 359 end if 360 361 362 ! Additional term for diagonal strains. 363 if(is2<=3) then 364 vxc10(:,:)=vxc10(:,:)+vxc(:,:) 365 if(n1xccc/=0) then 366 vxc10_core(:)=vxc10_core(:)+2.0_dp*vxc_core(:) 367 vxc1is_core(:)=vxc1is_core(:)+vxc_core(:) 368 end if 369 end if 370 371 ! For GGA, compute the contributions from the strain derivatives acting 372 ! on the gradient operators. 373 if(fgga==1) then 374 375 if (dtset%nspden==1) then 376 do ifft=1,nfft 377 ! Collect the needed derivatives of Exc. The factors introduced 378 ! deal with the difference between density as used here and 379 ! spin density as used with these kxc terms in other contexts. 380 dexdgs =half *kxc(ifft,2) 381 d2exdgs2=quarter*kxc(ifft,4) 382 ! Loop over 1st strain index 383 do is1=1,6 384 ! The notation here is .gs... for the derivatives of the squared- 385 ! gradient of (2X) each spin density, and .gst... for the total density. 386 dgsds10=zero;dgsds20=zero;d2gsds1ds2=zero 387 do jj=1,3 388 do ii=1,3 389 tmp0=rho0_redgr(ii,ifft,1)*rho0_redgr(jj,ifft,1) 390 dgsds10=dgsds10+dgm(ii,jj,is1)*tmp0 391 dgsds20=dgsds20+dgm(ii,jj,is2)*tmp0 392 d2gsds1ds2=d2gsds1ds2+d2gm(ii,jj,is1,is2)*tmp0 393 end do 394 end do 395 ! Volume derivative terms added 396 if(is1<=3) d2gsds1ds2=d2gsds1ds2+dgsds20 397 if(is2<=3) d2gsds1ds2=d2gsds1ds2+dgsds10 398 ! Add the gradient derivative terms to eltfrxc. 399 eltfrxc(is1,is2)=eltfrxc(is1,is2)+d2exdgs2*dgsds10*dgsds20+dexdgs*d2gsds1ds2 400 end do !is1 401 end do !ifft 402 403 else ! nspden==2 404 405 do ispden=1,dtset%nspden 406 ispden_c=dtset%nspden-ispden+1 407 408 do ifft=1,nfft 409 410 ! Collect the needed derivatives of Exc. The factors introduced 411 ! deal with the difference between density as used here and 412 ! spin density as used with these kxc terms in other contexts. 413 dexdgs =quarter *kxc(ifft,3+ispden) 414 d2exdgs2=quarter*eighth*kxc(ifft,7+ispden) 415 decdgs =eighth *kxc(ifft,10) 416 d2ecdgs2=eighth*eighth *kxc(ifft,13) 417 418 ! Loop over 1st strain index 419 do is1=1,6 420 421 ! The notation here is .gs... for the derivatives of the squared- 422 ! gradient of (2X) each spin density, and .gst... for the total 423 ! density. Note the hack that the the total density is given 424 ! by the same expression for either the non-polarized or spin- 425 ! polarized case, implemented with the "complementary" index ispden_c 426 ! in the expression for tmp0t below. 427 dgsds10=zero;dgsds20=zero;d2gsds1ds2=zero 428 dgstds10=zero;dgstds20=zero;d2gstds1ds2=zero 429 do jj=1,3 430 do ii=1,3 431 tmp0=rho0_redgr(ii,ifft,ispden)*rho0_redgr(jj,ifft,ispden) 432 tmp0t=(rho0_redgr(ii,ifft,ispden)+rho0_redgr(ii,ifft,ispden_c))& 433 & *(rho0_redgr(jj,ifft,ispden)+rho0_redgr(jj,ifft,ispden_c)) 434 dgsds10=dgsds10+dgm(ii,jj,is1)*tmp0 435 dgsds20=dgsds20+dgm(ii,jj,is2)*tmp0 436 dgstds10=dgstds10+dgm(ii,jj,is1)*tmp0t 437 dgstds20=dgstds20+dgm(ii,jj,is2)*tmp0t 438 d2gsds1ds2=d2gsds1ds2+d2gm(ii,jj,is1,is2)*tmp0 439 d2gstds1ds2=d2gstds1ds2+d2gm(ii,jj,is1,is2)*tmp0t 440 end do 441 end do 442 ! Volume derivative terms added 443 if(is1<=3) then 444 d2gsds1ds2=d2gsds1ds2+dgsds20 445 d2gstds1ds2=d2gstds1ds2+dgstds20 446 end if 447 if(is2<=3) then 448 d2gsds1ds2=d2gsds1ds2+dgsds10 449 d2gstds1ds2=d2gstds1ds2+dgstds10 450 end if 451 452 ! Add the gradient derivative terms to eltfrxc. 453 eltfrxc(is1,is2)=eltfrxc(is1,is2)+spnorm*& 454 & (d2exdgs2*(dgsds10 *dgsds20) + dexdgs*d2gsds1ds2& 455 & +d2ecdgs2*(dgstds10*dgstds20)+ decdgs*d2gstds1ds2) 456 457 end do !is1 458 end do !ifft 459 end do !ispden 460 461 end if ! nspden 462 463 end if !GGA 464 465 ! Compute valence electron 1st-order charge contributions. Recall that 466 ! the diagonal strain derivatives of the valence charge are minus the 467 ! zero-order density. The explicit symmetrization avoids the need 468 ! to store vxc10 for strain indices other than is2. 469 470 call dotprod_vn(1,rhor_,d2eacc,valuei,nfft,nfft,dtset%nspden,1,& 471 & vxc10,ucvol) 472 do is1=1,3 473 eltfrxc_tmp(is1,is2)=eltfrxc_tmp(is1,is2)-0.5_dp*d2eacc 474 eltfrxc_tmp(is2,is1)=eltfrxc_tmp(is2,is1)-0.5_dp*d2eacc 475 end do 476 477 ! Compute additional core contributions from is1 perturbation 478 ! Internal strain terms calculated here. 479 if(n1xccc/=0) then 480 481 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 482 ! Calculation in Reciprocal space for paw or NC with nc_xccc_gspace 483 optatm=0;optdyfr=0;optgr=0;optstr=0;optv=0;optn=n3xccc/nfft;optn2=1;opteltfr=1 484 ABI_ALLOCATE(vxc10_coreg,(2,nfft)) 485 ABI_ALLOCATE(vxc_coreg,(2,nfft)) 486 ABI_ALLOCATE(vxc1is_coreg,(2,nfft)) 487 488 vxc10_coreg(:,:)=zero;vxc10_coreg(:,:)=zero;vxc1is_coreg(:,:)=zero; 489 490 ! Fourier transform of Vxc_core/vxc10_core to use in atm2fft (reciprocal space calculation) 491 call fourdp(1,vxc10_coreg,vxc10_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 492 call fourdp(1,vxc_coreg,vxc_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 493 call fourdp(1,vxc1is_coreg,vxc1is_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 494 495 call atm2fft(atindx,dummy_out1,dummy_out2,dummy_out3,dummy_out4,eltfrxc_tmp2,dummy_in,gmet,gprimd,& 496 & dummy_out5,dummy_out6,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%ntypat,& 497 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,& 498 & dummy_in,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,vxc_coreg,vxc10_coreg,vxc1is_coreg,dtset%vprtrb,psps%vlspl,is2_in=is2,& 499 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 500 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 501 502 ! The indexing array atindx is used to reestablish the correct order of atoms 503 ABI_ALLOCATE(elt_work,(6+3*dtset%natom,6)) 504 elt_work(1:6,1:6)=eltfrxc_tmp2(1:6,1:6) 505 do ia=1,dtset%natom 506 ielt=7+3*(ia-1) 507 ieltx=7+3*(atindx(ia)-1) 508 elt_work(ielt:ielt+2,1:6)=eltfrxc_tmp2(ieltx:ieltx+2,1:6) 509 end do 510 eltfrxc_tmp2(:,:)=elt_work(:,:) 511 ABI_DEALLOCATE(elt_work) 512 513 514 ABI_DEALLOCATE(vxc10_coreg) 515 ABI_DEALLOCATE(vxc_coreg) 516 ABI_DEALLOCATE(vxc1is_coreg) 517 eltfrxc(:,:)= eltfrxc(:,:) + eltfrxc_tmp2(:,:) 518 519 else 520 521 call eltxccore(eltfrxc,is2,mpi_enreg%my_natom,dtset%natom,nfft,dtset%ntypat,& 522 & n1,n1xccc,n2,n3,rprimd,dtset%typat,ucvol,vxc_core,vxc10_core,vxc1is_core,& 523 & xcccrc,xccc1d,xred,mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 524 525 !DEBUG 526 ! TEST ZONE (DO NOT REMOVE) USE TO RECIPROCAL SPACE IN NC CASE 527 if (dtset%userid==567) then 528 eltfrxc_test1(:,is2)=zero;eltfrxc_test2(:,is2)=zero 529 call eltxccore(eltfrxc_test1,is2,mpi_enreg%my_natom,dtset%natom,nfft,dtset%ntypat,& 530 & n1,n1xccc,n2,n3,rprimd,dtset%typat,ucvol,vxc_core,vxc10_core,vxc1is_core,& 531 & xcccrc,xccc1d,xred,mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 532 ! if (is2==1) print*,"elt-frxc from eltxccore",is2,eltfrxc_test1(1,1)*ucvol/dble(nfft) 533 ABI_DATATYPE_ALLOCATE(pawtab_test,(dtset%ntypat)) 534 call pawtab_nullify(pawtab_test) 535 do jj=1,dtset%ntypat 536 pawtab_test(jj)%mqgrid=psps%mqgrid_vl 537 ABI_ALLOCATE(pawtab_test(jj)%tcorespl,(pawtab_test(jj)%mqgrid,2)) 538 rstep=xcccrc(jj)/dble(n1xccc-1) 539 call pawrad_init(mesh=core_mesh,mesh_size=n1xccc,mesh_type=1,rstep=rstep) 540 call pawpsp_cg(pawtab_test(jj)%dncdq0,pawtab_test(jj)%d2ncdq0,psps%mqgrid_vl,psps%qgrid_vl,& 541 & pawtab_test(jj)%tcorespl(:,1),core_mesh,xccc1d(:,1,jj),yp1,ypn) 542 call paw_spline(psps%qgrid_vl,pawtab_test(jj)%tcorespl(:,1),psps%mqgrid_vl,yp1,ypn,pawtab_test(jj)%tcorespl(:,2)) 543 ! if (is2==1) then 544 ! do ii=1,n1xccc;write(100+jj,*) (ii-1)*rstep,xccc1d(ii,1,jj);enddo 545 ! do ii=1,psps%mqgrid_vl;write(200+jj,*) psps%qgrid_vl(ii),pawtab_test(jj)%tcorespl(ii,1);enddo 546 ! end if 547 end do 548 ABI_ALLOCATE(vxc10_coreg,(2,nfft)) 549 ABI_ALLOCATE(vxc_coreg,(2,nfft)) 550 ABI_ALLOCATE(vxc1is_coreg,(2,nfft)) 551 vxc10_coreg(:,:)=zero;vxc10_coreg(:,:)=zero;vxc1is_coreg(:,:)=zero; 552 call fourdp(1,vxc10_coreg,vxc10_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 553 call fourdp(1,vxc_coreg,vxc_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 554 call fourdp(1,vxc1is_coreg,vxc1is_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 555 optatm=0;optdyfr=0;optgr=0;optstr=0;optv=0;optn=1;optn2=1;opteltfr=1;corstr=zero 556 call atm2fft(atindx,dummy_out1,dummy_out2,dummy_out3,dummy_out4,eltfrxc_test2,dummy_in,gmet,gprimd,& 557 & dummy_out5,dummy_out6,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%ntypat,& 558 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab_test,ph1d,psps%qgrid_vl,dtset%qprtrb,& 559 & dummy_in,corstr,dummy6,ucvol,psps%usepaw,vxc_coreg,vxc10_coreg,vxc1is_coreg,dtset%vprtrb,psps%vlspl,is2_in=is2) 560 ABI_DEALLOCATE(vxc10_coreg) 561 ABI_DEALLOCATE(vxc_coreg) 562 ABI_DEALLOCATE(vxc1is_coreg) 563 call pawrad_free(core_mesh) 564 call pawtab_free(pawtab_test) 565 ABI_DATATYPE_DEALLOCATE(pawtab_test) 566 eltfrxc(:,:)= eltfrxc(:,:)+eltfrxc_test2(:,:) 567 ! if (is2==1) print*,"cor-str from atm2fft",is2,corstr*ucvol 568 ! if (is2==1) print*,"elt-frxc from atm2fft ",is2,eltfrxc_test2(1,1) 569 end if 570 !DEBUG 571 572 end if 573 end if 574 575 ! Additional term for diagonal strains 576 if(is2<=3) then 577 do is1=1,3 578 eltfrxc_tmp(is1,is2)=eltfrxc_tmp(is1,is2)+enxc 579 end do 580 end if 581 end do !is2 outermost strain loop 582 583 !Accumulate eltfrxc accross processors 584 call timab(48,1,tsec) 585 call xmpi_sum(eltfrxc,mpi_enreg%comm_fft,ierr) 586 call timab(48,2,tsec) 587 588 !Normalize accumulated 2nd derivatives in NC case 589 if(psps%usepaw==1)then 590 eltfrxc(:,:)=eltfrxc_tmp(:,:)+eltfrxc 591 else 592 eltfrxc(:,:)=eltfrxc_tmp(:,:)+eltfrxc*ucvol/dble(nfft) 593 end if 594 595 ABI_DEALLOCATE(eltfrxc_tmp) 596 ABI_DEALLOCATE(eltfrxc_tmp2) 597 ABI_DEALLOCATE(vxc10) 598 ABI_DEALLOCATE(xccc3d1) 599 if(psps%usepaw==0)then 600 ABI_DEALLOCATE(xccc1d) 601 ABI_DEALLOCATE(xcccrc) 602 end if 603 if (usexcnhat==0.and.dtset%usepaw==1) then 604 ABI_DEALLOCATE(rhor_) 605 end if 606 607 if(n1xccc/=0) then 608 ABI_DEALLOCATE(vxc_core) 609 ABI_DEALLOCATE(vxc10_core) 610 ABI_DEALLOCATE(vxc1is_core) 611 end if 612 613 if(fgga==1) then 614 ABI_DEALLOCATE(rho0_redgr) 615 ABI_DEALLOCATE(dgm) 616 ABI_DEALLOCATE(d2gm) 617 end if 618 619 end subroutine dfpt_eltfrxc