TABLE OF CONTENTS


ABINIT/constrf [ Functions ]

[ Top ] [ Functions ]

NAME

 constrf

FUNCTION

 Computes projected forces, fpcart, which satisfy a set of
 constraint equations of the form
  Sum[mu,iatom]: wtatcon(mu,iatom,iconeq)*fpcart(mu,iatom) = 0 (iconeq=1,nconeq).
 These projected forces are returned in fcart and thus replace
 the original forces.

INPUTS

  iatfix(3,natom)=1 for frozen atom along each direction, 0 for unfrozen
  natom=number of atoms in cell
  nconeq=number of atomic constraint equations
  prtvol=control print volume and debugging
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  wtatcon(3,natom,nconeq)=weights for atomic constraints
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  diffor=maximum absolute change in component of projected fcart between present
          and previous SCF cycle
  gred(3,natom)=grads of Etot wrt reduced coordinates (hartree)
  maxfor=maximum absolute value of fcart

SIDE EFFECTS

  fcart(3,natom)=cartesian forces (hartree/bohr) on input, projected forces on output
  forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)

TODO

SOURCE

1494 subroutine constrf(diffor,fcart,forold,gred,iatfix,ionmov,maxfor,natom,&
1495 & nconeq,prtvol,rprimd,wtatcon,xred)
1496 
1497  use m_linalg_interfaces
1498 
1499 !Arguments ------------------------------------
1500 !scalars
1501  integer,intent(in) :: ionmov,natom,nconeq,prtvol
1502  real(dp),intent(out) :: diffor,maxfor
1503 !arrays
1504  integer,intent(in) :: iatfix(3,natom)
1505  real(dp),intent(in) :: rprimd(3,3),wtatcon(3,natom,nconeq)
1506  real(dp),intent(inout) :: fcart(3,natom),forold(3,natom),xred(3,natom)
1507  real(dp),intent(inout) :: gred(3,natom) !vz_i
1508 
1509 !Local variables -------------------------
1510 !scalars
1511  integer :: iatom,iconeq,iconeq1,iconeq2,index,info,mu,prtvel
1512  character(len=500) :: message
1513 !arrays
1514  real(dp),allocatable :: fcartvec(:),fpcart(:,:),fvector(:),vel_dummy(:,:)
1515  real(dp),allocatable :: wmatrix(:,:),wtatconvec(:,:),wtcoeffs(:),xcart(:,:)
1516 
1517 !************************************************************************
1518 
1519 !Allocate temporary variables
1520  ABI_MALLOC(fpcart,(3,natom))
1521  ABI_MALLOC(fcartvec,(3*natom))
1522  ABI_MALLOC(fvector,(nconeq))
1523  ABI_MALLOC(vel_dummy,(3,natom))
1524  ABI_MALLOC(wmatrix,(nconeq,nconeq))
1525  ABI_MALLOC(wtatconvec,(3*natom,nconeq))
1526  ABI_MALLOC(wtcoeffs,(nconeq))
1527  ABI_MALLOC(xcart,(3,natom))
1528 
1529 !If prtvol>10, output coordinates and forces prior to projecting
1530  if(prtvol>=10)then
1531    write(message,'(a)')' constrf - coordinates and forces prior to constraint projections:'
1532    call wrtout(std_out,message,'COLL')
1533    call xred2xcart(natom,rprimd,xcart,xred)
1534    prtvel=0
1535    call prtxvf(fcart,gred,iatfix,06,natom,prtvel,vel_dummy,xcart,xred)
1536  end if
1537 
1538 !Transfer fcart and wtatcon to flat vectors
1539  index=0
1540  do iatom=1,natom
1541    do mu=1,3
1542      index=index+1
1543      fcartvec(index)=fcart(mu,iatom)
1544      wtatconvec(index,:)=wtatcon(mu,iatom,:)
1545    end do
1546  end do
1547 
1548 !Compute a matrix (wmatrix) and vector (fvector) such that solving
1549 !the linear equations wmatrix*wcoeffs=fvector gives the coefficients
1550 !of wtatcon (wcoeffs) needed to compute the projected forces
1551  do iconeq2=1,nconeq
1552    fvector(iconeq2)=ddot(3*natom,fcartvec,1,wtatconvec(1,iconeq2),1)
1553    do iconeq1=1,nconeq
1554      wmatrix(iconeq1,iconeq2)=ddot(3*natom,wtatconvec(1,iconeq1),1,wtatconvec(1,iconeq2),1)
1555    end do
1556  end do
1557 
1558 !Solve the system of linear equations, wmatrix*wcoeffs=fvector
1559  call dposv('U',nconeq,1,wmatrix,nconeq,fvector,nconeq,info)
1560 
1561  if (info/=0) then
1562    write(message, '(a,a,a,a,a)' )&
1563 &   'Constraint matrix is not positive definite,',ch10,&
1564 &   'probably because constraints are linearly dependent.',ch10,&
1565 &   'Action: Check for linear dependence of constraints.'
1566    ABI_ERROR(message)
1567  end if
1568 
1569 !The solution vector is returned in fvector, so copy it to a more sensible location
1570  wtcoeffs(:)=fvector(:)
1571 
1572 !Compute the projected forces, which now satisfy all the constraint equations
1573  fpcart(:,:)=fcart(:,:)
1574  do iconeq=1,nconeq
1575    fpcart(:,:)=fpcart(:,:)-wtcoeffs(iconeq)*wtatcon(:,:,iconeq)
1576  end do
1577 
1578 !Reconvert constrained forces back from fpcart to gred
1579  do iatom=1,natom
1580    do mu=1,3
1581      gred(mu,iatom)= - (rprimd(1,mu)*fpcart(1,iatom)+&
1582 &     rprimd(2,mu)*fpcart(2,iatom)+&
1583 &     rprimd(3,mu)*fpcart(3,iatom))
1584    end do
1585  end do
1586 
1587 !If prtvol>=10, output coordinates and forces after projecting
1588  if(prtvol>=10)then
1589    write(message,'(a)')' constrf - coordinates and forces after constraint projections:'
1590    call wrtout(std_out,message,'COLL')
1591    prtvel=0
1592    call prtxvf(fpcart,gred,iatfix,06,natom,prtvel,vel_dummy,xcart,xred)
1593  end if
1594 
1595 !Copy the constrained forces, fpcart, back to fcart
1596  fcart(:,:)=fpcart(:,:)
1597 
1598 !Compute maximal force and maximal difference of the projected forces,
1599 !overriding the values already computed in forces
1600  maxfor=0.0_dp
1601  diffor=0.0_dp
1602  do iatom=1,natom
1603    do mu=1,3
1604      if (iatfix(mu,iatom) /= 1) then
1605        maxfor=max(maxfor,abs(fcart(mu,iatom)))
1606        diffor=max(diffor,abs(fcart(mu,iatom)-forold(mu,iatom)))
1607      else if (ionmov==4 .or. ionmov==5) then
1608 !      Make the force vanish on fixed atoms when ionmov=4 or 5
1609 !      This is because fixing of atom cannot be imposed at the
1610 !      level of a routine similar to brdmin or moldyn for these options.
1611        fcart(mu,iatom)=0.0_dp
1612      end if
1613    end do
1614  end do
1615  forold(:,:)=fcart(:,:)
1616 
1617 !Dellocate temporary variables
1618  ABI_FREE(fpcart)
1619  ABI_FREE(fcartvec)
1620  ABI_FREE(fvector)
1621  ABI_FREE(vel_dummy)
1622  ABI_FREE(wmatrix)
1623  ABI_FREE(wtatconvec)
1624  ABI_FREE(wtcoeffs)
1625  ABI_FREE(xcart)
1626 
1627 end subroutine constrf

ABINIT/forces [ Functions ]

[ Top ] [ Functions ]

NAME

 forces

FUNCTION

 Assemble gradients of various total energy terms with respect
 to reduced coordinates, including possible symmetrization,
 in order to produce forces.

     fcart(i,iat) = d(Etot)/(d(r(i,iat)))

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  dtefield <type(efield_type)> = variables related to Berry phase
  dtset <type(dataset_type)>=all input variables in this dataset
 berryopt    =  4/14: electric field is on -> add the contribution of the
                      -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j  terms to the total energy
     = 6/16, or 7/17: electric displacement field is on  -> add the contribution of the
                      Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j  terms to the total energy
   | efield = cartesian coordinates of the electric field in atomic units
   | dfield = cartesian coordinates of the electric displacement field in atomic units
   | iatfix(3,natom)=1 for frozen atom along specified direction, 0 for unfrozen
   | ionmov=governs the movement of atoms (see help file)
   | densfor_pred=governs the mixed electronic-atomic part of the preconditioner
   | natom=number of atoms in cell
   | nconeq=number of atomic constraint equations
   | nspden=number of spin-density components
   | nsym=number of symmetries in space group
   | prtvol=integer controlling volume of printed output
   | typat(natom)=type integer for each atom in cell
   | wtatcon(3,natom,nconeq)=weights for atomic constraints
  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)
  grnl(3*natom)=gradients of Etot due to nonlocal contributions
  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)=indirect indexing array for atom labels
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  n3xccc=dimension of the xccc3d array (0 or nfft).
  nattyp(ntypat)=number of atoms of each type
  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
  ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
  ntypat=number of types of atoms
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) array
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhog(2,nfft)=Fourier transform of charge density (bohr^-3)
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
  usefock=1 if fock operator is used; 0 otherwise.
  usekden= 1 is kinetic energy density has to be computed, 0 otherwise
  vresid(nfft,nspden)=potential residual (if non-collinear magn., only trace of it)
  vxc(nfft,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)
  xred(3,natom)=reduced dimensionless atomic coordinates
  xred_old(3,natom)=previous reduced dimensionless atomic coordinates

