mirror of https://github.com/lammps/lammps.git
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11954 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
a93cebcc2d
commit
153de89ab3
|
|
@ -18,7 +18,7 @@ PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \
|
|||
reax replica rigid shock srd voronoi xtc
|
||||
|
||||
PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
|
||||
user-cuda user-eff user-lb user-misc user-molfile \
|
||||
user-cuda user-eff user-fep user-lb user-misc user-molfile \
|
||||
user-omp user-phonon user-qmmm user-reaxc user-sph
|
||||
|
||||
PACKLIB = gpu kim meam poems reax voronoi \
|
||||
|
|
|
|||
|
|
@ -0,0 +1,48 @@
|
|||
# Install/unInstall package files in LAMMPS
|
||||
# mode = 0/1/2 for uninstall/install/update
|
||||
|
||||
# this is default Install.sh for all packages
|
||||
# if package has an auxiliary library or a file with a dependency,
|
||||
# then package dir has its own customized Install.sh
|
||||
|
||||
mode=$1
|
||||
|
||||
# arg1 = file, arg2 = file it depends on
|
||||
|
||||
action () {
|
||||
if (test $mode = 0) then
|
||||
rm -f ../$1
|
||||
elif (! cmp -s $1 ../$1) then
|
||||
if (test -z "$2" || test -e ../$2) then
|
||||
cp $1 ..
|
||||
if (test $mode = 2) then
|
||||
echo " updating src/$1"
|
||||
fi
|
||||
fi
|
||||
elif (test -n "$2") then
|
||||
if (test ! -e ../$2) then
|
||||
rm -f ../$1
|
||||
fi
|
||||
fi
|
||||
}
|
||||
|
||||
# all package files with dependencies
|
||||
|
||||
action compute_fep.cpp
|
||||
action compute_fep.h
|
||||
action fix_adapt_fep.cpp
|
||||
action fix_adapt_fep.h
|
||||
action pair_coul_cut_soft.cpp
|
||||
action pair_coul_cut_soft.h
|
||||
action pair_coul_long_soft.cpp pppm.cpp
|
||||
action pair_coul_long_soft.h pppm.cpp
|
||||
action pair_lj_cut_coul_cut_soft.cpp
|
||||
action pair_lj_cut_coul_cut_soft.h
|
||||
action pair_lj_cut_coul_long_soft.cpp pppm.cpp
|
||||
action pair_lj_cut_coul_long_soft.h pppm.cpp
|
||||
action pair_lj_cut_soft.cpp
|
||||
action pair_lj_cut_soft.h
|
||||
action pair_lj_cut_tip4p_long_soft.cpp pppm_tip4p.cpp
|
||||
action pair_lj_cut_tip4p_long_soft.h pppm_tip4p.cpp
|
||||
action pair_tip4p_long_soft.cpp pppm_tip4p.cpp
|
||||
action pair_tip4p_long_soft.h pppm_tip4p.cpp
|
||||
|
|
@ -0,0 +1,13 @@
|
|||
This package provides methods for performing free energy perturbation
|
||||
simulations with soft-core pair potentials in LAMMPS.
|
||||
|
||||
See the doc page for the fix adapt/fep command to get started.
|
||||
|
||||
There are example scripts for using this package in
|
||||
examples/USER/fep.
|
||||
|
||||
There are auxiliary tools for using this package in tools/fep.
|
||||
|
||||
The person who created this package is Agilio Padua at Université
|
||||
Blaise Pascal Clermont-Ferrand (agilio.padua at univ-bpclermont.fr)
|
||||
Contact him directly if you have questions.
|
||||
|
|
@ -0,0 +1,690 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "math.h"
|
||||
#include "mpi.h"
|
||||
#include "comm.h"
|
||||
#include "update.h"
|
||||
#include "atom.h"
|
||||
#include "domain.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "pair_hybrid.h"
|
||||
#include "kspace.h"
|
||||
#include "input.h"
|
||||
#include "variable.h"
|
||||
#include "timer.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "compute_fep.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
enum{PAIR,ATOM};
|
||||
enum{CHARGE};
|
||||
|
||||
#undef FEP_DEBUG
|
||||
#undef FEP_MAXDEBUG
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 5) error->all(FLERR,"Illegal number of arguments in compute fep");
|
||||
|
||||
scalar_flag = 0;
|
||||
vector_flag = 1;
|
||||
size_vector = 3;
|
||||
extvector = 0;
|
||||
|
||||
vector = new double[3];
|
||||
|
||||
fepinitflag = 0; // avoid init to run entirely when called by write_data
|
||||
|
||||
temp_fep = force->numeric(FLERR,arg[3]);
|
||||
|
||||
// count # of perturbations
|
||||
|
||||
npert = 0;
|
||||
int iarg = 4;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"pair") == 0) {
|
||||
if (iarg+6 > narg) error->all(FLERR,"Illegal pair attribute in compute fep");
|
||||
npert++;
|
||||
iarg += 6;
|
||||
} else if (strcmp(arg[iarg],"atom") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal atom attribute in compute fep");
|
||||
npert++;
|
||||
iarg += 4;
|
||||
} else break;
|
||||
}
|
||||
|
||||
if (npert == 0) error->all(FLERR,"Illegal syntax in compute fep");
|
||||
perturb = new Perturb[npert];
|
||||
|
||||
// parse keywords
|
||||
|
||||
npert = 0;
|
||||
chgflag = 0;
|
||||
|
||||
iarg = 4;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"pair") == 0) {
|
||||
perturb[npert].which = PAIR;
|
||||
int n = strlen(arg[iarg+1]) + 1;
|
||||
perturb[npert].pstyle = new char[n];
|
||||
strcpy(perturb[npert].pstyle,arg[iarg+1]);
|
||||
n = strlen(arg[iarg+2]) + 1;
|
||||
perturb[npert].pparam = new char[n];
|
||||
strcpy(perturb[npert].pparam,arg[iarg+2]);
|
||||
force->bounds(arg[iarg+3],atom->ntypes,
|
||||
perturb[npert].ilo,perturb[npert].ihi);
|
||||
force->bounds(arg[iarg+4],atom->ntypes,
|
||||
perturb[npert].jlo,perturb[npert].jhi);
|
||||
if (strstr(arg[iarg+5],"v_") == arg[iarg+5]) {
|
||||
n = strlen(&arg[iarg+5][2]) + 1;
|
||||
perturb[npert].var = new char[n];
|
||||
strcpy(perturb[npert].var,&arg[iarg+5][2]);
|
||||
} else error->all(FLERR,"Illegal variable in compute fep");
|
||||
npert++;
|
||||
iarg += 6;
|
||||
} else if (strcmp(arg[iarg],"atom") == 0) {
|
||||
perturb[npert].which = ATOM;
|
||||
if (strcmp(arg[iarg+1],"charge") == 0) {
|
||||
perturb[npert].aparam = CHARGE;
|
||||
chgflag = 1;
|
||||
} else error->all(FLERR,"Illegal atom argument in compute fep");
|
||||
force->bounds(arg[iarg+2],atom->ntypes,
|
||||
perturb[npert].ilo,perturb[npert].ihi);
|
||||
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
|
||||
int n = strlen(&arg[iarg+3][2]) + 1;
|
||||
perturb[npert].var = new char[n];
|
||||
strcpy(perturb[npert].var,&arg[iarg+3][2]);
|
||||
} else error->all(FLERR,"Illegal variable in compute fep");
|
||||
npert++;
|
||||
iarg += 4;
|
||||
} else break;
|
||||
}
|
||||
|
||||
// optional keywords
|
||||
|
||||
tailflag = 0;
|
||||
volumeflag = 0;
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"tail") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword "
|
||||
"in compute fep");
|
||||
if (strcmp(arg[iarg+1],"no") == 0) tailflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) tailflag = 1;
|
||||
else error->all(FLERR,"Illegal optional keyword in compute fep");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"volume") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword "
|
||||
"in compute fep");
|
||||
if (strcmp(arg[iarg+1],"no") == 0) volumeflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) volumeflag = 1;
|
||||
else error->all(FLERR,"Illegal optional keyword in compute fep");
|
||||
iarg += 2;
|
||||
} else
|
||||
error->all(FLERR,"Illegal optional keyword in compute fep");
|
||||
}
|
||||
|
||||
// allocate pair style arrays
|
||||
|
||||
int ntype = atom->ntypes;
|
||||
for (int m = 0; m < npert; m++) {
|
||||
if (perturb[m].which == PAIR)
|
||||
memory->create(perturb[m].array_orig,ntype+1,ntype+1,"fep:array_orig");
|
||||
}
|
||||
|
||||
// allocate space for charge, force, energy, virial arrays
|
||||
|
||||
allocate_storage();
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeFEP::~ComputeFEP()
|
||||
{
|
||||
delete [] vector;
|
||||
|
||||
for (int m = 0; m < npert; m++) {
|
||||
delete [] perturb[m].var;
|
||||
if (perturb[m].which == PAIR) {
|
||||
delete [] perturb[m].pstyle;
|
||||
delete [] perturb[m].pparam;
|
||||
memory->destroy(perturb[m].array_orig);
|
||||
}
|
||||
}
|
||||
delete [] perturb;
|
||||
|
||||
deallocate_storage();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeFEP::init()
|
||||
{
|
||||
int i,j;
|
||||
|
||||
if (!fepinitflag) // avoid init to run entirely when called by write_data
|
||||
fepinitflag = 1;
|
||||
else return;
|
||||
|
||||
// setup and error checks
|
||||
|
||||
pairflag = 0;
|
||||
|
||||
for (int m = 0; m < npert; m++) {
|
||||
Perturb *pert = &perturb[m];
|
||||
|
||||
pert->ivar = input->variable->find(pert->var);
|
||||
if (pert->ivar < 0)
|
||||
error->all(FLERR,"Variable name for compute fep does not exist");
|
||||
if (!input->variable->equalstyle(pert->ivar))
|
||||
error->all(FLERR,"Variable for compute fep is of invalid style");
|
||||
|
||||
if (force->pair == NULL)
|
||||
error->all(FLERR,"compute fep pair requires pair interactions");
|
||||
|
||||
if (pert->which == PAIR) {
|
||||
pairflag = 1;
|
||||
|
||||
Pair *pair = force->pair_match(pert->pstyle,1);
|
||||
if (pair == NULL) error->all(FLERR,"compute fep pair style "
|
||||
"does not exist");
|
||||
void *ptr = pair->extract(pert->pparam,pert->pdim);
|
||||
if (ptr == NULL)
|
||||
error->all(FLERR,"compute fep pair style param not supported");
|
||||
|
||||
pert->array = (double **) ptr;
|
||||
|
||||
// if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style
|
||||
|
||||
if ((strcmp(force->pair_style,"hybrid") == 0 ||
|
||||
strcmp(force->pair_style,"hybrid/overlay") == 0)) {
|
||||
PairHybrid *pair = (PairHybrid *) force->pair;
|
||||
for (i = pert->ilo; i <= pert->ihi; i++)
|
||||
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
|
||||
if (!pair->check_ijtype(i,j,pert->pstyle))
|
||||
error->all(FLERR,"compute fep type pair range is not valid for "
|
||||
"pair hybrid sub-style");
|
||||
}
|
||||
|
||||
} else if (pert->which == ATOM) {
|
||||
if (pert->aparam == CHARGE) {
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"compute fep requires atom attribute charge");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (tailflag) {
|
||||
if (force->pair->tail_flag == 0)
|
||||
error->all(FLERR,"Compute fep tail when pair style does not "
|
||||
"compute tail corrections");
|
||||
}
|
||||
|
||||
if (comm->me == 0) {
|
||||
if (screen) {
|
||||
fprintf(screen, "FEP settings ...\n");
|
||||
fprintf(screen, " temperature = %f\n", temp_fep);
|
||||
fprintf(screen, " tail %s\n", (tailflag ? "yes":"no"));
|
||||
for (int m = 0; m < npert; m++) {
|
||||
Perturb *pert = &perturb[m];
|
||||
if (pert->which == PAIR)
|
||||
fprintf(screen, " %s %s %d-%d %d-%d\n", pert->pstyle, pert->pparam,
|
||||
pert->ilo, pert->ihi, pert->jlo, pert->jhi);
|
||||
else if (pert->which == ATOM)
|
||||
fprintf(screen, " %d-%d charge\n", pert->ilo, pert->ihi);
|
||||
}
|
||||
}
|
||||
if (logfile) {
|
||||
fprintf(logfile, "FEP settings ...\n");
|
||||
fprintf(logfile, " temperature = %f\n", temp_fep);
|
||||
fprintf(logfile, " tail %s\n", (tailflag ? "yes":"no"));
|
||||
for (int m = 0; m < npert; m++) {
|
||||
Perturb *pert = &perturb[m];
|
||||
if (pert->which == PAIR)
|
||||
fprintf(logfile, " %s %s %d-%d %d-%d\n", pert->pstyle, pert->pparam,
|
||||
pert->ilo, pert->ihi, pert->jlo, pert->jhi);
|
||||
else if (pert->which == ATOM)
|
||||
fprintf(logfile, " %d-%d charge\n", pert->ilo, pert->ihi);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
|
||||
void ComputeFEP::compute_vector()
|
||||
{
|
||||
double pe0,pe1;
|
||||
|
||||
invoked_vector = update->ntimestep;
|
||||
|
||||
if (atom->nmax > nmax) { // reallocate working arrays if necessary
|
||||
deallocate_storage();
|
||||
allocate_storage();
|
||||
}
|
||||
|
||||
backup_qfev(); // backup charge, force, energy, virial array values
|
||||
backup_params(); // backup pair parameters
|
||||
|
||||
timer->stamp();
|
||||
if (force->pair && force->pair->compute_flag) {
|
||||
force->pair->compute(1,0);
|
||||
timer->stamp(TIME_PAIR);
|
||||
}
|
||||
if (force->kspace && force->kspace->compute_flag) {
|
||||
force->kspace->compute(1,0);
|
||||
timer->stamp(TIME_KSPACE);
|
||||
}
|
||||
pe0 = compute_epair();
|
||||
|
||||
perturb_params();
|
||||
|
||||
timer->stamp();
|
||||
if (force->pair && force->pair->compute_flag) {
|
||||
force->pair->compute(1,0);
|
||||
timer->stamp(TIME_PAIR);
|
||||
}
|
||||
if (force->kspace && force->kspace->compute_flag) {
|
||||
force->kspace->compute(1,0);
|
||||
timer->stamp(TIME_KSPACE);
|
||||
}
|
||||
pe1 = compute_epair();
|
||||
|
||||
restore_qfev(); // restore charge, force, energy, virial array values
|
||||
restore_params(); // restore pair parameters
|
||||
|
||||
vector[0] = pe1-pe0;
|
||||
vector[1] = exp(-(pe1-pe0)/(force->boltz*temp_fep));
|
||||
vector[2] = domain->xprd * domain->yprd * domain->zprd;
|
||||
if (volumeflag)
|
||||
vector[1] *= vector[2];
|
||||
|
||||
#ifdef FEP_DEBUG
|
||||
if (comm->me == 0 && screen)
|
||||
fprintf(screen, "FEP U0 = %f U1 = %f DU = %f exp(-DU/kT) = %f\n",
|
||||
pe0,pe1,vector[0],vector[1]);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
obtain pair energy from lammps accumulators
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double ComputeFEP::compute_epair()
|
||||
{
|
||||
double eng, eng_pair;
|
||||
|
||||
eng = 0.0;
|
||||
if (force->pair)
|
||||
eng = force->pair->eng_vdwl + force->pair->eng_coul;
|
||||
MPI_Allreduce(&eng,&eng_pair,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
if (tailflag) {
|
||||
double volume = domain->xprd * domain->yprd * domain->zprd;
|
||||
eng_pair += force->pair->etail / volume;
|
||||
}
|
||||
|
||||
if (force->kspace) eng_pair += force->kspace->energy;
|
||||
|
||||
return eng_pair;
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
apply perturbation to pair, atom parameters based on variable evaluation
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeFEP::perturb_params()
|
||||
{
|
||||
int i,j;
|
||||
|
||||
for (int m = 0; m < npert; m++) {
|
||||
Perturb *pert = &perturb[m];
|
||||
|
||||
double delta = input->variable->compute_equal(pert->ivar);
|
||||
|
||||
if (pert->which == PAIR) { // modify pair parameters
|
||||
for (i = pert->ilo; i <= pert->ihi; i++)
|
||||
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
|
||||
pert->array[i][j] = pert->array_orig[i][j] + delta;
|
||||
|
||||
#ifdef FEP_MAXDEBUG
|
||||
if (comm->me == 0 && screen) {
|
||||
fprintf(screen, "###FEP change %s %s, delta = %f\n",
|
||||
pert->pstyle, pert->pparam, delta);
|
||||
fprintf(screen, "###FEP I J old_param new_param\n");
|
||||
for (i = pert->ilo; i <= pert->ihi; i++)
|
||||
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
|
||||
fprintf(screen, "###FEP %2d %2d %9.5f %9.5f\n", i, j,
|
||||
pert->array_orig[i][j], pert->array[i][j]);
|
||||
}
|
||||
#endif
|
||||
|
||||
} else if (pert->which == ATOM) {
|
||||
|
||||
if (pert->aparam == CHARGE) { // modify charges
|
||||
int *atype = atom->type;
|
||||
double *q = atom->q;
|
||||
int *mask = atom->mask;
|
||||
int natom = atom->nlocal + atom->nghost;
|
||||
|
||||
for (i = 0; i < natom; i++)
|
||||
if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
|
||||
if (mask[i] & groupbit)
|
||||
q[i] += delta;
|
||||
|
||||
#ifdef FEP_MAXDEBUG
|
||||
if (comm->me == 0 && screen) {
|
||||
fprintf(screen, "###FEP change charge, delta = %f\n", delta);
|
||||
fprintf(screen, "###FEP atom I old_q new_q\n");
|
||||
for (i = 0; i < atom->nlocal; i++)
|
||||
if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
|
||||
if (mask[i] & groupbit)
|
||||
fprintf(screen, "###FEP %5d %2d %9.5f %9.5f\n", i, atype[i],
|
||||
q_orig[i], q[i]);
|
||||
}
|
||||
#endif
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// re-initialize pair styles if any PAIR settings were changed
|
||||
// this resets other coeffs that may depend on changed values,
|
||||
// and also offset and tail corrections
|
||||
|
||||
if (pairflag) force->pair->reinit();
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
backup pair parameters
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeFEP::backup_params()
|
||||
{
|
||||
int i,j;
|
||||
|
||||
for (int m = 0; m < npert; m++) {
|
||||
Perturb *pert = &perturb[m];
|
||||
if (pert->which == PAIR) {
|
||||
for (i = pert->ilo; i <= pert->ihi; i++)
|
||||
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
|
||||
pert->array_orig[i][j] = pert->array[i][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
restore pair parameters to original values
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeFEP::restore_params()
|
||||
{
|
||||
int i,j;
|
||||
|
||||
for (int m = 0; m < npert; m++) {
|
||||
Perturb *pert = &perturb[m];
|
||||
if (pert->which == PAIR) {
|
||||
for (i = pert->ilo; i <= pert->ihi; i++)
|
||||
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
|
||||
pert->array[i][j] = pert->array_orig[i][j];
|
||||
|
||||
#ifdef FEP_MAXDEBUG
|
||||
if (comm->me == 0 && screen) {
|
||||
fprintf(screen, "###FEP restore %s %s\n", pert->pstyle, pert->pparam);
|
||||
fprintf(screen, "###FEP I J param\n");
|
||||
for (i = pert->ilo; i <= pert->ihi; i++)
|
||||
for (j = MAX(pert->jlo,i); j <= pert->jhi; j++)
|
||||
fprintf(screen, "###FEP %2d %2d %9.5f\n", i, j, pert->array[i][j]);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
#ifdef FEP_MAXDEBUG
|
||||
if (pert->which == ATOM) {
|
||||
if (comm->me == 0 && screen) {
|
||||
int *atype = atom->type;
|
||||
double *q = atom->q;
|
||||
int *mask = atom->mask;
|
||||
int natom = atom->nlocal;
|
||||
fprintf(screen, "###FEP restore charge\n");
|
||||
fprintf(screen, "###FEP atom I q\n");
|
||||
for (i = 0; i < natom; i++)
|
||||
if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
|
||||
if (mask[i] & groupbit)
|
||||
fprintf(screen, "###FEP %5d %2d %9.5f\n", i, atype[i], q[i]);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
// re-initialize pair styles if any PAIR settings were changed
|
||||
// this resets other coeffs that may depend on changed values,
|
||||
// and also offset and tail corrections
|
||||
|
||||
if (pairflag) force->pair->reinit();
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
manage storage for charge, force, energy, virial arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeFEP::allocate_storage()
|
||||
{
|
||||
nmax = atom->nmax;
|
||||
if (chgflag)
|
||||
memory->create(q_orig,nmax,"fep:q_orig");
|
||||
memory->create(f_orig,nmax,3,"fep:f_orig");
|
||||
memory->create(peatom_orig,nmax,"fep:peatom_orig");
|
||||
memory->create(pvatom_orig,nmax,6,"fep:pvatom_orig");
|
||||
if (force->kspace) {
|
||||
memory->create(keatom_orig,nmax,"fep:keatom_orig");
|
||||
memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig");
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeFEP::deallocate_storage()
|
||||
{
|
||||
if (chgflag)
|
||||
memory->destroy(q_orig);
|
||||
memory->destroy(f_orig);
|
||||
memory->destroy(peatom_orig);
|
||||
memory->destroy(pvatom_orig);
|
||||
if (force->kspace) {
|
||||
memory->destroy(keatom_orig);
|
||||
memory->destroy(kvatom_orig);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
backup and restore arrays with charge, force, energy, virial
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeFEP::backup_qfev()
|
||||
{
|
||||
int i;
|
||||
|
||||
int nall = atom->nlocal + atom->nghost;
|
||||
int natom = atom->nlocal;
|
||||
if (force->newton || force->kspace->tip4pflag)
|
||||
natom += atom->nghost;
|
||||
|
||||
if (chgflag) {
|
||||
double *q = atom->q;
|
||||
for (i = 0; i < nall; i++)
|
||||
q_orig[i] = q[i];
|
||||
}
|
||||
|
||||
double **f = atom->f;
|
||||
for (i = 0; i < natom; i++) {
|
||||
f_orig[i][0] = f[i][0];
|
||||
f_orig[i][1] = f[i][1];
|
||||
f_orig[i][2] = f[i][2];
|
||||
}
|
||||
|
||||
eng_vdwl_orig = force->pair->eng_vdwl;
|
||||
eng_coul_orig = force->pair->eng_coul;
|
||||
|
||||
pvirial_orig[0] = force->pair->virial[0];
|
||||
pvirial_orig[1] = force->pair->virial[1];
|
||||
pvirial_orig[2] = force->pair->virial[2];
|
||||
pvirial_orig[3] = force->pair->virial[3];
|
||||
pvirial_orig[4] = force->pair->virial[4];
|
||||
pvirial_orig[5] = force->pair->virial[5];
|
||||
|
||||
if (update->eflag_atom) {
|
||||
double *peatom = force->pair->eatom;
|
||||
for (i = 0; i < natom; i++)
|
||||
peatom_orig[i] = peatom[i];
|
||||
}
|
||||
if (update->vflag_atom) {
|
||||
double **pvatom = force->pair->vatom;
|
||||
for (i = 0; i < natom; i++) {
|
||||
pvatom_orig[i][0] = pvatom[i][0];
|
||||
pvatom_orig[i][1] = pvatom[i][1];
|
||||
pvatom_orig[i][2] = pvatom[i][2];
|
||||
pvatom_orig[i][3] = pvatom[i][3];
|
||||
pvatom_orig[i][4] = pvatom[i][4];
|
||||
pvatom_orig[i][5] = pvatom[i][5];
|
||||
}
|
||||
}
|
||||
|
||||
if (force->kspace) {
|
||||
energy_orig = force->kspace->energy;
|
||||
kvirial_orig[0] = force->kspace->virial[0];
|
||||
kvirial_orig[1] = force->kspace->virial[1];
|
||||
kvirial_orig[2] = force->kspace->virial[2];
|
||||
kvirial_orig[3] = force->kspace->virial[3];
|
||||
kvirial_orig[4] = force->kspace->virial[4];
|
||||
kvirial_orig[5] = force->kspace->virial[5];
|
||||
|
||||
if (update->eflag_atom) {
|
||||
double *keatom = force->kspace->eatom;
|
||||
for (i = 0; i < natom; i++)
|
||||
keatom_orig[i] = keatom[i];
|
||||
}
|
||||
if (update->vflag_atom) {
|
||||
double **kvatom = force->kspace->vatom;
|
||||
for (i = 0; i < natom; i++) {
|
||||
kvatom_orig[i][0] = kvatom[i][0];
|
||||
kvatom_orig[i][1] = kvatom[i][1];
|
||||
kvatom_orig[i][2] = kvatom[i][2];
|
||||
kvatom_orig[i][3] = kvatom[i][3];
|
||||
kvatom_orig[i][4] = kvatom[i][4];
|
||||
kvatom_orig[i][5] = kvatom[i][5];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeFEP::restore_qfev()
|
||||
{
|
||||
int i;
|
||||
|
||||
int nall = atom->nlocal + atom->nghost;
|
||||
int natom = atom->nlocal;
|
||||
if (force->newton || force->kspace->tip4pflag)
|
||||
natom += atom->nghost;
|
||||
|
||||
if (chgflag) {
|
||||
double *q = atom->q;
|
||||
for (i = 0; i < nall; i++)
|
||||
q[i] = q_orig[i];
|
||||
}
|
||||
|
||||
double **f = atom->f;
|
||||
for (i = 0; i < natom; i++) {
|
||||
f[i][0] = f_orig[i][0];
|
||||
f[i][1] = f_orig[i][1];
|
||||
f[i][2] = f_orig[i][2];
|
||||
}
|
||||
|
||||
force->pair->eng_vdwl = eng_vdwl_orig;
|
||||
force->pair->eng_coul = eng_coul_orig;
|
||||
|
||||
force->pair->virial[0] = pvirial_orig[0];
|
||||
force->pair->virial[1] = pvirial_orig[1];
|
||||
force->pair->virial[2] = pvirial_orig[2];
|
||||
force->pair->virial[3] = pvirial_orig[3];
|
||||
force->pair->virial[4] = pvirial_orig[4];
|
||||
force->pair->virial[5] = pvirial_orig[5];
|
||||
|
||||
if (update->eflag_atom) {
|
||||
double *peatom = force->pair->eatom;
|
||||
for (i = 0; i < natom; i++)
|
||||
peatom[i] = peatom_orig[i];
|
||||
}
|
||||
if (update->vflag_atom) {
|
||||
double **pvatom = force->pair->vatom;
|
||||
for (i = 0; i < natom; i++) {
|
||||
pvatom[i][0] = pvatom_orig[i][0];
|
||||
pvatom[i][1] = pvatom_orig[i][1];
|
||||
pvatom[i][2] = pvatom_orig[i][2];
|
||||
pvatom[i][3] = pvatom_orig[i][3];
|
||||
pvatom[i][4] = pvatom_orig[i][4];
|
||||
pvatom[i][5] = pvatom_orig[i][5];
|
||||
}
|
||||
}
|
||||
|
||||
if (force->kspace) {
|
||||
force->kspace->energy = energy_orig;
|
||||
force->kspace->virial[0] = kvirial_orig[0];
|
||||
force->kspace->virial[1] = kvirial_orig[1];
|
||||
force->kspace->virial[2] = kvirial_orig[2];
|
||||
force->kspace->virial[3] = kvirial_orig[3];
|
||||
force->kspace->virial[4] = kvirial_orig[4];
|
||||
force->kspace->virial[5] = kvirial_orig[5];
|
||||
|
||||
if (update->eflag_atom) {
|
||||
double *keatom = force->kspace->eatom;
|
||||
for (i = 0; i < natom; i++)
|
||||
keatom[i] = keatom_orig[i];
|
||||
}
|
||||
if (update->vflag_atom) {
|
||||
double **kvatom = force->kspace->vatom;
|
||||
for (i = 0; i < natom; i++) {
|
||||
kvatom[i][0] = kvatom_orig[i][0];
|
||||
kvatom[i][1] = kvatom_orig[i][1];
|
||||
kvatom[i][2] = kvatom_orig[i][2];
|
||||
kvatom[i][3] = kvatom_orig[i][3];
|
||||
kvatom[i][4] = kvatom_orig[i][4];
|
||||
kvatom[i][5] = kvatom_orig[i][5];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -0,0 +1,109 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Agilio Padua (ICCF,UBP,CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
|
||||
ComputeStyle(fep,ComputeFEP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef COMPUTE_FEP_H
|
||||
#define COMPUTE_FEP_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeFEP : public Compute {
|
||||
public:
|
||||
ComputeFEP(class LAMMPS *, int, char **);
|
||||
~ComputeFEP();
|
||||
void init();
|
||||
void compute_vector();
|
||||
|
||||
private:
|
||||
int npert;
|
||||
int pairflag;
|
||||
int chgflag;
|
||||
int tailflag, volumeflag;
|
||||
int fepinitflag;
|
||||
double temp_fep;
|
||||
|
||||
int nmax;
|
||||
double *q_orig;
|
||||
double **f_orig;
|
||||
double eng_vdwl_orig,eng_coul_orig;
|
||||
double pvirial_orig[6];
|
||||
double *peatom_orig,**pvatom_orig;
|
||||
double energy_orig;
|
||||
double kvirial_orig[6];
|
||||
double *keatom_orig,**kvatom_orig;
|
||||
|
||||
struct Perturb {
|
||||
int which,ivar;
|
||||
char *var;
|
||||
char *pstyle,*pparam;
|
||||
int ilo,ihi,jlo,jhi;
|
||||
int pdim;
|
||||
double **array,**array_orig;
|
||||
int aparam;
|
||||
};
|
||||
|
||||
Perturb *perturb;
|
||||
|
||||
double compute_epair();
|
||||
void perturb_params();
|
||||
void backup_params();
|
||||
void restore_params();
|
||||
void allocate_storage();
|
||||
void deallocate_storage();
|
||||
void backup_qfev();
|
||||
void restore_qfev();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#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: Variable name for compute fep does not exist
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Variable for compute fep is invalid style
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Compute fep pair style does not exist
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Energy was not tallied on needed timestep
|
||||
|
||||
You are using a thermo keyword that requires potentials to
|
||||
have tallied energy, but they didn't on this timestep. See the
|
||||
variable doc page for ideas on how to make this work.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,460 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Charges by type and after option: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "string.h"
|
||||
#include "stdlib.h"
|
||||
#include "fix_adapt_fep.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "pair_hybrid.h"
|
||||
#include "kspace.h"
|
||||
#include "input.h"
|
||||
#include "variable.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
using namespace MathConst;
|
||||
|
||||
enum{PAIR,KSPACE,ATOM};
|
||||
enum{DIAMETER,CHARGE};
|
||||
|
||||
#undef ADAPT_DEBUG
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 5) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
nevery = force->inumeric(FLERR,arg[3]);
|
||||
if (nevery < 0) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
|
||||
// count # of adaptations
|
||||
|
||||
nadapt = 0;
|
||||
|
||||
int iarg = 4;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"pair") == 0) {
|
||||
if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
nadapt++;
|
||||
iarg += 6;
|
||||
} else if (strcmp(arg[iarg],"kspace") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
nadapt++;
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"atom") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
nadapt++;
|
||||
iarg += 4;
|
||||
} else break;
|
||||
}
|
||||
|
||||
if (nadapt == 0) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
adapt = new Adapt[nadapt];
|
||||
|
||||
// parse keywords
|
||||
|
||||
nadapt = 0;
|
||||
diamflag = 0;
|
||||
chgflag = 0;
|
||||
|
||||
iarg = 4;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"pair") == 0) {
|
||||
if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
adapt[nadapt].which = PAIR;
|
||||
int n = strlen(arg[iarg+1]) + 1;
|
||||
adapt[nadapt].pstyle = new char[n];
|
||||
strcpy(adapt[nadapt].pstyle,arg[iarg+1]);
|
||||
n = strlen(arg[iarg+2]) + 1;
|
||||
adapt[nadapt].pparam = new char[n];
|
||||
strcpy(adapt[nadapt].pparam,arg[iarg+2]);
|
||||
force->bounds(arg[iarg+3],atom->ntypes,
|
||||
adapt[nadapt].ilo,adapt[nadapt].ihi);
|
||||
force->bounds(arg[iarg+4],atom->ntypes,
|
||||
adapt[nadapt].jlo,adapt[nadapt].jhi);
|
||||
if (strstr(arg[iarg+5],"v_") == arg[iarg+5]) {
|
||||
n = strlen(&arg[iarg+5][2]) + 1;
|
||||
adapt[nadapt].var = new char[n];
|
||||
strcpy(adapt[nadapt].var,&arg[iarg+5][2]);
|
||||
} else error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
nadapt++;
|
||||
iarg += 6;
|
||||
} else if (strcmp(arg[iarg],"kspace") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
adapt[nadapt].which = KSPACE;
|
||||
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
|
||||
int n = strlen(&arg[iarg+1][2]) + 1;
|
||||
adapt[nadapt].var = new char[n];
|
||||
strcpy(adapt[nadapt].var,&arg[iarg+1][2]);
|
||||
} else error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
nadapt++;
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"atom") == 0) {
|
||||
if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
adapt[nadapt].which = ATOM;
|
||||
if (strcmp(arg[iarg+1],"diameter") == 0) {
|
||||
adapt[nadapt].aparam = DIAMETER;
|
||||
diamflag = 1;
|
||||
} else if (strcmp(arg[iarg+1],"charge") == 0) {
|
||||
adapt[nadapt].aparam = CHARGE;
|
||||
chgflag = 1;
|
||||
} else error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
force->bounds(arg[iarg+2],atom->ntypes,
|
||||
adapt[nadapt].ilo,adapt[nadapt].ihi);
|
||||
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
|
||||
int n = strlen(&arg[iarg+3][2]) + 1;
|
||||
adapt[nadapt].var = new char[n];
|
||||
strcpy(adapt[nadapt].var,&arg[iarg+3][2]);
|
||||
} else error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
nadapt++;
|
||||
iarg += 4;
|
||||
} else break;
|
||||
}
|
||||
|
||||
// optional keywords
|
||||
|
||||
resetflag = 0;
|
||||
scaleflag = 0;
|
||||
afterflag = 0;
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"reset") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
if (strcmp(arg[iarg+1],"no") == 0) resetflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) resetflag = 1;
|
||||
else error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"scale") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1;
|
||||
else error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"after") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
if (strcmp(arg[iarg+1],"no") == 0) afterflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) afterflag = 1;
|
||||
else error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
iarg += 2;
|
||||
} else error->all(FLERR,"Illegal fix adapt/fep command");
|
||||
}
|
||||
|
||||
// allocate pair style arrays
|
||||
|
||||
int n = atom->ntypes;
|
||||
for (int m = 0; m < nadapt; m++)
|
||||
if (adapt[m].which == PAIR)
|
||||
memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixAdaptFEP::~FixAdaptFEP()
|
||||
{
|
||||
for (int m = 0; m < nadapt; m++) {
|
||||
delete [] adapt[m].var;
|
||||
if (adapt[m].which == PAIR) {
|
||||
delete [] adapt[m].pstyle;
|
||||
delete [] adapt[m].pparam;
|
||||
memory->destroy(adapt[m].array_orig);
|
||||
}
|
||||
}
|
||||
delete [] adapt;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixAdaptFEP::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= PRE_FORCE;
|
||||
mask |= POST_RUN;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixAdaptFEP::init()
|
||||
{
|
||||
int i,j;
|
||||
|
||||
// setup and error checks
|
||||
|
||||
anypair = 0;
|
||||
|
||||
for (int m = 0; m < nadapt; m++) {
|
||||
Adapt *ad = &adapt[m];
|
||||
|
||||
ad->ivar = input->variable->find(ad->var);
|
||||
if (ad->ivar < 0)
|
||||
error->all(FLERR,"Variable name for fix adapt/fep does not exist");
|
||||
if (!input->variable->equalstyle(ad->ivar))
|
||||
error->all(FLERR,"Variable for fix adapt/fep is invalid style");
|
||||
|
||||
if (ad->which == PAIR) {
|
||||
anypair = 1;
|
||||
|
||||
Pair *pair = force->pair_match(ad->pstyle,1);
|
||||
if (pair == NULL) error->all(FLERR,"Fix adapt/fep pair style does not exist");
|
||||
void *ptr = pair->extract(ad->pparam,ad->pdim);
|
||||
if (ptr == NULL)
|
||||
error->all(FLERR,"Fix adapt/fep pair style param not supported");
|
||||
|
||||
ad->pdim = 2;
|
||||
if (ad->pdim == 0) ad->scalar = (double *) ptr;
|
||||
if (ad->pdim == 2) ad->array = (double **) ptr;
|
||||
|
||||
// if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style
|
||||
|
||||
if (ad->pdim == 2 && (strcmp(force->pair_style,"hybrid") == 0 ||
|
||||
strcmp(force->pair_style,"hybrid/overlay") == 0)) {
|
||||
PairHybrid *pair = (PairHybrid *) force->pair;
|
||||
for (i = ad->ilo; i <= ad->ihi; i++)
|
||||
for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
|
||||
if (!pair->check_ijtype(i,j,ad->pstyle))
|
||||
error->all(FLERR,"Fix adapt/fep type pair range is not valid for "
|
||||
"pair hybrid sub-style");
|
||||
}
|
||||
|
||||
} else if (ad->which == KSPACE) {
|
||||
if (force->kspace == NULL)
|
||||
error->all(FLERR,"Fix adapt/fep kspace style does not exist");
|
||||
kspace_scale = (double *) force->kspace->extract("scale");
|
||||
|
||||
} else if (ad->which == ATOM) {
|
||||
if (ad->aparam == DIAMETER) {
|
||||
if (!atom->radius_flag)
|
||||
error->all(FLERR,"Fix adapt/fep requires atom attribute diameter");
|
||||
}
|
||||
if (ad->aparam == CHARGE) {
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Fix adapt/fep requires atom attribute charge");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// make copy of original pair array values
|
||||
|
||||
for (int m = 0; m < nadapt; m++) {
|
||||
Adapt *ad = &adapt[m];
|
||||
if (ad->which == PAIR && ad->pdim == 2) {
|
||||
for (i = ad->ilo; i <= ad->ihi; i++)
|
||||
for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
|
||||
ad->array_orig[i][j] = ad->array[i][j];
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef ADAPT_DEBUG
|
||||
if (comm->me == 0 && screen) {
|
||||
for (int m = 0; m < nadapt; m++) {
|
||||
Adapt *ad = &adapt[m];
|
||||
if (ad->which == PAIR && ad->pdim == 2) {
|
||||
fprintf(screen, "###ADAPT original %s %s\n", ad->pstyle, ad->pparam);
|
||||
fprintf(screen, "###ADAPT I J old_param\n");
|
||||
for (i = ad->ilo; i <= ad->ihi; i++)
|
||||
for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
|
||||
fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j,
|
||||
ad->array_orig[i][j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixAdaptFEP::setup_pre_force(int vflag)
|
||||
{
|
||||
change_settings();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixAdaptFEP::pre_force(int vflag)
|
||||
{
|
||||
if (nevery == 0) return;
|
||||
|
||||
if (afterflag) { // update at n+1 (better with fix ave/time)
|
||||
if (nevery == 1 || update->ntimestep == 0)
|
||||
change_settings();
|
||||
else if (update->ntimestep > 1 && !((update->ntimestep - 1) % nevery))
|
||||
change_settings();
|
||||
} else { // original version: update at n
|
||||
if (update->ntimestep % nevery)
|
||||
return;
|
||||
change_settings();
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixAdaptFEP::post_run()
|
||||
{
|
||||
if (resetflag) restore_settings();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
change pair,kspace,atom parameters based on variable evaluation
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixAdaptFEP::change_settings()
|
||||
{
|
||||
int i,j;
|
||||
|
||||
// variable evaluation may invoke computes so wrap with clear/add
|
||||
|
||||
modify->clearstep_compute();
|
||||
|
||||
for (int m = 0; m < nadapt; m++) {
|
||||
Adapt *ad = &adapt[m];
|
||||
double value = input->variable->compute_equal(ad->ivar);
|
||||
|
||||
// set global scalar or type pair array values
|
||||
|
||||
if (ad->which == PAIR) {
|
||||
if (ad->pdim == 0) {
|
||||
if (scaleflag) *ad->scalar = value * ad->scalar_orig;
|
||||
else *ad->scalar = value;
|
||||
} else if (ad->pdim == 2) {
|
||||
if (scaleflag)
|
||||
for (i = ad->ilo; i <= ad->ihi; i++)
|
||||
for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
|
||||
ad->array[i][j] = value*ad->array_orig[i][j];
|
||||
else
|
||||
for (i = ad->ilo; i <= ad->ihi; i++)
|
||||
for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
|
||||
ad->array[i][j] = value;
|
||||
|
||||
#ifdef ADAPT_DEBUG
|
||||
if (comm->me == 0 && screen) {
|
||||
fprintf(screen, "###ADAPT change %s %s\n", ad->pstyle, ad->pparam);
|
||||
fprintf(screen, "###ADAPT I J param\n");
|
||||
for (i = ad->ilo; i <= ad->ihi; i++)
|
||||
for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
|
||||
fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j, ad->array[i][j]);
|
||||
}
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
// set kspace scale factor
|
||||
|
||||
} else if (ad->which == KSPACE) {
|
||||
*kspace_scale = value;
|
||||
|
||||
} else if (ad->which == ATOM) {
|
||||
|
||||
// set radius from diameter
|
||||
// also scale rmass to new value
|
||||
|
||||
if (ad->aparam == DIAMETER) {
|
||||
int mflag = 0;
|
||||
if (atom->rmass_flag) mflag = 1;
|
||||
double density;
|
||||
|
||||
int *atype = atom->type;
|
||||
double *radius = atom->radius;
|
||||
double *rmass = atom->rmass;
|
||||
int *mask = atom->mask;
|
||||
int natom = atom->nlocal + atom->nghost;
|
||||
|
||||
if (mflag == 0) {
|
||||
for (i = 0; i < natom; i++)
|
||||
if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
|
||||
if (mask[i] & groupbit)
|
||||
radius[i] = 0.5*value;
|
||||
} else {
|
||||
for (i = 0; i < natom; i++)
|
||||
if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
|
||||
if (mask[i] & groupbit) {
|
||||
density = rmass[i] / (4.0*MY_PI/3.0 *
|
||||
radius[i]*radius[i]*radius[i]);
|
||||
radius[i] = 0.5*value;
|
||||
rmass[i] = 4.0*MY_PI/3.0 *
|
||||
radius[i]*radius[i]*radius[i] * density;
|
||||
}
|
||||
}
|
||||
} else if (ad->aparam == CHARGE) {
|
||||
int *atype = atom->type;
|
||||
double *q = atom->q;
|
||||
int *mask = atom->mask;
|
||||
int natom = atom->nlocal + atom->nghost;
|
||||
|
||||
for (i = 0; i < natom; i++)
|
||||
if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
|
||||
if (mask[i] & groupbit) q[i] = value;
|
||||
|
||||
#ifdef ADAPT_DEBUG
|
||||
if (comm->me == 0 && screen) {
|
||||
fprintf(screen, "###ADAPT change charge\n");
|
||||
fprintf(screen, "###ADAPT atom I q\n");
|
||||
for (i = 0; i < atom->nlocal; i++)
|
||||
if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
|
||||
if (mask[i] & groupbit)
|
||||
fprintf(screen, "###ADAPT %5d %2d %9.5f\n", i, atype[i], q[i]);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
modify->addstep_compute(update->ntimestep + nevery);
|
||||
|
||||
// re-initialize pair styles if any PAIR settings were changed
|
||||
// this resets other coeffs that may depend on changed values,
|
||||
// and also offset and tail corrections
|
||||
|
||||
if (anypair) force->pair->reinit();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
restore pair,kspace.atom parameters to original values
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixAdaptFEP::restore_settings()
|
||||
{
|
||||
for (int m = 0; m < nadapt; m++) {
|
||||
Adapt *ad = &adapt[m];
|
||||
if (ad->which == PAIR) {
|
||||
if (ad->pdim == 0) *ad->scalar = ad->scalar_orig;
|
||||
else if (ad->pdim == 2) {
|
||||
for (int i = ad->ilo; i <= ad->ihi; i++)
|
||||
for (int j = MAX(ad->jlo,i); j <= ad->jhi; j++)
|
||||
ad->array[i][j] = ad->array_orig[i][j];
|
||||
}
|
||||
|
||||
} else if (ad->which == KSPACE) {
|
||||
*kspace_scale = 1.0;
|
||||
|
||||
} else if (ad->which == ATOM) {
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
if (anypair) force->pair->reinit();
|
||||
}
|
||||
|
|
@ -0,0 +1,107 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 FIX_CLASS
|
||||
|
||||
FixStyle(adapt/fep,FixAdaptFEP)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_ADAPT_FEP_H
|
||||
#define LMP_FIX_ADAPT_FEP_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixAdaptFEP : public Fix {
|
||||
public:
|
||||
int diamflag; // 1 if atom diameters will vary, for AtomVecGranular
|
||||
int chgflag;
|
||||
|
||||
FixAdaptFEP(class LAMMPS *, int, char **);
|
||||
~FixAdaptFEP();
|
||||
int setmask();
|
||||
void init();
|
||||
void setup_pre_force(int);
|
||||
void pre_force(int);
|
||||
void post_run();
|
||||
|
||||
private:
|
||||
int nadapt,resetflag,scaleflag,afterflag;
|
||||
int anypair;
|
||||
|
||||
struct Adapt {
|
||||
int which,ivar;
|
||||
char *var;
|
||||
char *pstyle,*pparam;
|
||||
int ilo,ihi,jlo,jhi;
|
||||
int pdim;
|
||||
double *scalar,scalar_orig;
|
||||
double **array,**array_orig;
|
||||
int aparam;
|
||||
};
|
||||
|
||||
Adapt *adapt;
|
||||
double *kspace_scale;
|
||||
|
||||
void change_settings();
|
||||
void restore_settings();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#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: Variable name for fix adapt does not exist
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Variable for fix adapt is invalid style
|
||||
|
||||
Only equal-style variables can be used.
|
||||
|
||||
E: Fix adapt pair style does not exist
|
||||
|
||||
Self-explanatory
|
||||
|
||||
E: Fix adapt pair style param not supported
|
||||
|
||||
The pair style does not know about the parameter you specified.
|
||||
|
||||
E: Fix adapt type pair range is not valid for pair hybrid sub-style
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix adapt kspace style does not exist
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix adapt requires atom attribute diameter
|
||||
|
||||
The atom style being used does not specify an atom diameter.
|
||||
|
||||
E: Fix adapt requires atom attribute charge
|
||||
|
||||
The atom style being used does not specify an atom charge.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,376 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_coul_cut_soft.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;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulCutSoft::PairCoulCutSoft(LAMMPS *lmp) : Pair(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulCutSoft::~PairCoulCutSoft()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(cut);
|
||||
memory->destroy(lambda);
|
||||
memory->destroy(lam1);
|
||||
memory->destroy(lam2);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
double rsq,forcecoul,factor_coul;
|
||||
double denc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
|
||||
denc = sqrt(lam2[itype][jtype] + rsq);
|
||||
forcecoul = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
|
||||
fpair = factor_coul*forcecoul;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag)
|
||||
ecoul = factor_coul * qqrd2e * lam1[itype][jtype] * qtmp*q[j] / denc;
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
0.0,ecoul,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::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(lambda,n+1,n+1,"pair:lambda");
|
||||
memory->create(lam1,n+1,n+1,"pair:lam1");
|
||||
memory->create(lam2,n+1,n+1,"pair:lam2");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
nlambda = force->numeric(FLERR,arg[0]);
|
||||
alphac = force->numeric(FLERR,arg[1]);
|
||||
|
||||
cut_global = force->numeric(FLERR,arg[2]);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i+1; j <= atom->ntypes; j++)
|
||||
if (setflag[i][j]) cut[i][j] = cut_global;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 3 || narg > 4)
|
||||
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 lambda_one = force->numeric(FLERR,arg[2]);
|
||||
|
||||
double cut_one = cut_global;
|
||||
if (narg == 4) cut_one = force->numeric(FLERR,arg[3]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
lambda[i][j] = lambda_one;
|
||||
cut[i][j] = cut_one;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::init_style()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Pair style coul/cut/soft requires atom attribute q");
|
||||
|
||||
neighbor->request(this);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulCutSoft::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) {
|
||||
if (lambda[i][i] != lambda[j][j])
|
||||
error->all(FLERR,"Pair coul/cut/soft different lambda values in mix");
|
||||
lambda[i][j] = lambda[i][i];
|
||||
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
|
||||
}
|
||||
|
||||
lam1[i][j] = pow(lambda[i][j], nlambda);
|
||||
lam2[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
|
||||
|
||||
cut[j][i] = cut[i][j];
|
||||
lambda[j][i] = lambda[i][j];
|
||||
lam1[j][i] = lam1[i][j];
|
||||
lam2[j][i] = lam2[i][j];
|
||||
|
||||
return cut[i][j];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::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(&lambda[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::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(&lambda[i][j],sizeof(double),1,fp);
|
||||
fread(&cut[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&nlambda,sizeof(double),1,fp);
|
||||
fwrite(&alphac,sizeof(double),1,fp);
|
||||
|
||||
fwrite(&cut_global,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&nlambda,sizeof(double),1,fp);
|
||||
fread(&alphac,sizeof(double),1,fp);
|
||||
|
||||
fread(&cut_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alphac,1,MPI_DOUBLE,0,world);
|
||||
|
||||
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
fprintf(fp,"%d %g\n",i,lambda[i][i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes all pairs to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCutSoft::write_data_all(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
for (int j = i; j <= atom->ntypes; j++)
|
||||
fprintf(fp,"%d %d %g\n",i,j,lambda[i][j]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulCutSoft::single(int i, int j, int itype, int jtype,
|
||||
double rsq, double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double forcecoul,phicoul;
|
||||
double denc;
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
denc = sqrt(lam2[itype][jtype] + rsq);
|
||||
forcecoul = force->qqrd2e * lam1[itype][jtype] * atom->q[i]*atom->q[j] /
|
||||
(denc*denc*denc);
|
||||
} else forcecoul = 0.0;
|
||||
fforce = factor_coul*forcecoul;
|
||||
|
||||
if (rsq < cutsq[itype][jtype])
|
||||
phicoul = force->qqrd2e * lam1[itype][jtype] * atom->q[i]*atom->q[j] / denc;
|
||||
else phicoul = 0.0;
|
||||
return factor_coul*phicoul;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairCoulCutSoft::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 2;
|
||||
if (strcmp(str,"lambda") == 0) return (void *) lambda;
|
||||
return NULL;
|
||||
}
|
||||
|
|
@ -0,0 +1,80 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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(coul/cut/soft,PairCoulCutSoft)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_COUL_CUT_SOFT_H
|
||||
#define LMP_PAIR_COUL_CUT_SOFT_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairCoulCutSoft : public Pair {
|
||||
public:
|
||||
PairCoulCutSoft(class LAMMPS *);
|
||||
virtual ~PairCoulCutSoft();
|
||||
virtual void compute(int, int);
|
||||
virtual void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
virtual void write_restart_settings(FILE *);
|
||||
virtual void read_restart_settings(FILE *);
|
||||
void write_data(FILE *);
|
||||
void write_data_all(FILE *);
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
void *extract(const char *, int &);
|
||||
|
||||
protected:
|
||||
double cut_global;
|
||||
double **cut;
|
||||
double **lambda;
|
||||
double nlambda, alphac;
|
||||
double **lam1, **lam2;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair style coul/cut/soft requires atom attribute q
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
E: Pair coul/cut/soft different lambda values in mix
|
||||
|
||||
The value of lambda has to be the same for I J pairs.
|
||||
|
||||
*/
|
||||
|
|
@ -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: Paul Crozier (SNL)
|
||||
Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_coul_long_soft.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulLongSoft::PairCoulLongSoft(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
ewaldflag = pppmflag = 1;
|
||||
qdist = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulLongSoft::~PairCoulLongSoft()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(scale);
|
||||
|
||||
memory->destroy(lambda);
|
||||
memory->destroy(lam1);
|
||||
memory->destroy(lam2);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulLongSoft::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
double r,rsq,forcecoul,factor_coul;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double denc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < 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;
|
||||
|
||||
denc = sqrt(lam2[itype][jtype] + rsq);
|
||||
prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
fpair = forcecoul;
|
||||
|
||||
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) {
|
||||
prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / denc;
|
||||
ecoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
0.0,ecoul,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulLongSoft::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(scale,n+1,n+1,"pair:scale");
|
||||
|
||||
memory->create(lambda,n+1,n+1,"pair:lambda");
|
||||
memory->create(lam1,n+1,n+1,"pair:lam1");
|
||||
memory->create(lam2,n+1,n+1,"pair:lam2");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulLongSoft::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
nlambda = force->numeric(FLERR,arg[0]);
|
||||
alphac = force->numeric(FLERR,arg[1]);
|
||||
|
||||
cut_coul = force->numeric(FLERR,arg[2]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulLongSoft::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3)
|
||||
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 lambda_one = force->numeric(FLERR,arg[2]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
lambda[i][j] = lambda_one;
|
||||
scale[i][j] = 1.0;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulLongSoft::init_style()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q");
|
||||
|
||||
neighbor->request(this);
|
||||
|
||||
cut_coulsq = cut_coul * cut_coul;
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulLongSoft::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) {
|
||||
if (lambda[i][i] != lambda[j][j])
|
||||
error->all(FLERR,"Pair coul/cut/soft different lambda values in mix");
|
||||
lambda[i][j] = lambda[i][i];
|
||||
}
|
||||
|
||||
lam1[i][j] = pow(lambda[i][j], nlambda);
|
||||
lam2[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
|
||||
|
||||
scale[j][i] = scale[i][j];
|
||||
lambda[j][i] = lambda[i][j];
|
||||
lam1[j][i] = lam1[i][j];
|
||||
lam2[j][i] = lam2[i][j];
|
||||
|
||||
return cut_coul+2.0*qdist;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulLongSoft::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(&lambda[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulLongSoft::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(&lambda[i][j],sizeof(double),1,fp);
|
||||
MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulLongSoft::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&nlambda,sizeof(double),1,fp);
|
||||
fwrite(&alphac,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 PairCoulLongSoft::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&nlambda,sizeof(double),1,fp);
|
||||
fread(&alphac,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(&nlambda,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alphac,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);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulLongSoft::single(int i, int j, int itype, int jtype,
|
||||
double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r,grij,expm2,t,erfc,prefactor;
|
||||
double forcecoul,phicoul;
|
||||
double denc;
|
||||
|
||||
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;
|
||||
|
||||
denc = sqrt(lam2[itype][jtype] + rsq);
|
||||
prefactor = force->qqrd2e * lam1[itype][jtype] * atom->q[i]*atom->q[j] /
|
||||
(denc*denc*denc);
|
||||
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
fforce = forcecoul;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
prefactor = force->qqrd2e * lam1[itype][jtype] * atom->q[i]*atom->q[j] / denc;
|
||||
phicoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
|
||||
} else phicoul = 0.0;
|
||||
|
||||
return phicoul;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairCoulLongSoft::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 0;
|
||||
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
|
||||
dim = 2;
|
||||
if (strcmp(str,"scale") == 0) return (void *) scale;
|
||||
if (strcmp(str,"lambda") == 0) return (void *) lambda;
|
||||
|
||||
return NULL;
|
||||
}
|
||||
|
|
@ -0,0 +1,80 @@
|
|||
/* -*- 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(coul/long/soft,PairCoulLongSoft)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_COUL_LONG_SOFT_H
|
||||
#define LMP_PAIR_COUL_LONG_SOFT_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairCoulLongSoft : public Pair {
|
||||
public:
|
||||
PairCoulLongSoft(class LAMMPS *);
|
||||
virtual ~PairCoulLongSoft();
|
||||
virtual void compute(int, int);
|
||||
virtual void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
virtual void init_style();
|
||||
double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
virtual void write_restart_settings(FILE *);
|
||||
virtual void read_restart_settings(FILE *);
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
virtual void *extract(const char *, int &);
|
||||
|
||||
protected:
|
||||
double cut_coul,cut_coulsq;
|
||||
double **scale;
|
||||
double **lambda;
|
||||
double nlambda, alphac;
|
||||
double **lam1, **lam2;
|
||||
double qdist; // TIP4P distance O to negative charge (compatibility of cutoffs)
|
||||
double g_ewald;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair style lj/cut/coul/long requires atom attribute q
|
||||
|
||||
The atom style defined does not have this attribute.
|
||||
|
||||
E: Pair style requires a KSpace style
|
||||
|
||||
No kspace style is defined.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,514 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_lj_cut_coul_cut_soft.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulCutSoft::PairLJCutCoulCutSoft(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
writedata = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulCutSoft::~PairLJCutCoulCutSoft()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(cut_lj);
|
||||
memory->destroy(cut_ljsq);
|
||||
memory->destroy(cut_coul);
|
||||
memory->destroy(cut_coulsq);
|
||||
memory->destroy(epsilon);
|
||||
memory->destroy(sigma);
|
||||
memory->destroy(lambda);
|
||||
memory->destroy(lj1);
|
||||
memory->destroy(lj2);
|
||||
memory->destroy(lj3);
|
||||
memory->destroy(lj4);
|
||||
memory->destroy(offset);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double rsq,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double denc, denlj, r4sig6;
|
||||
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;
|
||||
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;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_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]) {
|
||||
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
denc = sqrt(lj4[itype][jtype] + rsq);
|
||||
forcecoul = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
} else forcelj = 0.0;
|
||||
|
||||
fpair = factor_coul*forcecoul + factor_lj*forcelj;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq[itype][jtype])
|
||||
ecoul = factor_coul * qqrd2e * lj1[itype][jtype] * qtmp*q[j] / denc;
|
||||
else ecoul = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
evdwl,ecoul,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::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(cut_coul,n+1,n+1,"pair:cut_coul");
|
||||
memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq");
|
||||
memory->create(epsilon,n+1,n+1,"pair:epsilon");
|
||||
memory->create(sigma,n+1,n+1,"pair:sigma");
|
||||
memory->create(lambda,n+1,n+1,"pair:lambda");
|
||||
memory->create(lj1,n+1,n+1,"pair:lj1");
|
||||
memory->create(lj2,n+1,n+1,"pair:lj2");
|
||||
memory->create(lj3,n+1,n+1,"pair:lj3");
|
||||
memory->create(lj4,n+1,n+1,"pair:lj4");
|
||||
memory->create(offset,n+1,n+1,"pair:offset");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 4 || narg > 5) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
nlambda = force->numeric(FLERR,arg[0]);
|
||||
alphalj = force->numeric(FLERR,arg[1]);
|
||||
alphac = force->numeric(FLERR,arg[2]);
|
||||
|
||||
cut_lj_global = force->numeric(FLERR,arg[3]);
|
||||
if (narg == 4) cut_coul_global = cut_lj_global;
|
||||
else cut_coul_global = force->numeric(FLERR,arg[4]);
|
||||
|
||||
// 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;
|
||||
cut_coul[i][j] = cut_coul_global;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::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 epsilon_one = force->numeric(FLERR,arg[2]);
|
||||
double sigma_one = force->numeric(FLERR,arg[3]);
|
||||
double lambda_one = force->numeric(FLERR,arg[4]);
|
||||
|
||||
double cut_lj_one = cut_lj_global;
|
||||
double cut_coul_one = cut_coul_global;
|
||||
if (narg >= 6) cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[5]);
|
||||
if (narg == 7) cut_coul_one = force->numeric(FLERR,arg[6]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
epsilon[i][j] = epsilon_one;
|
||||
sigma[i][j] = sigma_one;
|
||||
lambda[i][j] = lambda_one;
|
||||
cut_lj[i][j] = cut_lj_one;
|
||||
cut_coul[i][j] = cut_coul_one;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::init_style()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Pair style lj/cut/coul/cut/soft requires atom attribute q");
|
||||
|
||||
neighbor->request(this);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulCutSoft::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]);
|
||||
if (lambda[i][i] != lambda[j][j])
|
||||
error->all(FLERR,"Pair lj/cut/coul/cut/soft different lambda values in mix");
|
||||
lambda[i][j] = lambda[i][i];
|
||||
cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
|
||||
cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]);
|
||||
}
|
||||
|
||||
double cut = MAX(cut_lj[i][j],cut_coul[i][j]);
|
||||
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
|
||||
cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j];
|
||||
|
||||
lj1[i][j] = pow(lambda[i][j], nlambda);
|
||||
lj2[i][j] = pow(sigma[i][j], 6.0);
|
||||
lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
|
||||
lj4[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
|
||||
|
||||
if (offset_flag) {
|
||||
double denlj = lj3[i][j] + pow(sigma[i][j] / cut_lj[i][j], 6.0);
|
||||
offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj);
|
||||
} else offset[i][j] = 0.0;
|
||||
|
||||
epsilon[j][i] = epsilon[i][j];
|
||||
sigma[j][i] = sigma[i][j];
|
||||
lambda[j][i] = lambda[i][j];
|
||||
cut_ljsq[j][i] = cut_ljsq[i][j];
|
||||
cut_coulsq[j][i] = cut_coulsq[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 sig2 = sigma[i][j]*sigma[i][j];
|
||||
double sig6 = sig2*sig2*sig2;
|
||||
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
|
||||
double rc6 = rc3*rc3;
|
||||
double rc9 = rc3*rc6;
|
||||
etail_ij = 8.0*MY_PI*all[0]*all[1]* lj1[i][j] * epsilon[i][j] *
|
||||
sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
|
||||
ptail_ij = 16.0*MY_PI*all[0]*all[1]* lj1[i][j] * epsilon[i][j] *
|
||||
sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
|
||||
}
|
||||
|
||||
return cut;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
fwrite(&setflag[i][j],sizeof(int),1,fp);
|
||||
if (setflag[i][j]) {
|
||||
fwrite(&epsilon[i][j],sizeof(double),1,fp);
|
||||
fwrite(&sigma[i][j],sizeof(double),1,fp);
|
||||
fwrite(&lambda[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_coul[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::read_restart(FILE *fp)
|
||||
{
|
||||
read_restart_settings(fp);
|
||||
allocate();
|
||||
|
||||
int i,j;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
|
||||
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
|
||||
if (setflag[i][j]) {
|
||||
if (me == 0) {
|
||||
fread(&epsilon[i][j],sizeof(double),1,fp);
|
||||
fread(&sigma[i][j],sizeof(double),1,fp);
|
||||
fread(&lambda[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_lj[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_coul[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&nlambda,sizeof(double),1,fp);
|
||||
fwrite(&alphalj,sizeof(double),1,fp);
|
||||
fwrite(&alphac,sizeof(double),1,fp);
|
||||
|
||||
fwrite(&cut_lj_global,sizeof(double),1,fp);
|
||||
fwrite(&cut_coul_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 PairLJCutCoulCutSoft::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&nlambda,sizeof(double),1,fp);
|
||||
fread(&alphalj,sizeof(double),1,fp);
|
||||
fread(&alphac,sizeof(double),1,fp);
|
||||
|
||||
fread(&cut_lj_global,sizeof(double),1,fp);
|
||||
fread(&cut_coul_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
fread(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alphalj,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alphac,1,MPI_DOUBLE,0,world);
|
||||
|
||||
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],lambda[i][i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes all pairs to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutSoft::write_data_all(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
for (int j = i; j <= atom->ntypes; j++)
|
||||
fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],
|
||||
lambda[i][j],cut_lj[i][j]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulCutSoft::single(int i, int j, int itype, int jtype,
|
||||
double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double forcecoul,forcelj,phicoul,philj;
|
||||
double denc, denlj, r4sig6;
|
||||
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
denc = sqrt(lj4[itype][jtype] + rsq);
|
||||
forcecoul = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] /
|
||||
(denc*denc*denc);
|
||||
} else forcecoul = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
} else forcelj = 0.0;
|
||||
fforce = factor_coul*forcecoul + factor_lj*forcelj;
|
||||
|
||||
double eng = 0.0;
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
phicoul = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] / denc;
|
||||
eng += factor_coul*phicoul;
|
||||
}
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
philj = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
eng += factor_lj*philj;
|
||||
}
|
||||
|
||||
return eng;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairLJCutCoulCutSoft::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 0;
|
||||
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
|
||||
dim = 2;
|
||||
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
|
||||
if (strcmp(str,"sigma") == 0) return (void *) sigma;
|
||||
if (strcmp(str,"lambda") == 0) return (void *) lambda;
|
||||
return NULL;
|
||||
}
|
||||
|
|
@ -0,0 +1,82 @@
|
|||
/* -*- 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/cut/coul/cut/soft,PairLJCutCoulCutSoft)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_LJ_CUT_COUL_CUT_SOFT_H
|
||||
#define LMP_PAIR_LJ_CUT_COUL_CUT_SOFT_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairLJCutCoulCutSoft : public Pair {
|
||||
public:
|
||||
PairLJCutCoulCutSoft(class LAMMPS *);
|
||||
virtual ~PairLJCutCoulCutSoft();
|
||||
virtual void compute(int, int);
|
||||
virtual void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
virtual void write_restart_settings(FILE *);
|
||||
virtual void read_restart_settings(FILE *);
|
||||
void write_data(FILE *);
|
||||
void write_data_all(FILE *);
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
void *extract(const char *, int &);
|
||||
|
||||
protected:
|
||||
double cut_lj_global,cut_coul_global;
|
||||
double **cut_lj,**cut_ljsq;
|
||||
double **cut_coul,**cut_coulsq;
|
||||
double **epsilon,**sigma, **lambda;
|
||||
double nlambda, alphalj, alphac;
|
||||
double **lj1,**lj2,**lj3,**lj4,**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: Pair style lj/cut/coul/cut/soft requires atom attribute q
|
||||
|
||||
The atom style defined does not have this attribute.
|
||||
|
||||
E: Pair lj/cut/coul/cut/soft different lambda values in mix
|
||||
|
||||
The value of lambda has to be the same for I J pairs.
|
||||
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,955 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Paul Crozier (SNL)
|
||||
Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_lj_cut_coul_long_soft.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulLongSoft::PairLJCutCoulLongSoft(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
ewaldflag = pppmflag = 1;
|
||||
respa_enable = 1;
|
||||
writedata = 1;
|
||||
qdist = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulLongSoft::~PairLJCutCoulLongSoft()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(cut_lj);
|
||||
memory->destroy(cut_ljsq);
|
||||
memory->destroy(epsilon);
|
||||
memory->destroy(sigma);
|
||||
memory->destroy(lambda);
|
||||
memory->destroy(lj1);
|
||||
memory->destroy(lj2);
|
||||
memory->destroy(lj3);
|
||||
memory->destroy(lj4);
|
||||
memory->destroy(offset);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,ii,j,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double r,rsq,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double denc, denlj, r4sig6;
|
||||
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;
|
||||
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;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_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]) {
|
||||
|
||||
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;
|
||||
|
||||
denc = sqrt(lj4[itype][jtype] + rsq);
|
||||
prefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
|
||||
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]) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
} else forcelj = 0.0;
|
||||
|
||||
fpair = forcecoul + factor_lj*forcelj;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq) {
|
||||
prefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / denc;
|
||||
ecoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
evdwl,ecoul,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::compute_inner()
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
|
||||
double rsq,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double rsw;
|
||||
double denc, denlj, r4sig6;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = listinner->inum;
|
||||
ilist = listinner->ilist;
|
||||
numneigh = listinner->numneigh;
|
||||
firstneigh = listinner->firstneigh;
|
||||
|
||||
double cut_out_on = cut_respa[0];
|
||||
double cut_out_off = cut_respa[1];
|
||||
|
||||
double cut_out_diff = cut_out_off - cut_out_on;
|
||||
double cut_out_on_sq = cut_out_on*cut_out_on;
|
||||
double cut_out_off_sq = cut_out_off*cut_out_off;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_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;
|
||||
|
||||
if (rsq < cut_out_off_sq) {
|
||||
jtype = type[j];
|
||||
|
||||
denc = sqrt(lj4[itype][jtype] + rsq);
|
||||
forcecoul = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
} else forcelj = 0.0;
|
||||
|
||||
fpair = forcecoul + factor_lj*forcelj;
|
||||
|
||||
if (rsq > cut_out_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
|
||||
fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::compute_middle()
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
|
||||
double rsq,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double rsw;
|
||||
double denc, denlj, r4sig6;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = listmiddle->inum;
|
||||
ilist = listmiddle->ilist;
|
||||
numneigh = listmiddle->numneigh;
|
||||
firstneigh = listmiddle->firstneigh;
|
||||
|
||||
double cut_in_off = cut_respa[0];
|
||||
double cut_in_on = cut_respa[1];
|
||||
double cut_out_on = cut_respa[2];
|
||||
double cut_out_off = cut_respa[3];
|
||||
|
||||
double cut_in_diff = cut_in_on - cut_in_off;
|
||||
double cut_out_diff = cut_out_off - cut_out_on;
|
||||
double cut_in_off_sq = cut_in_off*cut_in_off;
|
||||
double cut_in_on_sq = cut_in_on*cut_in_on;
|
||||
double cut_out_on_sq = cut_out_on*cut_out_on;
|
||||
double cut_out_off_sq = cut_out_off*cut_out_off;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_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;
|
||||
|
||||
if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
|
||||
jtype = type[j];
|
||||
|
||||
denc = sqrt(lj4[itype][jtype] + rsq);
|
||||
forcecoul = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
} else forcelj = 0.0;
|
||||
|
||||
fpair = forcecoul + factor_lj*forcelj;
|
||||
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
|
||||
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
|
||||
}
|
||||
if (rsq > cut_out_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
|
||||
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::compute_outer(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double r,rsq,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,fprefactor,eprefactor,t,erfc;
|
||||
double rsw;
|
||||
double denc, denlj, r4sig6;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = listouter->inum;
|
||||
ilist = listouter->ilist;
|
||||
numneigh = listouter->numneigh;
|
||||
firstneigh = listouter->firstneigh;
|
||||
|
||||
double cut_in_off = cut_respa[2];
|
||||
double cut_in_on = cut_respa[3];
|
||||
|
||||
double cut_in_diff = cut_in_on - cut_in_off;
|
||||
double cut_in_off_sq = cut_in_off*cut_in_off;
|
||||
double cut_in_on_sq = cut_in_on*cut_in_on;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_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]) {
|
||||
|
||||
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;
|
||||
|
||||
denc = sqrt(lj4[itype][jtype] + rsq);
|
||||
fprefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
|
||||
forcecoul = fprefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
|
||||
|
||||
if (rsq > cut_in_off_sq) {
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (r - cut_in_off)/cut_in_diff;
|
||||
forcecoul += fprefactor*rsw*rsw*(3.0 - 2.0*rsw);
|
||||
if (factor_coul < 1.0)
|
||||
forcecoul -=
|
||||
(1.0-factor_coul)*fprefactor*rsw*rsw*(3.0 - 2.0*rsw);
|
||||
} else {
|
||||
forcecoul += fprefactor;
|
||||
if (factor_coul < 1.0)
|
||||
forcecoul -= (1.0-factor_coul)*fprefactor;
|
||||
}
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
|
||||
forcelj *= rsw*rsw*(3.0 - 2.0*rsw);
|
||||
}
|
||||
} else forcelj = 0.0;
|
||||
|
||||
fpair = forcecoul + forcelj;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq) {
|
||||
eprefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / denc;
|
||||
ecoul = eprefactor*erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*eprefactor;
|
||||
} else ecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (vflag) {
|
||||
if (rsq < cut_coulsq) {
|
||||
forcecoul = fprefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*fprefactor;
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
if (rsq <= cut_in_off_sq) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
} else if (rsq < cut_in_on_sq) {
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
}
|
||||
fpair = forcecoul + factor_lj*forcelj;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
evdwl,ecoul,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::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(lambda,n+1,n+1,"pair:lambda");
|
||||
memory->create(lj1,n+1,n+1,"pair:lj1");
|
||||
memory->create(lj2,n+1,n+1,"pair:lj2");
|
||||
memory->create(lj3,n+1,n+1,"pair:lj3");
|
||||
memory->create(lj4,n+1,n+1,"pair:lj4");
|
||||
memory->create(offset,n+1,n+1,"pair:offset");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 4 || narg > 5) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
nlambda = force->numeric(FLERR,arg[0]);
|
||||
alphalj = force->numeric(FLERR,arg[1]);
|
||||
alphac = force->numeric(FLERR,arg[2]);
|
||||
|
||||
cut_lj_global = force->numeric(FLERR,arg[3]);
|
||||
if (narg == 4) cut_coul = cut_lj_global;
|
||||
else cut_coul = force->numeric(FLERR,arg[4]);
|
||||
|
||||
// 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 PairLJCutCoulLongSoft::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 5 || narg > 6)
|
||||
error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(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 lambda_one = force->numeric(FLERR,arg[4]);
|
||||
|
||||
double cut_lj_one = cut_lj_global;
|
||||
if (narg == 6) cut_lj_one = force->numeric(FLERR,arg[5]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
epsilon[i][j] = epsilon_one;
|
||||
sigma[i][j] = sigma_one;
|
||||
lambda[i][j] = lambda_one;
|
||||
cut_lj[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 PairLJCutCoulLongSoft::init_style()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Pair style lj/cut/coul/long/soft requires atom attribute q");
|
||||
|
||||
// request regular or rRESPA neighbor lists
|
||||
|
||||
int irequest;
|
||||
|
||||
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
|
||||
int respa = 0;
|
||||
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
|
||||
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
|
||||
|
||||
if (respa == 0) irequest = neighbor->request(this);
|
||||
else if (respa == 1) {
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 1;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respainner = 1;
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 3;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respaouter = 1;
|
||||
} else {
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 1;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respainner = 1;
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 2;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respamiddle = 1;
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 3;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respaouter = 1;
|
||||
}
|
||||
|
||||
} else irequest = neighbor->request(this);
|
||||
|
||||
cut_coulsq = cut_coul * cut_coul;
|
||||
|
||||
// set rRESPA cutoffs
|
||||
|
||||
if (strstr(update->integrate_style,"respa") &&
|
||||
((Respa *) update->integrate)->level_inner >= 0)
|
||||
cut_respa = ((Respa *) update->integrate)->cutoff;
|
||||
else cut_respa = NULL;
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
neighbor callback to inform pair style of neighbor list to use
|
||||
regular or rRESPA
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::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;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulLongSoft::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]);
|
||||
if (lambda[i][i] != lambda[j][j])
|
||||
error->all(FLERR,"Pair lj/cut/coul/long/soft different lambda values in mix");
|
||||
lambda[i][j] = lambda[i][i];
|
||||
cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
|
||||
}
|
||||
|
||||
// include TIP4P qdist in full cutoff, qdist = 0.0 if not TIP4P
|
||||
|
||||
double cut = MAX(cut_lj[i][j],cut_coul+2.0*qdist);
|
||||
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
|
||||
|
||||
lj1[i][j] = pow(lambda[i][j], nlambda);
|
||||
lj2[i][j] = pow(sigma[i][j], 6.0);
|
||||
lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
|
||||
lj4[i][j] = alphac * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
|
||||
|
||||
if (offset_flag) {
|
||||
double denlj = lj3[i][j] + pow(sigma[i][j] / cut_lj[i][j], 6.0);
|
||||
offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj);
|
||||
} else offset[i][j] = 0.0;
|
||||
|
||||
epsilon[j][i] = epsilon[i][j];
|
||||
sigma[j][i] = sigma[i][j];
|
||||
lambda[j][i] = lambda[i][j];
|
||||
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];
|
||||
|
||||
// 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");
|
||||
|
||||
// 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 sig2 = sigma[i][j]*sigma[i][j];
|
||||
double sig6 = sig2*sig2*sig2;
|
||||
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
|
||||
double rc6 = rc3*rc3;
|
||||
double rc9 = rc3*rc6;
|
||||
etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
|
||||
sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
|
||||
ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
|
||||
sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
|
||||
}
|
||||
|
||||
return cut;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
fwrite(&setflag[i][j],sizeof(int),1,fp);
|
||||
if (setflag[i][j]) {
|
||||
fwrite(&epsilon[i][j],sizeof(double),1,fp);
|
||||
fwrite(&sigma[i][j],sizeof(double),1,fp);
|
||||
fwrite(&lambda[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::read_restart(FILE *fp)
|
||||
{
|
||||
read_restart_settings(fp);
|
||||
|
||||
allocate();
|
||||
|
||||
int i,j;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
|
||||
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
|
||||
if (setflag[i][j]) {
|
||||
if (me == 0) {
|
||||
fread(&epsilon[i][j],sizeof(double),1,fp);
|
||||
fread(&sigma[i][j],sizeof(double),1,fp);
|
||||
fread(&lambda[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_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(&lambda[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 PairLJCutCoulLongSoft::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&nlambda,sizeof(double),1,fp);
|
||||
fwrite(&alphalj,sizeof(double),1,fp);
|
||||
fwrite(&alphac,sizeof(double),1,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(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&nlambda,sizeof(double),1,fp);
|
||||
fread(&alphalj,sizeof(double),1,fp);
|
||||
fread(&alphac,sizeof(double),1,fp);
|
||||
|
||||
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(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alphalj,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alphac,1,MPI_DOUBLE,0,world);
|
||||
|
||||
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(&tail_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],lambda[i][i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes all pairs to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongSoft::write_data_all(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
for (int j = i; j <= atom->ntypes; j++)
|
||||
fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],
|
||||
lambda[i][j],cut_lj[i][j]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulLongSoft::single(int i, int j, int itype, int jtype,
|
||||
double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r,grij,expm2,t,erfc,prefactor;
|
||||
double forcecoul,forcelj,phicoul,philj;
|
||||
double denc, denlj, r4sig6;
|
||||
|
||||
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;
|
||||
|
||||
denc = sqrt(lj4[itype][jtype] + rsq);
|
||||
prefactor = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] /
|
||||
(denc*denc*denc);
|
||||
|
||||
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]) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
} else forcelj = 0.0;
|
||||
|
||||
fforce = forcecoul + factor_lj*forcelj;
|
||||
|
||||
double eng = 0.0;
|
||||
if (rsq < cut_coulsq) {
|
||||
prefactor = force->qqrd2e * lj1[itype][jtype] * atom->q[i]*atom->q[j] / denc;
|
||||
phicoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
|
||||
eng += phicoul;
|
||||
}
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
philj = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
eng += factor_lj*philj;
|
||||
}
|
||||
|
||||
return eng;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairLJCutCoulLongSoft::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 0;
|
||||
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
|
||||
dim = 2;
|
||||
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
|
||||
if (strcmp(str,"sigma") == 0) return (void *) sigma;
|
||||
if (strcmp(str,"lambda") == 0) return (void *) lambda;
|
||||
return NULL;
|
||||
}
|
||||
|
|
@ -0,0 +1,94 @@
|
|||
/* -*- 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/cut/coul/long/soft,PairLJCutCoulLongSoft)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_LJ_CUT_COUL_LONG_SOFT_H
|
||||
#define LMP_PAIR_LJ_CUT_COUL_LONG_SOFT_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairLJCutCoulLongSoft : public Pair {
|
||||
public:
|
||||
PairLJCutCoulLongSoft(class LAMMPS *);
|
||||
virtual ~PairLJCutCoulLongSoft();
|
||||
virtual void compute(int, int);
|
||||
virtual void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
virtual void init_style();
|
||||
void init_list(int, class NeighList *);
|
||||
virtual double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
virtual void write_restart_settings(FILE *);
|
||||
virtual void read_restart_settings(FILE *);
|
||||
void write_data(FILE *);
|
||||
void write_data_all(FILE *);
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
virtual void *extract(const char *, int &);
|
||||
|
||||
void compute_inner();
|
||||
void compute_middle();
|
||||
virtual void compute_outer(int, int);
|
||||
|
||||
protected:
|
||||
double cut_lj_global;
|
||||
double **cut_lj,**cut_ljsq;
|
||||
double cut_coul,cut_coulsq;
|
||||
double **epsilon,**sigma, **lambda;
|
||||
double nlambda, alphalj, alphac;
|
||||
double **lj1,**lj2,**lj3,**lj4,**offset;
|
||||
double *cut_respa;
|
||||
double qdist; // TIP4P distance O to negative charge (compatibility of cutoffs)
|
||||
double g_ewald;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair style lj/cut/coul/long requires atom attribute q
|
||||
|
||||
The atom style defined does not have this attribute.
|
||||
|
||||
E: Pair style requires a KSpace style
|
||||
|
||||
No kspace style is defined.
|
||||
|
||||
E: Pair cutoff < Respa interior cutoff
|
||||
|
||||
One or more pairwise cutoffs are too short to use with the specified
|
||||
rRESPA cutoffs.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,779 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Paul Crozier (SNL)
|
||||
Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_lj_cut_soft.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutSoft::PairLJCutSoft(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
respa_enable = 1;
|
||||
writedata = 1;
|
||||
allocated = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutSoft::~PairLJCutSoft()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(cut);
|
||||
memory->destroy(epsilon);
|
||||
memory->destroy(sigma);
|
||||
memory->destroy(lambda);
|
||||
memory->destroy(lj1);
|
||||
memory->destroy(lj2);
|
||||
memory->destroy(lj3);
|
||||
memory->destroy(offset);
|
||||
allocated=0;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
|
||||
double rsq,forcelj,factor_lj;
|
||||
double denlj, r4sig6;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
|
||||
fpair = factor_lj*forcelj;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
evdwl,0.0,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::compute_inner()
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
|
||||
double rsq,forcelj,factor_lj,rsw;
|
||||
double denlj, r4sig6;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
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 = listinner->inum;
|
||||
ilist = listinner->ilist;
|
||||
numneigh = listinner->numneigh;
|
||||
firstneigh = listinner->firstneigh;
|
||||
|
||||
double cut_out_on = cut_respa[0];
|
||||
double cut_out_off = cut_respa[1];
|
||||
|
||||
double cut_out_diff = cut_out_off - cut_out_on;
|
||||
double cut_out_on_sq = cut_out_on*cut_out_on;
|
||||
double cut_out_off_sq = cut_out_off*cut_out_off;
|
||||
|
||||
// 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;
|
||||
|
||||
if (rsq < cut_out_off_sq) {
|
||||
jtype = type[j];
|
||||
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
|
||||
fpair = factor_lj*forcelj;
|
||||
|
||||
if (rsq > cut_out_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
|
||||
fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw);
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::compute_middle()
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
|
||||
double rsq,forcelj,factor_lj,rsw;
|
||||
double denlj, r4sig6;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
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 = listmiddle->inum;
|
||||
ilist = listmiddle->ilist;
|
||||
numneigh = listmiddle->numneigh;
|
||||
firstneigh = listmiddle->firstneigh;
|
||||
|
||||
double cut_in_off = cut_respa[0];
|
||||
double cut_in_on = cut_respa[1];
|
||||
double cut_out_on = cut_respa[2];
|
||||
double cut_out_off = cut_respa[3];
|
||||
|
||||
double cut_in_diff = cut_in_on - cut_in_off;
|
||||
double cut_out_diff = cut_out_off - cut_out_on;
|
||||
double cut_in_off_sq = cut_in_off*cut_in_off;
|
||||
double cut_in_on_sq = cut_in_on*cut_in_on;
|
||||
double cut_out_on_sq = cut_out_on*cut_out_on;
|
||||
double cut_out_off_sq = cut_out_off*cut_out_off;
|
||||
|
||||
// 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;
|
||||
|
||||
if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
|
||||
jtype = type[j];
|
||||
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
|
||||
fpair = factor_lj*forcelj;
|
||||
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
|
||||
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
|
||||
}
|
||||
if (rsq > cut_out_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
|
||||
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::compute_outer(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
|
||||
double rsq,forcelj,factor_lj,rsw;
|
||||
double denlj, r4sig6;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = 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 = listouter->inum;
|
||||
ilist = listouter->ilist;
|
||||
numneigh = listouter->numneigh;
|
||||
firstneigh = listouter->firstneigh;
|
||||
|
||||
double cut_in_off = cut_respa[2];
|
||||
double cut_in_on = cut_respa[3];
|
||||
|
||||
double cut_in_diff = cut_in_on - cut_in_off;
|
||||
double cut_in_off_sq = cut_in_off*cut_in_off;
|
||||
double cut_in_on_sq = cut_in_on*cut_in_on;
|
||||
|
||||
// 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]) {
|
||||
if (rsq > cut_in_off_sq) {
|
||||
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
|
||||
fpair = factor_lj*forcelj;
|
||||
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
|
||||
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
|
||||
}
|
||||
|
||||
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 = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
}
|
||||
|
||||
if (vflag) {
|
||||
if (rsq <= cut_in_off_sq) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
fpair = factor_lj*forcelj;
|
||||
} else if (rsq < cut_in_on_sq)
|
||||
fpair = factor_lj*forcelj;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
evdwl,0.0,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes;
|
||||
|
||||
memory->create(setflag,n+1,n+1,"pair:setflag");
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = i; j <= n; j++)
|
||||
setflag[i][j] = 0;
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
memory->create(cut,n+1,n+1,"pair:cut");
|
||||
memory->create(epsilon,n+1,n+1,"pair:epsilon");
|
||||
memory->create(sigma,n+1,n+1,"pair:sigma");
|
||||
memory->create(lambda,n+1,n+1,"pair:lambda");
|
||||
memory->create(lj1,n+1,n+1,"pair:lj1");
|
||||
memory->create(lj2,n+1,n+1,"pair:lj2");
|
||||
memory->create(lj3,n+1,n+1,"pair:lj3");
|
||||
memory->create(offset,n+1,n+1,"pair:offset");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
nlambda = force->numeric(FLERR,arg[0]);
|
||||
alphalj = force->numeric(FLERR,arg[1]);
|
||||
|
||||
cut_global = force->numeric(FLERR,arg[2]);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i+1; j <= atom->ntypes; j++)
|
||||
if (setflag[i][j]) cut[i][j] = cut_global;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 5 || narg > 6)
|
||||
error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(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 lambda_one = force->numeric(FLERR,arg[4]);
|
||||
|
||||
double cut_one = cut_global;
|
||||
if (narg == 6) cut_one = force->numeric(FLERR,arg[5]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
epsilon[i][j] = epsilon_one;
|
||||
sigma[i][j] = sigma_one;
|
||||
lambda[i][j] = lambda_one;
|
||||
cut[i][j] = cut_one;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::init_style()
|
||||
{
|
||||
// request regular or rRESPA neighbor lists
|
||||
|
||||
int irequest;
|
||||
|
||||
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
|
||||
int respa = 0;
|
||||
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
|
||||
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
|
||||
|
||||
if (respa == 0) irequest = neighbor->request(this);
|
||||
else if (respa == 1) {
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 1;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respainner = 1;
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 3;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respaouter = 1;
|
||||
} else {
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 1;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respainner = 1;
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 2;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respamiddle = 1;
|
||||
irequest = neighbor->request(this);
|
||||
neighbor->requests[irequest]->id = 3;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->respaouter = 1;
|
||||
}
|
||||
|
||||
} else irequest = neighbor->request(this);
|
||||
|
||||
// set rRESPA cutoffs
|
||||
|
||||
if (strstr(update->integrate_style,"respa") &&
|
||||
((Respa *) update->integrate)->level_inner >= 0)
|
||||
cut_respa = ((Respa *) update->integrate)->cutoff;
|
||||
else cut_respa = NULL;
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
neighbor callback to inform pair style of neighbor list to use
|
||||
regular or rRESPA
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::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;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutSoft::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[i][j] = mix_distance(cut[i][i],cut[j][j]);
|
||||
}
|
||||
|
||||
lj1[i][j] = pow(lambda[i][j], nlambda);
|
||||
lj2[i][j] = pow(sigma[i][j], 6.0);
|
||||
lj3[i][j] = alphalj * (1.0 - lambda[i][j])*(1.0 - lambda[i][j]);
|
||||
|
||||
if (offset_flag) {
|
||||
double denlj = lj3[i][j] + pow(sigma[i][j] / cut[i][j], 6.0);
|
||||
offset[i][j] = lj1[i][j] * 4.0 * epsilon[i][j] * (1.0/(denlj*denlj) - 1.0/denlj);
|
||||
} else offset[i][j] = 0.0;
|
||||
|
||||
epsilon[j][i] = epsilon[i][j];
|
||||
sigma[j][i] = sigma[i][j];
|
||||
lambda[j][i] = lambda[i][j];
|
||||
lj1[j][i] = lj1[i][j];
|
||||
lj2[j][i] = lj2[i][j];
|
||||
lj3[j][i] = lj3[i][j];
|
||||
offset[j][i] = offset[i][j];
|
||||
|
||||
// check interior rRESPA cutoff
|
||||
|
||||
if (cut_respa && cut[i][j] < cut_respa[3])
|
||||
error->all(FLERR,"Pair cutoff < Respa interior cutoff");
|
||||
|
||||
// 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 sig2 = sigma[i][j]*sigma[i][j];
|
||||
double sig6 = sig2*sig2*sig2;
|
||||
double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
|
||||
double rc6 = rc3*rc3;
|
||||
double rc9 = rc3*rc6;
|
||||
etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
|
||||
sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
|
||||
ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
|
||||
sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
|
||||
}
|
||||
|
||||
return cut[i][j];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
fwrite(&setflag[i][j],sizeof(int),1,fp);
|
||||
if (setflag[i][j]) {
|
||||
fwrite(&epsilon[i][j],sizeof(double),1,fp);
|
||||
fwrite(&sigma[i][j],sizeof(double),1,fp);
|
||||
fwrite(&lambda[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::read_restart(FILE *fp)
|
||||
{
|
||||
read_restart_settings(fp);
|
||||
allocate();
|
||||
|
||||
int i,j;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
|
||||
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
|
||||
if (setflag[i][j]) {
|
||||
if (me == 0) {
|
||||
fread(&epsilon[i][j],sizeof(double),1,fp);
|
||||
fread(&sigma[i][j],sizeof(double),1,fp);
|
||||
fread(&lambda[i][j],sizeof(double),1,fp);
|
||||
fread(&cut[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&lambda[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&nlambda,sizeof(double),1,fp);
|
||||
fwrite(&alphalj,sizeof(double),1,fp);
|
||||
|
||||
fwrite(&cut_global,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
fwrite(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::read_restart_settings(FILE *fp)
|
||||
{
|
||||
int me = comm->me;
|
||||
if (me == 0) {
|
||||
fread(&nlambda,sizeof(double),1,fp);
|
||||
fread(&alphalj,sizeof(double),1,fp);
|
||||
|
||||
fread(&cut_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
fread(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&nlambda,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alphalj,1,MPI_DOUBLE,0,world);
|
||||
|
||||
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],lambda[i][i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes all pairs to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutSoft::write_data_all(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
for (int j = i; j <= atom->ntypes; j++)
|
||||
fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],
|
||||
lambda[i][j],cut[i][j]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutSoft::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double forcelj,philj;
|
||||
double r4sig6, denlj;
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
} else forcelj = 0.0;
|
||||
fforce = factor_lj*forcelj;
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
philj = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
} else philj = 0.0;
|
||||
|
||||
return factor_lj*philj;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairLJCutSoft::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 2;
|
||||
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
|
||||
if (strcmp(str,"sigma") == 0) return (void *) sigma;
|
||||
if (strcmp(str,"lambda") == 0) return (void *) lambda;
|
||||
return NULL;
|
||||
}
|
||||
|
|
@ -0,0 +1,87 @@
|
|||
/* -*- 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/cut/soft,PairLJCutSoft)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_LJ_CUT_SOFT_H
|
||||
#define LMP_PAIR_LJ_CUT_SOFT_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairLJCutSoft : public Pair {
|
||||
public:
|
||||
PairLJCutSoft(class LAMMPS *);
|
||||
virtual ~PairLJCutSoft();
|
||||
virtual void compute(int, int);
|
||||
virtual void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
virtual void init_style();
|
||||
void init_list(int, class NeighList *);
|
||||
virtual double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
virtual void write_restart_settings(FILE *);
|
||||
virtual void read_restart_settings(FILE *);
|
||||
void write_data(FILE *);
|
||||
void write_data_all(FILE *);
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
void *extract(const char *, int &);
|
||||
|
||||
void compute_inner();
|
||||
void compute_middle();
|
||||
void compute_outer(int, int);
|
||||
|
||||
protected:
|
||||
double cut_global;
|
||||
double **cut;
|
||||
double **epsilon,**sigma, **lambda;
|
||||
double nlambda, alphalj;
|
||||
double **lj1,**lj2,**lj3,**offset;
|
||||
double *cut_respa;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair cutoff < Respa interior cutoff
|
||||
|
||||
One or more pairwise cutoffs are too short to use with the specified
|
||||
rRESPA cutoffs.
|
||||
|
||||
E: Pair lj/cut/soft different lambda values in mix
|
||||
|
||||
The value of lambda has to be the same for I J pairs.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,592 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 authors: Amalie Frischknecht and Ahmed Ismail (SNL)
|
||||
simpler force assignment added by Rolf Isele-Holder (Aachen University)
|
||||
Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_lj_cut_tip4p_long_soft.h"
|
||||
#include "angle.h"
|
||||
#include "atom.h"
|
||||
#include "bond.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "update.h"
|
||||
#include "respa.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutTIP4PLongSoft::PairLJCutTIP4PLongSoft(LAMMPS *lmp) :
|
||||
PairLJCutCoulLongSoft(lmp)
|
||||
{
|
||||
tip4pflag = 1;
|
||||
single_enable = 0;
|
||||
respa_enable = 0;
|
||||
|
||||
nmax = 0;
|
||||
hneigh = NULL;
|
||||
newsite = NULL;
|
||||
|
||||
// TIP4P cannot compute virial as F dot r
|
||||
// due to finding bonded H atoms which are not near O atom
|
||||
|
||||
no_virial_fdotr_compute = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutTIP4PLongSoft::~PairLJCutTIP4PLongSoft()
|
||||
{
|
||||
memory->destroy(hneigh);
|
||||
memory->destroy(newsite);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype,key;
|
||||
int n,vlist[6];
|
||||
int iH1,iH2,jH1,jH2;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
|
||||
double r,forcecoul,forcelj,cforce;
|
||||
double factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double denc, denlj, r4sig6;
|
||||
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
|
||||
double *x1,*x2;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
// reallocate hneigh & newsite if necessary
|
||||
// initialize hneigh[0] to -1 on steps when reneighboring occurred
|
||||
// initialize hneigh[2] to 0 every step
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
int nall = nlocal + atom->nghost;
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
nmax = atom->nmax;
|
||||
memory->destroy(hneigh);
|
||||
memory->create(hneigh,nmax,3,"pair:hneigh");
|
||||
memory->destroy(newsite);
|
||||
memory->create(newsite,nmax,3,"pair:newsite");
|
||||
}
|
||||
if (neighbor->ago == 0)
|
||||
for (i = 0; i < nall; i++) hneigh[i][0] = -1;
|
||||
for (i = 0; i < nall; i++) hneigh[i][2] = 0;
|
||||
|
||||
double **f = atom->f;
|
||||
double **x = atom->x;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
double *special_coul = force->special_coul;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
double cut_coulsqplus = (cut_coul+2.0*qdist) * (cut_coul+2.0*qdist);
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
|
||||
// if atom I = water O, set x1 = offset charge site
|
||||
// else x1 = x of atom I
|
||||
|
||||
if (itype == typeO) {
|
||||
if (hneigh[i][0] < 0) {
|
||||
hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1);
|
||||
hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2);
|
||||
hneigh[i][2] = 1;
|
||||
if (iH1 == -1 || iH2 == -1)
|
||||
error->one(FLERR,"TIP4P hydrogen is missing");
|
||||
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
|
||||
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
|
||||
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
|
||||
} else {
|
||||
iH1 = hneigh[i][0];
|
||||
iH2 = hneigh[i][1];
|
||||
if (hneigh[i][2] == 0) {
|
||||
hneigh[i][2] = 1;
|
||||
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
|
||||
}
|
||||
}
|
||||
x1 = newsite[i];
|
||||
} else x1 = x[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];
|
||||
|
||||
// LJ interaction based on true rsq
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
|
||||
r4sig6 = rsq*rsq / lj2[itype][jtype];
|
||||
denlj = lj3[itype][jtype] + rsq*r4sig6;
|
||||
forcelj = lj1[itype][jtype] * epsilon[itype][jtype] *
|
||||
(48.0*r4sig6/(denlj*denlj*denlj) - 24.0*r4sig6/(denlj*denlj));
|
||||
|
||||
forcelj *= factor_lj;
|
||||
|
||||
f[i][0] += delx*forcelj;
|
||||
f[i][1] += dely*forcelj;
|
||||
f[i][2] += delz*forcelj;
|
||||
f[j][0] -= delx*forcelj;
|
||||
f[j][1] -= dely*forcelj;
|
||||
f[j][2] -= delz*forcelj;
|
||||
|
||||
if (eflag) {
|
||||
evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
|
||||
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
evdwl,0.0,forcelj,delx,dely,delz);
|
||||
}
|
||||
|
||||
// adjust rsq and delxyz for off-site O charge(s) if necessary
|
||||
// but only if they are within reach
|
||||
|
||||
if (rsq < cut_coulsqplus) {
|
||||
if (itype == typeO || jtype == typeO) {
|
||||
|
||||
// if atom J = water O, set x2 = offset charge site
|
||||
// else x2 = x of atom J
|
||||
|
||||
if (jtype == typeO) {
|
||||
if (hneigh[j][0] < 0) {
|
||||
hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
|
||||
hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);
|
||||
hneigh[j][2] = 1;
|
||||
if (jH1 == -1 || jH2 == -1)
|
||||
error->one(FLERR,"TIP4P hydrogen is missing");
|
||||
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
|
||||
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
|
||||
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
|
||||
} else {
|
||||
jH1 = hneigh[j][0];
|
||||
jH2 = hneigh[j][1];
|
||||
if (hneigh[j][2] == 0) {
|
||||
hneigh[j][2] = 1;
|
||||
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
|
||||
}
|
||||
}
|
||||
x2 = newsite[j];
|
||||
} else x2 = x[j];
|
||||
|
||||
delx = x1[0] - x2[0];
|
||||
dely = x1[1] - x2[1];
|
||||
delz = x1[2] - x2[2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
}
|
||||
|
||||
// Coulombic interaction based on modified 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;
|
||||
|
||||
denc = sqrt(lj4[itype][jtype] + rsq);
|
||||
prefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
|
||||
cforce = forcecoul;
|
||||
|
||||
// if i,j are not O atoms, force is applied directly
|
||||
// if i or j are O atoms, force is on fictitious atom & partitioned
|
||||
// force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999)
|
||||
// f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f
|
||||
// preserves total force and torque on water molecule
|
||||
// virial = sum(r x F) where each water's atoms are near xi and xj
|
||||
// vlist stores 2,4,6 atoms whose forces contribute to virial
|
||||
|
||||
n = 0;
|
||||
key = 0;
|
||||
|
||||
if (itype != typeO) {
|
||||
f[i][0] += delx * cforce;
|
||||
f[i][1] += dely * cforce;
|
||||
f[i][2] += delz * cforce;
|
||||
|
||||
if (vflag) {
|
||||
v[0] = x[i][0] * delx * cforce;
|
||||
v[1] = x[i][1] * dely * cforce;
|
||||
v[2] = x[i][2] * delz * cforce;
|
||||
v[3] = x[i][0] * dely * cforce;
|
||||
v[4] = x[i][0] * delz * cforce;
|
||||
v[5] = x[i][1] * delz * cforce;
|
||||
}
|
||||
vlist[n++] = i;
|
||||
|
||||
} else {
|
||||
key++;
|
||||
|
||||
fd[0] = delx*cforce;
|
||||
fd[1] = dely*cforce;
|
||||
fd[2] = delz*cforce;
|
||||
|
||||
fO[0] = fd[0]*(1 - alpha);
|
||||
fO[1] = fd[1]*(1 - alpha);
|
||||
fO[2] = fd[2]*(1 - alpha);
|
||||
|
||||
fH[0] = 0.5 * alpha * fd[0];
|
||||
fH[1] = 0.5 * alpha * fd[1];
|
||||
fH[2] = 0.5 * alpha * fd[2];
|
||||
|
||||
f[i][0] += fO[0];
|
||||
f[i][1] += fO[1];
|
||||
f[i][2] += fO[2];
|
||||
|
||||
f[iH1][0] += fH[0];
|
||||
f[iH1][1] += fH[1];
|
||||
f[iH1][2] += fH[2];
|
||||
|
||||
f[iH2][0] += fH[0];
|
||||
f[iH2][1] += fH[1];
|
||||
f[iH2][2] += fH[2];
|
||||
|
||||
if (vflag) {
|
||||
domain->closest_image(x[i],x[iH1],xH1);
|
||||
domain->closest_image(x[i],x[iH2],xH2);
|
||||
|
||||
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
|
||||
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
|
||||
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
|
||||
v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
|
||||
v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
|
||||
v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
|
||||
}
|
||||
vlist[n++] = i;
|
||||
vlist[n++] = iH1;
|
||||
vlist[n++] = iH2;
|
||||
}
|
||||
|
||||
if (jtype != typeO) {
|
||||
f[j][0] -= delx * cforce;
|
||||
f[j][1] -= dely * cforce;
|
||||
f[j][2] -= delz * cforce;
|
||||
|
||||
if (vflag) {
|
||||
v[0] -= x[j][0] * delx * cforce;
|
||||
v[1] -= x[j][1] * dely * cforce;
|
||||
v[2] -= x[j][2] * delz * cforce;
|
||||
v[3] -= x[j][0] * dely * cforce;
|
||||
v[4] -= x[j][0] * delz * cforce;
|
||||
v[5] -= x[j][1] * delz * cforce;
|
||||
}
|
||||
vlist[n++] = j;
|
||||
|
||||
} else {
|
||||
key += 2;
|
||||
|
||||
fd[0] = -delx*cforce;
|
||||
fd[1] = -dely*cforce;
|
||||
fd[2] = -delz*cforce;
|
||||
|
||||
fO[0] = fd[0]*(1 - alpha);
|
||||
fO[1] = fd[1]*(1 - alpha);
|
||||
fO[2] = fd[2]*(1 - alpha);
|
||||
|
||||
fH[0] = 0.5 * alpha * fd[0];
|
||||
fH[1] = 0.5 * alpha * fd[1];
|
||||
fH[2] = 0.5 * alpha * fd[2];
|
||||
|
||||
f[j][0] += fO[0];
|
||||
f[j][1] += fO[1];
|
||||
f[j][2] += fO[2];
|
||||
|
||||
f[jH1][0] += fH[0];
|
||||
f[jH1][1] += fH[1];
|
||||
f[jH1][2] += fH[2];
|
||||
|
||||
f[jH2][0] += fH[0];
|
||||
f[jH2][1] += fH[1];
|
||||
f[jH2][2] += fH[2];
|
||||
|
||||
if (vflag) {
|
||||
domain->closest_image(x[j],x[jH1],xH1);
|
||||
domain->closest_image(x[j],x[jH2],xH2);
|
||||
|
||||
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
|
||||
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
|
||||
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
|
||||
v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
|
||||
v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
|
||||
v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
|
||||
}
|
||||
vlist[n++] = j;
|
||||
vlist[n++] = jH1;
|
||||
vlist[n++] = jH2;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
prefactor = qqrd2e * lj1[itype][jtype] * qtmp*q[j] / denc;
|
||||
ecoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
|
||||
if (evflag) ev_tally_tip4p(key,vlist,v,ecoul,alpha);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutTIP4PLongSoft::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 9 || narg > 10) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
typeO = force->inumeric(FLERR,arg[0]);
|
||||
typeH = force->inumeric(FLERR,arg[1]);
|
||||
typeB = force->inumeric(FLERR,arg[2]);
|
||||
typeA = force->inumeric(FLERR,arg[3]);
|
||||
qdist = force->numeric(FLERR,arg[4]);
|
||||
nlambda = force->numeric(FLERR,arg[5]);
|
||||
alphalj = force->numeric(FLERR,arg[6]);
|
||||
alphac = force->numeric(FLERR,arg[7]);
|
||||
|
||||
cut_lj_global = force->numeric(FLERR,arg[8]);
|
||||
if (narg == 9) cut_coul = cut_lj_global;
|
||||
else cut_coul = force->numeric(FLERR,arg[9]);
|
||||
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutTIP4PLongSoft::init_style()
|
||||
{
|
||||
if (atom->tag_enable == 0)
|
||||
error->all(FLERR,"Pair style lj/cut/tip4p/long/soft requires atom IDs");
|
||||
if (!force->newton_pair)
|
||||
error->all(FLERR,
|
||||
"Pair style lj/cut/tip4p/long/soft requires newton pair on");
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,
|
||||
"Pair style lj/cut/tip4p/long/soft requires atom attribute q");
|
||||
if (force->bond == NULL)
|
||||
error->all(FLERR,"Must use a bond style with TIP4P potential");
|
||||
if (force->angle == NULL)
|
||||
error->all(FLERR,"Must use an angle style with TIP4P potential");
|
||||
|
||||
PairLJCutCoulLongSoft::init_style();
|
||||
|
||||
// set alpha parameter
|
||||
|
||||
double theta = force->angle->equilibrium_angle(typeA);
|
||||
double blen = force->bond->equilibrium_distance(typeB);
|
||||
alpha = qdist / (cos(0.5*theta) * blen);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutTIP4PLongSoft::init_one(int i, int j)
|
||||
{
|
||||
double cut = PairLJCutCoulLongSoft::init_one(i,j);
|
||||
|
||||
// check that LJ epsilon = 0.0 for water H
|
||||
// set LJ cutoff to 0.0 for any interaction involving water H
|
||||
// so LJ term isn't calculated in compute()
|
||||
|
||||
if ((i == typeH && epsilon[i][i] != 0.0) ||
|
||||
(j == typeH && epsilon[j][j] != 0.0))
|
||||
error->all(FLERR,"Water H epsilon must be 0.0 for "
|
||||
"pair style lj/cut/tip4p/long/soft");
|
||||
|
||||
if (i == typeH || j == typeH)
|
||||
cut_ljsq[j][i] = cut_ljsq[i][j] = 0.0;
|
||||
|
||||
return cut;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutTIP4PLongSoft::write_restart_settings(FILE *fp)
|
||||
{
|
||||
PairLJCutCoulLongSoft::write_restart_settings(fp);
|
||||
|
||||
fwrite(&typeO,sizeof(int),1,fp);
|
||||
fwrite(&typeH,sizeof(int),1,fp);
|
||||
fwrite(&typeB,sizeof(int),1,fp);
|
||||
fwrite(&typeA,sizeof(int),1,fp);
|
||||
fwrite(&qdist,sizeof(double),1,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(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutTIP4PLongSoft::read_restart_settings(FILE *fp)
|
||||
{
|
||||
PairLJCutCoulLongSoft::read_restart_settings(fp);
|
||||
|
||||
if (comm->me == 0) {
|
||||
fread(&typeO,sizeof(int),1,fp);
|
||||
fread(&typeH,sizeof(int),1,fp);
|
||||
fread(&typeB,sizeof(int),1,fp);
|
||||
fread(&typeA,sizeof(int),1,fp);
|
||||
fread(&qdist,sizeof(double),1,fp);
|
||||
|
||||
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(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
MPI_Bcast(&typeO,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&typeH,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&typeB,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&typeA,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&qdist,1,MPI_DOUBLE,0,world);
|
||||
|
||||
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(&tail_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute position xM of fictitious charge site for O atom and 2 H atoms
|
||||
return it as xM
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutTIP4PLongSoft::compute_newsite(double *xO, double *xH1,
|
||||
double *xH2, double *xM)
|
||||
{
|
||||
double delx1 = xH1[0] - xO[0];
|
||||
double dely1 = xH1[1] - xO[1];
|
||||
double delz1 = xH1[2] - xO[2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
|
||||
double delx2 = xH2[0] - xO[0];
|
||||
double dely2 = xH2[1] - xO[1];
|
||||
double delz2 = xH2[2] - xO[2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
|
||||
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
|
||||
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
|
||||
xM[2] = xO[2] + alpha * 0.5 * (delz1 + delz2);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairLJCutTIP4PLongSoft::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 0;
|
||||
if (strcmp(str,"qdist") == 0) return (void *) &qdist;
|
||||
if (strcmp(str,"typeO") == 0) return (void *) &typeO;
|
||||
if (strcmp(str,"typeH") == 0) return (void *) &typeH;
|
||||
if (strcmp(str,"typeA") == 0) return (void *) &typeA;
|
||||
if (strcmp(str,"typeB") == 0) return (void *) &typeB;
|
||||
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
|
||||
dim = 2;
|
||||
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
|
||||
if (strcmp(str,"sigma") == 0) return (void *) sigma;
|
||||
if (strcmp(str,"lambda") == 0) return (void *) lambda;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of hneigh
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutTIP4PLongSoft::memory_usage()
|
||||
{
|
||||
double bytes = maxeatom * sizeof(double);
|
||||
bytes += maxvatom*6 * sizeof(double);
|
||||
bytes += 2 * nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,105 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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/tip4p/long/soft,PairLJCutTIP4PLongSoft)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_LJ_CUT_TIP4P_LONG_SOFT_H
|
||||
#define LMP_PAIR_LJ_CUT_TIP4P_LONG_SOFT_H
|
||||
|
||||
#include "pair_lj_cut_coul_long_soft.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairLJCutTIP4PLongSoft : public PairLJCutCoulLongSoft {
|
||||
public:
|
||||
PairLJCutTIP4PLongSoft(class LAMMPS *);
|
||||
virtual ~PairLJCutTIP4PLongSoft();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
void write_restart_settings(FILE *fp);
|
||||
void read_restart_settings(FILE *fp);
|
||||
void *extract(const char *, int &);
|
||||
virtual double memory_usage();
|
||||
|
||||
protected:
|
||||
int typeH,typeO; // atom types of TIP4P water H and O atoms
|
||||
int typeA,typeB; // angle and bond types of TIP4P water
|
||||
double alpha; // geometric constraint parameter for TIP4P
|
||||
|
||||
int nmax; // info on off-oxygen charge sites
|
||||
int **hneigh; // 0,1 = indices of 2 H associated with O
|
||||
// 2 = 0 if site loc not yet computed, 1 if yes
|
||||
double **newsite; // locations of charge sites
|
||||
|
||||
void compute_newsite(double *, double *, double *, double *);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: TIP4P hydrogen is missing
|
||||
|
||||
The TIP4P pairwise computation failed to find the correct H atom
|
||||
within a water molecule.
|
||||
|
||||
E: TIP4P hydrogen has incorrect atom type
|
||||
|
||||
The TIP4P pairwise computation found an H atom whose type does not
|
||||
agree with the specified H type.
|
||||
|
||||
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: Pair style lj/cut/tip4p/long/soft requires atom IDs
|
||||
|
||||
There are no atom IDs defined in the system and the TIP4P potential
|
||||
requires them to find O,H atoms with a water molecule.
|
||||
|
||||
E: Pair style lj/cut/tip4p/long/soft requires newton pair on
|
||||
|
||||
This is because the computation of constraint forces within a water
|
||||
molecule adds forces to atoms owned by other processors.
|
||||
|
||||
E: Pair style lj/cut/tip4p/long/soft requires atom attribute q
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
E: Must use a bond style with TIP4P potential
|
||||
|
||||
TIP4P potentials assume bond lengths in water are constrained
|
||||
by a fix shake command.
|
||||
|
||||
E: Must use an angle style with TIP4P potential
|
||||
|
||||
TIP4P potentials assume angles in water are constrained by a fix shake
|
||||
command.
|
||||
|
||||
E: Water H epsilon must be 0.0 for pair style lj/cut/tip4p/long/soft
|
||||
|
||||
This is because LAMMPS does not compute the Lennard-Jones interactions
|
||||
with these particles for efficiency reasons.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,514 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 authors: Amalie Frischknecht and Ahmed Ismail (SNL)
|
||||
simpler force assignment added by Rolf Isele-Holder (Aachen University)
|
||||
Soft-core version: Agilio Padua (Univ Blaise Pascal & CNRS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_tip4p_long_soft.h"
|
||||
#include "angle.h"
|
||||
#include "atom.h"
|
||||
#include "bond.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "update.h"
|
||||
#include "respa.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairTIP4PLongSoft::PairTIP4PLongSoft(LAMMPS *lmp) : PairCoulLongSoft(lmp)
|
||||
{
|
||||
tip4pflag = 1;
|
||||
single_enable = 0;
|
||||
respa_enable = 0;
|
||||
|
||||
nmax = 0;
|
||||
hneigh = NULL;
|
||||
newsite = NULL;
|
||||
|
||||
// TIP4P cannot compute virial as F dot r
|
||||
// due to finding bonded H atoms which are not near O atom
|
||||
|
||||
no_virial_fdotr_compute = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairTIP4PLongSoft::~PairTIP4PLongSoft()
|
||||
{
|
||||
memory->destroy(hneigh);
|
||||
memory->destroy(newsite);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairTIP4PLongSoft::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype,key;
|
||||
int n,vlist[6];
|
||||
int iH1,iH2,jH1,jH2;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul;
|
||||
double r,forcecoul,cforce;
|
||||
double factor_coul;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
double denc;
|
||||
double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
|
||||
double *x1,*x2;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
|
||||
ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
// reallocate hneigh & newsite if necessary
|
||||
// initialize hneigh[0] to -1 on steps when reneighboring occurred
|
||||
// initialize hneigh[2] to 0 every step
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
int nall = nlocal + atom->nghost;
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
nmax = atom->nmax;
|
||||
memory->destroy(hneigh);
|
||||
memory->create(hneigh,nmax,3,"pair:hneigh");
|
||||
memory->destroy(newsite);
|
||||
memory->create(newsite,nmax,3,"pair:newsite");
|
||||
}
|
||||
if (neighbor->ago == 0)
|
||||
for (i = 0; i < nall; i++) hneigh[i][0] = -1;
|
||||
for (i = 0; i < nall; i++) hneigh[i][2] = 0;
|
||||
|
||||
double **f = atom->f;
|
||||
double **x = atom->x;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
double *special_coul = force->special_coul;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
double cut_coulsqplus = (cut_coul+2.0*qdist) * (cut_coul+2.0*qdist);
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
|
||||
// if atom I = water O, set x1 = offset charge site
|
||||
// else x1 = x of atom I
|
||||
|
||||
if (itype == typeO) {
|
||||
if (hneigh[i][0] < 0) {
|
||||
hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1);
|
||||
hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2);
|
||||
hneigh[i][2] = 1;
|
||||
if (iH1 == -1 || iH2 == -1)
|
||||
error->one(FLERR,"TIP4P hydrogen is missing");
|
||||
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
|
||||
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
|
||||
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
|
||||
} else {
|
||||
iH1 = hneigh[i][0];
|
||||
iH2 = hneigh[i][1];
|
||||
if (hneigh[i][2] == 0) {
|
||||
hneigh[i][2] = 1;
|
||||
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
|
||||
}
|
||||
}
|
||||
x1 = newsite[i];
|
||||
} else x1 = x[i];
|
||||
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
// adjust rsq and delxyz for off-site O charge(s) if necessary
|
||||
// but only if they are within reach
|
||||
|
||||
if (rsq < cut_coulsqplus) {
|
||||
if (itype == typeO || jtype == typeO) {
|
||||
|
||||
// if atom J = water O, set x2 = offset charge site
|
||||
// else x2 = x of atom J
|
||||
|
||||
if (jtype == typeO) {
|
||||
if (hneigh[j][0] < 0) {
|
||||
hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
|
||||
hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);
|
||||
hneigh[j][2] = 1;
|
||||
if (jH1 == -1 || jH2 == -1)
|
||||
error->one(FLERR,"TIP4P hydrogen is missing");
|
||||
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
|
||||
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
|
||||
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
|
||||
} else {
|
||||
jH1 = hneigh[j][0];
|
||||
jH2 = hneigh[j][1];
|
||||
if (hneigh[j][2] == 0) {
|
||||
hneigh[j][2] = 1;
|
||||
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
|
||||
}
|
||||
}
|
||||
x2 = newsite[j];
|
||||
} else x2 = x[j];
|
||||
|
||||
delx = x1[0] - x2[0];
|
||||
dely = x1[1] - x2[1];
|
||||
delz = x1[2] - x2[2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
}
|
||||
|
||||
// Coulombic interaction based on modified 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;
|
||||
|
||||
denc = sqrt(lam2[itype][jtype] + rsq);
|
||||
prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
|
||||
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) {
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
|
||||
cforce = forcecoul;
|
||||
|
||||
// if i,j are not O atoms, force is applied directly
|
||||
// if i or j are O atoms, force is on fictitious atom & partitioned
|
||||
// force partitioning due to Feenstra, J Comp Chem, 20, 786 (1999)
|
||||
// f_f = fictitious force, fO = f_f (1 - 2 alpha), fH = alpha f_f
|
||||
// preserves total force and torque on water molecule
|
||||
// virial = sum(r x F) where each water's atoms are near xi and xj
|
||||
// vlist stores 2,4,6 atoms whose forces contribute to virial
|
||||
|
||||
n = 0;
|
||||
key = 0;
|
||||
|
||||
if (itype != typeO) {
|
||||
f[i][0] += delx * cforce;
|
||||
f[i][1] += dely * cforce;
|
||||
f[i][2] += delz * cforce;
|
||||
|
||||
if (vflag) {
|
||||
v[0] = x[i][0] * delx * cforce;
|
||||
v[1] = x[i][1] * dely * cforce;
|
||||
v[2] = x[i][2] * delz * cforce;
|
||||
v[3] = x[i][0] * dely * cforce;
|
||||
v[4] = x[i][0] * delz * cforce;
|
||||
v[5] = x[i][1] * delz * cforce;
|
||||
}
|
||||
vlist[n++] = i;
|
||||
|
||||
} else {
|
||||
key++;
|
||||
|
||||
fd[0] = delx*cforce;
|
||||
fd[1] = dely*cforce;
|
||||
fd[2] = delz*cforce;
|
||||
|
||||
fO[0] = fd[0]*(1 - alpha);
|
||||
fO[1] = fd[1]*(1 - alpha);
|
||||
fO[2] = fd[2]*(1 - alpha);
|
||||
|
||||
fH[0] = 0.5 * alpha * fd[0];
|
||||
fH[1] = 0.5 * alpha * fd[1];
|
||||
fH[2] = 0.5 * alpha * fd[2];
|
||||
|
||||
f[i][0] += fO[0];
|
||||
f[i][1] += fO[1];
|
||||
f[i][2] += fO[2];
|
||||
|
||||
f[iH1][0] += fH[0];
|
||||
f[iH1][1] += fH[1];
|
||||
f[iH1][2] += fH[2];
|
||||
|
||||
f[iH2][0] += fH[0];
|
||||
f[iH2][1] += fH[1];
|
||||
f[iH2][2] += fH[2];
|
||||
|
||||
if (vflag) {
|
||||
domain->closest_image(x[i],x[iH1],xH1);
|
||||
domain->closest_image(x[i],x[iH2],xH2);
|
||||
|
||||
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
|
||||
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
|
||||
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
|
||||
v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
|
||||
v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
|
||||
v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
|
||||
}
|
||||
vlist[n++] = i;
|
||||
vlist[n++] = iH1;
|
||||
vlist[n++] = iH2;
|
||||
}
|
||||
|
||||
if (jtype != typeO) {
|
||||
f[j][0] -= delx * cforce;
|
||||
f[j][1] -= dely * cforce;
|
||||
f[j][2] -= delz * cforce;
|
||||
|
||||
if (vflag) {
|
||||
v[0] -= x[j][0] * delx * cforce;
|
||||
v[1] -= x[j][1] * dely * cforce;
|
||||
v[2] -= x[j][2] * delz * cforce;
|
||||
v[3] -= x[j][0] * dely * cforce;
|
||||
v[4] -= x[j][0] * delz * cforce;
|
||||
v[5] -= x[j][1] * delz * cforce;
|
||||
}
|
||||
vlist[n++] = j;
|
||||
|
||||
} else {
|
||||
key += 2;
|
||||
|
||||
fd[0] = -delx*cforce;
|
||||
fd[1] = -dely*cforce;
|
||||
fd[2] = -delz*cforce;
|
||||
|
||||
fO[0] = fd[0]*(1 - alpha);
|
||||
fO[1] = fd[1]*(1 - alpha);
|
||||
fO[2] = fd[2]*(1 - alpha);
|
||||
|
||||
fH[0] = 0.5 * alpha * fd[0];
|
||||
fH[1] = 0.5 * alpha * fd[1];
|
||||
fH[2] = 0.5 * alpha * fd[2];
|
||||
|
||||
f[j][0] += fO[0];
|
||||
f[j][1] += fO[1];
|
||||
f[j][2] += fO[2];
|
||||
|
||||
f[jH1][0] += fH[0];
|
||||
f[jH1][1] += fH[1];
|
||||
f[jH1][2] += fH[2];
|
||||
|
||||
f[jH2][0] += fH[0];
|
||||
f[jH2][1] += fH[1];
|
||||
f[jH2][2] += fH[2];
|
||||
|
||||
if (vflag) {
|
||||
domain->closest_image(x[j],x[jH1],xH1);
|
||||
domain->closest_image(x[j],x[jH2],xH2);
|
||||
|
||||
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
|
||||
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
|
||||
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
|
||||
v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
|
||||
v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
|
||||
v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
|
||||
}
|
||||
vlist[n++] = j;
|
||||
vlist[n++] = jH1;
|
||||
vlist[n++] = jH2;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / denc;
|
||||
ecoul = prefactor*erfc;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
|
||||
if (evflag) ev_tally_tip4p(key,vlist,v,ecoul,alpha);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairTIP4PLongSoft::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 8) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
typeO = force->inumeric(FLERR,arg[0]);
|
||||
typeH = force->inumeric(FLERR,arg[1]);
|
||||
typeB = force->inumeric(FLERR,arg[2]);
|
||||
typeA = force->inumeric(FLERR,arg[3]);
|
||||
qdist = force->numeric(FLERR,arg[4]);
|
||||
|
||||
nlambda = force->numeric(FLERR,arg[5]);
|
||||
alphac = force->numeric(FLERR,arg[6]);
|
||||
|
||||
cut_coul = force->numeric(FLERR,arg[7]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairTIP4PLongSoft::init_style()
|
||||
{
|
||||
if (atom->tag_enable == 0)
|
||||
error->all(FLERR,"Pair style tip4p/long requires atom IDs");
|
||||
if (!force->newton_pair)
|
||||
error->all(FLERR,
|
||||
"Pair style tip4p/long requires newton pair on");
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,
|
||||
"Pair style tip4p/long requires atom attribute q");
|
||||
if (force->bond == NULL)
|
||||
error->all(FLERR,"Must use a bond style with TIP4P potential");
|
||||
if (force->angle == NULL)
|
||||
error->all(FLERR,"Must use an angle style with TIP4P potential");
|
||||
|
||||
PairCoulLongSoft::init_style();
|
||||
|
||||
// set alpha parameter
|
||||
|
||||
double theta = force->angle->equilibrium_angle(typeA);
|
||||
double blen = force->bond->equilibrium_distance(typeB);
|
||||
alpha = qdist / (cos(0.5*theta) * blen);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairTIP4PLongSoft::init_one(int i, int j)
|
||||
{
|
||||
return PairCoulLongSoft::init_one(i,j);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairTIP4PLongSoft::write_restart_settings(FILE *fp)
|
||||
{
|
||||
PairCoulLongSoft::write_restart_settings(fp);
|
||||
|
||||
fwrite(&typeO,sizeof(int),1,fp);
|
||||
fwrite(&typeH,sizeof(int),1,fp);
|
||||
fwrite(&typeB,sizeof(int),1,fp);
|
||||
fwrite(&typeA,sizeof(int),1,fp);
|
||||
fwrite(&qdist,sizeof(double),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairTIP4PLongSoft::read_restart_settings(FILE *fp)
|
||||
{
|
||||
PairCoulLongSoft::read_restart_settings(fp);
|
||||
|
||||
if (comm->me == 0) {
|
||||
fread(&typeO,sizeof(int),1,fp);
|
||||
fread(&typeH,sizeof(int),1,fp);
|
||||
fread(&typeB,sizeof(int),1,fp);
|
||||
fread(&typeA,sizeof(int),1,fp);
|
||||
fread(&qdist,sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&typeO,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&typeH,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&typeB,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&typeA,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&qdist,1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute position xM of fictitious charge site for O atom and 2 H atoms
|
||||
return it as xM
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairTIP4PLongSoft::compute_newsite(double *xO, double *xH1,
|
||||
double *xH2, double *xM)
|
||||
{
|
||||
double delx1 = xH1[0] - xO[0];
|
||||
double dely1 = xH1[1] - xO[1];
|
||||
double delz1 = xH1[2] - xO[2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
|
||||
double delx2 = xH2[0] - xO[0];
|
||||
double dely2 = xH2[1] - xO[1];
|
||||
double delz2 = xH2[2] - xO[2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
|
||||
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
|
||||
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
|
||||
xM[2] = xO[2] + alpha * 0.5 * (delz1 + delz2);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairTIP4PLongSoft::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 0;
|
||||
if (strcmp(str,"qdist") == 0) return (void *) &qdist;
|
||||
if (strcmp(str,"typeO") == 0) return (void *) &typeO;
|
||||
if (strcmp(str,"typeH") == 0) return (void *) &typeH;
|
||||
if (strcmp(str,"typeA") == 0) return (void *) &typeA;
|
||||
if (strcmp(str,"typeB") == 0) return (void *) &typeB;
|
||||
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
|
||||
dim = 2;
|
||||
if (strcmp(str,"lambda") == 0) return (void *) lambda;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of hneigh
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairTIP4PLongSoft::memory_usage()
|
||||
{
|
||||
double bytes = maxeatom * sizeof(double);
|
||||
bytes += maxvatom*6 * sizeof(double);
|
||||
bytes += 2 * nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,100 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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(tip4p/long/soft,PairTIP4PLongSoft)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_TIP4P_LONG_SOFT_H
|
||||
#define LMP_PAIR_TIP4P_LONG_SOFT_H
|
||||
|
||||
#include "pair_coul_long_soft.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairTIP4PLongSoft : public PairCoulLongSoft {
|
||||
public:
|
||||
PairTIP4PLongSoft(class LAMMPS *);
|
||||
virtual ~PairTIP4PLongSoft();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
void write_restart_settings(FILE *fp);
|
||||
void read_restart_settings(FILE *fp);
|
||||
void *extract(const char *, int &);
|
||||
virtual double memory_usage();
|
||||
|
||||
protected:
|
||||
int typeH,typeO; // atom types of TIP4P water H and O atoms
|
||||
int typeA,typeB; // angle and bond types of TIP4P water
|
||||
double alpha; // geometric constraint parameter for TIP4P
|
||||
|
||||
int nmax; // info on off-oxygen charge sites
|
||||
int **hneigh; // 0,1 = indices of 2 H associated with O
|
||||
// 2 = 0 if site loc not yet computed, 1 if yes
|
||||
double **newsite; // locations of charge sites
|
||||
|
||||
void compute_newsite(double *, double *, double *, double *);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: TIP4P hydrogen is missing
|
||||
|
||||
The TIP4P pairwise computation failed to find the correct H atom
|
||||
within a water molecule.
|
||||
|
||||
E: TIP4P hydrogen has incorrect atom type
|
||||
|
||||
The TIP4P pairwise computation found an H atom whose type does not
|
||||
agree with the specified H type.
|
||||
|
||||
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: Pair style tip4p/long/soft requires atom IDs
|
||||
|
||||
There are no atom IDs defined in the system and the TIP4P potential
|
||||
requires them to find O,H atoms with a water molecule.
|
||||
|
||||
E: Pair style tip4p/long/soft requires newton pair on
|
||||
|
||||
This is because the computation of constraint forces within a water
|
||||
molecule adds forces to atoms owned by other processors.
|
||||
|
||||
E: Pair style tip4p/long/soft requires atom attribute q
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
E: Must use a bond style with TIP4P potential
|
||||
|
||||
TIP4P potentials assume bond lengths in water are constrained
|
||||
by a fix shake command.
|
||||
|
||||
E: Must use an angle style with TIP4P potential
|
||||
|
||||
TIP4P potentials assume angles in water are constrained by a fix shake
|
||||
command.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue