ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10029 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2013-06-06 14:51:09 +00:00
parent d7cd69f288
commit 94d3ef2744
6 changed files with 1540 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
src/DIPOLE/pair_dipole_cut.cpp → src/DIPOLE/pair_lj_cut_dipole_cut.cpp Normal file → Executable file
View File

@ -13,7 +13,7 @@
#include "math.h"
#include "stdlib.h"
#include "pair_dipole_cut.h"
#include "pair_lj_cut_dipole_cut.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_list.h"
@ -21,19 +21,21 @@
#include "force.h"
#include "memory.h"
#include "error.h"
#include "update.h"
#include "string.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairDipoleCut::PairDipoleCut(LAMMPS *lmp) : Pair(lmp)
PairLJCutDipoleCut::PairLJCutDipoleCut(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
}
/* ---------------------------------------------------------------------- */
PairDipoleCut::~PairDipoleCut()
PairLJCutDipoleCut::~PairLJCutDipoleCut()
{
if (allocated) {
memory->destroy(setflag);
@ -55,7 +57,7 @@ PairDipoleCut::~PairDipoleCut()
/* ---------------------------------------------------------------------- */
void PairDipoleCut::compute(int eflag, int vflag)
void PairLJCutDipoleCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fx,fy,fz;
@ -259,7 +261,7 @@ void PairDipoleCut::compute(int eflag, int vflag)
allocate all arrays
------------------------------------------------------------------------- */
void PairDipoleCut::allocate()
void PairLJCutDipoleCut::allocate()
{
allocated = 1;
int n = atom->ntypes;
@ -288,11 +290,14 @@ void PairDipoleCut::allocate()
global settings
------------------------------------------------------------------------- */
void PairDipoleCut::settings(int narg, char **arg)
void PairLJCutDipoleCut::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2)
error->all(FLERR,"Incorrect args in pair_style command");
if (strcmp(update->unit_style,"electron") == 0)
error->all(FLERR,"Cannot (yet) use 'electron' units with dipoles");
cut_lj_global = force->numeric(arg[0]);
if (narg == 1) cut_coul_global = cut_lj_global;
else cut_coul_global = force->numeric(arg[1]);
@ -314,7 +319,7 @@ void PairDipoleCut::settings(int narg, char **arg)
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairDipoleCut::coeff(int narg, char **arg)
void PairLJCutDipoleCut::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 6)
error->all(FLERR,"Incorrect args for pair coefficients");
@ -351,7 +356,7 @@ void PairDipoleCut::coeff(int narg, char **arg)
init specific to this pair style
------------------------------------------------------------------------- */
void PairDipoleCut::init_style()
void PairLJCutDipoleCut::init_style()
{
if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag)
error->all(FLERR,"Pair dipole/cut requires atom attributes q, mu, torque");
@ -363,7 +368,7 @@ void PairDipoleCut::init_style()
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairDipoleCut::init_one(int i, int j)
double PairLJCutDipoleCut::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
@ -402,7 +407,7 @@ double PairDipoleCut::init_one(int i, int j)
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairDipoleCut::write_restart(FILE *fp)
void PairLJCutDipoleCut::write_restart(FILE *fp)
{
write_restart_settings(fp);
@ -423,7 +428,7 @@ void PairDipoleCut::write_restart(FILE *fp)
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairDipoleCut::read_restart(FILE *fp)
void PairLJCutDipoleCut::read_restart(FILE *fp)
{
read_restart_settings(fp);
@ -454,7 +459,7 @@ void PairDipoleCut::read_restart(FILE *fp)
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairDipoleCut::write_restart_settings(FILE *fp)
void PairLJCutDipoleCut::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul_global,sizeof(double),1,fp);
@ -466,7 +471,7 @@ void PairDipoleCut::write_restart_settings(FILE *fp)
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairDipoleCut::read_restart_settings(FILE *fp)
void PairLJCutDipoleCut::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_lj_global,sizeof(double),1,fp);

144
src/DIPOLE/pair_dipole_cut.h → src/DIPOLE/pair_lj_cut_dipole_cut.h Normal file → Executable file
View File

@ -1,70 +1,74 @@
/* ----------------------------------------------------------------------
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(dipole/cut,PairDipoleCut)
#else
#ifndef LMP_PAIR_DIPOLE_CUT_H
#define LMP_PAIR_DIPOLE_CUT_H
#include "pair.h"
namespace LAMMPS_NS {
class PairDipoleCut : public Pair {
public:
PairDipoleCut(class LAMMPS *);
virtual ~PairDipoleCut();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
protected:
double cut_lj_global,cut_coul_global;
double **cut_lj,**cut_ljsq;
double **cut_coul,**cut_coulsq;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**offset;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Incorrect args in pair_style command
Self-explanatory.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair dipole/cut requires atom attributes q, mu, torque
The atom style defined does not have these attributes.
*/
/* ----------------------------------------------------------------------
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/cut/dipole/cut,PairLJCutDipoleCut)
#else
#ifndef LMP_PAIR_LJ_CUT_DIPOLE_CUT_H
#define LMP_PAIR_LJ_CUT_DIPOLE_CUT_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJCutDipoleCut : public Pair {
public:
PairLJCutDipoleCut(class LAMMPS *);
virtual ~PairLJCutDipoleCut();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
protected:
double cut_lj_global,cut_coul_global;
double **cut_lj,**cut_ljsq;
double **cut_coul,**cut_coulsq;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**offset;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Incorrect args in pair_style command
Self-explanatory.
E: Cannot (yet) use 'electron' units with dipoles
This feature is not yet supported.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair dipole/cut requires atom attributes q, mu, torque
The atom style defined does not have these attributes.
*/

560
src/DIPOLE/pair_lj_cut_dipole_long.cpp Executable file
View File

@ -0,0 +1,560 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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 "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_lj_cut_dipole_long.h"
#include "atom.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "force.h"
#include "kspace.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "update.h"
#include "string.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
/* ---------------------------------------------------------------------- */
PairLJCutDipoleLong::PairLJCutDipoleLong(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
ewaldflag = dipoleflag = 1;
respa_enable = 0;
}
/* ----------------------------------------------------------------------
free all arrays
------------------------------------------------------------------------- */
PairLJCutDipoleLong::~PairLJCutDipoleLong()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut_lj);
memory->destroy(cut_ljsq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(offset);
}
}
/* ---------------------------------------------------------------------- */
void PairLJCutDipoleLong::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz;
double rsq,r,rinv,r2inv,r6inv;
double forcecoulx,forcecouly,forcecoulz,fforce;
double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul;
double fx,fy,fz,fdx,fdy,fdz,fax,fay,faz;
double pdotp,pidotr,pjdotr,pre1,pre2,pre3;
double grij,expm2,t,erfc;
double g0,g1,g2,b0,b1,b2,b3,d0,d1,d2,d3;
double zdix,zdiy,zdiz,zdjx,zdjy,zdjz,zaix,zaiy,zaiz,zajx,zajy,zajz;
double g0b1_g1b2_g2b3,g0d1_g1d2_g2d3;
double forcelj,factor_coul,factor_lj,facm1;
double evdwl,ecoul;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 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;
double **mu = atom->mu;
double **torque = atom->torque;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
pre1 = 2.0 * g_ewald / MY_PIS;
pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = atom->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_lj = special_lj[sbmask(j)];
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;
rinv = sqrt(r2inv);
if (rsq < cut_coulsq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];
pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
g0 = qtmp*q[j];
g1 = qtmp*pjdotr - q[j]*pidotr + pdotp;
g2 = -pidotr*pjdotr;
if (factor_coul > 0.0) {
b0 = erfc * rinv;
b1 = (b0 + pre1*expm2) * r2inv;
b2 = (3.0*b1 + pre2*expm2) * r2inv;
b3 = (5.0*b2 + pre3*expm2) * r2inv;
g0b1_g1b2_g2b3 = g0*b1 + g1*b2 + g2*b3;
fdx = delx * g0b1_g1b2_g2b3 -
b1 * (qtmp*mu[j][0] - q[j]*mu[i][0]) +
b2 * (pjdotr*mu[i][0] + pidotr*mu[j][0]);
fdy = dely * g0b1_g1b2_g2b3 -
b1 * (qtmp*mu[j][1] - q[j]*mu[i][1]) +
b2 * (pjdotr*mu[i][1] + pidotr*mu[j][1]);
fdz = delz * g0b1_g1b2_g2b3 -
b1 * (qtmp*mu[j][2] - q[j]*mu[i][2]) +
b2 * (pjdotr*mu[i][2] + pidotr*mu[j][2]);
zdix = delx * (q[j]*b1 + b2*pjdotr) - b1*mu[j][0];
zdiy = dely * (q[j]*b1 + b2*pjdotr) - b1*mu[j][1];
zdiz = delz * (q[j]*b1 + b2*pjdotr) - b1*mu[j][2];
zdjx = delx * (-qtmp*b1 + b2*pidotr) - b1*mu[i][0];
zdjy = dely * (-qtmp*b1 + b2*pidotr) - b1*mu[i][1];
zdjz = delz * (-qtmp*b1 + b2*pidotr) - b1*mu[i][2];
if (factor_coul < 1.0) {
fdx *= factor_coul;
fdy *= factor_coul;
fdz *= factor_coul;
zdix *= factor_coul;
zdiy *= factor_coul;
zdiz *= factor_coul;
zdjx *= factor_coul;
zdjy *= factor_coul;
zdjz *= factor_coul;
}
} else {
fdx = fdy = fdz = 0.0;
zdix = zdiy = zdiz = 0.0;
zdjx = zdjy = zdjz = 0.0;
}
if (factor_coul < 1.0) {
d0 = (erfc - 1.0) * rinv;
d1 = (d0 + pre1*expm2) * r2inv;
d2 = (3.0*d1 + pre2*expm2) * r2inv;
d3 = (5.0*d2 + pre3*expm2) * r2inv;
g0d1_g1d2_g2d3 = g0*d1 + g1*d2 + g2*d3;
fax = delx * g0d1_g1d2_g2d3 -
d1 * (qtmp*mu[j][0] - q[j]*mu[i][0]) +
d2 * (pjdotr*mu[i][0] + pidotr*mu[j][0]);
fay = dely * g0d1_g1d2_g2d3 -
d1 * (qtmp*mu[j][1] - q[j]*mu[i][1]) +
d2 * (pjdotr*mu[i][1] + pidotr*mu[j][1]);
faz = delz * g0d1_g1d2_g2d3 -
d1 * (qtmp*mu[j][2] - q[j]*mu[i][2]) +
d2 * (pjdotr*mu[i][2] + pidotr*mu[j][2]);
zaix = delx * (q[j]*d1 + d2*pjdotr) - d1*mu[j][0];
zaiy = dely * (q[j]*d1 + d2*pjdotr) - d1*mu[j][1];
zaiz = delz * (q[j]*d1 + d2*pjdotr) - d1*mu[j][2];
zajx = delx * (-qtmp*d1 + d2*pidotr) - d1*mu[i][0];
zajy = dely * (-qtmp*d1 + d2*pidotr) - d1*mu[i][1];
zajz = delz * (-qtmp*d1 + d2*pidotr) - d1*mu[i][2];
if (factor_coul > 0.0) {
facm1 = 1.0 - factor_coul;
fax *= facm1;
fay *= facm1;
faz *= facm1;
zaix *= facm1;
zaiy *= facm1;
zaiz *= facm1;
zajx *= facm1;
zajy *= facm1;
zajz *= facm1;
}
} else {
fax = fay = faz = 0.0;
zaix = zaiy = zaiz = 0.0;
zajx = zajy = zajz = 0.0;
}
forcecoulx = fdx + fax;
forcecouly = fdy + fay;
forcecoulz = fdz + faz;
tixcoul = mu[i][1]*(zdiz + zaiz) - mu[i][2]*(zdiy + zaiy);
tiycoul = mu[i][2]*(zdix + zaix) - mu[i][0]*(zdiz + zaiz);
tizcoul = mu[i][0]*(zdiy + zaiy) - mu[i][1]*(zdix + zaix);
tjxcoul = mu[j][1]*(zdjz + zajz) - mu[j][2]*(zdjy + zajy);
tjycoul = mu[j][2]*(zdjx + zajx) - mu[j][0]*(zdjz + zajz);
tjzcoul = mu[j][0]*(zdjy + zajy) - mu[j][1]*(zdjx + zajx);
} else {
forcecoulx = forcecouly = forcecoulz = 0.0;
tixcoul = tiycoul = tizcoul = 0.0;
tjxcoul = tjycoul = tjzcoul = 0.0;
}
// LJ interaction
if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
fforce = factor_lj * forcelj*r2inv;
} else fforce = 0.0;
// total force
fx = qqrd2e*forcecoulx + delx*fforce;
fy = qqrd2e*forcecouly + dely*fforce;
fz = qqrd2e*forcecoulz + delz*fforce;
// force & torque accumulation
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
torque[i][0] += qqrd2e*tixcoul;
torque[i][1] += qqrd2e*tiycoul;
torque[i][2] += qqrd2e*tizcoul;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
torque[j][0] += qqrd2e*tjxcoul;
torque[j][1] += qqrd2e*tjycoul;
torque[j][2] += qqrd2e*tjzcoul;
}
if (eflag) {
if (rsq < cut_coulsq) {
ecoul = qqrd2e*(b0*g0 + b1*g1 + b2*g2);
if (factor_coul < 1.0) {
ecoul *= factor_coul;
ecoul += (1-factor_coul) * qqrd2e * (d0*g0 + d1*g1 + d2*g2);
}
} else ecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl *= factor_lj;
} else evdwl = 0.0;
}
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fx,fy,fz,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJCutDipoleLong::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_lj,n+1,n+1,"pair:cut_lj");
memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
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");
memory->create(offset,n+1,n+1,"pair:offset");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJCutDipoleLong::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2)
error->all(FLERR,"Incorrect args in pair_style command");
cut_lj_global = atof(arg[0]);
if (narg == 1) cut_coul = cut_lj_global;
else cut_coul = atof(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_lj[i][j] = cut_lj_global;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJCutDipoleLong::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5)
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 = atof(arg[2]);
double sigma_one = atof(arg[3]);
double cut_lj_one = cut_lj_global;
if (narg == 5) cut_lj_one = atof(arg[4]);
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_lj[i][j] = cut_lj_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 PairLJCutDipoleLong::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_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
}
double cut = MAX(cut_lj[i][j],cut_coul);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[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);
if (offset_flag) {
double ratio = sigma[i][j] / cut_lj[i][j];
offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
} else offset[i][j] = 0.0;
cut_ljsq[j][i] = cut_ljsq[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];
offset[j][i] = offset[i][j];
return cut;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJCutDipoleLong::init_style()
{
if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag)
error->all(FLERR,"Pair dipole/long requires atom attributes q, mu, torque");
if (strcmp(update->unit_style,"electron") == 0)
error->all(FLERR,"Cannot (yet) use 'electron' units with dipoles");
// insure use of KSpace long-range solver, set g_ewald
if (force->kspace == NULL)
error->all(FLERR,"Pair style requires a KSpace style");
g_ewald = force->kspace->g_ewald;
cut_coulsq = cut_coul * cut_coul;
neighbor->request(this);
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJCutDipoleLong::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_lj[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJCutDipoleLong::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_lj[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_lj[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJCutDipoleLong::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul,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 PairLJCutDipoleLong::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_lj_global,sizeof(double),1,fp);
fread(&cut_coul,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
void *PairLJCutDipoleLong::extract(const char *str, int &dim)
{
if (strcmp(str,"cut_coul") == 0) {
dim = 0;
return (void *) &cut_coul;
} else if (strcmp(str,"ewald_order") == 0) {
ewald_order = 0;
ewald_order |= 1<<1;
ewald_order |= 1<<3;
dim = 0;
return (void *) &ewald_order;
} else if (strcmp(str,"ewald_mix") == 0) {
dim = 0;
return (void *) &mix_flag;
}
return NULL;
}

84
src/DIPOLE/pair_lj_cut_dipole_long.h Executable file
View File

@ -0,0 +1,84 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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/cut/dipole/long,PairLJCutDipoleLong)
#else
#ifndef LMP_PAIR_LJ_CUT_DIPOLE_LONG_H
#define LMP_PAIR_LJ_CUT_DIPOLE_LONG_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJCutDipoleLong : public Pair {
public:
double cut_coul;
double **sigma;
PairLJCutDipoleLong(class LAMMPS *);
~PairLJCutDipoleLong();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
private:
double cut_lj_global;
double **cut_lj,**cut_ljsq;
double cut_coulsq;
double **epsilon;
double **lj1,**lj2,**lj3,**lj4,**offset;
double g_ewald;
int ewald_order;
virtual void *extract(const char *, int &);
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Incorrect args in pair_style command
Self-explanatory.
E: Cannot (yet) use 'electron' units with dipoles
This feature is not yet supported.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair dipole/cut requires atom attributes q, mu, torque
The atom style defined does not have these attributes.
E: Pair style requires a KSpace style
No kspace style is defined.
*/

687
src/DIPOLE/pair_lj_long_dipole_long.cpp Executable file
View File

@ -0,0 +1,687 @@
/* ----------------------------------------------------------------------
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: Pieter J. in 't Veld and Stan Moore (Sandia)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "math_const.h"
#include "math_vector.h"
#include "pair_lj_long_dipole_long.h"
#include "atom.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "kspace.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
// ----------------------------------------------------------------------
PairLJLongDipoleLong::PairLJLongDipoleLong(LAMMPS *lmp) : Pair(lmp)
{
dispersionflag = ewaldflag = dipoleflag = 1;
respa_enable = 0;
single_enable = 0;
}
// ----------------------------------------------------------------------
// global settings
// ----------------------------------------------------------------------
#define PAIR_ILLEGAL "Illegal pair_style lj/long/dipole/long command"
#define PAIR_CUTOFF "Only one cut-off allowed when requesting all long"
#define PAIR_MISSING "Cut-offs missing in pair_style lj/long/dipole/long"
#define PAIR_COUL_CUT "Coulombic cut not supported in pair_style lj/long/dipole/long"
#define PAIR_LARGEST "Using largest cut-off for lj/long/dipole/long long long"
#define PAIR_MIX "Geometric mixing assumed for 1/r^6 coefficients"
void PairLJLongDipoleLong::options(char **arg, int order)
{
const char *option[] = {"long", "cut", "off", NULL};
int i;
if (!*arg) error->all(FLERR,PAIR_ILLEGAL);
for (i=0; option[i]&&strcmp(arg[0], option[i]); ++i);
switch (i) {
default: error->all(FLERR,PAIR_ILLEGAL);
case 0: ewald_order |= 1<<order; break; // set kspace r^-order
case 2: ewald_off |= 1<<order; // turn r^-order off
case 1: break;
}
}
void PairLJLongDipoleLong::settings(int narg, char **arg)
{
if (narg != 3 && narg != 4) error->all(FLERR,"Illegal pair_style command");
ewald_off = 0;
ewald_order = 0;
options(arg, 6);
options(++arg, 3);
options(arg, 1);
if (!comm->me && ewald_order&(1<<6))
error->warning(FLERR,PAIR_MIX);
if (!comm->me && ewald_order==((1<<3)|(1<<6)))
error->warning(FLERR,PAIR_LARGEST);
if (!*(++arg))
error->all(FLERR,PAIR_MISSING);
if (!((ewald_order^ewald_off)&(1<<3)))
error->all(FLERR,PAIR_COUL_CUT);
cut_lj_global = force->numeric(*(arg++));
if (narg == 4 && (ewald_order==74))
error->all(FLERR,PAIR_CUTOFF);
if (narg == 4) cut_coul = force->numeric(*(arg++));
else cut_coul = cut_lj_global;
if (allocated) { // reset explicit cuts
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) cut_lj[i][j] = cut_lj_global;
}
}
// ----------------------------------------------------------------------
// free all arrays
// ----------------------------------------------------------------------
PairLJLongDipoleLong::~PairLJLongDipoleLong()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut_lj_read);
memory->destroy(cut_lj);
memory->destroy(cut_ljsq);
memory->destroy(epsilon_read);
memory->destroy(epsilon);
memory->destroy(sigma_read);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(offset);
}
//if (ftable) free_tables();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::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_lj_read,n+1,n+1,"pair:cut_lj_read");
memory->create(cut_lj,n+1,n+1,"pair:cut_lj");
memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
memory->create(epsilon_read,n+1,n+1,"pair:epsilon_read");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma_read,n+1,n+1,"pair:sigma_read");
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");
memory->create(offset,n+1,n+1,"pair:offset");
}
/* ----------------------------------------------------------------------
extract protected data from object
------------------------------------------------------------------------- */
void *PairLJLongDipoleLong::extract(const char *id, int &dim)
{
const char *ids[] = {
"B", "sigma", "epsilon", "ewald_order", "ewald_cut", "ewald_mix",
"cut_coul", "cut_vdwl", NULL};
void *ptrs[] = {
lj4, sigma, epsilon, &ewald_order, &cut_coul, &mix_flag, &cut_coul,
&cut_lj_global, NULL};
int i;
for (i=0; ids[i]&&strcmp(ids[i], id); ++i);
if (i <= 2) dim = 2;
else dim = 0;
return ptrs[i];
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5) 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(arg[2]);
double sigma_one = force->numeric(arg[3]);
double cut_lj_one = cut_lj_global;
if (narg == 5) cut_lj_one = force->numeric(arg[4]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
epsilon_read[i][j] = epsilon_one;
sigma_read[i][j] = sigma_one;
cut_lj_read[i][j] = cut_lj_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::init_style()
{
const char *style3[] = {"ewald/disp", NULL};
const char *style6[] = {"ewald/disp", NULL};
int i;
if (strcmp(update->unit_style,"electron") == 0)
error->all(FLERR,"Cannot (yet) use 'electron' units with dipoles");
// require an atom style with charge defined
if (!atom->q_flag && (ewald_order&(1<<1)))
error->all(FLERR,
"Invoking coulombic in pair style lj/long/dipole/long requires atom attribute q");
if (!atom->mu && (ewald_order&(1<<3)))
error->all(FLERR,"Pair lj/long/dipole/long requires atom attributes mu, torque");
if (!atom->torque && (ewald_order&(1<<3)))
error->all(FLERR,"Pair lj/long/dipole/long requires atom attributes mu, torque");
neighbor->request(this);
cut_coulsq = cut_coul * cut_coul;
// ensure use of KSpace long-range solver, set g_ewald
if (ewald_order&(1<<3)) { // r^-1 kspace
if (force->kspace == NULL)
error->all(FLERR,"Pair style is incompatible with KSpace style");
for (i=0; style3[i]&&strcmp(force->kspace_style, style3[i]); ++i);
if (!style3[i])
error->all(FLERR,"Pair style is incompatible with KSpace style");
}
if (ewald_order&(1<<6)) { // r^-6 kspace
if (force->kspace == NULL)
error->all(FLERR,"Pair style is incompatible with KSpace style");
for (i=0; style6[i]&&strcmp(force->kspace_style, style6[i]); ++i);
if (!style6[i])
error->all(FLERR,"Pair style is incompatible with KSpace style");
}
if (force->kspace) g_ewald = force->kspace->g_ewald;
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
regular or rRESPA
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::init_list(int id, NeighList *ptr)
{
if (id == 0) list = ptr;
else if (id == 1) listinner = ptr;
else if (id == 2) listmiddle = ptr;
else if (id == 3) listouter = ptr;
if (id)
error->all(FLERR,"Pair style lj/long/dipole/long does not currently support respa");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJLongDipoleLong::init_one(int i, int j)
{
if ((ewald_order&(1<<6))||(setflag[i][j] == 0)) {
epsilon[i][j] = mix_energy(epsilon_read[i][i],epsilon_read[j][j],
sigma_read[i][i],sigma_read[j][j]);
sigma[i][j] = mix_distance(sigma_read[i][i],sigma_read[j][j]);
if (ewald_order&(1<<6))
cut_lj[i][j] = cut_lj_global;
else
cut_lj[i][j] = mix_distance(cut_lj_read[i][i],cut_lj_read[j][j]);
}
else {
sigma[i][j] = sigma_read[i][j];
epsilon[i][j] = epsilon_read[i][j];
cut_lj[i][j] = cut_lj_read[i][j];
}
double cut = MAX(cut_lj[i][j], cut_coul);
cutsq[i][j] = cut*cut;
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[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);
// check interior rRESPA cutoff
//if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3])
//error->all(FLERR,"Pair cutoff < Respa interior cutoff");
if (offset_flag) {
double ratio = sigma[i][j] / cut_lj[i][j];
offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
} else offset[i][j] = 0.0;
cutsq[j][i] = cutsq[i][j];
cut_ljsq[j][i] = cut_ljsq[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];
offset[j][i] = offset[i][j];
return cut;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::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_read[i][j],sizeof(double),1,fp);
fwrite(&sigma_read[i][j],sizeof(double),1,fp);
fwrite(&cut_lj_read[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::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_read[i][j],sizeof(double),1,fp);
fread(&sigma_read[i][j],sizeof(double),1,fp);
fread(&cut_lj_read[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&epsilon_read[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma_read[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_lj_read[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&ewald_order,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_lj_global,sizeof(double),1,fp);
fread(&cut_coul,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
fread(&ewald_order,sizeof(int),1,fp);
}
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul,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(&ewald_order,1,MPI_INT,0,world);
}
/* ----------------------------------------------------------------------
compute pair interactions
------------------------------------------------------------------------- */
void PairLJLongDipoleLong::compute(int eflag, int vflag)
{
double evdwl,ecoul,fpair;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x, *x0 = x[0];
double **mu = atom->mu, *mu0 = mu[0], *imu, *jmu;
double **tq = atom->torque, *tq0 = tq[0], *tqi;
double **f = atom->f, *f0 = f[0], *fi = f0, fx, fy, fz;
double *q = atom->q, qi = 0, qj;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
int i, j;
int order1 = ewald_order&(1<<1), order3 = ewald_order&(1<<3),
order6 = ewald_order&(1<<6);
int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni;
double *cutsqi, *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti;
double rsq, r2inv, force_coul, force_lj;
double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2;
double B0, B1, B2, B3, G0, G1, G2, mudi, mudj, muij;
vector force_d = VECTOR_NULL, ti = VECTOR_NULL, tj = VECTOR_NULL;
vector mui, muj, xi, d;
double C1 = 2.0 * g_ewald / MY_PIS;
double C2 = 2.0 * g2 * C1;
double C3 = 2.0 * g2 * C2;
ineighn = (ineigh = list->ilist)+list->inum;
for (; ineigh<ineighn; ++ineigh) { // loop over all neighs
i = *ineigh; fi = f0+3*i; tqi = tq0+3*i;
qi = q[i]; // initialize constants
offseti = offset[typei = type[i]];
lj1i = lj1[typei]; lj2i = lj2[typei]; lj3i = lj3[typei]; lj4i = lj4[typei];
cutsqi = cutsq[typei]; cut_ljsqi = cut_ljsq[typei];
memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
memcpy(mui, imu = mu0+(i<<2), sizeof(vector));
jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
for (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
ni = sbmask(j); // special index
j &= NEIGHMASK;
{ register double *xj = x0+(j+(j<<1));
d[0] = xi[0] - xj[0]; // pair vector
d[1] = xi[1] - xj[1];
d[2] = xi[2] - xj[2]; }
if ((rsq = vec_dot(d, d)) >= cutsqi[typej = type[j]]) continue;
r2inv = 1.0/rsq;
if (order3 && (rsq < cut_coulsq)) { // dipole
memcpy(muj, jmu = mu0+(j<<2), sizeof(vector));
{ // series real space
register double r = sqrt(rsq);
register double x = g_ewald*r;
register double f = exp(-x*x)*qqrd2e;
B0 = 1.0/(1.0+EWALD_P*x); // eqn 2.8
B0 *= ((((A5*B0+A4)*B0+A3)*B0+A2)*B0+A1)*f/r;
B1 = (B0 + C1 * f) * r2inv;
B2 = (3.0*B1 + C2 * f) * r2inv;
B3 = (5.0*B2 + C3 * f) * r2inv;
mudi = mui[0]*d[0]+mui[1]*d[1]+mui[2]*d[2];
mudj = muj[0]*d[0]+muj[1]*d[1]+muj[2]*d[2];
muij = mui[0]*muj[0]+mui[1]*muj[1]+mui[2]*muj[2];
G0 = qi*(qj = q[j]); // eqn 2.10
G1 = qi*mudj-qj*mudi+muij;
G2 = -mudi*mudj;
force_coul = G0*B1+G1*B2+G2*B3;
mudi *= B2; mudj *= B2; // torque contribs
ti[0] = mudj*d[0]+(qj*d[0]-muj[0])*B1;
ti[1] = mudj*d[1]+(qj*d[1]-muj[1])*B1;
ti[2] = mudj*d[2]+(qj*d[2]-muj[2])*B1;
if (newton_pair || j < nlocal) {
tj[0] = mudi*d[0]-(qi*d[0]+mui[0])*B1;
tj[1] = mudi*d[1]-(qi*d[1]+mui[1])*B1;
tj[2] = mudi*d[2]-(qi*d[2]+mui[2])*B1;
}
if (eflag) ecoul = G0*B0+G1*B1+G2*B2;
if (ni > 0) { // adj part, eqn 2.13
force_coul -= (f = qqrd2e*(1.0-special_coul[ni])/r)*(
(3.0*G1+15.0*G2*r2inv)*r2inv+G0)*r2inv;
if (eflag)
ecoul -= f*((G1+3.0*G2*r2inv)*r2inv+G0);
B1 -= f*r2inv;
}
B0 = mudj+qj*B1; B3 = -qi*B1+mudi; // position independent
if (ni > 0) B0 -= f*3.0*mudj*r2inv*r2inv/B2;
if (ni > 0) B3 -= f*3.0*mudi*r2inv*r2inv/B2;
force_d[0] = B0*mui[0]+B3*muj[0]; // force contribs
force_d[1] = B0*mui[1]+B3*muj[1];
force_d[2] = B0*mui[2]+B3*muj[2];
if (ni > 0) {
ti[0] -= f*(3.0*mudj*r2inv*r2inv*d[0]/B2+(qj*r2inv*d[0]-muj[0]*r2inv));
ti[1] -= f*(3.0*mudj*r2inv*r2inv*d[1]/B2+(qj*r2inv*d[1]-muj[1]*r2inv));
ti[2] -= f*(3.0*mudj*r2inv*r2inv*d[2]/B2+(qj*r2inv*d[2]-muj[2]*r2inv));
if (newton_pair || j < nlocal) {
tj[0] -= f*(3.0*mudi*r2inv*r2inv*d[0]/B2-(qi*r2inv*d[0]+mui[0]*r2inv));
tj[1] -= f*(3.0*mudi*r2inv*r2inv*d[1]/B2-(qi*r2inv*d[1]+mui[1]*r2inv));
tj[2] -= f*(3.0*mudi*r2inv*r2inv*d[2]/B2-(qi*r2inv*d[2]+mui[2]*r2inv));
}
}
} // table real space
} else {
force_coul = ecoul = 0.0;
memset(force_d, 0, 3*sizeof(double));
}
if (rsq < cut_ljsqi[typej]) { // lj
if (order6) { // long-range lj
register double rn = r2inv*r2inv*r2inv;
register double x2 = g2*rsq, a2 = 1.0/x2;
x2 = a2*exp(-x2)*lj4i[typej];
if (ni < 0) {
force_lj =
(rn*=rn)*lj1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
}
else { // special case
register double f = special_lj[ni], t = rn*(1.0-f);
force_lj = f*(rn *= rn)*lj1i[typej]-
g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
if (eflag) evdwl =
f*rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[typej];
}
}
else { // cut lj
register double rn = r2inv*r2inv*r2inv;
if (ni < 0) {
force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
}
else { // special case
register double f = special_lj[ni];
force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
if (eflag) evdwl = f*(
rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
}
}
force_lj *= r2inv;
}
else force_lj = evdwl = 0.0;
fpair = force_coul+force_lj; // force
if (newton_pair || j < nlocal) {
register double *fj = f0+(j+(j<<1));
fi[0] += fx = d[0]*fpair+force_d[0]; fj[0] -= fx;
fi[1] += fy = d[1]*fpair+force_d[1]; fj[1] -= fy;
fi[2] += fz = d[2]*fpair+force_d[2]; fj[2] -= fz;
tqi[0] += mui[1]*ti[2]-mui[2]*ti[1]; // torque
tqi[1] += mui[2]*ti[0]-mui[0]*ti[2];
tqi[2] += mui[0]*ti[1]-mui[1]*ti[0];
register double *tqj = tq0+(j+(j<<1));
tqj[0] += muj[1]*tj[2]-muj[2]*tj[1];
tqj[1] += muj[2]*tj[0]-muj[0]*tj[2];
tqj[2] += muj[0]*tj[1]-muj[1]*tj[0];
}
else {
fi[0] += fx = d[0]*fpair+force_d[0]; // force
fi[1] += fy = d[1]*fpair+force_d[1];
fi[2] += fz = d[2]*fpair+force_d[2];
tqi[0] += mui[1]*ti[2]-mui[2]*ti[1]; // torque
tqi[1] += mui[2]*ti[0]-mui[0]*ti[2];
tqi[2] += mui[0]*ti[1]-mui[1]*ti[0];
}
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fx,fy,fz,d[0],d[1],d[2]);
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
/*
double PairLJLongDipoleLong::single(int i, int j, int itype, int jtype,
double rsq, double factor_coul, double factor_lj,
double &fforce)
{
double r6inv, force_coul, force_lj;
double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2, *q = atom->q;
double eng = 0.0;
double r2inv = 1.0/rsq;
if ((ewald_order&(1<<3)) && (rsq < cut_coulsq)) { // coulombic
double *mui = atom->mu[i], *muj = atom->mu[j];
double *xi = atom->x[i], *xj = atom->x[j];
double qi = q[i], qj = q[j];
double G0, G1, G2, B0, B1, B2, B3, mudi, mudj, muij;
vector d = {xi[0]-xj[0], xi[1]-xj[1], xi[2]-xj[2]};
{ // series real space
register double r = sqrt(rsq);
register double x = g_ewald*r;
register double f = exp(-x*x)*qqrd2e;
B0 = 1.0/(1.0+EWALD_P*x); // eqn 2.8
B0 *= ((((A5*B0+A4)*B0+A3)*B0+A2)*B0+A1)*f/r;
B1 = (B0 + C1 * f) * r2inv;
B2 = (3.0*B1 + C2 * f) * r2inv;
B3 = (5.0*B2 + C3 * f) * r2inv;
mudi = mui[0]*d[0]+mui[1]*d[1]+mui[2]*d[2];
mudj = muj[0]*d[0]+muj[1]*d[1]+muj[2]*d[2];
muij = mui[0]*muj[0]+mui[1]*muj[1]+mui[2]*muj[2];
G0 = qi*(qj = q[j]); // eqn 2.10
G1 = qi*mudj-qj*mudi+muij;
G2 = -mudi*mudj;
force_coul = G0*B1+G1*B2+G2*B3;
eng += G0*B0+G1*B1+G2*B2;
if (factor_coul < 1.0) { // adj part, eqn 2.13
force_coul -= (f = force->qqrd2e*(1.0-factor_coul)/r)*(
(3.0*G1+6.0*muij+15.0*G2*r2inv)*r2inv+G0);
eng -= f*((G1+3.0*G2*r2inv)*r2inv+G0);
B1 -= f*r2inv;
}
B0 = mudj*B2-qj*B1; B3 = qi*B1+mudi*B2; // position independent
//force_d[0] = B0*mui[0]+B3*muj[0]; // force contributions
//force_d[1] = B0*mui[1]+B3*muj[1];
//force_d[2] = B0*mui[2]+B3*muj[2];
} // table real space
}
else force_coul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) { // lennard-jones
r6inv = r2inv*r2inv*r2inv;
if (ewald_order&0x40) { // long-range
register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
x2 = a2*exp(-x2)*lj4[itype][jtype];
force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]-
g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype];
eng += factor_lj*r6inv*lj3[itype][jtype]-
g6*((a2+1.0)*a2+0.5)*x2+t*lj4[itype][jtype];
}
else { // cut
force_lj = factor_lj*r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
eng += factor_lj*(r6inv*(r6inv*lj3[itype][jtype]-
lj4[itype][jtype])-offset[itype][jtype]);
}
}
else force_lj = 0.0;
fforce = (force_coul+force_lj)*r2inv;
return eng;
}
*/

117
src/DIPOLE/pair_lj_long_dipole_long.h Executable file
View File

@ -0,0 +1,117 @@
/* -*- 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/long/dipole/long,PairLJLongDipoleLong)
#else
#ifndef LMP_PAIR_LJ_LONG_DIPOLE_LONG_H
#define LMP_PAIR_LJ_LONG_DIPOLE_LONG_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJLongDipoleLong : public Pair {
public:
double cut_coul;
PairLJLongDipoleLong(class LAMMPS *);
virtual ~PairLJLongDipoleLong();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void *extract(const char *, int &);
protected:
double cut_lj_global;
double **cut_lj, **cut_lj_read, **cut_ljsq;
double cut_coulsq;
double **epsilon_read, **epsilon, **sigma_read, **sigma;
double **lj1, **lj2, **lj3, **lj4, **offset;
double *cut_respa;
double g_ewald;
int ewald_order, ewald_off;
void options(char **arg, int order);
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.
W: Geometric mixing assumed for 1/r^6 coefficients
Self-explanatory.
W: Using largest cutoff for lj/long/dipole/long
Self-exlanatory.
E: Cutoffs missing in pair_style lj/long/dipole/long
Self-exlanatory.
E: Coulombic cut not supported in pair_style lj/long/dipole/long
Must use long-range Coulombic interactions.
E: Only one cutoff allowed when requesting all long
Self-explanatory.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Cannot (yet) use 'electron' units with dipoles
This feature is not yet supported.
E: Invoking coulombic in pair style lj/long/dipole/long requires atom attribute q
The atom style defined does not have these attributes.
E: Pair lj/long/dipole/long requires atom attributes mu, torque
The atom style defined does not have these attributes.
E: Pair style is incompatible with KSpace style
Self-explanatory.
E: Pair style lj/long/dipole/long does not currently support respa
This feature is not yet supported.
*/