diff --git a/src/USER-FEP/pair_lj_class2_coul_cut_soft.cpp b/src/USER-FEP/pair_lj_class2_coul_cut_soft.cpp
index 4ce19aef6e41e9b61a9a39abc8afdc1a4c840705..7970097cfec1bb609ebe6fd4d8bb92079426f1d2 100644
--- a/src/USER-FEP/pair_lj_class2_coul_cut_soft.cpp
+++ b/src/USER-FEP/pair_lj_class2_coul_cut_soft.cpp
@@ -68,7 +68,7 @@ void PairLJClass2CoulCutSoft::compute(int eflag, int vflag)
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double rsq,forcecoul,forcelj;
double factor_coul,factor_lj;
- double denc, denlj, r4sig6;
+ double denc, denlj, r4sig6;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
@@ -124,8 +124,8 @@ void PairLJClass2CoulCutSoft::compute(int eflag, int vflag)
if (rsq < cut_ljsq[itype][jtype]) {
denlj = lj3[itype][jtype] + pow(rsq, 3) * pow(sigma[itype][jtype], -6.0);
r4sig6 = rsq*rsq / lj2[itype][jtype];
- forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
- (18.0*r4sig6*pow(denlj, -2.5) - 18.0*r4sig6*pow(denlj, -2));
+ forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
+ (18.0*r4sig6*pow(denlj, -2.5) - 18.0*r4sig6*pow(denlj, -2));
} else forcelj = 0.0;
fpair = factor_coul*forcecoul + factor_lj*forcelj;
@@ -477,7 +477,7 @@ double PairLJClass2CoulCutSoft::single(int i, int j, int itype, int jtype,
{
double forcecoul,forcelj,phicoul,philj;
double denc, denlj, r4sig6;
-
+
if (rsq < cut_coulsq[itype][jtype]) {
denc = sqrt(lj4[itype][jtype] + rsq);
forcecoul = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] /
@@ -486,7 +486,7 @@ double PairLJClass2CoulCutSoft::single(int i, int j, int itype, int jtype,
if (rsq < cut_ljsq[itype][jtype]) {
r4sig6 = rsq*rsq / lj2[itype][jtype];
denlj = lj3[itype][jtype] + rsq*r4sig6;
- forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
+ forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
(18.0*r4sig6/(denlj*denlj*sqrt(denlj)) - 18.0*r4sig6/(denlj*denlj));
} else forcelj = 0.0;
fforce = factor_coul*forcecoul + factor_lj*forcelj;
@@ -517,4 +517,3 @@ void *PairLJClass2CoulCutSoft::extract(const char *str, int &dim)
if (strcmp(str,"lambda") == 0) return (void *) lambda;
return NULL;
}
-
diff --git a/src/USER-FEP/pair_lj_class2_soft.cpp b/src/USER-FEP/pair_lj_class2_soft.cpp
index 523123d9831c789a8b76c88767b554dd8f5c3e20..b7f21fc59c67dc602264f891a2eb61dd71f1ded2 100644
--- a/src/USER-FEP/pair_lj_class2_soft.cpp
+++ b/src/USER-FEP/pair_lj_class2_soft.cpp
@@ -1,438 +1,438 @@
-/* ----------------------------------------------------------------------
- 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 <cmath>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-#include "pair_lj_class2_soft.h"
-#include "atom.h"
-#include "comm.h"
-#include "force.h"
-#include "neigh_list.h"
-#include "math_const.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-using namespace MathConst;
-
-/* ---------------------------------------------------------------------- */
-
-PairLJClass2Soft::PairLJClass2Soft(LAMMPS *lmp) : Pair(lmp)
-{
- writedata = 1;
- allocated = 0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-PairLJClass2Soft::~PairLJClass2Soft()
-{
- if (!copymode) {
- if (allocated) {
- memory->destroy(setflag);
- memory->destroy(cutsq);
-
- memory->destroy(cut);
- memory->destroy(epsilon);
- memory->destroy(sigma);
- memory->destroy(lambda);
- memory->destroy(lj1);
- memory->destroy(lj2);
- memory->destroy(lj3);
- memory->destroy(offset);
- allocated=0;
- }
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairLJClass2Soft::compute(int eflag, int vflag)
-{
- int i,j,ii,jj,inum,jnum,itype,jtype;
- double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
- double rsq,forcelj,factor_lj;
- double denlj, r4sig6;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- evdwl = 0.0;
- if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = vflag_fdotr = 0;
-
- double **x = atom->x;
- double **f = atom->f;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- double *special_lj = force->special_lj;
- int newton_pair = force->newton_pair;
-
- 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];
- 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_lj = special_lj[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]) {
- denlj = lj3[itype][jtype] + pow(rsq, 3) * pow(sigma[itype][jtype], -6.0);
- r4sig6 = rsq*rsq / lj2[itype][jtype];
- forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
- (18.0*r4sig6*pow(denlj, -2.5) - 18.0*r4sig6*pow(denlj, -2));
- fpair = factor_lj*forcelj;
-
- 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) {
- denlj = lj3[itype][jtype] + pow(rsq, 3) * pow(sigma[itype][jtype], -6.0);
- evdwl = lj1[itype][jtype] * epsilon[itype][jtype] * (2.0/(denlj*sqrt(denlj)) - 3.0/denlj) -
- offset[itype][jtype];
- evdwl *= factor_lj;
- }
-
- if (evflag) ev_tally(i,j,nlocal,newton_pair,
- evdwl,0.0,fpair,delx,dely,delz);
- }
- }
- }
-
- if (vflag_fdotr) virial_fdotr_compute();
-}
-
-/* ----------------------------------------------------------------------
- allocate all arrays
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::allocate()
-{
- allocated = 1;
- int n = atom->ntypes;
-
- memory->create(setflag,n+1,n+1,"pair:setflag");
- for (int i = 1; i <= n; i++)
- for (int j = i; j <= n; j++)
- setflag[i][j] = 0;
-
- memory->create(cutsq,n+1,n+1,"pair:cutsq");
-
- memory->create(cut,n+1,n+1,"pair:cut");
- memory->create(epsilon,n+1,n+1,"pair:epsilon");
- memory->create(sigma,n+1,n+1,"pair:sigma");
- memory->create(lambda,n+1,n+1,"pair:lambda");
- memory->create(lj1,n+1,n+1,"pair:lj1");
- memory->create(lj2,n+1,n+1,"pair:lj2");
- memory->create(lj3,n+1,n+1,"pair:lj3");
- memory->create(offset,n+1,n+1,"pair:offset");
-}
-
-/* ----------------------------------------------------------------------
- global settings
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::settings(int narg, char **arg)
-{
- if (narg != 3) error->all(FLERR,"Illegal pair_style command");
-
- nlambda = force->numeric(FLERR,arg[0]);
- alphalj = force->numeric(FLERR,arg[1]);
-
- cut_global = force->numeric(FLERR,arg[2]);
-
- // 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;
- }
-}
-
-/* ----------------------------------------------------------------------
- set coeffs for one or more type pairs
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::coeff(int narg, char **arg)
-{
- if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
- if (!allocated) allocate();
-
- int ilo,ihi,jlo,jhi;
- force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
- force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
-
- double epsilon_one = force->numeric(FLERR,arg[2]);
- double sigma_one = force->numeric(FLERR,arg[3]);
- double lambda_one = force->numeric(FLERR,arg[4]);
- if (sigma_one <= 0.0) error->all(FLERR,"Incorrect args for pair coefficients");
-
- double cut_one = cut_global;
- if (narg == 6) cut_one = force->numeric(FLERR,arg[5]);
-
- int count = 0;
- for (int i = ilo; i <= ihi; i++) {
- for (int j = MAX(jlo,i); j <= jhi; j++) {
- epsilon[i][j] = epsilon_one;
- sigma[i][j] = sigma_one;
- lambda[i][j] = lambda_one;
- cut[i][j] = cut_one;
- setflag[i][j] = 1;
- count++;
- }
- }
-
- if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
-}
-
-/* ----------------------------------------------------------------------
- init for one type pair i,j and corresponding j,i
-------------------------------------------------------------------------- */
-
-double PairLJClass2Soft::init_one(int i, int j)
-{
- // always mix epsilon,sigma via sixthpower rules
- // mix distance via user-defined rule
-
- if (setflag[i][j] == 0) {
- epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) *
- pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) /
- (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0));
- sigma[i][j] =
- pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0);
- if (lambda[i][i] != lambda[j][j])
- error->all(FLERR,"Pair lj/class2/coul/cut/soft different lambda values in mix");
- lambda[i][j] = lambda[i][i];
- cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
- }
-
- lj1[i][j] = pow(lambda[i][j], nlambda);
- lj2[i][j] = pow(sigma[i][j], 6.0);
- lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
-
- if (offset_flag && (cut[i][j] > 0.0)) {
- double denlj = lj3[i][j] + pow(cut[i][j] / sigma[i][j], 6.0);
- offset[i][j] = lj1[i][j] * epsilon[i][j] * (2.0/(denlj*sqrt(denlj)) - 3.0/denlj);
- } else offset[i][j] = 0.0;
-
- epsilon[j][i] = epsilon[i][j];
- sigma[j][i] = sigma[i][j];
- lambda[j][i] = lambda[i][j];
- lj1[j][i] = lj1[i][j];
- lj2[j][i] = lj2[i][j];
- lj3[j][i] = lj3[i][j];
- offset[j][i] = offset[i][j];
-
- // compute I,J contribution to long-range tail correction
- // count total # of atoms of type I and J via Allreduce
-
- if (tail_flag) {
- int *type = atom->type;
- int nlocal = atom->nlocal;
-
- double count[2],all[2];
- count[0] = count[1] = 0.0;
- for (int k = 0; k < nlocal; k++) {
- if (type[k] == i) count[0] += 1.0;
- if (type[k] == j) count[1] += 1.0;
- }
- MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
-
- double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
- double sig6 = sig3*sig3;
- double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
- double rc6 = rc3*rc3;
- // check the following expressions for etail_lj & ptail_ij they are not correct !
- etail_ij = 2.0*MY_PI*all[0]*all[1]*lj1[i][j] *epsilon[i][j] *
- sig6 * (sig3 - 3.0*rc3) / (3.0*rc6);
- ptail_ij = 2.0*MY_PI*all[0]*all[1]*lj1[i][j] *epsilon[i][j] *
- sig6 * (sig3 - 2.0*rc3) / rc6;
- }
-
- return cut[i][j];
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes to restart file
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::write_restart(FILE *fp)
-{
- write_restart_settings(fp);
-
- int i,j;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i; j <= atom->ntypes; j++) {
- fwrite(&setflag[i][j],sizeof(int),1,fp);
- if (setflag[i][j]) {
- fwrite(&epsilon[i][j],sizeof(double),1,fp);
- fwrite(&sigma[i][j],sizeof(double),1,fp);
- fwrite(&lambda[i][j],sizeof(double),1,fp);
- fwrite(&cut[i][j],sizeof(double),1,fp);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- proc 0 reads from restart file, bcasts
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::read_restart(FILE *fp)
-{
- read_restart_settings(fp);
- allocate();
-
- int i,j;
- int me = comm->me;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i; j <= atom->ntypes; j++) {
- if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
- MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
- if (setflag[i][j]) {
- if (me == 0) {
- fread(&epsilon[i][j],sizeof(double),1,fp);
- fread(&sigma[i][j],sizeof(double),1,fp);
- fread(&lambda[i][j],sizeof(double),1,fp);
- fread(&cut[i][j],sizeof(double),1,fp);
- }
- MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes to restart file
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::write_restart_settings(FILE *fp)
-{
- fwrite(&nlambda,sizeof(double),1,fp);
- fwrite(&alphalj,sizeof(double),1,fp);
- fwrite(&cut_global,sizeof(double),1,fp);
- fwrite(&offset_flag,sizeof(int),1,fp);
- fwrite(&mix_flag,sizeof(int),1,fp);
- fwrite(&tail_flag,sizeof(int),1,fp);
-}
-
-/* ----------------------------------------------------------------------
- proc 0 reads from restart file, bcasts
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::read_restart_settings(FILE *fp)
-{
- int me = comm->me;
- if (me == 0) {
- fread(&nlambda,sizeof(double),1,fp);
- fread(&alphalj,sizeof(double),1,fp);
- fread(&cut_global,sizeof(double),1,fp);
- fread(&offset_flag,sizeof(int),1,fp);
- fread(&mix_flag,sizeof(int),1,fp);
- fread(&tail_flag,sizeof(int),1,fp);
- }
- MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&alphalj,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
- MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
- MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
-}
-
-
-/* ----------------------------------------------------------------------
- proc 0 writes to data file
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::write_data(FILE *fp)
-{
- for (int i = 1; i <= atom->ntypes; i++)
- fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],lambda[i][i]);
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes all pairs to data file
-------------------------------------------------------------------------- */
-
-void PairLJClass2Soft::write_data_all(FILE *fp)
-{
- for (int i = 1; i <= atom->ntypes; i++)
- for (int j = i; j <= atom->ntypes; j++)
- fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],
- lambda[i][j],cut[i][j]);
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairLJClass2Soft::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
- double /*factor_coul*/, double factor_lj,
- double &fforce)
-{
- double forcelj,philj;
- double r4sig6, denlj;
-
- if (rsq < cutsq[itype][jtype]) {
- r4sig6 = rsq*rsq / lj2[itype][jtype];
- denlj = lj3[itype][jtype] + rsq*r4sig6;
- forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
- (18.0*r4sig6/(denlj*denlj*sqrt(denlj)) - 18.0*r4sig6/(denlj*denlj));
- } else forcelj = 0.0;
- fforce = factor_lj*forcelj;
-
- if (rsq < cutsq[itype][jtype]) {
- denlj = lj3[itype][jtype] + rsq*rsq*rsq / lj2[itype][jtype];
- philj = lj1[itype][jtype] * epsilon[itype][jtype] * (2.0/(denlj*sqrt(denlj)) - 3.0/denlj) -
- offset[itype][jtype];
- } else philj = 0.0;
-
- return factor_lj*philj;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void *PairLJClass2Soft::extract(const char *str, int &dim)
-{
- dim = 2;
- if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
- if (strcmp(str,"sigma") == 0) return (void *) sigma;
- if (strcmp(str,"lambda") == 0) return (void *) lambda;
- return NULL;
-}
+/* ----------------------------------------------------------------------
+ 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 <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include "pair_lj_class2_soft.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neigh_list.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairLJClass2Soft::PairLJClass2Soft(LAMMPS *lmp) : Pair(lmp)
+{
+ writedata = 1;
+ allocated = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairLJClass2Soft::~PairLJClass2Soft()
+{
+ if (!copymode) {
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(cut);
+ memory->destroy(epsilon);
+ memory->destroy(sigma);
+ memory->destroy(lambda);
+ memory->destroy(lj1);
+ memory->destroy(lj2);
+ memory->destroy(lj3);
+ memory->destroy(offset);
+ allocated=0;
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJClass2Soft::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+ double rsq,forcelj,factor_lj;
+ double denlj, r4sig6;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+
+ 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];
+ 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_lj = special_lj[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]) {
+ denlj = lj3[itype][jtype] + pow(rsq, 3) * pow(sigma[itype][jtype], -6.0);
+ r4sig6 = rsq*rsq / lj2[itype][jtype];
+ forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
+ (18.0*r4sig6*pow(denlj, -2.5) - 18.0*r4sig6*pow(denlj, -2));
+ fpair = factor_lj*forcelj;
+
+ 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) {
+ denlj = lj3[itype][jtype] + pow(rsq, 3) * pow(sigma[itype][jtype], -6.0);
+ evdwl = lj1[itype][jtype] * epsilon[itype][jtype] * (2.0/(denlj*sqrt(denlj)) - 3.0/denlj) -
+ offset[itype][jtype];
+ evdwl *= factor_lj;
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,0.0,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(cut,n+1,n+1,"pair:cut");
+ memory->create(epsilon,n+1,n+1,"pair:epsilon");
+ memory->create(sigma,n+1,n+1,"pair:sigma");
+ memory->create(lambda,n+1,n+1,"pair:lambda");
+ memory->create(lj1,n+1,n+1,"pair:lj1");
+ memory->create(lj2,n+1,n+1,"pair:lj2");
+ memory->create(lj3,n+1,n+1,"pair:lj3");
+ memory->create(offset,n+1,n+1,"pair:offset");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::settings(int narg, char **arg)
+{
+ if (narg != 3) error->all(FLERR,"Illegal pair_style command");
+
+ nlambda = force->numeric(FLERR,arg[0]);
+ alphalj = force->numeric(FLERR,arg[1]);
+
+ cut_global = force->numeric(FLERR,arg[2]);
+
+ // 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;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::coeff(int narg, char **arg)
+{
+ if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+ double epsilon_one = force->numeric(FLERR,arg[2]);
+ double sigma_one = force->numeric(FLERR,arg[3]);
+ double lambda_one = force->numeric(FLERR,arg[4]);
+ if (sigma_one <= 0.0) error->all(FLERR,"Incorrect args for pair coefficients");
+
+ double cut_one = cut_global;
+ if (narg == 6) cut_one = force->numeric(FLERR,arg[5]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ epsilon[i][j] = epsilon_one;
+ sigma[i][j] = sigma_one;
+ lambda[i][j] = lambda_one;
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairLJClass2Soft::init_one(int i, int j)
+{
+ // always mix epsilon,sigma via sixthpower rules
+ // mix distance via user-defined rule
+
+ if (setflag[i][j] == 0) {
+ epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) *
+ pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) /
+ (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0));
+ sigma[i][j] =
+ pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0);
+ if (lambda[i][i] != lambda[j][j])
+ error->all(FLERR,"Pair lj/class2/coul/cut/soft different lambda values in mix");
+ lambda[i][j] = lambda[i][i];
+ cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
+ }
+
+ lj1[i][j] = pow(lambda[i][j], nlambda);
+ lj2[i][j] = pow(sigma[i][j], 6.0);
+ lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
+
+ if (offset_flag && (cut[i][j] > 0.0)) {
+ double denlj = lj3[i][j] + pow(cut[i][j] / sigma[i][j], 6.0);
+ offset[i][j] = lj1[i][j] * epsilon[i][j] * (2.0/(denlj*sqrt(denlj)) - 3.0/denlj);
+ } else offset[i][j] = 0.0;
+
+ epsilon[j][i] = epsilon[i][j];
+ sigma[j][i] = sigma[i][j];
+ lambda[j][i] = lambda[i][j];
+ lj1[j][i] = lj1[i][j];
+ lj2[j][i] = lj2[i][j];
+ lj3[j][i] = lj3[i][j];
+ offset[j][i] = offset[i][j];
+
+ // compute I,J contribution to long-range tail correction
+ // count total # of atoms of type I and J via Allreduce
+
+ if (tail_flag) {
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ double count[2],all[2];
+ count[0] = count[1] = 0.0;
+ for (int k = 0; k < nlocal; k++) {
+ if (type[k] == i) count[0] += 1.0;
+ if (type[k] == j) count[1] += 1.0;
+ }
+ MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
+
+ double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
+ double sig6 = sig3*sig3;
+ double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
+ double rc6 = rc3*rc3;
+ // check the following expressions for etail_lj & ptail_ij they are not correct !
+ etail_ij = 2.0*MY_PI*all[0]*all[1]*lj1[i][j] *epsilon[i][j] *
+ sig6 * (sig3 - 3.0*rc3) / (3.0*rc6);
+ ptail_ij = 2.0*MY_PI*all[0]*all[1]*lj1[i][j] *epsilon[i][j] *
+ sig6 * (sig3 - 2.0*rc3) / rc6;
+ }
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::write_restart(FILE *fp)
+{
+ write_restart_settings(fp);
+
+ int i,j;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ fwrite(&setflag[i][j],sizeof(int),1,fp);
+ if (setflag[i][j]) {
+ fwrite(&epsilon[i][j],sizeof(double),1,fp);
+ fwrite(&sigma[i][j],sizeof(double),1,fp);
+ fwrite(&lambda[i][j],sizeof(double),1,fp);
+ fwrite(&cut[i][j],sizeof(double),1,fp);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::read_restart(FILE *fp)
+{
+ read_restart_settings(fp);
+ allocate();
+
+ int i,j;
+ int me = comm->me;
+ for (i = 1; i <= atom->ntypes; i++)
+ for (j = i; j <= atom->ntypes; j++) {
+ if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+ MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+ if (setflag[i][j]) {
+ if (me == 0) {
+ fread(&epsilon[i][j],sizeof(double),1,fp);
+ fread(&sigma[i][j],sizeof(double),1,fp);
+ fread(&lambda[i][j],sizeof(double),1,fp);
+ fread(&cut[i][j],sizeof(double),1,fp);
+ }
+ MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::write_restart_settings(FILE *fp)
+{
+ fwrite(&nlambda,sizeof(double),1,fp);
+ fwrite(&alphalj,sizeof(double),1,fp);
+ fwrite(&cut_global,sizeof(double),1,fp);
+ fwrite(&offset_flag,sizeof(int),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+ fwrite(&tail_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::read_restart_settings(FILE *fp)
+{
+ int me = comm->me;
+ if (me == 0) {
+ fread(&nlambda,sizeof(double),1,fp);
+ fread(&alphalj,sizeof(double),1,fp);
+ fread(&cut_global,sizeof(double),1,fp);
+ fread(&offset_flag,sizeof(int),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ fread(&tail_flag,sizeof(int),1,fp);
+ }
+ MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&alphalj,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+ MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
+}
+
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::write_data(FILE *fp)
+{
+ for (int i = 1; i <= atom->ntypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],lambda[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairLJClass2Soft::write_data_all(FILE *fp)
+{
+ for (int i = 1; i <= atom->ntypes; i++)
+ for (int j = i; j <= atom->ntypes; j++)
+ fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],
+ lambda[i][j],cut[i][j]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJClass2Soft::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
+ double /*factor_coul*/, double factor_lj,
+ double &fforce)
+{
+ double forcelj,philj;
+ double r4sig6, denlj;
+
+ if (rsq < cutsq[itype][jtype]) {
+ r4sig6 = rsq*rsq / lj2[itype][jtype];
+ denlj = lj3[itype][jtype] + rsq*r4sig6;
+ forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
+ (18.0*r4sig6/(denlj*denlj*sqrt(denlj)) - 18.0*r4sig6/(denlj*denlj));
+ } else forcelj = 0.0;
+ fforce = factor_lj*forcelj;
+
+ if (rsq < cutsq[itype][jtype]) {
+ denlj = lj3[itype][jtype] + rsq*rsq*rsq / lj2[itype][jtype];
+ philj = lj1[itype][jtype] * epsilon[itype][jtype] * (2.0/(denlj*sqrt(denlj)) - 3.0/denlj) -
+ offset[itype][jtype];
+ } else philj = 0.0;
+
+ return factor_lj*philj;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairLJClass2Soft::extract(const char *str, int &dim)
+{
+ dim = 2;
+ if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
+ if (strcmp(str,"sigma") == 0) return (void *) sigma;
+ if (strcmp(str,"lambda") == 0) return (void *) lambda;
+ return NULL;
+}
diff --git a/src/USER-FEP/pair_lj_class2_soft.h b/src/USER-FEP/pair_lj_class2_soft.h
index 24cb2527d545f0c44486683b174219fa887fe6a1..ef031a1911a2d7ebbfd33a831d4e8a57671f6fec 100644
--- a/src/USER-FEP/pair_lj_class2_soft.h
+++ b/src/USER-FEP/pair_lj_class2_soft.h
@@ -1,69 +1,69 @@
-/* -*- c++ -*- ----------------------------------------------------------
- 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.
-------------------------------------------------------------------------- */
-
-#ifdef PAIR_CLASS
-
-PairStyle(lj/class2/soft,PairLJClass2Soft)
-
-#else
-
-#ifndef LMP_PAIR_LJ_CLASS2_SOFT_H
-#define LMP_PAIR_LJ_CLASS2_SOFT_H
-
-#include "pair.h"
-
-namespace LAMMPS_NS {
-
-class PairLJClass2Soft : public Pair {
- public:
- PairLJClass2Soft(class LAMMPS *);
- virtual ~PairLJClass2Soft();
- virtual void compute(int, int);
- virtual void settings(int, char **);
- void coeff(int, char **);
- virtual double init_one(int, int);
- void write_restart(FILE *);
- void read_restart(FILE *);
- virtual void write_restart_settings(FILE *);
- virtual void read_restart_settings(FILE *);
- void write_data(FILE *);
- void write_data_all(FILE *);
- virtual double single(int, int, int, int, double, double, double, double &);
- void *extract(const char *, int &);
-
- protected:
- double cut_global;
- double **cut;
- double **epsilon, **sigma, **lambda;
- double nlambda, alphalj;
- double **lj1,**lj2,**lj3,**offset;
-
- virtual void allocate();
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-E: Illegal ... command
-Self-explanatory. Check the input script syntax and compare to the
-documentation for the command. You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-E: Incorrect args for pair coefficients
-Self-explanatory. Check the input script or data file.
-
-E: Pair lj/cut/soft different lambda values in mix
-The value of lambda has to be the same for I J pairs.
-
-*/
\ No newline at end of file
+/* -*- c++ -*- ----------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(lj/class2/soft,PairLJClass2Soft)
+
+#else
+
+#ifndef LMP_PAIR_LJ_CLASS2_SOFT_H
+#define LMP_PAIR_LJ_CLASS2_SOFT_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairLJClass2Soft : public Pair {
+ public:
+ PairLJClass2Soft(class LAMMPS *);
+ virtual ~PairLJClass2Soft();
+ virtual void compute(int, int);
+ virtual void settings(int, char **);
+ void coeff(int, char **);
+ virtual double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ virtual void write_restart_settings(FILE *);
+ virtual void read_restart_settings(FILE *);
+ void write_data(FILE *);
+ void write_data_all(FILE *);
+ virtual double single(int, int, int, int, double, double, double, double &);
+ void *extract(const char *, int &);
+
+ protected:
+ double cut_global;
+ double **cut;
+ double **epsilon, **sigma, **lambda;
+ double nlambda, alphalj;
+ double **lj1,**lj2,**lj3,**offset;
+
+ virtual void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+E: Illegal ... command
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+Self-explanatory. Check the input script or data file.
+
+E: Pair lj/cut/soft different lambda values in mix
+The value of lambda has to be the same for I J pairs.
+
+*/