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

This commit is contained in:
sjplimp 2007-08-03 19:36:35 +00:00
parent a2d5aa2e1f
commit 72ec2d9204
5 changed files with 67 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/fix.h
View File

@ -93,7 +93,7 @@ class Fix : protected Pointers {
virtual double thermo(int) {return 0.0;}
virtual int dof(int) {return 0;}
virtual void dilate(int, double, double, double, double) {}
virtual void deform(int) {}
virtual int modify_param(int, char **) {return 0;}

4
src/fix_nph.h
View File

@ -40,7 +40,7 @@ class FixNPH : public Fix {
double boltz,nktv2p;
double vol0,nkt;
int press_couple,dilate_partial;
int press_couple,allremap;
int p_flag[3]; // 1 if control P on this dim, 0 if not
double p_start[3],p_stop[3];
double p_freq[3],p_target[3];
@ -60,7 +60,7 @@ class FixNPH : public Fix {
int tflag,pflag;
void couple();
void box_dilate(int);
void remap(int);
};
}

132
src/fix_npt.cpp
View File

@ -22,7 +22,6 @@
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "output.h"
#include "modify.h"
#include "compute.h"
#include "kspace.h"
@ -117,7 +116,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
// process extra keywords
drag = 0.0;
dilate_partial = 0;
allremap = 1;
int iarg;
if (press_couple == 0) iarg = 10;
@ -130,8 +129,8 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"dilate") == 0) {
if (iarg+2 > narg) error->all("Illegal fix npt command");
if (strcmp(arg[iarg+1],"all") == 0) dilate_partial = 0;
else if (strcmp(arg[iarg+1],"partial") == 0) dilate_partial = 1;
if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0;
else error->all("Illegal fix npt command");
iarg += 2;
} else error->all("Illegal fix npt command");
@ -284,7 +283,7 @@ void FixNPT::init()
step_respa = ((Respa *) update->integrate)->step;
}
// detect if any rigid fixes exist so rigid bodies move when box is dilated
// detect if any rigid fixes exist so rigid bodies move when box is remapped
// rfix[] = indices to each fix rigid
delete [] rfix;
@ -381,9 +380,9 @@ void FixNPT::initial_integrate()
}
}
// rescale simulation box and all owned atoms by 1/2 step
// remap simulation box and all owned atoms by 1/2 step
box_dilate(0);
remap(0);
// x update by full step only for atoms in NPT group
@ -395,10 +394,10 @@ void FixNPT::initial_integrate()
}
}
// rescale simulation box and all owned atoms by 1/2 step
// remap simulation box and all owned atoms by 1/2 step
// redo KSpace coeffs since volume has changed
box_dilate(0);
remap(0);
if (kspace_flag) force->kspace->setup();
}
@ -468,11 +467,11 @@ void FixNPT::final_integrate()
void FixNPT::initial_integrate_respa(int ilevel, int flag)
{
// if flag = 1, then is 2nd call at outermost level from rRESPA
// perform 2nd half of box dilate on own + ghost atoms and return
// perform 2nd half of box remap on own + ghost atoms and return
// redo KSpace coeffs since volume has changed
if (flag == 1) {
box_dilate(1);
remap(1);
if (kspace_flag) force->kspace->setup();
return;
}
@ -494,7 +493,7 @@ void FixNPT::initial_integrate_respa(int ilevel, int flag)
int *mask = atom->mask;
int nlocal = atom->nlocal;
// outermost level - update eta_dot and omega_dot, apply to v, dilate box
// outermost level - update eta_dot and omega_dot, apply to v, remap box
// all other levels - NVE update of v
// x,v updates only performed for atoms in NPT group
@ -540,9 +539,9 @@ void FixNPT::initial_integrate_respa(int ilevel, int flag)
}
}
// rescale simulation box and all owned atoms by 1/2 step
// remap simulation box and all owned atoms by 1/2 step
box_dilate(0);
remap(0);
} else {
@ -638,89 +637,62 @@ void FixNPT::couple()
}
/* ----------------------------------------------------------------------
dilate the box around center of box
change box size
remap owned or owned+ghost atoms depending on flag
remap all atoms or fix group atoms depending on allremap flag
if rigid bodies exist, scale rigid body centers-of-mass
------------------------------------------------------------------------- */
void FixNPT::box_dilate(int flag)
void FixNPT::remap(int flag)
{
int i,n;
// ctr = geometric center of box in a dimension
// scale owned or owned+ghost atoms depending on flag
// re-define simulation box via xprd/yprd/zprd
// scale atom coords for all atoms or only for fix group atoms
// if fix rigid exists, scale rigid body centers-of-mass
// don't do anything if non-periodic or press style is constant volume
double oldlo,oldhi,ctr;
double **x = atom->x;
int *mask = atom->mask;
if (flag) n = atom->nlocal + atom->nghost;
else n = atom->nlocal;
double oldlo,oldhi,ctr;
// convert pertinent atoms and rigid bodies to lamda coords
if (domain->xperiodic && p_flag[0]) {
oldlo = domain->boxlo[0];
oldhi = domain->boxhi[0];
ctr = 0.5 * (oldlo + oldhi);
domain->boxlo[0] = (oldlo-ctr)*dilation[0] + ctr;
domain->boxhi[0] = (oldhi-ctr)*dilation[0] + ctr;
domain->prd[0] = domain->xprd = domain->boxhi[0] - domain->boxlo[0];
if (dilate_partial) {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
x[i][0] = ctr + (x[i][0]-ctr)*dilation[0];
} else {
for (i = 0; i < n; i++)
x[i][0] = ctr + (x[i][0]-ctr)*dilation[0];
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->
dilate(0,oldlo,oldhi,domain->boxlo[0],domain->boxhi[0]);
if (allremap) domain->x2lamda(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->x2lamda(x[i],x[i]);
}
if (domain->yperiodic && p_flag[1]) {
oldlo = domain->boxlo[1];
oldhi = domain->boxhi[1];
ctr = 0.5 * (oldlo + oldhi);
domain->boxlo[1] = (oldlo-ctr)*dilation[1] + ctr;
domain->boxhi[1] = (oldhi-ctr)*dilation[1] + ctr;
domain->prd[1] = domain->yprd = domain->boxhi[1] - domain->boxlo[1];
if (dilate_partial) {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
x[i][1] = ctr + (x[i][1]-ctr)*dilation[1];
} else {
for (i = 0; i < n; i++)
x[i][1] = ctr + (x[i][1]-ctr)*dilation[1];
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(0);
// reset global and local box to new size/shape
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
oldlo = domain->boxlo[i];
oldhi = domain->boxhi[i];
ctr = 0.5 * (oldlo + oldhi);
domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr;
domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr;
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->
dilate(1,oldlo,oldhi,domain->boxlo[1],domain->boxhi[1]);
}
if (domain->zperiodic && p_flag[2]) {
oldlo = domain->boxlo[2];
oldhi = domain->boxhi[2];
ctr = 0.5 * (oldlo + oldhi);
domain->boxlo[2] = (oldlo-ctr)*dilation[2] + ctr;
domain->boxhi[2] = (oldhi-ctr)*dilation[2] + ctr;
domain->prd[2] = domain->zprd = domain->boxhi[2] - domain->boxlo[2];
if (dilate_partial) {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
x[i][2] = ctr + (x[i][2]-ctr)*dilation[2];
} else {
for (i = 0; i < n; i++)
x[i][2] = ctr + (x[i][2]-ctr)*dilation[2];
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->
dilate(2,oldlo,oldhi,domain->boxlo[2],domain->boxhi[2]);
domain->set_global_box();
domain->set_local_box();
// convert pertinent atoms and rigid bodies back to box coords
if (allremap) domain->lamda2x(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->lamda2x(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(1);
}
/* ----------------------------------------------------------------------

4
src/fix_npt.h
View File

@ -45,7 +45,7 @@ class FixNPT : public Fix {
double t_freq;
double f_eta,eta_dot,eta;
int press_couple,dilate_partial;
int press_couple,allremap;
int p_flag[3]; // 1 if control P on this dim, 0 if not
double p_start[3],p_stop[3];
double p_freq[3],p_target[3];
@ -65,7 +65,7 @@ class FixNPT : public Fix {
int tflag,pflag;
void couple();
void box_dilate(int);
void remap(int);
};
}

19
src/fix_rigid.cpp
View File

@ -847,18 +847,19 @@ int FixRigid::dof(int igroup)
}
/* ----------------------------------------------------------------------
adjust xcm of each rigid body due to box dilation in idim
called by various fixes that change box
adjust xcm of each rigid body due to box deformation
called by various fixes that change box size/shape
flag = 0/1 means map from box to lamda coords or vice versa
------------------------------------------------------------------------- */
void FixRigid::dilate(int idim,
double oldlo, double oldhi, double newlo, double newhi)
void FixRigid::deform(int flag)
{
double ratio;
for (int ibody = 0; ibody < nbody; ibody++) {
ratio = (xcm[ibody][idim] - oldlo) / (oldhi - oldlo);
xcm[ibody][idim] = newlo + ratio*(newhi - newlo);
}
if (flag == 0)
for (int ibody = 0; ibody < nbody; ibody++)
domain->x2lamda(xcm[ibody],xcm[ibody]);
else
for (int ibody = 0; ibody < nbody; ibody++)
domain->lamda2x(xcm[ibody],xcm[ibody]);
}
/* ----------------------------------------------------------------------