jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2016-06-17 23:48:15 +00:00
parent b161fbb52a
commit d55f968432
112 changed files with 1525 additions and 1350 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/KOKKOS/verlet_kokkos.cpp
View File

@ -68,9 +68,10 @@ void VerletKokkos::setup()
{
if (comm->me == 0 && screen) {
fprintf(screen,"Setting up Verlet run ...\n");
fprintf(screen," Unit style : %s\n", update->unit_style);
fprintf(screen," Current step: " BIGINT_FORMAT "\n", update->ntimestep);
fprintf(screen," Time step : %g\n", update->dt);
fprintf(screen," Unit style : %s\n", update->unit_style);
fprintf(screen," Current step : " BIGINT_FORMAT "\n", update->ntimestep);
fprintf(screen," Time step : %g\n", update->dt);
timer->print_timeout(screen);
}
update->setupflag = 1;
@ -294,12 +295,16 @@ void VerletKokkos::run(int n)
f_merge_copy = DAT::t_f_array("VerletKokkos::f_merge_copy",atomKK->k_f.dimension_0());
static double time = 0.0;
static int count = 0;
atomKK->sync(Device,ALL_MASK);
Kokkos::Impl::Timer ktimer;
timer->init_timeout();
for (int i = 0; i < n; i++) {
if (timer->check_timeout(i)) {
update->nsteps = i;
break;
}
ntimestep = ++update->ntimestep;
ev_set(ntimestep);
@ -383,7 +388,6 @@ void VerletKokkos::run(int n)
unsigned int datamask_read_device = 0;
unsigned int datamask_modify_device = 0;
unsigned int datamask_read_host = 0;
unsigned int datamask_modify_host = 0;
if ( pair_compute_flag ) {
if (force->pair->execution_space==Host) {

16
src/MANYBODY/fix_qeq_comb.cpp
View File

@ -45,6 +45,8 @@ FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
peratom_flag = 1;
size_peratom_cols = 0;
peratom_freq = 1;
respa_level_support = 1;
ilevel_respa = 0;
nevery = force->inumeric(FLERR,arg[3]);
precision = force->numeric(FLERR,arg[4]);
@ -125,8 +127,10 @@ void FixQEQComb::init()
if (comb == NULL && comb3 == NULL)
error->all(FLERR,"Must use pair_style comb or comb3 with fix qeq/comb");
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");
@ -140,9 +144,9 @@ void FixQEQComb::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
firstflag = 0;
}
@ -275,7 +279,7 @@ void FixQEQComb::post_force(int vflag)
void FixQEQComb::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ----------------------------------------------------------------------

2
src/MANYBODY/fix_qeq_comb.h
View File

@ -43,7 +43,7 @@ class FixQEQComb : public Fix {
protected:
int me,firstflag;
double precision;
int nlevels_respa;
int ilevel_respa;
bigint ngroup;
FILE *fp;

16
src/MISC/fix_efield.cpp
View File

@ -52,6 +52,8 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extvector = 1;
extscalar = 1;
respa_level_support = 1;
ilevel_respa = 0;
qe2f = force->qe2f;
xstr = ystr = zstr = NULL;
@ -216,8 +218,10 @@ void FixEfield::init()
update->whichflag == 2 && estyle == NONE)
error->all(FLERR,"Must use variable energy with fix efield");
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -227,9 +231,9 @@ void FixEfield::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -393,7 +397,7 @@ void FixEfield::post_force(int vflag)
void FixEfield::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/MISC/fix_efield.h
View File

@ -45,7 +45,7 @@ class FixEfield : public Fix {
char *xstr,*ystr,*zstr,*estr;
char *idregion;
int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
int nlevels_respa;
int ilevel_respa;
double qe2f;
int qflag,muflag;

16
src/MISC/fix_orient_bcc.cpp
View File

@ -72,6 +72,8 @@ FixOrientBCC::FixOrientBCC(LAMMPS *lmp, int narg, char **arg) :
peratom_flag = 1;
size_peratom_cols = 2;
peratom_freq = 1;
respa_level_support = 1;
ilevel_respa = 0;
nstats = force->inumeric(FLERR,arg[3]);
direction_of_motion = force->inumeric(FLERR,arg[4]);
@ -210,8 +212,10 @@ int FixOrientBCC::setmask()
void FixOrientBCC::init()
{
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
// need a full neighbor list, built whenever re-neighboring occurs
@ -236,9 +240,9 @@ void FixOrientBCC::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -467,7 +471,7 @@ void FixOrientBCC::post_force(int vflag)
void FixOrientBCC::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/MISC/fix_orient_bcc.h
View File

@ -58,7 +58,7 @@ class FixOrientBCC : public Fix {
private:
int me;
int nlevels_respa;
int ilevel_respa;
int direction_of_motion; // 1 = center shrinks, 0 = center grows
int nstats; // stats output every this many steps

16
src/MISC/fix_orient_fcc.cpp
View File

@ -69,6 +69,8 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
peratom_flag = 1;
size_peratom_cols = 2;
peratom_freq = 1;
respa_level_support = 1;
ilevel_respa = 0;
nstats = force->inumeric(FLERR,arg[3]);
direction_of_motion = force->inumeric(FLERR,arg[4]);
@ -207,8 +209,10 @@ int FixOrientFCC::setmask()
void FixOrientFCC::init()
{
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
// need a full neighbor list, built whenever re-neighboring occurs
@ -233,9 +237,9 @@ void FixOrientFCC::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -464,7 +468,7 @@ void FixOrientFCC::post_force(int vflag)
void FixOrientFCC::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/MISC/fix_orient_fcc.h
View File

@ -58,7 +58,7 @@ class FixOrientFCC : public Fix {
private:
int me;
int nlevels_respa;
int ilevel_respa;
int direction_of_motion; // 1 = center shrinks, 0 = center grows
int nstats; // stats output every this many steps

8
src/RIGID/fix_ehex.cpp
View File

@ -244,13 +244,11 @@ void FixEHEX::end_of_step() {
------------------------------------------------------------------------- */
void FixEHEX::rescale() {
double heat,ke,massone;
double Kr, Ke, escale;
double vsub[3],vcm[3], sfr[3];
double dvcmdt[3];
double mi;
double dt;
double F, mr, Fo2Kr, epsr_ik, sfvr, eta_ik;
double F, mr, epsr_ik, sfvr, eta_ik;
dt = update->dt;
@ -266,8 +264,6 @@ void FixEHEX::rescale() {
mr = masstotal;
Fo2Kr = F / (2.*Kr);
// energy scaling factor
escale = 1. + (F*dt)/Kr;
@ -406,7 +402,7 @@ void FixEHEX::update_scalingmask() {
inside the region.
------------------------------------------------------------------------- */
bool FixEHEX::check_cluster(int *shake_atom, int n, Region * region) {
bool FixEHEX::check_cluster(tagint *shake_atom, int n, Region * region) {
// IMPORTANT NOTE: If any site of the cluster belongs to a group
// which should not be rescaled than all of the sites

2
src/RIGID/fix_ehex.h
View File

@ -43,7 +43,7 @@ class FixEHEX : public Fix {
void com_properties(double *, double *, double *, double*, double *, double*);
bool rescale_atom(int i, Region*region);
virtual void grow_arrays(int nmax);
bool check_cluster(int *shake_atom, int n, Region * region);
bool check_cluster(tagint *shake_atom, int n, Region * region);
private:
int iregion;

11
src/RIGID/fix_rattle.cpp
View File

@ -240,7 +240,6 @@ void FixRattle::final_integrate_respa(int ilevel, int iloop)
void FixRattle::vrattle3angle(int m)
{
tagint i0,i1,i2;
int nlist,list[3];
double c[3], l[3], a[3][3], r01[3], imass[3],
r02[3], r12[3], vp01[3], vp02[3], vp12[3];
@ -323,7 +322,6 @@ void FixRattle::vrattle3angle(int m)
void FixRattle::vrattle2(int m)
{
tagint i0, i1;
int nlist,list[2];
double imass[2], r01[3], vp01[3];
// local atom IDs and constraint distances
@ -372,7 +370,6 @@ void FixRattle::vrattle2(int m)
void FixRattle::vrattle3(int m)
{
tagint i0,i1,i2;
int nlist,list[3];
double imass[3], r01[3], r02[3], vp01[3], vp02[3],
a[2][2],c[2],l[2];
@ -442,7 +439,6 @@ void FixRattle::vrattle3(int m)
void FixRattle::vrattle4(int m)
{
tagint i0,i1,i2,i3;
int nlist,list[4];
double imass[4], c[3], l[3], a[3][3],
r01[3], r02[3], r03[3], vp01[3], vp02[3], vp03[3];
@ -592,7 +588,6 @@ void FixRattle::solve3x3exactly(const double a[][3],
void FixRattle::reset_dt()
{
FixShake::reset_dt();
dtfv = 0.5 * update->dt * force->ftm2v;
}
/* ----------------------------------------------------------------------
@ -601,8 +596,8 @@ void FixRattle::reset_dt()
void FixRattle::update_v_half_nocons()
{
const double dtfv = 0.5 * update->dt * force->ftm2v;
double dtfvinvm;
dtfv = 0.5 * update->dt * force->ftm2v;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
@ -632,10 +627,6 @@ void FixRattle::update_v_half_nocons()
void FixRattle::update_v_half_nocons_respa(int ilevel)
{
// select timestep for current level
dtfv = 0.5 * step_respa[ilevel] * force->ftm2v;
// carry out unconstrained velocity update
update_v_half_nocons();

1
src/RIGID/fix_rattle.h
View File

@ -28,7 +28,6 @@ namespace LAMMPS_NS {
class FixRattle : public FixShake {
public:
double **vp; // array for unconstrained velocities
double dtfv; // timestep for velocity update
int comm_mode; // mode for communication pack/unpack
double derr_max; // distance error
double verr_max; // velocity error

111
src/USER-DPD/atom_vec_dpd.cpp
View File

@ -76,9 +76,8 @@ void AtomVecDPD::grow(int n)
uChem = memory->grow(atom->uChem,nmax,"atom:uChem");
uCG = memory->grow(atom->uCG,nmax,"atom:uCG");
uCGnew = memory->grow(atom->uCGnew,nmax,"atom:uCGnew");
duCond = memory->grow(atom->duCond,nmax,"atom:duCond");
duMech = memory->grow(atom->duMech,nmax,"atom:duMech");
duChem = memory->grow(atom->duChem,nmax,"atom:duChem");
ssaAIR = memory->grow(atom->ssaAIR,nmax,"atom:ssaAIR");
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -101,9 +100,8 @@ void AtomVecDPD::grow_reset()
uChem = atom->uChem;
uCG = atom->uCG;
uCGnew = atom->uCGnew;
duCond = atom->duCond;
duMech = atom->duMech;
duChem = atom->duChem;
ssaAIR = atom->ssaAIR;
}
/* ----------------------------------------------------------------------
@ -128,6 +126,7 @@ void AtomVecDPD::copy(int i, int j, int delflag)
uChem[j] = uChem[i];
uCG[j] = uCG[i];
uCGnew[j] = uCGnew[i];
ssaAIR[j] = ssaAIR[i];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -220,10 +219,10 @@ int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
buf[m++] = dpdTheta[j];
buf[m++] = uCond[j];
buf[m++] = uMech[j];
buf[m++] = uChem[j];
buf[m++] = dpdTheta[j];
buf[m++] = uCond[j];
buf[m++] = uMech[j];
buf[m++] = uChem[j];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
@ -243,10 +242,10 @@ int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
buf[m++] = dpdTheta[j];
buf[m++] = uCond[j];
buf[m++] = uMech[j];
buf[m++] = uChem[j];
buf[m++] = dpdTheta[j];
buf[m++] = uCond[j];
buf[m++] = uMech[j];
buf[m++] = uChem[j];
}
}
}
@ -428,18 +427,18 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf,
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
buf[m++] = dpdTheta[j];
buf[m++] = uCond[j];
buf[m++] = uMech[j];
buf[m++] = uChem[j];
buf[m++] = uCG[j];
buf[m++] = uCGnew[j];
buf[m++] = dpdTheta[j];
buf[m++] = uCond[j];
buf[m++] = uMech[j];
buf[m++] = uChem[j];
buf[m++] = uCG[j];
buf[m++] = uCGnew[j];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
@ -450,9 +449,9 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf,
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
@ -462,12 +461,12 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf,
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
buf[m++] = dpdTheta[j];
buf[m++] = uCond[j];
buf[m++] = uMech[j];
buf[m++] = uChem[j];
buf[m++] = uCG[j];
buf[m++] = uCGnew[j];
buf[m++] = dpdTheta[j];
buf[m++] = uCond[j];
buf[m++] = uMech[j];
buf[m++] = uChem[j];
buf[m++] = uCG[j];
buf[m++] = uCGnew[j];
}
}
}
@ -517,6 +516,41 @@ int AtomVecDPD::pack_border_hybrid(int n, int *list, double *buf)
return m;
}
/* ----------------------------------------------------------------------
convert atom coords into the ssa active interaction region number
------------------------------------------------------------------------- */
int AtomVecDPD::coord2ssaAIR(double *x)
{
int ix, iy, iz;
ix = iy = iz = 0;
if (x[2] < domain->sublo[2]) iz = -1;
if (x[2] > domain->subhi[2]) iz = 1;
if (x[1] < domain->sublo[1]) iy = -1;
if (x[1] > domain->subhi[1]) iy = 1;
if (x[0] < domain->sublo[0]) ix = -1;
if (x[0] > domain->subhi[0]) ix = 1;
if(iz < 0){
return 0;
} else if(iz == 0){
if( iy<0 ) return 0; // bottom left/middle/right
if( (iy==0) && (ix<0) ) return 0; // left atoms
if( (iy==0) && (ix==0) ) return 1; // Locally owned atoms
if( (iy==0) && (ix>0) ) return 3; // Right atoms
if( (iy>0) && (ix==0) ) return 2; // Top-middle atoms
if( (iy>0) && (ix!=0) ) return 4; // Top-right and top-left atoms
} else { // iz > 0
if((ix==0) && (iy==0)) return 5; // Back atoms
if((ix==0) && (iy!=0)) return 6; // Top-back and bottom-back atoms
if((ix!=0) && (iy==0)) return 7; // Left-back and right-back atoms
if((ix!=0) && (iy!=0)) return 8; // Back corner atoms
}
return 0;
}
/* ---------------------------------------------------------------------- */
void AtomVecDPD::unpack_border(int n, int first, double *buf)
@ -539,6 +573,7 @@ void AtomVecDPD::unpack_border(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
ssaAIR[i] = coord2ssaAIR(x[i]);
}
if (atom->nextra_border)
@ -572,6 +607,7 @@ void AtomVecDPD::unpack_border_vel(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
ssaAIR[i] = coord2ssaAIR(x[i]);
}
if (atom->nextra_border)
@ -612,6 +648,7 @@ int AtomVecDPD::unpack_border_hybrid(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
ssaAIR[i] = coord2ssaAIR(x[i]);
}
return m;
}
@ -673,6 +710,7 @@ int AtomVecDPD::unpack_exchange(double *buf)
uChem[nlocal] = buf[m++];
uCG[nlocal] = buf[m++];
uCGnew[nlocal] = buf[m++];
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -763,6 +801,7 @@ int AtomVecDPD::unpack_restart(double *buf)
uCond[nlocal] = buf[m++];
uMech[nlocal] = buf[m++];
uChem[nlocal] = buf[m++];
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
double **extra = atom->extra;
if (atom->nextra_store) {
@ -802,9 +841,8 @@ void AtomVecDPD::create_atom(int itype, double *coord)
uChem[nlocal] = 0.0;
uCG[nlocal] = 0.0;
uCGnew[nlocal] = 0.0;
duCond[nlocal] = 0.0;
duMech[nlocal] = 0.0;
duChem[nlocal] = 0.0;
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
atom->nlocal++;
}
@ -845,6 +883,7 @@ void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values)
uChem[nlocal] = 0.0;
uCG[nlocal] = 0.0;
uCGnew[nlocal] = 0.0;
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
atom->nlocal++;
}
@ -857,6 +896,7 @@ void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values)
int AtomVecDPD::data_atom_hybrid(int nlocal, char **values)
{
dpdTheta[nlocal] = atof(values[0]);
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
return 1;
}
@ -900,9 +940,9 @@ void AtomVecDPD::write_data(FILE *fp, int n, double **buf)
for (int i = 0; i < n; i++)
fprintf(fp,TAGINT_FORMAT " %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
(tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
buf[i][2],buf[i][3],buf[i][4],buf[i][5],
buf[i][2],buf[i][3],buf[i][4],buf[i][5],
(int) ubuf(buf[i][6]).i,(int) ubuf(buf[i][7]).i,
(int) ubuf(buf[i][8]).i);
(int) ubuf(buf[i][8]).i);
}
/* ----------------------------------------------------------------------
@ -937,9 +977,8 @@ bigint AtomVecDPD::memory_usage()
if (atom->memcheck("uChem")) bytes += memory->usage(uChem,nmax);
if (atom->memcheck("uCG")) bytes += memory->usage(uCG,nmax);
if (atom->memcheck("uCGnew")) bytes += memory->usage(uCGnew,nmax);
if (atom->memcheck("duCond")) bytes += memory->usage(duCond,nmax);
if (atom->memcheck("duMech")) bytes += memory->usage(duMech,nmax);
if (atom->memcheck("duChem")) bytes += memory->usage(duChem,nmax);
if (atom->memcheck("ssaAIR")) bytes += memory->usage(ssaAIR,nmax);
return bytes;
}

7
src/USER-DPD/atom_vec_dpd.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -59,13 +59,16 @@ class AtomVecDPD : public AtomVec {
int write_data_hybrid(FILE *, double *);
bigint memory_usage();
double *uCond,*uMech,*uChem,*uCG,*uCGnew,*rho,*dpdTheta;
double *duCond,*duMech,*duChem;
double *duChem;
int *ssaAIR; // Shardlow Splitting Algorithm Active Interaction Region number
protected:
tagint *tag;
int *type,*mask;
imageint *image;
double **x,**v,**f;
int coord2ssaAIR(double *); // map atom coord to an AIR number
};
}

8
src/USER-DPD/compute_dpd.cpp
View File

@ -36,9 +36,9 @@ ComputeDpd::ComputeDpd(LAMMPS *lmp, int narg, char **arg) :
vector_flag = 1;
size_vector = 5;
extvector = 0;
vector = new double[size_vector];
if (atom->dpd_flag != 1) error->all(FLERR,"compute dpd requires atom_style with internal temperature and energies (e.g. dpd)");
}
@ -67,7 +67,7 @@ void ComputeDpd::compute_vector()
dpdU = new double[size_vector];
for (int i = 0; i < size_vector; i++) dpdU[i] = double(0.0);
for (int i = 0; i < size_vector; i++) dpdU[i] = 0.0;
for (int i = 0; i < nlocal; i++){
if (mask[i] & groupbit){
@ -78,7 +78,7 @@ void ComputeDpd::compute_vector()
dpdU[4]++;
}
}
MPI_Allreduce(dpdU,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
natoms = vector[4];

4
src/USER-DPD/compute_dpd.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -51,6 +51,6 @@ command-line option when running LAMMPS to see the offending line.
E: compute dpd requires atom_style with internal temperature and energies (e.g. dpd)
Self-explanatory.
Self-explanatory.
*/

2
src/USER-DPD/compute_dpd_atom.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov

20
src/USER-DPD/fix_eos_cv.cpp
View File

@ -32,7 +32,7 @@ FixEOScv::FixEOScv(LAMMPS *lmp, int narg, char **arg) :
{
if (narg != 4) error->all(FLERR,"Illegal fix eos/cv command");
cvEOS = force->numeric(FLERR,arg[3]);
if(cvEOS <= double(0.0)) error->all(FLERR,"EOS cv must be > 0.0");
if(cvEOS <= 0.0) error->all(FLERR,"EOS cv must be > 0.0");
restart_peratom = 1;
nevery = 1;
@ -61,14 +61,14 @@ void FixEOScv::init()
if(this->restart_reset){
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
uCond[i] = double(0.5)*cvEOS*dpdTheta[i];
uMech[i] = double(0.5)*cvEOS*dpdTheta[i];
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
uCond[i] = 0.5*cvEOS*dpdTheta[i];
uMech[i] = 0.5*cvEOS*dpdTheta[i];
}
}
}
@ -86,8 +86,8 @@ void FixEOScv::post_integrate()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
}
}
@ -104,7 +104,7 @@ void FixEOScv::end_of_step()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
dpdTheta[i] = (uCond[i]+uMech[i])/cvEOS;
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
}
}

2
src/USER-DPD/fix_eos_cv.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov

26
src/USER-DPD/fix_eos_table.cpp
View File

@ -112,15 +112,15 @@ void FixEOStable::init()
if(this->restart_reset){
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
energy_lookup(dpdTheta[i],tmp);
uCond[i] = tmp / double(2.0);
uMech[i] = tmp / double(2.0);
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
energy_lookup(dpdTheta[i],tmp);
uCond[i] = tmp / 2.0;
uMech[i] = tmp / 2.0;
}
}
}
@ -138,8 +138,8 @@ void FixEOStable::post_integrate()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
}
}
@ -156,8 +156,8 @@ void FixEOStable::end_of_step()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
temperature_lookup(uCond[i]+uMech[i],dpdTheta[i]);
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
}
}
@ -195,13 +195,13 @@ void FixEOStable::read_table(Table *tb, Table *tb2, char *file, char *keyword)
// open file
FILE *fp = fopen(file,"r");
FILE *fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
// loop until section found with matching keyword
while (1) {
@ -426,7 +426,7 @@ void FixEOStable::temperature_lookup(double u, double &t)
double fraction;
Table *tb = &tables[1];
if(u < tb->lo || u > tb->hi){
if(u < tb->lo || u > tb->hi){
printf("Energy=%lf TableMin=%lf TableMax=%lf\n",u,tb->lo,tb->hi);
error->one(FLERR,"Energy is not within table cutoffs");
}

8
src/USER-DPD/fix_eos_table.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -99,7 +99,7 @@ Self-explanatory. EOS may not be valid under current simulation conditions.
E: Cannot open file %s
The specified file cannot be opened. Check that the path and name are
correct.
correct.
E: Did not find keyword in table file
@ -115,12 +115,12 @@ List of fix eos/table parameters must include N setting.
E: Temperature is not within table cutoffs
The internal temperature does not lie with the minimum
The internal temperature does not lie with the minimum
and maximum temperature cutoffs of the table
E: Energy is not within table cutoffs
The internal energy does not lie with the minimum
The internal energy does not lie with the minimum
and maximum energy cutoffs of the table
*/

88
src/USER-DPD/fix_eos_table_rx.cpp
View File

@ -57,6 +57,8 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
ntables = 0;
tables = NULL;
tables2 = NULL;
eosSpecies = NULL;
int me;
MPI_Comm_rank(world,&me);
for (int ii=0;ii<nspecies;ii++){
@ -70,7 +72,7 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
null_table(tb);
null_table(tb2);
ntables++;
}
@ -79,7 +81,7 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
Table *tb2 = &tables2[ntables];
if (me == 0) read_table(tb,tb2,arg[4],arg[6]);
for (int ii=0;ii<nspecies;ii++){
Table *tb = &tables[ntables];
Table *tb2 = &tables2[ntables];
@ -88,21 +90,21 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
bcast_table(tb2);
// error check on table parameters
if (tb->ninput <= 1) error->one(FLERR,"Invalid eos/table/rx length");
tb->lo = tb->rfile[0];
tb->hi = tb->rfile[tb->ninput-1];
if (tb->lo >= tb->hi) error->all(FLERR,"eos/table/rx values are not increasing");
if (tb2->ninput <= 1) error->one(FLERR,"Invalid eos/table/rx length");
tb2->lo = tb2->rfile[0];
tb2->hi = tb2->rfile[tb2->ninput-1];
if (tb2->lo >= tb2->hi) error->all(FLERR,"eos/table/rx values are not increasing");
// spline read-in and compute r,e,f vectors within table
spline_table(tb);
compute_table(tb);
spline_table(tb2);
@ -129,6 +131,7 @@ FixEOStableRX::~FixEOStableRX()
memory->sfree(tables2);
delete [] dHf;
delete [] eosSpecies;
}
/* ---------------------------------------------------------------------- */
@ -159,8 +162,8 @@ void FixEOStableRX::setup(int vflag)
if (mask[i] & groupbit){
duChem = uCG[i] - uCGnew[i];
uChem[i] += duChem;
uCG[i] = double(0.0);
uCGnew[i] = double(0.0);
uCG[i] = 0.0;
uCGnew[i] = 0.0;
}
// Communicate the updated momenta and velocities to all nodes
@ -182,20 +185,20 @@ void FixEOStableRX::init()
double *uChem = atom->uChem;
double *dpdTheta = atom->dpdTheta;
double tmp;
if(this->restart_reset){
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
temperature_lookup(i,uCond[i]+uMech[i]+uChem[i],dpdTheta[i]);
temperature_lookup(i,uCond[i]+uMech[i]+uChem[i],dpdTheta[i]);
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
energy_lookup(i,dpdTheta[i],tmp);
uCond[i] = tmp / double(2.0);
uMech[i] = tmp / double(2.0);
uChem[i] = double(0.0);
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
energy_lookup(i,dpdTheta[i],tmp);
uCond[i] = tmp / 2.0;
uMech[i] = tmp / 2.0;
uChem[i] = 0.0;
}
}
}
@ -215,8 +218,8 @@ void FixEOStableRX::post_integrate()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
temperature_lookup(i,uCond[i]+uMech[i]+uChem[i],dpdTheta[i]);
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
}
}
@ -241,8 +244,8 @@ void FixEOStableRX::end_of_step()
if (mask[i] & groupbit){
duChem = uCG[i] - uCGnew[i];
uChem[i] += duChem;
uCG[i] = double(0.0);
uCGnew[i] = double(0.0);
uCG[i] = 0.0;
uCGnew[i] = 0.0;
}
// Communicate the updated momenta and velocities to all nodes
@ -251,8 +254,8 @@ void FixEOStableRX::end_of_step()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
temperature_lookup(i,uCond[i]+uMech[i]+uChem[i],dpdTheta[i]);
if(dpdTheta[i] <= double(0.0))
error->one(FLERR,"Internal temperature <= zero");
if(dpdTheta[i] <= 0.0)
error->one(FLERR,"Internal temperature <= zero");
}
}
@ -332,9 +335,9 @@ void FixEOStableRX::read_file(char *file)
for (ispecies = 0; ispecies < nspecies; ispecies++)
if (strcmp(words[0],&atom->dname[ispecies][0]) == 0) break;
if (ispecies == nspecies) continue;
dHf[ispecies] = atof(words[1]);
if (ispecies < nspecies)
dHf[ispecies] = atof(words[1]);
}
delete [] words;
@ -380,7 +383,7 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
// loop until section found with matching keyword
while (1) {
@ -406,13 +409,13 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
memory->create(tb->efile,tb->ninput,"eos:efile");
memory->create(tb2->rfile,tb2->ninput,"eos:rfile");
memory->create(tb2->efile,tb2->ninput,"eos:efile");
for (int ispecies=1;ispecies<nspecies;ispecies++){
Table *tbl = &tables[ispecies];
Table *tbl2 = &tables2[ispecies];
tbl->ninput = tb->ninput;
tbl2->ninput = tb2->ninput;
memory->create(tbl->rfile,tbl->ninput,"eos:rfile");
memory->create(tbl->efile,tbl->ninput,"eos:efile");
memory->create(tbl2->rfile,tbl2->ninput,"eos:rfile");
@ -421,7 +424,6 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
// read r,e table values from file
int itmp;
double rtmp, tmpE;
int nwords;
char * word;
@ -431,7 +433,7 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
fgets(line,MAXLINE,fp);
for (int i = 0; i < ninputs; i++) {
fgets(line,MAXLINE,fp);
nwords = atom->count_words(line);
if(nwords != nspecies+2){
printf("nwords=%d nspecies=%d\n",nwords,nspecies);
@ -439,7 +441,6 @@ void FixEOStableRX::read_table(Table *tb, Table *tb2, char *file, char *keyword)
}
nwords = 0;
word = strtok(line," \t\n\r\f");
itmp = atoi(word);
word = strtok(NULL," \t\n\r\f");
rtmp = atof(word);
@ -523,10 +524,11 @@ void FixEOStableRX::param_extract(Table *tb, char *line)
int ispecies;
ncolumn = 0;
eosSpecies = new int[nspecies];
if (!eosSpecies)
eosSpecies = new int[nspecies];
for (ispecies = 0; ispecies < nspecies; ispecies++)
eosSpecies[ispecies] = -1;
tb->ninput = 0;
char *word = strtok(line," \t\n\r\f");
@ -540,9 +542,9 @@ void FixEOStableRX::param_extract(Table *tb, char *line)
while (word) {
for (ispecies = 0; ispecies < nspecies; ispecies++)
if (strcmp(word,&atom->dname[ispecies][0]) == 0){
eosSpecies[ncolumn] = ispecies;
ncolumn++;
break;
eosSpecies[ncolumn] = ispecies;
ncolumn++;
break;
}
if (ispecies == nspecies){
printf("name=%s not found in species list\n",word);
@ -552,7 +554,7 @@ void FixEOStableRX::param_extract(Table *tb, char *line)
}
for (int icolumn = 0; icolumn < ncolumn; icolumn++)
if(eosSpecies[icolumn]==-1)
if(eosSpecies[icolumn]==-1)
error->one(FLERR,"EOS data is missing from fix eos/table/rx tabe");
if(ncolumn != nspecies){
printf("ncolumns=%d nspecies=%d\n",ncolumn,nspecies);
@ -648,8 +650,8 @@ void FixEOStableRX::energy_lookup(int id, double thetai, double &ui)
int itable;
double fraction, uTmp, nTotal;
ui = double(0.0);
nTotal = double(0.0);
ui = 0.0;
nTotal = 0.0;
for(int ispecies=0;ispecies<nspecies;ispecies++){
Table *tb = &tables[ispecies];
thetai = MAX(thetai,tb->lo);
@ -659,7 +661,7 @@ void FixEOStableRX::energy_lookup(int id, double thetai, double &ui)
itable = static_cast<int> ((thetai - tb->lo) * tb->invdelta);
fraction = (thetai - tb->r[itable]) * tb->invdelta;
uTmp = tb->e[itable] + fraction*tb->de[itable];
uTmp += dHf[ispecies];
// mol fraction form:
ui += atom->dvector[ispecies][id]*uTmp;
@ -691,11 +693,11 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai)
f1 = u1 - ui;
// Compute guess of t2
t2 = (double(1.0) + double(0.001))*t1;
t2 = (1.0 + 0.001)*t1;
// Compute u2 at t2
energy_lookup(id,t2,u2);
// Compute f1
f2 = u2 - ui;

6
src/USER-DPD/fix_eos_table_rx.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -66,7 +66,7 @@ class FixEOStableRX : public Fix {
int nspecies;
void read_file(char *);
double *dHf;
int pack_reverse_comm(int, int, double *);
@ -116,7 +116,7 @@ Self-explanatory.
E: Incorrect format in eos table/rx file
Self-explanatory.
Self-explanatory.
E: Cannot open file %s

207
src/USER-DPD/fix_rx.cpp
View File

@ -50,6 +50,8 @@ FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
params = NULL;
mol2param = NULL;
pairDPDE = NULL;
id_fix_species = NULL;
id_fix_species_old = NULL;
// Keep track of the argument list.
int iarg = 3;
@ -110,11 +112,19 @@ FixRX::~FixRX()
delete [] stoichReactants[ii];
delete [] stoichProducts[ii];
}
delete [] Arr;
delete [] nArr;
delete [] Ea;
delete [] tempExp;
delete [] stoich;
delete [] stoichReactants;
delete [] stoichProducts;
delete [] kR;
}
delete [] id_fix_species;
delete [] id_fix_species_old;
}
/* ---------------------------------------------------------------------- */
@ -125,10 +135,10 @@ void FixRX::post_constructor()
bool match;
for (int i = 0; i < modify->nfix; i++)
if (strncmp(modify->fix[i]->style,"property/atom",13) == 0)
error->all(FLERR,"fix rx cannot be combined with fix property/atom");
if (strncmp(modify->fix[i]->style,"property/atom",13) == 0)
error->all(FLERR,"fix rx cannot be combined with fix property/atom");
char **tmpspecies = new char*[maxspecies];
char **tmpspecies = new char*[maxspecies];
for(int jj=0; jj < maxspecies; jj++)
tmpspecies[jj] = NULL;
@ -137,7 +147,7 @@ void FixRX::post_constructor()
FILE *fp;
fp = NULL;
if (comm->me == 0) {
fp = fopen(kineticsFile,"r");
fp = force->open_potential(kineticsFile);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open rx file %s",kineticsFile);
@ -147,7 +157,7 @@ void FixRX::post_constructor()
// Assign species names to tmpspecies array and determine the number of unique species
int n,nwords,ispecies;
int n,nwords;
char line[MAXLINE],*ptr;
int eof = 0;
char * word;
@ -179,17 +189,17 @@ void FixRX::post_constructor()
word = strtok(NULL, " \t\n\r\f");
match=false;
for(int jj=0;jj<nUniqueSpecies;jj++){
if(strcmp(word,tmpspecies[jj])==0){
match=true;
break;
}
if(strcmp(word,tmpspecies[jj])==0){
match=true;
break;
}
}
if(!match){
if(nUniqueSpecies+1>=maxspecies)
error->all(FLERR,"Exceeded the maximum number of species permitted in fix rx.");
tmpspecies[nUniqueSpecies] = new char[strlen(word)];
strcpy(tmpspecies[nUniqueSpecies],word);
nUniqueSpecies++;
if(!match){
if(nUniqueSpecies+1>=maxspecies)
error->all(FLERR,"Exceeded the maximum number of species permitted in fix rx.");
tmpspecies[nUniqueSpecies] = new char[strlen(word)+1];
strcpy(tmpspecies[nUniqueSpecies],word);
nUniqueSpecies++;
}
word = strtok(NULL, " \t\n\r\f");
if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0) break;
@ -202,9 +212,6 @@ void FixRX::post_constructor()
// new id = fix-ID + FIX_STORE_ATTRIBUTE
// new fix group = group for this fix
id_fix_species = NULL;
id_fix_species_old = NULL;
n = strlen(id) + strlen("_SPECIES") + 1;
id_fix_species = new char[n];
n = strlen(id) + strlen("_SPECIES_OLD") + 1;
@ -231,8 +238,8 @@ void FixRX::post_constructor()
strncat(str1,tmpspecies[ii],strlen(tmpspecies[ii]));
strncat(str2,tmpspecies[ii],strlen(tmpspecies[ii]));
strncat(str2,"Old",3);
newarg[ii+3] = new char[strlen(str1)];
newarg2[ii+3] = new char[strlen(str2)];
newarg[ii+3] = new char[strlen(str1)+1];
newarg2[ii+3] = new char[strlen(str2)+1];
strcpy(newarg[ii+3],str1);
strcpy(newarg2[ii+3],str2);
}
@ -249,11 +256,15 @@ void FixRX::post_constructor()
if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
for(int jj=0;jj<maxspecies;jj++) delete tmpspecies[jj];
delete tmpspecies;
for(int jj=0;jj<nspecies;jj++) {
delete[] tmpspecies[jj];
delete[] newarg[jj+3];
delete[] newarg2[jj+3];
}
delete newarg;
delete newarg2;
delete[] newarg;
delete[] newarg2;
delete[] tmpspecies;
read_file( kineticsFile );
@ -282,7 +293,7 @@ void FixRX::init()
bool eos_flag = false;
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"eos/table/rx") == 0) eos_flag = true;
if(!eos_flag) error->all(FLERR,"fix rx requires fix eos/table/rx to be specified");
if(!eos_flag) error->all(FLERR,"fix rx requires fix eos/table/rx to be specified");
}
/* ---------------------------------------------------------------------- */
@ -292,20 +303,19 @@ void FixRX::setup_pre_force(int vflag)
int nlocal = atom->nlocal;
int nghost = atom->nghost;
int *mask = atom->mask;
double *dpdTheta = atom->dpdTheta;
int newton_pair = force->newton_pair;
double tmp, theta;
double tmp;
int ii;
if(localTempFlag){
if (newton_pair) {
dpdThetaLocal = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
dpdThetaLocal[ii] = double(0.0);
dpdThetaLocal[ii] = 0.0;
} else {
dpdThetaLocal = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
dpdThetaLocal[ii] = double(0.0);
dpdThetaLocal[ii] = 0.0;
}
computeLocalTemperature();
}
@ -317,12 +327,10 @@ void FixRX::setup_pre_force(int vflag)
}
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
if(localTempFlag) theta = dpdThetaLocal[i];
else theta=dpdTheta[i];
// Set the reaction rate constants to zero: no reactions occur at step 0
for(int irxn=0;irxn<nreactions;irxn++)
kR[irxn] = double(0.0);
kR[irxn] = 0.0;
if(odeIntegrationFlag==ODE_LAMMPS_RK4) rk4(i);
}
@ -347,14 +355,13 @@ void FixRX::pre_force(int vflag)
if (newton_pair) {
dpdThetaLocal = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
dpdThetaLocal[ii] = double(0.0);
dpdThetaLocal[ii] = 0.0;
} else {
dpdThetaLocal = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
dpdThetaLocal[ii] = double(0.0);
dpdThetaLocal[ii] = 0.0;
}
computeLocalTemperature();
double *dpdTheta = this->dpdThetaLocal;
}
for (int i = 0; i < nlocal; i++)
@ -364,7 +371,7 @@ void FixRX::pre_force(int vflag)
//Compute the reaction rate constants
for(int irxn=0;irxn<nreactions;irxn++)
kR[irxn] = Arr[irxn]*pow(theta,nArr[irxn])*exp(-Ea[irxn]/force->boltz/theta);
kR[irxn] = Arr[irxn]*pow(theta,nArr[irxn])*exp(-Ea[irxn]/force->boltz/theta);
if(odeIntegrationFlag==ODE_LAMMPS_RK4) rk4(i);
}
@ -384,7 +391,7 @@ void FixRX::read_file(char *file)
FILE *fp;
fp = NULL;
if (comm->me == 0) {
fp = fopen(file,"r");
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open rx file %s",file);
@ -421,7 +428,7 @@ void FixRX::read_file(char *file)
}
// open file on proc 0
if (comm->me == 0) fp = fopen(file,"r");
if (comm->me == 0) fp = force->open_potential(file);
// read each reaction from kinetics file
eof=0;
@ -444,14 +451,14 @@ void FixRX::read_file(char *file)
kR = new double[nreactions];
for (int ii=0;ii<nreactions;ii++){
for (int jj=0;jj<nspecies;jj++){
stoich[ii][jj] = double(0.0);
stoichReactants[ii][jj] = double(0.0);
stoichProducts[ii][jj] = double(0.0);
stoich[ii][jj] = 0.0;
stoichReactants[ii][jj] = 0.0;
stoichProducts[ii][jj] = 0.0;
}
}
nreactions=0;
sign = double(-1.0);
sign = -1.0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
@ -479,31 +486,37 @@ void FixRX::read_file(char *file)
tmpStoich = atof(word);
word = strtok(NULL, " \t\n\r\f");
for (ispecies = 0; ispecies < nspecies; ispecies++){
if (strcmp(word,&atom->dname[ispecies][0]) == 0){
stoich[nreactions][ispecies] += sign*tmpStoich;
if(sign<double(0.0)) stoichReactants[nreactions][ispecies] += tmpStoich;
else stoichProducts[nreactions][ispecies] += tmpStoich;
break;
}
if (strcmp(word,&atom->dname[ispecies][0]) == 0){
stoich[nreactions][ispecies] += sign*tmpStoich;
if(sign<0.0)
stoichReactants[nreactions][ispecies] += tmpStoich;
else stoichProducts[nreactions][ispecies] += tmpStoich;
break;
}
}
if(ispecies==nspecies){
printf("%s mol fraction is not found in data file\n",word);
printf("nspecies=%d ispecies=%d\n",nspecies,ispecies);
error->all(FLERR,"Illegal fix rx command");
if (comm->me) {
fprintf(stderr,"%s mol fraction is not found in data file\n",word);
fprintf(stderr,"nspecies=%d ispecies=%d\n",nspecies,ispecies);
}
error->all(FLERR,"Illegal fix rx command");
}
word = strtok(NULL, " \t\n\r\f");
if(strcmp(word,"=") == 0) sign = double(1.0);
if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0){
if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
Arr[nreactions] = atof(word);
word = strtok(NULL, " \t\n\r\f");
if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
nArr[nreactions] = atof(word);
word = strtok(NULL, " \t\n\r\f");
if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
Ea[nreactions] = atof(word);
sign = double(-1.0);
break;
if(word==NULL)
error->all(FLERR,"Missing parameters in reaction kinetic equation");
if(strcmp(word,"=") == 0) sign = 1.0;
if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0) {
Arr[nreactions] = atof(word);
word = strtok(NULL, " \t\n\r\f");
if(word==NULL)
error->all(FLERR,"Missing parameters in reaction kinetic equation");
nArr[nreactions] = atof(word);
word = strtok(NULL, " \t\n\r\f");
if(word==NULL)
error->all(FLERR,"Missing parameters in reaction kinetic equation");
Ea[nreactions] = atof(word);
sign = -1.0;
break;
}
word = strtok(NULL, " \t\n\r\f");
}
@ -527,8 +540,8 @@ void FixRX::setupParams()
n = -1;
for (j = 0; j < nreactions; j++) {
if (i == params[j].ispecies) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = j;
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = j;
}
}
mol2param[i] = n;
@ -568,13 +581,13 @@ void FixRX::rk4(int id)
// k2
for (int ispecies = 0; ispecies < nspecies; ispecies++)
yp[ispecies] = y[ispecies] + double(0.5)*h*k1[ispecies];
yp[ispecies] = y[ispecies] + 0.5*h*k1[ispecies];
rhs(0.0,yp,k2,dummyArray);
// k3
for (int ispecies = 0; ispecies < nspecies; ispecies++)
yp[ispecies] = y[ispecies] + double(0.5)*h*k2[ispecies];
yp[ispecies] = y[ispecies] + 0.5*h*k2[ispecies];
rhs(0.0,yp,k3,dummyArray);
@ -585,17 +598,17 @@ void FixRX::rk4(int id)
rhs(0.0,yp,k4,dummyArray);
for (int ispecies = 0; ispecies < nspecies; ispecies++)
y[ispecies] += h*(k1[ispecies]/6.0 + k2[ispecies]/3.0 +
k3[ispecies]/3.0 + k4[ispecies]/6.0);
y[ispecies] += h*(k1[ispecies]/6.0 + k2[ispecies]/3.0 +
k3[ispecies]/3.0 + k4[ispecies]/6.0);
} // end for (int step...
// Store the solution back in atom->dvector.
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if(y[ispecies] < double(-1.0e-10))
if(y[ispecies] < -1.0e-10)
error->one(FLERR,"Computed concentration in RK4 solver is < -1.0e-10");
else if(y[ispecies] < double(0.0))
y[ispecies] = double(0.0);
else if(y[ispecies] < 0.0)
y[ispecies] = 0.0;
atom->dvector[ispecies][id] = y[ispecies];
}
delete [] k1;
@ -612,7 +625,7 @@ int FixRX::rhs(double t, const double *y, double *dydt, void *params)
int nspecies = atom->nspecies_dpd;
for(int ispecies=0; ispecies<nspecies; ispecies++)
dydt[ispecies] = double(0.0);
dydt[ispecies] = 0.0;
// Construct the reaction rate laws
for(int jrxn=0; jrxn<nreactions; jrxn++){
@ -624,7 +637,7 @@ int FixRX::rhs(double t, const double *y, double *dydt, void *params)
}
rxnRateLaw[jrxn] = rxnRateLawForward;
}
// Construct the reaction rates for each species
for(int ispecies=0; ispecies<nspecies; ispecies++)
for(int jrxn=0; jrxn<nreactions; jrxn++)
@ -649,18 +662,18 @@ void FixRX::computeLocalTemperature()
int newton_pair = force->newton_pair;
// local temperature variables
double wij, fr, fr4;
double wij = 0.0;
double *dpdTheta = atom->dpdTheta;
// Initialize the local density and local temperature arrays
if (newton_pair) {
sumWeights = new double[nlocal+nghost];
for (ii = 0; ii < nlocal+nghost; ii++)
sumWeights[ii] = double(0.0);
sumWeights[ii] = 0.0;
} else {
sumWeights = new double[nlocal];
for (ii = 0; ii < nlocal; ii++)
dpdThetaLocal[ii] = double(0.0);
dpdThetaLocal[ii] = 0.0;
}
inum = pairDPDE->list->inum;
@ -690,24 +703,22 @@ void FixRX::computeLocalTemperature()
if (rsq < pairDPDE->cutsq[itype][jtype]) {
double rcut = sqrt(pairDPDE->cutsq[itype][jtype]);
double rij = sqrt(rsq);
double tmpFactor = double(1.0)-rij/rcut;
double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
double ratio = rij/rcut;
double rij = sqrt(rsq);
double ratio = rij/rcut;
// Lucy's Weight Function
if(wtFlag==LUCY){
wij = (double(1.0)+double(3.0)*ratio) * (double(1.0)-ratio)*(double(1.0)-ratio)*(double(1.0)-ratio);
dpdThetaLocal[i] += wij/dpdTheta[j];
if (newton_pair || j < nlocal)
dpdThetaLocal[j] += wij/dpdTheta[i];
}
// Lucy's Weight Function
if(wtFlag==LUCY){
wij = (1.0+3.0*ratio) * (1.0-ratio)*(1.0-ratio)*(1.0-ratio);
dpdThetaLocal[i] += wij/dpdTheta[j];
if (newton_pair || j < nlocal)
dpdThetaLocal[j] += wij/dpdTheta[i];
}
sumWeights[i] += wij;
sumWeights[i] += wij;
if (newton_pair || j < nlocal) {
sumWeights[j] += wij;
}
sumWeights[j] += wij;
}
}
}
}
@ -718,17 +729,17 @@ void FixRX::computeLocalTemperature()
// Lucy Weight Function
if(wtFlag==LUCY){
wij = double(1.0);
wij = 1.0;
dpdThetaLocal[i] += wij / dpdTheta[i];
}
sumWeights[i] += wij;
// Normalized local temperature
dpdThetaLocal[i] = dpdThetaLocal[i] / sumWeights[i];
if(localTempFlag == HARMONIC)
dpdThetaLocal[i] = double(1.0) / dpdThetaLocal[i];
dpdThetaLocal[i] = 1.0 / dpdThetaLocal[i];
}
delete [] sumWeights;

4
src/USER-DPD/fix_rx.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -129,6 +129,6 @@ Self-explanatory
E: Computed concentration in RK4 solver is < -1.0e-10.
Self-explanatory: Adjust settings for the RK4 solver.
Self-explanatory: Adjust settings for the RK4 solver.
*/

499
src/USER-DPD/fix_shardlow.cpp
View File

@ -12,11 +12,11 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Contributing authors:
James Larentzos (U.S. Army Research Laboratory)
and Timothy I. Mattox (Engility Corporation)
Martin Lisal (Institute of Chemical Process Fundamentals
Martin Lisal (Institute of Chemical Process Fundamentals
of the Czech Academy of Sciences and J. E. Purkinje University)
John Brennan, Joshua Moore and William Mattson (Army Research Lab)
@ -24,7 +24,7 @@
Please cite the related publications:
J. P. Larentzos, J. K. Brennan, J. D. Moore, M. Lisal, W. D. Mattson,
"Parallel implementation of isothermal and isoenergetic Dissipative
Particle Dynamics using Shardlow-like splitting algorithms",
Particle Dynamics using Shardlow-like splitting algorithms",
Computer Physics Communications, 2014, 185, pp 1987--1998.
M. Lisal, J. K. Brennan, J. Bonet Avalos, "Dissipative particle dynamics
@ -96,15 +96,17 @@ FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) :
pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy",1);
if(pairDPDE){
comm_forward = 10;
comm_forward = 3;
comm_reverse = 5;
} else {
comm_forward = 6;
comm_forward = 3;
comm_reverse = 3;
}
if(pairDPD == NULL && pairDPDE == NULL)
error->all(FLERR,"Must use pair_style dpd/fdt or dpd/fdt/energy with fix shardlow");
if(!(atom->dpd_flag))
error->all(FLERR,"Must use atom_style dpd with fix shardlow");
}
/* ---------------------------------------------------------------------- */
@ -113,7 +115,6 @@ int FixShardlow::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= PRE_NEIGHBOR;
return mask;
}
@ -144,13 +145,6 @@ void FixShardlow::setup(int vflag)
}
}
/* ---------------------------------------------------------------------- */
void FixShardlow::setup_pre_force(int vflag)
{
neighbor->build_one(list);
}
/* ----------------------------------------------------------------------
Perform the stochastic integration and Shardlow update
Allow for both per-type and per-atom mass
@ -176,14 +170,12 @@ void FixShardlow::initial_integrate(int vflag)
int *type = atom->type;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
int nall = nlocal + nghost;
int newton_pair = force->newton_pair;
double randPair;
int *ssaAIR = atom->ssaAIR;
double *uCond = atom->uCond;
double *uMech = atom->uMech;
double *duCond = atom->duCond;
double *duMech = atom->duMech;
double *dpdTheta = atom->dpdTheta;
double kappa_ij, alpha_ij, theta_ij, gamma_ij, sigma_ij;
double vxi, vyi, vzi, vxj, vyj, vzj;
@ -203,7 +195,7 @@ void FixShardlow::initial_integrate(int vflag)
double bby = domain->subhi[1] - domain->sublo[1];
double bbz = domain->subhi[2] - domain->sublo[2];
double rcut = double(2.0)*neighbor->cutneighmax;
double rcut = 2.0*neighbor->cutneighmax;
if (domain->triclinic)
error->all(FLERR,"Fix shardlow does not yet support triclinic geometries");
@ -211,25 +203,8 @@ void FixShardlow::initial_integrate(int vflag)
if(rcut >= bbx || rcut >= bby || rcut>= bbz )
error->all(FLERR,"Shardlow algorithm requires sub-domain length > 2*(rcut+skin). Either reduce the number of processors requested, or change the cutoff/skin\n");
// Allocate memory for the dvSSA arrays
dvSSA = new double*[nall];
for (ii = 0; ii < nall; ii++) {
dvSSA[ii] = new double[3];
}
// Zero the momenta
for (ii = 0; ii < nlocal; ii++) {
dvSSA[ii][0] = double(0.0);
dvSSA[ii][1] = double(0.0);
dvSSA[ii][2] = double(0.0);
if(pairDPDE){
duCond[ii] = double(0.0);
duMech[ii] = double(0.0);
}
}
// Communicate the updated momenta and velocities to all nodes
comm->forward_comm_fix(this);
// Allocate memory for v_t0 to hold the initial velocities for the ghosts
v_t0 = (double (*)[3]) memory->smalloc(sizeof(double)*3*nghost, "FixShardlow:v_t0");
// Define pointers to access the neighbor list
if(pairDPDE){
@ -244,9 +219,20 @@ void FixShardlow::initial_integrate(int vflag)
firstneigh = pairDPD->list->firstneigh;
}
//Loop over all 14 directions (8 stages)
//Loop over all 14 directions (8 stages)
for (int idir = 1; idir <=8; idir++){
if (idir > 1) {
// Communicate the updated velocities to all nodes
comm->forward_comm_fix(this);
if(pairDPDE){
// Zero out the ghosts' uCond & uMech to be used as delta accumulators
memset(&(uCond[nlocal]), 0, sizeof(double)*nghost);
memset(&(uMech[nlocal]), 0, sizeof(double)*nghost);
}
}
// Loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
@ -254,221 +240,190 @@ void FixShardlow::initial_integrate(int vflag)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
// load velocity for i from memory
vxi = v[i][0];
vyi = v[i][1];
vzi = v[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
// Loop over Directional Neighbors only
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
if (neighbor->ssa_airnum[j] != idir) continue;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
j = jlist[jj];
j &= NEIGHMASK;
if (ssaAIR[j] != idir) continue;
jtype = type[j];
if(pairDPDE){
cut2 = pairDPDE->cutsq[itype][jtype];
cut = pairDPDE->cut[itype][jtype];
} else {
cut2 = pairDPD->cutsq[itype][jtype];
cut = pairDPD->cut[itype][jtype];
}
// if (rsq < pairDPD->cutsq[itype][jtype]) {
if (rsq < cut2) {
r = sqrt(rsq);
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = double(1.0)/r;
// Store the velocities from previous Shardlow step
vx0i = v[i][0] + dvSSA[i][0];
vy0i = v[i][1] + dvSSA[i][1];
vz0i = v[i][2] + dvSSA[i][2];
vx0j = v[j][0] + dvSSA[j][0];
vy0j = v[j][1] + dvSSA[j][1];
vz0j = v[j][2] + dvSSA[j][2];
// Compute the velocity difference between atom i and atom j
delvx = vx0i - vx0j;
delvy = vy0i - vy0j;
delvz = vz0i - vz0j;
dot = (delx*delvx + dely*delvy + delz*delvz);
// wr = double(1.0) - r/pairDPD->cut[itype][jtype];
wr = double(1.0) - r/cut;
wd = wr*wr;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if(pairDPDE){
// Compute the current temperature
theta_ij = double(0.5)*(double(1.0)/dpdTheta[i] + double(1.0)/dpdTheta[j]);
theta_ij = double(1.0)/theta_ij;
sigma_ij = pairDPDE->sigma[itype][jtype];
randnum = pairDPDE->random->gaussian();
} else {
theta_ij = pairDPD->temperature;
sigma_ij = pairDPD->sigma[itype][jtype];
randnum = pairDPD->random->gaussian();
}
if(pairDPDE){
cut2 = pairDPDE->cutsq[itype][jtype];
cut = pairDPDE->cut[itype][jtype];
} else {
cut2 = pairDPD->cutsq[itype][jtype];
cut = pairDPD->cut[itype][jtype];
}
gamma_ij = sigma_ij*sigma_ij / (2.0*force->boltz*theta_ij);
randPair = sigma_ij*wr*randnum*dtsqrt;
factor_dpd = -dt*gamma_ij*wd*dot*rinv;
factor_dpd += randPair;
factor_dpd *= double(0.5);
// Compute momentum change between t and t+dt
dpx = factor_dpd*delx*rinv;
dpy = factor_dpd*dely*rinv;
dpz = factor_dpd*delz*rinv;
if (rmass) {
mass_i = rmass[i];
mass_j = rmass[j];
} else {
mass_i = mass[itype];
mass_j = mass[jtype];
}
massinv_i = double(1.0) / mass_i;
massinv_j = double(1.0) / mass_j;
// Update the delta velocity on i
dvSSA[i][0] += dpx*force->ftm2v*massinv_i;
dvSSA[i][1] += dpy*force->ftm2v*massinv_i;
dvSSA[i][2] += dpz*force->ftm2v*massinv_i;
if (newton_pair || j < nlocal) {
// Update the delta velocity on j
dvSSA[j][0] -= dpx*force->ftm2v*massinv_j;
dvSSA[j][1] -= dpy*force->ftm2v*massinv_j;
dvSSA[j][2] -= dpz*force->ftm2v*massinv_j;
}
//ii. Compute the velocity diff
delvx = v[i][0] + dvSSA[i][0] - v[j][0] - dvSSA[j][0];
delvy = v[i][1] + dvSSA[i][1] - v[j][1] - dvSSA[j][1];
delvz = v[i][2] + dvSSA[i][2] - v[j][2] - dvSSA[j][2];
dot = delx*delvx + dely*delvy + delz*delvz;
//iii. Compute dpi again
mu_ij = massinv_i + massinv_j;
denom = double(1.0) + double(0.5)*mu_ij*gamma_ij*wd*dt*force->ftm2v;
factor_dpd = -double(0.5)*dt*gamma_ij*wd*force->ftm2v/denom;
factor_dpd1 = factor_dpd*(mu_ij*randPair);
factor_dpd1 += randPair;
factor_dpd1 *= double(0.5);
// Compute the momentum change between t and t+dt
dpx = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*delx*rinv;
dpy = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*dely*rinv;
dpz = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*delz*rinv;
//Update the velocity change on i
dvSSA[i][0] += dpx*force->ftm2v*massinv_i;
dvSSA[i][1] += dpy*force->ftm2v*massinv_i;
dvSSA[i][2] += dpz*force->ftm2v*massinv_i;
if (newton_pair || j < nlocal) {
//Update the velocity change on j
dvSSA[j][0] -= dpx*force->ftm2v*massinv_j;
dvSSA[j][1] -= dpy*force->ftm2v*massinv_j;
dvSSA[j][2] -= dpz*force->ftm2v*massinv_j;
}
// if (rsq < pairDPD->cutsq[itype][jtype])
if (rsq < cut2) {
r = sqrt(rsq);
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
if(pairDPDE){
// Compute uCond
randnum = pairDPDE->random->gaussian();
kappa_ij = pairDPDE->kappa[itype][jtype];
alpha_ij = sqrt(2.0*force->boltz*kappa_ij);
randPair = alpha_ij*wr*randnum*dtsqrt;
factor_dpd = kappa_ij*(double(1.0)/dpdTheta[i] - double(1.0)/dpdTheta[j])*wd*dt;
factor_dpd += randPair;
duCond[i] += factor_dpd;
if (newton_pair || j < nlocal) {
duCond[j] -= factor_dpd;
}
// Compute uMech
vxi = v[i][0] + dvSSA[i][0];
vyi = v[i][1] + dvSSA[i][1];
vzi = v[i][2] + dvSSA[i][2];
vxj = v[j][0] + dvSSA[j][0];
vyj = v[j][1] + dvSSA[j][1];
vzj = v[j][2] + dvSSA[j][2];
dot1 = vxi*vxi + vyi*vyi + vzi*vzi;
dot2 = vxj*vxj + vyj*vyj + vzj*vzj;
dot3 = vx0i*vx0i + vy0i*vy0i + vz0i*vz0i;
dot4 = vx0j*vx0j + vy0j*vy0j + vz0j*vz0j;
dot1 = dot1*mass_i;
dot2 = dot2*mass_j;
dot3 = dot3*mass_i;
dot4 = dot4*mass_j;
factor_dpd = double(0.25)*(dot1+dot2-dot3-dot4)/force->ftm2v;
duMech[i] -= factor_dpd;
if (newton_pair || j < nlocal) {
duMech[j] -= factor_dpd;
}
}
}
// Keep a copy of the velocities from previous Shardlow step
vx0i = vxi;
vy0i = vyi;
vz0i = vzi;
vx0j = vxj = v[j][0];
vy0j = vyj = v[j][1];
vz0j = vzj = v[j][2];
// Compute the velocity difference between atom i and atom j
delvx = vx0i - vx0j;
delvy = vy0i - vy0j;
delvz = vz0i - vz0j;
dot = (delx*delvx + dely*delvy + delz*delvz);
// wr = 1.0 - r/pairDPD->cut[itype][jtype];
wr = 1.0 - r/cut;
wd = wr*wr;
if(pairDPDE){
// Compute the current temperature
theta_ij = 0.5*(1.0/dpdTheta[i] + 1.0/dpdTheta[j]);
theta_ij = 1.0/theta_ij;
sigma_ij = pairDPDE->sigma[itype][jtype];
randnum = pairDPDE->random->gaussian();
} else {
theta_ij = pairDPD->temperature;
sigma_ij = pairDPD->sigma[itype][jtype];
randnum = pairDPD->random->gaussian();
}
gamma_ij = sigma_ij*sigma_ij / (2.0*force->boltz*theta_ij);
randPair = sigma_ij*wr*randnum*dtsqrt;
factor_dpd = -dt*gamma_ij*wd*dot*rinv;
factor_dpd += randPair;
factor_dpd *= 0.5;
// Compute momentum change between t and t+dt
dpx = factor_dpd*delx*rinv;
dpy = factor_dpd*dely*rinv;
dpz = factor_dpd*delz*rinv;
if (rmass) {
mass_i = rmass[i];
mass_j = rmass[j];
} else {
mass_i = mass[itype];
mass_j = mass[jtype];
}
massinv_i = 1.0 / mass_i;
massinv_j = 1.0 / mass_j;
// Update the velocity on i
vxi += dpx*force->ftm2v*massinv_i;
vyi += dpy*force->ftm2v*massinv_i;
vzi += dpz*force->ftm2v*massinv_i;
if (newton_pair || j < nlocal) {
// Update the velocity on j
vxj -= dpx*force->ftm2v*massinv_j;
vyj -= dpy*force->ftm2v*massinv_j;
vzj -= dpz*force->ftm2v*massinv_j;
}
//ii. Compute the velocity diff
delvx = vxi - vxj;
delvy = vyi - vyj;
delvz = vzi - vzj;
dot = delx*delvx + dely*delvy + delz*delvz;
//iii. Compute dpi again
mu_ij = massinv_i + massinv_j;
denom = 1.0 + 0.5*mu_ij*gamma_ij*wd*dt*force->ftm2v;
factor_dpd = -0.5*dt*gamma_ij*wd*force->ftm2v/denom;
factor_dpd1 = factor_dpd*(mu_ij*randPair);
factor_dpd1 += randPair;
factor_dpd1 *= 0.5;
// Compute the momentum change between t and t+dt
dpx = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*delx*rinv;
dpy = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*dely*rinv;
dpz = (factor_dpd*dot*rinv/force->ftm2v + factor_dpd1)*delz*rinv;
// Update the velocity on i
vxi += dpx*force->ftm2v*massinv_i;
vyi += dpy*force->ftm2v*massinv_i;
vzi += dpz*force->ftm2v*massinv_i;
if (newton_pair || j < nlocal) {
// Update the velocity on j
vxj -= dpx*force->ftm2v*massinv_j;
vyj -= dpy*force->ftm2v*massinv_j;
vzj -= dpz*force->ftm2v*massinv_j;
// Store updated velocity for j
v[j][0] = vxj;
v[j][1] = vyj;
v[j][2] = vzj;
}
if(pairDPDE){
// Compute uCond
randnum = pairDPDE->random->gaussian();
kappa_ij = pairDPDE->kappa[itype][jtype];
alpha_ij = sqrt(2.0*force->boltz*kappa_ij);
randPair = alpha_ij*wr*randnum*dtsqrt;
factor_dpd = kappa_ij*(1.0/dpdTheta[i] - 1.0/dpdTheta[j])*wd*dt;
factor_dpd += randPair;
uCond[i] += factor_dpd;
if (newton_pair || j < nlocal) {
uCond[j] -= factor_dpd;
}
// Compute uMech
dot1 = vxi*vxi + vyi*vyi + vzi*vzi;
dot2 = vxj*vxj + vyj*vyj + vzj*vzj;
dot3 = vx0i*vx0i + vy0i*vy0i + vz0i*vz0i;
dot4 = vx0j*vx0j + vy0j*vy0j + vz0j*vz0j;
dot1 = dot1*mass_i;
dot2 = dot2*mass_j;
dot3 = dot3*mass_i;
dot4 = dot4*mass_j;
factor_dpd = 0.25*(dot1+dot2-dot3-dot4)/force->ftm2v;
uMech[i] -= factor_dpd;
if (newton_pair || j < nlocal) {
uMech[j] -= factor_dpd;
}
}
}
}
// store updated velocity for i
v[i][0] = vxi;
v[i][1] = vyi;
v[i][2] = vzi;
}
// Communicate the ghost delta velocities to the locally owned atoms
comm->reverse_comm_fix(this);
// Zero the dv
for (ii = 0; ii < nlocal; ii++) {
// Shardlow update
v[ii][0] += dvSSA[ii][0];
v[ii][1] += dvSSA[ii][1];
v[ii][2] += dvSSA[ii][2];
dvSSA[ii][0] = double(0.0);
dvSSA[ii][1] = double(0.0);
dvSSA[ii][2] = double(0.0);
if(pairDPDE){
uCond[ii] += duCond[ii];
uMech[ii] += duMech[ii];
duCond[ii] = double(0.0);
duMech[ii] = double(0.0);
}
}
// Communicate the updated momenta and velocities to all nodes
comm->forward_comm_fix(this);
// Communicate the ghost deltas to the atom owners
if (idir > 1) comm->reverse_comm_fix(this);
} //End Loop over all directions For idir = Top, Top-Right, Right, Bottom-Right, Back
for (ii = 0; ii < nall; ii++) {
delete dvSSA[ii];
}
delete [] dvSSA;
}
/* ----------------------------------------------------------------------
* assign owned and ghost atoms their ssa active interaction region numbers
------------------------------------------------------------------------- */
void FixShardlow::setup_pre_neighbor()
{
neighbor->assign_ssa_airnums();
}
void FixShardlow::pre_neighbor()
{
neighbor->assign_ssa_airnums();
memory->sfree(v_t0);
v_t0 = NULL;
}
/* ---------------------------------------------------------------------- */
@ -477,26 +432,13 @@ int FixShardlow::pack_forward_comm(int n, int *list, double *buf, int pbc_flag,
{
int ii,jj,m;
double **v = atom->v;
double *duCond = atom->duCond;
double *duMech = atom->duMech;
double *uCond = atom->uCond;
double *uMech = atom->uMech;
m = 0;
for (ii = 0; ii < n; ii++) {
jj = list[ii];
buf[m++] = dvSSA[jj][0];
buf[m++] = dvSSA[jj][1];
buf[m++] = dvSSA[jj][2];
buf[m++] = v[jj][0];
buf[m++] = v[jj][1];
buf[m++] = v[jj][2];
if(pairDPDE){
buf[m++] = duCond[jj];
buf[m++] = duMech[jj];
buf[m++] = uCond[jj];
buf[m++] = uMech[jj];
}
}
return m;
}
@ -506,27 +448,15 @@ int FixShardlow::pack_forward_comm(int n, int *list, double *buf, int pbc_flag,
void FixShardlow::unpack_forward_comm(int n, int first, double *buf)
{
int ii,m,last;
int nlocal = atom->nlocal;
double **v = atom->v;
double *duCond = atom->duCond;
double *duMech = atom->duMech;
double *uCond = atom->uCond;
double *uMech = atom->uMech;
m = 0;
last = first + n ;
for (ii = first; ii < last; ii++) {
dvSSA[ii][0] = buf[m++];
dvSSA[ii][1] = buf[m++];
dvSSA[ii][2] = buf[m++];
v[ii][0] = buf[m++];
v[ii][1] = buf[m++];
v[ii][2] = buf[m++];
if(pairDPDE){
duCond[ii] = buf[m++];
duMech[ii] = buf[m++];
uCond[ii] = buf[m++];
uMech[ii] = buf[m++];
}
v_t0[ii - nlocal][0] = v[ii][0] = buf[m++];
v_t0[ii - nlocal][1] = v[ii][1] = buf[m++];
v_t0[ii - nlocal][2] = v[ii][2] = buf[m++];
}
}
@ -535,18 +465,20 @@ void FixShardlow::unpack_forward_comm(int n, int first, double *buf)
int FixShardlow::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
double *duCond = atom->duCond;
double *duMech = atom->duMech;
int nlocal = atom->nlocal;
double **v = atom->v;
double *uCond = atom->uCond;
double *uMech = atom->uMech;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = dvSSA[i][0];
buf[m++] = dvSSA[i][1];
buf[m++] = dvSSA[i][2];
buf[m++] = v[i][0] - v_t0[i - nlocal][0];
buf[m++] = v[i][1] - v_t0[i - nlocal][1];
buf[m++] = v[i][2] - v_t0[i - nlocal][2];
if(pairDPDE){
buf[m++] = duCond[i];
buf[m++] = duMech[i];
buf[m++] = uCond[i]; // for ghosts, this is an accumulated delta
buf[m++] = uMech[i]; // for ghosts, this is an accumulated delta
}
}
return m;
@ -557,19 +489,20 @@ int FixShardlow::pack_reverse_comm(int n, int first, double *buf)
void FixShardlow::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
double *duCond = atom->duCond;
double *duMech = atom->duMech;
double **v = atom->v;
double *uCond = atom->uCond;
double *uMech = atom->uMech;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
dvSSA[j][0] += buf[m++];
dvSSA[j][1] += buf[m++];
dvSSA[j][2] += buf[m++];
v[j][0] += buf[m++];
v[j][1] += buf[m++];
v[j][2] += buf[m++];
if(pairDPDE){
duCond[j] += buf[m++];
duMech[j] += buf[m++];
uCond[j] += buf[m++]; // add in the accumulated delta
uMech[j] += buf[m++]; // add in the accumulated delta
}
}
}

