TABLE OF CONTENTS


ABINIT/forstr [ Functions ]

[ Top ] [ Functions ]

NAME

 forstr

FUNCTION

 Drives the computation of forces and/or stress tensor

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  cg(2,mcg)=wavefunctions (may be read from disk instead of input)
  cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
  dtefield <type(efield_type)> = variables related to Berry phase
  dtset <type(dataset_type)>=all input variables in this dataset
   | berryopt  = 4: electric field is on -> add the contribution of the
   |                - \Omega E.P term to the total energy
   |          /= 4: electric field is off
   |  from Etot(npw) data (at fixed geometry), used for making
   |  Pulay correction to stress tensor (hartree).  Should be <=0.
   | ecut=cut-off energy for plane wave basis sphere (Ha)
   | ecutsm=smearing energy for plane wave kinetic energy (Ha)
   | effmass_free=effective mass for electrons (1. in common case)
   | efield = cartesian coordinates of the electric field in atomic units
   | ionmov=governs the movement of atoms (see help file)
   | densfor_pred=governs the mixed electronic-atomic part of the preconditioner
   | istwfk(nkpt)=input option parameter that describes the storage of wfs
   | kptns(3,nkpt)=reduced coordinates of k points in Brillouin zone
   | mband=maximum number of bands
   | mgfft=maximum size of 1D FFTs
   | mkmem=maximum number of k points in core memory
   | mpw = maximum number of plane waves
   | natom=number of atoms in cell
   | nband(nkpt*nsppol)=number of bands to be included in summation at each k point
   | nfft=(effective) number of FFT grid points (for this processor)
   | ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
   | nkpt=number of k points in Brillouin zone
   | nloalg(3)=governs the choice of the algorithm for non-local operator.
   | nspden=number of spin-density components
   | nsppol=1 for unpolarized, 2 for spin-polarized
   | pawprtvol=control print volume and debugging output for PAW
   | prtvol=integer controlling volume of printed output
   | symafm(nsym)=(anti)ferromagnetic part of symmetry operations
   | tfkinfunc=1 if use of Thomas-Fermi kinetic functional
   |          =2 if use of recursion method
   | typat(natom)=type integer for each atom in cell
   | wtk(nkpt)=weights associated with various k points
   | nsym=number of symmetries in space group
  energies <type(energies_type)>=all part of total energy.
   | e_localpsp(IN)=local psp energy (hartree)
   | e_hartree(IN)=Hartree part of total energy (hartree units)
   | e_corepsp(IN)=psp core-core energy
   | e_kinetic(IN)=kinetic energy part of total energy.
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  fock <type(fock_type)>= quantities to calculate Fock exact exchange
  grchempottn(3,natom)=d(E_chemical potential)/d(xred) (hartree)
  grcondft(3,natom)=d(E_constrainedDFT)/d(xred) (hartree)
  grewtn(3,natom)=d(Ewald)/d(xred) (hartree)
  grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree)
  gsqcut=cutoff value on G**2 for (large) sphere inside FFT box.
                       gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2)
  indsym(4,nsym,natom)=index showing transformation of atom labels
                       under symmetry operations (computed in symatm)
  kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere
  kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mgfftf= -PAW ONLY- maximum size of 1D FFTs for the fine grid
         (mgfftf=mgfft for norm-conserving potential runs)
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  n3xccc=dimension of the xccc3d array (0 or nfftf).
  nattyp(ntypat)=number of atoms of each type
  nfftf= -PAW ONLY- number of FFT grid points for the fine grid
         (nfftf=nfft for norm-conserving potential runs)
  ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid
              (ngfftf=ngfft for norm-conserving potential runs)
  ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
  nhat(nfftf,nspden*psps%usepaw)= -PAW only- compensation density
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  npwarr(nkpt)=number of planewaves in basis and on boundary for each k
  ntypat=number of types of atoms
  nvresid(nfftf,nspden)=array for the residual of the density/potential
  occ(mband*nkpt*nsppol)=occupancies of bands at various k points
  optfor=1 if computation of forces is required
  optres=0 if the potential residual has to be used for forces corrections
        =1 if the density residual has to be used for forces corrections
  paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
  pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases
  ph1df(2,3*(2*mgfftf+1)*natom)=-PAW only- 1-dim structure factor phases for the fine grid
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
   | mpsang= 1+maximum angular momentum
  rhog(2,nfftf)=Fourier transform of charge density (bohr^-3)
  rhor(nfftf,nspden)=array for electron density in electrons/bohr**3.
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  strscondft(6)=cDFT correction to stress
  strsxc(6)=xc correction to stress
  stress_needed=1 if computation of stress tensor is required
  symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
  ucvol=unit cell volume in bohr**3
  usecprj=1 if cprj datastructure is stored in memory
  vhartr(nfftf)=array for holding Hartree potential
  vpsp(nfftf)=array for holding local psp
  vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space
  vxctau(nfft,nspden,4*usekden)=(only for meta-GGA): derivative of XC energy density
                                wrt kinetic energy density (depsxcdtau)
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
  xcctau3d(n3xccc*usekden)=(only for meta-GGA): 3D core electron kinetic energy density for XC core correction
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics

OUTPUT

  ==== if (optfor==1) ====
   diffor=maximal absolute value of changes in the components of
          force between the input and the output.
   favg(3)=mean of the forces before correction for translational symmetry
   fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
     at input, previous value of forces,
     at output, new value.
     Note : unlike gred, this array has been corrected by enforcing
     the translational symmetry, namely that the sum of force
     on all atoms is zero.
   forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
   gred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
   gresid(3,natom)=forces due to the residual of the density/potential
   grhf(3,natom)=Hellman-Feynman derivatives of the total energy
   grxc(9+3*natom)=d(Exc)/d(xred) if core charges are used
   maxfor=maximal absolute value of the output array force.
   synlgr(3,natom)=symmetrized gradients of energy due to nonlocal contributions
  ==== if (stress_needed==1) ====
   strten(6)=components of the stress tensor (hartree/bohr^3) for the
    6 unique components of this symmetric 3x3 tensor:
    Given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).

SIDE EFFECTS

  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
  forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
  ===== if psps%usepaw==1
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
    (gradients of rhoij for each atom with respect to atomic positions are computed here)
  wvl <type(wvl_data)>=all wavelets data.

NOTES

  Be careful to the meaning of nfft (size of FFT grids):
   - In case of norm-conserving calculations the FFT grid is the usual FFT grid.
   - In case of PAW calculations:
     Two FFT grids are used; one with nfft points (coarse grid) for
     the computation of wave functions ; one with nfftf points
     (fine grid) for the computation of total density.

SOURCE

