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
  fred(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

PARENTS

      forces

CHILDREN

      dposv,prtxvf,wrtout,xred2xcart

SOURCE

1518 subroutine constrf(diffor,fcart,forold,fred,iatfix,ionmov,maxfor,natom,&
1519 & nconeq,prtvol,rprimd,wtatcon,xred)
1520 
1521  use m_linalg_interfaces
1522 
1523 !This section has been created automatically by the script Abilint (TD).
1524 !Do not modify the following lines by hand.
1525 #undef ABI_FUNC
1526 #define ABI_FUNC 'constrf'
1527 !End of the abilint section
1528 
1529  implicit none
1530 
1531 !Arguments ------------------------------------
1532 !scalars
1533  integer,intent(in) :: ionmov,natom,nconeq,prtvol
1534  real(dp),intent(out) :: diffor,maxfor
1535 !arrays
1536  integer,intent(in) :: iatfix(3,natom)
1537  real(dp),intent(in) :: rprimd(3,3),wtatcon(3,natom,nconeq)
1538  real(dp),intent(inout) :: fcart(3,natom),forold(3,natom),xred(3,natom)
1539  real(dp),intent(inout) :: fred(3,natom) !vz_i
1540 
1541 !Local variables -------------------------
1542 !scalars
1543  integer :: iatom,iconeq,iconeq1,iconeq2,index,info,mu,prtvel
1544  character(len=500) :: message
1545 !arrays
1546  real(dp),allocatable :: fcartvec(:),fpcart(:,:),fvector(:),vel_dummy(:,:)
1547  real(dp),allocatable :: wmatrix(:,:),wtatconvec(:,:),wtcoeffs(:),xcart(:,:)
1548 
1549 !************************************************************************
1550 
1551 !Allocate temporary variables
1552  ABI_ALLOCATE(fpcart,(3,natom))
1553  ABI_ALLOCATE(fcartvec,(3*natom))
1554  ABI_ALLOCATE(fvector,(nconeq))
1555  ABI_ALLOCATE(vel_dummy,(3,natom))
1556  ABI_ALLOCATE(wmatrix,(nconeq,nconeq))
1557  ABI_ALLOCATE(wtatconvec,(3*natom,nconeq))
1558  ABI_ALLOCATE(wtcoeffs,(nconeq))
1559  ABI_ALLOCATE(xcart,(3,natom))
1560 
1561 !If prtvol>10, output coordinates and forces prior to projecting
1562  if(prtvol>=10)then
1563    write(message,'(a)')' constrf - coordinates and forces prior to constraint projections:'
1564    call wrtout(std_out,message,'COLL')
1565    call xred2xcart(natom,rprimd,xcart,xred)
1566    prtvel=0
1567    call prtxvf(fcart,fred,iatfix,06,natom,prtvel,vel_dummy,xcart,xred)
1568  end if
1569 
1570 !Transfer fcart and wtatcon to flat vectors
1571  index=0
1572  do iatom=1,natom
1573    do mu=1,3
1574      index=index+1
1575      fcartvec(index)=fcart(mu,iatom)
1576      wtatconvec(index,:)=wtatcon(mu,iatom,:)
1577    end do
1578  end do
1579 
1580 !Compute a matrix (wmatrix) and vector (fvector) such that solving
1581 !the linear equations wmatrix*wcoeffs=fvector gives the coefficients
1582 !of wtatcon (wcoeffs) needed to compute the projected forces
1583  do iconeq2=1,nconeq
1584    fvector(iconeq2)=ddot(3*natom,fcartvec,1,wtatconvec(1,iconeq2),1)
1585    do iconeq1=1,nconeq
1586      wmatrix(iconeq1,iconeq2)=ddot(3*natom,wtatconvec(1,iconeq1),1,wtatconvec(1,iconeq2),1)
1587    end do
1588  end do
1589 
1590 !Solve the system of linear equations, wmatrix*wcoeffs=fvector
1591  call dposv('U',nconeq,1,wmatrix,nconeq,fvector,nconeq,info)
1592 
1593  if (info/=0) then
1594    write(message, '(a,a,a,a,a)' )&
1595 &   'Constraint matrix is not positive definite,',ch10,&
1596 &   'probably because constraints are linearly dependent.',ch10,&
1597 &   'Action: Check for linear dependence of constraints.'
1598    MSG_ERROR(message)
1599  end if
1600 
1601 !The solution vector is returned in fvector, so copy it to a more sensible location
1602  wtcoeffs(:)=fvector(:)
1603 
1604 !Compute the projected forces, which now satisfy all the constraint equations
1605  fpcart(:,:)=fcart(:,:)
1606  do iconeq=1,nconeq
1607    fpcart(:,:)=fpcart(:,:)-wtcoeffs(iconeq)*wtatcon(:,:,iconeq)
1608  end do
1609 
1610 !Reconvert constrained forces back from fpcart to fred
1611  do iatom=1,natom
1612    do mu=1,3
1613      fred(mu,iatom)= - (rprimd(1,mu)*fpcart(1,iatom)+&
1614 &     rprimd(2,mu)*fpcart(2,iatom)+&
1615 &     rprimd(3,mu)*fpcart(3,iatom))
1616    end do
1617  end do
1618 
1619 !If prtvol>=10, output coordinates and forces after projecting
1620  if(prtvol>=10)then
1621    write(message,'(a)')' constrf - coordinates and forces after constraint projections:'
1622    call wrtout(std_out,message,'COLL')
1623    prtvel=0
1624    call prtxvf(fpcart,fred,iatfix,06,natom,prtvel,vel_dummy,xcart,xred)
1625  end if
1626 
1627 !Copy the constrained forces, fpcart, back to fcart
1628  fcart(:,:)=fpcart(:,:)
1629 
1630 !Compute maximal force and maximal difference of the projected forces,
1631 !overriding the values already computed in forces
1632  maxfor=0.0_dp
1633  diffor=0.0_dp
1634  do iatom=1,natom
1635    do mu=1,3
1636      if (iatfix(mu,iatom) /= 1) then
1637        maxfor=max(maxfor,abs(fcart(mu,iatom)))
1638        diffor=max(diffor,abs(fcart(mu,iatom)-forold(mu,iatom)))
1639      else if (ionmov==4 .or. ionmov==5) then
1640 !      Make the force vanish on fixed atoms when ionmov=4 or 5
1641 !      This is because fixing of atom cannot be imposed at the
1642 !      level of a routine similar to brdmin or moldyn for these options.
1643        fcart(mu,iatom)=0.0_dp
1644      end if
1645    end do
1646  end do
1647  forold(:,:)=fcart(:,:)
1648 
1649 !Dellocate temporary variables
1650  ABI_DEALLOCATE(fpcart)
1651  ABI_DEALLOCATE(fcartvec)
1652  ABI_DEALLOCATE(fvector)
1653  ABI_DEALLOCATE(vel_dummy)
1654  ABI_DEALLOCATE(wmatrix)
1655  ABI_DEALLOCATE(wtatconvec)
1656  ABI_DEALLOCATE(wtcoeffs)
1657  ABI_DEALLOCATE(xcart)
1658 
1659 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)
  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.
  vresid(nfft,nspden)=potential residual (if non-collinear magn., only trace of it)
  vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space
  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)
  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 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 fred, 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.

