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

This commit is contained in:
sjplimp
2015-03-31 22:02:12 +00:00
parent 8428bd3d9d
commit 536ae8000d
13 changed files with 54 additions and 49 deletions

View File

@ -179,6 +179,8 @@ void ComputeTempAsphere::dof_compute()
}
dof -= extra_dof + fix_dof;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -267,8 +269,6 @@ double ComputeTempAsphere::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic || tempbias == 2) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -214,6 +214,8 @@ void ComputeTempCS::dof_compute()
dof = nper * natoms;
dof -= nper * nshells;
dof -= extra_dof + fix_dof;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -256,8 +258,6 @@ double ComputeTempCS::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -508,7 +508,10 @@ void FixRigidSmall::init()
/* ----------------------------------------------------------------------
setup static/dynamic properties of rigid bodies, using current atom info
allows resetting of atom properties like mass between runs
only do static initialization once if using an infile
only do static initialization must be done once
do not do again if using an infile or mol keyword
this is b/c the infile and Molecule class define the COM, etc
so should not reset those values
cannot do this until now, b/c requires comm->setup() to have setup stencil
invoke pre_neighbor() to insure body xcmimage flags are reset
needed if Verlet::setup::pbc() has remapped/migrated atoms for 2nd run
@ -517,9 +520,14 @@ void FixRigidSmall::init()
void FixRigidSmall::setup_pre_neighbor()
{
if (!staticflag || !infile) setup_bodies_static();
if (!staticflag || (!infile && !onemols)) setup_bodies_static();
else pre_neighbor();
setup_bodies_dynamic();
//if (!staticflag || (!infile)) setup_bodies_static();
//else pre_neighbor();
//setup_bodies_dynamic();
if (!onemols) setup_bodies_dynamic();
}
/* ----------------------------------------------------------------------
@ -880,24 +888,27 @@ void FixRigidSmall::final_integrate_respa(int ilevel, int iloop)
due to first-time definition of rigid body in setup_bodies_static()
or due to box flip
also adjust imagebody = rigid body image flags, due to xcm remap
also adjust rgdimage flags of all atoms in bodies for two effects
then communicate bodies so other procs will know of changes to body xcm
then adjust xcmimage flags of all atoms in bodies via image_shift()
for two effects
(1) change in true image flags due to pbc() call during exchange
(2) change in imagebody due to xcm remap
rgdimage flags are always -1,0,-1 so that body can be unwrapped
xcmimage flags are always -1,0,-1 so that body can be unwrapped
around in-box xcm and stay close to simulation box
if don't do this, then a body could end up very far away
when unwrapped by true image flags, and set_xv() will
compute huge displacements every step to reset coords of
if just inferred unwrapped from atom image flags,
then a body could end up very far away
when unwrapped by true image flags
then set_xv() will compute huge displacements every step to reset coords of
all the body atoms to be back inside the box, ditto for triclinic box flip
note: so just want to avoid that numeric probem?
------------------------------------------------------------------------- */
void FixRigidSmall::pre_neighbor()
{
// acquire ghost bodies via forward comm
// also gets new remapflags needed for adjusting atom image flags
// reset atom2body for owned atoms
// forward comm sets atom2body for ghost atoms
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
domain->remap(b->xcm,b->image);
}
nghost_body = 0;
commflag = FULL_BODY;
@ -905,10 +916,6 @@ void FixRigidSmall::pre_neighbor()
reset_atom2body();
//check(4);
for (int ibody = 0; ibody < nlocal_body; ibody++) {
Body *b = &body[ibody];
domain->remap(b->xcm,b->image);
}
image_shift();
}
@ -1654,7 +1661,8 @@ void FixRigidSmall::ring_farthest(int n, char *cbuf)
/* ----------------------------------------------------------------------
initialization of rigid body attributes
called at setup, so body/atom properties can be changed between runs
unless reading from infile, in which case only called once
unless reading from infile or defined by Molecule class,
in which case only called once
sets extended flags, masstotal, center-of-mass
sets Cartesian and diagonalized inertia tensor
sets body image flags, but only on first call
@ -2180,7 +2188,7 @@ void FixRigidSmall::setup_bodies_static()
one-time initialization of dynamic rigid body attributes
vcm and angmom, computed explicitly from constituent particles
even if wrong for overlapping particles, is OK,
since is just setting initial time=0 Vcm and angmom of the body
since is just setting initial time=0 vcm and angmom of the body
which can be estimated value
------------------------------------------------------------------------- */