245 subroutine forstr(atindx1,cg,cprj,diffor,dtefield,dtset,eigen,electronpositron,energies,favg,fcart,fock,&
246 &                 forold,gred,grchempottn,grcondft,gresid,grewtn,grhf,grvdw,grxc,gsqcut,extfpmd,indsym,&
247 &                 kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,n3xccc,nattyp,&
248 &                 nfftf,ngfftf,ngrvdw,nhat,nkxc,npwarr,&
249 &                 ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,&
250 &                 pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,psps,rhog,rhor,rprimd,stress_needed,&
251 &                 strscondft,strsxc,strten,symrec,synlgr,ucvol,usecprj,vhartr,vpsp,&
252 &                 vxc,vxctau,wvl,xccc3d,xcctau3d,xred,ylm,ylmgr,qvpotzero)
253 
254 !Arguments ------------------------------------
255 !scalars
256  integer,intent(in) :: mcg,mcprj,mgfftf,my_natom,n3xccc,nfftf,ngrvdw,nkxc,ntypat,optfor,optres
257  integer,intent(in) :: stress_needed,usecprj
258  real(dp),intent(in) :: gsqcut,qvpotzero,ucvol
259  real(dp),intent(inout) :: diffor,maxfor
260  type(electronpositron_type),pointer :: electronpositron
261  type(MPI_type),intent(inout) :: mpi_enreg
262  type(efield_type),intent(in) :: dtefield
263  type(dataset_type),intent(in) :: dtset
264  type(energies_type),intent(in) :: energies
265  type(extfpmd_type),pointer,intent(inout) :: extfpmd
266  type(pawang_type),intent(in) :: pawang
267  type(pawfgr_type),intent(in) :: pawfgr
268  type(pseudopotential_type),intent(in) :: psps
269  type(wvl_data),intent(inout) :: wvl
270  type(fock_type),pointer, intent(inout) :: fock
271 !arrays
272  integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
273  integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfftf(18)
274  integer,intent(in) :: npwarr(dtset%nkpt),symrec(3,3,dtset%nsym)
275  real(dp),intent(in) :: cg(2,mcg)
276  real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
277  real(dp),intent(in) :: grchempottn(3,dtset%natom),grcondft(3,dtset%natom),grewtn(3,dtset%natom)
278  real(dp),intent(in) :: grvdw(3,ngrvdw),kxc(dtset%nfft,nkxc)
279  real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
280  real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
281  real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom)
282  real(dp),intent(in) :: rhog(2,nfftf),strscondft(6),strsxc(6),vhartr(nfftf)
283  real(dp),intent(in) :: vpsp(nfftf),vxc(nfftf,dtset%nspden),vxctau(nfftf,dtset%nspden,4*dtset%usekden)
284  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
285  real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
286  real(dp),intent(inout) :: forold(3,dtset%natom)
287  real(dp),intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw),rhor(nfftf,dtset%nspden),rprimd(3,3)
288  real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*dtset%usekden),xred(3,dtset%natom)
289  real(dp),intent(inout),target :: nvresid(nfftf,dtset%nspden)
290  real(dp),intent(out) :: favg(3)
291  real(dp),intent(inout) :: fcart(3,dtset%natom),gred(3,dtset%natom)
292  real(dp),intent(inout) :: gresid(3,dtset%natom),grhf(3,dtset%natom)
293  real(dp),intent(inout) :: grxc(3,dtset%natom),strten(6),synlgr(3,dtset%natom)
294  type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj)
295  type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
296  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
297  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
298  type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw)
299  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
300 
301 !Local variables-------------------------------
302 !scalars
303 
304  integer :: comm_grid,ifft,ispden,ncpgr,occopt_,optgr,optgr2,option,optnc,optstr,optstr2,iorder_cprj,ctocprj_choice
305  integer :: idir,iatom,unpaw,mcgbz
306  integer,allocatable :: dimcprj(:)
307  real(dp) ::dum,dum1,ucvol_
308  logical :: apply_residual
309 !arrays
310  real(dp),parameter :: k0(3)=(/zero,zero,zero/)
311  real(dp) :: kinstr(6),nlstr(6),tsec(2),strdum(6),gmet(3,3),gprimd(3,3),rmet(3,3)
312  real(dp) :: dummy(0)
313  real(dp),allocatable :: grnl(:),vlocal(:,:),vxc_hf(:,:),xcart(:,:),ylmbz(:,:),ylmgrbz(:,:,:)
314  real(dp), ABI_CONTIGUOUS pointer :: resid(:,:)
315 
316 ! *************************************************************************
317 
318  call timab(910,1,tsec)
319  call timab(911,1,tsec)
320 
321 !Do nothing if nothing is required
322  if (optfor==0.and.stress_needed==0) return
323 
324 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
325  if (dtset%usewvl==0) then
326    if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
327      ABI_BUG(' wrong values for nfft, nfftf !')
328    end if
329    if ((psps%usepaw==1.and.pawfgr%mgfft/=mgfftf).or.(psps%usepaw==0.and.dtset%mgfft/=mgfftf)) then
330      ABI_BUG('wrong values for mgfft, mgfftf!')
331    end if
332  end if
333 
334 !==========================================================================
335 !Here compute terms common to forces and stresses
336 !==========================================================================
337 
338  !output only if (optfor==1) but we have to allocate it
339  ABI_MALLOC(grnl,(3*dtset%natom*optfor))
340  grnl(:)=zero
341 
342 !Compute nonlocal psp + potential Fock ACE parts of forces and stress tensor
343 !-involves summations over wavefunctions at all k points
344  if (dtset%tfkinfunc>0.and.stress_needed==1) then
345    kinstr(1:3)=-two/three*energies%e_kinetic/ucvol ; kinstr(4:6)=zero
346    nlstr(1:6)=zero
347  else if (dtset%usewvl==0) then
348    occopt_=0 ! This means that occ are now fixed
349    if(dtset%usefock==1 .and. associated(fock)) then
350 !     if((dtset%optstress/=0).and.(psps%usepaw==1)) then
351      if((psps%usepaw==1).and.((dtset%optstress/=0).or.(dtset%optforces==2))) then
352        if(dtset%optstress==0) then
353          ctocprj_choice=2
354          ncpgr=3
355        end if
356        if(dtset%optstress/=0) then
357          ctocprj_choice=20*optfor+3*dtset%optstress
358          ncpgr=6*dtset%optstress+3*optfor
359        end if
360        if (allocated(fock%fock_BZ%cwaveocc_prj)) then
361          call pawcprj_free(fock%fock_BZ%cwaveocc_prj)
362          ABI_FREE(fock%fock_BZ%cwaveocc_prj)
363          ABI_MALLOC(fock%fock_BZ%cwaveocc_prj,(dtset%natom,fock%fock_BZ%mcprj))
364          ABI_MALLOC(dimcprj,(dtset%natom))
365          call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
366          call pawcprj_alloc(fock%fock_BZ%cwaveocc_prj,ncpgr,dimcprj)
367          ABI_FREE(dimcprj)
368        end if
369        iatom=-1;idir=0;iorder_cprj=0;unpaw=26
370        call metric(gmet,gprimd,-1,rmet,rprimd,dum)
371        if (fock%fock_BZ%mkpt/=dtset%mkmem.or.(fock%fock_BZ%mpi_enreg%paral_hf ==1)) then
372          ABI_MALLOC(ylmbz,(dtset%mpw*fock%fock_BZ%mkpt,psps%mpsang*psps%mpsang*psps%useylm))
373          ABI_MALLOC(ylmgrbz,(dtset%mpw*fock%fock_BZ%mkpt,3,psps%mpsang*psps%mpsang*psps%useylm))
374          option=1; mcgbz=dtset%mpw*fock%fock_BZ%mkptband*fock%fock_common%my_nsppol
375          call initylmg(gprimd,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,fock%fock_BZ%mkpt,fock%fock_BZ%mpi_enreg,&
376 &         psps%mpsang,dtset%mpw,fock%fock_BZ%nbandocc_bz,fock%fock_BZ%mkpt,&
377 &         fock%fock_BZ%npwarr,dtset%nsppol,option,rprimd,ylmbz,ylmgrbz)
378          call ctocprj(fock%fock_common%atindx,fock%fock_BZ%cgocc,ctocprj_choice,fock%fock_BZ%cwaveocc_prj,gmet,gprimd,iatom,idir,&
379 &         iorder_cprj,fock%fock_BZ%istwfk_bz,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,mcgbz,&
380 &         fock%fock_BZ%mcprj,dtset%mgfft,fock%fock_BZ%mkpt,fock%fock_BZ%mpi_enreg,psps%mpsang,&
381 &         dtset%mpw,dtset%natom,nattyp,fock%fock_BZ%nbandocc_bz,dtset%natom,dtset%ngfft,fock%fock_BZ%mkpt,&
382 &         dtset%nloalg,fock%fock_BZ%npwarr,dtset%nspinor,&
383 &         dtset%nsppol,fock%fock_common%my_nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,unpaw,&
384 &         xred,ylmbz,ylmgrbz)
385          ABI_FREE(ylmbz)
386          ABI_FREE(ylmgrbz)
387        else
388          call ctocprj(fock%fock_common%atindx,fock%fock_BZ%cgocc,ctocprj_choice,fock%fock_BZ%cwaveocc_prj,gmet,gprimd,iatom,idir,&
389 &         iorder_cprj,fock%fock_BZ%istwfk_bz,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,mcg,&
390 &         fock%fock_BZ%mcprj,dtset%mgfft,fock%fock_BZ%mkpt,mpi_enreg,psps%mpsang,&
391 &         dtset%mpw,dtset%natom,nattyp,fock%fock_BZ%nbandocc_bz,dtset%natom,dtset%ngfft,fock%fock_BZ%mkpt,&
392 &         dtset%nloalg,fock%fock_BZ%npwarr,dtset%nspinor,&
393 &         dtset%nsppol,fock%fock_common%my_nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,unpaw,&
394 &         xred,ylm,ylmgr)
395        end if
396      end if
397    end if
398    call forstrnps(cg,cprj,dtset%ecut,dtset%ecutsm,dtset%effmass_free,eigen,electronpositron,fock,grnl,&
399 &   dtset%istwfk,kg,kinstr,nlstr,dtset%kptns,dtset%mband,mcg,mcprj,dtset%mgfft,dtset%mkmem,&
400 &   mpi_enreg,psps%mpsang,dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,dtset%ngfft,&
401 &   dtset%nkpt,dtset%nloalg,npwarr,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%nsym,ntypat,&
402 &   dtset%nucdipmom,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,stress_needed,symrec,dtset%typat,&
403 &   usecprj,dtset%usefock,dtset%use_gpu_cuda,dtset%wtk,xred,ylm,ylmgr)
404 !DEBUG
405 !   write(6,*)' after forstrnps, nlstr=',nlstr(1:6)
406 !ENDDEBUG
407  else if (optfor>0) then !WVL
408    ABI_MALLOC(xcart,(3, dtset%natom))
409    call xred2xcart(dtset%natom, rprimd, xcart, xred)
410    call wvl_nl_gradient(grnl, mpi_enreg, dtset%natom, rprimd, wvl, xcart)
411    ABI_FREE(xcart)
412  end if
413 
414  call timab(911,2,tsec)
415  call timab(912,1,tsec)
416 
417 !PAW: add gradients due to Dij derivatives to non-local term
418  if (psps%usepaw==1) then
419    ABI_MALLOC(vlocal,(nfftf,dtset%nspden))
420 
421 !$OMP PARALLEL DO COLLAPSE(2)
422    do ispden=1,min(dtset%nspden,2)
423      do ifft=1,nfftf
424        vlocal(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft)
425      end do
426    end do
427    if (dtset%nspden==4) then
428 !$OMP PARALLEL DO COLLAPSE(2)
429      do ispden=3,4
430        do ifft=1,nfftf
431          vlocal(ifft,ispden)=vxc(ifft,ispden)
432        end do
433      end do
434    end if
435    ucvol_=ucvol
436 #if defined HAVE_BIGDFT
437    if (dtset%usewvl==1) ucvol_=product(wvl%den%denspot%dpbox%hgrids)*real(product(wvl%den%denspot%dpbox%ndims),dp)
438 #endif
439    optgr=optfor;optgr2=0;optstr=stress_needed;optstr2=0
440    comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl
441    call pawgrnl(atindx1,dtset%nspden,dummy,1,dummy,grnl,gsqcut,mgfftf,my_natom,dtset%natom,&
442 &   nattyp,nfftf,ngfftf,nhat,nlstr,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,&
443 &   pawang,pawfgrtab,pawrhoij,pawtab,ph1df,psps,k0,rprimd,symrec,dtset%typat,ucvol_,vlocal,vxc,xred,&
444 &   mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom,mpi_comm_grid=comm_grid)
445 !DEBUG
446 !   write(6,*)' after pawgrnl, nlstr=',nlstr(1:6)
447 !ENDDEBUG
448    ABI_FREE(vlocal)
449 
450  end if
451  call timab(912,2,tsec)
452  call timab(913,1,tsec)
453 
454 !==========================================================================
455 !Here compute forces (if required)
456 !==========================================================================
457  if (optfor==1) then
458    apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. &
459 &   abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5)
460 
461 !  If residual is a density residual (and forces from residual asked),
462 !  has to convert it into a potential residual before calling forces routine
463    if (apply_residual) then
464      ABI_MALLOC(resid,(nfftf,dtset%nspden))
465      option=0; if (dtset%densfor_pred<0) option=1
466      optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2
467      call nres2vres(dtset,gsqcut,psps%usepaw,kxc,mpi_enreg,my_natom,nfftf,ngfftf,nhat,&
468 &     nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,&
469 &     rhor,rprimd,psps%usepaw,resid,xccc3d,xred,vxc)
470    else
471      resid => nvresid
472    end if
473 
474    call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,gred,grchempottn,grcondft,gresid,grewtn,&
475 &   grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfftf,&
476 &   mpi_enreg,psps%n1xccc,n3xccc,nattyp,&
477 &   nfftf,ngfftf,ngrvdw,ntypat,pawrad,pawtab,ph1df,psps,rhog,&
478 &   rhor,rprimd,symrec,synlgr,dtset%usefock,resid,vxc,vxctau,wvl%descr,wvl%den,xred,&
479 &   electronpositron=electronpositron)
480 
481    if (apply_residual) then
482      ABI_FREE(resid)
483    end if
484  end if
485 
486  call timab(913,2,tsec)
487  call timab(914,1,tsec)
488 
489 !==========================================================================
490 !Here compute stress tensor (if required)
491 !==========================================================================
492 
493  if (stress_needed==1.and.dtset%usewvl==0) then
494 !   if (dtset%usefock==1 .and. associated(fock).and.fock%fock_common%optstr.and.psps%usepaw==0) then
495    if (dtset%usefock==1 .and. associated(fock)) then
496      if (fock%fock_common%optstr) then
497        fock%fock_common%stress(1:3)=fock%fock_common%stress(1:3)-(two*energies%e_fock-energies%e_fock0)/ucvol
498        if (n3xccc>0.and.psps%usepaw==0 .and. &
499 &        (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) then
500          ABI_MALLOC(vxc_hf,(nfftf,dtset%nspden))
501 !compute Vxc^GGA(rho_val)
502          call xchybrid_ncpp_cc(dtset,dum,mpi_enreg,nfftf,ngfftf,n3xccc,rhor,rprimd,strdum,dum1,xccc3d,vxc=vxc_hf,optstr=1)
503        end if
504      end if
505    end if
506    call stress(atindx1,dtset%berryopt,dtefield,energies%e_localpsp,dtset%efield,&
507 &   energies%e_hartree,energies%e_corepsp,fock,gsqcut,extfpmd,dtset%ixc,kinstr,mgfftf,&
508 &   mpi_enreg,psps%mqgrid_vl,psps%n1xccc,n3xccc,dtset%natom,nattyp,&
509 &   nfftf,ngfftf,nlstr,dtset%nspden,dtset%nsym,ntypat,psps,pawrad,pawtab,ph1df,&
510 &   dtset%prtvol,psps%qgrid_vl,dtset%red_efieldbar,rhog,rprimd,strten,strscondft,strsxc,symrec,&
511 &   dtset%typat,dtset%usefock,dtset%usekden,psps%usepaw,&
512 &   dtset%vdw_tol,dtset%vdw_tol_3bt,dtset%vdw_xc,psps%vlspl,vxc,vxctau,vxc_hf,&
513 &   psps%xccc1d,xccc3d,xcctau3d,psps%xcccrc,xred,psps%ziontypat,psps%znucltypat,qvpotzero,&
514 &   electronpositron=electronpositron)
515  end if
516 
517 !Memory deallocation
518  ABI_FREE(grnl)
519  if (allocated(vxc_hf)) then
520    ABI_FREE(vxc_hf)
521  end if
522 
523 
524  call timab(914,2,tsec)
525  call timab(910,2,tsec)
526 
527 end subroutine forstr

ABINIT/forstrnps [ Functions ]

[ Top ] [ Functions ]

NAME

 forstrnps

FUNCTION

 Compute nonlocal pseudopotential energy contribution to forces and/or stress tensor
 as well as kinetic energy contribution to stress tensor.

INPUTS

  cg(2,mcg)=wavefunctions (may be read from disk file)
  cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
  ecut=cut-off energy for plane wave basis sphere (Ha)
  ecutsm=smearing energy for plane wave kinetic energy (Ha)
  effmass_free=effective mass for electrons (1. in common case)
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  fock <type(fock_type)>= quantities to calculate Fock exact exchange
  istwfk(nkpt)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis
  kpt(3,nkpt)=k points in reduced coordinates
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mgfft=maximum size of 1D FFTs
  mkmem=number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mpw= maximum number of plane waves
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nband(nkpt)=number of bands at each k point
  nfft=number of FFT grid points
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt=number of k points in Brillouin zone
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npwarr(nkpt)=number of planewaves in basis and boundary at each k
  nspden=Number of spin Density components
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym=number of elements in symmetry group
  ntypat=number of types of atoms
  nucdipmom(3,my_natom)= nuclear dipole moments
  occ(mband*nkpt*nsppol)=occupation numbers for each band over all k points
  optfor=1 if computation of forces is required
  paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  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 in real space (bohr)
  stress_needed=1 if computation of stress tensor is required
  symrec(3,3,nsym)=symmetries in reciprocal space (dimensionless)
  typat(natom)=type of each atom
  usecprj=1 if cprj datastructure has been allocated
  use_gpu_cuda= 0 or 1 to know if we use cuda for nonlop call
  wtk(nkpt)=weight associated with each k point
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics

OUTPUT

  if (optfor==1)
   grnl(3*natom*optfor)=stores grads of nonlocal energy wrt atomic coordinates
  if (stress_needed==1)
   kinstr(6)=kinetic energy part of stress tensor (hartree/bohr^3)
   Store 6 unique components of symmetric 3x3 tensor in the order
   11, 22, 33, 32, 31, 21
   npsstr(6)=nonlocal pseudopotential energy part of stress tensor
    (hartree/bohr^3)

SOURCE

 601 subroutine forstrnps(cg,cprj,ecut,ecutsm,effmass_free,eigen,electronpositron,fock,&
 602 &  grnl,istwfk,kg,kinstr,npsstr,kpt,mband,mcg,mcprj,mgfft,mkmem,mpi_enreg,mpsang,&
 603 &  mpw,my_natom,natom,nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,nsym,&
 604 &  ntypat,nucdipmom,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,&
 605 &  stress_needed,symrec,typat,usecprj,usefock,use_gpu_cuda,wtk,xred,ylm,ylmgr)
 606 
 607 !Arguments ------------------------------------
 608 !scalars
 609  integer,intent(in) :: mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,nfft,nkpt
 610  integer,intent(in) :: nspden,nsppol,nspinor,nsym,ntypat,optfor,stress_needed
 611  integer,intent(in) :: usecprj,usefock,use_gpu_cuda
 612  real(dp),intent(in) :: ecut,ecutsm,effmass_free
 613  type(electronpositron_type),pointer :: electronpositron
 614  type(MPI_type),intent(inout) :: mpi_enreg
 615  type(pseudopotential_type),intent(in) :: psps
 616 !arrays
 617  integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol)
 618  integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt)
 619  integer,intent(in) :: symrec(3,3,nsym),typat(natom)
 620  real(dp),intent(in) :: cg(2,mcg)
 621  real(dp),intent(in) :: eigen(mband*nkpt*nsppol),kpt(3,nkpt),nucdipmom(3,my_natom)
 622  real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom)
 623  real(dp),intent(in) :: rprimd(3,3),wtk(nkpt),xred(3,natom)
 624  real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
 625  real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm)
 626  real(dp),intent(out) :: grnl(3*natom*optfor),kinstr(6),npsstr(6)
 627  type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj)
 628  type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
 629  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
 630  type(fock_type),pointer, intent(inout) :: fock
 631 !Local variables-------------------------------
 632 !scalars
 633  integer,parameter :: tim_rwwf=7
 634  integer :: bandpp,bdtot_index,choice,cpopt,dimffnl,dimffnl_str,iband,iband_cprj,iband_last,ibg,icg,ider,ider_str
 635  integer :: idir,idir_str,ierr,ii,ikg,ikpt,ilm,ipositron,ipw,ishift,isppol,istwf_k
 636  integer :: mband_cprj,me_distrb,my_ikpt,my_nspinor,nband_k,nband_cprj_k,ndat,nkpg
 637  integer :: nnlout,npw_k,paw_opt,signs,spaceComm
 638  integer :: tim_nonlop,tim_nonlop_prep,usecprj_local,use_ACE_old
 639  integer :: blocksize,iblock,iblocksize,ibs,nblockbd
 640  real(dp) :: ar,renorm_factor,dfsm,ecutsm_inv,fact_kin,fsm,htpisq,kgc1
 641  real(dp) :: kgc2,kgc3,kin,xx
 642  type(gs_hamiltonian_type) :: gs_hamk
 643  logical :: compute_gbound,usefock_loc
 644  character(len=500) :: msg
 645  type(fock_common_type),pointer :: fockcommon
 646 !arrays
 647  integer,allocatable :: kg_k(:,:)
 648  real(dp) :: kpoint(3),nonlop_dum(1,1),rmet(3,3),tsec(2)
 649  real(dp),allocatable :: cwavef(:,:),enlout(:),ffnl_sav(:,:,:,:),ffnl_str(:,:,:,:)
 650  real(dp),allocatable :: ghc_dum(:,:),gprimd(:,:),kpg_k(:,:),kpg_k_sav(:,:)
 651  real(dp),allocatable :: kstr1(:),kstr2(:),kstr3(:),kstr4(:),kstr5(:),kstr6(:)
 652  real(dp),allocatable :: lambda(:),occblock(:),ph3d(:,:,:),ph3d_sav(:,:,:)
 653  real(dp),allocatable :: weight(:),ylm_k(:,:),ylmgr_k(:,:,:)
 654  real(dp),allocatable,target :: ffnl(:,:,:,:)
 655  type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
 656  type(pawcprj_type),target,allocatable :: cwaveprj(:,:)
 657  type(pawcprj_type),pointer :: cwaveprj_idat(:,:)
 658 
 659 
 660 !*************************************************************************
 661 
 662  call timab(920,1,tsec) ; call timab(921,-1,tsec)
 663 
 664 !Init mpicomm and me
 665  if(mpi_enreg%paral_kgb==1)then
 666    spaceComm=mpi_enreg%comm_kpt
 667    me_distrb=mpi_enreg%me_kpt
 668  else
 669 !* In case of HF calculation
 670    if (mpi_enreg%paral_hf==1) then
 671      spaceComm=mpi_enreg%comm_kpt
 672      me_distrb=mpi_enreg%me_kpt
 673    else
 674      spaceComm=mpi_enreg%comm_cell
 675      me_distrb=mpi_enreg%me_cell
 676    end if
 677  end if
 678 
 679 !Some constants
 680  ipositron=abs(electronpositron_calctype(electronpositron))
 681  my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
 682 !Smearing of plane wave kinetic energy
 683  ecutsm_inv=zero;if( ecutsm>1.0d-20) ecutsm_inv=1/ecutsm
 684 !htpisq is (1/2) (2 Pi) **2:
 685  htpisq=0.5_dp*(two_pi)**2
 686 
 687 !Check that fock is present if want to use fock option
 688  compute_gbound=.false.
 689  usefock_loc = (usefock==1 .and. associated(fock))
 690 !Arrays initializations
 691  grnl(:)=zero
 692  if (usefock_loc) then
 693    fockcommon =>fock%fock_common
 694    fockcommon%optfor=.false.
 695    fockcommon%optstr=.false.
 696    use_ACE_old=fockcommon%use_ACE
 697    fockcommon%use_ACE=0
 698    if (fockcommon%optfor) compute_gbound=.true.
 699  end if
 700  if (stress_needed==1) then
 701    kinstr(:)=zero;npsstr(:)=zero
 702    if (usefock_loc) then
 703      fockcommon%optstr=.TRUE.
 704      fockcommon%stress=zero
 705      compute_gbound=.true.
 706    end if
 707  end if
 708 
 709  usecprj_local=usecprj
 710 
 711  if ((usefock_loc).and.(psps%usepaw==1)) then
 712    usecprj_local=1
 713    if(optfor==1)then
 714      fockcommon%optfor=.true.
 715      if (.not.allocated(fockcommon%forces_ikpt)) then
 716        ABI_MALLOC(fockcommon%forces_ikpt,(3,natom,mband))
 717      end if
 718      if (.not.allocated(fockcommon%forces)) then
 719        ABI_MALLOC(fockcommon%forces,(3,natom))
 720      end if
 721      fockcommon%forces=zero
 722      compute_gbound=.true.
 723    end if
 724  end if
 725 
 726 !Initialize Hamiltonian (k-independent terms)
 727 
 728  call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,nsppol,nspden,natom,&
 729 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,usecprj=usecprj_local,&
 730 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
 731 & paw_ij=paw_ij,ph1d=ph1d,electronpositron=electronpositron,fock=fock,&
 732 & nucdipmom=nucdipmom,use_gpu_cuda=use_gpu_cuda)
 733  rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
 734 
 735 !need to reorder cprj=<p_lmn|Cnk> (from unsorted to atom-sorted)
 736  if (psps%usepaw==1.and.usecprj_local==1) then
 737    call pawcprj_reorder(cprj,gs_hamk%atindx)
 738  end if
 739 
 740 !Common data for "nonlop" routine
 741  signs=1 ; idir=0  ; ishift=0 ; tim_nonlop=4 ; tim_nonlop_prep=12
 742  choice=2*optfor;if (stress_needed==1) choice=10*choice+3
 743  if (optfor==1.and.stress_needed==1)  ishift=6
 744  nnlout=max(1,6*stress_needed+3*natom*optfor)
 745  if (psps%usepaw==0) then
 746    paw_opt=0 ; cpopt=-1
 747  else
 748    paw_opt=2 ; cpopt=-1+3*usecprj_local
 749  end if
 750 
 751  call timab(921,2,tsec)
 752 
 753 !LOOP OVER SPINS
 754  bdtot_index=0;ibg=0;icg=0
 755  do isppol=1,nsppol
 756 
 757 !  Continue to initialize the Hamiltonian (PAW DIJ coefficients)
 758    call gs_hamk%load_spin(isppol,with_nonlocal=.true.)
 759    if (usefock_loc) fockcommon%isppol=isppol
 760 
 761 !  Loop over k points
 762    ikg=0
 763    do ikpt=1,nkpt
 764      if (usefock_loc) fockcommon%ikpt=ikpt
 765      nband_k=nband(ikpt+(isppol-1)*nkpt)
 766      istwf_k=istwfk(ikpt)
 767      npw_k=npwarr(ikpt)
 768      kpoint(:)=kpt(:,ikpt)
 769 
 770      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
 771        bdtot_index=bdtot_index+nband_k
 772        cycle
 773      end if
 774 
 775      call timab(922,1,tsec)
 776 
 777 !    Parallelism over FFT and/or bands: define sizes and tabs
 778      if (mpi_enreg%paral_kgb==1) then
 779        my_ikpt=mpi_enreg%my_kpttab(ikpt)
 780        nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
 781        bandpp=mpi_enreg%bandpp
 782        my_bandfft_kpt => bandfft_kpt(my_ikpt)
 783      else
 784        my_ikpt=ikpt
 785        bandpp=mpi_enreg%bandpp
 786        nblockbd=nband_k/bandpp
 787      end if
 788      blocksize=nband_k/nblockbd
 789      mband_cprj=mband/mpi_enreg%nproc_band
 790      nband_cprj_k=nband_k/mpi_enreg%nproc_band
 791 
 792      ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize))
 793      if (psps%usepaw==1.and.usecprj_local==1) then
 794        ABI_MALLOC(cwaveprj,(natom,my_nspinor*bandpp))
 795        call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
 796      else
 797        ABI_MALLOC(cwaveprj,(0,0))
 798      end if
 799 
 800      if (stress_needed==1) then
 801        ABI_MALLOC(kstr1,(npw_k))
 802        ABI_MALLOC(kstr2,(npw_k))
 803        ABI_MALLOC(kstr3,(npw_k))
 804        ABI_MALLOC(kstr4,(npw_k))
 805        ABI_MALLOC(kstr5,(npw_k))
 806        ABI_MALLOC(kstr6,(npw_k))
 807      end if
 808 
 809      ABI_MALLOC(kg_k,(3,mpw))
 810 !$OMP PARALLEL DO
 811      do ipw=1,npw_k
 812        kg_k(:,ipw)=kg(:,ipw+ikg)
 813      end do
 814 
 815      ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
 816      if (stress_needed==1) then
 817        ABI_MALLOC(ylmgr_k,(npw_k,3,mpsang*mpsang*psps%useylm))
 818      else
 819        ABI_MALLOC(ylmgr_k,(0,0,0))
 820      end if
 821      if (psps%useylm==1) then
 822 !$OMP PARALLEL DO COLLAPSE(2)
 823        do ilm=1,mpsang*mpsang
 824          do ipw=1,npw_k
 825            ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
 826          end do
 827        end do
 828        if (stress_needed==1) then
 829 !$OMP PARALLEL DO COLLAPSE(2)
 830          do ilm=1,mpsang*mpsang
 831            do ii=1,3
 832              do ipw=1,npw_k
 833                ylmgr_k(ipw,ii,ilm)=ylmgr(ipw+ikg,ii,ilm)
 834              end do
 835            end do
 836          end do
 837        end if
 838      end if
 839 
 840 !    Prepare kinetic contribution to stress tensor (Warning : the symmetry
 841 !    has not been broken, like in mkkin.f or kpg3.f . It should be, in order to be coherent).
 842      if (stress_needed==1) then
 843        ABI_MALLOC(gprimd,(3,3))
 844        gprimd=gs_hamk%gprimd
 845 !$OMP PARALLEL DO PRIVATE(fact_kin,ipw,kgc1,kgc2,kgc3,kin,xx,fsm,dfsm) &
 846 !$OMP&SHARED(ecut,ecutsm,ecutsm_inv,gs_hamk,htpisq,kg_k,kpoint,kstr1,kstr2,kstr3,kstr4,kstr5,kstr6,npw_k)
 847        do ipw=1,npw_k
 848 !        Compute Cartesian coordinates of (k+G)
 849          kgc1=gprimd(1,1)*(kpoint(1)+kg_k(1,ipw))+&
 850 &         gprimd(1,2)*(kpoint(2)+kg_k(2,ipw))+&
 851 &         gprimd(1,3)*(kpoint(3)+kg_k(3,ipw))
 852          kgc2=gprimd(2,1)*(kpoint(1)+kg_k(1,ipw))+&
 853 &         gprimd(2,2)*(kpoint(2)+kg_k(2,ipw))+&
 854 &         gprimd(2,3)*(kpoint(3)+kg_k(3,ipw))
 855          kgc3=gprimd(3,1)*(kpoint(1)+kg_k(1,ipw))+&
 856 &         gprimd(3,2)*(kpoint(2)+kg_k(2,ipw))+&
 857 &         gprimd(3,3)*(kpoint(3)+kg_k(3,ipw))
 858          kin=htpisq* ( kgc1**2 + kgc2**2 + kgc3**2 )
 859          fact_kin=1.0_dp
 860          if(kin>ecut-ecutsm)then
 861            if(kin>ecut)then
 862              fact_kin=0.0_dp
 863            else
 864 !            See the routine mkkin.f, for the smearing procedure
 865              xx=(ecut-kin)*ecutsm_inv
 866 !            This kinetic cutoff smoothing function and its xx derivatives
 867 !            were produced with Mathematica and the fortran code has been
 868 !            numerically checked against Mathematica.
 869              fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
 870              dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
 871 !            d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
 872 !            &                         (-144+45*xx))))))*fsm**3
 873              fact_kin=fsm+kin*(-ecutsm_inv)*dfsm
 874            end if
 875          end if
 876          kstr1(ipw)=fact_kin*kgc1*kgc1
 877          kstr2(ipw)=fact_kin*kgc2*kgc2
 878          kstr3(ipw)=fact_kin*kgc3*kgc3
 879          kstr4(ipw)=fact_kin*kgc3*kgc2
 880          kstr5(ipw)=fact_kin*kgc3*kgc1
 881          kstr6(ipw)=fact_kin*kgc2*kgc1
 882        end do ! ipw
 883        ABI_FREE(gprimd)
 884      end if
 885 
 886 !    Compute (k+G) vectors
 887      nkpg=3*nloalg(3)
 888      ABI_MALLOC(kpg_k,(npw_k,nkpg))
 889      if (nkpg>0) then
 890        call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
 891      end if
 892 
 893 !    Compute nonlocal form factors ffnl at all (k+G)
 894      ider=0;idir=0;dimffnl=1
 895      if (stress_needed==1) then
 896        ider=1;dimffnl=2+2*psps%useylm
 897      end if
 898      ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
 899      call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
 900 &     ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
 901 &     nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
 902      if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
 903        ider_str=1; dimffnl_str=7;idir_str=-7
 904        ABI_MALLOC(ffnl_str,(npw_k,dimffnl_str,psps%lmnmax,ntypat))
 905        call mkffnl(psps%dimekb,dimffnl_str,psps%ekb,ffnl_str,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
 906 &       ider_str,idir_str,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
 907 &       nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
 908      end if
 909 
 910 !    Load k-dependent part in the Hamiltonian datastructure
 911 !     - Compute 3D phase factors
 912 !     - Prepare various tabs in case of band-FFT parallelism
 913 !     - Load k-dependent quantities in the Hamiltonian
 914      ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk))
 915      call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
 916 &     kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_gbound=compute_gbound,compute_ph3d=.true.)
 917 
 918 !    Load band-FFT tabs (transposed k-dependent arrays)
 919      if (mpi_enreg%paral_kgb==1) then
 920        call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
 921        call prep_bandfft_tabs(gs_hamk,ikpt,mkmem,mpi_enreg)
 922        call gs_hamk%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, &
 923 &       kg_k     =my_bandfft_kpt%kg_k_gather, &
 924 &       kpg_k    =my_bandfft_kpt%kpg_k_gather, &
 925        ffnl_k   =my_bandfft_kpt%ffnl_gather, &
 926        ph3d_k   =my_bandfft_kpt%ph3d_gather,compute_gbound=compute_gbound)
 927      end if
 928 
 929 !    Setup gemm_nonlop
 930      if (gemm_nonlop_use_gemm) then
 931        gemm_nonlop_ikpt_this_proc_being_treated = my_ikpt
 932        call make_gemm_nonlop(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%lmnmax,gs_hamk%ntypat, &
 933 &            gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%istwf_k, gs_hamk%ucvol, &
 934 &            gs_hamk%ffnl_k, gs_hamk%ph3d_k, gs_hamk%kpt_k, gs_hamk%kg_k, gs_hamk%kpg_k, &
 935 &            compute_grad_strain=(stress_needed>0),compute_grad_atom=(optfor>0))
 936        end if
 937 
 938 !    Loop over (blocks of) bands; accumulate forces and/or stresses
 939 !    The following is now wrong. In sequential, nblockbd=nband_k/bandpp
 940 !    blocksize= bandpp (JB 2016/04/16)
 941 !    Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1
 942      ABI_MALLOC(lambda,(blocksize))
 943      ABI_MALLOC(occblock,(blocksize))
 944      ABI_MALLOC(weight,(blocksize))
 945      ABI_MALLOC(enlout,(nnlout*blocksize))
 946      occblock=zero;weight=zero;enlout(:)=zero
 947      if (usefock_loc) then
 948        if (fockcommon%optstr) then
 949          ABI_MALLOC(fockcommon%stress_ikpt,(6,nband_k))
 950          fockcommon%stress_ikpt=zero
 951        end if
 952      end if
 953      if ((usefock_loc).and.(psps%usepaw==1)) then
 954        if (fockcommon%optfor) then
 955          fockcommon%forces_ikpt=zero
 956        end if
 957      end if
 958 
 959      call timab(922,2,tsec)
 960 
 961      do iblock=1,nblockbd
 962 
 963        iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
 964        iband_cprj=(iblock-1)*bandpp+1
 965        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
 966 
 967 !      Select occupied bandsddk
 968        occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
 969        if( abs(maxval(occblock))>=tol8 ) then
 970          call timab(923,1,tsec)
 971          weight(:)=wtk(ikpt)*occblock(:)
 972 
 973 !        gs_hamk%ffnl_k is changed in fock_getghc, so that it is necessary to restore it when stresses are to be calculated.
 974          if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
 975            call gs_hamk%load_k(ffnl_k=ffnl)
 976          end if
 977 
 978 !        Load contribution from n,k
 979          cwavef(:,1:npw_k*my_nspinor*blocksize)=&
 980 &         cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
 981          if (psps%usepaw==1.and.usecprj_local==1) then
 982            call pawcprj_get(gs_hamk%atindx1,cwaveprj,cprj,natom,iband_cprj,ibg,ikpt,0,isppol,&
 983 &           mband_cprj,mkmem,natom,bandpp,nband_cprj_k,my_nspinor,nsppol,0,&
 984 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
 985          end if
 986 
 987          call timab(923,2,tsec)
 988          call timab(924,-1,tsec)
 989 
 990          lambda(1:blocksize)= eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
 991          if (mpi_enreg%paral_kgb/=1) then
 992            call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,mpi_enreg,blocksize,nnlout,&
 993 &           paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
 994          else
 995            call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,blocksize,&
 996 &           mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef,cwavef)
 997          end if
 998          if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
 999            call gs_hamk%load_k(ffnl_k=ffnl_str)
