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@11515 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2014-02-10 20:44:43 +00:00
parent 63c9df46a9
commit 92d0ee3c32
5 changed files with 164 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/Makefile.lib
View File

File diff suppressed because one or more lines are too long

39
src/compute_pressure.cpp
View File

@ -49,16 +49,19 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
// store temperature ID used by pressure computation
// insure it is valid for temperature computation
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
if (strcmp(arg[3],"NULL") == 0) id_temp = NULL;
else {
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find compute pressure temperature ID");
if (modify->compute[icompute]->tempflag == 0)
error->all(FLERR,
"Compute pressure temperature ID does not compute temperature");
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find compute pressure temperature ID");
if (modify->compute[icompute]->tempflag == 0)
error->all(FLERR,
"Compute pressure temperature ID does not compute temperature");
}
// process optional args
@ -91,6 +94,12 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
}
}
// error check
if (keflag && id_temp == NULL)
error->all(FLERR,"Compute pressure requires temperature ID "
"to include kinetic energy");
vector = new double[6];
nvirial = 0;
vptr = NULL;
@ -116,10 +125,12 @@ void ComputePressure::init()
// set temperature compute, must be done in init()
// fixes could have changed or compute_modify could have changed it
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find compute pressure temperature ID");
temperature = modify->compute[icompute];
if (keflag) {
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find compute pressure temperature ID");
temperature = modify->compute[icompute];
}
// detect contributions to virial
// vptr points to all virial[6] contributions
@ -169,7 +180,7 @@ double ComputePressure::compute_scalar()
if (update->vflag_global != invoked_scalar)
error->all(FLERR,"Virial was not tallied on needed timestep");
// invoke temperature it it hasn't been already
// invoke temperature if it hasn't been already
double t;
if (keflag) {

127
src/compute_stress_atom.cpp
View File

@ -31,12 +31,14 @@
using namespace LAMMPS_NS;
enum{NOBIAS,BIAS};
/* ---------------------------------------------------------------------- */
ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg < 3) error->all(FLERR,"Illegal compute stress/atom command");
if (narg < 4) error->all(FLERR,"Illegal compute stress/atom command");
peratom_flag = 1;
size_peratom_cols = 6;
@ -44,7 +46,26 @@ ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) :
timeflag = 1;
comm_reverse = 6;
if (narg == 3) {
// store temperature ID used by stress computation
// insure it is valid for temperature computation
if (strcmp(arg[3],"NULL") == 0) id_temp = NULL;
else {
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find compute stress/atom temperature ID");
if (modify->compute[icompute]->tempflag == 0)
error->all(FLERR,
"Compute stress/atom temperature ID does not compute temperature");
}
// process optional args
if (narg == 4) {
keflag = 1;
pairflag = 1;
bondflag = angleflag = dihedralflag = improperflag = 1;
@ -56,7 +77,7 @@ ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) :
bondflag = angleflag = dihedralflag = improperflag = 0;
kspaceflag = 0;
fixflag = 0;
int iarg = 3;
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"ke") == 0) keflag = 1;
else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
@ -83,11 +104,29 @@ ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) :
ComputeStressAtom::~ComputeStressAtom()
{
delete [] id_temp;
memory->destroy(stress);
}
/* ---------------------------------------------------------------------- */
void ComputeStressAtom::init()
{
// set temperature compute, must be done in init()
// fixes could have changed or compute_modify could have changed it
if (id_temp) {
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
error->all(FLERR,"Could not find compute stress/atom temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempbias) biasflag = BIAS;
else biasflag = NOBIAS;
} else biasflag = NOBIAS;
}
/* ---------------------------------------------------------------------- */
void ComputeStressAtom::compute_peratom()
{
int i,j;
@ -211,6 +250,7 @@ void ComputeStressAtom::compute_peratom()
}
// include kinetic energy term for each atom in group
// apply temperature bias is applicable
// mvv2e converts mv^2 to energy
if (keflag) {
@ -220,29 +260,68 @@ void ComputeStressAtom::compute_peratom()
int *type = atom->type;
double mvv2e = force->mvv2e;
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
onemass = mvv2e * rmass[i];
stress[i][0] += onemass*v[i][0]*v[i][0];
stress[i][1] += onemass*v[i][1]*v[i][1];
stress[i][2] += onemass*v[i][2]*v[i][2];
stress[i][3] += onemass*v[i][0]*v[i][1];
stress[i][4] += onemass*v[i][0]*v[i][2];
stress[i][5] += onemass*v[i][1]*v[i][2];
}
if (biasflag == NOBIAS) {
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
onemass = mvv2e * rmass[i];
stress[i][0] += onemass*v[i][0]*v[i][0];
stress[i][1] += onemass*v[i][1]*v[i][1];
stress[i][2] += onemass*v[i][2]*v[i][2];
stress[i][3] += onemass*v[i][0]*v[i][1];
stress[i][4] += onemass*v[i][0]*v[i][2];
stress[i][5] += onemass*v[i][1]*v[i][2];
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
onemass = mvv2e * mass[type[i]];
stress[i][0] += onemass*v[i][0]*v[i][0];
stress[i][1] += onemass*v[i][1]*v[i][1];
stress[i][2] += onemass*v[i][2]*v[i][2];
stress[i][3] += onemass*v[i][0]*v[i][1];
stress[i][4] += onemass*v[i][0]*v[i][2];
stress[i][5] += onemass*v[i][1]*v[i][2];
}
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
onemass = mvv2e * mass[type[i]];
stress[i][0] += onemass*v[i][0]*v[i][0];
stress[i][1] += onemass*v[i][1]*v[i][1];
stress[i][2] += onemass*v[i][2]*v[i][2];
stress[i][3] += onemass*v[i][0]*v[i][1];
stress[i][4] += onemass*v[i][0]*v[i][2];
stress[i][5] += onemass*v[i][1]*v[i][2];
}
// invoke temperature if it hasn't been already
// this insures bias factor is pre-computed
if (keflag && temperature->invoked_scalar != update->ntimestep)
temperature->compute_scalar();
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
onemass = mvv2e * rmass[i];
stress[i][0] += onemass*v[i][0]*v[i][0];
stress[i][1] += onemass*v[i][1]*v[i][1];
stress[i][2] += onemass*v[i][2]*v[i][2];
stress[i][3] += onemass*v[i][0]*v[i][1];
stress[i][4] += onemass*v[i][0]*v[i][2];
stress[i][5] += onemass*v[i][1]*v[i][2];
temperature->restore_bias(i,v[i]);
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
onemass = mvv2e * mass[type[i]];
stress[i][0] += onemass*v[i][0]*v[i][0];
stress[i][1] += onemass*v[i][1]*v[i][1];
stress[i][2] += onemass*v[i][2]*v[i][2];
stress[i][3] += onemass*v[i][0]*v[i][1];
stress[i][4] += onemass*v[i][0]*v[i][2];
stress[i][5] += onemass*v[i][1]*v[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
}

