Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F68526675
pair_lj_class2_coul_long.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jun 27, 19:50
Size
14 KB
Mime Type
text/x-c
Expires
Sat, Jun 29, 19:50 (2 d)
Engine
blob
Format
Raw Data
Handle
18596326
Attached To
rLAMMPS lammps
pair_lj_class2_coul_long.cpp
View Options
/* ----------------------------------------------------------------------
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_class2_coul_long.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "kspace.h"
#include "update.h"
#include "memory.h"
#include "neighbor.h"
#include "error.h"
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#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
/* ----------------------------------------------------------------------
free all arrays
------------------------------------------------------------------------- */
PairLJClass2CoulLong::~PairLJClass2CoulLong()
{
if (allocated) {
memory->destroy_2d_int_array(setflag);
memory->destroy_2d_double_array(cutsq);
memory->destroy_2d_double_array(cut_lj);
memory->destroy_2d_double_array(cut_ljsq);
memory->destroy_2d_double_array(epsilon);
memory->destroy_2d_double_array(sigma);
memory->destroy_2d_double_array(lj1);
memory->destroy_2d_double_array(lj2);
memory->destroy_2d_double_array(lj3);
memory->destroy_2d_double_array(lj4);
memory->destroy_2d_double_array(offset);
}
}
/* ---------------------------------------------------------------------- */
void PairLJClass2CoulLong::compute(int eflag, int vflag)
{
int i,j,k,numneigh,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz;
double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,fforce;
double grij,expm2,prefactor,t,erfc;
double factor_coul,factor_lj,factor,phicoul,philj;
int *neighs;
double **f;
eng_vdwl = eng_coul = 0.0;
if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0;
if (vflag == 2) f = update->f_pair;
else f = atom->f;
double **x = atom->x;
double *q = atom->q;
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;
// loop over neighbors of my atoms
for (i = 0; i < nlocal; i++) {
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
neighs = neighbor->firstneigh[i];
numneigh = neighbor->numneigh[i];
for (k = 0; k < numneigh; k++) {
j = neighs[k];
if (j < nall) factor_coul = factor_lj = 1.0;
else {
factor_coul = special_coul[j/nall];
factor_lj = special_lj[j/nall];
j %= nall;
}
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;
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;
prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
} else forcelj = 0.0;
fforce = (forcecoul + factor_lj*forcelj) * r2inv;
f[i][0] += delx*fforce;
f[i][1] += dely*fforce;
f[i][2] += delz*fforce;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fforce;
f[j][1] -= dely*fforce;
f[j][2] -= delz*fforce;
}
if (eflag) {
if (newton_pair || j < nlocal) factor = 1.0;
else factor = 0.5;
if (rsq < cut_coulsq) {
phicoul = prefactor*erfc;
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
eng_coul += factor*phicoul;
}
if (rsq < cut_ljsq[itype][jtype]) {
philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
offset[itype][jtype];
eng_vdwl += factor*factor_lj*philj;
}
}
if (vflag == 1) {
if (newton_pair || j < nlocal) {
virial[0] += delx*delx*fforce;
virial[1] += dely*dely*fforce;
virial[2] += delz*delz*fforce;
virial[3] += delx*dely*fforce;
virial[4] += delx*delz*fforce;
virial[5] += dely*delz*fforce;
} else {
virial[0] += 0.5*delx*delx*fforce;
virial[1] += 0.5*dely*dely*fforce;
virial[2] += 0.5*delz*delz*fforce;
virial[3] += 0.5*delx*dely*fforce;
virial[4] += 0.5*delx*delz*fforce;
virial[5] += 0.5*dely*delz*fforce;
}
}
}
}
}
if (vflag == 2) virial_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJClass2CoulLong::allocate()
{
allocated = 1;
int n = atom->ntypes;
setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
cut_lj = memory->create_2d_double_array(n+1,n+1,"pair:cut_lj");
cut_ljsq = memory->create_2d_double_array(n+1,n+1,"pair:cut_ljsq");
epsilon = memory->create_2d_double_array(n+1,n+1,"pair:epsilon");
sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma");
lj1 = memory->create_2d_double_array(n+1,n+1,"pair:lj1");
lj2 = memory->create_2d_double_array(n+1,n+1,"pair:lj2");
lj3 = memory->create_2d_double_array(n+1,n+1,"pair:lj3");
lj4 = memory->create_2d_double_array(n+1,n+1,"pair:lj4");
offset = memory->create_2d_double_array(n+1,n+1,"pair:offset");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJClass2CoulLong::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all("Illegal 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 PairLJClass2CoulLong::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 6) error->all("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("Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJClass2CoulLong::init_one(int i, int j)
{
// always mix epsilon,sigma via sixthpower rules
// mix distance via user-defined rule
if (setflag[i][j] == 0) {
epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) *
pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) /
(pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0));
sigma[i][j] =
pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0);
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] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0);
lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0);
lj4[i][j] = 3.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] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.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];
// 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 PI = 4.0*atan(1.0);
double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
double sig6 = sig3*sig3;
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
double rc6 = rc3*rc3;
etail_ij = 2.0*PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig3 - 3.0*rc3) / (3.0*rc6);
ptail_ij = 2.0*PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig3 - 2.0*rc3) / rc6;
}
return cut;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJClass2CoulLong::init_style()
{
// require an atom style with charge defined
if (atom->charge_allow == 0)
error->all("Must use charged atom style with this pair style");
cut_coulsq = cut_coul * cut_coul;
// insure use of KSpace long-range solver, set g_ewald
if (force->kspace == NULL)
error->all("Pair style is incompatible with KSpace style");
else if (strcmp(force->kspace_style,"ewald") == 0)
g_ewald = force->kspace->g_ewald;
else if (strcmp(force->kspace_style,"pppm") == 0)
g_ewald = force->kspace->g_ewald;
else error->all("Pair style is incompatible with KSpace style");
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJClass2CoulLong::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 PairLJClass2CoulLong::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 PairLJClass2CoulLong::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 PairLJClass2CoulLong::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 PairLJClass2CoulLong::single(int i, int j, int itype, int jtype,
double rsq, double factor_coul,
double factor_lj, int eflag, One &one)
{
double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor;
double forcecoul,forcelj,phicoul,philj;
r2inv = 1.0/rsq;
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;
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
r6inv = r3inv*r3inv;
forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
} else forcelj = 0.0;
one.fforce = (forcecoul + factor_lj*forcelj) * r2inv;
if (eflag) {
if (rsq < cut_coulsq) {
phicoul = prefactor*erfc;
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
one.eng_coul = phicoul;
} else one.eng_coul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
offset[itype][jtype];
one.eng_vdwl = factor_lj*philj;
} else one.eng_vdwl = 0.0;
}
}
Event Timeline
Log In to Comment