TABLE OF CONTENTS
ABINIT/m_pead_nl_loop [ Modules ]
NAME
m_pead_nl_loop
FUNCTION
COPYRIGHT
Copyright (C) 2002-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_pead_nl_loop 23 24 use defs_basis 25 use defs_wvltypes 26 use m_wffile 27 use m_abicore 28 use m_xmpi 29 use m_hdr 30 use m_dtset 31 use m_dtfil 32 #if defined HAVE_MPI2 33 use mpi 34 #endif 35 36 use defs_datatypes, only : pseudopotential_type 37 use defs_abitypes, only : MPI_type 38 use m_time, only : timab 39 use m_kg, only : getph, mkkpg 40 use m_cgtools, only : dotprod_vn, dotprod_g 41 use m_fft, only : fourdp, fftpac, fourwf 42 use m_ioarr, only : read_rhor 43 use m_pawtab, only : pawtab_type 44 use m_pawrhoij, only : pawrhoij_type 45 use m_pawcprj, only : pawcprj_type 46 use m_inwffil, only : inwffil 47 use m_spacepar, only : hartre 48 use m_initylmg, only : initylmg 49 use m_dfpt_mkvxc, only : dfpt_mkvxc 50 use m_mkcore, only : dfpt_mkcore 51 use m_mklocl, only : dfpt_vlocal 52 use m_hamiltonian,only : init_hamiltonian, gs_hamiltonian_type 53 use m_mkffnl, only : mkffnl 54 use m_mpinfo, only : proc_distrb_cycle 55 use m_nonlop, only : nonlop 56 use m_dfptnl_pert, only : dfptnl_exc3 57 58 implicit none 59 60 private 61 62 #if defined HAVE_MPI1 63 include 'mpif.h' 64 #endif
ABINIT/pead_nl_loop [ Functions ]
NAME
pead_nl_loop
FUNCTION
Loop over the perturbations j1, j2 and j3
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 rhohxc.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
SOURCE
142 subroutine pead_nl_loop(blkflg,cg,cgindex,dtfil,dtset,d3lo,& 143 & gmet,gprimd,gsqcut, & 144 & hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,mband,mgfft,mkmem,mkmem_max,mk1mem,& 145 & mpert,mpi_enreg,mpw,mvwtk,natom,nfft,nkpt,nkpt3,nkxc,nk3xc,nneigh,nspinor,nsppol,& 146 & npwarr,occ,psps,pwind,& 147 & rfpert,rprimd,ucvol,xred) 148 149 !Arguments ------------------------------------ 150 !scalars 151 integer,intent(in) :: mband,mgfft,mk1mem,mkmem,mkmem_max,mpert,mpw,natom,nfft 152 integer,intent(in) :: nk3xc,nkpt,nkpt3,nkxc,nneigh,nspinor,nsppol 153 real(dp),intent(in) :: gsqcut,ucvol 154 type(MPI_type),intent(inout) :: mpi_enreg 155 type(datafiles_type),intent(in) :: dtfil 156 type(dataset_type),intent(inout) :: dtset 157 type(hdr_type),intent(inout) :: hdr 158 type(pseudopotential_type),intent(in) :: psps 159 !arrays 160 integer,intent(in) :: cgindex(nkpt,nsppol),kg(3,mk1mem*mpw),kneigh(30,nkpt) 161 integer,intent(in) :: kg_neigh(30,nkpt,3) 162 integer,intent(in) :: kptindex(2,nkpt3),npwarr(nkpt),pwind(mpw,nneigh,mkmem) 163 integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert) 164 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) !vz_i 165 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3) 166 real(dp),intent(in) :: gprimd(3,3),k3xc(nfft,nk3xc),kpt3(3,nkpt3) 167 real(dp),intent(in) :: kxc(nfft,nkxc),mvwtk(30,nkpt),rprimd(3,3) 168 real(dp),intent(in) :: xred(3,natom) 169 real(dp),intent(inout) :: occ(mband*nkpt*nsppol) 170 real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert) !vz_i 171 172 !Local variables------------------------------- 173 !scalars 174 integer,parameter :: level=51 175 integer :: ask_accurate,counter,cplex,formeig,i1dir 176 integer :: i1pert,i2dir,i2pert,i3dir,i3pert,iatom,ierr,index,ir 177 integer :: ireadwf,itypat,mcg,mpsang,n1,n2,n3,n3xccc,nfftot,nspden,option,optorth 178 integer :: pert1case,pert2case,pert3case,rdwrpaw,timrev,comm_cell 179 logical :: nmxc 180 real(dp) :: ecut_eff,exc3(2) 181 character(len=500) :: message 182 character(len=fnlen) :: fiden1i,fiwf1i,fiwf3i 183 type(wffile_type) :: wff1,wff2,wfft1,wfft2 184 type(wvl_data) :: wvl 185 type(hdr_type) :: hdr_den 186 !arrays 187 integer,allocatable :: atindx(:),atindx1(:),nattyp(:) 188 real(dp) :: d3_berry(2,3),rho_dum(1),tsec(2),ylmgr_dum(1) 189 real(dp),allocatable :: cg1(:,:),cg3(:,:),eigen1(:),ph1d(:,:),rho1r1(:,:) 190 real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:),rho3r1(:,:),vhartr1(:) 191 real(dp),allocatable :: vpsp1(:),vtrial1(:,:),vxc1(:,:),work(:) 192 real(dp),allocatable :: xccc3d1(:),xccc3d2(:),xccc3d3(:),ylm(:,:,:) 193 type(pawrhoij_type),allocatable :: rhoij_dum(:) 194 195 ! *********************************************************************** 196 197 call timab(502,1,tsec) 198 199 comm_cell = mpi_enreg%comm_cell 200 201 timrev = 1 202 cplex = 2 - timrev 203 nspden = dtset%nspden 204 ecut_eff = (dtset%ecut)*(dtset%dilatmx)**2 205 mpsang = psps%mpsang 206 optorth=1;if (psps%usepaw==1) optorth=0 207 208 ABI_MALLOC(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 209 ABI_MALLOC(cg3,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 210 ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol)) 211 ABI_MALLOC(rho1r1,(cplex*nfft,dtset%nspden)) 212 ABI_MALLOC(rho2r1,(cplex*nfft,dtset%nspden)) 213 ABI_MALLOC(rho2g1,(2,nfft)) 214 ABI_MALLOC(rho3r1,(cplex*nfft,dtset%nspden)) 215 ABI_MALLOC(ylm,(2,dtset%mpw*dtset%mkmem,mpsang*mpsang*psps%useylm)) 216 217 ask_accurate=1 ; formeig = 1 ; ireadwf = 1 218 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 219 nfftot=n1*n2*n3 220 221 !Generate an index table of atoms, in order for them to be used 222 !type after type. 223 ABI_MALLOC(atindx,(natom)) 224 ABI_MALLOC(atindx1,(natom)) 225 ABI_MALLOC(nattyp,(psps%ntypat)) 226 index=1 227 do itypat=1,psps%ntypat 228 nattyp(itypat)=0 229 do iatom=1,natom 230 if(dtset%typat(iatom)==itypat)then 231 atindx(iatom)=index 232 atindx1(index)=iatom 233 index=index+1 234 nattyp(itypat)=nattyp(itypat)+1 235 end if 236 end do 237 end do 238 239 !Generate the 1-dimensional phases 240 ABI_MALLOC(ph1d,(2,3*(2*mgfft+1)*natom)) 241 call getph(atindx,natom,n1,n2,n3,ph1d,xred) 242 243 !Set up the Ylm for each k point 244 if (psps%useylm==1) then 245 call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,& 246 & dtset%mpw,dtset%nband,dtset%nkpt,& 247 & npwarr,dtset%nsppol,0,rprimd,ylm,ylmgr_dum) 248 end if 249 250 ABI_MALLOC(vpsp1,(cplex*nfft)) 251 ABI_MALLOC(xccc3d1,(cplex*nfft)) 252 ABI_MALLOC(xccc3d2,(cplex*nfft)) 253 ABI_MALLOC(xccc3d3,(cplex*nfft)) 254 ABI_MALLOC(vhartr1,(cplex*nfft)) 255 ABI_MALLOC(vxc1,(cplex*nfft,dtset%nspden)) 256 ABI_MALLOC(vtrial1,(cplex*nfft,dtset%nspden)) 257 258 !Loop over the perturbations j1, j2, j3 259 260 pert1case = 0 ; pert2case = 0 ; pert3case = 0 261 262 do i1pert = 1, mpert 263 do i1dir = 1, 3 264 265 if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then 266 267 pert1case = i1dir + (i1pert-1)*3 268 counter = pert1case 269 call appdig(pert1case,dtfil%fnamewff1,fiwf1i) 270 271 mcg=mpw*nspinor*mband*mkmem*nsppol 272 call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,& 273 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 274 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 275 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 276 & dtset%nsppol,dtset%nsym,& 277 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 278 & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl) 279 280 if (ireadwf==1) then 281 call WffClose (wff1,ierr) 282 end if 283 284 rho1r1(:,:) = 0._dp 285 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 286 rdwrpaw=0 287 call appdig(pert1case,dtfil%fildens1in,fiden1i) 288 289 call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho1r1, & 290 hdr_den, rhoij_dum, comm_cell, check_hdr=hdr) 291 call hdr_den%free() 292 end if 293 294 xccc3d1(:) = 0._dp 295 if ((psps%n1xccc/=0).and.(i1pert <= natom)) then 296 call dfpt_mkcore(cplex,i1dir,i1pert,natom,psps%ntypat,n1,psps%n1xccc,& 297 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,& 298 & psps%xcccrc,psps%xccc1d,xccc3d1,xred) 299 end if ! psps%n1xccc/=0 300 301 do i3pert = 1, mpert 302 do i3dir = 1, 3 303 304 if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then 305 306 pert3case = i3dir + (i3pert-1)*3 307 counter = 100*pert3case + pert1case 308 call appdig(pert3case,dtfil%fnamewff1,fiwf3i) 309 310 mcg=mpw*nspinor*mband*mkmem*nsppol 311 call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,& 312 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 313 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 314 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 315 & dtset%nsppol,dtset%nsym,& 316 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 317 & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,& 318 & fiwf3i,wvl) 319 if (ireadwf==1) then 320 call WffClose (wff2,ierr) 321 end if 322 323 rho3r1(:,:) = 0._dp 324 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 325 rdwrpaw=0 326 call appdig(pert3case,dtfil%fildens1in,fiden1i) 327 328 call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho3r1, & 329 hdr_den, rhoij_dum, comm_cell, check_hdr=hdr) 330 call hdr_den%free() 331 end if 332 333 xccc3d3(:) = 0._dp 334 if ((psps%n1xccc/=0).and.(i3pert <= natom)) then 335 call dfpt_mkcore(cplex,i3dir,i3pert,natom,psps%ntypat,n1,psps%n1xccc,& 336 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,& 337 & psps%xcccrc,psps%xccc1d,xccc3d3,xred) 338 end if ! psps%n1xccc/=0 339 340 do i2pert = 1, mpert 341 342 ! In case of electric field perturbation, evaluate the ddk 343 ! using the finite difference expression of 344 ! Marzari and Vanderbilt PRB 56, 12847 (1997) [[cite:Marzari1997]]. 345 346 d3_berry(:,:) = 0._dp 347 348 if ((i2pert==dtset%natom+2).and.& 349 & (maxval(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert)) == 1)) then 350 351 call timab(511,1,tsec) 352 call pead_nl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,& 353 & i1pert,i3pert,i1dir,i3dir,& 354 & kneigh,kg_neigh,kptindex,kpt3,mband,mkmem,mkmem_max,mk1mem,& 355 & mpi_enreg,mpw,mvwtk,natom,nkpt,nkpt3,nneigh,npwarr,nspinor,nsppol,pwind) 356 call timab(511,2,tsec) 357 358 end if 359 360 if (mpi_enreg%me == 0) then 361 362 if(sum(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert))>0)then 363 write(message,'(a,a,a,a,a,a)')ch10,ch10,& 364 & ' Decomposition of the third-order energy for the set of perturbations',ch10 365 call wrtout(std_out,message,'COLL') 366 call wrtout(ab_out,message,'COLL') 367 if (i1pert < natom + 1) then 368 write(message,'(a,i3,a,i3)') & 369 & ' j1 : displacement of atom ',i1pert,' along direction ', i1dir 370 end if 371 if (i1pert == dtset%natom + 2) then 372 write(message,'(a,i4)')' j1 : homogeneous electric field along direction ',i1dir 373 end if 374 call wrtout(std_out,message,'COLL') 375 call wrtout(ab_out,message,'COLL') 376 if (i3pert < natom + 1) then 377 write(message,'(a,i3,a,i3,a)') & 378 & ' j3 : displacement of atom ',i3pert,' along direction ', i3dir,ch10 379 end if 380 if (i3pert == dtset%natom + 2) then 381 write(message,'(a,i4,a)')' j3 : homogeneous electric field along direction ',i3dir,ch10 382 end if 383 call wrtout(std_out,message,'COLL') 384 call wrtout(ab_out,message,'COLL') 385 end if 386 387 end if ! mpi_enreg%me == 0 388 389 do i2dir = 1, 3 390 391 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 392 pert2case = i2dir + (i2pert-1)*3 393 394 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 395 396 ! Read the first-order densities from disk-files 397 rho2r1(:,:) = 0._dp ; rho2g1(:,:) = 0._dp 398 399 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 400 rdwrpaw=0 401 call appdig(pert2case,dtfil%fildens1in,fiden1i) 402 403 call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho2r1, & 404 hdr_den, rhoij_dum, comm_cell, check_hdr=hdr) 405 call hdr_den%free() 406 407 ! Compute up+down rho1(G) by fft 408 ABI_MALLOC(work,(cplex*nfft)) 409 work(:)=rho2r1(:,1) 410 call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfft,1,dtset%ngfft,0) 411 ABI_FREE(work) 412 413 end if 414 415 ! Compute first-order local potentials 416 ! (hartree, xc and pseudopotential) 417 418 n3xccc=0; if(psps%n1xccc/=0)n3xccc=nfft 419 xccc3d2(:)=0._dp ; vpsp1(:)=0._dp 420 421 if (i2pert <= natom) then 422 423 call dfpt_vlocal(atindx,cplex,gmet,gsqcut,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,natom,& 424 & nattyp,nfft,dtset%ngfft,psps%ntypat,n1,n2,n3,ph1d,psps%qgrid_vl,& 425 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 426 427 if (psps%n1xccc/=0) then 428 call dfpt_mkcore(cplex,i2dir,i2pert,natom,psps%ntypat,n1,psps%n1xccc,& 429 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,& 430 & psps%xcccrc,psps%xccc1d,xccc3d2,xred) 431 end if ! psps%n1xccc/=0 432 433 end if ! i2pert <= natom 434 435 call hartre(cplex,gsqcut,3,0,mpi_enreg,nfft,dtset%ngfft,dtset%nkpt,& 436 &dtset%rcut,rho2g1,rprimd,dtset%vcutgeo,vhartr1) 437 option=1 ; nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 438 call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,dtset%ngfft,& 439 & rho_dum,0,rho_dum,0,nkxc,nmxc,dtset%nspden,n3xccc,option,& 440 & dtset%qptn,rho2r1,rprimd,0,vxc1,xccc3d2) 441 442 if(dtset%nsppol==1)then 443 if(cplex==1)then 444 do ir=1,nfft 445 vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1) 446 end do 447 else 448 do ir=1,nfft 449 vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1) 450 vtrial1(2*ir ,1)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,1) 451 end do 452 end if 453 else 454 if(cplex==1)then 455 do ir=1,nfft 456 vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1) 457 vtrial1(ir,2)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,2) 458 end do 459 else 460 ! fab: I think there was an error in the definition of vtrial1(2*ir-1,2); I have corrected it... 461 do ir=1,nfft 462 vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1) 463 vtrial1(2*ir ,1)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,1) 464 vtrial1(2*ir-1,2)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1 ,2) 465 vtrial1(2*ir ,2)=vpsp1(2*ir )+vhartr1(2*ir )+vxc1(2*ir ,2) 466 end do 467 end if 468 end if 469 470 ! Compute the third-order xc energy 471 call dfptnl_exc3(cplex,exc3,k3xc,mpi_enreg,nk3xc,nfft,nfftot,dtset%nspden,& 472 & rho1r1,rho2r1,rho3r1,ucvol,xccc3d1,xccc3d2,xccc3d3) 473 474 ! Perform DFPT part of the 3dte calculation 475 476 call timab(512,1,tsec) 477 call pead_nl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,& 478 & kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,& 479 & nspinor,nsppol,npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm) 480 call timab(512,2,tsec) 481 482 483 ! Describe the perturbation and write out the result 484 if (mpi_enreg%me == 0) then 485 if (i2pert < natom + 1) then 486 write(message,'(a,i3,a,i3)') & 487 & ' j2 : displacement of atom ',i2pert,& 488 & ' along direction ', i2dir 489 end if 490 if (i2pert == dtset%natom + 2) then 491 write(message,'(a,i4)') & 492 & ' j2 : homogeneous electric field along direction ',& 493 & i2dir 494 end if 495 call wrtout(std_out,message,'COLL') 496 call wrtout(ab_out,message,'COLL') 497 write(ab_out,'(20x,a,13x,a)')'real part','imaginary part' 498 write(ab_out,'(5x,a2,1x,f22.10,3x,f22.10)')'xc',exc3(1)*sixth,zero 499 if (i2pert == natom + 2) then 500 write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',& 501 & d3_berry(1,i2dir),d3_berry(2,i2dir) 502 end if 503 write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'dft',& 504 & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),& 505 & d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 506 write(ab_out,*) 507 write(std_out,'(18x,a,11x,a)')'real part','imaginary part' 508 write(std_out,'(5x,a2,1x,f20.10,3x,f20.10)')'xc',exc3(1)*sixth,zero 509 write(std_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',& 510 & d3_berry(1,i2dir),d3_berry(2,i2dir) 511 write(std_out,'(5x,a3,f22.10,3x,f22.10)')'dft',& 512 & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),& 513 & d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 514 write(std_out,*) 515 end if ! mpi_enreg%me == 0 516 517 d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = & 518 & d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + exc3(1)*sixth 519 d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = & 520 & d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + d3_berry(:,i2dir) 521 522 end if !rfpert 523 end do !i2dir 524 end do ! i2pert 525 526 end if ! rfpert 527 end do ! i3dir 528 end do ! i3pert 529 530 end if ! rfpert 531 end do ! i1dir 532 end do ! i1pert 533 534 535 ABI_FREE(cg1) 536 ABI_FREE(cg3) 537 ABI_FREE(eigen1) 538 ABI_FREE(rho1r1) 539 ABI_FREE(rho2r1) 540 ABI_FREE(rho2g1) 541 ABI_FREE(rho3r1) 542 ABI_FREE(atindx1) 543 ABI_FREE(atindx) 544 ABI_FREE(nattyp) 545 ABI_FREE(ph1d) 546 ABI_FREE(ylm) 547 ABI_FREE(vtrial1) 548 ABI_FREE(vxc1) 549 ABI_FREE(vhartr1) 550 ABI_FREE(vpsp1) 551 ABI_FREE(xccc3d1) 552 ABI_FREE(xccc3d2) 553 ABI_FREE(xccc3d3) 554 555 call timab(502,2,tsec) 556 557 end subroutine pead_nl_loop
ABINIT/pead_nl_mv [ Functions ]
NAME
pead_nl_mv
FUNCTION
Compute the finite difference expression of the k-point derivative using the PEAD formulation of the third-order energy (see Nunes and Gonze PRB 63, 155107 (2001) [[cite:Nunes2001]] Eq. 102) and the finite difference formula of Marzari and Vanderbilt (see Marzari and Vanderbilt, PRB 56, 12847 (1997) [[cite:Marzari1997]], Appendix B)
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions cgindex(nkpt2,nsppol) = for each k-point, cgindex stores the location of the WF in the cg array cg1 = first-order wavefunction relative to the perturbations i1pert cg3 = first-order wavefunction relative to the perturbations i3pert 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 i1pert,i3pert = type of perturbation that has to be computed i1dir,i3dir=directions of the corresponding perturbations kneigh(30,nkpt2) = index of the neighbours of each k-point kg_neigh(30,nkpt2,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. See getshell.F90 for details 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 mband = maximum number of bands mkmem = maximum number of k points which can fit in core memory mkmem_max = maximal number of k-points on each processor (MPI //) mk1mem = maximum number of k points for first-order WF which can fit in core memory 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 nkpt2 = number of k-points in the reduced part of the BZ nkpt2 = nkpt/2 in case of time-reversal symmetry (kptopt = 2) nkpt3 = number of k-points in the full BZ nneigh = total number of neighbours required to evaluate the finite difference formula npwarr(nkpt) = array holding npw for each k point nspinor = number of spinorial components of the wavefunctions nsppol = number of channels for spin-polarization (1 or 2) pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat between k-points
OUTPUT
d3_berry(2,3) = Berry-phase part of the third-order energy
SIDE EFFECTS
mpi_enreg=MPI-parallelisation information
NOTES
For a given set of values of i1pert,i3pert,i1dir and i3dir, the routine computes the k-point derivatives for 12dir = 1,2,3
SOURCE
920 subroutine pead_nl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,& 921 & i1pert,i3pert,i1dir,i3dir,kneigh,kg_neigh,kptindex,& 922 & kpt3,mband,mkmem,mkmem_max,mk1mem,mpi_enreg,& 923 & mpw,mvwtk,natom,nkpt2,nkpt3,nneigh,npwarr,nspinor,& 924 & nsppol,pwind) 925 926 use m_hide_lapack, only : dzgedi, dzgefa 927 928 !Arguments ------------------------------------ 929 ! 930 !--- Arguments : integer scalars 931 integer, intent(in) :: i1dir,i1pert,i3dir,i3pert,mband,mk1mem 932 integer, intent(in) :: mkmem,mkmem_max,mpw,natom 933 integer, intent(in) :: nkpt2,nkpt3,nneigh,nspinor,nsppol 934 ! 935 !--- Arguments : integer arrays 936 integer, intent(in) :: cgindex(nkpt2,nsppol) 937 integer, intent(in) :: kneigh(30,nkpt2),kg_neigh(30,nkpt2,3),kptindex(2,nkpt3) 938 integer, intent(in) :: npwarr(nkpt2),pwind(mpw,nneigh,mkmem) 939 ! 940 !--- Arguments : real(dp) scalars 941 ! 942 !--- Arguments : real(dp) arrays 943 real(dp), intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 944 real(dp), intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol) 945 real(dp), intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol) 946 real(dp), intent(in) :: gmet(3,3),kpt3(3,nkpt3) 947 real(dp), intent(in) :: mvwtk(30,nkpt2) 948 real(dp), intent(out) :: d3_berry(2,3) 949 ! 950 !--- Arguments : structured datatypes 951 type(MPI_type), intent(in) :: mpi_enreg 952 type(datafiles_type), intent(in) :: dtfil 953 type(dataset_type), intent(in) :: dtset 954 955 !Local variables------------------------------- 956 ! 957 !---- Local variables : integer scalars 958 integer :: count,counter,count1,iband,icg 959 integer :: ierr,ii,ikpt,ikpt_loc,ikpt2 960 integer :: ikpt_rbz,ineigh,info,ipw,isppol,jband,jcg,jj,jkpt,job,jpw, jkpt2, jkpt_rbz 961 integer :: lband,lpband,nband_occ,npw_k,npw_k1,my_source,his_source,dest,tag 962 integer :: spaceComm 963 integer,parameter :: level=52 964 integer :: bdtot_index 965 ! 966 !---- Local variables : integer arrays 967 integer,allocatable :: ipvt(:) 968 integer, allocatable :: bd_index(:,:) 969 ! 970 !---- Local variables : real(dp) scalars 971 real(dp) :: dotnegi,dotnegr,dotposi,dotposr 972 ! real(dp) :: c1,c2 ! appear commented out below 973 ! 974 !---- Local variables : real(dp) arrays 975 real(dp) :: d3_aux(2,3),det(2,2),dk(3),dk_(3) 976 real(dp) :: z1(2),z2(2) 977 real(dp),allocatable :: buffer(:,:),cgq(:,:),cg1q(:,:),cg3q(:,:) 978 real(dp),allocatable :: qmat(:,:,:),s13mat(:,:,:),s1mat(:,:,:),s3mat(:,:,:) 979 real(dp),allocatable :: smat(:,:,:),zgwork(:,:) 980 ! 981 !---- Local variables : character variables 982 character(len=500) :: message 983 ! 984 !---- Local variables : structured datatypes 985 986 987 #if defined HAVE_MPI 988 integer :: status1(MPI_STATUS_SIZE) 989 spaceComm=mpi_enreg%comm_cell 990 #endif 991 992 ABI_UNUSED(dtfil%ireadwf) 993 994 ! *********************************************************************** 995 996 write(message,'(8a)') ch10,& 997 & ' pead_nl_mv : finite difference expression of the k-point derivative',ch10,& 998 & ' is performed using the PEAD formulation of ',& 999 & 'the third-order energy',ch10,& 1000 & ' (see Nunes and Gonze PRB 63, 155107 (2001) [[cite:Nunes2001]] Eq. 102)',ch10 1001 !call wrtout(ab_out,message,'COLL') 1002 call wrtout(std_out, message,'COLL') 1003 1004 1005 !fab: I think that the following restriction must be eliminated: 1006 !isppol = 1 1007 1008 ikpt_loc = 0 1009 d3_aux(:,:) = 0_dp 1010 1011 ABI_MALLOC(s13mat,(2,mband,mband)) 1012 ABI_MALLOC(smat,(2,mband,mband)) 1013 ABI_MALLOC(s1mat,(2,mband,mband)) 1014 ABI_MALLOC(qmat,(2,mband,mband)) 1015 ABI_MALLOC(ipvt,(mband)) 1016 ABI_MALLOC(s3mat,(2,mband,mband)) 1017 ABI_MALLOC(zgwork,(2,mband)) 1018 ABI_MALLOC(bd_index, (nkpt2, nsppol)) 1019 1020 bdtot_index = 0 1021 do isppol = 1, nsppol 1022 do ikpt_rbz = 1, nkpt2 1023 bd_index(ikpt_rbz,isppol) = bdtot_index 1024 bdtot_index = bdtot_index + dtset%nband(ikpt_rbz+nkpt2*(isppol-1)) 1025 end do 1026 end do 1027 1028 !fab: I think here I have to add the loop over spin 1029 1030 do isppol = 1, nsppol 1031 1032 ! Loop over k-points 1033 ! COMMENT: Every processor has to make mkmem_max iterations 1034 ! even if mkmem < mkemem_max. This is due to the fact 1035 ! that it still has to communicate its wavefunctions 1036 ! to other processors even if it has no more overlap 1037 ! matrices to compute. 1038 1039 ikpt_loc = 0 ; ikpt = 0 1040 1041 do while (ikpt_loc < mkmem_max) 1042 1043 if (ikpt_loc < mkmem) ikpt = ikpt + 1 1044 1045 if (xmpi_paral == 1) then 1046 ! if ((minval(abs(mpi_enreg%proc_distrb(ikpt,1:mband,1:dtset%nsppol) & 1047 ! & - mpi_enreg%me)) /= 0).and.(ikpt_loc < mkmem)) cycle 1048 if(ikpt>nkpt2)then 1049 ikpt_loc=mkmem_max 1050 cycle 1051 end if 1052 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me)) then 1053 if(ikpt==nkpt2) ikpt_loc=mkmem_max 1054 cycle 1055 end if 1056 end if 1057 1058 ikpt_loc = ikpt_loc + 1 1059 npw_k = npwarr(ikpt) 1060 counter = 100*ikpt 1061 1062 ii = cgindex(ikpt,isppol) 1063 1064 ! Loop on the neighbours 1065 1066 do ineigh = 1,nneigh 1067 1068 s13mat(:,:,:) = zero 1069 smat(:,:,:) = zero 1070 s1mat(:,:,:) = zero 1071 s3mat(:,:,:) = zero 1072 qmat(:,:,:) = zero 1073 1074 ikpt2 = kneigh(ineigh,ikpt) 1075 ikpt_rbz = kptindex(1,ikpt2) ! index of the k-point in the reduced BZ 1076 jj = cgindex(ikpt_rbz,isppol) 1077 ! previous fixed value for nband_k now called nband_occ: 1078 !nband_occ = dtset%nband(ikpt_rbz+nkpt2*(isppol-1)) 1079 ! TODO: check if all these bands are occupied in nsppol = 2 case 1080 nband_occ = 0 1081 do iband = 1, dtset%nband(ikpt_rbz+nkpt2*(isppol-1)) 1082 !Note, only one image is allowed here (or occ_orig should be the same or all images) 1083 if (dtset%occ_orig(bd_index(ikpt_rbz,isppol) + iband,1) > tol10) nband_occ = nband_occ + 1 1084 end do 1085 npw_k1 = npwarr(ikpt_rbz) 1086 dk_(:) = kpt3(:,ikpt2) - dtset%kptns(:,ikpt) 1087 dk(:) = dk_(:) - nint(dk_(:)) + real(kg_neigh(ineigh,ikpt,:),dp) 1088 1089 count = nspinor*mband*npw_k1 1090 ABI_MALLOC(cgq,(2,count)) 1091 ABI_MALLOC(cg1q,(2,count)) 1092 ABI_MALLOC(cg3q,(2,count)) 1093 1094 #if defined HAVE_MPI 1095 1096 my_source = mpi_enreg%proc_distrb(ikpt_rbz,1,1) 1097 1098 ! do dest = 0, mpi_enreg%nproc-1 1099 do dest = 0, maxval(mpi_enreg%proc_distrb(1:nkpt2,1:mband,1:dtset%nsppol)) 1100 1101 if ((dest==mpi_enreg%me).and.(ikpt_loc <= mkmem)) then 1102 ! I am dest and have something to do 1103 1104 if (my_source == mpi_enreg%me) then 1105 ! I am destination and source 1106 jcg = cgindex(ikpt_rbz,isppol) 1107 1108 cgq(:,1:count) = cg(:,jcg+1:jcg+count) 1109 cg1q(:,1:count) = cg1(:,jcg+1:jcg+count) 1110 cg3q(:,1:count) = cg3(:,jcg+1:jcg+count) 1111 1112 else 1113 ! I am the destination but not the source -> receive 1114 1115 tag = ikpt_rbz 1116 1117 ABI_MALLOC(buffer,(2,3*count)) 1118 1119 call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,my_source,tag,spaceComm,status1,ierr) 1120 1121 cgq(:,1:count) = buffer(:,1:count) 1122 cg1q(:,1:count) = buffer(:,count+1:2*count) 1123 cg3q(:,1:count) = buffer(:,2*count+1:3*count) 1124 ABI_FREE(buffer) 1125 1126 end if 1127 1128 else if (ikpt_loc <= mpi_enreg%mkmem(dest)) then ! dest != me and the dest has a k-point to treat 1129 1130 jkpt=mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,1) 1131 jkpt2 = kneigh(ineigh,jkpt) 1132 jkpt_rbz = kptindex(1,jkpt2) ! index of the k-point in the reduced BZ 1133 1134 his_source = mpi_enreg%proc_distrb(jkpt_rbz,1,1) 1135 1136 if (his_source == mpi_enreg%me) then 1137 1138 jcg = cgindex(jkpt_rbz,isppol) 1139 1140 tag = jkpt_rbz 1141 count1 = npwarr(jkpt_rbz)*mband*nspinor 1142 ABI_MALLOC(buffer,(2,3*count1)) 1143 buffer(:,1:count1) = cg(:,jcg+1:jcg+count1) 1144 buffer(:,count1+1:2*count1) = cg1(:,jcg+1:jcg+count1) 1145 buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1) 1146 1147 call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,dest,tag,spaceComm,ierr) 1148 1149 ABI_FREE(buffer) 1150 1151 end if 1152 1153 end if 1154 1155 end do 1156 ! 1157 ! do jkpt = 1, nkpt2 1158 ! 1159 ! if ((jkpt == ikpt_rbz).and.(source /= mpi_enreg%me).and.& 1160 ! & (ikpt_loc <= mkmem)) then 1161 ! 1162 ! tag = jkpt 1163 ! 1164 ! allocate(buffer(2,3*count)) 1165 ! call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,& 1166 ! source,tag,spaceComm,status1,ierr) 1167 ! 1168 ! cgq(:,1:count) = buffer(:,1:count) 1169 ! cg1q(:,1:count) = buffer(:,count+1:2*count) 1170 ! cg3q(:,1:count) = buffer(:,2*count+1:3*count) 1171 ! deallocate(buffer) 1172 ! 1173 ! end if 1174 ! 1175 ! ! ---------------------------------------------------------------------------- 1176 ! ! --------------- Here: send the WF to all the cpus that need it ------------- 1177 ! ! ---------------------------------------------------------------------------- 1178 ! 1179 ! do dest = 1, mpi_enreg%nproc 1180 ! 1181 ! if ((minval(abs(mpi_enreg%proc_distrb(jkpt,1:mband,1:dtset%nsppol) & 1182 ! & - mpi_enreg%me)) == 0).and.& 1183 ! & (mpi_enreg%kptdstrb(dest,ineigh,ikpt_loc) == jkpt)) then 1184 ! 1185 ! 1186 ! 1187 ! jcg = cgindex(jkpt,isppol) 1188 ! 1189 ! if (((dest-1) == mpi_enreg%me)) then 1190 ! 1191 ! cgq(:,1:count) = cg(:,jcg+1:jcg+count) 1192 ! cg1q(:,1:count) = cg1(:,jcg+1:jcg+count) 1193 ! cg3q(:,1:count) = cg3(:,jcg+1:jcg+count) 1194 ! 1195 ! else 1196 ! 1197 ! tag = jkpt 1198 ! count1 = npwarr(jkpt)*mband*nspinor 1199 ! allocate(buffer(2,3*count1)) 1200 ! buffer(:,1:count1) = cg(:,jcg+1:jcg+count1) 1201 ! buffer(:,count1+1:2*count1) = cg1(:,jcg+1:jcg+count1) 1202 ! buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1) 1203 ! 1204 ! call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,(dest-1),tag,spaceComm,ierr) 1205 ! 1206 ! deallocate(buffer) 1207 ! 1208 ! end if 1209 ! 1210 ! end if 1211 ! 1212 ! end do ! loop over dest 1213 ! 1214 ! end do ! loop over jkpt 1215 1216 if (ikpt_loc > mkmem) then 1217 ABI_FREE(cgq) 1218 ABI_FREE(cg1q) 1219 ABI_FREE(cg3q) 1220 cycle 1221 end if 1222 1223 #else 1224 ! no // over k-points 1225 1226 cgq(:,1:count) = cg(:,jj+1:jj+count) 1227 cg1q(:,1:count) = cg1(:,jj+1:jj+count) 1228 cg3q(:,1:count) = cg3(:,jj+1:jj+count) 1229 1230 #endif 1231 1232 ! Compute overlap matrices 1233 1234 if (kptindex(2,ikpt2) == 0) then ! no time-reversal symmetry 1235 1236 do ipw = 1, npw_k 1237 1238 jpw = pwind(ipw,ineigh,ikpt_loc) 1239 if (jpw /= 0) then 1240 1241 do iband = 1, nband_occ 1242 do jband = 1, nband_occ 1243 1244 icg = ii + (iband-1)*npw_k + ipw 1245 jcg = (jband-1)*npw_k1 + jpw 1246 1247 smat(1,iband,jband) = smat(1,iband,jband) + & 1248 & cg(1,icg)*cgq(1,jcg) + cg(2,icg)*cgq(2,jcg) 1249 smat(2,iband,jband) = smat(2,iband,jband) + & 1250 & cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg) 1251 1252 s13mat(1,iband,jband) = s13mat(1,iband,jband) + & 1253 & cg1(1,icg)*cg3q(1,jcg) + cg1(2,icg)*cg3q(2,jcg) 1254 s13mat(2,iband,jband) = s13mat(2,iband,jband) + & 1255 & cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg) 1256 1257 s1mat(1,iband,jband) = s1mat(1,iband,jband) + & 1258 & cg1(1,icg)*cgq(1,jcg) + cg1(2,icg)*cgq(2,jcg) + & 1259 & cg(1,icg)*cg1q(1,jcg) + cg(2,icg)*cg1q(2,jcg) 1260 s1mat(2,iband,jband) = s1mat(2,iband,jband) + & 1261 & cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) + & 1262 & cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg) 1263 1264 s3mat(1,iband,jband) = s3mat(1,iband,jband) + & 1265 & cg3(1,icg)*cgq(1,jcg) + cg3(2,icg)*cgq(2,jcg) + & 1266 & cg(1,icg)*cg3q(1,jcg) + cg(2,icg)*cg3q(2,jcg) 1267 s3mat(2,iband,jband) = s3mat(2,iband,jband) + & 1268 & cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) + & 1269 & cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg) 1270 1271 end do 1272 end do 1273 1274 end if 1275 1276 end do ! ipw 1277 1278 else ! use time-reversal symmetry 1279 1280 do ipw = 1,npw_k 1281 1282 jpw = pwind(ipw,ineigh,ikpt_loc) 1283 if (jpw /= 0) then 1284 1285 do iband = 1, nband_occ 1286 do jband = 1, nband_occ 1287 1288 icg = ii + (iband-1)*npw_k + ipw 1289 jcg = (jband-1)*npw_k1 + jpw 1290 1291 smat(1,iband,jband) = smat(1,iband,jband) + & 1292 & cg(1,icg)*cgq(1,jcg) - cg(2,icg)*cgq(2,jcg) 1293 smat(2,iband,jband) = smat(2,iband,jband) - & 1294 & cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg) 1295 1296 s13mat(1,iband,jband) = s13mat(1,iband,jband) + & 1297 & cg1(1,icg)*cg3q(1,jcg) - cg1(2,icg)*cg3q(2,jcg) 1298 s13mat(2,iband,jband) = s13mat(2,iband,jband) - & 1299 & cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg) 1300 1301 s1mat(1,iband,jband) = s1mat(1,iband,jband) + & 1302 & cg1(1,icg)*cgq(1,jcg) - cg1(2,icg)*cgq(2,jcg) + & 1303 & cg(1,icg)*cg1q(1,jcg) - cg(2,icg)*cg1q(2,jcg) 1304 s1mat(2,iband,jband) = s1mat(2,iband,jband) - & 1305 & cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) - & 1306 & cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg) 1307 1308 s3mat(1,iband,jband) = s3mat(1,iband,jband) + & 1309 & cg3(1,icg)*cgq(1,jcg) - cg3(2,icg)*cgq(2,jcg) + & 1310 & cg(1,icg)*cg3q(1,jcg) - cg(2,icg)*cg3q(2,jcg) 1311 s3mat(2,iband,jband) = s3mat(2,iband,jband) - & 1312 & cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) - & 1313 & cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg) 1314 1315 end do 1316 end do 1317 1318 end if 1319 1320 end do ! ipw 1321 1322 end if 1323 1324 ABI_FREE(cgq) 1325 ABI_FREE(cg1q) 1326 ABI_FREE(cg3q) 1327 1328 ! Compute qmat, the inverse of smat 1329 1330 job = 1 ! compute inverse only 1331 qmat(:,:,:) = smat(:,:,:) 1332 1333 call dzgefa(qmat,mband,nband_occ,ipvt,info) 1334 call dzgedi(qmat,mband,nband_occ,ipvt,det,zgwork,job) 1335 1336 ! DEBUG 1337 ! write(100,*) 1338 ! write(100,*)'ikpt = ',ikpt,'ineigh = ',ineigh 1339 ! do iband = 1,nband_occ 1340 ! do jband = 1,nband_occ 1341 ! c1 = 0_dp ; c2 = 0_dp 1342 ! do lband = 1,nband_occ 1343 ! c1 = c1 + smat(1,iband,lband)*qmat(1,lband,jband) - & 1344 ! & smat(2,iband,lband)*qmat(2,lband,jband) 1345 ! c2 = c2 + smat(1,iband,lband)*qmat(2,lband,jband) + & 1346 ! & smat(2,iband,lband)*qmat(1,lband,jband) 1347 ! end do 1348 ! write(100,'(2(2x,i2),2(2x,f16.9))')iband,jband,& 1349 ! & c1,c2 1350 ! end do 1351 ! end do 1352 ! ENDDEBUG 1353 1354 1355 1356 ! Accumulate sum over bands 1357 1358 dotposr = 0_dp ; dotposi = 0_dp 1359 dotnegr = 0_dp ; dotnegi = 0_dp 1360 do iband = 1, nband_occ 1361 do jband = 1, nband_occ 1362 1363 dotposr = dotposr + & 1364 & s13mat(1,iband,jband)*qmat(1,jband,iband) - & 1365 & s13mat(2,iband,jband)*qmat(2,jband,iband) 1366 dotposi = dotposi + & 1367 & s13mat(1,iband,jband)*qmat(2,jband,iband) + & 1368 & s13mat(2,iband,jband)*qmat(1,jband,iband) 1369 1370 1371 do lband = 1, nband_occ 1372 do lpband= 1, nband_occ 1373 1374 z1(1) = s1mat(1,iband,jband)*qmat(1,jband,lband) - & 1375 & s1mat(2,iband,jband)*qmat(2,jband,lband) 1376 z1(2) = s1mat(1,iband,jband)*qmat(2,jband,lband) + & 1377 & s1mat(2,iband,jband)*qmat(1,jband,lband) 1378 1379 z2(1) = s3mat(1,lband,lpband)*qmat(1,lpband,iband) - & 1380 & s3mat(2,lband,lpband)*qmat(2,lpband,iband) 1381 z2(2) = s3mat(1,lband,lpband)*qmat(2,lpband,iband) + & 1382 & s3mat(2,lband,lpband)*qmat(1,lpband,iband) 1383 1384 dotnegr = dotnegr + & 1385 & z1(1)*z2(1) - z1(2)*z2(2) 1386 dotnegi = dotnegi + & 1387 & z1(1)*z2(2) + z1(2)*z2(1) 1388 1389 end do ! lpband 1390 end do ! lband 1391 1392 end do ! jband 1393 end do ! iband 1394 1395 d3_aux(1,:) = d3_aux(1,:) + & 1396 & dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposr-dotnegr) 1397 d3_aux(2,:) = d3_aux(2,:) + & 1398 & dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposi-dotnegi) 1399 1400 end do ! End loop over neighbours 1401 1402 1403 end do ! End loop over k-points 1404 1405 end do ! fab: end loop over spin 1406 1407 1408 1409 1410 call xmpi_sum(d3_aux,spaceComm,ierr) 1411 1412 1413 ABI_FREE(s13mat) 1414 ABI_FREE(smat) 1415 ABI_FREE(s1mat) 1416 ABI_FREE(qmat) 1417 ABI_FREE(ipvt) 1418 ABI_FREE(s3mat) 1419 ABI_FREE(zgwork) 1420 ABI_FREE(bd_index) 1421 1422 1423 !fab: I think that in the following we have to make a distinction: 1424 !for the spin unpolarized case we leave the PEAD expression as it is, while 1425 !in the spin polarized case we have simply to divide by 2 1426 !(see eq.19 di PRB 63,155107 [[cite:Nunes2001]], eq. 7 di PRB 71,125107 [[cite:Veithen2005]] 1427 ! and eq 13 di PRB 71, 125107 [[cite:Veithen2005]] ... 1428 !in this latter equation the 2 must be simply replaced by the sum over the spin components... 1429 !and indeed we have inserted the loop over the spin, 1430 !but there was a factor 2 already present in the routine due to spin degenracy that had to be removed) 1431 1432 1433 if (nsppol==1) then 1434 1435 ! Take minus the imaginary part 1436 1437 d3_berry(1,:) = -1_dp*d3_aux(2,:) 1438 d3_berry(2,:) = d3_aux(1,:) 1439 1440 d3_berry(2,:) = 0_dp 1441 1442 else 1443 1444 d3_berry(1,:) = -1_dp*d3_aux(2,:)/2._dp 1445 d3_berry(2,:) = d3_aux(1,:)/2._dp 1446 1447 d3_berry(2,:) = 0_dp/2._dp 1448 1449 end if 1450 1451 !DEBUG 1452 !write(100,*)'pead_nl_mv.f : d3_berry' 1453 !write(100,*)'Perturbation',i1dir,i3dir 1454 !write(100,*) 1455 !write(100,*)'before transformation' 1456 !write(100,*)'real part' 1457 !write(100,'(3(2x,f20.9))')d3_berry(1,:) 1458 !write(100,*) 1459 !write(100,*)'imaginary part' 1460 !write(100,'(3(2x,f20.9))')d3_berry(2,:) 1461 !write(100,*) 1462 !write(100,*)'after transformation' 1463 !ENDDEBUG 1464 1465 !Compute the projection on the basis vectors of 1466 !reciprocal space 1467 1468 d3_aux(:,:) = 0_dp 1469 do ii = 1,3 1470 do jj = 1,3 1471 d3_aux(:,ii) = d3_aux(:,ii) + gmet(ii,jj)*d3_berry(:,jj) 1472 end do 1473 end do 1474 d3_berry(:,:) = d3_aux(:,:) 1475 1476 !Write out the berryphase part of the third order energy 1477 1478 if (mpi_enreg%me == 0) then 1479 1480 write(message,'(a,a,a)')ch10,' Berryphase part of the third-order energy:',ch10 1481 call wrtout(std_out, message,'COLL') 1482 1483 if (i1pert < natom + 1) then 1484 write(message,'(a,i3,a,i3)')& 1485 & ' j1: Displacement of atom ',i1pert,& 1486 & ' along direction ',i1dir 1487 else if (i1pert == natom + 2) then 1488 write(message,'(a,i3)')& 1489 & ' j1: homogenous electric field along direction ',i1dir 1490 end if 1491 call wrtout(std_out, message,'COLL') 1492 1493 write(message,'(a)')& 1494 & ' j2: k-point derivative along direction i2dir ' 1495 call wrtout(std_out, message,'COLL') 1496 1497 if (i3pert < natom + 1) then 1498 write(message,'(a,i3,a,i3,a)')& 1499 & ' j3: Displacement of atom ',i3pert,& 1500 & ' along direction ',i3dir,ch10 1501 else if (i3pert == natom + 2) then 1502 write(message,'(a,i3,a)')& 1503 & ' j3: homogenous electric field along direction ',i3dir,ch10 1504 end if 1505 call wrtout(std_out, message,'COLL') 1506 1507 ! write(ab_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part' 1508 write(std_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part' 1509 do ii = 1,3 1510 write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,& 1511 & d3_berry(1,ii),d3_berry(2,ii) 1512 write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,& 1513 & d3_berry(1,ii),d3_berry(2,ii) 1514 end do 1515 1516 end if ! mpi_enreg%me == 0 1517 1518 !DEBUG 1519 !write(100,*)'real part' 1520 !write(100,'(3(2x,f20.9))')d3_berry(1,:) 1521 !write(100,*) 1522 !write(100,*)'imaginary part' 1523 !write(100,'(3(2x,f20.9))')d3_berry(2,:) 1524 !ENDDEBUG 1525 1526 end subroutine pead_nl_mv
ABINIT/pead_nl_resp [ Functions ]
NAME
pead_nl_resp
FUNCTION
Compute the linear response part to the 3dte
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions cg1 = first-order wavefunction relative to the perturbations i1pert cg3 = first-order wavefunction relative to the perturbations 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 i1dir,i2dir,i3dir=directions of the corresponding perturbations i1pert,i2pert,i3pert = type of perturbation that has to be computed kg(3,mpw*mkmem)=reduced planewave coordinates mband = maximum number of bands mgfft=maximum size of 1D FFTs mkmem = 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 nfft = (effective) number of FFT grid points (for this processor) nkpt = number of k points nspden = number of spin-density components 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 ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information psps <type(pseudopotential_type)> = variables related to pseudopotentials rprimd(3,3) = dimensional primitive translations (bohr) vtrial1(cplex*nfft,nspden)=firs-order local potential xred(3,natom) = reduced atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= spherical harmonics for each G and k point
OUTPUT
d3lo(2,3,mpert,3,mpert,3,mpert) = matrix of the 3DTEs
SOURCE
609 subroutine pead_nl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,& 610 & i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,& 611 & kg,mband,mgfft,mkmem,mk1mem,& 612 & mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,nspinor,nsppol,& 613 & npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm) 614 615 !Arguments ------------------------------------ 616 !scalars 617 integer,intent(in) :: cplex,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband,mgfft 618 integer,intent(in) :: mk1mem,mkmem,mpert,mpsang,mpw,natom,nfft,nkpt,nspden 619 integer,intent(in) :: nspinor,nsppol 620 type(MPI_type),intent(in) :: mpi_enreg 621 type(datafiles_type),intent(in) :: dtfil 622 type(dataset_type),intent(in) :: dtset 623 type(pseudopotential_type),intent(in) :: psps 624 !arrays 625 integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt) 626 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 627 real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol) 628 real(dp),intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol) 629 real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom),rprimd(3,3) 630 real(dp),intent(in) :: xred(3,natom),ylm(mpw*mkmem,mpsang*mpsang*psps%useylm) 631 real(dp),intent(inout) :: vtrial1(cplex*nfft,nspden) 632 real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert) 633 634 !Local variables------------------------------- 635 !scalars 636 integer,parameter :: level=52 637 integer :: bantot,choice,counter,cpopt,dimffnl,iband,icg0,ider,ierr 638 integer :: ii,ikg,ikpt,ilm,ipw,isppol,istwf_k,jband,jj 639 integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nnlout,npw_k 640 integer :: option,paw_opt,signs,spaceComm,tim_fourwf,tim_nonlop 641 real(dp) :: dot1i,dot1r,dot2i,dot2r,doti,dotr,lagi,lagr,sumi,sumr,weight 642 type(gs_hamiltonian_type) :: gs_hamk 643 !arrays 644 integer,allocatable :: kg_k(:,:) 645 real(dp) :: buffer(2),enlout(3),kpq(3),kpt(3) 646 real(dp) :: dum_svectout(1,1),dum(1),rmet(3,3),ylmgr_dum(1,1,1) 647 real(dp),allocatable :: cwave0(:,:),cwavef3(:,:),ffnlk(:,:,:,:) 648 real(dp),allocatable :: gh0(:,:),gh1(:,:),gvnl(:,:),kpg_k(:,:) 649 real(dp),allocatable :: vlocal1(:,:,:),wfraug(:,:,:,:),ylm_k(:,:) 650 type(pawcprj_type) :: cprj_dum(1,1) 651 type(pawtab_type) :: pawtab_dum(0) 652 653 !*********************************************************************** 654 655 ABI_UNUSED(dtfil%ireadwf) 656 657 me = mpi_enreg%me 658 spaceComm=mpi_enreg%comm_cell 659 660 bantot = 0 661 icg0 = 0 662 option = 2 663 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 664 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 665 666 ABI_MALLOC(vlocal1,(cplex*n4,n5,n6)) 667 ABI_MALLOC(wfraug,(2,n4,n5,n6)) 668 669 !Initialize Hamiltonian (k-independent terms) - NCPP only 670 call init_hamiltonian(gs_hamk,psps,pawtab_dum,nspinor,nsppol,nspden,natom,& 671 & dtset%typat,xred,nfft,mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 672 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 673 !& paw_ij=paw_ij) 674 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 675 676 sumr = zero ; sumi = zero 677 678 !Loop over spins 679 680 do isppol = 1, nsppol 681 682 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,vtrial1,vlocal1,option) 683 684 ! Loop over k-points 685 686 ikg = 0 687 do ikpt = 1, nkpt 688 689 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me))cycle 690 691 counter = 100*ikpt 692 693 nband_k = dtset%nband(ikpt+(isppol-1)*nkpt) 694 npw_k = npwarr(ikpt) 695 istwf_k = dtset%istwfk(ikpt) 696 697 kpt(:) = dtset%kptns(:,ikpt) 698 kpq(:) = dtset%kptns(:,ikpt) ! In case of non zero q, kpt = kpt + q 699 700 ABI_MALLOC(cwave0,(2,npw_k*dtset%nspinor)) 701 ABI_MALLOC(cwavef3,(2,npw_k*dtset%nspinor)) 702 ABI_MALLOC(gh0,(2,npw_k*dtset%nspinor)) 703 ABI_MALLOC(gvnl,(2,npw_k*dtset%nspinor)) 704 ABI_MALLOC(gh1,(2,npw_k*dtset%nspinor)) 705 706 ABI_MALLOC(kg_k,(3,npw_k)) 707 ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm)) 708 kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg) 709 if (psps%useylm==1) then 710 do ilm=1,mpsang*mpsang 711 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 712 end do 713 end if 714 715 ! Compute (k+G) and (k+q+G) vectors (only if useylm=1) 716 nkpg=0;if (i2pert<natom+1) nkpg=3*dtset%nloalg(3) 717 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 718 if (nkpg>0) then 719 call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k) 720 end if 721 722 ! Compute nonlocal form factors ffnl at (k+G), for all atoms 723 dimffnl=1 724 ABI_MALLOC(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 725 if (i2pert<natom+1) then 726 ider=0 727 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,& 728 & ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,& 729 & psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,& 730 & psps%usepaw,psps%useylm,ylm_k,ylmgr_dum) 731 end if 732 733 ! Load k-dependent part in the Hamiltonian datastructure 734 call gs_hamk%load_k(kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,& 735 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk,compute_gbound=.true.) 736 ! Load k+q-dependent part in the Hamiltonian datastructure 737 ! call load_kprime_hamiltonian... !! To be activated when q/=0 738 739 ! Loop over bands 740 741 do iband = 1,nband_k 742 743 cwave0(:,:)=cg(:,1+(iband - 1)*npw_k*dtset%nspinor+icg0:& 744 & iband*npw_k*dtset%nspinor+icg0) 745 cwavef3(:,:)=cg3(:,1+(iband-1)*npw_k*dtset%nspinor+icg0:& 746 & iband*npw_k*dtset%nspinor+icg0) 747 748 ! Compute vtrial1 | cwafef3 > 749 tim_fourwf = 0 ; weight = one 750 call fourwf(cplex,vlocal1,cwavef3,gh1,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,& 751 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,& 752 & tim_fourwf,weight,weight,& 753 & gpu_option=dtset%gpu_option) 754 755 ! In case i2pert = phonon-type perturbation 756 ! add first-order change in the nonlocal potential 757 if (i2pert<natom+1) then 758 signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1 759 call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,& 760 & signs,dum_svectout,tim_nonlop,cwavef3,gvnl,iatom_only=i2pert) 761 gh1(:,:) = gh1(:,:) + gvnl(:,:) 762 end if 763 764 ii = (iband-1)*npw_k*dtset%nspinor + icg0 765 call dotprod_g(dotr,doti,istwf_k,npw_k,2,cg1(:,ii+1:ii+npw_k),gh1,mpi_enreg%me_g0,xmpi_comm_self) 766 767 ! Compute vtrial1 | cwave0 > 768 tim_fourwf = 0 ; weight = one 769 call fourwf(cplex,vlocal1,cwave0,gh0,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,& 770 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,& 771 & tim_fourwf,weight,weight,gpu_option=dtset%gpu_option) 772 773 ! In case i2pert = phonon-type perturbation 774 ! add first-order change in the nonlocal potential 775 if (i2pert<natom+1) then 776 signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1 777 call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,& 778 & signs,dum_svectout,tim_nonlop,cwave0,gvnl,iatom_only=i2pert) 779 gh0(:,:) = gh0(:,:) + gvnl(:,:) 780 end if 781 782 ! Compute the dft contribution to the Lagrange multiplier 783 ! cwavef3 and cwave0 have been transferred to gh1 and gh0 784 ! these vectors will be used to store the wavefunctions of band iband 785 ! cg1 and gh0 contain the wavefunctions of band jband 786 787 lagr = zero ; lagi = zero 788 do jband = 1, nband_k 789 790 ii = (jband - 1)*npw_k*dtset%nspinor + icg0 791 jj = (iband - 1)*npw_k*dtset%nspinor + icg0 792 793 ! dot1r and dot1i contain < u_mk | v^(1) | u_nk > 794 ! dot2r and dot2i contain < u_nk^(1) | u_mk^(1) > 795 ! m -> jband and n -> iband 796 797 dot1r = zero ; dot1i = zero 798 dot2r = zero ; dot2i = zero 799 do ipw = 1, npw_k 800 ii = ii + 1 ; jj = jj + 1 801 dot1r = dot1r + cg(1,ii)*gh0(1,ipw) + cg(2,ii)*gh0(2,ipw) 802 dot1i = dot1i + cg(1,ii)*gh0(2,ipw) - cg(2,ii)*gh0(1,ipw) 803 dot2r = dot2r + cg1(1,jj)*cg3(1,ii) + & 804 & cg1(2,jj)*cg3(2,ii) 805 dot2i = dot2i + cg1(1,jj)*cg3(2,ii) - & 806 & cg1(2,jj)*cg3(1,ii) 807 end do ! ipw 808 809 lagr = lagr + dot1r*dot2r - dot1i*dot2i 810 lagi = lagi + dot1r*dot2i + dot1i*dot2r 811 812 end do ! jband 813 814 sumr = sumr + & 815 & dtset%wtk(ikpt)*occ(bantot+iband)*(dotr-lagr) 816 sumi = sumi + & 817 & dtset%wtk(ikpt)*occ(bantot+iband)*(doti-lagi) 818 819 end do ! end loop over bands 820 821 bantot = bantot + nband_k 822 icg0 = icg0 + npw_k*dtset%nspinor*nband_k 823 ikg = ikg + npw_k 824 825 ABI_FREE(cwave0) 826 ABI_FREE(cwavef3) 827 ABI_FREE(gh0) 828 ABI_FREE(gh1) 829 ABI_FREE(gvnl) 830 ABI_FREE(kg_k) 831 ABI_FREE(ylm_k) 832 ABI_FREE(ffnlk) 833 ABI_FREE(kpg_k) 834 835 end do ! end loop over k-points 836 837 end do ! end loop over spins 838 839 if (xmpi_paral == 1) then 840 buffer(1) = sumr ; buffer(2) = sumi 841 call xmpi_sum(buffer,spaceComm,ierr) 842 sumr = buffer(1) ; sumi = buffer(2) 843 end if 844 845 846 d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr 847 !d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi 848 849 !In some cases, the imaginary part is /= 0 because of the 850 !use of time reversal symmetry 851 d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero 852 853 call gs_hamk%free() 854 855 ABI_FREE(vlocal1) 856 ABI_FREE(wfraug) 857 858 end subroutine pead_nl_resp