1000          end if
1001 
1002 !        Accumulate non-local contributions from n,k
1003          if (optfor==1) then
1004            do iblocksize=1,blocksize
1005              ibs=nnlout*(iblocksize-1)
1006              grnl(1:3*natom)=grnl(1:3*natom)+weight(iblocksize)*enlout(ibs+1+ishift:ibs+3*natom+ishift)
1007            end do
1008          end if
1009          if (stress_needed==1) then
1010            do iblocksize=1,blocksize
1011              ibs=nnlout*(iblocksize-1)
1012              npsstr(1:6)=npsstr(1:6) + weight(iblocksize)*enlout(ibs+1:ibs+6)
1013            end do
1014          end if
1015 
1016          call timab(924,2,tsec)
1017 
1018 !        Accumulate stress tensor kinetic contributions
1019          if (stress_needed==1) then
1020            call timab(925,1,tsec)
1021            do iblocksize=1,blocksize
1022              call meanvalue_g(ar,kstr1,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1023 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1024 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1025              kinstr(1)=kinstr(1)+weight(iblocksize)*ar
1026              call meanvalue_g(ar,kstr2,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1027 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1028 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1029              kinstr(2)=kinstr(2)+weight(iblocksize)*ar
1030              call meanvalue_g(ar,kstr3,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1031 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1032 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1033              kinstr(3)=kinstr(3)+weight(iblocksize)*ar
1034              call meanvalue_g(ar,kstr4,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1035 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1036 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1037              kinstr(4)=kinstr(4)+weight(iblocksize)*ar
1038              call meanvalue_g(ar,kstr5,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1039 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1040 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1041              kinstr(5)=kinstr(5)+weight(iblocksize)*ar
1042              call meanvalue_g(ar,kstr6,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
1043 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
1044 &             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
1045              kinstr(6)=kinstr(6)+weight(iblocksize)*ar
1046            end do
1047            call timab(925,2,tsec)
1048          end if
1049 
1050 !        Accumulate stress tensor and forces for the Fock part
1051          if (usefock_loc) then
1052            if(fockcommon%optstr.or.fockcommon%optfor) then
1053              call timab(926,1,tsec)
1054              if (mpi_enreg%paral_kgb==1) then
1055                msg='forsrtnps: Paral_kgb is not yet implemented for fock stresses'
1056                ABI_BUG(msg)
1057              end if
1058              ndat=mpi_enreg%bandpp
1059              if (gs_hamk%usepaw==0) cwaveprj_idat => cwaveprj
1060              ABI_MALLOC(ghc_dum,(0,0))
1061              do iblocksize=1,blocksize
1062                fockcommon%ieigen=(iblock-1)*blocksize+iblocksize
1063                fockcommon%iband=(iblock-1)*blocksize+iblocksize
1064                if (gs_hamk%usepaw==1) then
1065                  cwaveprj_idat => cwaveprj(:,(iblocksize-1)*my_nspinor+1:iblocksize*my_nspinor)
1066                end if
1067                call fock_getghc(cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),cwaveprj_idat,&
1068 &               ghc_dum,gs_hamk,mpi_enreg)
1069                if (fockcommon%optstr) then
1070                  fockcommon%stress(:)=fockcommon%stress(:)+weight(iblocksize)*fockcommon%stress_ikpt(:,fockcommon%ieigen)
1071                end if
1072                if (fockcommon%optfor) then
1073                  fockcommon%forces(:,:)=fockcommon%forces(:,:)+weight(iblocksize)*fockcommon%forces_ikpt(:,:,fockcommon%ieigen)
1074                end if
1075              end do
1076              ABI_FREE(ghc_dum)
1077              call timab(926,2,tsec)
1078            end if
1079          end if ! usefock_loc
1080        end if
1081 
1082      end do ! End of loop on block of bands
1083 
1084      call timab(927,1,tsec)
1085 
1086 !    Restore the bandfft tabs
1087      if (mpi_enreg%paral_kgb==1) then
1088        call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
1089      end if
1090 
1091 !    Increment indexes
1092      bdtot_index=bdtot_index+nband_k
1093      if (mkmem/=0) then
1094        ibg=ibg+my_nspinor*nband_cprj_k
1095        icg=icg+npw_k*my_nspinor*nband_k
1096        ikg=ikg+npw_k
1097      end if
1098 
1099      if (usefock_loc) then
1100        if (fockcommon%optstr) then
1101          ABI_FREE(fockcommon%stress_ikpt)
1102        end if
1103      end if
1104 
1105      if (psps%usepaw==1) then
1106        call pawcprj_free(cwaveprj)
1107      end if
1108      ABI_FREE(cwaveprj)
1109      ABI_FREE(cwavef)
1110 
1111      ABI_FREE(lambda)
1112      ABI_FREE(occblock)
1113      ABI_FREE(weight)
1114      ABI_FREE(enlout)
1115      ABI_FREE(ffnl)
1116      ABI_FREE(kg_k)
1117      ABI_FREE(kpg_k)
1118      ABI_FREE(ylm_k)
1119      ABI_FREE(ylmgr_k)
1120      ABI_FREE(ph3d)
1121      if (stress_needed==1) then
1122        ABI_FREE(kstr1)
1123        ABI_FREE(kstr2)
1124        ABI_FREE(kstr3)
1125        ABI_FREE(kstr4)
1126        ABI_FREE(kstr5)
1127        ABI_FREE(kstr6)
1128      end if
1129      if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then
1130        ABI_FREE(ffnl_str)
1131      end if
1132 
1133      call timab(927,2,tsec)
1134 
1135    end do ! End k point loop
1136  end do ! End loop over spins
1137 
1138  call timab(928,1,tsec)
1139 
1140 !Stress is equal to dE/d_strain * (1/ucvol)
1141  npsstr(:)=npsstr(:)/gs_hamk%ucvol
1142 
1143 !Parallel case: accumulate (n,k) contributions
1144  if (xmpi_paral==1) then
1145 !  Forces
1146    if (optfor==1) then
1147      call timab(65,1,tsec)
1148      call xmpi_sum(grnl,spaceComm,ierr)
1149      call timab(65,2,tsec)
1150      if ((usefock_loc).and.(psps%usepaw==1)) then
1151        call xmpi_sum(fockcommon%forces,spaceComm,ierr)
1152      end if
1153    end if
1154 !  Stresses
1155    if (stress_needed==1) then
1156      call timab(65,1,tsec)
1157      call xmpi_sum(kinstr,spaceComm,ierr)
1158      call xmpi_sum(npsstr,spaceComm,ierr)
1159      if (usefock_loc) then
1160        if (fockcommon%optstr) then
1161          call xmpi_sum(fockcommon%stress,spaceComm,ierr)
1162        end if
1163      end if
1164      call timab(65,2,tsec)
1165    end if
1166  end if
1167 
1168 !Do final normalizations and symmetrizations of stress tensor contributions
1169  if (stress_needed==1) then
1170    renorm_factor=-(two_pi**2)/effmass_free/gs_hamk%ucvol
1171    kinstr(:)=kinstr(:)*renorm_factor
1172    if (nsym>1) then
1173      call stresssym(gs_hamk%gprimd,nsym,kinstr,symrec)
1174      call stresssym(gs_hamk%gprimd,nsym,npsstr,symrec)
1175      if (usefock_loc) then
1176        if (fockcommon%optstr) then
1177          call stresssym(gs_hamk%gprimd,nsym,fockcommon%stress,symrec)
1178        end if
1179      end if
1180    end if
1181  end if
1182 
1183 !Need to reorder cprj=<p_lmn|Cnk> (from atom-sorted to unsorted)
1184  if (psps%usepaw==1.and.usecprj_local==1) then
1185    call pawcprj_reorder(cprj,gs_hamk%atindx1)
1186  end if
1187 
1188 !Deallocate temporary space
1189  call gs_hamk%free()
1190  if (usefock_loc) then
1191    fockcommon%use_ACE=use_ACE_old
1192  end if
1193  call timab(928,2,tsec) ; call timab(920,-2,tsec)
1194 
1195 end subroutine forstrnps

ABINIT/m_forstr [ Modules ]

[ Top ] [ Modules ]

NAME

  m_forstr

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, AF, AR, MB, MT)
  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_forstr
23 
24  use defs_basis
25  use defs_wvltypes
26  use m_abicore
27  use m_efield
28  use m_errors
29  use m_xmpi
30  use m_fock
31  use m_cgtools
32  use m_xcdata
33  use m_dtset
34  use m_extfpmd
35 
36  use defs_datatypes,     only : pseudopotential_type
37  use defs_abitypes,      only : MPI_type
38  use m_time,             only : timab
39  use m_geometry,         only : xred2xcart, metric, stresssym
40  use m_energies,         only : energies_type
41  use m_pawang,           only : pawang_type
42  use m_pawrad,           only : pawrad_type
43  use m_pawtab,           only : pawtab_type
44  use m_paw_ij,           only : paw_ij_type
45  use m_pawfgrtab,        only : pawfgrtab_type
46  use m_pawrhoij,         only : pawrhoij_type
47  use m_pawfgr,           only : pawfgr_type
48  use m_pawcprj,          only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get, pawcprj_reorder, pawcprj_getdim
49  use m_paw_dfpt,         only : pawgrnl
50  use libxc_functionals,  only : libxc_functionals_is_hybrid
51  use m_stress,           only : stress
52  use m_forces,           only : forces
53  use m_initylmg,         only : initylmg
54  use m_xchybrid,         only : xchybrid_ncpp_cc
55  use m_kg,               only : mkkpg
56  use m_hamiltonian,      only : init_hamiltonian, gs_hamiltonian_type !,K_H_KPRIME
57  use m_electronpositron, only : electronpositron_type, electronpositron_calctype
58  use m_bandfft_kpt,      only : bandfft_kpt, bandfft_kpt_type, prep_bandfft_tabs, &
59 &                               bandfft_kpt_savetabs, bandfft_kpt_restoretabs
60  use m_spacepar,         only : meanvalue_g, hartre
61  use m_mkffnl,           only : mkffnl
62  use m_mpinfo,           only : proc_distrb_cycle
63  use m_nonlop,           only : nonlop
64  use m_gemm_nonlop,      only : make_gemm_nonlop,gemm_nonlop_use_gemm, &
65 &                               gemm_nonlop_ikpt_this_proc_being_treated
66  use m_fock_getghc,      only : fock_getghc
67  use m_prep_kgb,         only : prep_nonlop
68  use m_paw_nhat,         only : pawmknhat
69  use m_rhotoxc,          only : rhotoxc
70  use m_dfpt_mkvxc,       only : dfpt_mkvxc, dfpt_mkvxc_noncoll
71  use m_cgprj,            only : ctocprj
72  use m_psolver,          only : psolver_hartree
73  use m_wvl_psi,          only : wvl_nl_gradient
74  use m_fft,              only : fourdp
75  use, intrinsic :: iso_c_binding,      only : c_loc,c_f_pointer
76 
77  implicit none
78 
79  private

ABINIT/nres2vres [ Functions ]

[ Top ] [ Functions ]

NAME

 nres2vres

FUNCTION

 Convert a density residual into a potential residual
 using a first order formula:
     V^res(r)=dV/dn.n^res(r)
             =V_hartree(n^res)(r) + Kxc.n^res(r)

INPUTS

 dtset <type(dataset_type)>=all input variables in this dataset
  | icoulomb=0 periodic treatment of Hartree potential, 1 use of Poisson solver
  | natom= number of atoms in cell
  | nspden=number of spin-density components
  | ntypat=number of atom types
  | typat(natom)=type (integer) for each atom
 gsqcut=cutoff value on G**2 for sphere inside fft box
 izero=if 1, unbalanced components of Vhartree(g) are set to zero
 kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
 mpi_enreg=information about MPI parallelization
 my_natom=number of atoms treated by current processor
 nfft=(effective) number of FFT grid points (for this processor)
 ngfft(18)=contain all needed information about 3D FFT
 nhat(nfft,nspden*usepaw)= -PAW only- compensation density
 nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description
 nresid(nfft,nspden)= the input density residual
 n3xccc=dimension of the xccc3d array (0 or nfft).
 optnc=option for non-collinear magnetism (nspden=4):
       1: the whole 2x2 Vres matrix is computed
       2: only Vres^{11} and Vres^{22} are computed
 optxc=0 if LDA part of XC kernel has only to be taken into account (even for GGA)
       1 if XC kernel has to be fully taken into
      -1 if XC kernel does not have to be taken into account
 pawang <type(pawang_type)>=paw angular mesh and related data
 pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
 pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
 pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
 rhor(nfft,nspden)=electron density in real space
                   (used only if Kxc was not computed before)
 rprimd(3,3)=dimensional primitive translation vectors (bohr)
 usepaw= 0 for non paw calculation; =1 for paw calculation
 xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
 xred(3,natom)=reduced dimensionless atomic coordinates

 === optional inputs ===
 vxc(cplex*nfft,nspden)=XC GS potential

OUTPUT

 vresid(nfft,nspden)= the output potential residual

SOURCE

1251 subroutine nres2vres(dtset,gsqcut,izero,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,&
1252 &                 nkxc,nresid,n3xccc,optnc,optxc,pawang,pawfgrtab,pawrhoij,pawtab,&
1253 &                 rhor,rprimd,usepaw,vresid,xccc3d,xred,vxc)
1254 
1255 !Arguments ------------------------------------
1256 !scalars
1257  integer,intent(in) :: izero,my_natom,n3xccc,nfft,nkxc,optnc,optxc,usepaw
1258  real(dp),intent(in) :: gsqcut
1259  type(MPI_type),intent(in) :: mpi_enreg
1260  type(dataset_type),intent(in) :: dtset
1261  type(pawang_type),intent(in) :: pawang
1262 !arrays
1263  integer,intent(in) :: ngfft(18)
1264  real(dp),intent(in) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden)
1265  real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3),xccc3d(n3xccc),xred(3,dtset%natom)
1266  real(dp),intent(inout) :: nhat(nfft,dtset%nspden*usepaw)
1267  real(dp),intent(out) :: vresid(nfft,dtset%nspden)
1268  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*usepaw)
1269  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*usepaw)
1270  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw)
1271  real(dp),intent(in) :: vxc(nfft,dtset%nspden) !FR TODO:cplex?
1272 
1273 !Local variables-------------------------------
1274 !scalars
1275  integer :: cplex,ider,idir,ipert,ispden,nhatgrdim,nkxc_cur,option,me,nproc,comm,usexcnhat
1276  logical :: has_nkxc_gga,non_magnetic_xc
1277  real(dp) :: dum,energy,m_norm_min,ucvol,vxcavg
1278  character(len=500) :: message
1279  type(xcdata_type) :: xcdata
1280 !arrays
1281  integer :: nk3xc
1282  real(dp) :: dummy6(6),gmet(3,3),gprimd(3,3),qq(3),rmet(3,3)
1283  real(dp),allocatable :: dummy(:),kxc_cur(:,:),nhatgr(:,:,:)
1284  real(dp),allocatable :: nresg(:,:),rhor0(:,:),vhres(:)
1285 
1286 ! *************************************************************************
1287 
1288 !Compatibility tests:
1289  has_nkxc_gga=(nkxc==7.or.nkxc==19)
1290 
1291  if(optxc<-1.or.optxc>1)then
1292    write(message,'(a,i0)')' Wrong value for optxc ',optxc
1293    ABI_BUG(message)
1294  end if
1295 
1296  if((optnc/=1.and.optnc/=2).or.(dtset%nspden/=4.and.optnc/=1))then
1297    write(message,'(a,i0)')' Wrong value for optnc ',optnc
1298    ABI_BUG(message)
1299  end if
1300 
1301  if(dtset%icoulomb==1.and.optxc/=-1)then
1302    write(message,'(a)')' This routine is not compatible with icoulomb==1 and optxc/=-1 !'
1303    ABI_BUG(message)
1304  end if
1305 
1306  if(dtset%nspden==4.and.dtset%xclevel==2.and.optxc==1.and.(.not.has_nkxc_gga))then
1307    ABI_ERROR(' Wrong values for optxc and nkxc !')
1308  end if
1309 
1310  qq=zero
1311  nkxc_cur=0
1312  m_norm_min=EPSILON(0.0_dp)**2
1313  usexcnhat=0;if (usepaw==1) usexcnhat=maxval(pawtab(1:dtset%ntypat)%usexcnhat)
1314  non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
1315  if (dtset%xclevel==1.or.optxc==0) nkxc_cur= 2*min(dtset%nspden,2)-1 ! LDA: nkxc=1,3
1316  if (dtset%xclevel==2.and.optxc==1)nkxc_cur=12*min(dtset%nspden,2)-5 ! GGA: nkxc=7,19
1317  ABI_MALLOC(vhres,(nfft))
1318 
1319 !Compute different geometric tensor, as well as ucvol, from rprimd
1320  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1321 
1322 !Compute density residual in reciprocal space
1323  if (dtset%icoulomb==0) then
1324    ABI_MALLOC(nresg,(2,nfft))
1325    ABI_MALLOC(dummy,(nfft))
1326    dummy(:)=nresid(:,1)
1327    call fourdp(1,nresg,dummy,-1,mpi_enreg,nfft,1,ngfft,0)
1328    ABI_FREE(dummy)
1329  end if
1330 
1331 !For GGA, has to recompute gradients of nhat
1332  nhatgrdim=0
1333  if ((nkxc==nkxc_cur.and.has_nkxc_gga).or.(optxc==-1.and.has_nkxc_gga).or.&
1334 & (optxc/=-1.and.nkxc/=nkxc_cur)) then
1335    if (usepaw==1.and.dtset%xclevel==2.and.usexcnhat>0.and.dtset%pawnhatxc>0) then
1336      nhatgrdim=1
1337      ABI_MALLOC(nhatgr,(nfft,dtset%nspden,3))
1338      ider=1;cplex=1;ipert=0;idir=0
1339      call pawmknhat(dum,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
1340 &     nfft,ngfft,nhatgrdim,dtset%nspden,dtset%ntypat,pawang,pawfgrtab,&
1341 &     nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qq,rprimd,ucvol,dtset%usewvl,xred,&
1342 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1343 &     comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
1344 &     distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
1345    else
1346      ABI_MALLOC(nhatgr,(0,0,0))
1347    end if
1348  else
1349    ABI_MALLOC(nhatgr,(0,0,0))
1350  end if
1351 
1352  ABI_MALLOC(dummy,(0))
1353 !First case: Kxc has already been computed
1354 !-----------------------------------------
1355  if (nkxc==nkxc_cur.or.optxc==-1) then
1356 
1357 !  Compute VH(n^res)(r)
1358    if (dtset%icoulomb == 0) then
1359      call hartre(1,gsqcut,dtset%icutcoul,izero,mpi_enreg,nfft,ngfft,&
1360                  &dtset%nkpt,dtset%rcut,nresg,rprimd,dtset%vcutgeo,vhres)
1361    else
1362      comm=mpi_enreg%comm_cell
1363      nproc=xmpi_comm_size(comm)
1364      me=xmpi_comm_rank(comm)
1365      call psolver_hartree(energy, (/ rprimd(1,1) / dtset%ngfft(1), &
1366 &     rprimd(2,2) / dtset%ngfft(2), rprimd(3,3) / dtset%ngfft(3) /), dtset%icoulomb, &
1367 &     me, comm, dtset%nfft, dtset%ngfft(1:3), nproc, dtset%nscforder, dtset%nspden, &
1368 &     nresid(:,1), vhres, dtset%usewvl)
1369    end if
1370 
1371 !  Compute Kxc(r).n^res(r)
1372    if (optxc/=-1) then
1373 
1374 !    Collinear magnetism or non-polarized
1375      if (dtset%nspden/=4) then
1376        !Note: imposing usexcnhat=1 avoid nhat to be substracted
1377        call dfpt_mkvxc(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,&
1378 &       nkxc,non_magnetic_xc,dtset%nspden,0,2,qq,nresid,rprimd,1,vresid,dummy)
1379      else
1380 !FR    call routine for Non-collinear magnetism
1381        ABI_MALLOC(rhor0,(nfft,dtset%nspden))
1382        rhor0(:,:)=rhor(:,:)-nresid(:,:)
1383        !Note: imposing usexcnhat=1 avoid nhat to be substracted
1384        call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,&
1385 &       nkxc,non_magnetic_xc,dtset%nspden,0,2,2,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d)
1386        ABI_FREE(rhor0)
1387      end if
1388 
1389    else
1390      vresid=zero
1391    end if
1392 
1393  end if
1394 
1395 !2nd case: Kxc has to be computed
1396 !--------------------------------
1397  if (nkxc/=nkxc_cur.and.optxc/=-1) then
1398 
1399 !  Has to use the "initial" density to compute Kxc
1400    ABI_MALLOC(rhor0,(nfft,dtset%nspden))
1401    rhor0(:,:)=rhor(:,:)-nresid(:,:)
1402 
1403 !  Compute VH(n^res) and XC kernel (Kxc) together
1404    ABI_MALLOC(kxc_cur,(nfft,nkxc_cur))
1405 
1406    option=2;if (dtset%xclevel==2.and.optxc==0) option=12
1407 
1408    call hartre(1,gsqcut,dtset%icutcoul,izero,mpi_enreg,nfft,ngfft,&
1409                &dtset%nkpt,dtset%rcut,nresg,rprimd,dtset%vcutgeo,vhres)
1410    call xcdata_init(xcdata,dtset=dtset)
1411 
1412 !  To be adjusted for the call to rhotoxc
1413    nk3xc=1
1414    call rhotoxc(energy,kxc_cur,mpi_enreg,nfft,ngfft,&
1415 &   nhat,usepaw,nhatgr,nhatgrdim,nkxc_cur,nk3xc,non_magnetic_xc,n3xccc,option,&
1416 &   rhor0,rprimd,dummy6,usexcnhat,vresid,vxcavg,xccc3d,xcdata,vhartr=vhres)  !vresid=work space
1417    if (dtset%nspden/=4)  then
1418      ABI_FREE(rhor0)
1419    end if
1420 
1421 !  Compute Kxc(r).n^res(r)
1422 
1423    if (dtset%nspden/=4) then
1424 !    Collinear magnetism or non-polarized
1425      call dfpt_mkvxc(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,&
1426 &     nkxc_cur,non_magnetic_xc,dtset%nspden,0,2,qq,nresid,rprimd,1,vresid,dummy)
1427    else
1428 !    Non-collinear magnetism
1429      ABI_MALLOC(rhor0,(nfft,dtset%nspden))
1430      rhor0(:,:)=rhor(:,:)-nresid(:,:)
1431      call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,&
1432 &     nkxc,non_magnetic_xc,dtset%nspden,0,2,2,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d)
1433      ABI_FREE(rhor0)
1434    end if
1435 
1436    ABI_FREE(kxc_cur)
1437  end if
1438 
1439  !if (nhatgrdim>0)  then
1440  ABI_FREE(nhatgr)
1441  !end if
1442 
1443 !Assemble potential residual: V^res(r)=VH(n^res)(r) + Kxc(r).n^res(r)
1444 !--------------------------------------------------------------------
1445  do ispden=1,dtset%nspden/optnc
1446    vresid(:,ispden)=vresid(:,ispden)+vhres(:)
1447  end do
1448 
1449  if (dtset%icoulomb==0)  then
1450    ABI_FREE(nresg)
1451  end if
1452  ABI_FREE(vhres)
1453  ABI_FREE(dummy)
1454 
1455 end subroutine nres2vres