7
src/compute_stress_atom.h
View File

@ -28,7 +28,7 @@ class ComputeStressAtom : public Compute {
public:
ComputeStressAtom(class LAMMPS *, int, char **);
~ComputeStressAtom();
void init() {}
void init();
void compute_peratom();
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
@ -36,7 +36,10 @@ class ComputeStressAtom : public Compute {
private:
int keflag,pairflag,bondflag,angleflag,dihedralflag,improperflag;
int kspaceflag,fixflag;
int kspaceflag,fixflag,biasflag;
Compute *temperature;
char *id_temp;
int nmax;
double **stress;
};

34
src/pair_hybrid.cpp
View File

@ -100,9 +100,10 @@ void PairHybrid::compute(int eflag, int vflag)
for (m = 0; m < nstyles; m++) {
// invoke compute() unless
// invoke compute() unless compute flag is turned off or
// outerflag is set and sub-style has a compute_outer() method
if (styles[m]->compute_flag == 0) continue;
if (outerflag && styles[m]->respa_enable)
styles[m]->compute_outer(eflag,vflag_substyle);
else styles[m]->compute(eflag,vflag_substyle);
@ -280,8 +281,10 @@ void PairHybrid::flags()
// ghostneigh = 1 if any sub-style is set
// ewaldflag, pppmflag, msmflag, dispersionflag, tip4pflag = 1
// if any sub-style is set
// compute_flag = 1 if any sub-style is set
single_enable = 0;
compute_flag = 0;
for (m = 0; m < nstyles; m++) {
if (styles[m]->single_enable) single_enable = 1;
if (styles[m]->respa_enable) respa_enable = 1;
@ -293,6 +296,7 @@ void PairHybrid::flags()
if (styles[m]->msmflag) msmflag = 1;
if (styles[m]->dispersionflag) dispersionflag = 1;
if (styles[m]->tip4pflag) tip4pflag = 1;
if (styles[m]->compute_flag) compute_flag = 1;
}
// single_extra = min of all sub-style single_extra
@ -694,14 +698,34 @@ double PairHybrid::single(int i, int j, int itype, int jtype,
/* ----------------------------------------------------------------------
modify parameters of the pair style
call modify_params of PairHybrid
also pass command args to each sub-style of hybrid
if 1st keyword is pair, then applies to one sub-style
else pass command args to every sub-style of hybrid
------------------------------------------------------------------------- */
void PairHybrid::modify_params(int narg, char **arg)
{
Pair::modify_params(narg,arg);
for (int m = 0; m < nstyles; m++) styles[m]->modify_params(narg,arg);
if (narg == 0) error->all(FLERR,"Illegal pair_modify command");
if (strcmp(arg[0],"pair") == 0) {
if (narg < 2) error->all(FLERR,"Illegal pair_modify command");
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
if (m == nstyles) error->all(FLERR,"Unknown pair_modify hybrid sub-style");
if (multiple[m] == 0)
styles[m]->modify_params(narg-2,&arg[2]);
else {
if (narg < 3) error->all(FLERR,"Illegal pair_modify command");
int multiflag = force->inumeric(FLERR,arg[2]);
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0 && multiflag == multiple[m]) break;
if (m == nstyles)
error->all(FLERR,"Unknown pair_modify hybrid sub-style");
styles[m]->modify_params(narg-3,&arg[3]);
}
} else
for (int m = 0; m < nstyles; m++) styles[m]->modify_params(narg,arg);
}
/* ----------------------------------------------------------------------