OUTPUT

  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
  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 d(enl)/d(xred)

SIDE EFFECTS

  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
  fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
    Note : unlike gred, this array has been corrected by enforcing
    the translational symmetry, namely that the sum of force
    on all atoms is zero.

NOTES

 * Symmetrization of gradients with respect to reduced
   coordinates xred is conducted according to the expression
   [d(e)/d(t(n,a))]_symmetrized = (1/Nsym) Sum(S) symrec(n,m,S)*
                [d(e)/d(t(m,b))]_unsymmetrized
   where t(m,b)= (symrel^-1)(m,n)*(t(n,a)-tnons(n)) and tnons
   is a possible nonsymmorphic translation.  The label "b" here
   refers to the atom which gets rotated into "a" under symmetry "S".
   symrel is the symmetry matrix in real space, which is the inverse
   transpose of symrec.  symrec is the symmetry matrix in reciprocal
   space.  sym_cartesian = R * symrel * R^-1 = G * symrec * G^-1
   where the columns of R and G are the dimensional primitive translations
   in real and reciprocal space respectively.
 * Note the use of "symrec" in the symmetrization expression above.

SOURCE

165 subroutine forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,&
166 &                  forold,gred,grchempottn,grcondft,gresid,grewtn,&
167 &                  grhf,grnl,grvdw,grxc,gsqcut,indsym,&
168 &                  maxfor,mgfft,mpi_enreg,n1xccc,n3xccc,&
169 &                  nattyp,nfft,ngfft,ngrvdw,ntypat,&
170 &                  pawrad,pawtab,ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,usefock,&
171 &                  vresid,vxc,vxctau,wvl,wvl_den,xred,&
172 &                  electronpositron) ! optional argument
173 
174 !Arguments ------------------------------------
175 !scalars
176  integer,intent(in) :: mgfft,n1xccc,n3xccc,nfft,ngrvdw,ntypat,usefock
177  real(dp),intent(in) :: gsqcut
178  real(dp),intent(out) :: diffor,maxfor
179  type(MPI_type),intent(in) :: mpi_enreg
180  type(efield_type),intent(in) :: dtefield
181  type(dataset_type),intent(in) :: dtset
182  type(electronpositron_type),pointer,optional :: electronpositron
183  type(pseudopotential_type),intent(in) :: psps
184  type(wvl_internal_type), intent(in) :: wvl
185  type(wvl_denspot_type), intent(inout) :: wvl_den
186  type(fock_type),pointer, intent(inout) :: fock
187 !arrays
188  integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
189  integer,intent(in) :: nattyp(ntypat),ngfft(18),symrec(3,3,dtset%nsym)
190  real(dp),intent(in) :: grchempottn(3,dtset%natom),grcondft(3,dtset%natom),grewtn(3,dtset%natom)
191  real(dp),intent(in) :: grvdw(3,ngrvdw),grnl(3*dtset%natom)
192  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom)
193  real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,dtset%nspden)
194  real(dp),intent(in) :: vxc(nfft,dtset%nspden),vxctau(nfft,dtset%nspden,4*dtset%usekden)
195  real(dp),intent(inout) :: fcart(3,dtset%natom),forold(3,dtset%natom)
196  real(dp),intent(inout) :: vresid(nfft,dtset%nspden),xred(3,dtset%natom)
197  real(dp),intent(out) :: favg(3),gred(3,dtset%natom),gresid(3,dtset%natom)
198  real(dp),intent(out) :: grhf(3,dtset%natom),rprimd(3,3)
199  real(dp),intent(inout) :: grxc(3,dtset%natom)
200  real(dp),intent(out) :: synlgr(3,dtset%natom)
201  type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw)
202  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
203 
204 !Local variables-------------------------------
205 !scalars
206  integer :: coredens_method,coretau_method,fdir,iatom,idir,indx,ipositron,itypat,mu
207  integer :: optatm,optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method
208  real(dp) :: eei_dum1,eei_dum2,ucvol,ucvol_local,vol_element
209  logical :: calc_epaw3_forces, efield_flag
210  logical :: is_hybrid_ncpp
211 !arrays
212  integer :: qprtrb_dum(3)
213  real(dp) :: dummy6(6),ep3(3),fioncart(3),gmet(3,3),gprimd(3,3)
214  real(dp) :: rmet(3,3),strn_dummy6(6),strv_dummy6(6),tsec(2),vprtrb_dum(2)
215  real(dp),allocatable :: atmrho_dum(:),atmvloc_dum(:),dyfrlo_dum(:,:,:)
216  real(dp),allocatable :: dyfrn_dum(:,:,:),dyfrv_dum(:,:,:)
217  real(dp),allocatable :: dyfrx2_dum(:,:,:),eltfrn_dum(:,:),gauss_dum(:,:)
218  real(dp),allocatable :: epawf3red(:,:),fin(:,:),fionred(:,:),grl(:,:),grl_dum(:,:)
219  real(dp),allocatable :: grnl_tmp(:,:),grtn(:,:),grtn_indx(:,:),grxctau(:,:),v_dum(:),vxctotg(:,:)
220  real(dp),allocatable :: xccc3d_dum(:)
221 
222 ! *************************************************************************
223 
224  call timab(69,1,tsec)
225 
226 !Save input value of forces
227  ABI_MALLOC(fin,(3,dtset%natom))
228  fin(:,:)=fcart(:,:)
229 
230 !Compute different geometric tensor, as well as ucvol, from rprimd
231  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
232 
233 !Check if we're in hybrid norm conserving pseudopotential
234  is_hybrid_ncpp=(psps%usepaw==0 .and. &
235 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid()))
236 
237 !=======================================================================
238 !========= Local pseudopotential and core charge contributions =========
239 !=======================================================================
240 
241  ABI_MALLOC(grl,(3,dtset%natom))
242 
243 !Determine by which method the local ionic potential and/or the pseudo core
244 !  charge density contributions have to be computed
245 !Local ionic potential:
246 ! Method 1: PAW
247 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
248  vloc_method=1;if (psps%usepaw==0) vloc_method=2
249  if (dtset%icoulomb>0) vloc_method=2
250  if (psps%usewvl==1) vloc_method=2
251 !Pseudo core charge density:
252 ! Method 1: PAW, nc_xccc_gspace
253 ! Method 2: Norm-conserving PP, wavelets
254  coredens_method=1;if (psps%usepaw==0) coredens_method=2
255  if (psps%nc_xccc_gspace==1) coredens_method=1
256  if (psps%nc_xccc_gspace==0) coredens_method=2
257  if (psps%usewvl==1) coredens_method=2
258  coretau_method=0
259  if (dtset%usekden==1.and.psps%usepaw==1) then
260    coretau_method=1;if (psps%nc_xccc_gspace==0) coretau_method=2
261  end if
262 
263 !Local ionic potential and/or pseudo core charge by method 1
264  if (vloc_method==1.or.coredens_method==1.or.coretau_method==1) then
265    if (psps%nc_xccc_gspace==1.and.psps%usepaw==0.and.is_hybrid_ncpp) then
266      ABI_BUG(' Not yet implemented !')
267    end if
268    call timab(550,1,tsec)
269 !  Allocate (unused) dummy variables, otherwise some compilers complain
270    ABI_MALLOC(gauss_dum,(0,0))
271    ABI_MALLOC(atmrho_dum,(0))
272    ABI_MALLOC(atmvloc_dum,(0))
273    ABI_MALLOC(dyfrn_dum,(0,0,0))
274    ABI_MALLOC(dyfrv_dum,(0,0,0))
275    ABI_MALLOC(eltfrn_dum,(0,0))
276 !  Compute Vxc in reciprocal space
277    if (coredens_method==1.and.n3xccc>0) then
278      ABI_MALLOC(v_dum,(nfft))
279      ABI_MALLOC(vxctotg,(2,nfft))
280      v_dum(:)=vxc(:,1);if (dtset%nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2))
281      call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0)
282      call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
283 &     comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
284      ABI_FREE(v_dum)
285    else
286      ABI_MALLOC(vxctotg,(0,0))
287    end if
288 !  Compute contribution to forces from Vloc and/or pseudo core density
289    optv=0;if (vloc_method==1) optv=1
290    optn=0;if (coredens_method==1) optn=n3xccc/nfft
291    optatm=0;optdyfr=0;optgr=1;optstr=0;optn2=1;opteltfr=0
292    if (vloc_method==1.or.coredens_method==1) then
293      call atm2fft(atindx1,atmrho_dum,atmvloc_dum,dyfrn_dum,dyfrv_dum,&
294 &     eltfrn_dum,gauss_dum,gmet,gprimd,&
295 &     grxc,grl,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
296 &     optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb_dum,&
297 &     dtset%rcut,rhog,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,psps%vlspl,&
298 &     comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
299 &     paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
300    end if
301    if (n3xccc==0.and.coredens_method==1) grxc=zero
302    ABI_FREE(vxctotg)
303    if (dtset%usekden==1.and.coretau_method==1..and.n3xccc>0) then
304 !    Compute contribution to forces from pseudo kinetic energy core density
305      optv=0;optn=1;optn2=4
306      ABI_MALLOC(grxctau,(3,dtset%natom))
307      ABI_MALLOC(grl_dum,(0,0))
308      ABI_MALLOC(v_dum,(nfft))
309      ABI_MALLOC(vxctotg,(2,nfft))
310      v_dum(:)=vxctau(:,1,1);if (dtset%nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxctau(:,2,1))
311      call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0)
312      call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
313 &     comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
314      ABI_FREE(v_dum)
315      call atm2fft(atindx1,atmrho_dum,atmvloc_dum,dyfrn_dum,dyfrv_dum,&
316 &     eltfrn_dum,gauss_dum,gmet,gprimd,&
317 &     grxctau,grl_dum,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
318 &     optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb_dum,&
319 &     dtset%rcut,rhog,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,psps%vlspl,&
320 &     comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
321 &     paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
322      grxc(:,:)=grxc(:,:)+grxctau(:,:)
323      ABI_FREE(grl_dum)
324      ABI_FREE(grxctau)
325      ABI_FREE(vxctotg)
326    end if
327 !  Deallocate temporary arrays
328    ABI_FREE(gauss_dum)
329    ABI_FREE(atmrho_dum)
330    ABI_FREE(atmvloc_dum)
331    ABI_FREE(dyfrn_dum)
332    ABI_FREE(dyfrv_dum)
333    ABI_FREE(eltfrn_dum)
334    call timab(550,2,tsec)
335  end if
336 
337 !Local ionic potential by method 2
338  if (vloc_method==2) then
339    option=2
340    ABI_MALLOC(dyfrlo_dum,(3,3,dtset%natom))
341    ABI_MALLOC(grtn_indx,(3,dtset%natom))
342    ABI_MALLOC(v_dum,(nfft))
343    call mklocl(dtset,dyfrlo_dum,eei_dum1,gmet,gprimd,grtn_indx,gsqcut,dummy6,mgfft,&
344 &   mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,ntypat,option,pawtab,ph1d,psps,&
345 &   qprtrb_dum,rhog,rhor,rprimd,ucvol,vprtrb_dum,v_dum,wvl,wvl_den,xred)
346    do iatom=1,dtset%natom
347 !    Has to use the indexing array atindx1
348      grl(1:3,atindx1(iatom))=grtn_indx(1:3,iatom)
349    end do
350    ABI_FREE(dyfrlo_dum)
351    ABI_FREE(grtn_indx)
352    ABI_FREE(v_dum)
353 !  If gradients are computed in real space, we need to symmetrize the system before summing.
354 !  Rshaltaf: I changed the following line to include surfaces BC
355    if (dtset%icoulomb == 1 .or. dtset%icoulomb == 2) then
356      ABI_MALLOC(grnl_tmp,(3,dtset%natom))
357      call sygrad(grnl_tmp,dtset%natom,grl,dtset%nsym,symrec,indsym)
358      grl(:, :) = grnl_tmp(:, :)
359      ABI_FREE(grnl_tmp)
360    end if
361  end if
362 
363 !Pseudo core electron density by method 2
364  if (coredens_method==2.or.coretau_method==2) then
365    if (n1xccc/=0) then
366      call timab(53,1,tsec)
367      option=2
368      ABI_MALLOC(dyfrx2_dum,(3,3,dtset%natom))
369      ABI_MALLOC(xccc3d_dum,(n3xccc))
370      if (coredens_method==2) then
371        if (is_hybrid_ncpp) then
372          call xchybrid_ncpp_cc(dtset,eei_dum1,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
373 &         dummy6,eei_dum2,xccc3d_dum,grxc=grxc,xcccrc=psps%xcccrc,xccc1d=psps%xccc1d,xred=xred,n1xccc=n1xccc)
374        else
375          if (psps%usewvl==0.and.psps%usepaw==0.and.dtset%icoulomb==0) then
376            call mkcore(dummy6,dyfrx2_dum,grxc,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
377 &           ngfft(1),n1xccc, ngfft(2),ngfft(3),option,rprimd,dtset%typat,ucvol,vxc,&
378 &           psps%xcccrc,psps%xccc1d,xccc3d_dum,xred)
379          else if (psps%usewvl==0.and.(psps%usepaw==1.or.dtset%icoulomb==1)) then
380            call mkcore_alt(atindx1,dummy6,dyfrx2_dum,grxc,dtset%icoulomb,mpi_enreg,dtset%natom,nfft,&
381 &           dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
382 &           ucvol,vxc,psps%xcccrc,psps%xccc1d,xccc3d_dum,xred,pawrad,pawtab,psps%usepaw)
383          else if (psps%usewvl==1.and.psps%usepaw==1) then
384            ucvol_local=ucvol
385 #if defined HAVE_BIGDFT
386 !          ucvol_local=product(wvl_den%denspot%dpbox%hgrids)*real(product(wvl_den%denspot%dpbox%ndims),dp)
387 !          call mkcore_wvl_old(atindx1,dummy6,dyfrx2_dum,wvl%atoms%astruct%geocode,grxc,wvl%h,dtset%natom,&
388 ! &           nattyp,nfft,wvl_den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,:),dtset%nspden,ntypat,&
389 ! &           wvl%Glr%d%n1,wvl%Glr%d%n1i,wvl%Glr%d%n2,wvl%Glr%d%n2i,wvl%Glr%d%n3,wvl_den%denspot%dpbox%n3pi,&
390 ! &           n3xccc,option,pawrad,pawtab,psps%gth_params%psppar,rprimd,ucvol_local,vxc,xccc3d_dum,xred,&
391 ! &           mpi_comm_wvl=mpi_enreg%comm_wvl)
392            call mkcore_wvl(atindx1,dummy6,grxc,dtset%natom,nattyp,nfft,dtset%nspden,ntypat,&
393 &           n1xccc,n3xccc,option,pawrad,pawtab,rprimd,vxc,psps%xccc1d,xccc3d_dum,&
394 &           psps%xcccrc,xred,wvl_den,wvl,mpi_comm_wvl=mpi_enreg%comm_wvl)
395 #endif
396          end if
397        end if
398      end if
399      if (dtset%usekden==1.and.coretau_method==2) then
400        ABI_MALLOC(grxctau,(3,dtset%natom))
401        call mkcore_alt(atindx1,dummy6,dyfrx2_dum,grxctau,dtset%icoulomb,mpi_enreg,dtset%natom,nfft,&
402 &       dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
403 &       ucvol,vxctau(:,:,1),psps%xcccrc,psps%xccc1d,xccc3d_dum,xred,pawrad,pawtab,psps%usepaw,&
404 &       usekden=.true.)
405        grxc(:,:)=grxc(:,:)+grxctau(:,:)
406        ABI_FREE(grxctau)
407      end if
408      ABI_FREE(xccc3d_dum)
409      ABI_FREE(dyfrx2_dum)
410      call timab(53,2,tsec)
411    else
412      grxc(:,:)=zero
413    end if
414  end if
415 
416 !=======================================================================
417 !===================== Nonlocal contributions ==========================
418 !=======================================================================
419 
420 !Only has to apply symmetries
421  ABI_MALLOC(grnl_tmp,(3,dtset%natom))
422  do iatom=1,dtset%natom
423    indx=3*(iatom-1);grnl_tmp(1:3,atindx1(iatom))=grnl(indx+1:indx+3)
424  end do
425  if (dtset%usewvl == 0) then
426    call sygrad(synlgr,dtset%natom,grnl_tmp,dtset%nsym,symrec,indsym)
427  else
428    synlgr = grnl_tmp
429  end if
430  ABI_FREE(grnl_tmp)
431 
432 !=======================================================================
433 !============ Density/potential residual contributions =================
434 !=======================================================================
435 
436  if (dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1.and.abs(dtset%densfor_pred)<=3) then
437    call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,ntypat,1,&
438 &   pawtab,rhor,rprimd,ucvol,vresid,xred,xred,psps%znuclpsp)
439  else if (dtset%usewvl==0.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) then
440    call fresidrsp(atindx1,dtset,gmet,gprimd,gresid,gsqcut,mgfft,&
441 &   mpi_enreg,psps%mqgrid_vl,nattyp,nfft,ngfft,ntypat,psps,pawtab,ph1d,&
442 &   psps%qgrid_vl,rprimd,ucvol,psps%usepaw,vresid,psps%zionpsp,psps%znuclpsp)
443  else
444    gresid(:,:)=zero
445  end if
446 
447 !=======================================================================
448 !======================= Other contributions ===========================
449 !=======================================================================
450 
451 !Ewald energy contribution to forces as already been computed in "ewald"
452 
453 !Potential residual contribution to forces as already been computed (forstr)
454 
455 !Add Berry phase contributions (berryopt == 4/6/7/14/16/17)
456 !(compute the electric field force on the ion cores)
457  efield_flag = (dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. &
458 & dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17)
459  calc_epaw3_forces = (efield_flag .and. dtset%optforces /= 0 .and. psps%usepaw == 1)
460  if ( efield_flag ) then
461    ABI_MALLOC(fionred,(3,dtset%natom))
462    fionred(:,:)=zero
463    do iatom=1,dtset%natom
464      itypat=dtset%typat(iatom)
465 ! force on ion due to electric field, cartesian representation
466      fioncart(:)=psps%ziontypat(itypat)*dtset%efield(:)
467 ! form fionred = rprimd^T * fioncart, note that forces transform
468 ! oppositely to coordinates, because they are derivative with respect to
469 ! coordinates
470      call dgemv('T',3,3,one,rprimd,3,fioncart,1,zero,fionred(1:3,iatom),1)
471 !     do mu=1,3
472 !       fionred(mu,iatom)=rprimd(1,mu)*fioncart(1) &
473 !&       +rprimd(2,mu)*fioncart(2) &
474 !&       +rprimd(3,mu)*fioncart(3)
475 !     end do
476    end do
477  end if
478 
479 !(compute additional F3-type force due to projectors for electric field with PAW)
480  if ( efield_flag .and. calc_epaw3_forces ) then
481    ABI_MALLOC(epawf3red,(3,dtset%natom))
482 ! dtefield%epawf3(iatom,idir,fdir) contains
483    epawf3red(:,:)=zero
484    do iatom=1,dtset%natom
485      do fdir = 1, 3
486        do idir = 1, 3
487 ! vol_element is volume/pt for integration of epawf3. volume is BZ volume
488 ! so 1/ucvol, and number of kpts is nstr(idir)*nkstr(idir)
489          vol_element=one/(ucvol*dtefield%nstr(idir)*dtefield%nkstr(idir))
490          ep3(idir) = vol_element*dtefield%epawf3(iatom,idir,fdir)
491        end do
492        epawf3red(fdir,iatom) = -ucvol*dot_product(dtset%red_efieldbar(1:3),ep3(1:3))
493      end do
494    end do ! end loop over iatom
495  end if
496 
497 !This was incorrect coding. Bug found by Jiawang Hong
498 !if (dtset%berryopt==4) then
499 !allocate(fionred(3,dtset%natom));fionred(:,:)=zero
500 !iatom = 0
501 !do itypat=1,ntypat
502 !do iattyp=1,nattyp(itypat)
503 !iatom=iatom+1
504 !fioncart(:)=psps%ziontypat(itypat)*dtset%efield(:)
505 !do mu=1,3
506 !fionred(mu,iatom)=rprimd(1,mu)*fioncart(1) &
507 !&         +rprimd(2,mu)*fioncart(2) &
508 !&         +rprimd(3,mu)*fioncart(3)
509 !end do
510 !end do
511 !end do
512 !end if
513 
514 !=======================================================================
515 !======= Assemble the various contributions to the forces ==============
516 !=======================================================================
517 
518 !Collect grads of etot wrt reduced coordinates
519 !This gives non-symmetrized Hellman-Feynman reduced gradients
520  ABI_MALLOC(grtn,(3,dtset%natom))
521  grtn(:,:)=grl(:,:)+grchempottn(:,:)+grcondft(:,:)+grewtn(:,:)+synlgr(:,:)+grxc(:,:)
522 
523  if (usefock==1 .and. associated(fock)) then
524    if (fock%fock_common%optfor) then
525      grtn(:,:)=grtn(:,:)+fock%fock_common%forces(:,:)
526    end if
527  end if
528 
529  if (ngrvdw==dtset%natom) grtn(:,:)=grtn(:,:)+grvdw(:,:)
530 ! note that fionred is subtracted, because it really is a force and we need to
531 ! turn it back into a gradient. The gred2fcart routine below includes the minus
532 ! sign to convert gradients back to forces
533  if ( efield_flag ) grtn(:,:)=grtn(:,:)-fionred(:,:)
534 ! epawf3red is added, because it actually is a gradient, not a force
535  if ( efield_flag .and. calc_epaw3_forces ) grtn(:,:) = grtn(:,:) + epawf3red(:,:)
536 
537 !Symmetrize explicitly for given space group and store in grhf :
538  call sygrad(grhf,dtset%natom,grtn,dtset%nsym,symrec,indsym)
539 
540 !If residual forces are too large, there must be a problem: cancel them !
541  if (dtset%usewvl==0.and.abs(dtset%densfor_pred)>0.and.abs(dtset%densfor_pred)/=5) then
542    do iatom=1,dtset%natom
543      do mu=1,3
544        if (abs(gresid(mu,iatom))>10000._dp*abs(grtn(mu,iatom))) gresid(mu,iatom)=zero
545      end do
546    end do
547  end if
548 
549 !Add residual potential correction
550  grtn(:,:)=grtn(:,:)+gresid(:,:)
551 
552 !Additional stuff for electron-positron
553  ipositron=0
554  if (present(electronpositron)) then
555    if (associated(electronpositron)) then
556      if (allocated(electronpositron%gred_ep)) ipositron=electronpositron_calctype(electronpositron)
557    end if
558  end if
559  if (abs(ipositron)==1) then
560    grtn(:,:)=grtn(:,:)-grxc(:,:)-grchempottn(:,:)-grcondft(:,:)-grewtn(:,:)-gresid(:,:)-two*grl(:,:)
561 !  grtn(:,:)=grtn(:,:)-grxc(:,:)-grewtn(:,:)-gresid(:,:)-two*grl(:,:)
562    grl(:,:)=-grl(:,:);grxc(:,:)=zero;gresid(:,:)=zero
563    if (ngrvdw==dtset%natom) grtn(:,:)=grtn(:,:)-grvdw(:,:)
564    if ( dtset%berryopt== 4 .or. dtset%berryopt== 6 .or. dtset%berryopt== 7 .or. &
565 &   dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17)  then
566      grtn(:,:)=grtn(:,:)+fionred(:,:)
567      fionred(:,:)=zero
568    end if
569  end if
570  if (ipositron>0) grtn(:,:)=grtn(:,:)+electronpositron%gred_ep(:,:)
571 
572 !Symmetrize all grads explicitly for given space group:
573  if (dtset%usewvl == 0) then
574    call sygrad(gred,dtset%natom,grtn,dtset%nsym,symrec,indsym)
575  else
576    gred = grtn
577  end if
578 
579 !Conversion to cartesian coordinates (bohr) AND
580 !Subtract off average force from each force component
581 !to avoid spurious drifting of atoms across cell.
582 ! notice that gred2fcart multiplies gred by -1 to convert it
583 ! from a gradient (input) to a force (output)
584 
585  call gred2fcart(favg,(dtset%jellslab==0 .and. dtset%nzchempot==0),fcart,gred,gprimd,dtset%natom)
586 
587 !Compute maximal force and maximal difference
588  maxfor=zero;diffor=zero
589  do iatom=1,dtset%natom
590    do mu=1,3
591      if (dtset%iatfix(mu,iatom) /= 1) then
592        maxfor=max(maxfor,abs(fcart(mu,iatom)))
593        diffor=max(diffor,abs(fcart(mu,iatom)-fin(mu,iatom)))
594      else if (dtset%ionmov==4 .or. dtset%ionmov==5) then
595 !      Make the force vanish on fixed atoms when ionmov=4 or 5
596 !      This is because fixing of atom cannot be imposed at the
597 !      level of a routine similar to brdmin or moldyn for these options.
598        fcart(mu,iatom)=zero
599      end if
600    end do
601  end do
602 
603 !Apply any generalized constraints to the forces
604  if (dtset%nconeq>0) call constrf(diffor,fcart,forold,gred,dtset%iatfix,dtset%ionmov,maxfor,&
605 & dtset%natom,dtset%nconeq,dtset%prtvol,rprimd,dtset%wtatcon,xred)
606 
607 !=======================================================================
608 !Memory deallocations
609  ABI_FREE(grl)
610  ABI_FREE(grtn)
611  ABI_FREE(fin)
612  if ( efield_flag )  then
613    ABI_FREE(fionred)
614    if ( calc_epaw3_forces ) then
615      ABI_FREE(epawf3red)
616    end if
617  end if
618 
619  call timab(69,2,tsec)
620 
621 end subroutine forces

