diff --git a/src/USER-MISC/README b/src/USER-MISC/README index c8ca571b47..9f98bc6a02 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -48,12 +48,16 @@ fix ttm/mod, Sergey Starikov and Vasily Pisarev (JIHT), pisarevvv at gmail.com, improper_style cossq, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12 improper_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 improper_style ring, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12 +improper_style distance, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15 +pair_style buck/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15 pair_style coul/diel, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11 pair_style dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11 pair_style edip, Luca Ferraro, luca.ferraro at caspur.it, 15 Sep 11 pair_style eam/cd, Alexander Stukowski, stukowski at mm.tu-darmstadt.de, 7 Nov 09 pair_style gauss/cut, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11 +pair_style lennard/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15 pair_style list, Axel Kohlmeyer (Temple U), akohlmey at gmail.com, 1 Jun 13 +pair_style lj/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15 pair_style lj/sf, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12 pair_style meam/sw/spline, Robert Rudd (LLNL), robert.rudd at llnl.gov, 1 Oct 12 diff --git a/src/USER-MISC/improper_distance.cpp b/src/USER-MISC/improper_distance.cpp new file mode 100644 index 0000000000..a98f404d52 --- /dev/null +++ b/src/USER-MISC/improper_distance.cpp @@ -0,0 +1,262 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Paolo Raiteri (Curtin University) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "improper_distance.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "domain.h" +#include "force.h" +#include "update.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define TOLERANCE 0.05 +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +ImproperDistance::ImproperDistance(LAMMPS *lmp) : Improper(lmp) {} + +/* ---------------------------------------------------------------------- */ + +ImproperDistance::~ImproperDistance() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(chi); + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperDistance::compute(int eflag, int vflag) +{ + int i1,i2,i3,i4,n,type; + double xab, yab, zab; // bond 1-2 + double xac, yac, zac; // bond 1-3 + double xad, yad, zad; // bond 1-4 + double xbc, ybc, zbc; // bond 2-3 + double xbd, ybd, zbd; // bond 2-4 + double xna, yna, zna, rna; // normal + double da; + + double eimproper,f1[3],f2[3],f3[3],f4[3]; +// double ss1,ss2,ss3,r1,r2,r3,c0,c1,c2,s1,s2; +// double s12,c,s,domega,a,a11,a22,a33,a12,a13,a23; + double domega,a; + double sx2,sy2,sz2; + + eimproper = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **improperlist = neighbor->improperlist; + int nimproperlist = neighbor->nimproperlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nimproperlist; n++) { + i1 = improperlist[n][0]; + i2 = improperlist[n][1]; + i3 = improperlist[n][2]; + i4 = improperlist[n][3]; + type = improperlist[n][4]; + + // geometry of 4-body + // 1 is the central atom + // 2-3-4 are ment to be equivalent + // I need the bonds between 2-3 and 2-4 to get the plane normal + // Then I need the bond 1-2 to project it onto the normal to the plane + + // bond 1->2 + xab = x[i2][0] - x[i1][0]; + yab = x[i2][1] - x[i1][1]; + zab = x[i2][2] - x[i1][2]; + domain->minimum_image(xab,yab,zab); + + // bond 1->3 + xac = x[i3][0] - x[i1][0]; + yac = x[i3][1] - x[i1][1]; + zac = x[i3][2] - x[i1][2]; + domain->minimum_image(xac,yac,zac); + + // bond 1->4 + xad = x[i4][0] - x[i1][0]; + yad = x[i4][1] - x[i1][1]; + zad = x[i4][2] - x[i1][2]; + domain->minimum_image(xad,yad,zad); + + // bond 2-3 + xbc = x[i3][0] - x[i2][0]; + ybc = x[i3][1] - x[i2][1]; + zbc = x[i3][2] - x[i2][2]; + domain->minimum_image(xbc,ybc,zbc); + + // bond 2-4 + xbd = x[i4][0] - x[i2][0]; + ybd = x[i4][1] - x[i2][1]; + zbd = x[i4][2] - x[i2][2]; + domain->minimum_image(xbd,ybd,zbd); + + xna = ybc*zbd - zbc*ybd; + yna = -(xbc*zbd - zbc*xbd); + zna = xbc*ybd - ybc*xbd; + rna = 1.0 / sqrt(xna*xna+yna*yna+zna*zna); + xna *= rna; + yna *= rna; + zna *= rna; + + da = xna*xab + yna*yab + zna*zab; + + domega = k[type]*da*da + chi[type]*da*da*da*da; + //printf("%3i %3i %3i %3i %10.5f %10.5f \n",i1,i2,i3,i4,da,domega); + a = 2.0* (k[type]*da + 2.0*chi[type]*da*da*da); + + if (eflag) eimproper = domega; + + f1[0] = a*( xna); + f1[1] = a*( yna); + f1[2] = a*( zna); + + f2[0] = a*( -xna -yab*(zbd-zbc)*rna +zab*(ybd-ybc)*rna -da*( -yna*(zbd-zbc) + zna*(ybd-ybc) )*rna); + f2[1] = a*( +xab*(zbd-zbc)*rna -yna +zab*(xbc-xbd)*rna -da*( +xna*(zbd-zbc) + zna*(xbc-xbd) )*rna); + f2[2] = a*( -xab*(ybd-ybc)*rna -yab*(xbc-xbd)*rna -zna -da*( +xna*(ybc-ybd) - yna*(xbc-xbd) )*rna); + + f3[0] = a*( ( yab*zbd -zab*ybd ) *rna +da*( -yna*zbd +zna*ybd )*rna); + f3[1] = a*( ( -xab*zbd +zab*xbd ) *rna +da*( +xna*zbd -zna*xbd )*rna); + f3[2] = a*( ( +xab*ybd -yab*xbd ) *rna +da*( -xna*ybd +yna*xbd )*rna); + + f4[0] = a*( ( -yab*zbc +zab*ybc ) *rna -da*( -yna*zbc +zna*ybc )*rna); + f4[1] = a*( ( +xab*zbc -zab*xbc ) *rna -da*( +xna*zbc -zna*xbc )*rna); + f4[2] = a*( ( -xab*ybc +yab*xbc ) *rna -da*( -xna*ybc +yna*xbc )*rna); + //printf("%10.5f %10.5f %10.5f \n",f1[0],f1[1],f1[2]); + //printf("%10.5f %10.5f %10.5f \n",f2[0],f2[1],f2[2]); + //printf("%10.5f %10.5f %10.5f \n",f3[0],f3[1],f3[2]); + //printf("%10.5f %10.5f %10.5f \n",f4[0],f4[1],f4[2]); + + // apply force to each of 4 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] += f2[0]; + f[i2][1] += f2[1]; + f[i2][2] += f2[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (newton_bond || i4 < nlocal) { + f[i4][0] += f4[0]; + f[i4][1] += f4[1]; + f[i4][2] += f4[2]; + } + + if (evflag) + ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f2,f3,f4, + xab,yab,zab,xac,yac,zac,xad-xac,yad-yac,zad-zac); + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperDistance::allocate() +{ + allocated = 1; + int n = atom->nimpropertypes; + + memory->create(k,n+1,"improper:k"); + memory->create(chi,n+1,"improper:chi"); + + memory->create(setflag,n+1,"improper:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void ImproperDistance::coeff(int narg, char **arg) +{ +// if (which > 0) return; + if (narg != 3) error->all(FLERR,"Incorrect args for improper coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(arg[0],atom->nimpropertypes,ilo,ihi); + + double k_one = force->numeric(FLERR,arg[1]); + double chi_one = force->numeric(FLERR,arg[2]); + + // convert chi from degrees to radians + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + //chi[i] = chi_one/180.0 * PI; + chi[i] = chi_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void ImproperDistance::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp); + fwrite(&chi[1],sizeof(double),atom->nimpropertypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void ImproperDistance::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k[1],sizeof(double),atom->nimpropertypes,fp); + fread(&chi[1],sizeof(double),atom->nimpropertypes,fp); + } + MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + MPI_Bcast(&chi[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1; +} diff --git a/src/USER-MISC/improper_distance.h b/src/USER-MISC/improper_distance.h new file mode 100644 index 0000000000..6301e6897a --- /dev/null +++ b/src/USER-MISC/improper_distance.h @@ -0,0 +1,47 @@ +/* ---------------------------------------------------------------------- + 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 IMPROPER_CLASS + +ImproperStyle(distance,ImproperDistance) + +#else + +#ifndef LMP_IMPROPER_DISTANCE_H +#define LMP_IMPROPER_DISTANCE_H + +#include +#include "improper.h" + +namespace LAMMPS_NS { + +class ImproperDistance : public Improper { + public: + ImproperDistance(class LAMMPS *); + ~ImproperDistance(); + void compute(int, int); + void coeff(int, char **); + void write_restart(FILE *); + void read_restart(FILE *); + + private: + double *k,*chi; + + void allocate(); +}; + +} + +#endif +#endif + diff --git a/src/USER-MISC/pair_buck_mdf.cpp b/src/USER-MISC/pair_buck_mdf.cpp new file mode 100644 index 0000000000..96f49e42dd --- /dev/null +++ b/src/USER-MISC/pair_buck_mdf.cpp @@ -0,0 +1,433 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Paolo Raiteri (Curtin University) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_buck_mdf.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; + +/* ---------------------------------------------------------------------- */ + +PairBuckMDF::PairBuckMDF(LAMMPS *lmp) : Pair(lmp) {} + +/* ---------------------------------------------------------------------- */ + +PairBuckMDF::~PairBuckMDF() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(cut_inner); + memory->destroy(cut_inner_sq); + memory->destroy(a); + memory->destroy(rho); + memory->destroy(c); + memory->destroy(rhoinv); + memory->destroy(buck1); + memory->destroy(buck2); + memory->destroy(offset); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairBuckMDF::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,r2inv,r6inv,forcebuck,factor_lj; + double r,rexp,phibuck; + int *ilist,*jlist,*numneigh,**firstneigh; + double dp, d, tt, dt, dd; + + 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]) { + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + r = sqrt(rsq); + rexp = exp(-r*rhoinv[itype][jtype]); + forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; + if (rsq > cut_inner_sq[itype][jtype]) { + phibuck = a[itype][jtype]*rexp - c[itype][jtype]*r6inv; + dp = (cut[itype][jtype] - cut_inner[itype][jtype]); + d = (r-cut_inner[itype][jtype]) / dp; + dd = 1.-d; + // taperig function - mdf style + tt = (1. + 3.*d + 6.*d*d)*dd*dd*dd; + // minus the derivative of the tapering function + dt = 30.* d*d * dd*dd * r / dp; + + forcebuck = forcebuck*tt + phibuck*dt; + } else { + tt = 1; + } + + fpair = factor_lj*forcebuck*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) { + evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv; + if (rsq > cut_inner_sq[itype][jtype]) evdwl *= tt; + + 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 PairBuckMDF::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_lj"); + memory->create(cut_inner,n+1,n+1,"pair:cut_inner"); + memory->create(cut_inner_sq,n+1,n+1,"pair:cut_inner_sq"); + memory->create(a,n+1,n+1,"pair:a"); + memory->create(rho,n+1,n+1,"pair:rho"); + memory->create(c,n+1,n+1,"pair:c"); + memory->create(rhoinv,n+1,n+1,"pair:rhoinv"); + memory->create(buck1,n+1,n+1,"pair:buck1"); + memory->create(buck2,n+1,n+1,"pair:buck2"); + memory->create(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairBuckMDF::settings(int narg, char **arg) +{ + if (narg != 2) error->all(FLERR,"Illegal pair_style command"); + + cut_inner_global = 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+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairBuckMDF::coeff(int narg, char **arg) +{ + if (narg != 5 && narg != 7) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double a_one = force->numeric(FLERR,arg[2]); + double rho_one = force->numeric(FLERR,arg[3]); + if (rho_one <= 0) error->all(FLERR,"Incorrect args for pair coefficients"); + double c_one = force->numeric(FLERR,arg[4]); + + double cut_inner_one = cut_inner_global; + double cut_one = cut_global; + if (narg == 7) { + cut_inner_one = force->numeric(FLERR,arg[5]); + cut_one = force->numeric(FLERR,arg[6]); + } + if (cut_inner_global <= 0.0 || cut_inner_global > cut_global) + error->all(FLERR,"Illegal pair_style command"); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + a[i][j] = a_one; + rho[i][j] = rho_one; + c[i][j] = c_one; + cut[i][j] = cut_one; + cut[j][i] = cut_one; + cut_inner[i][j] = cut_inner_one; + cut_inner[j][i] = cut_inner_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 PairBuckMDF::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + rhoinv[i][j] = 1.0/rho[i][j]; + buck1[i][j] = a[i][j]/rho[i][j]; + buck2[i][j] = 6.0*c[i][j]; + + if (offset_flag) { + double rexp = exp(-cut[i][j]/rho[i][j]); + offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut[i][j],6.0); + } else offset[i][j] = 0.0; + + cut_inner[j][i] = cut_inner[i][j]; + cut_inner_sq[i][j] = cut_inner[i][j]*cut_inner[i][j]; + cut_inner_sq[j][i] = cut_inner_sq[i][j]; + + a[j][i] = a[i][j]; + c[j][i] = c[i][j]; + rhoinv[j][i] = rhoinv[i][j]; + buck1[j][i] = buck1[i][j]; + buck2[j][i] = buck2[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 rho1 = rho[i][j]; + double rho2 = rho1*rho1; + double rho3 = rho2*rho1; + double rc = cut[i][j]; + double rc2 = rc*rc; + double rc3 = rc2*rc; + etail_ij = 2.0*MY_PI*all[0]*all[1]* + (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) - + c[i][j]/(3.0*rc3)); + ptail_ij = (-1/3.0)*2.0*MY_PI*all[0]*all[1]* + (-a[i][j]*exp(-rc/rho1)* + (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3); + } + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairBuckMDF::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(&a[i][j],sizeof(double),1,fp); + fwrite(&rho[i][j],sizeof(double),1,fp); + fwrite(&c[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairBuckMDF::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(&a[i][j],sizeof(double),1,fp); + fread(&rho[i][j],sizeof(double),1,fp); + fread(&c[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&a[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&rho[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&c[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairBuckMDF::write_restart_settings(FILE *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 PairBuckMDF::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + 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(&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); +} + +/* ---------------------------------------------------------------------- */ + +double PairBuckMDF::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,rexp,forcebuck,phibuck; + double dp, d, tt, dt, dd; + + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + r = sqrt(rsq); + rexp = exp(-r*rhoinv[itype][jtype]); + forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; + phibuck = a[itype][jtype]*rexp - c[itype][jtype]*r6inv; + + if (rsq > cut_inner_sq[itype][jtype]) { + + dp = (cut[itype][jtype] - cut_inner[itype][jtype]); + d = (r - cut_inner[itype][jtype]) / dp; + dd = 1-d; + tt = (1. + 3.*d + 6.*d*d)* dd*dd*dd; + dt = 30.* d*d * dd*dd * r / dp; + + forcebuck = forcebuck*tt + phibuck*dt; + phibuck *= tt; + } + + fforce = factor_lj*forcebuck*r2inv; + + return factor_lj*phibuck; +} + +/* ---------------------------------------------------------------------- */ + +void *PairBuckMDF::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"aparm") == 0) return (void *) a; + if (strcmp(str,"cparm") == 0) return (void *) c; + return NULL; +} diff --git a/src/USER-MISC/pair_buck_mdf.h b/src/USER-MISC/pair_buck_mdf.h new file mode 100644 index 0000000000..c0ecceeea8 --- /dev/null +++ b/src/USER-MISC/pair_buck_mdf.h @@ -0,0 +1,73 @@ +/* ---------------------------------------------------------------------- + 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(buck/mdf,PairBuckMDF) + +#else + +#ifndef LMP_PAIR_BUCK_MDF_H +#define LMP_PAIR_BUCK_MDF_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairBuckMDF : public Pair { + public: + PairBuckMDF(class LAMMPS *); + virtual ~PairBuckMDF(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double cut_global,cut_inner_global; + double **cut,**cut_inner,**cut_inner_sq; + double **a,**rho,**c; + double **rhoinv,**buck1,**buck2,**offset; + + 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: All pair coeffs are not set + +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. + +*/ diff --git a/src/USER-MISC/pair_lennard_mdf.cpp b/src/USER-MISC/pair_lennard_mdf.cpp new file mode 100644 index 0000000000..81c4435963 --- /dev/null +++ b/src/USER-MISC/pair_lennard_mdf.cpp @@ -0,0 +1,395 @@ + +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Paolo Raiteri (Curtin University) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_lennard_mdf.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; + +/* ---------------------------------------------------------------------- */ + +PairLJ_AB_MDF::PairLJ_AB_MDF(LAMMPS *lmp) : Pair(lmp) {} + +/* ---------------------------------------------------------------------- */ + +PairLJ_AB_MDF::~PairLJ_AB_MDF() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(cut_inner); + memory->destroy(cut_inner_sq); + memory->destroy(aparm); + memory->destroy(bparm); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJ_AB_MDF::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,r2inv,r6inv,forcelj,factor_lj; + 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; + double rr, d, dd, tt, dt, dp, philj; + + 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]) { + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + + if (rsq > cut_inner_sq[itype][jtype]) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + + rr = sqrt(rsq); + dp = (cut[itype][jtype] - cut_inner[itype][jtype]); + d = (rr-cut_inner[itype][jtype]) / dp; + dd = 1.-d; +// taperig function - mdf style + tt = (1. + 3.*d + 6.*d*d)*dd*dd*dd; +// minus the derivative of the tapering function + dt = 30.* d*d * dd*dd * rr / dp; + + forcelj = forcelj*tt + philj*dt; + } else { + tt = 1; + } + fpair = factor_lj*forcelj*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) { + evdwl = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); + if (rsq > cut_inner_sq[itype][jtype]) evdwl *= tt; + + 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 PairLJ_AB_MDF::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(cut_inner,n+1,n+1,"pair:cut_inner"); + memory->create(cut_inner_sq,n+1,n+1,"pair:cut_inner_sq"); + memory->create(aparm,n+1,n+1,"pair:aparm"); + memory->create(bparm,n+1,n+1,"pair:bparm"); + 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(lj4,n+1,n+1,"pair:lj4"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJ_AB_MDF::settings(int narg, char **arg) +{ + if (narg != 2) error->all(FLERR,"Illegal pair_style command"); + + cut_inner_global = force->numeric(FLERR,arg[0]); + cut_global = force->numeric(FLERR,arg[1]); + + if (cut_inner_global <= 0.0 || cut_inner_global > cut_global) + error->all(FLERR,"Illegal pair_style command"); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_inner[i][j] = cut_inner_global; + cut[i][j] = cut_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJ_AB_MDF::coeff(int narg, char **arg) +{ + if (narg != 4 && narg != 6) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double aparm_one = force->numeric(FLERR,arg[2]); + double bparm_one = force->numeric(FLERR,arg[3]); + + double cut_inner_one = cut_inner_global; + double cut_one = cut_global; + if (narg == 6) { + cut_inner_one = force->numeric(FLERR,arg[4]); + cut_one = force->numeric(FLERR,arg[5]); + } + if (cut_inner_global <= 0.0 || cut_inner_global > cut_global) + error->all(FLERR,"Illegal pair_style command"); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + aparm[i][j] = aparm_one; + bparm[i][j] = bparm_one; + cut_inner[i][j] = cut_inner_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 PairLJ_AB_MDF::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + aparm[i][j] = mix_energy(aparm[i][i],aparm[j][j], + bparm[i][i],bparm[j][j]); + bparm[i][j] = mix_distance(bparm[i][i],bparm[j][j]); + cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + cut_inner_sq[i][j] = cut_inner[i][j]*cut_inner[i][j]; + + lj1[i][j] = 12.0 * aparm[i][j]; + lj2[i][j] = 6.0 * bparm[i][j]; + lj3[i][j] = aparm[i][j]; + lj4[i][j] = bparm[i][j]; + + cut_inner[j][i] = cut_inner[i][j]; + cut_inner_sq[j][i] = cut_inner_sq[i][j]; + lj1[j][i] = lj1[i][j]; + lj2[j][i] = lj2[i][j]; + lj3[j][i] = lj3[i][j]; + lj4[j][i] = lj4[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJ_AB_MDF::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(&aparm[i][j],sizeof(double),1,fp); + fwrite(&bparm[i][j],sizeof(double),1,fp); + fwrite(&cut_inner[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJ_AB_MDF::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(&aparm[i][j],sizeof(double),1,fp); + fread(&bparm[i][j],sizeof(double),1,fp); + fread(&cut_inner[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&aparm[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&bparm[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJ_AB_MDF::write_restart_settings(FILE *fp) +{ + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJ_AB_MDF::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJ_AB_MDF::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,forcelj,philj; + double rr, dp, d, tt, dt, dd; + + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + + if (rsq > cut_inner_sq[itype][jtype]) { + + rr = sqrt(rsq); + dp = (cut[itype][jtype] - cut_inner[itype][jtype]); + d = (rr - cut_inner[itype][jtype]) / dp; + dd = 1-d; + tt = (1. + 3.*d + 6.*d*d)* dd*dd*dd; + dt = 30.* d*d * dd*dd * rr / dp; + + forcelj = forcelj*tt + philj*dt; + philj *= dt; + } + + fforce = factor_lj*forcelj*r2inv; + + return factor_lj*philj; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJ_AB_MDF::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"aparm") == 0) return (void *) aparm; + if (strcmp(str,"bparm") == 0) return (void *) bparm; + return NULL; +} diff --git a/src/USER-MISC/pair_lennard_mdf.h b/src/USER-MISC/pair_lennard_mdf.h new file mode 100644 index 0000000000..6fefe6fc3f --- /dev/null +++ b/src/USER-MISC/pair_lennard_mdf.h @@ -0,0 +1,69 @@ +/* ---------------------------------------------------------------------- + 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(lennard/mdf,PairLJ_AB_MDF) + +#else + +#ifndef LMP_PAIR_LJ_AB_MDF_H +#define LMP_PAIR_LJ_AB_MDF_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJ_AB_MDF : public Pair { + public: + PairLJ_AB_MDF(class LAMMPS *); + virtual ~PairLJ_AB_MDF(); + + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double cut_global,cut_inner_global; + double **cut,**cut_inner,**cut_inner_sq; + double **aparm,**bparm; + double **lj1,**lj2,**lj3,**lj4; + + 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. + +*/ diff --git a/src/USER-MISC/pair_lj_mdf.cpp b/src/USER-MISC/pair_lj_mdf.cpp new file mode 100644 index 0000000000..7c3b8d586e --- /dev/null +++ b/src/USER-MISC/pair_lj_mdf.cpp @@ -0,0 +1,395 @@ + +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Paolo Raiteri (Curtin University) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_lj_mdf.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; + +/* ---------------------------------------------------------------------- */ + +PairLJMDF::PairLJMDF(LAMMPS *lmp) : Pair(lmp) {} + +/* ---------------------------------------------------------------------- */ + +PairLJMDF::~PairLJMDF() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(cut_inner); + memory->destroy(cut_inner_sq); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJMDF::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,r2inv,r6inv,forcelj,factor_lj; + 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; + double rr, d, dd, tt, dt, dp, philj; + + 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]) { + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + + if (rsq > cut_inner_sq[itype][jtype]) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + + rr = sqrt(rsq); + dp = (cut[itype][jtype] - cut_inner[itype][jtype]); + d = (rr-cut_inner[itype][jtype]) / dp; + dd = 1.-d; +// taperig function - mdf style + tt = (1. + 3.*d + 6.*d*d)*dd*dd*dd; +// minus the derivative of the tapering function + dt = 30.* d*d * dd*dd * rr / dp; + + forcelj = forcelj*tt + philj*dt; + } else { + tt = 1; + } + fpair = factor_lj*forcelj*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) { + evdwl = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); + if (rsq > cut_inner_sq[itype][jtype]) evdwl *= tt; + + 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 PairLJMDF::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(cut_inner,n+1,n+1,"pair:cut_inner"); + memory->create(cut_inner_sq,n+1,n+1,"pair:cut_inner_sq"); + memory->create(epsilon,n+1,n+1,"pair:epsilon"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + 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(lj4,n+1,n+1,"pair:lj4"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJMDF::settings(int narg, char **arg) +{ + if (narg != 2) error->all(FLERR,"Illegal pair_style command"); + + cut_inner_global = force->numeric(FLERR,arg[0]); + cut_global = force->numeric(FLERR,arg[1]); + + if (cut_inner_global <= 0.0 || cut_inner_global > cut_global) + error->all(FLERR,"Illegal pair_style command"); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_inner[i][j] = cut_inner_global; + cut[i][j] = cut_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJMDF::coeff(int narg, char **arg) +{ + if (narg != 4 && narg != 6) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = force->numeric(FLERR,arg[2]); + double sigma_one = force->numeric(FLERR,arg[3]); + + double cut_inner_one = cut_inner_global; + double cut_one = cut_global; + if (narg == 6) { + cut_inner_one = force->numeric(FLERR,arg[4]); + cut_one = force->numeric(FLERR,arg[5]); + } + if (cut_inner_global <= 0.0 || cut_inner_global > cut_global) + error->all(FLERR,"Illegal pair_style command"); + + 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; + cut_inner[i][j] = cut_inner_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 PairLJMDF::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], + sigma[i][i],sigma[j][j]); + sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); + cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + cut_inner_sq[i][j] = cut_inner[i][j]*cut_inner[i][j]; + lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + + cut[j][i] = cut[i][j]; // BUG FIX + cut_inner[j][i] = cut_inner[i][j]; + cut_inner_sq[j][i] = cut_inner_sq[i][j]; + lj1[j][i] = lj1[i][j]; + lj2[j][i] = lj2[i][j]; + lj3[j][i] = lj3[i][j]; + lj4[j][i] = lj4[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJMDF::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(&cut_inner[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJMDF::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(&cut_inner[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(&cut_inner[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJMDF::write_restart_settings(FILE *fp) +{ + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJMDF::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJMDF::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,forcelj,philj; + double rr, dp, d, tt, dt, dd; + + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + + if (rsq > cut_inner_sq[itype][jtype]) { + + rr = sqrt(rsq); + dp = (cut[itype][jtype] - cut_inner[itype][jtype]); + d = (rr - cut_inner[itype][jtype]) / dp; + dd = 1-d; + tt = (1. + 3.*d + 6.*d*d)* dd*dd*dd; + dt = 30.* d*d * dd*dd * rr / dp; + + forcelj = forcelj*tt + philj*dt; + philj *= tt; + } + + fforce = factor_lj*forcelj*r2inv; + + return factor_lj*philj; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJMDF::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + return NULL; +} diff --git a/src/USER-MISC/pair_lj_mdf.h b/src/USER-MISC/pair_lj_mdf.h new file mode 100644 index 0000000000..f49020bf2d --- /dev/null +++ b/src/USER-MISC/pair_lj_mdf.h @@ -0,0 +1,69 @@ +/* ---------------------------------------------------------------------- + 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/mdf,PairLJMDF) + +#else + +#ifndef LMP_PAIR_LJ_MDF_H +#define LMP_PAIR_LJ_MDF_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJMDF : public Pair { + public: + PairLJMDF(class LAMMPS *); + virtual ~PairLJMDF(); + + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + double single(int, int, int, int, double, double, double, double &); + void *extract(const char *, int &); + + protected: + double cut_global,cut_inner_global; + double **cut,**cut_inner,**cut_inner_sq; + double **epsilon,**sigma; + double **lj1,**lj2,**lj3,**lj4; + + 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. + +*/