16
src/USER-DPD/fix_shardlow.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -31,12 +31,8 @@ class FixShardlow : public Fix {
int setmask();
virtual void init_list(int,class NeighList *);
virtual void setup(int);
virtual void setup_pre_force(int);
virtual void initial_integrate(int);
void setup_pre_neighbor();
void pre_neighbor();
protected:
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
@ -45,7 +41,7 @@ class FixShardlow : public Fix {
class PairDPDfdt *pairDPD;
class PairDPDfdtEnergy *pairDPDE;
double **dvSSA;
double (*v_t0)[3];
private:
class NeighList *list;
@ -71,7 +67,7 @@ Self-explanatory.
E: Must use pair_style dpd/fdt or dpd/fdt/energy with fix shardlow
E: A deterministic integrator must be specified after fix shardlow in input
E: A deterministic integrator must be specified after fix shardlow in input
file (e.g. fix nve or fix nph).
Self-explanatory.
@ -84,11 +80,11 @@ E: Fix shardlow does not yet support triclinic geometries
Self-explanatory.
E: Shardlow algorithm requires sub-domain length > 2*(rcut+skin). Either
E: Shardlow algorithm requires sub-domain length > 2*(rcut+skin). Either
reduce the number of processors requested, or change the cutoff/skin
The Shardlow splitting algorithm requires the size of the sub-domain lengths
to be are larger than twice the cutoff+skin. Generally, the domain decomposition
The Shardlow splitting algorithm requires the size of the sub-domain lengths
to be are larger than twice the cutoff+skin. Generally, the domain decomposition
is dependant on the number of processors requested.
*/

10
src/USER-DPD/pair_dpd_fdt.cpp
View File

@ -114,12 +114,12 @@ void PairDPDfdt::compute(int eflag, int vflag)
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
wr = 1.0 - r/cut[itype][jtype];
wd = wr*wr;
wd = wr*wr;
// conservative force = a0 * wr
fpair = a0[itype][jtype]*wr;
fpair *= factor_dpd*rinv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
@ -179,7 +179,9 @@ void PairDPDfdt::settings(int narg, char **arg)
temperature = force->numeric(FLERR,arg[0]);
cut_global = force->numeric(FLERR,arg[1]);
seed = force->inumeric(FLERR,arg[2]);
if (atom->dpd_flag != 1)
error->all(FLERR,"pair_style dpd/fdt requires atom_style dpd");
// initialize Marsaglia RNG with processor-unique seed
if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
@ -212,7 +214,7 @@ void PairDPDfdt::coeff(int narg, char **arg)
double a0_one = force->numeric(FLERR,arg[2]);
double sigma_one = force->numeric(FLERR,arg[3]);
double cut_one = cut_global;
if (narg == 5) cut_one = force->numeric(FLERR,arg[4]);
int count = 0;

2
src/USER-DPD/pair_dpd_fdt.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov

10
src/USER-DPD/pair_dpd_fdt_energy.cpp
View File

@ -115,12 +115,12 @@ void PairDPDfdtEnergy::compute(int eflag, int vflag)
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
wr = 1.0 - r/cut[itype][jtype];
wd = wr*wr;
wd = wr*wr;
// conservative force = a0 * wr
fpair = a0[itype][jtype]*wr;
fpair *= factor_dpd*rinv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
@ -180,9 +180,9 @@ void PairDPDfdtEnergy::settings(int narg, char **arg)
cut_global = force->numeric(FLERR,arg[0]);
seed = force->inumeric(FLERR,arg[1]);
if (atom->dpd_flag != 1)
if (atom->dpd_flag != 1)
error->all(FLERR,"pair_style dpd/fdt/energy requires atom_style with internal temperature and energies (e.g. dpd)");
// initialize Marsaglia RNG with processor-unique seed
if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
@ -259,7 +259,7 @@ void PairDPDfdtEnergy::init_style()
bool eos_flag = false;
for (int i = 0; i < modify->nfix; i++)
if (strncmp(modify->fix[i]->style,"eos",3) == 0) eos_flag = true;
if(!eos_flag) error->all(FLERR,"pair_style dpd/fdt/energy requires an EOS to be specified");
if(!eos_flag) error->all(FLERR,"pair_style dpd/fdt/energy requires an EOS to be specified");
}
/* ----------------------------------------------------------------------

2
src/USER-DPD/pair_dpd_fdt_energy.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov

510
src/USER-DPD/pair_exp6_rx.cpp
View File

@ -43,12 +43,17 @@ PairExp6rx::PairExp6rx(LAMMPS *lmp) : Pair(lmp)
nparams = maxparam = 0;
params = NULL;
mol2param = NULL;
site1 = site2 = NULL;
}
/* ---------------------------------------------------------------------- */
PairExp6rx::~PairExp6rx()
{
for (int i=0; i < nparams; ++i) {
delete[] params[i].name;
delete[] params[i].potential;
}
memory->destroy(params);
memory->destroy(mol2param);
@ -58,6 +63,8 @@ PairExp6rx::~PairExp6rx()
memory->destroy(cut);
}
delete[] site1;
delete[] site2;
}
/* ---------------------------------------------------------------------- */
@ -66,7 +73,7 @@ void PairExp6rx::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwlOld,fpair;
double rsq,rinv,r2inv,r6inv,forceExp6,factor_lj;
double rsq,r2inv,r6inv,forceExp6,factor_lj;
double rCut,rCutInv,rCut2inv,rCut6inv,rCutExp,urc,durc;
double rm2ij,rm6ij;
double r,rexp;
@ -106,8 +113,8 @@ void PairExp6rx::compute(int eflag, int vflag)
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
const double nRep = double(12.0);
const double shift = double(1.05);
const double nRep = 12.0;
const double shift = 1.05;
double rin1, aRep, uin1, win1, uin1rep, rin1exp, rin6, rin6inv;
inum = list->inum;
@ -125,7 +132,7 @@ void PairExp6rx::compute(int eflag, int vflag)
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
getParamsEXP6(i,epsilon1_i,alpha1_i,rm1_i,fraction1_i,epsilon2_i,alpha2_i,rm2_i,fraction2_i,epsilonOld1_i,alphaOld1_i,rmOld1_i,fractionOld1_i,epsilonOld2_i,alphaOld2_i,rmOld2_i,fractionOld2_i);
for (jj = 0; jj < jnum; jj++) {
@ -141,236 +148,233 @@ void PairExp6rx::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = double(1.0)/rsq;
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
r = sqrt(rsq);
rinv = double(1.0)/r;
rCut2inv = double(1.0)/cutsq[itype][jtype];
rCut6inv = rCut2inv*rCut2inv*rCut2inv;
rCut = sqrt(cutsq[itype][jtype]);
rCutInv = double(1.0)/rCut;
rCut2inv = 1.0/cutsq[itype][jtype];
rCut6inv = rCut2inv*rCut2inv*rCut2inv;
rCut = sqrt(cutsq[itype][jtype]);
rCutInv = 1.0/rCut;
//
// A. Compute the exp-6 potential
//
//
// A. Compute the exp-6 potential
//
// A1. Get alpha, epsilon and rm for particle j
getParamsEXP6(j,epsilon1_j,alpha1_j,rm1_j,fraction1_j,epsilon2_j,alpha2_j,rm2_j,fraction2_j,epsilonOld1_j,alphaOld1_j,rmOld1_j,fractionOld1_j,epsilonOld2_j,alphaOld2_j,rmOld2_j,fractionOld2_j);
// A1. Get alpha, epsilon and rm for particle j
getParamsEXP6(j,epsilon1_j,alpha1_j,rm1_j,fraction1_j,epsilon2_j,alpha2_j,rm2_j,fraction2_j,epsilonOld1_j,alphaOld1_j,rmOld1_j,fractionOld1_j,epsilonOld2_j,alphaOld2_j,rmOld2_j,fractionOld2_j);
// A2. Apply Lorentz-Berthelot mixing rules for the i-j pair
alphaOld12_ij = sqrt(alphaOld1_i*alphaOld2_j);
rmOld12_ij = 0.5*(rmOld1_i + rmOld2_j);
epsilonOld12_ij = sqrt(epsilonOld1_i*epsilonOld2_j);
alphaOld21_ij = sqrt(alphaOld2_i*alphaOld1_j);
rmOld21_ij = 0.5*(rmOld2_i + rmOld1_j);
epsilonOld21_ij = sqrt(epsilonOld2_i*epsilonOld1_j);
// A2. Apply Lorentz-Berthelot mixing rules for the i-j pair
alphaOld12_ij = sqrt(alphaOld1_i*alphaOld2_j);
rmOld12_ij = 0.5*(rmOld1_i + rmOld2_j);
epsilonOld12_ij = sqrt(epsilonOld1_i*epsilonOld2_j);
alphaOld21_ij = sqrt(alphaOld2_i*alphaOld1_j);
rmOld21_ij = 0.5*(rmOld2_i + rmOld1_j);
epsilonOld21_ij = sqrt(epsilonOld2_i*epsilonOld1_j);
alpha12_ij = sqrt(alpha1_i*alpha2_j);
rm12_ij = 0.5*(rm1_i + rm2_j);
epsilon12_ij = sqrt(epsilon1_i*epsilon2_j);
alpha21_ij = sqrt(alpha2_i*alpha1_j);
rm21_ij = 0.5*(rm2_i + rm1_j);
epsilon21_ij = sqrt(epsilon2_i*epsilon1_j);
alpha12_ij = sqrt(alpha1_i*alpha2_j);
rm12_ij = 0.5*(rm1_i + rm2_j);
epsilon12_ij = sqrt(epsilon1_i*epsilon2_j);
alpha21_ij = sqrt(alpha2_i*alpha1_j);
rm21_ij = 0.5*(rm2_i + rm1_j);
epsilon21_ij = sqrt(epsilon2_i*epsilon1_j);
if(rmOld12_ij!=double(0.0) && rmOld21_ij!=double(0.0)){
if(alphaOld21_ij == double(6.0) || alphaOld12_ij == double(6.0))
error->all(FLERR,"alpha_ij is 6.0 in pair exp6");
if(rmOld12_ij!=0.0 && rmOld21_ij!=0.0){
if(alphaOld21_ij == 6.0 || alphaOld12_ij == 6.0)
error->all(FLERR,"alpha_ij is 6.0 in pair exp6");
// A3. Compute some convenient quantities for evaluating the force
rminv = 1.0/rmOld12_ij;
buck1 = epsilonOld12_ij / (alphaOld12_ij - 6.0);
rexp = expValue(alphaOld12_ij*(1.0-r*rminv));
rm2ij = rmOld12_ij*rmOld12_ij;
rm6ij = rm2ij*rm2ij*rm2ij;
// A3. Compute some convenient quantities for evaluating the force
rminv = 1.0/rmOld12_ij;
buck1 = epsilonOld12_ij / (alphaOld12_ij - 6.0);
rexp = expValue(alphaOld12_ij*(1.0-r*rminv));
rm2ij = rmOld12_ij*rmOld12_ij;
rm6ij = rm2ij*rm2ij*rm2ij;
// Compute the shifted potential
rCutExp = expValue(alphaOld12_ij*(1.0-rCut*rminv));
buck2 = 6.0*alphaOld12_ij;
urc = buck1*(6.0*rCutExp - alphaOld12_ij*rm6ij*rCut6inv);
durc = -buck1*buck2*(rCutExp* rminv - rCutInv*rm6ij*rCut6inv);
rin1 = shift*rmOld12_ij*func_rin(alphaOld12_ij);
if(r < rin1){
rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
rin6inv = double(1.0)/rin6;
// Compute the shifted potential
rCutExp = expValue(alphaOld12_ij*(1.0-rCut*rminv));
buck2 = 6.0*alphaOld12_ij;
urc = buck1*(6.0*rCutExp - alphaOld12_ij*rm6ij*rCut6inv);
durc = -buck1*buck2*(rCutExp* rminv - rCutInv*rm6ij*rCut6inv);
rin1 = shift*rmOld12_ij*func_rin(alphaOld12_ij);
if(r < rin1){
rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
rin6inv = 1.0/rin6;
rin1exp = expValue(alphaOld12_ij*(1.0-rin1*rminv));
rin1exp = expValue(alphaOld12_ij*(1.0-rin1*rminv));
uin1 = buck1*(6.0*rin1exp - alphaOld12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
uin1 = buck1*(6.0*rin1exp - alphaOld12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
aRep = double(-1.0)*win1*pow(rin1,nRep)/nRep;
aRep = -1.0*win1*pow(rin1,nRep)/nRep;
uin1rep = aRep/pow(rin1,nRep);
uin1rep = aRep/pow(rin1,nRep);
evdwlOldEXP6_12 = uin1 - uin1rep + aRep/pow(r,nRep);
evdwlOldEXP6_12 = uin1 - uin1rep + aRep/pow(r,nRep);
} else {
evdwlOldEXP6_12 = buck1*(6.0*rexp - alphaOld12_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
}
} else {
evdwlOldEXP6_12 = buck1*(6.0*rexp - alphaOld12_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
}
// A3. Compute some convenient quantities for evaluating the force
rminv = 1.0/rmOld21_ij;
buck1 = epsilonOld21_ij / (alphaOld21_ij - 6.0);
buck2 = 6.0*alphaOld21_ij;
rexp = expValue(alphaOld21_ij*(1.0-r*rminv));
rm2ij = rmOld21_ij*rmOld21_ij;
rm6ij = rm2ij*rm2ij*rm2ij;
// Compute the shifted potential
rCutExp = expValue(alphaOld21_ij*(1.0-rCut*rminv));
buck2 = 6.0*alphaOld21_ij;
urc = buck1*(6.0*rCutExp - alphaOld21_ij*rm6ij*rCut6inv);
durc = -buck1*buck2*(rCutExp* rminv - rCutInv*rm6ij*rCut6inv);
rin1 = shift*rmOld21_ij*func_rin(alphaOld21_ij);
// A3. Compute some convenient quantities for evaluating the force
rminv = 1.0/rmOld21_ij;
buck1 = epsilonOld21_ij / (alphaOld21_ij - 6.0);
buck2 = 6.0*alphaOld21_ij;
rexp = expValue(alphaOld21_ij*(1.0-r*rminv));
rm2ij = rmOld21_ij*rmOld21_ij;
rm6ij = rm2ij*rm2ij*rm2ij;
if(r < rin1){
rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
rin6inv = double(1.0)/rin6;
// Compute the shifted potential
rCutExp = expValue(alphaOld21_ij*(1.0-rCut*rminv));
buck2 = 6.0*alphaOld21_ij;
urc = buck1*(6.0*rCutExp - alphaOld21_ij*rm6ij*rCut6inv);
durc = -buck1*buck2*(rCutExp* rminv - rCutInv*rm6ij*rCut6inv);
rin1 = shift*rmOld21_ij*func_rin(alphaOld21_ij);
rin1exp = expValue(alphaOld21_ij*(1.0-rin1*rminv));
if(r < rin1){
rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
rin6inv = 1.0/rin6;
uin1 = buck1*(6.0*rin1exp - alphaOld21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
rin1exp = expValue(alphaOld21_ij*(1.0-rin1*rminv));
win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
uin1 = buck1*(6.0*rin1exp - alphaOld21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
aRep = double(-1.0)*win1*pow(rin1,nRep)/nRep;
win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
uin1rep = aRep/pow(rin1,nRep);
aRep = -1.0*win1*pow(rin1,nRep)/nRep;
evdwlOldEXP6_21 = uin1 - uin1rep + aRep/pow(r,nRep);
uin1rep = aRep/pow(rin1,nRep);
} else {
evdwlOldEXP6_21 = buck1*(6.0*rexp - alphaOld21_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
}
evdwlOldEXP6_21 = uin1 - uin1rep + aRep/pow(r,nRep);
if (strcmp(site1,site2) == 0)
evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12;
else
evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*evdwlOldEXP6_21;
evdwlOld *= factor_lj;
uCG[i] += double(0.5)*evdwlOld;
uCG[j] += double(0.5)*evdwlOld;
}
} else {
evdwlOldEXP6_21 = buck1*(6.0*rexp - alphaOld21_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
}
if(rm12_ij!=double(0.0) && rm21_ij!=double(0.0)){
if(alpha21_ij == double(6.0) || alpha12_ij == double(6.0))
error->all(FLERR,"alpha_ij is 6.0 in pair exp6");
if (strcmp(site1,site2) == 0)
evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12;
else
evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*evdwlOldEXP6_21;
// A3. Compute some convenient quantities for evaluating the force
rminv = 1.0/rm12_ij;
buck1 = epsilon12_ij / (alpha12_ij - 6.0);
buck2 = 6.0*alpha12_ij;
rexp = expValue(alpha12_ij*(1.0-r*rminv));
rm2ij = rm12_ij*rm12_ij;
rm6ij = rm2ij*rm2ij*rm2ij;
// Compute the shifted potential
rCutExp = expValue(alpha12_ij*(1.0-rCut*rminv));
urc = buck1*(6.0*rCutExp - alpha12_ij*rm6ij*rCut6inv);
durc = -buck1*buck2*(rCutExp*rminv - rCutInv*rm6ij*rCut6inv);
rin1 = shift*rm12_ij*func_rin(alpha12_ij);
evdwlOld *= factor_lj;
if(r < rin1){
rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
rin6inv = double(1.0)/rin6;
uCG[i] += 0.5*evdwlOld;
uCG[j] += 0.5*evdwlOld;
}
rin1exp = expValue(alpha12_ij*(1.0-rin1*rminv));
if(rm12_ij!=0.0 && rm21_ij!=0.0){
if(alpha21_ij == 6.0 || alpha12_ij == 6.0)
error->all(FLERR,"alpha_ij is 6.0 in pair exp6");
uin1 = buck1*(6.0*rin1exp - alpha12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
// A3. Compute some convenient quantities for evaluating the force
rminv = 1.0/rm12_ij;
buck1 = epsilon12_ij / (alpha12_ij - 6.0);
buck2 = 6.0*alpha12_ij;
rexp = expValue(alpha12_ij*(1.0-r*rminv));
rm2ij = rm12_ij*rm12_ij;
rm6ij = rm2ij*rm2ij*rm2ij;
win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
// Compute the shifted potential
rCutExp = expValue(alpha12_ij*(1.0-rCut*rminv));
urc = buck1*(6.0*rCutExp - alpha12_ij*rm6ij*rCut6inv);
durc = -buck1*buck2*(rCutExp*rminv - rCutInv*rm6ij*rCut6inv);
rin1 = shift*rm12_ij*func_rin(alpha12_ij);
aRep = double(-1.0)*win1*pow(rin1,nRep)/nRep;
if(r < rin1){
rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
rin6inv = 1.0/rin6;
uin1rep = aRep/pow(rin1,nRep);
rin1exp = expValue(alpha12_ij*(1.0-rin1*rminv));
evdwlEXP6_12 = uin1 - uin1rep + aRep/pow(r,nRep);
uin1 = buck1*(6.0*rin1exp - alpha12_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
forceExp6 = double(-1.0)*nRep*aRep/pow(r,nRep);
fpairEXP6_12 = factor_lj*forceExp6*r2inv;
win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
} else {
aRep = -1.0*win1*pow(rin1,nRep)/nRep;
// A4. Compute the exp-6 force and energy
forceExp6 = buck1*buck2*(r*rexp*rminv - rm6ij*r6inv) + r*durc;
fpairEXP6_12 = factor_lj*forceExp6*r2inv;
evdwlEXP6_12 = buck1*(6.0*rexp - alpha12_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
}
rminv = 1.0/rm21_ij;
buck1 = epsilon21_ij / (alpha21_ij - 6.0);
buck2 = 6.0*alpha21_ij;
rexp = expValue(alpha21_ij*(1.0-r*rminv));
rm2ij = rm21_ij*rm21_ij;
rm6ij = rm2ij*rm2ij*rm2ij;
uin1rep = aRep/pow(rin1,nRep);
// Compute the shifted potential
rCutExp = expValue(alpha21_ij*(1.0-rCut*rminv));
urc = buck1*(6.0*rCutExp - alpha21_ij*rm6ij*rCut6inv);
durc = -buck1*buck2*(rCutExp*rminv - rCutInv*rm6ij*rCut6inv);
rin1 = shift*rm21_ij*func_rin(alpha21_ij);
evdwlEXP6_12 = uin1 - uin1rep + aRep/pow(r,nRep);
if(r < rin1){
rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
rin6inv = double(1.0)/rin6;
forceExp6 = -1.0*nRep*aRep/pow(r,nRep);
fpairEXP6_12 = factor_lj*forceExp6*r2inv;
rin1exp = expValue(alpha21_ij*(1.0-rin1*rminv));
} else {
uin1 = buck1*(6.0*rin1exp - alpha21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
// A4. Compute the exp-6 force and energy
forceExp6 = buck1*buck2*(r*rexp*rminv - rm6ij*r6inv) + r*durc;
fpairEXP6_12 = factor_lj*forceExp6*r2inv;
evdwlEXP6_12 = buck1*(6.0*rexp - alpha12_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
}
win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
rminv = 1.0/rm21_ij;
buck1 = epsilon21_ij / (alpha21_ij - 6.0);
buck2 = 6.0*alpha21_ij;
rexp = expValue(alpha21_ij*(1.0-r*rminv));
rm2ij = rm21_ij*rm21_ij;
rm6ij = rm2ij*rm2ij*rm2ij;
aRep = double(-1.0)*win1*pow(rin1,nRep)/nRep;
// Compute the shifted potential
rCutExp = expValue(alpha21_ij*(1.0-rCut*rminv));
urc = buck1*(6.0*rCutExp - alpha21_ij*rm6ij*rCut6inv);
durc = -buck1*buck2*(rCutExp*rminv - rCutInv*rm6ij*rCut6inv);
rin1 = shift*rm21_ij*func_rin(alpha21_ij);
uin1rep = aRep/pow(rin1,nRep);
if(r < rin1){
rin6 = rin1*rin1*rin1*rin1*rin1*rin1;
rin6inv = 1.0/rin6;
evdwlEXP6_21 = uin1 - uin1rep + aRep/pow(r,nRep);
rin1exp = expValue(alpha21_ij*(1.0-rin1*rminv));
forceExp6 = double(-1.0)*nRep*aRep/pow(r,nRep);
fpairEXP6_21 = factor_lj*forceExp6*r2inv;
uin1 = buck1*(6.0*rin1exp - alpha21_ij*rm6ij*rin6inv) - urc - durc*(rin1-rCut);
} else {
win1 = buck1*buck2*(rin1*rin1exp*rminv - rm6ij*rin6inv) - rin1*durc;
// A4. Compute the exp-6 force and energy
forceExp6 = buck1*buck2*(r*rexp*rminv - rm6ij*r6inv) + r*durc;
fpairEXP6_21 = factor_lj*forceExp6*r2inv;
evdwlEXP6_21 = buck1*(6.0*rexp - alpha21_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
}
aRep = -1.0*win1*pow(rin1,nRep)/nRep;
//
// Apply Mixing Rule to get the overall force for the CG pair
//
if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairEXP6_12;
else fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*fpairEXP6_21;
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 (strcmp(site1,site2) == 0)
evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12;
else
evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12 + sqrt(fraction2_i*fraction1_j)*evdwlEXP6_21;
evdwl *= factor_lj;
uCGnew[i] += double(0.5)*evdwl;
uCGnew[j] += double(0.5)*evdwl;
evdwl = evdwlOld;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
uin1rep = aRep/pow(rin1,nRep);
evdwlEXP6_21 = uin1 - uin1rep + aRep/pow(r,nRep);
forceExp6 = -1.0*nRep*aRep/pow(r,nRep);
fpairEXP6_21 = factor_lj*forceExp6*r2inv;
} else {
// A4. Compute the exp-6 force and energy
forceExp6 = buck1*buck2*(r*rexp*rminv - rm6ij*r6inv) + r*durc;
fpairEXP6_21 = factor_lj*forceExp6*r2inv;
evdwlEXP6_21 = buck1*(6.0*rexp - alpha21_ij*rm6ij*r6inv) - urc - durc*(r-rCut);
}
//
// Apply Mixing Rule to get the overall force for the CG pair
//
if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairEXP6_12;
else fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*fpairEXP6_21;
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 (strcmp(site1,site2) == 0)
evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12;
else
evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12 + sqrt(fraction2_i*fraction1_j)*evdwlEXP6_21;
evdwl *= factor_lj;
uCGnew[i] += 0.5*evdwl;
uCGnew[j] += 0.5*evdwl;
evdwl = evdwlOld;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}
}
@ -451,11 +455,11 @@ void PairExp6rx::coeff(int narg, char **arg)
}
if (ispecies == nspecies && strcmp(site1,"1fluid") != 0)
error->all(FLERR,"Site1 name not recognized in pair coefficients");
n = strlen(arg[4]) + 1;
site2 = new char[n];
strcpy(site2,arg[4]);
for (ispecies = 0; ispecies < nspecies; ispecies++){
if (strcmp(site2,&atom->dname[ispecies][0]) == 0) break;
}
@ -597,8 +601,8 @@ void PairExp6rx::read_file(char *file)
params[nparams].epsilon = atof(words[3]);
params[nparams].rm = atof(words[4]);
if (params[nparams].epsilon <= 0.0 || params[nparams].rm <= 0.0 ||
params[nparams].alpha < 0.0)
error->all(FLERR,"Illegal exp6/rx parameters. Rm and Epsilon must be greater than zero. Alpha cannot be negative.");
params[nparams].alpha < 0.0)
error->all(FLERR,"Illegal exp6/rx parameters. Rm and Epsilon must be greater than zero. Alpha cannot be negative.");
} else {
error->all(FLERR,"Illegal exp6/rx parameters. Interaction potential does not exist.");
}
@ -624,8 +628,8 @@ void PairExp6rx::setup()
n = -1;
for (j = 0; j < nparams; j++) {
if (i == params[j].ispecies) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = j;
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = j;
}
}
mol2param[i] = n;
@ -713,31 +717,29 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
double rmi, rmj, rmij, rm3ij;
double epsiloni, epsilonj, epsilonij;
double alphai, alphaj, alphaij;
double epsilon_old, rm3_old, alpha_old, fraction_old;
double epsilon, rm3, alpha, fraction;
double epsilon_old, rm3_old, alpha_old;
double epsilon, rm3, alpha;
double fractionOFA, fractionOFA_old;
double nTotalOFA, nTotalOFA_old;
double nTotal, nTotal_old;
double xMolei, xMolej, xMolei_old, xMolej_old;
rm3 = double(0.0);
epsilon = double(0.0);
alpha = double(0.0);
epsilon_old = double(0.0);
rm3_old = double(0.0);
alpha_old = double(0.0);
fraction_old = double(0.0);
fraction = double(0.0);
fractionOFA = double(0.0);
fractionOFA_old = double(0.0);
nTotalOFA = double(0.0);
nTotalOFA_old = double(0.0);
nTotal = double(0.0);
nTotal_old = double(0.0);
rm3 = 0.0;
epsilon = 0.0;
alpha = 0.0;
epsilon_old = 0.0;
rm3_old = 0.0;
alpha_old = 0.0;
fractionOFA = 0.0;
fractionOFA_old = 0.0;
nTotalOFA = 0.0;
nTotalOFA_old = 0.0;
nTotal = 0.0;
nTotal_old = 0.0;
// Compute the total number of molecules in the old and new CG particle as well as the total number of molecules in the fluid portion of the old and new CG particle
for (int ispecies = 0; ispecies < nspecies; ispecies++){
nTotal += atom->dvector[ispecies][id];
nTotal += atom->dvector[ispecies][id];
nTotal_old += atom->dvector[ispecies+nspecies][id];
iparam = mol2param[ispecies];
@ -748,8 +750,8 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
nTotalOFA += atom->dvector[ispecies][id];
}
}
if(nTotal < 1e-8 || nTotal_old < 1e-8)
error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
if(nTotal < 1e-8 || nTotal_old < 1e-8)
error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
// Compute the mole fraction of molecules within the fluid portion of the particle (One Fluid Approximation)
fractionOFA_old = nTotalOFA_old / nTotal_old;
@ -758,7 +760,7 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
iparam = mol2param[ispecies];
if (iparam < 0 || strcmp(params[iparam].potential,"exp6") != 0) continue;
// If Site1 matches a pure species, then grab the parameters
if (strcmp(site1,params[iparam].name) == 0){
rm1_old = params[iparam].rm;
@ -771,7 +773,7 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
// Compute the mole fraction of Site1
fraction1_old = atom->dvector[ispecies+nspecies][id]/nTotal_old;
fraction1 = atom->dvector[ispecies][id]/nTotal;
}
}
// If Site2 matches a pure species, then grab the parameters
if (strcmp(site2,params[iparam].name) == 0){
@ -785,7 +787,7 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
// Compute the mole fraction of Site2
fraction2_old = atom->dvector[ispecies+nspecies][id]/nTotal_old;
fraction2 = atom->dvector[ispecies][id]/nTotal;
}
}
// If Site1 or Site2 matches is a fluid, then compute the paramters
if (strcmp(site1,"1fluid") == 0 || strcmp(site2,"1fluid") == 0) {
@ -797,30 +799,30 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
xMolei_old = atom->dvector[ispecies+nspecies][id]/nTotalOFA_old;
for (int jspecies = 0; jspecies < nspecies; jspecies++) {
jparam = mol2param[jspecies];
if (jparam < 0 || strcmp(params[jparam].potential,"exp6") != 0) continue;
if (strcmp(site1,params[jparam].name) == 0 || strcmp(site2,params[jparam].name) == 0) continue;
rmj = params[jparam].rm;
epsilonj = params[jparam].epsilon;
alphaj = params[jparam].alpha;
xMolej = atom->dvector[jspecies][id]/nTotalOFA;
xMolej_old = atom->dvector[jspecies+nspecies][id]/nTotalOFA_old;
jparam = mol2param[jspecies];
if (jparam < 0 || strcmp(params[jparam].potential,"exp6") != 0) continue;
if (strcmp(site1,params[jparam].name) == 0 || strcmp(site2,params[jparam].name) == 0) continue;
rmj = params[jparam].rm;
epsilonj = params[jparam].epsilon;
alphaj = params[jparam].alpha;
xMolej = atom->dvector[jspecies][id]/nTotalOFA;
xMolej_old = atom->dvector[jspecies+nspecies][id]/nTotalOFA_old;
rmij = (rmi+rmj)/2.0;
rm3ij = rmij*rmij*rmij;
epsilonij = sqrt(epsiloni*epsilonj);
alphaij = sqrt(alphai*alphaj);
if(fractionOFA_old > double(0.0)){
rm3_old += xMolei_old*xMolej_old*rm3ij;
epsilon_old += xMolei_old*xMolej_old*rm3ij*epsilonij;
alpha_old += xMolei_old*xMolej_old*rm3ij*epsilonij*alphaij;
}
if(fractionOFA > double(0.0)){
rm3 += xMolei*xMolej*rm3ij;
epsilon += xMolei*xMolej*rm3ij*epsilonij;
alpha += xMolei*xMolej*rm3ij*epsilonij*alphaij;
}
rmij = (rmi+rmj)/2.0;
rm3ij = rmij*rmij*rmij;
epsilonij = sqrt(epsiloni*epsilonj);
alphaij = sqrt(alphai*alphaj);
if(fractionOFA_old > 0.0){
rm3_old += xMolei_old*xMolej_old*rm3ij;
epsilon_old += xMolei_old*xMolej_old*rm3ij*epsilonij;
alpha_old += xMolei_old*xMolej_old*rm3ij*epsilonij*alphaij;
}
if(fractionOFA > 0.0){
rm3 += xMolei*xMolej*rm3ij;
epsilon += xMolei*xMolej*rm3ij*epsilonij;
alpha += xMolei*xMolej*rm3ij*epsilonij*alphaij;
}
}
}
}
@ -828,9 +830,9 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
if(strcmp(site1,"1fluid") == 0){
rm1 = cbrt(rm3);
if(rm1 < 1e-16) {
rm1 = double(0.0);
epsilon1 = double(0.0);
alpha1 = double(0.0);
rm1 = 0.0;
epsilon1 = 0.0;
alpha1 = 0.0;
} else {
epsilon1 = epsilon / rm3;
alpha1 = alpha / epsilon1 / rm3;
@ -840,9 +842,9 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rm1_old = cbrt(rm3_old);
if(rm1_old < 1e-16) {
rm1_old = double(0.0);
epsilon1_old = double(0.0);
alpha1_old = double(0.0);
rm1_old = 0.0;
epsilon1_old = 0.0;
alpha1_old = 0.0;
} else {
epsilon1_old = epsilon_old / rm3_old;
alpha1_old = alpha_old / epsilon1_old / rm3_old;
@ -880,14 +882,14 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rm1 *= pow(nTotalOFA,fuchslinR);
rm1_old *= pow(nTotalOFA_old,fuchslinR);
}
}
}
if(strcmp(site2,"1fluid") == 0){
rm2 = cbrt(rm3);
if(rm2 < 1e-16) {
rm2 = double(0.0);
epsilon2 = double(0.0);
alpha2 = double(0.0);
rm2 = 0.0;
epsilon2 = 0.0;
alpha2 = 0.0;
} else {
epsilon2 = epsilon / rm3;
alpha2 = alpha / epsilon2 / rm3;
@ -896,9 +898,9 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rm2_old = cbrt(rm3_old);
if(rm2_old < 1e-16) {
rm2_old = double(0.0);
epsilon2_old = double(0.0);
alpha2_old = double(0.0);
rm2_old = 0.0;
epsilon2_old = 0.0;
alpha2_old = 0.0;
} else {
epsilon2_old = epsilon_old / rm3_old;
alpha2_old = alpha_old / epsilon2_old / rm3_old;
@ -945,13 +947,11 @@ double PairExp6rx::func_rin(double &alpha)
{
double function;
const double alpha_min = double(11.0);
const double alpha_max = double(16.0);
const double a = double(3.7682065);
const double b = double(-1.4308614);
function = a+b*pow(alpha,0.5);
function = expValue(function);
const double a = 3.7682065;
const double b = -1.4308614;
function = a+b*sqrt(alpha);
function = expValue(function);
return function;
}
@ -961,8 +961,8 @@ double PairExp6rx::func_rin(double &alpha)
double PairExp6rx::expValue(double value)
{
double returnValue;
if(value < DBL_MIN_EXP) returnValue = double(0.0);
if(value < DBL_MIN_EXP) returnValue = 0.0;
else returnValue = exp(value);
return returnValue;
}

4
src/USER-DPD/pair_exp6_rx.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -130,7 +130,7 @@ Self-explanatory
E: The number of molecules in CG particle is less than 1e-8.
Self-explanatory. Check the species concentrations have been properly set
and check the reaction kinetic solver parameters in fix rx to more for
and check the reaction kinetic solver parameters in fix rx to more for
sufficient accuracy.

50
src/USER-DPD/pair_multi_lucy.cpp
View File

@ -12,11 +12,11 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------------------------
Contributing authors:
Contributing authors:
James Larentzos and Joshua Moore (U.S. Army Research Laboratory)
Please cite the related publications:
J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
"A coarse-grain force field for RDX: Density dependent and energy conserving"
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
@ -87,12 +87,10 @@ PairMultiLucy::~PairMultiLucy()
void PairMultiLucy::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype,itable;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,factor_lj,fraction,value,a,b;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
Table *tb;
union_int_float_t rsq_lookup;
int tlm1 = tablength - 1;
evdwl = 0.0;
@ -103,14 +101,12 @@ void PairMultiLucy::compute(int eflag, int vflag)
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double pi = MathConst::MY_PI;
double A_i;
double A_j;
double fraction_i,fraction_j;
double a_i,a_j,b_i,b_j;
int jtable;
double *rho = atom->rho;
@ -134,7 +130,6 @@ void PairMultiLucy::compute(int eflag, int vflag)
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
@ -151,19 +146,19 @@ void PairMultiLucy::compute(int eflag, int vflag)
if (tabstyle == LOOKUP) {
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
if (itable >= tlm1 || jtable >= tlm1)
error->one(FLERR,"Density > table outer cutoff");
if (itable >= tlm1 || jtable >= tlm1)
error->one(FLERR,"Density > table outer cutoff");
A_i = tb->f[itable];
A_j = tb->f[jtable];
fpair = double(0.5)*(A_i + A_j)*(double(1.0)+double(3.0)*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair/sqrt(rsq);
} else if (tabstyle == LINEAR) {
itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
if (itable >= tlm1 || jtable >= tlm1)
error->one(FLERR,"Density > table outer cutoff");
if (itable >= tlm1 || jtable >= tlm1)
error->one(FLERR,"Density > table outer cutoff");
fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
fraction_j = (((rho[j]*rho[j]) - tb->rsq[jtable]) * tb->invdelta);
@ -171,7 +166,7 @@ void PairMultiLucy::compute(int eflag, int vflag)
A_i = tb->f[itable] + fraction_i*tb->df[itable];
A_j = tb->f[jtable] + fraction_j*tb->df[jtable];
fpair = double(0.5)*(A_i + A_j)*(double(1.0)+double(3.0)*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair / sqrt(rsq);
@ -185,8 +180,8 @@ void PairMultiLucy::compute(int eflag, int vflag)
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,0.0,fpair,delx,dely,delz);
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,0.0,fpair,delx,dely,delz);
}
}
@ -195,17 +190,17 @@ void PairMultiLucy::compute(int eflag, int vflag)
error->one(FLERR,"Density < table inner cutoff");
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
if (tabstyle == LOOKUP) evdwl = tb->e[itable];
else if (tabstyle == LINEAR){
else if (tabstyle == LINEAR){
if (itable >= tlm1) error->one(FLERR,"Density > table outer cutoff");
if(itable==0) fraction_i=double(0.0);
if(itable==0) fraction_i=0.0;
else fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/double(84.0);
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
if (evflag) ev_tally(0,0,nlocal,newton_pair,
evdwl,0.0,0.0,0.0,0.0,0.0);
evdwl,0.0,0.0,0.0,0.0,0.0);
}
if (vflag_fdotr) virial_fdotr_compute();
@ -293,13 +288,11 @@ void PairMultiLucy::coeff(int narg, char **arg)
// insure cutoff is within table
if (tb->ninput <= 1) error->one(FLERR,"Invalid pair table length");
double rlo,rhi;
double rlo;
if (tb->rflag == 0) {
rlo = tb->rfile[0];
rhi = tb->rfile[tb->ninput-1];
} else {
rlo = tb->rlo;
rhi = tb->rhi;
}
rho_0 = rlo;
@ -390,7 +383,6 @@ void PairMultiLucy::read_table(Table *tb, char *file, char *keyword)
int itmp;
double rtmp;
union_int_float_t rsq_lookup;
fgets(line,MAXLINE,fp);
for (int i = 0; i < tb->ninput; i++) {
@ -733,7 +725,7 @@ void PairMultiLucy::computeLocalDensity()
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double factor, volume;
double factor;
inum = list->inum;
ilist = list->ilist;
@ -742,9 +734,9 @@ void PairMultiLucy::computeLocalDensity()
double pi = MathConst::MY_PI;
double *rho = atom->rho;
// zero out density
if (newton_pair) {
m = nlocal + atom->nghost;
for (i = 0; i < m; i++) rho[i] = 0.0;
@ -773,7 +765,7 @@ void PairMultiLucy::computeLocalDensity()
if (rsq < cutsq[itype][jtype]) {
double rcut = sqrt(cutsq[itype][jtype]);
factor= (double(84.0)/(double(5.0)*pi*rcut*rcut*rcut))*(double(1.0)+double(3.0)*sqrt(rsq)/(double(2.0)*rcut))*(double(1.0)-sqrt(rsq)/rcut)*(double(1.0)-sqrt(rsq)/rcut)*(double(1.0)-sqrt(rsq)/rcut)*(double(1.0)-sqrt(rsq)/rcut);
factor= (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut);
rho[i] += factor;
if (newton_pair || j < nlocal) {
rho[j] += factor;

2
src/USER-DPD/pair_multi_lucy.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov

140
src/USER-DPD/pair_multi_lucy_rx.cpp
View File

@ -12,15 +12,15 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------------------------
Contributing authors:
Contributing authors:
James Larentzos and Joshua Moore (U.S. Army Research Laboratory)
Please cite the related publications:
J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
"A coarse-grain force field for RDX: Density dependent and energy conserving"
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
#include "mpi.h"
#include <math.h>
#include "math_const.h"
@ -90,7 +90,7 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype,itable;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,evdwlOld,fpair;
double rsq,factor_lj;
double rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
Table *tb;
@ -105,7 +105,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double fractionOld1_i,fractionOld1_j;
@ -143,7 +142,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
@ -153,63 +151,63 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
fpair = double(0.0);
getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
fpair = 0.0;
getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
tb = &tables[tabindex[itype][jtype]];
if (rho[i]*rho[i] < tb->innersq || rho[j]*rho[j] < tb->innersq){
printf("Table inner cutoff = %lf\n",sqrt(tb->innersq));
printf("rho[%d]=%lf\n",i,rho[i]);
printf("rho[%d]=%lf\n",j,rho[j]);
printf("Table inner cutoff = %lf\n",sqrt(tb->innersq));
printf("rho[%d]=%lf\n",i,rho[i]);
printf("rho[%d]=%lf\n",j,rho[j]);
error->one(FLERR,"Density < table inner cutoff");
}
}
if (tabstyle == LOOKUP) {
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
if (itable >= tlm1 || jtable >= tlm1){
printf("Table outer index = %d\n",tlm1);
printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
printf("jtableIndex=%d rho[%d]=%lf\n",jtable,j,rho[j]);
error->one(FLERR,"Density > table outer cutoff");
}
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
if (itable >= tlm1 || jtable >= tlm1){
printf("Table outer index = %d\n",tlm1);
printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
printf("jtableIndex=%d rho[%d]=%lf\n",jtable,j,rho[j]);
error->one(FLERR,"Density > table outer cutoff");
}
A_i = tb->f[itable];
A_j = tb->f[jtable];
fpair = double(0.5)*(A_i + A_j)*(double(1.0)+double(3.0)*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair/sqrt(rsq);
} else if (tabstyle == LINEAR) {
itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
if (itable >= tlm1 || jtable >= tlm1){
printf("Table outer index = %d\n",tlm1);
printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
printf("jtableIndex=%d rho[%d]=%lf\n",jtable,j,rho[j]);
error->one(FLERR,"Density > table outer cutoff");
}
if(itable<0) itable=0;
if(itable>=tlm1) itable=tlm1;
if(jtable<0) jtable=0;
if(jtable>=tlm1)jtable=tlm1;
itable = static_cast<int> ((rho[i]*rho[i] - tb->innersq) * tb->invdelta);
jtable = static_cast<int> (((rho[j]*rho[j]) - tb->innersq) * tb->invdelta);
if (itable >= tlm1 || jtable >= tlm1){
printf("Table outer index = %d\n",tlm1);
printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
printf("jtableIndex=%d rho[%d]=%lf\n",jtable,j,rho[j]);
error->one(FLERR,"Density > table outer cutoff");
}
if(itable<0) itable=0;
if(itable>=tlm1) itable=tlm1;
if(jtable<0) jtable=0;
if(jtable>=tlm1)jtable=tlm1;
fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
fraction_j = (((rho[j]*rho[j]) - tb->rsq[jtable]) * tb->invdelta);
if(itable==0) fraction_i=double(0.0);
if(itable==tlm1) fraction_i=double(0.0);
if(jtable==0) fraction_j=double(0.0);
if(jtable==tlm1) fraction_j=double(0.0);
if(itable==0) fraction_i=0.0;
if(itable==tlm1) fraction_i=0.0;
if(jtable==0) fraction_j=0.0;
if(jtable==tlm1) fraction_j=0.0;
A_i = tb->f[itable] + fraction_i*tb->df[itable];
A_j = tb->f[jtable] + fraction_j*tb->df[jtable];
fpair = double(0.5)*(A_i + A_j)*(double(1.0)+double(3.0)*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(double(1.0) - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = 0.5*(A_i + A_j)*(1.0+3.0*sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]))*(1.0 - sqrt(rsq)/sqrt(cutsq[itype][jtype]));
fpair = fpair / sqrt(rsq);
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
@ -219,35 +217,35 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,0.0,fpair,delx,dely,delz);
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,0.0,fpair,delx,dely,delz);
}
}
tb = &tables[tabindex[itype][itype]];
itable = static_cast<int> (((rho[i]*rho[i]) - tb->innersq) * tb->invdelta);
if (tabstyle == LOOKUP) evdwl = tb->e[itable];
else if (tabstyle == LINEAR){
else if (tabstyle == LINEAR){
if (itable >= tlm1){
printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
error->one(FLERR,"Density > table outer cutoff");
printf("itableIndex=%d rho[%d]=%lf\n",itable,i,rho[i]);
error->one(FLERR,"Density > table outer cutoff");
}
if(itable==0) fraction_i=double(0.0);
if(itable==0) fraction_i=0.0;
else fraction_i = (((rho[i]*rho[i]) - tb->rsq[itable]) * tb->invdelta);
evdwl = tb->e[itable] + fraction_i*tb->de[itable];
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/double(84.0);
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
evdwlOld = fractionOld1_i*evdwl;
evdwl = fraction1_i*evdwl;
uCG[i] += evdwlOld;
uCGnew[i] += evdwl;
evdwl = evdwlOld;
if (evflag) ev_tally(0,0,nlocal,newton_pair,
evdwl,0.0,0.0,0.0,0.0,0.0);
evdwl,0.0,0.0,0.0,0.0,0.0);
}
if (vflag_fdotr) virial_fdotr_compute();
@ -337,7 +335,7 @@ void PairMultiLucyRX::coeff(int narg, char **arg)
n = strlen(arg[3]) + 1;
site1 = new char[n];
strcpy(site1,arg[4]);
n = strlen(arg[4]) + 1;
site2 = new char[n];
strcpy(site2,arg[5]);
@ -352,13 +350,11 @@ void PairMultiLucyRX::coeff(int narg, char **arg)
// insure cutoff is within table
if (tb->ninput <= 1) error->one(FLERR,"Invalid pair table length");
double rlo,rhi;
double rlo;
if (tb->rflag == 0) {
rlo = tb->rfile[0];
rhi = tb->rfile[tb->ninput-1];
} else {
rlo = tb->rlo;
rhi = tb->rhi;
}
rho_0 = rlo;
@ -412,7 +408,7 @@ void PairMultiLucyRX::read_table(Table *tb, char *file, char *keyword)
// open file
FILE *fp = fopen(file,"r");
FILE *fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
@ -461,7 +457,7 @@ void PairMultiLucyRX::read_table(Table *tb, char *file, char *keyword)
rtmp = tb->rlo*tb->rlo +
(tb->rhi*tb->rhi - tb->rlo*tb->rlo)*i/(tb->ninput-1);
rtmp = sqrt(rtmp);
}
}
tb->rfile[i] = rtmp;
}
@ -800,9 +796,9 @@ void PairMultiLucyRX::computeLocalDensity()
double pi = MathConst::MY_PI;
double *rho = atom->rho;
// zero out density
if (newton_pair) {
m = nlocal + atom->nghost;
for (i = 0; i < m; i++) rho[i] = 0.0;
@ -831,9 +827,9 @@ void PairMultiLucyRX::computeLocalDensity()
if (rsq < cutsq[itype][jtype]) {
double rcut = sqrt(cutsq[itype][jtype]);
double tmpFactor = double(1.0)-sqrt(rsq)/rcut;
double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
factor = (double(84.0)/(double(5.0)*pi*rcut*rcut*rcut))*(double(1.0)+double(3.0)*sqrt(rsq)/(double(2.0)*rcut))*tmpFactor4;
double tmpFactor = 1.0-sqrt(rsq)/rcut;
double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor;
factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4;
rho[i] += factor;
if (newton_pair || j < nlocal) {
rho[j] += factor;
@ -854,17 +850,17 @@ void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOl
double fractionOld, fraction;
double nTotal, nTotalOld;
fractionOld = double(0.0);
fraction = double(0.0);
fractionOld1 = double(0.0);
fractionOld2 = double(0.0);
fraction1 = double(0.0);
fraction2 = double(0.0);
fractionOld = 0.0;
fraction = 0.0;
fractionOld1 = 0.0;
fractionOld2 = 0.0;
fraction1 = 0.0;
fraction2 = 0.0;
nTotal = double(0.0);
nTotalOld = double(0.0);
nTotal = 0.0;
nTotalOld = 0.0;
for(int ispecies=0;ispecies<nspecies;ispecies++){
nTotal += atom->dvector[ispecies][id];
nTotal += atom->dvector[ispecies][id];
nTotalOld += atom->dvector[ispecies+nspecies][id];
}
@ -880,7 +876,7 @@ void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOl
}
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
fractionOld1 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction1 = atom->dvector[ispecies][id]/nTotal;
}

2
src/USER-DPD/pair_multi_lucy_rx.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov

100
src/USER-DPD/pair_table_rx.cpp
View File

@ -70,9 +70,9 @@ void PairTableRX::compute(int eflag, int vflag)
union_int_float_t rsq_lookup;
int tlm1 = tablength - 1;
fraction = double(0.0);
a = double(0.0);
b = double(0.0);
fraction = 0.0;
a = 0.0;
b = 0.0;
evdwlOld = 0.0;
evdwl = 0.0;
@ -122,7 +122,7 @@ void PairTableRX::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
tb = &tables[tabindex[itype][jtype]];
if (rsq < tb->innersq)
@ -158,8 +158,8 @@ void PairTableRX::compute(int eflag, int vflag)
value = tb->f[itable] + fraction*tb->df[itable];
fpair = factor_lj * value;
}
if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
if (strcmp(site1,site2) == 0) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
@ -170,29 +170,29 @@ void PairTableRX::compute(int eflag, int vflag)
f[j][2] -= delz*fpair;
}
if (tabstyle == LOOKUP)
evdwl = tb->e[itable];
else if (tabstyle == LINEAR || tabstyle == BITMAP){
evdwl = tb->e[itable] + fraction*tb->de[itable];
}
else
evdwl = a * tb->e[itable] + b * tb->e[itable+1] +
((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
tb->deltasq6;
if (strcmp(site1,site2) == 0){
evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwl;
evdwl = sqrt(fraction1_i*fraction2_j)*evdwl;
} else {
evdwlOld = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*evdwl;
evdwl = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*evdwl;
}
evdwlOld *= factor_lj;
evdwl *= factor_lj;
uCG[i] += double(0.5)*evdwlOld;
uCG[j] += double(0.5)*evdwlOld;
uCGnew[i] += double(0.5)*evdwl;
uCGnew[j] += double(0.5)*evdwl;
evdwl = evdwlOld;
if (tabstyle == LOOKUP)
evdwl = tb->e[itable];
else if (tabstyle == LINEAR || tabstyle == BITMAP){
evdwl = tb->e[itable] + fraction*tb->de[itable];
}
else
evdwl = a * tb->e[itable] + b * tb->e[itable+1] +
((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
tb->deltasq6;
if (strcmp(site1,site2) == 0){
evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwl;
evdwl = sqrt(fraction1_i*fraction2_j)*evdwl;
} else {
evdwlOld = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*evdwl;
evdwl = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*evdwl;
}
evdwlOld *= factor_lj;
evdwl *= factor_lj;
uCG[i] += 0.5*evdwlOld;
uCG[j] += 0.5*evdwlOld;
uCGnew[i] += 0.5*evdwl;
uCGnew[j] += 0.5*evdwl;
evdwl = evdwlOld;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
@ -309,7 +309,7 @@ void PairTableRX::coeff(int narg, char **arg)
}
if (ispecies == nspecies && strcmp(site1,"1fluid") != 0)
error->all(FLERR,"Site1 name not recognized in pair coefficients");
n = strlen(arg[4]) + 1;
site2 = new char[n];
strcpy(site2,arg[5]);
@ -599,7 +599,7 @@ void PairTableRX::compute_table(Table *tb)
else inner = tb->rfile[0];
tb->innersq = double(inner)*double(inner);
tb->delta = double(tb->cut*tb->cut - double(tb->innersq)) / double(tlm1);
tb->invdelta = double(1.0)/double(tb->delta);
tb->invdelta = 1.0/double(tb->delta);
// direct lookup tables
// N-1 evenly spaced bins in rsq from inner to cut
@ -988,9 +988,9 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
double fractionOld1_i, fractionOld1_j;
double fractionOld2_i, fractionOld2_j;
fraction = double(0.0);
a = double(0.0);
b = double(0.0);
fraction = 0.0;
a = 0.0;
b = 0.0;
getParams(i,fractionOld1_i,fractionOld2_i,fraction1_i,fraction2_i);
getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
@ -1026,7 +1026,7 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
fforce = factor_lj * value;
}
if (strcmp(site1,site2) == 0) fforce = sqrt(fraction1_i*fraction2_j)*fforce;
if (strcmp(site1,site2) == 0) fforce = sqrt(fraction1_i*fraction2_j)*fforce;
else fforce = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*fforce;
if (tabstyle == LOOKUP)
@ -1071,39 +1071,39 @@ void PairTableRX::getParams(int id, double &fractionOld1, double &fractionOld2,
double fractionOld, fraction;
double nTotal, nTotalOld;
fractionOld = double(0.0);
fraction = double(0.0);
fractionOld1 = double(0.0);
fractionOld2 = double(0.0);
fraction1 = double(0.0);
fraction2 = double(0.0);
fractionOld = 0.0;
fraction = 0.0;
fractionOld1 = 0.0;
fractionOld2 = 0.0;
fraction1 = 0.0;
fraction2 = 0.0;
nTotal = double(0.0);
nTotalOld = double(0.0);
nTotal = 0.0;
nTotalOld = 0.0;
for(int ispecies=0;ispecies<nspecies;ispecies++){
nTotal += atom->dvector[ispecies][id];
nTotal += atom->dvector[ispecies][id];
nTotalOld += atom->dvector[ispecies+nspecies][id];
}
if(nTotal < 1e-8 || nTotalOld < 1e-8)
error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
if(nTotal < 1e-8 || nTotalOld < 1e-8)
error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
fractionOld1 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction1 = atom->dvector[ispecies][id]/nTotal;
}
if (strcmp(site2,&atom->dname[ispecies][0]) == 0){
if (strcmp(site2,&atom->dname[ispecies][0]) == 0){
fractionOld2 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction2 = atom->dvector[ispecies][id]/nTotal;
}
}
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
if (strcmp(site1,&atom->dname[ispecies][0]) == 0){
fractionOld1 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction1 = atom->dvector[ispecies][id]/nTotal;
}
if (strcmp(site2,&atom->dname[ispecies][0]) == 0){
if (strcmp(site2,&atom->dname[ispecies][0]) == 0){
fractionOld2 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction2 = atom->dvector[ispecies][id]/nTotal;
}

4
src/USER-DPD/pair_table_rx.h
View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -166,7 +166,7 @@ long-range solver starts at that cutoff.
E: The number of molecules in CG particle is less than 1e-8
Self-explanatory. Check the species concentrations have been properly set
and check the reaction kinetic solver parameters in fix rx to more for
and check the reaction kinetic solver parameters in fix rx to more for
sufficient accuracy.
*/

4
src/USER-INTEL/improper_cvff_intel.cpp
View File

@ -364,10 +364,8 @@ void ImproperCvffIntel::eval(const int vflag,
/* ---------------------------------------------------------------------- */
void ImproperCvffIntel::init()
void ImproperCvffIntel::init_style()
{
ImproperCvff::init();
int ifix = modify->find_fix("package_intel");
if (ifix < 0)
error->all(FLERR,

2
src/USER-INTEL/improper_cvff_intel.h
View File

@ -35,7 +35,7 @@ class ImproperCvffIntel : public ImproperCvff {
ImproperCvffIntel(class LAMMPS *);
virtual ~ImproperCvffIntel();
virtual void compute(int, int);
virtual void init();
virtual void init_style();
protected:
FixIntel *fix;

4
src/USER-INTEL/improper_harmonic_intel.cpp
View File

@ -325,10 +325,8 @@ void ImproperHarmonicIntel::eval(const int vflag,
/* ---------------------------------------------------------------------- */
void ImproperHarmonicIntel::init()
void ImproperHarmonicIntel::init_style()
{
ImproperHarmonic::init();
int ifix = modify->find_fix("package_intel");
if (ifix < 0)
error->all(FLERR,

2
src/USER-INTEL/improper_harmonic_intel.h
View File

@ -35,7 +35,7 @@ class ImproperHarmonicIntel : public ImproperHarmonic {
ImproperHarmonicIntel(class LAMMPS *);
virtual ~ImproperHarmonicIntel();
virtual void compute(int, int);
virtual void init();
virtual void init_style();
protected:
FixIntel *fix;

16
src/USER-MISC/fix_addtorque.cpp
View File

@ -48,6 +48,8 @@ FixAddTorque::FixAddTorque(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
xstr = ystr = zstr = NULL;
@ -133,8 +135,10 @@ void FixAddTorque::init()
varflag = EQUAL;
else varflag = CONSTANT;
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -144,9 +148,9 @@ void FixAddTorque::setup(int vflag)
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -249,7 +253,7 @@ void FixAddTorque::post_force(int vflag)
void FixAddTorque::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/USER-MISC/fix_addtorque.h
View File

@ -45,7 +45,7 @@ class FixAddTorque : public Fix {
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
double foriginal[4],foriginal_all[4];
int force_flag;
int nlevels_respa;
int ilevel_respa;
};
}

35
src/USER-MISC/fix_smd.cpp
View File

@ -56,6 +56,8 @@ FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) :
size_vector = 7;
global_freq = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
int argoffs=3;
if (strcmp(arg[argoffs],"cvel") == 0) {
@ -156,8 +158,10 @@ void FixSMD::init()
zn = dz/r_old;
}
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -167,9 +171,9 @@ void FixSMD::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -180,7 +184,12 @@ void FixSMD::post_force(int vflag)
if (styleflag & SMD_TETHER) smd_tether();
else smd_couple();
if (styleflag & SMD_CVEL) r_old += v_smd * update->dt;
if (styleflag & SMD_CVEL) {
if (strstr(update->integrate_style,"verlet"))
r_old += v_smd * update->dt;
else
r_old += v_smd * ((Respa *) update->integrate)->step[ilevel_respa];
}
}
/* ---------------------------------------------------------------------- */
@ -190,6 +199,10 @@ void FixSMD::smd_tether()
double xcm[3];
group->xcm(igroup,masstotal,xcm);
double dt = update->dt;
if (strstr(update->integrate_style,"respa"))
dt = ((Respa *) update->integrate)->step[ilevel_respa];
// fx,fy,fz = components of k * (r-r0)
double dx,dy,dz,fx,fy,fz,r,dr;
@ -209,7 +222,7 @@ void FixSMD::smd_tether()
fx = k_smd*dx*dr/r;
fy = k_smd*dy*dr/r;
fz = k_smd*dz*dr/r;
pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt;
pmf += (fx*xn + fy*yn + fz*zn) * v_smd * dt;
} else {
fx = 0.0;
fy = 0.0;
@ -269,6 +282,10 @@ void FixSMD::smd_couple()
group->xcm(igroup,masstotal,xcm);
group->xcm(igroup2,masstotal2,xcm2);
double dt = update->dt;
if (strstr(update->integrate_style,"respa"))
dt = ((Respa *) update->integrate)->step[ilevel_respa];
// renormalize direction of spring
double dx,dy,dz,r,dr;
if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0];
@ -309,7 +326,7 @@ void FixSMD::smd_couple()
fx = k_smd*dx*dr/r;
fy = k_smd*dy*dr/r;
fz = k_smd*dz*dr/r;
pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * update->dt;
pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * dt;
} else {
fx = 0.0;
fy = 0.0;
@ -417,7 +434,7 @@ void FixSMD::restart(char *buf)
void FixSMD::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ----------------------------------------------------------------------

2
src/USER-MISC/fix_smd.h
View File

@ -46,7 +46,7 @@ class FixSMD : public Fix {
int igroup2,group2bit;
double masstotal,masstotal2;
int nlevels_respa;
int ilevel_respa;
double ftotal[3],ftotal_all[7];
int force_flag;

69
src/USER-MISC/fix_ttm_mod.cpp
View File

@ -78,24 +78,33 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) :
nevery = 1;
restart_peratom = 1;
restart_global = 1;
seed = atoi(arg[3]);
fpr_2 = fopen(arg[4],"r");
nxnodes = atoi(arg[5]);
nynodes = atoi(arg[6]);
nznodes = atoi(arg[7]);
fpr = fopen(arg[8],"r");
if (fpr == NULL) {
seed = force->inumeric(FLERR,arg[3]);
if (seed <= 0)
error->all(FLERR,"Invalid random number seed in fix ttm/mod command");
FILE *fpr_2 = force->open_potential(arg[4]);
if (fpr_2 == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",arg[7]);
error->one(FLERR,str);
sprintf(str,"Cannot open file %s",arg[4]);
error->all(FLERR,str);
}
nxnodes = force->inumeric(FLERR,arg[5]);
nynodes = force->inumeric(FLERR,arg[6]);
nznodes = force->inumeric(FLERR,arg[7]);
if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
error->all(FLERR,"Fix ttm/mod number of nodes must be > 0");
FILE *fpr = force->open_potential(arg[8]);
if (fpr == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",arg[8]);
error->one(FLERR,str);
error->all(FLERR,str);
}
nfileevery = atoi(arg[9]);
if (nfileevery) {
nfileevery = force->inumeric(FLERR,arg[9]);
if (nfileevery > 0) {
if (narg != 11) error->all(FLERR,"Illegal fix ttm/mod command");
MPI_Comm_rank(world,&me);
if (me == 0) {
@ -226,19 +235,15 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) :
mult_factor = intensity;
duration = 0.0;
v_0_sq = v_0*v_0;
// error checks
if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
error->all(FLERR,"Fix ttm number of nodes must be > 0");
surface_double = double(t_surface_l)*(domain->xprd/nxnodes);
if ((C_limit+esheat_0) < 0.0)
error->all(FLERR,"Fix ttm electronic_specific_heat must be >= 0.0");
error->all(FLERR,"Fix ttm/mod electronic_specific_heat must be >= 0.0");
if (electronic_density <= 0.0)
error->all(FLERR,"Fix ttm electronic_density must be > 0.0");
if (gamma_p < 0.0) error->all(FLERR,"Fix ttm gamma_p must be >= 0.0");
if (gamma_s < 0.0) error->all(FLERR,"Fix ttm gamma_s must be >= 0.0");
if (v_0 < 0.0) error->all(FLERR,"Fix ttm v_0 must be >= 0.0");
if (ionic_density <= 0.0) error->all(FLERR,"Fix ttm ionic_density must be > 0.0");
if (seed <= 0) error->all(FLERR,"Invalid random number seed in fix ttm command");
error->all(FLERR,"Fix ttm/mod electronic_density must be > 0.0");
if (gamma_p < 0.0) error->all(FLERR,"Fix ttm/mod gamma_p must be >= 0.0");
if (gamma_s < 0.0) error->all(FLERR,"Fix ttm/mod gamma_s must be >= 0.0");
if (v_0 < 0.0) error->all(FLERR,"Fix ttm/mod v_0 must be >= 0.0");
if (ionic_density <= 0.0) error->all(FLERR,"Fix ttm/mod ionic_density must be > 0.0");
if (surface_l < 0) error->all(FLERR,"Surface coordinates must be >= 0");
if (surface_l >= surface_r) error->all(FLERR, "Left surface coordinate must be less than right surface coordinate");
// initialize Marsaglia RNG with processor-unique seed
@ -274,8 +279,10 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) :
atom->add_callback(0);
atom->add_callback(1);
// set initial electron temperatures from user input file
if (me == 0) read_initial_electron_temperatures();
if (me == 0) read_initial_electron_temperatures(fpr);
MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
fclose(fpr);
fclose(fpr_2);
}
/* ---------------------------------------------------------------------- */
@ -317,11 +324,11 @@ int FixTTMMod::setmask()
void FixTTMMod::init()
{
if (domain->dimension == 2)
error->all(FLERR,"Cannot use fix ttm with 2d simulation");
error->all(FLERR,"Cannot use fix ttm/mod with 2d simulation");
if (domain->nonperiodic != 0)
error->all(FLERR,"Cannot use nonperiodic boundares with fix ttm");
error->all(FLERR,"Cannot use nonperiodic boundares with fix ttm/mod");
if (domain->triclinic)
error->all(FLERR,"Cannot use fix ttm with triclinic box");
error->all(FLERR,"Cannot use fix ttm/mod with triclinic box");
// set force prefactors
for (int i = 1; i <= atom->ntypes; i++) {
gfactor1[i] = - gamma_p / force->ftm2v;
@ -488,7 +495,7 @@ void FixTTMMod::reset_dt()
only called by proc 0
------------------------------------------------------------------------- */
void FixTTMMod::read_initial_electron_temperatures()
void FixTTMMod::read_initial_electron_temperatures(FILE *fpr)
{
char line[MAXLINE];
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
@ -501,7 +508,7 @@ void FixTTMMod::read_initial_electron_temperatures()
while (1) {
if (fgets(line,MAXLINE,fpr) == NULL) break;
sscanf(line,"%d %d %d %lg",&ixnode,&iynode,&iznode,&T_tmp);
if (T_tmp < 0.0) error->one(FLERR,"Fix ttm electron temperatures must be >= 0.0");
if (T_tmp < 0.0) error->one(FLERR,"Fix ttm/mod electron temperatures must be >= 0.0");
T_electron[ixnode][iynode][iznode] = T_tmp;
T_initial_set[ixnode][iynode][iznode] = 1;
}
@ -509,9 +516,7 @@ void FixTTMMod::read_initial_electron_temperatures()
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
if (T_initial_set[ixnode][iynode][iznode] == 0)
error->one(FLERR,"Initial temperatures not all set in fix ttm");
// close file
fclose(fpr);
error->one(FLERR,"Initial temperatures not all set in fix ttm/mod");
}
/* ---------------------------------------------------------------------- */
@ -634,7 +639,7 @@ void FixTTMMod::end_of_step()
num_inner_timesteps = static_cast<unsigned int>(update->dt/inner_dt) + 1;
inner_dt = update->dt/double(num_inner_timesteps);
if (num_inner_timesteps > 1000000)
error->warning(FLERR,"Too many inner timesteps in fix ttm",0);
error->warning(FLERR,"Too many inner timesteps in fix ttm/mod",0);
for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
ith_inner_timestep++) {
for (int ixnode = 0; ixnode < nxnodes; ixnode++)

4
src/USER-MISC/fix_ttm_mod.h
View File

@ -58,7 +58,7 @@ class FixTTMMod : public Fix {
int nlevels_respa;
int seed;
class RanMars *random;
FILE *fp,*fpr,*fpr_2;
FILE *fp;
int nxnodes,nynodes,nznodes,total_nnodes;
int ***nsum;
int ***nsum_all,***T_initial_set;
@ -79,7 +79,7 @@ class FixTTMMod : public Fix {
double electron_temperature_min;
el_heat_capacity_thermal_conductivity el_properties(double);
double el_sp_heat_integral(double);
void read_initial_electron_temperatures();
void read_initial_electron_temperatures(FILE *);
};
}

2
src/USER-OMP/fix_gravity_omp.cpp
View File

@ -108,5 +108,5 @@ void FixGravityOMP::post_force(int vflag)
void FixGravityOMP::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}

6
src/USER-OMP/fix_qeq_comb_omp.cpp
View File

@ -63,8 +63,10 @@ void FixQEQCombOMP::init()
error->all(FLERR,"Must use pair_style comb or "
"comb/omp with fix qeq/comb/omp");
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");

13
src/USER-OMP/respa_omp.cpp
View File

@ -73,7 +73,18 @@ void RespaOMP::setup()
fprintf(screen,"Setting up r-RESPA/omp run ...\n");
fprintf(screen," Unit style : %s\n", update->unit_style);
fprintf(screen," Current step : " BIGINT_FORMAT "\n", update->ntimestep);
fprintf(screen," OuterTime step: %g\n", update->dt);
fprintf(screen," Time steps :");
for (int ilevel=0; ilevel < nlevels; ++ilevel)
fprintf(screen," %d:%g",ilevel+1, step[ilevel]);
fprintf(screen,"\n r-RESPA fixes :");
for (int l=0; l < modify->n_post_force_respa; ++l) {
Fix *f = modify->fix[modify->list_post_force_respa[l]];
if (f->respa_level >= 0)
fprintf(screen," %d:%s[%s]",
MIN(f->respa_level+1,nlevels),f->style,f->id);
}
fprintf(screen,"\n");
timer->print_timeout(screen);
}
update->setupflag = 1;

10
src/USER-SMD/fix_smd_setvel.cpp
View File

@ -31,7 +31,6 @@
#include "modify.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
@ -192,9 +191,6 @@ void FixSMDSetVel::init() {
else
varflag = CONSTANT;
if (strstr(update->integrate_style, "respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
// cannot use non-zero forces for a minimization since no energy is integrated
// use fix addforce instead
@ -223,11 +219,7 @@ void FixSMDSetVel::setup(int vflag) {
if (strstr(update->integrate_style, "verlet"))
post_force(vflag);
else
for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
((Respa *) update->integrate)->copy_flevel_f(ilevel);
post_force_respa(vflag, ilevel, 0);
((Respa *) update->integrate)->copy_f_flevel(ilevel);
}
error->all(FLERR,"Fix smd/setvel does not support RESPA");
}
/* ---------------------------------------------------------------------- */

6
src/atom.cpp
View File

@ -94,8 +94,9 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
// USER-DPD
uCond = uMech = uChem = uCG = uCGnew = NULL;
duCond = duMech = duChem = NULL;
duChem = NULL;
dpdTheta = NULL;
ssaAIR = NULL;
// USER-SMD
@ -284,9 +285,8 @@ Atom::~Atom()
memory->destroy(uChem);
memory->destroy(uCG);
memory->destroy(uCGnew);
memory->destroy(duCond);
memory->destroy(duMech);
memory->destroy(duChem);
memory->destroy(ssaAIR);
memory->destroy(nspecial);
memory->destroy(special);

3
src/atom.h
View File

@ -88,9 +88,10 @@ class Atom : protected Pointers {
// USER-DPD package
double *uCond,*uMech,*uChem,*uCGnew,*uCG;
double *duCond,*duMech,*duChem;
double *duChem;
double *dpdTheta;
int nspecies_dpd;
int *ssaAIR; // Shardlow Splitting Algorithm Active Interaction Region number
// molecular info

8
src/balance.cpp
View File

@ -325,10 +325,6 @@ void Balance::command(int narg, char **arg)
bisection(1);
}
// output of final result
if (outflag) dumpout(update->ntimestep,fp);
// reset proc sub-domains
// for either brick or tiled comm style
@ -344,6 +340,10 @@ void Balance::command(int narg, char **arg)
delete irregular;
if (domain->triclinic) domain->lamda2x(atom->nlocal);
// output of final result
if (outflag) dumpout(update->ntimestep,fp);
// check if any atoms were lost
bigint natoms;

1
src/dump_image.cpp
View File

@ -893,7 +893,6 @@ void DumpImage::create_image()
int *molatom = atom->molatom;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
int newton_bond = force->newton_bond;
int molecular = atom->molecular;
Molecule **onemols = atom->avec->onemols;

10
src/fix.cpp
View File

@ -16,6 +16,7 @@
#include "fix.h"
#include "atom.h"
#include "group.h"
#include "force.h"
#include "atom_masks.h"
#include "memory.h"
#include "error.h"
@ -68,6 +69,8 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
dof_flag = 0;
special_alter_flag = 0;
enforce2d_flag = 0;
respa_level_support = 0;
respa_level = -1;
scalar_flag = vector_flag = array_flag = 0;
peratom_flag = local_flag = 0;
@ -130,6 +133,13 @@ void Fix::modify_params(int narg, char **arg)
else if (strcmp(arg[iarg+1],"yes") == 0) thermo_energy = 1;
else error->all(FLERR,"Illegal fix_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"respa") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
if (!respa_level_support) error->all(FLERR,"Illegal fix_modify command");
int lvl = force->inumeric(FLERR,arg[1]);
if (lvl < 0) error->all(FLERR,"Illegal fix_modify command");
respa_level = lvl-1;
iarg += 2;
} else {
int n = modify_param(narg-iarg,&arg[iarg]);
if (n == 0) error->all(FLERR,"Illegal fix_modify command");

2
src/fix.h
View File

@ -52,6 +52,8 @@ class Fix : protected Pointers {
int dof_flag; // 1 if has dof() method (not min_dof())
int special_alter_flag; // 1 if has special_alter() meth for spec lists
int enforce2d_flag; // 1 if has enforce2d method
int respa_level_support; // 1 if fix supports fix_modify respa
int respa_level; // which respa level to apply fix (1-Nrespa)
int scalar_flag; // 0/1 if compute_scalar() function exists
int vector_flag; // 0/1 if compute_vector() function exists

17
src/fix_addforce.cpp
View File

@ -47,6 +47,8 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
xstr = ystr = zstr = NULL;
@ -201,8 +203,13 @@ void FixAddForce::init()
update->whichflag == 2 && estyle == NONE)
error->all(FLERR,"Must use variable energy with fix addforce");
int max_respa = 0;
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
max_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0)
ilevel_respa = MIN(respa_level,max_respa);
}
/* ---------------------------------------------------------------------- */
@ -212,9 +219,9 @@ void FixAddForce::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -324,7 +331,7 @@ void FixAddForce::post_force(int vflag)
void FixAddForce::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_addforce.h
View File

@ -47,7 +47,7 @@ class FixAddForce : public Fix {
int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
double foriginal[4],foriginal_all[4];
int force_flag;
int nlevels_respa;
int ilevel_respa;
int maxatom;
double **sforce;

13
src/fix_aveforce.cpp
View File

@ -43,6 +43,8 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) :
size_vector = 3;
global_freq = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = nlevels_respa = 0;
xstr = ystr = zstr = NULL;
@ -161,8 +163,11 @@ void FixAveForce::init()
if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL) varflag = EQUAL;
else varflag = CONSTANT;
if (strstr(update->integrate_style,"respa"))
if (strstr(update->integrate_style,"respa")) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,nlevels_respa-1);
else ilevel_respa = nlevels_respa-1;
}
}
/* ---------------------------------------------------------------------- */
@ -255,10 +260,10 @@ void FixAveForce::post_force(int vflag)
void FixAveForce::post_force_respa(int vflag, int ilevel, int iloop)
{
// ave + extra force on outermost level
// just ave on inner levels
// ave + extra force on selected RESPA level
// just ave on all other levels
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
else {
Region *region = NULL;
if (iregion >= 0) {

2
src/fix_aveforce.h
View File

@ -45,7 +45,7 @@ class FixAveForce : public Fix {
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
int iregion;
double foriginal_all[4];
int nlevels_respa;
int nlevels_respa,ilevel_respa;
};
}

16
src/fix_drag.cpp
View File

@ -36,6 +36,8 @@ FixDrag::FixDrag(LAMMPS *lmp, int narg, char **arg) :
size_vector = 3;
global_freq = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
xflag = yflag = zflag = 1;
@ -67,8 +69,10 @@ int FixDrag::setmask()
void FixDrag::init()
{
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -78,9 +82,9 @@ void FixDrag::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -130,7 +134,7 @@ void FixDrag::post_force(int vflag)
void FixDrag::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ----------------------------------------------------------------------

2
src/fix_drag.h
View File

@ -39,7 +39,7 @@ class FixDrag : public Fix {
double f_mag;
int xflag,yflag,zflag;
double delta;
int nlevels_respa;
int ilevel_respa;
double ftotal[3],ftotal_all[3];
int force_flag;
};

16
src/fix_gravity.cpp
View File

@ -45,6 +45,8 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
respa_level_support = 1;
ilevel_respa = 0;
mstr = vstr = pstr = tstr = xstr = ystr = zstr = NULL;
mstyle = vstyle = pstyle = tstyle = xstyle = ystyle = zstyle = CONSTANT;
@ -162,8 +164,10 @@ int FixGravity::setmask()
void FixGravity::init()
{
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
// check variables
@ -234,9 +238,9 @@ void FixGravity::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -297,7 +301,7 @@ void FixGravity::post_force(int vflag)
void FixGravity::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_gravity.h
View File

@ -44,7 +44,7 @@ class FixGravity : public Fix {
double xdir,ydir,zdir;
double xgrav,ygrav,zgrav,xacc,yacc,zacc;
double degree2rad;
int nlevels_respa;
int ilevel_respa;
int time_origin;
int eflag;
double egrav,egrav_all;

16
src/fix_indent.cpp
View File

@ -50,6 +50,8 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
k = force->numeric(FLERR,arg[3]);
k3 = k/3.0;
@ -151,8 +153,10 @@ void FixIndent::init()
error->all(FLERR,"Variable for fix indent is not equal style");
}
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -162,9 +166,9 @@ void FixIndent::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -354,7 +358,7 @@ void FixIndent::post_force(int vflag)
void FixIndent::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_indent.h
View File

@ -47,7 +47,7 @@ class FixIndent : public Fix {
int indenter_flag,planeside;
double indenter[4],indenter_all[4];
int cdim,varflag;
int nlevels_respa;
int ilevel_respa;
void options(int, char **);
};

16
src/fix_restrain.cpp
View File

@ -51,6 +51,8 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
respa_level_support = 1;
ilevel_respa = 0;
// parse args
@ -147,8 +149,10 @@ int FixRestrain::setmask()
void FixRestrain::init()
{
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -158,9 +162,9 @@ void FixRestrain::setup(int vflag)
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -187,7 +191,7 @@ void FixRestrain::post_force(int vflag)
void FixRestrain::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_restrain.h
View File

@ -38,7 +38,7 @@ class FixRestrain : public Fix {
double compute_scalar();
private:
int nlevels_respa;
int ilevel_respa;
int nrestrain,maxrestrain;
int *rstyle;
int **ids;

12
src/fix_setforce.cpp
View File

@ -43,7 +43,8 @@ FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) :
size_vector = 3;
global_freq = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = nlevels_respa = 0;
xstr = ystr = zstr = NULL;
if (strstr(arg[3],"v_") == arg[3]) {
@ -172,8 +173,11 @@ void FixSetForce::init()
varflag = EQUAL;
else varflag = CONSTANT;
if (strstr(update->integrate_style,"respa"))
if (strstr(update->integrate_style,"respa")) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,nlevels_respa-1);
else ilevel_respa = nlevels_respa-1;
}
// cannot use non-zero forces for a minimization since no energy is integrated
// use fix addforce instead
@ -290,9 +294,9 @@ void FixSetForce::post_force(int vflag)
void FixSetForce::post_force_respa(int vflag, int ilevel, int iloop)
{
// set force to desired value on outermost level, 0.0 on other levels
// set force to desired value on requested level, 0.0 on other levels
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
else {
Region *region = NULL;
if (iregion >= 0) {

3
src/fix_setforce.h
View File

@ -36,6 +36,7 @@ class FixSetForce : public Fix {
void post_force_respa(int, int, int);
void min_post_force(int);
double compute_vector(int);
double memory_usage();
protected:
@ -46,7 +47,7 @@ class FixSetForce : public Fix {
int xvar,yvar,zvar,xstyle,ystyle,zstyle;
double foriginal[3],foriginal_all[3];
int force_flag;
int nlevels_respa;
int nlevels_respa,ilevel_respa;
int maxatom;
double **sforce;

16
src/fix_spring.cpp
View File

@ -48,6 +48,8 @@ FixSpring::FixSpring(LAMMPS *lmp, int narg, char **arg) :
extscalar = 1;
extvector = 1;
dynamic_group_allow = 1;
respa_level_support = 1;
ilevel_respa = 0;
group2 = NULL;
@ -130,8 +132,10 @@ void FixSpring::init()
masstotal = group->mass(igroup);
if (styleflag == COUPLE) masstotal2 = group->mass(igroup2);
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -141,9 +145,9 @@ void FixSpring::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -334,7 +338,7 @@ void FixSpring::spring_couple()
void FixSpring::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_spring.h
View File

@ -46,7 +46,7 @@ class FixSpring : public Fix {
char *group2;
int igroup2,group2bit;
double masstotal,masstotal2;
int nlevels_respa;
int ilevel_respa;
double espring,ftotal[4];
void spring_tether();

16
src/fix_spring_chunk.cpp
View File

@ -41,6 +41,8 @@ FixSpringChunk::FixSpringChunk(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
respa_level_support = 1;
ilevel_respa = 0;
k_spring = force->numeric(FLERR,arg[3]);
@ -114,8 +116,10 @@ void FixSpringChunk::init()
if (strcmp(idchunk,ccom->idchunk) != 0)
error->all(FLERR,"Fix spring chunk chunkID not same as comID chunkID");
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -125,9 +129,9 @@ void FixSpringChunk::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -229,7 +233,7 @@ void FixSpringChunk::post_force(int vflag)
void FixSpringChunk::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_spring_chunk.h
View File

@ -38,7 +38,7 @@ class FixSpringChunk : public Fix {
double compute_scalar();
private:
int nlevels_respa;
int ilevel_respa;
double k_spring;
double esprings;
char *idchunk,*idcom;

16
src/fix_spring_rg.cpp
View File

@ -44,6 +44,8 @@ FixSpringRG::FixSpringRG(LAMMPS *lmp, int narg, char **arg) :
else rg0 = force->numeric(FLERR,arg[4]);
dynamic_group_allow = 1;
respa_level_support = 1;
ilevel_respa = 0;
}
/* ---------------------------------------------------------------------- */
@ -72,8 +74,10 @@ void FixSpringRG::init()
rg0_flag = 0;
}
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -83,9 +87,9 @@ void FixSpringRG::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -140,5 +144,5 @@ void FixSpringRG::post_force(int vflag)
void FixSpringRG::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}

2
src/fix_spring_rg.h
View File

@ -34,7 +34,7 @@ class FixSpringRG : public Fix {
void post_force_respa(int, int, int);
private:
int nlevels_respa,rg0_flag;
int ilevel_respa,rg0_flag;
double rg0,k,masstotal;
};

15
src/fix_spring_self.cpp
View File

@ -41,6 +41,7 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
respa_level_support = 1;
k = force->numeric(FLERR,arg[3]);
if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command");
@ -118,8 +119,10 @@ int FixSpringSelf::setmask()
void FixSpringSelf::init()
{
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -129,9 +132,9 @@ void FixSpringSelf::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -179,7 +182,7 @@ void FixSpringSelf::post_force(int vflag)
void FixSpringSelf::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_spring_self.h
View File

@ -51,7 +51,7 @@ class FixSpringSelf : public Fix {
double k,espring;
double **xoriginal; // original coords of atoms
int xflag, yflag, zflag;
int nlevels_respa;
int ilevel_respa;
};
}

19
src/fix_viscous.cpp
View File

@ -49,6 +49,9 @@ FixViscous::FixViscous(LAMMPS *lmp, int narg, char **arg) :
iarg += 3;
} else error->all(FLERR,"Illegal fix viscous command");
}
respa_level_support = 1;
ilevel_respa = 0;
}
/* ---------------------------------------------------------------------- */
@ -73,8 +76,12 @@ int FixViscous::setmask()
void FixViscous::init()
{
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
int max_respa = 0;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = max_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,max_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -84,9 +91,9 @@ void FixViscous::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -126,7 +133,7 @@ void FixViscous::post_force(int vflag)
void FixViscous::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_viscous.h
View File

@ -38,7 +38,7 @@ class FixViscous : public Fix {
protected:
double *gamma;
int nlevels_respa;
int ilevel_respa;
};
}

18
src/fix_wall.cpp
View File

@ -42,6 +42,8 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
// parse args
@ -227,8 +229,6 @@ int FixWall::setmask()
void FixWall::init()
{
dt = update->dt;
for (int m = 0; m < nwall; m++) {
if (xstyle[m] == VARIABLE) {
xindex[m] = input->variable->find(xstr[m]);
@ -257,8 +257,10 @@ void FixWall::init()
for (int m = 0; m < nwall; m++) precompute(m);
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -268,9 +270,9 @@ void FixWall::setup(int vflag)
if (strstr(update->integrate_style,"verlet")) {
if (!fldflag) post_force(vflag);
} else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -335,7 +337,7 @@ void FixWall::post_force(int vflag)
void FixWall::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

3
src/fix_wall.h
View File

@ -53,8 +53,7 @@ class FixWall : public Fix {
char *estr[6],*sstr[6];
int varflag; // 1 if any wall position,epsilon,sigma is a var
int eflag; // per-wall flag for energy summation
int nlevels_respa;
double dt;
int ilevel_respa;
int fldflag;
};

16
src/fix_wall_region.cpp
View File

@ -44,6 +44,8 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extscalar = 1;
extvector = 1;
respa_level_support = 1;
ilevel_respa = 0;
// parse args
@ -151,8 +153,10 @@ void FixWallRegion::init()
offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
}
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
@ -162,9 +166,9 @@ void FixWallRegion::setup(int vflag)
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -245,7 +249,7 @@ void FixWallRegion::post_force(int vflag)
void FixWallRegion::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */

2
src/fix_wall_region.h
View File

@ -43,7 +43,7 @@ class FixWallRegion : public Fix {
double epsilon,sigma,cutoff;
int eflag;
double ewall[4],ewall_all[4];
int nlevels_respa;
int ilevel_respa;
char *idregion;
double coeff1,coeff2,coeff3,coeff4,offset;

2
src/improper.cpp
View File

@ -68,6 +68,8 @@ void Improper::init()
error->all(FLERR,"Improper coeffs are not set");
for (int i = 1; i <= atom->nimpropertypes; i++)
if (setflag[i] == 0) error->all(FLERR,"All improper coeffs are not set");
init_style();
}
/* ----------------------------------------------------------------------

1
src/improper.h
View File

@ -40,6 +40,7 @@ class Improper : protected Pointers {
Improper(class LAMMPS *);
virtual ~Improper();
virtual void init();
virtual void init_style() {}
virtual void compute(int, int) = 0;
virtual void settings(int, char **) {}
virtual void coeff(int, char **) = 0;

9
src/improper_hybrid.cpp
View File

@ -156,6 +156,15 @@ void ImproperHybrid::allocate()
for (int m = 0; m < nstyles; m++) improperlist[m] = NULL;
}
/* ---------------------------------------------------------------------- */
void ImproperHybrid::init_style()
{
for (int i = 0; i < nstyles; i++)
styles[i]->init_style();
}
/* ----------------------------------------------------------------------
create one improper style for each arg in list
------------------------------------------------------------------------- */

1
src/improper_hybrid.h
View File

@ -33,6 +33,7 @@ class ImproperHybrid : public Improper {
ImproperHybrid(class LAMMPS *);
~ImproperHybrid();
void init_style();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);

8
src/min.cpp
View File

@ -185,7 +185,8 @@ void Min::setup()
if (comm->me == 0 && screen) {
fprintf(screen,"Setting up %s style minimization ...\n",
update->minimize_style);
fprintf(screen," Unit style: %s\n", update->unit_style);
fprintf(screen," Unit style : %s\n", update->unit_style);
timer->print_timeout(screen);
}
update->setupflag = 1;
@ -397,7 +398,7 @@ void Min::run(int n)
// add ntimestep to all computes that store invocation times
// since are hardwiring call to thermo/dumps and computes may not be ready
if (stop_condition) {
if (stop_condition != MAXITER) {
update->nsteps = niter;
if (update->restrict_output == 0) {
@ -798,6 +799,7 @@ char *Min::stopstrings(int n)
"forces are zero",
"quadratic factors are zero",
"trust region too small",
"HFTN minimizer error"};
"HFTN minimizer error",
"walltime limit reached"};
return (char *) strings[n];
}

4
src/min.h
View File

@ -46,6 +46,10 @@ class Min : protected Pointers {
virtual void reset_vectors() = 0;
virtual int iterate(int) = 0;
// possible return values of iterate() method
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,
ZEROQUAD,TRSMALL,INTERROR,TIMEOUT};
protected:
int eflag,vflag; // flags for energy/virial computation
int virial_style; // compute virial explicitly or implicitly

7
src/min_cg.cpp
View File

@ -26,10 +26,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
/* ---------------------------------------------------------------------- */
MinCG::MinCG(LAMMPS *lmp) : MinLineSearch(lmp) {}
@ -67,6 +63,9 @@ int MinCG::iterate(int maxiter)
for (int iter = 0; iter < maxiter; iter++) {
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;

7
src/min_fire.cpp
View File

@ -27,10 +27,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
#define DELAYSTEP 5
#define DT_GROW 1.1
#define DT_SHRINK 0.5
@ -93,6 +89,9 @@ int MinFire::iterate(int maxiter)
for (int iter = 0; iter < maxiter; iter++) {
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;

16
src/min_hftn.cpp
View File

@ -42,12 +42,12 @@ using namespace LAMMPS_NS;
//---- CONSTANTS MAP TO stopstrings DECLARED IN Min.run (min.cpp).
static const int STOP_MAX_ITER = 0; //-- MAX ITERATIONS EXCEEDED
static const int STOP_MAX_FORCE_EVALS = 1; //-- MAX FORCE EVALUATIONS EXCEEDED
static const int STOP_ENERGY_TOL = 2; //-- STEP DID NOT CHANGE ENERGY
static const int STOP_FORCE_TOL = 3; //-- CONVERGED TO DESIRED FORCE TOL
static const int STOP_TR_TOO_SMALL = 8; //-- TRUST REGION TOO SMALL
static const int STOP_ERROR = 9; //-- INTERNAL ERROR
static const int STOP_MAX_ITER = Min::MAXITER; //-- MAX ITERATIONS EXCEEDED
static const int STOP_MAX_FORCE_EVALS = Min::MAXEVAL; //-- MAX FORCE EVALUATIONS EXCEEDED
static const int STOP_ENERGY_TOL = Min::ETOL; //-- STEP DID NOT CHANGE ENERGY
static const int STOP_FORCE_TOL = Min::FTOL; //-- CONVERGED TO DESIRED FORCE TOL
static const int STOP_TR_TOO_SMALL = Min::TRSMALL; //-- TRUST REGION TOO SMALL
static const int STOP_ERROR = Min::INTERROR; //-- INTERNAL ERROR
static const int NO_CGSTEP_BECAUSE_F_TOL_SATISFIED = 0;
static const int CGSTEP_NEWTON = 1;
@ -292,6 +292,10 @@ int MinHFTN::execute_hftn_(const bool bPrintProgress,
double dCurrentEnergy = dInitialEnergy;
double dCurrentForce2 = dInitialForce2;
for (niter = 0; niter < update->nsteps; niter++) {
if (timer->check_timeout(niter))
return(Min::TIMEOUT);
(update->ntimestep)++;
//---- CALL THE INNER LOOP TO GET THE NEXT TRUST REGION STEP.

4
src/min_linesearch.cpp
View File

@ -52,10 +52,6 @@ using namespace LAMMPS_NS;
#define EMACH 1.0e-8
#define EPS_QUAD 1.0e-28
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
/* ---------------------------------------------------------------------- */
MinLineSearch::MinLineSearch(LAMMPS *lmp) : Min(lmp)

7
src/min_quickmin.cpp
View File

@ -28,10 +28,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
#define DELAYSTEP 5
/* ---------------------------------------------------------------------- */
@ -88,6 +84,9 @@ int MinQuickMin::iterate(int maxiter)
for (int iter = 0; iter < maxiter; iter++) {
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;

7
src/min_sd.cpp
View File

@ -25,10 +25,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
/* ---------------------------------------------------------------------- */
MinSD::MinSD(LAMMPS *lmp) : MinLineSearch(lmp) {}
@ -58,6 +54,9 @@ int MinSD::iterate(int maxiter)
for (int iter = 0; iter < maxiter; iter++) {
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;

Some files were not shown because too many files have changed in this diff Show More