ABINIT/fresid [ Functions ]

[ Top ] [ Functions ]

NAME

 fresid

FUNCTION

 If option=1, compute the forces due to the residual of the potential
 If option=2, generate approximate new density from old one,
              old atomic positions, and new atomic positions

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
  | typat(natom)=integer type for each atom in cell
  | usepaw= 0 for non paw calculation; =1 for paw calculation
  | xclevel= level of the XC functional
 mpi_enreg=information about MPI parallelization
 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
 ntypat=number of types of atoms in cell.
 option=see below
 pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
 rhor(nfft,nspden)=electron density in electrons/bohr**3 (slices of it if FTT parallelism).
 rprimd(3,3)=dimensional primitive translation vectors (bohr)
 ucvol=unit cell volume (bohr**3).
 xred_new(3,natom)=new reduced coordinates for atoms in unit cell
 xred_old(3,natom)=old reduced coordinates for atoms in unit cell
 znucl(ntypat)=real(dp), atomic number of atom type

OUTPUT

 gresid(3,natom)=forces due to the residual of the potential

SIDE EFFECTS

 work(nfft,nspden)=functions on the fft grid (slices of it if FTT parallelism):
  if option==1, the POTENTIAL residual is input
  if option==2, the interpolated density is output

NOTES

 FFT parallelism:
 At the beginning of this routine, the plane-waves are ditributed over all the processors.
 In the main part, all the processors perform the same calculations over the whole FFT grid.
 At the end, each processor gets its part of the whole FFT grid.
 These modifications are not efficient when large FFT grids are used.
 So they have to be considered as a first step before a comprehensive parallelization of this routine.

