TABLE OF CONTENTS
ABINIT/constrf [ 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 ]
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 ]
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 ]
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 ]
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 ]
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