PARENTS

      etotfor,forstr

CHILDREN

      atm2fft,constrf,dgemv,fourdp,fred2fcart,fresid,fresidrsp,metric,mkcore
      mkcore_alt,mkcore_wvl,mklocl,sygrad,timab,xchybrid_ncpp_cc,zerosym

SOURCE

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

PARENTS

      forces,prcref,prcref_PMA,scfcv

CHILDREN

      atomdata_from_znucl,mean_fftr,pre_gather,pre_scatter

SOURCE

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

PARENTS

      forces

CHILDREN

      atm2fft,fourdp,wrtout

SOURCE

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

ABINIT/m_forces [ Modules ]

[ Top ] [ Modules ]

NAME

  m_forces

FUNCTION

COPYRIGHT

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

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_forces
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_atomdata
37 
38  use m_time,             only : timab
39  use m_geometry,         only : fred2fcart, metric, xred2xcart
40  use m_fock,             only : fock_type
41  use m_pawrad,           only : pawrad_type
42  use m_pawtab,           only : pawtab_type
43  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
44  use libxc_functionals,  only : libxc_functionals_is_hybrid
45  use m_fft,              only : zerosym, fourdp
46  use m_cgtools,          only : mean_fftr
47  use m_mpinfo,           only : pre_gather, pre_scatter
48  use m_atm2fft,          only : atm2fft
49  use m_mklocl,           only : mklocl
50  use m_predtk,           only : prtxvf
51  use m_xchybrid,         only : xchybrid_ncpp_cc
52  use m_mkcore,           only : mkcore, mkcore_alt
53  use m_mkcore_wvl,       only : mkcore_wvl
54 
55  implicit none
56 
57  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 fred.
 If nsym=1 simply copy dedt into fred (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

 fred(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.

PARENTS

      forces

CHILDREN

SOURCE

638 subroutine sygrad(fred,natom,dedt,nsym,symrec,indsym)
639 
640 
641 !This section has been created automatically by the script Abilint (TD).
642 !Do not modify the following lines by hand.
643 #undef ABI_FUNC
644 #define ABI_FUNC 'sygrad'
645 !End of the abilint section
646 
647  implicit none
648 
649 !Arguments ------------------------------------
650 !scalars
651  integer,intent(in) :: natom,nsym
652 !arrays
653  integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym)
654  real(dp),intent(in) :: dedt(3,natom)
655  real(dp),intent(out) :: fred(3,natom)
656 
657 !Local variables-------------------------------
658 !scalars
659  integer :: ia,ind,isym,mu
660  real(dp),parameter :: tol=1.0d-30
661  real(dp) :: summ
662 
663 ! *************************************************************************
664 !
665  if (nsym==1) then
666 !  only symmetry is identity so simply copy
667    fred(:,:)=dedt(:,:)
668  else
669 !  actually conduct symmetrization
670    do ia=1,natom
671      do mu=1,3
672        summ=0._dp
673        do isym=1,nsym
674          ind=indsym(4,isym,ia)
675          summ=summ+dble(symrec(mu,1,isym))*dedt(1,ind)+&
676 &         dble(symrec(mu,2,isym))*dedt(2,ind)+&
677 &         dble(symrec(mu,3,isym))*dedt(3,ind)
678        end do
679        fred(mu,ia)=summ/dble(nsym)
680        if(abs(fred(mu,ia))<tol)fred(mu,ia)=0.0_dp
681      end do
682    end do
683  end if
684 
685 end subroutine sygrad