SOURCE

 879 subroutine fresid(dtset,gresid,mpi_enreg,nfft,ngfft,ntypat,option,&
 880 &                 pawtab,rhor,rprimd,ucvol,work,xred_new,xred_old,znucl)
 881 
 882 !Arguments ------------------------------------
 883 !scalars
 884  integer,intent(in) :: nfft,ntypat,option
 885  real(dp),intent(in) :: ucvol
 886  type(MPI_type),intent(in) :: mpi_enreg
 887  type(dataset_type),intent(in) :: dtset
 888 !arrays
 889  integer,intent(in) :: ngfft(18)
 890  real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3)
 891  real(dp),intent(in) :: xred_new(3,dtset%natom),xred_old(3,dtset%natom)
 892  real(dp),intent(in) :: znucl(ntypat)
 893  real(dp),intent(inout) :: work(nfft,dtset%nspden)
 894  real(dp),intent(out) :: gresid(3,dtset%natom)
 895  type(pawtab_type),intent(in) :: pawtab(ntypat*dtset%usepaw)
 896 
 897 !Local variables-------------------------------
 898 !real(dp), parameter :: app_remain=0.001_dp
 899 !scalars
 900  integer,parameter :: natnum=110
 901  integer :: atmove,i1,i1_new,i1m,i1p,i2,i2_new,i2m,i2p,i3,i3_new,i3m,i3p
 902  integer :: iatom,ifft,ifft_new,iloop,ind2m,ind2m3m,ind2p,ind2p3p,ind3m,ind3p
 903  integer :: index,index_new,ishift,ishift1,ishift2,ishift3,ispden,ixp,mshift,mu
 904  integer :: n1,n2,n3,n4,nfft_tmp,nfftot,nu,quit
 905  real(dp),parameter :: app_remain=0.01_dp
 906  real(dp) :: diff_rem1,diff_rem2,diff_rem3,difmag,difmag2
 907  real(dp) :: difmag2_fact,difmag2_part,drho1,drho11,drho12,drho13,drho14
 908  real(dp) :: drho1dn,drho1mx,drho1my,drho1mz,drho1tot,drho1up,drho2,drho21
 909  real(dp) :: drho22,drho23,drho24,drho2dn,drho2mx,drho2my,drho2mz,drho2tot
 910  real(dp) :: drho2up,drho3,drho31,drho32,drho33,drho34,drho3dn,drho3mx,drho3my
 911  real(dp) :: drho3mz,drho3tot,drho3up,drhox00,drhox01,drhox10,drhox11,drhoxy0
 912  real(dp) :: drhoxy1,drhoxyz,fact,range,range2,rcov,rcov2,rcovm1,rdiff1
 913  real(dp) :: rdiff2,rdiff3,vresid1,vresid2,vresid3,vresid4,xx
 914  type(atomdata_t) :: atom
 915 !arrays
 916  integer :: diff_igrid(3),igrid(3),irange(3)
 917  integer,allocatable :: ii(:,:)
 918  real(dp) :: diff_grid(3),diff_rem(3),diff_tau(3),diff_xred(3),lencp(3)
 919  real(dp) :: rho_tot(4),rhosum(4),rmet(3,3),scale(3),tau(3)
 920  real(dp),allocatable :: approp(:),atmrho(:,:),rhor_tot(:,:),rrdiff(:,:)
 921  real(dp),allocatable :: work_tot(:,:)
 922  logical,allocatable :: my_sphere(:)
 923 
 924 ! *************************************************************************
 925 
 926 !Compute lengths of cross products for pairs of primitive
 927 !translation vectors (used in setting index search range below)
 928  lencp(1)=cross_fr(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
 929 & rprimd(1,3),rprimd(2,3),rprimd(3,3))
 930  lencp(2)=cross_fr(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
 931 & rprimd(1,1),rprimd(2,1),rprimd(3,1))
 932  lencp(3)=cross_fr(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
 933 & rprimd(1,2),rprimd(2,2),rprimd(3,2))
 934 
 935 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
 936 !(recall ucvol=R1.(R2xR3))
 937  scale(:)=ucvol/lencp(:)
 938 
 939 !initialize diff_igrid, otherwise valgrind complains
 940  diff_igrid=0
 941 
 942 !Compute metric tensor in real space rmet
 943  do nu=1,3
 944    rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu)
 945  end do
 946 
 947 !FFT parallelization: Starting from now, calculations are performed on the whole FFT grid
 948 !and no more on slices. The nfft variable becomes nfft_tmp until the end
 949  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
 950  n4=n3/mpi_enreg%nproc_fft
 951  nfftot=PRODUCT(ngfft(1:3));nfft_tmp=nfftot
 952  if(mpi_enreg%paral_kgb==1) then
 953    ABI_MALLOC(rhor_tot,(nfftot,dtset%nspden))
 954    ABI_MALLOC(work_tot,(nfftot,dtset%nspden))
 955    do ispden=1,dtset%nspden
 956      call pre_gather(rhor(:,ispden),rhor_tot(:,ispden),n1,n2,n3,n4,mpi_enreg)
 957      call pre_gather(work(:,ispden),work_tot(:,ispden),n1,n2,n3,n4,mpi_enreg)
 958    end do
 959  end if
 960 
 961  gresid(1:3,1:dtset%natom)=0.0_dp
 962  quit=0
 963 
 964 !Initialize appropriation function
 965  ABI_MALLOC(approp,(nfft_tmp))
 966  ABI_MALLOC(atmrho,(nfft_tmp,dtset%nspden))
 967  ABI_MALLOC(my_sphere,(nfft_tmp))
 968 
 969  approp(:)=app_remain
 970 !First loop over atoms in unit cell : build appropriation function
 971 !Second loop : compute forces
 972  do iloop=1,2
 973 
 974 !  Take into account the remaining density
 975    if(option==2 .and. iloop==2)then
 976      if(mpi_enreg%paral_kgb==1) then
 977 !      FFT parallelization: All the processors perform the same calculation.
 978 !      We divided by nproc_fft in order to "remove" the xmpi_sum made in mean_fftr
 979        do ispden=1,dtset%nspden
 980          do ifft=1,nfft_tmp
 981            work_tot(ifft,ispden)=rhor_tot(ifft,ispden)*approp(ifft)*app_remain
 982          end do
 983        end do
 984        call mean_fftr(work_tot,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
 985        rhosum(1:dtset%nspden)=rhosum(1:dtset%nspden)/mpi_enreg%nproc_fft
 986      else
 987        do ispden=1,dtset%nspden
 988          do ifft=1,nfft_tmp
 989            work(ifft,ispden)=rhor(ifft,ispden)*approp(ifft)*app_remain
 990          end do
 991        end do
 992        call mean_fftr(work,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
 993      end if
 994 
 995 !    This will be used to restore proper normalization of density
 996      rho_tot(1:dtset%nspden)=rhosum(1:dtset%nspden)*nfftot
 997    end if
 998 
 999    do iatom=1,dtset%natom
1000 
1001 !    Get the covalent radius
1002      call atomdata_from_znucl(atom,znucl(dtset%typat(iatom)))
1003      rcov = atom%rcov
1004 !    PAW choose PAW radius instead...
1005      if (dtset%usepaw==1) rcov=max(rcov,pawtab(dtset%typat(iatom))%rpaw)
1006 
1007 !    Set search range
1008      rcov2=rcov**2
1009      range=2._dp*rcov
1010      range2=range**2
1011      rcovm1=1.0_dp/rcov
1012 
1013 !    Use range to compute an index range along R(1:3)
1014 !    (add 1 to make sure it covers full range)
1015      irange(1)=1+nint((range/scale(1))*dble(n1))
1016      irange(2)=1+nint((range/scale(2))*dble(n2))
1017      irange(3)=1+nint((range/scale(3))*dble(n3))
1018 
1019 !    Allocate ii and rrdiff
1020      mshift=2*maxval(irange(1:3))+1
1021      ABI_MALLOC(ii,(mshift,3))
1022      ABI_MALLOC(rrdiff,(mshift,3))
1023 
1024 !    Consider each component in turn
1025      do mu=1,3
1026 
1027 !      Convert reduced coord of given atom to [0,1)
1028        tau(mu)=mod(xred_old(mu,iatom)+1._dp-aint(xred_old(mu,iatom)),1._dp)
1029 
1030 !      Use tau to find nearest grid point along R(mu)
1031 !      (igrid=0 is the origin; shift by 1 to agree with usual index)
1032        igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
1033 
1034 !      Set up a counter that explore the relevant range
1035 !      of points around the atom
1036        ishift=0
1037        do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
1038          ishift=ishift+1
1039          ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
1040          rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu)
1041        end do
1042 
1043 !      If option 2, set up quantities related with the change of atomic coordinates
1044        if(option==2 .and. iloop==2)then
1045          diff_xred(mu)=xred_new(mu,iatom)-xred_old(mu,iatom)
1046 !        Convert to [0,1)
1047          diff_tau(mu)=mod(diff_xred(mu)+1._dp-aint(diff_xred(mu)),1._dp)
1048 !        Convert to [0,ngfft)
1049          diff_grid(mu)=diff_tau(mu)*dble(ngfft(mu))
1050 !        Integer part
1051          diff_igrid(mu)=int(diff_grid(mu))
1052 !        Compute remainder
1053          diff_rem(mu)=diff_grid(mu)-diff_igrid(mu)
1054 
1055 !        DEBUG
1056 !        write(std_out,*)' mu,diff',mu,diff_igrid(mu),diff_rem(mu)
1057 !        ENDDEBUG
1058 
1059        end if
1060 
1061 !      End loop on mu
1062      end do
1063 
1064 !    May be the atom is fixed
1065      atmove=1
1066      if(option==2 .and. iloop==2)then
1067        if(diff_xred(1)**2+diff_xred(2)**2+diff_xred(3)**2 < 1.0d-24)then
1068          atmove=0
1069        else
1070          diff_rem1=diff_rem(1)
1071          diff_rem2=diff_rem(2)
1072          diff_rem3=diff_rem(3)
1073        end if
1074      end if
1075 
1076 !    If second loop, initialize atomic density, and the variable
1077 !    that says whether a fft point belongs to the sphere of the atom
1078      if(iloop==2) then
1079        atmrho(:,:)=0.0_dp
1080        my_sphere(:)=.false.
1081      end if
1082 
1083 !    Conduct triple loop over restricted range of grid points for iatom
1084 
1085      do ishift3=1,1+2*irange(3)
1086 !      map back to [1,ngfft(3)] for usual fortran index in unit cell
1087        i3=ii(ishift3,3)
1088        i3m=i3-1 ; if(i3==1)i3m=n3
1089        i3p=i3+1 ; if(i3==n3)i3p=1
1090 
1091 !      find vector from atom location to grid point (reduced)
1092        rdiff3=rrdiff(ishift3,3)
1093 
1094        do ishift2=1,1+2*irange(2)
1095          i2=ii(ishift2,2)
1096          i2m=i2-1 ; if(i2==1)i2m=n2
1097          i2p=i2+1 ; if(i2==n2)i2p=1
1098          index=n1*(i2-1+n2*(i3-1))
1099          ind3m=n1*(i2-1+n2*(i3m-1))
1100          ind3p=n1*(i2-1+n2*(i3p-1))
1101          ind2m=n1*(i2m-1+n2*(i3-1))
1102          ind2p=n1*(i2p-1+n2*(i3-1))
1103          ind2p3p=n1*(i2p-1+n2*(i3p-1))
1104 
1105          rdiff2=rrdiff(ishift2,2)
1106 !        Prepare the computation of difmag2
1107          difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
1108 &         +2.0_dp*rmet(3,2)*rdiff3*rdiff2
1109          difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
1110 
1111          do ishift1=1,1+2*irange(1)
1112            rdiff1=rrdiff(ishift1,1)
1113 
1114 !          Compute (rgrid-tau-Rprim)**2
1115            difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
1116 
1117 !          Only accept contribution inside defined range
1118 !          This condition means that x, calculated below, cannot exceed 2.0_dp
1119            if (difmag2<range2) then
1120 
1121 !            Will compute contribution to appropriation function based on
1122 !            rcov2, range2 and difmag2
1123              i1=ii(ishift1,1)
1124              ifft=i1+index
1125 
1126              if(iloop==1)then
1127 
1128 !              Build appropriation function
1129                if (difmag2<rcov2)then
1130                  approp(ifft)=approp(ifft)+1.0_dp
1131                else
1132                  difmag=sqrt(difmag2)
1133                  xx=difmag*rcovm1
1134 !                The following function is 1. at xx=1, 0. at xx=2, with vanishing
1135 !                derivatives at these points.
1136                  approp(ifft)=approp(ifft)+((2.0_dp*xx-9.0_dp)*xx+12.0_dp)*xx-4.0_dp
1137                end if
1138 
1139              else
1140 
1141                if (difmag2<rcov2) then
1142                  fact=one
1143                else
1144                  difmag=sqrt(difmag2)
1145                  xx=difmag*rcovm1
1146                  fact=((2.0_dp*xx-9.0_dp)*xx+12.0_dp)*xx-4.0_dp
1147                end if
1148 
1149 !              Build atomic density
1150                if(mpi_enreg%paral_kgb==1) then
1151                  atmrho(ifft,1:dtset%nspden)=atmrho(ifft,1:dtset%nspden) &
1152 &                 +rhor_tot(ifft,1:dtset%nspden)*fact*approp(ifft)
1153                else
1154                  atmrho(ifft,1:dtset%nspden)=atmrho(ifft,1:dtset%nspden) &
1155 &                 +rhor(ifft,1:dtset%nspden)*fact*approp(ifft)
1156                end if
1157 
1158 !              Compute the sphere of the atom : it is different for
1159 !              option 1 and for option 2
1160                i1p=i1+1 ; if(i1==n1)i1p=1
1161                if(option==1)then
1162                  i1m=i1-1 ; if(i1==1)i1m=n1
1163                  my_sphere(ifft)=.true.
1164                  my_sphere(i1p+index)=.true. ; my_sphere(i1m+index)=.true.
1165                  my_sphere(i1+ind2p)=.true. ; my_sphere(i1+ind2m)=.true.
1166                  my_sphere(i1+ind3p)=.true. ; my_sphere(i1+ind3m)=.true.
1167                else
1168                  my_sphere(ifft)=.true. ; my_sphere(i1p+index)=.true.
1169                  my_sphere(i1+ind2p)=.true. ; my_sphere(i1p+ind2p)=.true.
1170                  my_sphere(i1+ind3p)=.true. ; my_sphere(i1p+ind3p)=.true.
1171                  my_sphere(i1+ind2p3p)=.true. ; my_sphere(i1p+ind2p3p)=.true.
1172                end if
1173 
1174              end if
1175 
1176 !            End of condition on the range
1177            end if
1178 
1179 !          End loop on ishift1
1180          end do
1181 
1182 !        End loop on ishift2
1183        end do
1184 
1185 !      End loop on ishift3
1186      end do
1187 !    At the end of the second loop for each atom, compute the force
1188 !    from the atomic densities, or translate density.
1189 !    In the first case, use a two-point finite-difference approximation,
1190 !    since this calculation serves only to decrease the error,
1191 !    and should not be very accurate, but fast.
1192 !    In the second case, using a crude trilinear interpolation scheme
1193 !    for the same reason.
1194 !
1195 !    The section is skipped if option==2 and the atom is fixed
1196      if(iloop==2 .and. (option==1 .or. atmove==1) )then
1197 
1198        do i3=1,n3
1199          i3m=i3-1 ; if(i3==1)i3m=n3
1200          i3p=i3+1 ; if(i3==n3)i3p=1
1201 !        note: diff_igrid is only set  if(option==2 .and. iloop==2)
1202          i3_new=i3+diff_igrid(3) ; if(i3_new > n3)i3_new=i3_new-n3
1203          do i2=1,n2
1204            i2m=i2-1 ; if(i2==1)i2m=n2
1205            i2p=i2+1 ; if(i2==n2)i2p=1
1206            i2_new=i2+diff_igrid(2) ; if(i2_new > n2)i2_new=i2_new-n2
1207            index=n1*(i2-1+n2*(i3-1))
1208            index_new=n1*(i2_new-1+n2*(i3_new-1))
1209            ind3m=n1*(i2-1+n2*(i3m-1))
1210            ind3p=n1*(i2-1+n2*(i3p-1))
1211            ind2m=n1*(i2m-1+n2*(i3-1))
1212            ind2p=n1*(i2p-1+n2*(i3-1))
1213            ind2m3m=n1*(i2m-1+n2*(i3m-1))
1214            do i1=1,n1
1215              ifft=i1+index
1216              if(my_sphere(ifft))then
1217 
1218                i1m=i1-1 ; if(i1==1)i1m=n1
1219 
1220                if(option==1)then
1221 !                Treat option 1 : computation of residual forces
1222                  i1p=i1+1 ; if(i1==n1)i1p=1
1223 !                Distinguish spin-unpolarized and spin-polarized
1224                  if(dtset%nspden==1)then ! Non magnetic
1225 !                  Note that the factor needed to obtain a true finite difference
1226 !                  estimation of the derivative will be applied afterwards, for speed
1227                    drho1=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1228                    drho2=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1229                    drho3=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1230                    if(mpi_enreg%paral_kgb==1) then
1231                      vresid1=work_tot(ifft,1)
1232                    else
1233                      vresid1=work(ifft,1)
1234                    end if
1235                    gresid(1,iatom)=gresid(1,iatom)+drho1*vresid1
1236                    gresid(2,iatom)=gresid(2,iatom)+drho2*vresid1
1237                    gresid(3,iatom)=gresid(3,iatom)+drho3*vresid1
1238                  else if(dtset%nspden==2) then ! Collinear magnetism
1239                    drho1tot=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1240                    drho2tot=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1241                    drho3tot=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1242                    drho1up=atmrho(i1p+index,2)-atmrho(i1m+index,2)
1243                    drho2up=atmrho(i1+ind2p,2) -atmrho(i1+ind2m,2)
1244                    drho3up=atmrho(i1+ind3p,2) -atmrho(i1+ind3m,2)
1245                    drho1dn=drho1tot-drho1up
1246                    drho2dn=drho2tot-drho2up
1247                    drho3dn=drho3tot-drho3up
1248                    if(mpi_enreg%paral_kgb==1) then
1249                      vresid1=work_tot(ifft,1)
1250                      vresid2=work_tot(ifft,2)
1251                    else
1252                      vresid1=work(ifft,1)
1253                      vresid2=work(ifft,2)
1254                    end if
1255                    gresid(1,iatom)=gresid(1,iatom)+drho1up*vresid1+drho1dn*vresid2
1256                    gresid(2,iatom)=gresid(2,iatom)+drho2up*vresid1+drho2dn*vresid2
1257                    gresid(3,iatom)=gresid(3,iatom)+drho3up*vresid1+drho3dn*vresid2
1258                  else ! Non-collinear magnetism
1259                    drho1tot=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1260                    drho1mx =atmrho(i1p+index,2)-atmrho(i1m+index,2)
1261                    drho1my =atmrho(i1p+index,3)-atmrho(i1m+index,3)
1262                    drho1mz =atmrho(i1p+index,4)-atmrho(i1m+index,4)
1263                    drho2tot=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1264                    drho2mx =atmrho(i1+ind2p,2) -atmrho(i1+ind2m,2)
1265                    drho2my =atmrho(i1+ind2p,3) -atmrho(i1+ind2m,3)
1266                    drho2mz =atmrho(i1+ind2p,4) -atmrho(i1+ind2m,4)
1267                    drho3tot=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1268                    drho3mx =atmrho(i1+ind3p,2) -atmrho(i1+ind3m,2)
1269                    drho3my =atmrho(i1+ind3p,3) -atmrho(i1+ind3m,3)
1270                    drho3mz =atmrho(i1+ind3p,4) -atmrho(i1+ind3m,4)
1271                    drho11=half*(drho1tot+drho1mz)
1272                    drho12=half*(drho1tot-drho1mz)
1273                    drho13= half*drho1mx
1274                    drho14=-half*drho1my
1275                    drho21=half*(drho2tot+drho2mz)
1276                    drho22=half*(drho2tot-drho2mz)
1277                    drho23= half*drho2mx
1278                    drho24=-half*drho2my
1279                    drho31=half*(drho3tot+drho3mz)
1280                    drho32=half*(drho3tot-drho3mz)
1281                    drho33= half*drho3mx
1282                    drho34=-half*drho3my
1283                    if(mpi_enreg%paral_kgb==1) then
1284                      vresid1=work_tot(ifft,1)
1285                      vresid2=work_tot(ifft,2)
1286                      vresid3=work_tot(ifft,3)
1287                      vresid4=work_tot(ifft,4)
1288                    else
1289                      vresid1=work(ifft,1)
1290                      vresid2=work(ifft,2)
1291                      vresid3=work(ifft,3)
1292                      vresid4=work(ifft,4)
1293                    end if
1294                    gresid(1,iatom)=gresid(1,iatom)+drho11*vresid1+drho12*vresid2+two*(drho13*vresid3+drho14*vresid4)
1295                    gresid(2,iatom)=gresid(2,iatom)+drho21*vresid1+drho22*vresid2+two*(drho23*vresid3+drho24*vresid4)
1296                    gresid(3,iatom)=gresid(3,iatom)+drho31*vresid1+drho32*vresid2+two*(drho33*vresid3+drho34*vresid4)
1297                  end if
1298 !                Treat the case option==2 now : trilinear interpolation of the density
1299                else
1300                  i1_new=i1+diff_igrid(1) ; if(i1_new > n1)i1_new=i1_new-n1
1301                  ifft_new=i1_new+index_new
1302                  do ispden=1,dtset%nspden
1303                    drhox00=(atmrho(i1m+index,ispden)-atmrho(i1+index,ispden))*diff_rem1 &
1304 &                   +atmrho(i1+index,ispden)
1305                    drhox10=(atmrho(i1m+ind2m,ispden)-atmrho(i1+ind2m,ispden))*diff_rem1 &
1306 &                   +atmrho(i1+ind2m,ispden)
1307                    drhox01=(atmrho(i1m+ind3m,ispden)-atmrho(i1+ind3m,ispden))*diff_rem1 &
1308 &                   +atmrho(i1+ind3m,ispden)
1309                    drhox11=(atmrho(i1m+ind2m3m,ispden)-atmrho(i1+ind2m3m,ispden))*diff_rem1 &
1310 &                   +atmrho(i1+ind2m3m,ispden)
1311                    drhoxy0=(drhox10-drhox00)*diff_rem2+drhox00
1312                    drhoxy1=(drhox11-drhox01)*diff_rem2+drhox01
1313                    drhoxyz=(drhoxy1-drhoxy0)*diff_rem3+drhoxy0
1314                    if(mpi_enreg%paral_kgb==1) then
1315                      work_tot(ifft_new,ispden)=work_tot(ifft_new,ispden)+drhoxyz
1316                    else
1317                      work(ifft_new,ispden)=work(ifft_new,ispden)+drhoxyz
1318                    end if
1319                    rho_tot(ispden)=rho_tot(ispden)+drhoxyz
1320                  end do
1321                end if
1322 
1323 !              End condition of belonging to the sphere of influence of the atom
1324              end if
1325            end do
1326          end do
1327        end do
1328 !      The finite-difference factor applied here also take
1329 !      into account diverse factors
1330        fact=-ucvol/dble(nfftot)
1331        gresid(1,iatom)=gresid(1,iatom)*dble(n1)*.5_dp*fact
1332        gresid(2,iatom)=gresid(2,iatom)*dble(n2)*.5_dp*fact
1333        gresid(3,iatom)=gresid(3,iatom)*dble(n3)*.5_dp*fact
1334      end if
1335 
1336 !    Update work if the atom is fixed.
1337      if(iloop==2 .and. option==2 .and. atmove==0)then
1338        if(mpi_enreg%paral_kgb==1) then
1339 !        FFT parallelization: All the processors perform the same calculation.
1340 !        We divided by nproc_fft in order to "remove" the xmpi_sum made in mean_fftr
1341          do ispden=1,dtset%nspden
1342            do ifft=1,nfft_tmp
1343              work_tot(ifft,ispden)=work_tot(ifft,ispden)+atmrho(ifft,ispden)
1344            end do
1345          end do
1346          call mean_fftr(atmrho,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1347          rhosum(1:dtset%nspden)=rhosum(1:dtset%nspden)/mpi_enreg%nproc_fft
1348        else
1349          do ispden=1,dtset%nspden
1350            do ifft=1,nfft_tmp
1351              work(ifft,ispden)=work(ifft,ispden)+atmrho(ifft,ispden)
1352            end do
1353          end do
1354          call mean_fftr(atmrho,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1355        end if
1356 
1357        rho_tot(1:dtset%nspden)=rho_tot(1:dtset%nspden)+rhosum(1:dtset%nspden)*nfftot
1358      end if
1359 
1360      ABI_FREE(ii)
1361      ABI_FREE(rrdiff)
1362 
1363 !    End loop on atoms
1364    end do
1365 
1366 !  DEBUG
1367 !  if(option==2)then
1368 !  if(iloop==1)then
1369 !  write(std_out,*)' fresid : rhor, approp'
1370 !  do ifft=1,n1
1371 !  write(std_out,*)ifft,rhor(ifft,1),approp(ifft)
1372 !  end do
1373 !  end if
1374 !  if(iloop==2)then
1375 !  write(std_out,*)' fresid : rhor, approp, work(:,:)'
1376 !  do ifft=1,n1
1377 !  write(std_out,'(i4,3es18.8)' )ifft,rhor(ifft,1),approp(ifft),work(ifft,1)
1378 !  end do
1379 !  do ifft=1,nfft_tmp
1380 !  if(work(ifft,1)<0.0_dp)then
1381 !  write(std_out,*)' f_fft negative value',work(ifft,1),' for ifft=',ifft
1382 !  end if
1383 !  if(rhor(ifft,1)<0.0_dp)then
1384 !  write(std_out,*)' rhor  negative value',rhor(ifft,1),' for ifft=',ifft
1385 !  end if
1386 !  end do
1387 !  end if
1388 !  end if
1389 !  ENDDEBUG
1390 
1391    if(quit==1)exit
1392 
1393 !  At the end of the first loop, where the appropriation function is generated,
1394 !  invert it, to save cpu time later.
1395    if(iloop==1)approp(:)=1.0_dp/approp(:)
1396 
1397 !  End first or second pass through atoms
1398  end do
1399 
1400 !Restore proper normalisation of density
1401 !(Non-collinear magnetism: n, mx,my,mz integral conservation)
1402  if(option==2)then
1403    if(mpi_enreg%paral_kgb==1) then
1404 !    FFT parallelization: All the processors perform the same calculation.
1405 !    We divided by nproc_fft in order to "remove" the xmpi_sum made in mean_fftr
1406 !    Trangel: mpicomm now is optional in mean_fftr, no need to divide over nproc_fft
1407      call mean_fftr(rhor_tot,rhosum,nfft_tmp,nfftot,dtset%nspden)
1408 !    rhosum(1:dtset%nspden)=rhosum(1:dtset%nspden)/mpi_enreg%nproc_fft
1409    else
1410      call mean_fftr(rhor,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1411    end if
1412 !  "!OCL NOPREEX" to avoid zero division after optimization (-Of) by MM
1413 !  (Even if nspden=1, "1.0/rho_tot" will appear on vpp fujitsu
1414 !  OCL NOPREEX
1415    if(mpi_enreg%paral_kgb==1) then
1416      do ispden=1,dtset%nspden
1417        fact=rhosum(ispden)*dble(nfftot)/rho_tot(ispden)
1418        work_tot(:,ispden)=fact*work_tot(:,ispden)
1419        call pre_scatter(work(:,ispden),work_tot(:,ispden),n1,n2,n3,n4,mpi_enreg)
1420      end do
1421    else
1422      do ispden=1,dtset%nspden
1423        fact=rhosum(ispden)*dble(nfftot)/rho_tot(ispden)
1424        work(:,ispden)=fact*work(:,ispden)
1425      end do
1426    end if
1427 !  DEBUG
1428 !  Here, zero all the hard work, for checking purposes !
1429 !  work(:,:)=rhor(:,:)
1430 !  ENDDEBUG
1431  end if
1432 
1433  ABI_FREE(approp)
1434  ABI_FREE(atmrho)
1435  ABI_FREE(my_sphere)
1436  if(mpi_enreg%paral_kgb==1) then
1437    ABI_FREE(rhor_tot)
1438    ABI_FREE(work_tot)
1439  end if
1440 
1441 !DEBUG
1442 !write(std_out,*)' fresid : exit '
1443 !do iatom=1,dtset%natom
1444 !write(std_out,*)iatom,gresid(1:3,iatom)
1445 !enddo
1446 !ENDDEBUG
1447 
1448  contains
1449 
1450    function cross_fr(xx,yy,zz,aa,bb,cc)
1451 !Define magnitude of cross product of two vectors
1452    real(dp) :: cross_fr
1453    real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
1454    cross_fr=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
1455  end function cross_fr
1456 
1457 end subroutine fresid

ABINIT/fresidrsp [ Functions ]

[ Top ] [ Functions ]

NAME

 fresidrsp

FUNCTION

 Compute the forces due to the residual of the potential (or density)
 in RECIPROCAL SPACE, using
  - the atomic density read in psp file (PAW or NC with nctval_spl e.g. psp8 format)
  - a gaussian atomic density (norm-conserving psps if nctval_spl is not available)

INPUTS

 atindx1(natom)=index table for atoms, inverse of atindx
 dtset <type(dataset_type)>=all input variables in this dataset
  | densty(ntypat,4)=parameters for initialisation of the density of each atom type
  | icoulomb=0 periodic treatment of Hartree potential, 1 use of Poisson solver
  | ixc= choice of exchange-correlation scheme
  | natom=number of atoms in cell.
  | nspden=number of spin-density components
  | typat(natom)=integer type for each atom in cell
 gmet(3,3)=reciprocal space metric
 gprimd(3,3)=reciprocal space dimensional primitive translations
 gsqcut=cutoff value on G**2 for sphere inside fft box
 mgfft=maximum size of 1D FFTs
 mpi_enreg=information about MPI parallelization
 mqgrid=number of grid pts in q array for atomic density spline n^AT(q)
 nattyp(ntypat)=number of atoms of each type in cell
 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
 ntypat=number of types of atoms in cell.
 psps <type(pseudopotential_type)>=variables related to pseudopotentials
 pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
 ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase information for given atom coordinates.
 qgrid(mqgrid)=q grid for spline atomic valence density n^AT(q) from 0 to qmax
 ucvol=unit cell volume (bohr**3).
 usepaw= 0 for non paw calculation; =1 for paw calculation
 vresid(nfft,nspden)=potential residual - (non-collinear magn. : only V11 and V22 are used)
 zion(ntypat)=charge on each type of atom (real number)
 znucl(ntypat)=atomic number, for each type of atom

OUTPUT

 gresid(3,natom)=forces due to the residual of the potential

SOURCE

751 subroutine fresidrsp(atindx1,dtset,gmet,gprimd,gresid,gsqcut,mgfft,mpi_enreg,mqgrid,nattyp,nfft,&
752 &          ngfft,ntypat,psps,pawtab,ph1d,qgrid,rprimd,ucvol,usepaw,vresid,zion,znucl)
753 
754 !Arguments ------------------------------------
755 !scalars
756  integer,intent(in) :: mgfft,mqgrid,nfft,ntypat,usepaw
757  real(dp),intent(in) :: gsqcut,ucvol
758  type(pseudopotential_type),intent(in) :: psps
759  type(MPI_type),intent(in) :: mpi_enreg
760  type(dataset_type),intent(in) :: dtset
761 !arrays
762  integer,intent(in) :: atindx1(dtset%natom),nattyp(ntypat),ngfft(18)
763  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom)
764  real(dp),intent(in) :: qgrid(mqgrid),vresid(nfft,dtset%nspden),zion(ntypat)
765  real(dp),intent(inout) :: rprimd(3,3)
766  real(dp),intent(in) :: znucl(ntypat)
767  real(dp),intent(out) :: gresid(3,dtset%natom)
768  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
769 
770 !Local variables-------------------------------
771 !scalars
772  integer :: itypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv
773  logical :: usegauss
774 !arrays
775  integer :: dummy3(3)
776  real(dp) :: dummy2(2)
777  real(dp) :: dummy_in1(0),dummy_in2(0)
778  real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0)
779  real(dp) :: strn_dummy6(6),strv_dummy6(6)
780  real(dp),allocatable :: gauss(:,:),vresg(:,:),work(:)
781 
782 ! *************************************************************************
783 
784 !Inits
785  optatm=0;optdyfr=0;opteltfr=0;optgr=1;optstr=0;optv=0;optn=1
786  ABI_MALLOC(vresg,(2,nfft))
787 
788 !Transfer potential residual to reciprocal space
789 !Use only Vres=Vres11+Vres22=Vres_up+Vres_dn
790  ABI_MALLOC(work,(nfft))
791  work(:)=vresid(:,1)
792  if (dtset%nspden>=2) work(:)=work(:)+vresid(:,2)
793  call fourdp(1,vresg,work,-1,mpi_enreg,nfft,1,ngfft,0)
794  ABI_FREE(work)
795 
796 !Determine wether a gaussan atomic density has to be used or not
797  usegauss=.true.
798  if (usepaw==0) usegauss = any(.not.psps%nctab(1:ntypat)%has_tvale)
799  if (usepaw==1) usegauss=(minval(pawtab(1:ntypat)%has_tvale)==0)
800  if (usegauss) then
801    optn2=3
802    ABI_MALLOC(gauss,(2,ntypat))
803    do itypat=1,ntypat
804      gauss(1,itypat)=zion(itypat)
805      gauss(2,itypat) = atom_length(dtset%densty(itypat,1),zion(itypat),znucl(itypat))
806    end do
807    call wrtout(std_out," Computing residual forces using gaussian functions as atomic densities", "COLL")
808  else
809    optn2=2
810    ABI_MALLOC(gauss,(2,0))
811    call wrtout(std_out," Computing residual forces using atomic densities taken from pseudos", "COLL")
812  end if
813 
814 !Compute forces due to residual
815  call atm2fft(atindx1,dummy_out1,dummy_out2,dummy_out3,dummy_out4,&
816 & dummy_out5,gauss,gmet,gprimd,gresid,dummy_out6,gsqcut,mgfft,&
817 & mqgrid,dtset%natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
818 & psps,pawtab,ph1d,qgrid,dummy3,dtset%rcut,dummy_in1,rprimd,strn_dummy6,strv_dummy6,ucvol,usepaw,&
819 & vresg,vresg,vresg,dummy2,dummy_in2,comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
820 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
821 
822 !In case of nspden>=2, has to apply 1/2 factor
823  if (dtset%nspden>=2) gresid=gresid*half
824 
825  ABI_FREE(gauss)
826  ABI_FREE(vresg)
827 
828 end subroutine fresidrsp

ABINIT/m_forces [ Modules ]

[ Top ] [ Modules ]

NAME

  m_forces

FUNCTION

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, FJ, MM, MT, SCE)
  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_forces
23 
24  use defs_basis
25  use defs_wvltypes
26  use m_abicore
27  use m_efield
28  use m_errors
29  use m_atomdata
30  use m_dtset
31 
32  use defs_datatypes,     only : pseudopotential_type
33  use defs_abitypes,      only : MPI_type
34  use m_time,             only : timab
35  use m_geometry,         only : gred2fcart, metric, xred2xcart
36  use m_fock,             only : fock_type
37  use m_pawrad,           only : pawrad_type
38  use m_pawtab,           only : pawtab_type
39  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
40  use libxc_functionals,  only : libxc_functionals_is_hybrid
41  use m_fft,              only : zerosym, fourdp
42  use m_cgtools,          only : mean_fftr
43  use m_mpinfo,           only : pre_gather, pre_scatter
44  use m_atm2fft,          only : atm2fft
45  use m_mklocl,           only : mklocl
46  use m_predtk,           only : prtxvf
47  use m_xchybrid,         only : xchybrid_ncpp_cc
48  use m_mkcore,           only : mkcore, mkcore_alt
49  use m_mkcore_wvl,       only : mkcore_wvl
50 
51  implicit none
52 
53  private

ABINIT/sygrad [ Functions ]

[ Top ] [ Functions ]

NAME

 sygrad

FUNCTION

 Symmetrize derivatives of energy with respect to coordinates.
 Unsymmetrized gradients are input as dedt; symmetrized grads are then placed in gred.
 If nsym=1 simply copy dedt into gred (only symmetry is identity).

INPUTS

  natom=number of atoms in cell
  dedt(3,natom)=unsymmetrized gradients wrt dimensionless tn (hartree)
  nsym=number of symmetry operators in group
  symrec(3,3,nsym)=symmetries of group in terms of operations on
    reciprocal space primitive translations--see comments below
  indsym(4,nsym,natom)=label given by subroutine symatm, indicating atom
   label which gets rotated into given atom by given symmetry
   (first three elements are related primitive translation--
   see symatm where this is computed)

OUTPUT

 gred(3,3,natom)=symmetrized gradients wrt reduced coordinates (hartree)

NOTES

 symmetrization of gradients with respect to reduced
 coordinates tn is conducted according to the expression
 $[d(e)/d(t(n,a))]_{symmetrized} = (1/Nsym)*Sum(S)*symrec(n,m,S)*
              [d(e)/d(t(m,b))]_{unsymmetrized}$
 where $t(m,b)= (symrel^{-1})(m,n)*(t(n,a)-tnons(n))$ and tnons
 is a possible nonsymmorphic translation.  The label "b" here
 refers to the atom which gets rotated into "a" under symmetry "S".
 symrel is the symmetry matrix in real space, which is the inverse
 transpose of symrec.  symrec is the symmetry matrix in reciprocal
 space.  $sym_{cartesian} = R * symrel * R^{-1} = G * symrec * G^{-1}$
 where the columns of R and G are the dimensional primitive translations
 in real and reciprocal space respectively.
 Note the use of "symrec" in the symmetrization expression above.

SOURCE

665 subroutine sygrad(gred,natom,dedt,nsym,symrec,indsym)
666 
667 !Arguments ------------------------------------
668 !scalars
669  integer,intent(in) :: natom,nsym
670 !arrays
671  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym)
672  real(dp),intent(in) :: dedt(3,natom)
673  real(dp),intent(out) :: gred(3,natom)
674 
675 !Local variables-------------------------------
676 !scalars
677  integer :: ia,ind,isym,mu
678  real(dp),parameter :: tol=1.0d-30
679  real(dp) :: summ
680 
681 ! *************************************************************************
682 !
683  if (nsym==1) then
684 !  only symmetry is identity so simply copy
685    gred(:,:)=dedt(:,:)
686  else
687 !  actually conduct symmetrization
688    do ia=1,natom
689      do mu=1,3
690        summ=0._dp
691        do isym=1,nsym
692          ind=indsym(4,isym,ia)
693          summ=summ+dble(symrec(mu,1,isym))*dedt(1,ind)+&
694 &         dble(symrec(mu,2,isym))*dedt(2,ind)+&
695 &         dble(symrec(mu,3,isym))*dedt(3,ind)
696        end do
697        gred(mu,ia)=summ/dble(nsym)
698        if(abs(gred(mu,ia))<tol)gred(mu,ia)=0.0_dp
699      end do
700    end do
701  end if
702 
703 end subroutine sygrad