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)
  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)
  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
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
  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 fred, 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)
   fred(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.

PARENTS

      afterscfloop,setup_positron

CHILDREN

      ctocprj,forces,forstrnps,initylmg,metric,nres2vres,pawcprj_alloc
      pawcprj_free,pawcprj_getdim,pawgrnl,stress,timab,wvl_nl_gradient
      xchybrid_ncpp_cc,xred2xcart

SOURCE

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

PARENTS

      forstr

CHILDREN

      bandfft_kpt_restoretabs,bandfft_kpt_savetabs,destroy_hamiltonian
      fock_getghc,init_hamiltonian,load_k_hamiltonian,load_spin_hamiltonian
      meanvalue_g,mkffnl,mkkpg,nonlop,pawcprj_alloc,pawcprj_free,pawcprj_get
      pawcprj_reorder,prep_bandfft_tabs,prep_nonlop,stresssym,timab,xmpi_sum

SOURCE

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

ABINIT/m_forstr [ Modules ]

[ Top ] [ Modules ]

NAME

  m_forstr

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_forstr
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_abicore
34  use m_efield
35  use m_errors
36  use m_xmpi
37  use m_fock
38  use m_cgtools
39  use m_xcdata
40 
41  use m_time,             only : timab
42  use m_geometry,         only : xred2xcart, metric, stresssym
43  use m_energies,         only : energies_type
44  use m_pawang,           only : pawang_type
45  use m_pawrad,           only : pawrad_type
46  use m_pawtab,           only : pawtab_type
47  use m_paw_ij,           only : paw_ij_type
48  use m_pawfgrtab,        only : pawfgrtab_type
49  use m_pawrhoij,         only : pawrhoij_type
50  use m_pawfgr,           only : pawfgr_type
51  use m_pawcprj,          only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get, pawcprj_reorder, pawcprj_getdim
52  use m_paw_dfpt,         only : pawgrnl
53  use libxc_functionals,  only : libxc_functionals_is_hybrid
54  use m_stress,           only : stress
55  use m_forces,           only : forces
56  use m_initylmg,         only : initylmg
57  use m_xchybrid,         only : xchybrid_ncpp_cc
58  use m_kg,               only : mkkpg
59  use m_hamiltonian,      only : init_hamiltonian, destroy_hamiltonian, load_spin_hamiltonian,&
60 &                               load_k_hamiltonian, gs_hamiltonian_type, load_kprime_hamiltonian!,K_H_KPRIME
61  use m_electronpositron, only : electronpositron_type, electronpositron_calctype
62  use m_bandfft_kpt,      only : bandfft_kpt, bandfft_kpt_type, prep_bandfft_tabs, &
63 &                               bandfft_kpt_savetabs, bandfft_kpt_restoretabs
64  use m_spacepar,         only : meanvalue_g, hartre
65  use m_mkffnl,           only : mkffnl
66  use m_mpinfo,           only : proc_distrb_cycle
67  use m_nonlop,           only : nonlop
68  use m_fock_getghc,      only : fock_getghc
69  use m_prep_kgb,         only : prep_nonlop
70  use m_paw_nhat,         only : pawmknhat
71  use m_rhotoxc,          only : rhotoxc
72  use m_dfpt_mkvxc,       only : dfpt_mkvxc, dfpt_mkvxc_noncoll
73  use m_cgprj,            only : ctocprj
74  use m_psolver,          only : psolver_hartree
75  use m_wvl_psi,          only : wvl_nl_gradient
76  use m_fft,              only : fourdp
77 
78  implicit none
79 
80  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

PARENTS

      etotfor,forstr

CHILDREN

      dfpt_mkvxc,dfpt_mkvxc_noncoll,fourdp,hartre,metric,pawmknhat
      psolver_hartree,rhotoxc,xcdata_init

SOURCE

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