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

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

Browse Source
This commit is contained in:
sjplimp
2008-04-10 20:53:37 +00:00
parent b729a79670
commit 943216c548
6 changed files with 44 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
src/COLLOID/pair_lubricate.cpp
View File

@ -86,6 +86,7 @@ void PairLubricate::compute(int eflag, int vflag)
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int newton_pair = force->newton_pair; int newton_pair = force->newton_pair;
int omega_flag = atom->omega_flag; int omega_flag = atom->omega_flag;
double vxmu2f = force->vxmu2f;
inum = list->inum; inum = list->inum;
ilist = list->ilist; ilist = list->ilist;
@ -198,16 +199,17 @@ void PairLubricate::compute(int eflag, int vflag)
h_sep = r - 2.0*radi; h_sep = r - 2.0*radi;
if (flag1) if (flag1)
a_squeeze = (3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep); a_squeeze = vxmu2f *
(3.0*PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep);
if (flag2) if (flag2)
a_shear = (PI*mu*2.*radi/2.0) * a_shear = vxmu2f * (PI*mu*2.*radi/2.0) *
log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0; log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0;
if (flag3) if (flag3)
a_pump = (PI*mu*pow(2.0*radi,4)/8.0) * a_pump = vxmu2f * (PI*mu*pow(2.0*radi,4)/8.0) *
((3.0/20.0) * log(2.0*radi/2.0/h_sep) + ((3.0/20.0) * log(2.0*radi/2.0/h_sep) +
(63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep)); (63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep));
if (flag4) if (flag4)
a_twist = (PI*mu*pow(2.0*radi,4)/4.0) * a_twist = vxmu2f * (PI*mu*pow(2.0*radi,4)/4.0) *
(h_sep/2.0/radi) * log(2.0/(2.0*h_sep)); (h_sep/2.0/radi) * log(2.0/(2.0*h_sep));
if (h_sep >= cut_inner[itype][jtype]) { if (h_sep >= cut_inner[itype][jtype]) {

32
src/fix_nph.cpp
View File

@ -37,6 +37,8 @@ using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B) #define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B) #define MAX(A,B) ((A) > (B)) ? (A) : (B)
enum{XYZ,XY,YZ,XZ,ANISO};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
@ -55,7 +57,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[3],"xyz") == 0) { if (strcmp(arg[3],"xyz") == 0) {
if (narg < 7) error->all("Illegal fix nph command"); if (narg < 7) error->all("Illegal fix nph command");
press_couple = 0; press_couple = XYZ;
p_start[0] = p_start[1] = p_start[2] = atof(arg[4]); p_start[0] = p_start[1] = p_start[2] = atof(arg[4]);
p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[5]); p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[5]);
p_period[0] = p_period[1] = p_period[2] = atof(arg[6]); p_period[0] = p_period[1] = p_period[2] = atof(arg[6]);
@ -66,16 +68,16 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
} }
} else { } else {
if (strcmp(arg[3],"xy") == 0) press_couple = 1; if (strcmp(arg[3],"xy") == 0) press_couple = XY;
else if (strcmp(arg[3],"yz") == 0) press_couple = 2; else if (strcmp(arg[3],"yz") == 0) press_couple = YZ;
else if (strcmp(arg[3],"xz") == 0) press_couple = 3; else if (strcmp(arg[3],"xz") == 0) press_couple = XZ;
else if (strcmp(arg[3],"aniso") == 0) press_couple = 4; else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO;
else error->all("Illegal fix nph command"); else error->all("Illegal fix nph command");
if (narg < 11) error->all("Illegal fix nph command"); if (narg < 11) error->all("Illegal fix nph command");
if (domain->dimension == 2 && if (domain->dimension == 2 &&
(press_couple == 1 || press_couple == 2 || press_couple == 3)) (press_couple == XY || press_couple == YZ || press_couple == XZ))
error->all("Invalid fix nph command for a 2d simulation"); error->all("Invalid fix nph command for a 2d simulation");
if (strcmp(arg[4],"NULL") == 0) { if (strcmp(arg[4],"NULL") == 0) {
@ -117,7 +119,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
allremap = 1; allremap = 1;
int iarg; int iarg;
if (press_couple == 0) iarg = 7; if (press_couple == XYZ) iarg = 7;
else iarg = 11; else iarg = 11;
while (iarg < narg) { while (iarg < narg) {
@ -134,7 +136,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
} else error->all("Illegal fix nph command"); } else error->all("Illegal fix nph command");
} }
// error checks // check for periodicity in controlled dimensions
if (p_flag[0] && domain->xperiodic == 0) if (p_flag[0] && domain->xperiodic == 0)
error->all("Cannot use fix nph on a non-periodic dimension"); error->all("Cannot use fix nph on a non-periodic dimension");
@ -313,7 +315,7 @@ void FixNPH::setup(int vflag)
p_target[1] = p_start[1]; p_target[1] = p_start[1];
p_target[2] = p_start[2]; p_target[2] = p_start[2];
if (press_couple == 0) { if (press_couple == XYZ) {
double tmp = temperature->compute_scalar(); double tmp = temperature->compute_scalar();
tmp = pressure->compute_scalar(); tmp = pressure->compute_scalar();
} else { } else {
@ -426,7 +428,7 @@ void FixNPH::final_integrate()
// compute new pressure // compute new pressure
if (press_couple == 0) { if (press_couple == XYZ) {
double tmp = temperature->compute_scalar(); double tmp = temperature->compute_scalar();
tmp = pressure->compute_scalar(); tmp = pressure->compute_scalar();
} else { } else {
@ -598,21 +600,21 @@ void FixNPH::couple()
{ {
double *tensor = pressure->vector; double *tensor = pressure->vector;
if (press_couple == 0) if (press_couple == XYZ)
p_current[0] = p_current[1] = p_current[2] = pressure->scalar; p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
else if (press_couple == 1) { else if (press_couple == XY) {
double ave = 0.5 * (tensor[0] + tensor[1]); double ave = 0.5 * (tensor[0] + tensor[1]);
p_current[0] = p_current[1] = ave; p_current[0] = p_current[1] = ave;
p_current[2] = tensor[2]; p_current[2] = tensor[2];
} else if (press_couple == 2) { } else if (press_couple == YZ) {
double ave = 0.5 * (tensor[1] + tensor[2]); double ave = 0.5 * (tensor[1] + tensor[2]);
p_current[1] = p_current[2] = ave; p_current[1] = p_current[2] = ave;
p_current[0] = tensor[0]; p_current[0] = tensor[0];
} else if (press_couple == 3) { } else if (press_couple == XZ) {
double ave = 0.5 * (tensor[0] + tensor[2]); double ave = 0.5 * (tensor[0] + tensor[2]);
p_current[0] = p_current[2] = ave; p_current[0] = p_current[2] = ave;
p_current[1] = tensor[1]; p_current[1] = tensor[1];
} if (press_couple == 4) { } else if (press_couple == ANISO) {
p_current[0] = tensor[0]; p_current[0] = tensor[0];
p_current[1] = tensor[1]; p_current[1] = tensor[1];
p_current[2] = tensor[2]; p_current[2] = tensor[2];

31
src/fix_npt.cpp
View File

@ -37,6 +37,7 @@ using namespace LAMMPS_NS;
#define MAX(A,B) ((A) > (B)) ? (A) : (B) #define MAX(A,B) ((A) > (B)) ? (A) : (B)
enum{NOBIAS,BIAS}; enum{NOBIAS,BIAS};
enum{XYZ,XY,YZ,XZ,ANISO};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -63,7 +64,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[6],"xyz") == 0) { if (strcmp(arg[6],"xyz") == 0) {
if (narg < 10) error->all("Illegal fix npt command"); if (narg < 10) error->all("Illegal fix npt command");
press_couple = 0; press_couple = XYZ;
p_start[0] = p_start[1] = p_start[2] = atof(arg[7]); p_start[0] = p_start[1] = p_start[2] = atof(arg[7]);
p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[8]); p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[8]);
p_period[0] = p_period[1] = p_period[2] = atof(arg[9]); p_period[0] = p_period[1] = p_period[2] = atof(arg[9]);
@ -74,16 +75,16 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
} }
} else { } else {
if (strcmp(arg[6],"xy") == 0) press_couple = 1; if (strcmp(arg[6],"xy") == 0) press_couple = XY;
else if (strcmp(arg[6],"yz") == 0) press_couple = 2; else if (strcmp(arg[6],"yz") == 0) press_couple = YZ;
else if (strcmp(arg[6],"xz") == 0) press_couple = 3; else if (strcmp(arg[6],"xz") == 0) press_couple = XZ;
else if (strcmp(arg[6],"aniso") == 0) press_couple = 4; else if (strcmp(arg[6],"aniso") == 0) press_couple = ANISO;
else error->all("Illegal fix npt command"); else error->all("Illegal fix npt command");
if (narg < 14) error->all("Illegal fix npt command"); if (narg < 14) error->all("Illegal fix npt command");
if (domain->dimension == 2 && if (domain->dimension == 2 &&
(press_couple == 1 || press_couple == 2 || press_couple == 3)) (press_couple == XY || press_couple == YZ || press_couple == XZ))
error->all("Invalid fix npt command for a 2d simulation"); error->all("Invalid fix npt command for a 2d simulation");
if (strcmp(arg[7],"NULL") == 0) { if (strcmp(arg[7],"NULL") == 0) {
@ -125,7 +126,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
allremap = 1; allremap = 1;
int iarg; int iarg;
if (press_couple == 0) iarg = 10; if (press_couple == XYZ) iarg = 10;
else iarg = 14; else iarg = 14;
while (iarg < narg) { while (iarg < narg) {
@ -142,7 +143,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
} else error->all("Illegal fix npt command"); } else error->all("Illegal fix npt command");
} }
// error checks // check for periodicity in controlled dimensions
if (p_flag[0] && domain->xperiodic == 0) if (p_flag[0] && domain->xperiodic == 0)
error->all("Cannot use fix npt on a non-periodic dimension"); error->all("Cannot use fix npt on a non-periodic dimension");
@ -322,7 +323,7 @@ void FixNPT::setup(int vflag)
p_target[2] = p_start[2]; p_target[2] = p_start[2];
t_current = temperature->compute_scalar(); t_current = temperature->compute_scalar();
if (press_couple == 0) { if (press_couple == XYZ) {
double tmp = pressure->compute_scalar(); double tmp = pressure->compute_scalar();
} else { } else {
temperature->compute_vector(); temperature->compute_vector();
@ -473,7 +474,7 @@ void FixNPT::final_integrate()
// compute new T,P // compute new T,P
t_current = temperature->compute_scalar(); t_current = temperature->compute_scalar();
if (press_couple == 0) { if (press_couple == XYZ) {
double tmp = pressure->compute_scalar(); double tmp = pressure->compute_scalar();
} else { } else {
temperature->compute_vector(); temperature->compute_vector();
@ -687,21 +688,21 @@ void FixNPT::couple()
{ {
double *tensor = pressure->vector; double *tensor = pressure->vector;
if (press_couple == 0) if (press_couple == XYZ)
p_current[0] = p_current[1] = p_current[2] = pressure->scalar; p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
else if (press_couple == 1) { else if (press_couple == XY) {
double ave = 0.5 * (tensor[0] + tensor[1]); double ave = 0.5 * (tensor[0] + tensor[1]);
p_current[0] = p_current[1] = ave; p_current[0] = p_current[1] = ave;
p_current[2] = tensor[2]; p_current[2] = tensor[2];
} else if (press_couple == 2) { } else if (press_couple == YZ) {
double ave = 0.5 * (tensor[1] + tensor[2]); double ave = 0.5 * (tensor[1] + tensor[2]);
p_current[1] = p_current[2] = ave; p_current[1] = p_current[2] = ave;
p_current[0] = tensor[0]; p_current[0] = tensor[0];
} else if (press_couple == 3) { } else if (press_couple == XZ) {
double ave = 0.5 * (tensor[0] + tensor[2]); double ave = 0.5 * (tensor[0] + tensor[2]);
p_current[0] = p_current[2] = ave; p_current[0] = p_current[2] = ave;
p_current[1] = tensor[1]; p_current[1] = tensor[1];
} if (press_couple == 4) { } else if (press_couple == ANISO) {
p_current[0] = tensor[0]; p_current[0] = tensor[0];
p_current[1] = tensor[1]; p_current[1] = tensor[1];
p_current[2] = tensor[2]; p_current[2] = tensor[2];

1
src/force.h
View File

@ -26,6 +26,7 @@ class Force : protected Pointers {
double nktv2p; // conversion of NkT/V to pressure double nktv2p; // conversion of NkT/V to pressure
double qqr2e; // conversion of q^2/r to energy double qqr2e; // conversion of q^2/r to energy
double qe2f; // conversion of qE to force double qe2f; // conversion of qE to force
double vxmu2f; // conversion of vx mu to force
double dielectric; // dielectric constant double dielectric; // dielectric constant
double qqrd2e; // q^2/r to energy w/ dielectric constant double qqrd2e; // q^2/r to energy w/ dielectric constant

1
src/min_cg.cpp
View File

@ -186,6 +186,7 @@ void MinCG::run()
eng_force(&ntmp,&xtmp,&htmp,&etmp); eng_force(&ntmp,&xtmp,&htmp,&etmp);
output->write(update->ntimestep); output->write(update->ntimestep);
} }
timer->barrier_stop(TIME_LOOP); timer->barrier_stop(TIME_LOOP);
// delete fix at end of run, so its atom arrays won't persist // delete fix at end of run, so its atom arrays won't persist

3
src/update.cpp
View File

@ -102,6 +102,7 @@ void Update::set_units(const char *style)
force->nktv2p = 1.0; force->nktv2p = 1.0;
force->qqr2e = 1.0; force->qqr2e = 1.0;
force->qe2f = 1.0; force->qe2f = 1.0;
force->vxmu2f = 1.0;
dt = 0.005; dt = 0.005;
neighbor->skin = 0.3; neighbor->skin = 0.3;
@ -112,6 +113,7 @@ void Update::set_units(const char *style)
force->nktv2p = 68568.415; force->nktv2p = 68568.415;
force->qqr2e = 332.06371; force->qqr2e = 332.06371;
force->qe2f = 23.060549; force->qe2f = 23.060549;
force->vxmu2f = 1.4393264316e4;
dt = 1.0; dt = 1.0;
neighbor->skin = 2.0; neighbor->skin = 2.0;
@ -122,6 +124,7 @@ void Update::set_units(const char *style)
force->nktv2p = 1.6021765e6; force->nktv2p = 1.6021765e6;
force->qqr2e = 14.399645; force->qqr2e = 14.399645;
force->qe2f = 1.0; force->qe2f = 1.0;
force->vxmu2f = 0.6241509647;
dt = 0.001; dt = 0.001;
neighbor->skin = 2.0; neighbor->skin = 2.0;
} else error->all("Illegal units command"); } else error->all("Illegal units command");