View File

@ -255,7 +255,7 @@ NPT/NPH fix must be defined in input script after all rigid fixes,
else the rigid fix contribution to the pressure virial is
incorrect.
W: Cannot count rigid body degrees-of-freedom before bodies are fully initialized
W: Cannot count rigid body degrees-of-freedom before bodies are fully initialized h
This means the temperature associated with the rigid bodies may be
incorrect on this timestep.

View File

@ -123,6 +123,8 @@ void ComputeTempDeformEff::dof_compute()
dof -= domain->dimension * nelectrons;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -174,8 +176,6 @@ double ComputeTempDeformEff::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -87,6 +87,8 @@ void ComputeTempEff::dof_compute()
dof -= domain->dimension * nelectrons;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -120,8 +122,6 @@ double ComputeTempEff::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -64,6 +64,8 @@ void ComputeTemp::dof_compute()
double natoms = group->count(igroup);
dof = domain->dimension * natoms;
dof -= extra_dof + fix_dof;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -96,8 +98,6 @@ double ComputeTemp::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -288,6 +288,8 @@ double ComputeTempChunk::compute_scalar()
MPI_Allreduce(&rcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world);
double dof = nchunk*cdof + adof*allcount;
if (dof < 0.0 && allcount > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
double tfactor = 0.0;
if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
scalar *= tfactor;

View File

@ -71,9 +71,10 @@ void ComputeTempCOM::dof_compute()
{
adjust_dof_fix();
double natoms = group->count(igroup);
int nper = domain->dimension;
dof = nper * natoms;
dof = domain->dimension * natoms;
dof -= extra_dof + fix_dof;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -112,8 +113,6 @@ double ComputeTempCOM::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -99,6 +99,8 @@ void ComputeTempDeform::dof_compute()
double natoms = group->count(igroup);
dof = domain->dimension * natoms;
dof -= extra_dof + fix_dof;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -148,8 +150,6 @@ double ComputeTempDeform::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -196,9 +196,10 @@ void ComputeTempProfile::dof_compute()
{
adjust_dof_fix();
double natoms = group->count(igroup);
int nper = domain->dimension;
dof = nper * natoms;
dof = domain->dimension * natoms;
dof -= extra_dof + fix_dof;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -242,8 +243,6 @@ double ComputeTempProfile::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -127,9 +127,10 @@ void ComputeTempRamp::dof_compute()
{
adjust_dof_fix();
double natoms = group->count(igroup);
int nper = domain->dimension;
dof = nper * natoms;
dof = domain->dimension * natoms;
dof -= extra_dof + fix_dof;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -171,8 +172,6 @@ double ComputeTempRamp::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}

View File

@ -125,6 +125,7 @@ void ComputeTempSphere::dof_compute()
int count,count_all;
adjust_dof_fix();
double natoms = group->count(igroup);
// 6 or 3 dof for extended/point particles for 3d
// 3 or 2 dof for extended/point particles for 2d
@ -165,10 +166,7 @@ void ComputeTempSphere::dof_compute()
// additional adjustments to dof
if (tempbias == 1) {
if (mode == ALL) {
double natoms = group->count(igroup);
dof -= tbias->dof_remove(-1) * natoms;
}
if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms;
} else if (tempbias == 2) {
int *mask = atom->mask;
@ -208,6 +206,8 @@ void ComputeTempSphere::dof_compute()
}
dof -= extra_dof + fix_dof;
if (dof < 0.0 && natoms > 0.0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
else tfactor = 0.0;
}
@ -252,8 +252,6 @@ double ComputeTempSphere::compute_scalar()
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic || tempbias == 2) dof_compute();
if (tfactor == 0.0 && atom->natoms != 0)
error->all(FLERR,"Temperature compute degrees of freedom < 0");
scalar *= tfactor;
return scalar;
}