TABLE OF CONTENTS
ABINIT/dfptnl_loop [ Functions ]
NAME
dfptnl_loop
FUNCTION
Loop over the perturbations j1, j2 and j3
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (MVeithen,MB) 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
cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions
cgindex(nkpt,nsppol) = for each k-point, cgindex tores the location
of the WF in the cg array
dtfil <type(datafiles_type)>=variables related to files
dtset <type(dataset_type)>=all input variables for this dataset
gmet(3,3)=reciprocal space metric tensor in bohr**-2
gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double
that of the basis sphere--appropriate for charge density rho(G),
Hartree potential, and pseudopotentials
kg(3,mpw*mkmem)=reduced planewave coordinates
kneigh(30,nkpt) = index of the neighbours of each k-point
kg_neigh(30,nkpt,3) = necessary to construct the vector joining a k-point
to its nearest neighbour in case of a single k-point,
a line of k-points or a plane of k-points.
kptindex(2,nkpt3)= index of the k-points in the reduced BZ
related to a k-point in the full BZ
kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ
kxc(nfft,nkxc)=exchange-correlation kernel
k3xc(nfft,nk3xc)=third-order exchange-correlation kernel
mband = maximum number of bands
mgfft = maximum single fft dimension
mkmem = Number of k points treated by this node.
mkmem_max = maximal number of k-points on each processor (MPI //)
mk1mem = Number of k points for first-order WF treated by this node.
mpert =maximum number of ipert
mpi_enreg=MPI-parallelisation information
mpw = maximum number of planewaves in basis sphere (large number)
mvwtk(30,nkpt) = weights to compute the finite difference ddk
natom = number of atoms in unit cell
nfft = (effective) number of FFT grid points (for this processor)
nkpt = number of k points
nkpt3 = number of k-points in the full BZ
nkxc=second dimension of the array kxc, see rhotoxc.f for a description
nneigh = total number of neighbours required to evaluate the finite
difference formula
nspinor = number of spinorial components of the wavefunctions
nsppol = number of channels for spin-polarization (1 or 2)
npwarr(nkpt) = array holding npw for each k point
occ(mband*nkpt*nsppol) = occupation number for each band and k
psps <type(pseudopotential_type)> = variables related to pseudopotentials
pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat
between k-points
rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
that have to be computed
1 -> element has to be computed explicitely
-1 -> use symmetry operations to obtain the corresponding element
rprimd(3,3)=dimensional primitive translations (bohr)
ucvol = unit cell volume (bohr^3)
xred(3,natom) = reduced atomic coordinates
OUTPUT
blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE
(=1 if computed)
d3lo(2,3,mpert,3,mpert,3,mpert) = matrix of the 3DTEs
SIDE EFFECTS
hdr <type(hdr_type)>=the header of wf, den and pot files
PARENTS
nonlinear
CHILDREN
appdig,dfpt_mkcore,dfpt_mkvxc,dfpt_vlocal,dfptnl_mv,dfptnl_resp dotprod_vn,fourdp,getph,hartre,hdr_free,initylmg,inwffil,read_rhor status,timab,wffclose,wrtout
SOURCE
86 #if defined HAVE_CONFIG_H 87 #include "config.h" 88 #endif 89 90 #include "abi_common.h" 91 92 93 subroutine dfptnl_loop(blkflg,cg,cgindex,dtfil,dtset,d3lo,& 94 & gmet,gprimd,gsqcut, & 95 & hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,mband,mgfft,mkmem,mkmem_max,mk1mem,& 96 & mpert,mpi_enreg,mpw,mvwtk,natom,nfft,nkpt,nkpt3,nkxc,nk3xc,nneigh,nspinor,nsppol,& 97 & npwarr,occ,psps,pwind,& 98 & rfpert,rprimd,ucvol,xred) 99 100 use defs_basis 101 use defs_datatypes 102 use defs_abitypes 103 use defs_wvltypes 104 use m_wffile 105 use m_profiling_abi 106 use m_hdr 107 108 use m_ioarr, only : read_rhor 109 use m_kg, only : getph 110 use m_cgtools, only : dotprod_vn 111 use m_pawrhoij, only : pawrhoij_type 112 113 !This section has been created automatically by the script Abilint (TD). 114 !Do not modify the following lines by hand. 115 #undef ABI_FUNC 116 #define ABI_FUNC 'dfptnl_loop' 117 use interfaces_14_hidewrite 118 use interfaces_18_timing 119 use interfaces_32_util 120 use interfaces_53_ffts 121 use interfaces_56_recipspace 122 use interfaces_56_xc 123 use interfaces_72_response 124 use interfaces_79_seqpar_mpi 125 !End of the abilint section 126 127 implicit none 128 129 !Arguments ------------------------------------ 130 !scalars 131 integer,intent(in) :: mband,mgfft,mk1mem,mkmem,mkmem_max,mpert,mpw,natom,nfft 132 integer,intent(in) :: nk3xc,nkpt,nkpt3,nkxc,nneigh,nspinor,nsppol 133 real(dp),intent(in) :: gsqcut,ucvol 134 type(MPI_type),intent(inout) :: mpi_enreg 135 type(datafiles_type),intent(in) :: dtfil 136 type(dataset_type),intent(in) :: dtset 137 type(hdr_type),intent(inout) :: hdr 138 type(pseudopotential_type),intent(in) :: psps 139 !arrays 140 integer,intent(in) :: cgindex(nkpt,nsppol),kg(3,mk1mem*mpw),kneigh(30,nkpt) 141 integer,intent(in) :: kg_neigh(30,nkpt,3) 142 integer,intent(in) :: kptindex(2,nkpt3),npwarr(nkpt),pwind(mpw,nneigh,mkmem) 143 integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert) 144 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) !vz_i 145 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3) 146 real(dp),intent(in) :: gprimd(3,3),k3xc(nfft,nk3xc),kpt3(3,nkpt3) 147 real(dp),intent(in) :: kxc(nfft,nkxc),mvwtk(30,nkpt),rprimd(3,3) 148 real(dp),intent(in) :: xred(3,natom) 149 real(dp),intent(inout) :: occ(mband*nkpt*nsppol) 150 real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert) !vz_i 151 152 !Local variables------------------------------- 153 !scalars 154 integer,parameter :: level=51 155 integer :: ask_accurate,counter,cplex,formeig,i1dir 156 integer :: i1pert,i2dir,i2pert,i3dir,i3pert,iatom,ierr,iexit,ifft,index,ir 157 integer :: ireadwf,itypat,mcg,mpsang,n1,n2,n3,n3xccc,nfftot,nspden,option,optorth 158 integer :: pert1case,pert2case,pert3case,rdwrpaw,timrev,comm_cell 159 real(dp) :: ecut_eff,exc3,valuei 160 character(len=500) :: message 161 character(len=fnlen) :: fiden1i,fiwf1i,fiwf3i 162 type(wffile_type) :: wff1,wff2,wfft1,wfft2 163 type(wvl_data) :: wvl 164 type(hdr_type) :: hdr_den 165 !arrays 166 integer,allocatable :: atindx(:),atindx1(:),nattyp(:) 167 real(dp) :: d3_berry(2,3),rho_dum(1),tsec(2),ylmgr_dum(1) 168 real(dp),allocatable :: cg1(:,:),cg3(:,:),eigen1(:),ph1d(:,:),rho1r1(:,:) 169 real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:),rho3r1(:,:),vhartr1(:) 170 real(dp),allocatable :: vpsp1(:),vtrial1(:,:),vxc1(:,:),work(:),xc_tmp(:,:) 171 real(dp),allocatable :: xccc3d1(:),xccc3d2(:),xccc3d3(:),ylm(:,:,:) 172 type(pawrhoij_type),allocatable :: rhoij_dum(:) 173 174 ! *********************************************************************** 175 176 call timab(502,1,tsec) 177 call status(0,dtfil%filstat,iexit,level,'enter ') 178 179 comm_cell = mpi_enreg%comm_cell 180 181 timrev = 1 182 cplex = 2 - timrev 183 nspden = dtset%nspden 184 ecut_eff = (dtset%ecut)*(dtset%dilatmx)**2 185 mpsang = psps%mpsang 186 optorth=1;if (psps%usepaw==1) optorth=0 187 188 ABI_ALLOCATE(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 189 ABI_ALLOCATE(cg3,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 190 ABI_ALLOCATE(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol)) 191 ABI_ALLOCATE(rho1r1,(cplex*nfft,dtset%nspden)) 192 ABI_ALLOCATE(rho2r1,(cplex*nfft,dtset%nspden)) 193 ABI_ALLOCATE(rho2g1,(2,nfft)) 194 ABI_ALLOCATE(rho3r1,(cplex*nfft,dtset%nspden)) 195 ABI_ALLOCATE(ylm,(2,dtset%mpw*dtset%mkmem,mpsang*mpsang*psps%useylm)) 196 197 ask_accurate=1 ; formeig = 1 ; ireadwf = 1 198 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 199 nfftot=n1*n2*n3 200 201 !Generate an index table of atoms, in order for them to be used 202 !type after type. 203 ABI_ALLOCATE(atindx,(natom)) 204 ABI_ALLOCATE(atindx1,(natom)) 205 ABI_ALLOCATE(nattyp,(psps%ntypat)) 206 index=1 207 do itypat=1,psps%ntypat 208 nattyp(itypat)=0 209 do iatom=1,natom 210 if(dtset%typat(iatom)==itypat)then 211 atindx(iatom)=index 212 atindx1(index)=iatom 213 index=index+1 214 nattyp(itypat)=nattyp(itypat)+1 215 end if 216 end do 217 end do 218 219 !Generate the 1-dimensional phases 220 ABI_ALLOCATE(ph1d,(2,3*(2*mgfft+1)*natom)) 221 call getph(atindx,natom,n1,n2,n3,ph1d,xred) 222 223 !Set up the Ylm for each k point 224 if (psps%useylm==1) then 225 call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,& 226 & dtset%mpw,dtset%nband,dtset%nkpt,& 227 & npwarr,dtset%nsppol,0,rprimd,ylm,ylmgr_dum) 228 end if 229 230 ABI_ALLOCATE(vpsp1,(cplex*nfft)) 231 ABI_ALLOCATE(xccc3d1,(cplex*nfft)) 232 ABI_ALLOCATE(xccc3d2,(cplex*nfft)) 233 ABI_ALLOCATE(xccc3d3,(cplex*nfft)) 234 ABI_ALLOCATE(vhartr1,(cplex*nfft)) 235 ABI_ALLOCATE(vxc1,(cplex*nfft,dtset%nspden)) 236 ABI_ALLOCATE(vtrial1,(cplex*nfft,dtset%nspden)) 237 238 !Loop over the perturbations j1, j2, j3 239 240 pert1case = 0 ; pert2case = 0 ; pert3case = 0 241 242 do i1pert = 1, mpert 243 do i1dir = 1, 3 244 245 if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then 246 247 pert1case = i1dir + (i1pert-1)*3 248 counter = pert1case 249 call appdig(pert1case,dtfil%fnamewff1,fiwf1i) 250 251 call status(counter,dtfil%filstat,iexit,level,'call inwffil ') 252 253 mcg=mpw*nspinor*mband*mkmem*nsppol 254 call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,& 255 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 256 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 257 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 258 & dtset%nsppol,dtset%nsym,& 259 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 260 & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl) 261 262 if (ireadwf==1) then 263 call WffClose (wff1,ierr) 264 end if 265 266 rho1r1(:,:) = 0._dp 267 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 268 rdwrpaw=0 269 call appdig(pert1case,dtfil%fildens1in,fiden1i) 270 call status(counter,dtfil%filstat,iexit,level,'call ioarr ') 271 272 call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho1r1, & 273 hdr_den, rhoij_dum, comm_cell, check_hdr=hdr) 274 call hdr_free(hdr_den) 275 end if 276 277 xccc3d1(:) = 0._dp 278 if ((psps%n1xccc/=0).and.(i1pert <= natom)) then 279 call status(counter,dtfil%filstat,iexit,level,'call dfpt_mkcore ') 280 call dfpt_mkcore(cplex,i1dir,i1pert,natom,psps%ntypat,n1,psps%n1xccc,& 281 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,& 282 & psps%xcccrc,psps%xccc1d,xccc3d1,xred) 283 end if ! psps%n1xccc/=0 284 285 do i3pert = 1, mpert 286 do i3dir = 1, 3 287 288 if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then 289 290 pert3case = i3dir + (i3pert-1)*3 291 counter = 100*pert3case + pert1case 292 call appdig(pert3case,dtfil%fnamewff1,fiwf3i) 293 294 call status(counter,dtfil%filstat,iexit,level,'call inwffil ') 295 mcg=mpw*nspinor*mband*mkmem*nsppol 296 call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,& 297 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 298 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 299 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 300 & dtset%nsppol,dtset%nsym,& 301 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 302 & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,& 303 & fiwf3i,wvl) 304 if (ireadwf==1) then 305 call WffClose (wff2,ierr) 306 end if 307 308 rho3r1(:,:) = 0._dp 309 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 310 rdwrpaw=0 311 call appdig(pert3case,dtfil%fildens1in,fiden1i) 312 call status(counter,dtfil%filstat,iexit,level,'call ioarr ') 313 314 call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho3r1, & 315 hdr_den, rhoij_dum, comm_cell, check_hdr=hdr) 316 call hdr_free(hdr_den) 317 end if 318 319 xccc3d3(:) = 0._dp 320 if ((psps%n1xccc/=0).and.(i3pert <= natom)) then 321 call status(counter,dtfil%filstat,iexit,level,'call dfpt_mkcore ') 322 call dfpt_mkcore(cplex,i3dir,i3pert,natom,psps%ntypat,n1,psps%n1xccc,& 323 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,& 324 & psps%xcccrc,psps%xccc1d,xccc3d3,xred) 325 end if ! psps%n1xccc/=0 326 327 do i2pert = 1, mpert 328 329 ! In case of electric field perturbation, evaluate the ddk 330 ! using the finite difference expression of 331 ! Marzari and Vanderbilt PRB 56, 12847 (1997). 332 333 d3_berry(:,:) = 0._dp 334 335 if ((i2pert==dtset%natom+2).and.& 336 & (maxval(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert)) == 1)) then 337 338 call timab(511,1,tsec) 339 call status(counter,dtfil%filstat,iexit,level,'call dfptnl_mv ') 340 call dfptnl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,& 341 & i1pert,i3pert,i1dir,i3dir,& 342 & kneigh,kg_neigh,kptindex,kpt3,mband,mkmem,mkmem_max,mk1mem,& 343 & mpi_enreg,mpw,mvwtk,natom,nkpt,nkpt3,nneigh,npwarr,nspinor,nsppol,pwind) 344 call timab(511,2,tsec) 345 346 end if 347 348 if (mpi_enreg%me == 0) then 349 350 if(sum(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert))>0)then 351 write(message,'(a,a,a,a,a,a)')ch10,ch10,& 352 & ' Decomposition of the third-order energy for the set of perturbations',ch10 353 call wrtout(std_out,message,'COLL') 354 call wrtout(ab_out,message,'COLL') 355 if (i1pert < natom + 1) then 356 write(message,'(a,i3,a,i3)') & 357 & ' j1 : displacement of atom ',i1pert,' along direction ', i1dir 358 end if 359 if (i1pert == dtset%natom + 2) then 360 write(message,'(a,i4)')' j1 : homogeneous electric field along direction ',i1dir 361 end if 362 call wrtout(std_out,message,'COLL') 363 call wrtout(ab_out,message,'COLL') 364 if (i3pert < natom + 1) then 365 write(message,'(a,i3,a,i3,a)') & 366 & ' j3 : displacement of atom ',i3pert,' along direction ', i3dir,ch10 367 end if 368 if (i3pert == dtset%natom + 2) then 369 write(message,'(a,i4,a)')' j3 : homogeneous electric field along direction ',i3dir,ch10 370 end if 371 call wrtout(std_out,message,'COLL') 372 call wrtout(ab_out,message,'COLL') 373 end if 374 375 end if ! mpi_enreg%me == 0 376 377 do i2dir = 1, 3 378 379 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 380 pert2case = i2dir + (i2pert-1)*3 381 382 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 383 384 ! Read the first-order densities from disk-files 385 rho2r1(:,:) = 0._dp ; rho2g1(:,:) = 0._dp 386 387 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 388 rdwrpaw=0 389 call appdig(pert2case,dtfil%fildens1in,fiden1i) 390 call status(counter,dtfil%filstat,iexit,level,'call ioarr ') 391 392 call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho2r1, & 393 hdr_den, rhoij_dum, comm_cell, check_hdr=hdr) 394 call hdr_free(hdr_den) 395 396 ! Compute up+down rho1(G) by fft 397 ABI_ALLOCATE(work,(cplex*nfft)) 398 work(:)=rho2r1(:,1) 399 call status(counter,dtfil%filstat,iexit,level,'call fourdp ') 400 call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfft,dtset%ngfft,dtset%paral_kgb,0) 401 ABI_DEALLOCATE(work) 402 403 end if 404 405 ! Compute first-order local potentials 406 ! (hartree, xc and pseudopotential) 407 408 n3xccc=0; if(psps%n1xccc/=0)n3xccc=nfft 409 xccc3d2(:)=0._dp ; vpsp1(:)=0._dp 410 411 if (i2pert <= natom) then 412 413 call status(counter,dtfil%filstat,iexit,level,'call dfpt_vlocal ') 414 call dfpt_vlocal(atindx,cplex,gmet,gsqcut,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,natom,& 415 & nattyp,nfft,dtset%ngfft,psps%ntypat,n1,n2,n3,dtset%paral_kgb,ph1d,psps%qgrid_vl,& 416 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 417 418 if (psps%n1xccc/=0) then 419 call status(counter,dtfil%filstat,iexit,level,'call dfpt_mkcore ') 420 call dfpt_mkcore(cplex,i2dir,i2pert,natom,psps%ntypat,n1,psps%n1xccc,& 421 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,& 422 & psps%xcccrc,psps%xccc1d,xccc3d2,xred) 423 end if ! psps%n1xccc/=0 424 425 end if ! i2pert <= natom 426 427 call status(counter,dtfil%filstat,iexit,level,'call hartre ') 428 call hartre(cplex,gsqcut,0,mpi_enreg,nfft,dtset%ngfft,dtset%paral_kgb,rho2g1,rprimd,vhartr1,qpt=dtset%qptn) 429 option=1 430 call status(counter,dtfil%filstat,iexit,level,'call dfpt_mkvxc ') 431 call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,dtset%ngfft,& 432 & rho_dum,0,rho_dum,0,nkxc,dtset%nspden,n3xccc,option,& 433 & dtset%paral_kgb,dtset%qptn,rho2r1,rprimd,0,vxc1,xccc3d2) 434 435 if(dtset%nsppol==1)then 436 if(cplex==1)then 437 do ir=1,nfft 438 vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1) 439 end do 440 else 441 do ir=1,nfft 442 vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1) 443 vtrial1(2*ir ,1)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,1) 444 end do 445 end if 446 else 447 if(cplex==1)then 448 do ir=1,nfft 449 vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1) 450 vtrial1(ir,2)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,2) 451 end do 452 else 453 ! fab: I think there was an error in the definition of vtrial1(2*ir-1,2); I have corrected it... 454 do ir=1,nfft 455 vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1) 456 vtrial1(2*ir ,1)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,1) 457 vtrial1(2*ir-1,2)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1 ,2) 458 vtrial1(2*ir ,2)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,2) 459 end do 460 end if 461 end if 462 463 ! Compute the third-order xc energy 464 ! take into account the contribution of the term 465 !$ 466 ! \frac{d}{d \lambda} 467 ! \frac{\delta^2 E_{Hxc}}{\delta n(r) \delta n(r\prim)} 468 !$ 469 ! (seventh term of Eq. (110) of X. Gonze, PRA 52, 1096 (1995)). 470 471 ! the following are essentially the 4th and the 3rd terms of PRB 71,125107, but the 472 ! multiplication for rho1 will be done by dotprod_vn later 473 474 ! in the non spin polarized case xc_tmp has only 1 component 475 if (nspden==1)then 476 477 ABI_ALLOCATE(xc_tmp,(cplex*nfft,1)) 478 479 if (cplex==1) then 480 ! This, and the next lines, have to be changed in case cplex=2 481 do ifft=1,nfft 482 xc_tmp(ifft,1)= k3xc(ifft,1)*(rho2r1(ifft,1)+3*xccc3d2(ifft))*rho3r1(ifft,1) 483 end do 484 else 485 do ifft=1,nfft ! 2*ifft-1 denotes the real part, 2*ifft the imaginary part 486 xc_tmp(2*ifft-1,1)= k3xc(ifft,1)*( (rho2r1(2*ifft-1,1)+3*xccc3d2(2*ifft-1))*rho3r1(2*ifft-1,1) & 487 & -( rho2r1(2*ifft,1)+3*xccc3d2(2*ifft))*rho3r1(2*ifft,1)) 488 489 xc_tmp(2*ifft,1)= k3xc(ifft,1)*( (rho2r1(2*ifft-1,1)+3*xccc3d2(2*ifft-1))*rho3r1(2*ifft,1) & 490 & +( rho2r1(2*ifft,1)+3*xccc3d2(2*ifft))*rho3r1(2*ifft-1,1)) 491 end do 492 493 end if 494 495 end if 496 497 ! fab: modifications for the spin polarized raman part: 498 ! in the spin polarized case xc_tmp has 2 components 499 ! note that now the non linear core correction is divided by 2 500 if (nspden==2) then 501 502 ABI_ALLOCATE(xc_tmp,(cplex*nfft,2)) 503 504 if (cplex==1) then 505 do ifft=1,nfft 506 xc_tmp(ifft,1)= k3xc(ifft,1)*(rho2r1(ifft,2)+(3._dp/2._dp)*xccc3d2(ifft))*rho3r1(ifft,2)+ & 507 & k3xc(ifft,2)*(rho2r1(ifft,2)+(3._dp/2._dp)*xccc3d2(ifft))*(rho3r1(ifft,1)-rho3r1(ifft,2))+ & 508 & k3xc(ifft,2)*((rho2r1(ifft,1)-rho2r1(ifft,2))+(3._dp/2._dp)*xccc3d2(ifft))*rho3r1(ifft,2)+ & 509 & k3xc(ifft,3)*((rho2r1(ifft,1)-rho2r1(ifft,2))+(3._dp/2._dp)*xccc3d2(ifft))*(rho3r1(ifft,1)-rho3r1(ifft,2)) 510 xc_tmp(ifft,2)= k3xc(ifft,2)*(rho2r1(ifft,2)+(3._dp/2._dp)*xccc3d2(ifft))*rho3r1(ifft,2)+ & 511 & k3xc(ifft,3)*(rho2r1(ifft,2)+(3._dp/2._dp)*xccc3d2(ifft))*(rho3r1(ifft,1)-rho3r1(ifft,2))+ & 512 & k3xc(ifft,3)*((rho2r1(ifft,1)-rho2r1(ifft,2))+(3._dp/2._dp)*xccc3d2(ifft))*rho3r1(ifft,2)+ & 513 & k3xc(ifft,4)*((rho2r1(ifft,1)-rho2r1(ifft,2))+(3._dp/2._dp)*xccc3d2(ifft))*(rho3r1(ifft,1)-rho3r1(ifft,2)) 514 end do 515 516 else 517 do ifft=1,nfft 518 ! These sections should be rewritten, to be easier to read ... (defining intermediate scalars) 519 xc_tmp(2*ifft-1,1)= k3xc(ifft,1)*& 520 & ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft-1,2)- & 521 & (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft,2))+ & 522 & k3xc(ifft,2)*& 523 & ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*(rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))- & 524 & (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*(rho3r1(2*ifft,1)-rho3r1(2*ifft,2)))+ & 525 & k3xc(ifft,2)*& 526 & ( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft-1,2)- & 527 & ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft,2))+ & 528 & k3xc(ifft,3)*& 529 & ( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 530 & (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))- & 531 & ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*& 532 & (rho3r1(2*ifft,1)-rho3r1(2*ifft,2))) 533 xc_tmp(2*ifft,1)=k3xc(ifft,1)*& 534 & ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft,2)+ & 535 & (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft-1,2))+ & 536 & k3xc(ifft,2)*& 537 & ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*(rho3r1(2*ifft,1)-rho3r1(2*ifft,2))+ & 538 & (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*(rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2)))+ & 539 & k3xc(ifft,2)*& 540 & ( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft,2)+ & 541 & ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft-1,2))+ & 542 & k3xc(ifft,3)*& 543 & ( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 544 & (rho3r1(2*ifft,1)-rho3r1(2*ifft,2))+ & 545 & ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*& 546 & (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))) 547 ! fab: now the spin down component 548 xc_tmp(2*ifft-1,2)= k3xc(ifft,2)*& 549 & ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft-1,2)- & 550 & (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft,2))+ & 551 & k3xc(ifft,3)*( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 552 & (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))- & 553 & (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*(rho3r1(2*ifft,1)-rho3r1(2*ifft,2)))+ & 554 & k3xc(ifft,3)*( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 555 & rho3r1(2*ifft-1,2)- & 556 & ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft,2))+ & 557 & k3xc(ifft,4)*( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 558 & (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))- & 559 ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*& 560 & (rho3r1(2*ifft,1)-rho3r1(2*ifft,2))) 561 xc_tmp(2*ifft,2)=k3xc(ifft,1)*( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 562 & rho3r1(2*ifft,2)+ & 563 & (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft-1,2))+ & 564 & k3xc(ifft,3)*( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 565 & (rho3r1(2*ifft,1)-rho3r1(2*ifft,2))+ & 566 & (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*(rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2)))+ & 567 & k3xc(ifft,3)*( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 568 & rho3r1(2*ifft,2)+ & 569 & ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft-1,2))+ & 570 & k3xc(ifft,4)*( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*& 571 & (rho3r1(2*ifft,1)-rho3r1(2*ifft,2))+ & 572 & ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*& 573 & (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))) 574 end do 575 576 ! fab: this is the end if over cplex 577 end if 578 ! fab: this is the enf if over nspden 579 end if 580 581 call dotprod_vn(1,rho1r1,exc3,valuei,nfft,nfftot,nspden,1,xc_tmp,ucvol,mpi_comm_sphgrid=mpi_enreg%comm_fft) 582 ABI_DEALLOCATE(xc_tmp) 583 584 ! Perform DFPT part of the 3dte calculation 585 586 call timab(512,1,tsec) 587 call status(counter,dtfil%filstat,iexit,level,'call dfptnl_resp ') 588 call dfptnl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,& 589 & kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,& 590 & nspinor,nsppol,npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm) 591 call timab(512,2,tsec) 592 593 call status(counter,dtfil%filstat,iexit,level,'after dfptnl_resp') 594 595 ! Describe the perturbation and write out the result 596 if (mpi_enreg%me == 0) then 597 if (i2pert < natom + 1) then 598 write(message,'(a,i3,a,i3)') & 599 & ' j2 : displacement of atom ',i2pert,& 600 & ' along direction ', i2dir 601 end if 602 if (i2pert == dtset%natom + 2) then 603 write(message,'(a,i4)') & 604 & ' j2 : homogeneous electric field along direction ',& 605 & i2dir 606 end if 607 call wrtout(std_out,message,'COLL') 608 call wrtout(ab_out,message,'COLL') 609 write(ab_out,'(20x,a,13x,a)')'real part','imaginary part' 610 write(ab_out,'(5x,a2,1x,f22.10,3x,f22.10)')'xc',exc3*sixth,zero 611 if (i2pert == natom + 2) then 612 write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',& 613 & d3_berry(1,i2dir),d3_berry(2,i2dir) 614 end if 615 write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'dft',& 616 & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),& 617 & d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 618 write(ab_out,*) 619 write(std_out,'(18x,a,11x,a)')'real part','imaginary part' 620 write(std_out,'(5x,a2,1x,f20.10,3x,f20.10)')'xc',exc3*sixth,zero 621 write(std_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',& 622 & d3_berry(1,i2dir),d3_berry(2,i2dir) 623 write(std_out,'(5x,a3,f22.10,3x,f22.10)')'dft',& 624 & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),& 625 & d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 626 write(std_out,*) 627 end if ! mpi_enreg%me == 0 628 629 d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = & 630 & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + exc3*sixth 631 d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = & 632 & d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + d3_berry(:,i2dir) 633 634 end if !rfpert 635 end do !i2dir 636 end do ! i2pert 637 638 end if ! rfpert 639 end do ! i3dir 640 end do ! i3pert 641 642 end if ! rfpert 643 end do ! i1dir 644 end do ! i1pert 645 646 call status(0,dtfil%filstat,iexit,level,'exit ') 647 648 ABI_DEALLOCATE(cg1) 649 ABI_DEALLOCATE(cg3) 650 ABI_DEALLOCATE(eigen1) 651 ABI_DEALLOCATE(rho1r1) 652 ABI_DEALLOCATE(rho2r1) 653 ABI_DEALLOCATE(rho2g1) 654 ABI_DEALLOCATE(rho3r1) 655 ABI_DEALLOCATE(atindx1) 656 ABI_DEALLOCATE(atindx) 657 ABI_DEALLOCATE(nattyp) 658 ABI_DEALLOCATE(ph1d) 659 ABI_DEALLOCATE(ylm) 660 ABI_DEALLOCATE(vtrial1) 661 ABI_DEALLOCATE(vxc1) 662 ABI_DEALLOCATE(vhartr1) 663 ABI_DEALLOCATE(vpsp1) 664 ABI_DEALLOCATE(xccc3d1) 665 ABI_DEALLOCATE(xccc3d2) 666 ABI_DEALLOCATE(xccc3d3) 667 668 call timab(502,2,tsec) 669 670 end subroutine dfptnl_loop