diff --git a/src/pair_coul_slater_cut.cpp b/src/pair_coul_slater_cut.cpp deleted file mode 100644 index ba4da22c80..0000000000 --- a/src/pair_coul_slater_cut.cpp +++ /dev/null @@ -1,183 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include -#include -#include -#include -#include "pair_coul_slater_cut.h" -#include "atom.h" -#include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -PairCoulSlaterCut::PairCoulSlaterCut(LAMMPS *lmp) : PairCoulCut(lmp) {} - -/* ---------------------------------------------------------------------- */ - -void PairCoulSlaterCut::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; - double rsq,r2inv,r,rinv,forcecoul,factor_coul,bracket_term; - int *ilist,*jlist,*numneigh,**firstneigh; - - ecoul = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - double **x = atom->x; - double **f = atom->f; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - qtmp = q[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_coul = special_coul[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; - r = sqrt(rsq); - rinv = 1.0/r; - bracket_term = 1 - exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda))); - forcecoul = qqrd2e * scale[itype][jtype] * - qtmp*q[j] * bracket_term * rinv; - fpair = factor_coul*forcecoul * r2inv; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; - } - - if (eflag) ecoul = factor_coul * qqrd2e * - scale[itype][jtype] * qtmp*q[j] * rinv * - (1 - (1 + r/lamda)*exp(-2*r/lamda)); - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - 0.0,ecoul,fpair,delx,dely,delz); - } - } - } - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairCoulSlaterCut::settings(int narg, char **arg) -{ - if (narg != 2) error->all(FLERR,"Illegal pair_style command"); - - lamda = force->numeric(FLERR,arg[0]); - cut_global = force->numeric(FLERR,arg[1]); - - // reset cutoffs that have been explicitly set - - if (allocated) { - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - if (setflag[i][j]) cut[i][j] = cut_global; - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairCoulSlaterCut::write_restart_settings(FILE *fp) -{ - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&lamda,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairCoulSlaterCut::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&cut_global,sizeof(double),1,fp); - fread(&lamda,sizeof(double),1,fp); - fread(&offset_flag,sizeof(int),1,fp); - fread(&mix_flag,sizeof(int),1,fp); - } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&lamda,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); -} - -/* ---------------------------------------------------------------------- */ - -double PairCoulSlaterCut::single(int i, int j, int /*itype*/, int /*jtype*/, - double rsq, double factor_coul, double /*factor_lj*/, - double &fforce) -{ - double r2inv,r,rinv,forcecoul,phicoul,bracket_term; - - r2inv = 1.0/rsq; - r = sqrt(rsq); - rinv = 1.0/r; - bracket_term = 1 - exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda))); - forcecoul = force->qqrd2e * atom->q[i]*atom->q[j] * - bracket_term * rinv; - fforce = factor_coul*forcecoul * r2inv; - - phicoul = force->qqrd2e * atom->q[i]*atom->q[j] * rinv * (1 - (1 + r/lamda)*exp(-2*r/lamda)); - return factor_coul*phicoul; -}