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

Browse Source
This commit is contained in:
sjplimp
2009-12-03 17:45:22 +00:00
parent 2ac56e521f
commit 69d350af34
25 changed files with 1065 additions and 270 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/GRANULAR/fix_wall_gran.cpp
View File

@ -46,8 +46,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 10) error->all("Illegal fix wall/gran command");
time_depend = 1;
if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag)
error->all("Fix wall/gran requires atom attributes radius, omega, torque");
@ -166,7 +164,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
if (wiggle) {
double PI = 4.0 * atan(1.0);
omega = 2.0*PI / period;
time_origin = update->ntimestep;
}
// perform initial allocation of atom-based arrays
@ -254,7 +251,7 @@ void FixWallGran::post_force(int vflag)
double whi = hi;
vwall[0] = vwall[1] = vwall[2] = 0.0;
if (wiggle) {
double arg = omega * (update->ntimestep - time_origin) * dt;
double arg = omega * (update->ntimestep - update->beginstep) * dt;
wlo = lo + amplitude - amplitude*cos(arg);
whi = hi + amplitude - amplitude*cos(arg);
vwall[axis] = amplitude*omega*sin(arg);

2
src/PERI/compute_damage_atom.cpp
View File

@ -81,7 +81,7 @@ void ComputeDamageAtom::compute_peratom()
memory->sfree(damage);
nmax = atom->nmax;
damage = (double *) memory->smalloc(nmax*sizeof(double),
"compute/damage/atom:damage");
"damage/atom:damage");
scalar_atom = damage;
}

6
src/compute_centro_atom.cpp
View File

@ -100,7 +100,7 @@ void ComputeCentroAtom::compute_peratom()
memory->sfree(centro);
nmax = atom->nmax;
centro = (double *)
memory->smalloc(nmax*sizeof(double),"compute/centro:centro");
memory->smalloc(nmax*sizeof(double),"centro/atom:centro");
scalar_atom = centro;
}
@ -137,9 +137,9 @@ void ComputeCentroAtom::compute_peratom()
memory->sfree(nearest);
maxneigh = jnum;
distsq = (double *) memory->smalloc(maxneigh*sizeof(double),
"compute/centro:distsq");
"centro/atom:distsq");
nearest = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/centro:nearest");
"centro/atom:nearest");
}
// loop over list of all neighbors within force cutoff

2
src/compute_coord_atom.cpp
View File

@ -101,7 +101,7 @@ void ComputeCoordAtom::compute_peratom()
memory->sfree(coordination);
nmax = atom->nmax;
coordination = (double *)
memory->smalloc(nmax*sizeof(double),"compute/coord/atom:coordination");
memory->smalloc(nmax*sizeof(double),"coord/atom:coordination");
scalar_atom = coordination;
}

2
src/compute_displace_atom.cpp
View File

@ -82,7 +82,7 @@ void ComputeDisplaceAtom::compute_peratom()
memory->destroy_2d_double_array(displace);
nmax = atom->nmax;
displace =
memory->create_2d_double_array(nmax,4,"compute/displace/atom:displace");
memory->create_2d_double_array(nmax,4,"displace/atom:displace");
vector_atom = displace;
}

2
src/compute_ke_atom.cpp
View File

@ -66,7 +66,7 @@ void ComputeKEAtom::compute_peratom()
if (atom->nlocal > nmax) {
memory->sfree(ke);
nmax = atom->nmax;
ke = (double *) memory->smalloc(nmax*sizeof(double),"compute/ke/atom:ke");
ke = (double *) memory->smalloc(nmax*sizeof(double),"ke/atom:ke");
scalar_atom = ke;
}

2
src/compute_pe_atom.cpp
View File

@ -86,7 +86,7 @@ void ComputePEAtom::compute_peratom()
memory->sfree(energy);
nmax = atom->nmax;
energy = (double *)
memory->smalloc(nmax*sizeof(double),"compute/pe/atom:energy");
memory->smalloc(nmax*sizeof(double),"pe/atom:energy");
scalar_atom = energy;
}

2
src/compute_reduce.cpp
View File

@ -336,7 +336,7 @@ double ComputeReduce::compute_one(int m)
maxatom = atom->nmax;
memory->sfree(varatom);
varatom = (double *)
memory->smalloc(maxatom*sizeof(double),"compute/reduce:varatom");
memory->smalloc(maxatom*sizeof(double),"reduce:varatom");
}
input->variable->compute_atom(n,igroup,varatom,1,0);

2
src/compute_reduce_region.cpp
View File

@ -125,7 +125,7 @@ double ComputeReduceRegion::compute_one(int m)
maxatom = atom->nmax;
memory->sfree(varatom);
varatom = (double *)
memory->smalloc(maxatom*sizeof(double),"compute/reduce:varatom");
memory->smalloc(maxatom*sizeof(double),"reduce/region:varatom");
}
input->variable->compute_atom(n,igroup,varatom,1,0);

2
src/compute_stress_atom.cpp
View File

@ -100,7 +100,7 @@ void ComputeStressAtom::compute_peratom()
memory->destroy_2d_double_array(stress);
nmax = atom->nmax;
stress =
memory->create_2d_double_array(nmax,6,"compute/stress/atom:stress");
memory->create_2d_double_array(nmax,6,"stress/atom:stress");
vector_atom = stress;
}

2
src/compute_temp_deform.cpp
View File

@ -227,7 +227,7 @@ void ComputeTempDeform::remove_bias_all()
memory->destroy_2d_double_array(vbiasall);
maxbias = atom->nmax;
vbiasall = memory->create_2d_double_array(maxbias,3,
"compute/temp:vbiasall");
"temp/deform:vbiasall");
}
double lamda[3];

2
src/compute_temp_partial.cpp
View File

@ -190,7 +190,7 @@ void ComputeTempPartial::remove_bias_all()
memory->destroy_2d_double_array(vbiasall);
maxbias = atom->nmax;
vbiasall = memory->create_2d_double_array(maxbias,3,
"compute/temp:vbiasall");
"temp/partial:vbiasall");
}
if (!xflag) {

2
src/compute_temp_profile.cpp
View File

@ -286,7 +286,7 @@ void ComputeTempProfile::remove_bias_all()
memory->destroy_2d_double_array(vbiasall);
maxbias = atom->nmax;
vbiasall = memory->create_2d_double_array(maxbias,3,
"compute/temp:vbiasall");
"temp/profile:vbiasall");
}
int ibin;

2
src/compute_temp_ramp.cpp
View File

@ -254,7 +254,7 @@ void ComputeTempRamp::remove_bias_all()
memory->destroy_2d_double_array(vbiasall);
maxbias = atom->nmax;
vbiasall = memory->create_2d_double_array(maxbias,3,
"compute/temp:vbiasall");
"temp/ramp:vbiasall");
}
double fraction;

2
src/compute_temp_region.cpp
View File

@ -182,7 +182,7 @@ void ComputeTempRegion::remove_bias_all()
memory->destroy_2d_double_array(vbiasall);
maxbias = atom->nmax;
vbiasall = memory->create_2d_double_array(maxbias,3,
"compute/temp:vbiasall");
"temp/region:vbiasall");
}
Region *region = domain->regions[iregion];

86
src/domain.cpp
View File

@ -78,6 +78,7 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
h_ratelo[0] = h_ratelo[1] = h_ratelo[2] = 0.0;
prd_lamda[0] = prd_lamda[1] = prd_lamda[2] = 1.0;
prd_half_lamda[0] = prd_half_lamda[1] = prd_half_lamda[2] = 0.5;
boxlo_lamda[0] = boxlo_lamda[1] = boxlo_lamda[2] = 0.0;
boxhi_lamda[0] = boxhi_lamda[1] = boxhi_lamda[2] = 1.0;
@ -186,9 +187,9 @@ void Domain::set_global_box()
h_inv[1] = 1.0/h[1];
h_inv[2] = 1.0/h[2];
xprd_half = 0.5*xprd;
yprd_half = 0.5*yprd;
zprd_half = 0.5*zprd;
prd_half[0] = xprd_half = 0.5*xprd;
prd_half[1] = yprd_half = 0.5*yprd;
prd_half[2] = zprd_half = 0.5*zprd;
if (triclinic) {
h[3] = yz;
@ -680,7 +681,7 @@ void Domain::remap(double *x, int &image)
remap the point into the periodic box no matter how far away
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
for triclinic, point is converted to lamda coords (0-1) before doing remap
for triclinic, point is converted to lamda coords (0-1) before remap
------------------------------------------------------------------------- */
void Domain::remap(double *x)
@ -722,6 +723,83 @@ void Domain::remap(double *x)
if (triclinic) lamda2x(coord,x);
}
/* ----------------------------------------------------------------------
remap xnew to be within half box length of xold
do it directly, not iteratively, in case is far away
for triclinic, both points are converted to lamda coords (0-1) before remap
------------------------------------------------------------------------- */
void Domain::remap_near(double *xnew, double *xold)
{
int n;
double *coordnew,*coordold,*period,*half;
double lamdanew[3],lamdaold[3];
if (triclinic == 0) {
period = prd;
half = prd_half;
coordnew = xnew;
coordold = xold;
} else {
period = prd_lamda;
half = prd_half_lamda;
x2lamda(xnew,lamdanew);
coordnew = lamdanew;
x2lamda(xold,lamdaold);
coordold = lamdaold;
}
// iterative form
// if (xperiodic) {
// while (xnew[0]-xold[0] > half[0]) xnew[0] -= period[0];
// while (xold[0]-xnew[0] > half[0]) xnew[0] += period[0];
// }
if (xperiodic) {
if (coordnew[0]-coordold[0] > period[0]) {
n = static_cast<int> ((coordnew[0]-coordold[0])/period[0]);
xnew[0] -= n*period[0];
}
while (xnew[0]-xold[0] > half[0]) xnew[0] -= period[0];
if (xold[0]-xnew[0] > period[0]) {
n = static_cast<int> ((xold[0]-xnew[0])/period[0]);
xnew[0] += n*period[0];
}
while (xold[0]-xnew[0] > half[0]) xnew[0] += period[0];
}
if (yperiodic) {
if (coordnew[1]-coordold[1] > period[1]) {
n = static_cast<int> ((coordnew[1]-coordold[1])/period[1]);
xnew[1] -= n*period[1];
}
while (xnew[1]-xold[1] > half[1]) xnew[1] -= period[1];
if (xold[1]-xnew[1] > period[1]) {
n = static_cast<int> ((xold[1]-xnew[1])/period[1]);
xnew[1] += n*period[1];
}
while (xold[1]-xnew[1] > half[1]) xnew[1] += period[1];
}
if (zperiodic) {
if (coordnew[2]-coordold[2] > period[2]) {
n = static_cast<int> ((coordnew[2]-coordold[2])/period[2]);
xnew[2] -= n*period[2];
}
while (xnew[2]-xold[2] > half[2]) xnew[2] -= period[2];
if (xold[2]-xnew[2] > period[2]) {
n = static_cast<int> ((xold[2]-xnew[2])/period[2]);
xnew[2] += n*period[2];
}
while (xold[2]-xnew[2] > half[2]) xnew[2] += period[2];
}
if (triclinic) {
lamda2x(coordnew,xnew);
lamda2x(coordold,xold);
}
}
/* ----------------------------------------------------------------------
unmap the point via image flags
x overwritten with result, don't reset image flag

8
src/domain.h
View File

@ -40,15 +40,18 @@ class Domain : protected Pointers {
double xprd,yprd,zprd; // global box dimensions
double xprd_half,yprd_half,zprd_half; // half dimensions
double prd[3]; // array form of dimensions
double prd_half[3]; // array form of half dimensions
// triclinic box
// xprd,half,prd = same (as if untilt)
// xprd,xprd_half,prd,prd_half =
// same as if untilted
double prd_lamda[3]; // lamda box = (1,1,1)
double prd_half_lamda[3]; // lamda half box = (0.5,0.5,0.5)
double boxlo[3],boxhi[3]; // orthogonal box global bounds
// triclinic box
// boxlo/hi = same (as if untilted)
// boxlo/hi = same as if untilted
double boxlo_lamda[3],boxhi_lamda[3]; // lamda box = (0,1)
double boxlo_bound[3],boxhi_bound[3]; // bounding box of tilted domain
@ -91,6 +94,7 @@ class Domain : protected Pointers {
void pbc();
void remap(double *, int &);
void remap(double *);
void remap_near(double *, double *);
void unmap(double *, int);
void unmap(double *, int, double *);
void minimum_image(double &, double &, double &);

10
src/fix_gravity.cpp
View File

@ -33,8 +33,6 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 5) error->all("Illegal fix gravity command");
time_depend = 1;
magnitude = atof(arg[3]);
if (strcmp(arg[4],"chute") == 0) {
@ -88,8 +86,6 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
zgrav = 0.0;
}
}
time_initial = update->ntimestep;
}
/* ---------------------------------------------------------------------- */
@ -138,15 +134,15 @@ void FixGravity::post_force(int vflag)
if (style == GRADIENT) {
if (domain->dimension == 3) {
double phi_current = degree2rad *
(phi + (update->ntimestep-time_initial)*dt*phigrad*360.0);
(phi + (update->ntimestep-update->beginstep)*dt*phigrad*360.0);
double theta_current = degree2rad *
(theta + (update->ntimestep-time_initial)*dt*thetagrad*360.0);
(theta + (update->ntimestep-update->beginstep)*dt*thetagrad*360.0);
xgrav = sin(theta_current) * cos(phi_current);
ygrav = sin(theta_current) * sin(phi_current);
zgrav = cos(theta_current);
} else {
double theta_current = degree2rad *
(theta + (update->ntimestep-time_initial)*dt*thetagrad*360.0);
(theta + (update->ntimestep-update->beginstep)*dt*thetagrad*360.0);
xgrav = sin(theta_current);
ygrav = cos(theta_current);
}

2
src/fix_gravity.h
View File

@ -30,7 +30,7 @@ class FixGravity : public Fix {
void post_force_respa(int, int, int);
private:
int style,time_initial;
int style;
double magnitude,dt;
double phi,theta,phigrad,thetagrad;
double xdir,ydir,zdir;

890
src/fix_move.cpp Normal file
View File

@ -0,0 +1,890 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "fix_move.h"
#include "atom.h"
#include "group.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "lattice.h"
#include "comm.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{LINEAR,WIGGLE,ROTATE,VARIABLE};
enum{EQUAL,ATOM};
/* ---------------------------------------------------------------------- */
FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 4) error->all("Illegal fix move command");
restart_peratom = 1;
peratom_flag = 1;
size_peratom = 3;
peratom_freq = 1;
time_integrate = 1;
// parse args
int iarg;
xvarstr = yvarstr = zvarstr = NULL;
vxvarstr = vyvarstr = vzvarstr = NULL;
if (strcmp(arg[3],"linear") == 0) {
if (narg < 7) error->all("Illegal fix move command");
iarg = 7;
mstyle = LINEAR;
if (strcmp(arg[4],"NULL") == 0) vxflag = 0;
else {
vxflag = 1;
vx = atof(arg[4]);
}
if (strcmp(arg[5],"NULL") == 0) vyflag = 0;
else {
vyflag = 1;
vy = atof(arg[5]);
}
if (strcmp(arg[6],"NULL") == 0) vzflag = 0;
else {
vzflag = 1;
vz = atof(arg[6]);
}
} else if (strcmp(arg[3],"wiggle") == 0) {
if (narg < 8) error->all("Illegal fix move command");
iarg = 8;
mstyle = WIGGLE;
if (strcmp(arg[4],"NULL") == 0) axflag = 0;
else {
axflag = 1;
ax = atof(arg[4]);
}
if (strcmp(arg[5],"NULL") == 0) ayflag = 0;
else {
ayflag = 1;
ay = atof(arg[5]);
}
if (strcmp(arg[6],"NULL") == 0) azflag = 0;
else {
azflag = 1;
az = atof(arg[6]);
}
period = atof(arg[7]);
} else if (strcmp(arg[3],"rotate") == 0) {
if (narg < 11) error->all("Illegal fix move command");
iarg = 11;
mstyle = ROTATE;
point[0] = atof(arg[4]);
point[1] = atof(arg[5]);
point[2] = atof(arg[6]);
axis[0] = atof(arg[7]);
axis[1] = atof(arg[8]);
axis[2] = atof(arg[9]);
period = atof(arg[10]);
} else if (strcmp(arg[3],"variable") == 0) {
if (narg < 10) error->all("Illegal fix move command");
iarg = 10;
mstyle = VARIABLE;
if (strcmp(arg[4],"NULL") == 0) xvarstr = NULL;
else if (strstr(arg[4],"v_") == arg[4]) {
int n = strlen(&arg[4][2]) + 1;
xvarstr = new char[n];
strcpy(xvarstr,&arg[4][2]);
} else error->all("Illegal fix move command");
if (strcmp(arg[5],"NULL") == 0) yvarstr = NULL;
else if (strstr(arg[5],"v_") == arg[5]) {
int n = strlen(&arg[5][2]) + 1;
yvarstr = new char[n];
strcpy(yvarstr,&arg[5][2]);
} else error->all("Illegal fix move command");
if (strcmp(arg[6],"NULL") == 0) zvarstr = NULL;
else if (strstr(arg[6],"v_") == arg[6]) {
int n = strlen(&arg[6][2]) + 1;
zvarstr = new char[n];
strcpy(zvarstr,&arg[6][2]);
} else error->all("Illegal fix move command");
if (strcmp(arg[7],"NULL") == 0) vxvarstr = NULL;
else if (strstr(arg[7],"v_") == arg[7]) {
int n = strlen(&arg[7][2]) + 1;
vxvarstr = new char[n];
strcpy(vxvarstr,&arg[7][2]);
} else error->all("Illegal fix move command");
if (strcmp(arg[8],"NULL") == 0) vyvarstr = NULL;
else if (strstr(arg[8],"v_") == arg[8]) {
int n = strlen(&arg[8][2]) + 1;
vyvarstr = new char[n];
strcpy(vyvarstr,&arg[8][2]);
} else error->all("Illegal fix move command");
if (strcmp(arg[9],"NULL") == 0) vzvarstr = NULL;
else if (strstr(arg[9],"v_") == arg[9]) {
int n = strlen(&arg[9][2]) + 1;
vzvarstr = new char[n];
strcpy(vzvarstr,&arg[9][2]);
} else error->all("Illegal fix move command");
} else error->all("Illegal fix move command");
// optional args
int scaleflag = 1;
while (iarg < narg) {
if (strcmp(arg[iarg],"units") == 0) {
if (iarg+2 > narg) error->all("Illegal fix move command");
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
else error->all("Illegal fix move command");
iarg += 2;
} else error->all("Illegal fix move command");
}
// error checks and warnings
if (domain->dimension == 2) {
if (mstyle == LINEAR && vzflag && vz != 0.0)
error->all("Fix move cannot set linear z motion for 2d problem");
if (mstyle == WIGGLE && azflag && az != 0.0)
error->all("Fix move cannot set wiggle z motion for 2d problem");
if (mstyle == ROTATE && (axis[0] != 0.0 || axis[1] != 0.0))
error->all("Fix move cannot rotate aroung non z-axis for 2d problem");
if (mstyle == VARIABLE && (zvarstr || vzvarstr))
error->all("Fix move cannot define z or vz variable for 2d problem");
}
if (atom->angmom_flag && comm->me == 0)
error->warning("Fix move does not update angular momentum");
if (atom->quat_flag && comm->me == 0)
error->warning("Fix move does not update quaternions");
// setup scaling and apply scaling factors to velocity & amplitude
if ((mstyle == LINEAR || mstyle == WIGGLE || mstyle == ROTATE) &&
scaleflag) {
if (domain->lattice == NULL)
error->all("Use of fix move with undefined lattice");
double xscale,yscale,zscale;
if (scaleflag) {
xscale = domain->lattice->xlattice;
yscale = domain->lattice->ylattice;
zscale = domain->lattice->zlattice;
}
else xscale = yscale = zscale = 1.0;
if (mstyle == LINEAR) {
if (vxflag) vx *= xscale;
if (vyflag) vy *= yscale;
if (vzflag) vz *= zscale;
} else if (mstyle == WIGGLE) {
if (axflag) ax *= xscale;
if (ayflag) ay *= yscale;
if (azflag) az *= zscale;
} else if (mstyle == ROTATE) {
point[0] *= xscale;
point[1] *= yscale;
point[2] *= zscale;
}
}
// set omega_rotate from period
if (mstyle == WIGGLE || mstyle == ROTATE) {
double PI = 4.0 * atan(1.0);
omega_rotate = 2.0*PI / period;
}
// runit = unit vector along rotation axis
if (mstyle == ROTATE) {
double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
if (len == 0.0)
error->all("Fix move cannot have 0 length rotation vector");
runit[0] = axis[0]/len;
runit[1] = axis[1]/len;
runit[2] = axis[2]/len;
}
// set omega_flag if particles store omega
omega_flag = atom->omega_flag;
// perform initial allocation of atom-based array
// register with Atom class
xoriginal = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
atom->add_callback(1);
maxatom = 0;
displace = velocity = NULL;
// xoriginal = initial unwrapped positions of atoms
double **x = atom->x;
int *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
FixMove::~FixMove()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
atom->delete_callback(id,1);
// delete locally stored arrays
memory->destroy_2d_double_array(xoriginal);
memory->destroy_2d_double_array(displace);
memory->destroy_2d_double_array(velocity);
delete [] xvarstr;
delete [] yvarstr;
delete [] zvarstr;
delete [] vxvarstr;
delete [] vyvarstr;
delete [] vzvarstr;
}
/* ---------------------------------------------------------------------- */
int FixMove::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixMove::init()
{
dt = update->dt;
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
// set indices and style of all variables
if (mstyle == VARIABLE) {
if (xvarstr) {
xvar = input->variable->find(xvarstr);
if (xvar < 0) error->all("Variable name for fix move does not exist");
if (input->variable->equalstyle(xvar)) xvarstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xvarstyle = ATOM;
else error->all("Variable for fix move is invalid style");
}
if (yvarstr) {
yvar = input->variable->find(yvarstr);
if (yvar < 0) error->all("Variable name for fix move does not exist");
if (input->variable->equalstyle(yvar)) yvarstyle = EQUAL;
else if (input->variable->atomstyle(yvar)) yvarstyle = ATOM;
else error->all("Variable for fix move is invalid style");
}
if (zvarstr) {
zvar = input->variable->find(zvarstr);
if (zvar < 0) error->all("Variable name for fix move does not exist");
if (input->variable->equalstyle(zvar)) zvarstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zvarstyle = ATOM;
else error->all("Variable for fix move is invalid style");
}
if (vxvarstr) {
vxvar = input->variable->find(vxvarstr);
if (vxvar < 0) error->all("Variable name for fix move does not exist");
if (input->variable->equalstyle(vxvar)) vxvarstyle = EQUAL;
else if (input->variable->atomstyle(vxvar)) vxvarstyle = ATOM;
else error->all("Variable for fix move is invalid style");
}
if (vyvarstr) {
vyvar = input->variable->find(vyvarstr);
if (vyvar < 0) error->all("Variable name for fix move does not exist");
if (input->variable->equalstyle(vyvar)) vyvarstyle = EQUAL;
else if (input->variable->atomstyle(vyvar)) vyvarstyle = ATOM;
else error->all("Variable for fix move is invalid style");
}
if (vzvarstr) {
vzvar = input->variable->find(vzvarstr);
if (vzvar < 0) error->all("Variable name for fix move does not exist");
if (input->variable->equalstyle(vzvar)) vzvarstyle = EQUAL;
else if (input->variable->atomstyle(vzvar)) vzvarstyle = ATOM;
else error->all("Variable for fix move is invalid style");
}
displaceflag = velocityflag = 0;
if (xvarstr && xvarstyle == ATOM) displaceflag = 1;
if (yvarstr && yvarstyle == ATOM) displaceflag = 1;
if (zvarstr && zvarstyle == ATOM) displaceflag = 1;
if (vxvarstr && vxvarstyle == ATOM) velocityflag = 1;
if (vyvarstr && vyvarstyle == ATOM) velocityflag = 1;
if (vzvarstr && vzvarstyle == ATOM) velocityflag = 1;
}
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ----------------------------------------------------------------------
set x,v of particles
------------------------------------------------------------------------- */
void FixMove::initial_integrate(int vflag)
{
double dtfm;
double xold[3],a[3],b[3],c[3],d[3],disp[3];
double x0dotr,dx,dy,dz;
double delta = (update->ntimestep - update->beginstep) * dt;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **omega = atom->omega;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// for linear: X = X0 + V*dt
if (mstyle == LINEAR) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
xold[0] = x[i][0];
xold[1] = x[i][1];
xold[2] = x[i][2];
if (vxflag) {
v[i][0] = vx;
x[i][0] = xoriginal[i][0] + vx*delta;
} else if (rmass) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
x[i][0] += dtv * v[i][0];
} else {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
x[i][0] += dtv * v[i][0];
}
if (vyflag) {
v[i][1] = vy;
x[i][1] = xoriginal[i][1] + vy*delta;
} else if (rmass) {
dtfm = dtf / rmass[i];
v[i][1] += dtfm * f[i][1];
x[i][1] += dtv * v[i][1];
} else {
dtfm = dtf / mass[type[i]];
v[i][1] += dtfm * f[i][1];
x[i][1] += dtv * v[i][1];
}
if (vzflag) {
v[i][2] = vz;
x[i][2] = xoriginal[i][2] + vz*delta;
} else if (rmass) {
dtfm = dtf / rmass[i];
v[i][2] += dtfm * f[i][2];
x[i][2] += dtv * v[i][2];
} else {
dtfm = dtf / mass[type[i]];
v[i][2] += dtfm * f[i][2];
x[i][2] += dtv * v[i][2];
}
domain->remap_near(x[i],xold);
}
}
// for wiggle: X = X0 + A sin(w*dt)
} else if (mstyle == WIGGLE) {
double arg = omega_rotate * delta;
double sine = sin(arg);
double cosine = cos(arg);
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
xold[0] = x[i][0];
xold[1] = x[i][1];
xold[2] = x[i][2];
if (axflag) {
v[i][0] = ax*omega_rotate*cosine;
x[i][0] = xoriginal[i][0] + ax*sine;
} else if (rmass) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
x[i][0] += dtv * v[i][0];
} else {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
x[i][0] += dtv * v[i][0];
}
if (ayflag) {
v[i][1] = ay*omega_rotate*cosine;
x[i][1] = xoriginal[i][1] + ay*sine;
} else if (rmass) {
dtfm = dtf / rmass[i];
v[i][1] += dtfm * f[i][1];
x[i][1] += dtv * v[i][1];
} else {
dtfm = dtf / mass[type[i]];
v[i][1] += dtfm * f[i][1];
x[i][1] += dtv * v[i][1];
}
if (azflag) {
v[i][2] = az*omega_rotate*cosine;
x[i][2] = xoriginal[i][2] + az*sine;
} else if (rmass) {
dtfm = dtf / rmass[i];
v[i][2] += dtfm * f[i][2];
x[i][2] += dtv * v[i][2];
} else {
dtfm = dtf / mass[type[i]];
v[i][2] += dtfm * f[i][2];
x[i][2] += dtv * v[i][2];
}
domain->remap_near(x[i],xold);
}
}
// for rotate:
// P = vector = point of rotation
// R = vector = axis of rotation
// w = omega of rotation (from period)
// X0 = initial coord of atom
// R0 = unit vector for R
// C = (X0 dot R0) R0 = projection of atom coord onto R
// D = X0 - P = vector from P to X0
// A = D - C = vector from R line to X0
// B = R0 cross A = vector perp to A in plane of rotation
// A,B define plane of circular rotation around R line
// X = P + C + A cos(w*dt) + B sin(w*dt)
// V = w R0 cross (A cos(w*dt) + B sin(w*dt))
} else if (mstyle == ROTATE) {
double arg = omega_rotate * delta;
double sine = sin(arg);
double cosine = cos(arg);
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
xold[0] = x[i][0];
xold[1] = x[i][1];
xold[2] = x[i][2];
x0dotr = xoriginal[i][0]*runit[0] +
xoriginal[i][1]*runit[1] + xoriginal[i][2]*runit[2];
c[0] = x0dotr * runit[0];
c[1] = x0dotr * runit[1];
c[2] = x0dotr * runit[2];
d[0] = xoriginal[i][0] - point[0];
d[1] = xoriginal[i][1] - point[1];
d[2] = xoriginal[i][2] - point[2];
a[0] = d[0] - c[0];
a[1] = d[1] - c[1];
a[2] = d[2] - c[2];
b[0] = runit[1]*a[2] - runit[2]*a[1];
b[1] = runit[2]*a[0] - runit[0]*a[2];
b[2] = runit[0]*a[1] - runit[1]*a[0];
disp[0] = a[0]*cosine + b[0]*sine;
disp[1] = a[1]*cosine + b[1]*sine;
disp[2] = a[2]*cosine + b[2]*sine;
x[i][0] = point[0] + c[0] + disp[0];
x[i][1] = point[1] + c[1] + disp[1];
x[i][2] = point[2] + c[2] + disp[2];
v[i][0] = omega_rotate * (runit[1]*disp[2] - runit[2]*disp[1]);
v[i][1] = omega_rotate * (runit[2]*disp[0] - runit[0]*disp[2]);
v[i][2] = omega_rotate * (runit[0]*disp[1] - runit[1]*disp[0]);
if (omega_flag) {
omega[i][0] = omega_rotate*runit[0];
omega[i][1] = omega_rotate*runit[1];
omega[i][2] = omega_rotate*runit[2];
}
domain->remap_near(x[i],xold);
}
}
// for variable: compute x,v from variables
} else if (mstyle == VARIABLE) {
// reallocate displace and velocity arrays as necessary
if ((displaceflag || velocityflag) && nlocal > maxatom) {
maxatom = atom->nmax;
if (displaceflag) {
memory->destroy_2d_double_array(displace);
displace = memory->create_2d_double_array(maxatom,3,"move:displace");
}
if (velocityflag) {
memory->destroy_2d_double_array(velocity);
velocity = memory->create_2d_double_array(maxatom,3,"move:velocity");
}
}
// pre-compute variable values
if (xvarstr) {
if (xvarstyle == EQUAL) dx = input->variable->compute_equal(xvar);
else input->variable->compute_atom(xvar,igroup,&displace[0][0],3,0);
}
if (yvarstr) {
if (yvarstyle == EQUAL) dy = input->variable->compute_equal(yvar);
else input->variable->compute_atom(yvar,igroup,&displace[0][1],3,0);
}
if (zvarstr) {
if (zvarstyle == EQUAL) dz = input->variable->compute_equal(zvar);
else input->variable->compute_atom(zvar,igroup,&displace[0][2],3,0);
}
if (vxvarstr) {
if (vxvarstyle == EQUAL) vx = input->variable->compute_equal(vxvar);
else input->variable->compute_atom(vxvar,igroup,&velocity[0][0],3,0);
}
if (vyvarstr) {
if (vyvarstyle == EQUAL) vy = input->variable->compute_equal(vyvar);
else input->variable->compute_atom(vyvar,igroup,&velocity[0][1],3,0);
}
if (vzvarstr) {
if (vzvarstyle == EQUAL) vz = input->variable->compute_equal(vzvar);
else input->variable->compute_atom(vzvar,igroup,&velocity[0][2],3,0);
}
// update x,v
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
xold[0] = x[i][0];
xold[1] = x[i][1];
xold[2] = x[i][2];
if (xvarstr && vxvarstr) {
if (vxvarstyle == EQUAL) v[i][0] = vx;
else v[i][0] = velocity[i][0];
if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx;
else x[i][0] = xoriginal[i][0] + displace[i][0];
} else if (xvarstr) {
if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx;
else x[i][0] = xoriginal[i][0] + displace[i][0];
} else if (vxvarstr) {
if (vxvarstyle == EQUAL) v[i][0] = vx;
else v[i][0] = velocity[i][0];
if (rmass) {
dtfm = dtf / rmass[i];
x[i][0] += dtv * v[i][0];
} else {
dtfm = dtf / mass[type[i]];
x[i][0] += dtv * v[i][0];
}
} else {
if (rmass) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
x[i][0] += dtv * v[i][0];
} else {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
x[i][0] += dtv * v[i][0];
}
}
if (yvarstr && vyvarstr) {
if (vyvarstyle == EQUAL) v[i][1] = vy;
else v[i][1] = velocity[i][1];
if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy;
else x[i][1] = xoriginal[i][1] + displace[i][1];
} else if (yvarstr) {
if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy;
else x[i][1] = xoriginal[i][1] + displace[i][1];
} else if (vyvarstr) {
if (vyvarstyle == EQUAL) v[i][1] = vy;
else v[i][1] = velocity[i][1];
if (rmass) {
dtfm = dtf / rmass[i];
x[i][1] += dtv * v[i][1];
} else {
dtfm = dtf / mass[type[i]];
x[i][1] += dtv * v[i][1];
}
} else {
if (rmass) {
dtfm = dtf / rmass[i];
v[i][1] += dtfm * f[i][1];
x[i][1] += dtv * v[i][1];
} else {
dtfm = dtf / mass[type[i]];
v[i][1] += dtfm * f[i][1];
x[i][1] += dtv * v[i][1];
}
}
if (zvarstr && vzvarstr) {
if (vzvarstyle == EQUAL) v[i][2] = vz;
else v[i][2] = velocity[i][2];
if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz;
else x[i][2] = xoriginal[i][2] + displace[i][2];
} else if (zvarstr) {
if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz;
else x[i][2] = xoriginal[i][2] + displace[i][2];
} else if (vzvarstr) {
if (vzvarstyle == EQUAL) v[i][2] = vz;
else v[i][2] = velocity[i][2];
if (rmass) {
dtfm = dtf / rmass[i];
x[i][2] += dtv * v[i][2];
} else {
dtfm = dtf / mass[type[i]];
x[i][2] += dtv * v[i][2];
}
} else {
if (rmass) {
dtfm = dtf / rmass[i];
v[i][2] += dtfm * f[i][2];
x[i][2] += dtv * v[i][2];
} else {
dtfm = dtf / mass[type[i]];
v[i][2] += dtfm * f[i][2];
x[i][2] += dtv * v[i][2];
}
}
domain->remap_near(x[i],xold);
}
}
}
}
/* ----------------------------------------------------------------------
final NVE of particles with NULL components
------------------------------------------------------------------------- */
void FixMove::final_integrate()
{
double dtfm;
int xflag = 1;
if (mstyle == LINEAR && vxflag) xflag = 0;
else if (mstyle == WIGGLE && axflag) xflag = 0;
else if (mstyle == ROTATE) xflag = 0;
else if (mstyle == VARIABLE && (xvarstr || vxvarstr)) xflag = 0;
int yflag = 1;
if (mstyle == LINEAR && vyflag) yflag = 0;
else if (mstyle == WIGGLE && ayflag) yflag = 0;
else if (mstyle == ROTATE) yflag = 0;
else if (mstyle == VARIABLE && (yvarstr || vyvarstr)) yflag = 0;
int zflag = 1;
if (mstyle == LINEAR && vzflag) zflag = 0;
else if (mstyle == WIGGLE && azflag) zflag = 0;
else if (mstyle == ROTATE) zflag = 0;
else if (mstyle == VARIABLE && (zvarstr || vzvarstr)) zflag = 0;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **omega = atom->omega;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (xflag)
if (rmass) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
} else {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
}
if (yflag)
if (rmass) {
dtfm = dtf / rmass[i];
v[i][1] += dtfm * f[i][1];
} else {
dtfm = dtf / mass[type[i]];
v[i][1] += dtfm * f[i][1];
}
if (zflag)
if (rmass) {
dtfm = dtf / rmass[i];
v[i][2] += dtfm * f[i][2];
} else {
dtfm = dtf / mass[type[i]];
v[i][2] += dtfm * f[i][2];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixMove::initial_integrate_respa(int vflag, int ilevel, int flag)
{
if (flag) return; // only used by NPT,NPH
// outermost level - update v and x
// all other levels - nothing
if (ilevel == nlevels_respa-1) initial_integrate(vflag);
}
/* ---------------------------------------------------------------------- */
void FixMove::final_integrate_respa(int ilevel)
{
if (ilevel == nlevels_respa-1) final_integrate();
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixMove::memory_usage()
{
double bytes = atom->nmax*3 * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate atom-based array
------------------------------------------------------------------------- */
void FixMove::grow_arrays(int nmax)
{
xoriginal =
memory->grow_2d_double_array(xoriginal,nmax,3,"move:xoriginal");
vector_atom = xoriginal;
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
void FixMove::copy_arrays(int i, int j)
{
xoriginal[j][0] = xoriginal[i][0];
xoriginal[j][1] = xoriginal[i][1];
xoriginal[j][2] = xoriginal[i][2];
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixMove::pack_exchange(int i, double *buf)
{
buf[0] = xoriginal[i][0];
buf[1] = xoriginal[i][1];
buf[2] = xoriginal[i][2];
return 3;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixMove::unpack_exchange(int nlocal, double *buf)
{
xoriginal[nlocal][0] = buf[0];
xoriginal[nlocal][1] = buf[1];
xoriginal[nlocal][2] = buf[2];
return 3;
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */
int FixMove::pack_restart(int i, double *buf)
{
buf[0] = 4;
buf[1] = xoriginal[i][0];
buf[2] = xoriginal[i][1];
buf[3] = xoriginal[i][2];
return 4;
}
/* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */
void FixMove::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
// skip to Nth set of extra values
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;
xoriginal[nlocal][0] = extra[nlocal][m++];
xoriginal[nlocal][1] = extra[nlocal][m++];
xoriginal[nlocal][2] = extra[nlocal][m++];
}
/* ----------------------------------------------------------------------
maxsize of any atom's restart data
------------------------------------------------------------------------- */
int FixMove::maxsize_restart()
{
return 4;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
int FixMove::size_restart(int nlocal)
{
return 4;
}

62
src/fix_move.h Normal file
View File

@ -0,0 +1,62 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef FIX_MOVE_H
#define FIX_MOVE_H
#include "fix.h"
namespace LAMMPS_NS {
class FixMove : public Fix {
public:
FixMove(class LAMMPS *, int, char **);
~FixMove();
int setmask();
void init();
void initial_integrate(int);
void final_integrate();
void initial_integrate_respa(int, int, int);
void final_integrate_respa(int);
double memory_usage();
void grow_arrays(int);
void copy_arrays(int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
int pack_restart(int, double *);
void unpack_restart(int, int);
int maxsize_restart();
int size_restart(int);
private:
char *xvarstr,*yvarstr,*zvarstr,*vxvarstr,*vyvarstr,*vzvarstr;
int mstyle;
int vxflag,vyflag,vzflag,axflag,ayflag,azflag;
double vx,vy,vz,ax,ay,az;
double period,omega_rotate;
double point[3],axis[3],runit[3];
double dt,dtv,dtf;
int xvar,yvar,zvar,vxvar,vyvar,vzvar;
int xvarstyle,yvarstyle,zvarstyle,vxvarstyle,vyvarstyle,vzvarstyle;
int omega_flag,nlevels_respa;
double **xoriginal; // original coords of atoms
int displaceflag,velocityflag;
int maxatom;
double **displace,**velocity;
};
}
#endif

186
src/fix_wiggle.cpp
View File

@ -1,186 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_wiggle.h"
#include "atom.h"
#include "update.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixWiggle::FixWiggle(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg != 6) error->all("Illegal fix wiggle command");
time_depend = 1;
// parse command-line args
if (strcmp(arg[3],"x") == 0) axis = 0;
else if (strcmp(arg[3],"y") == 0) axis = 1;
else if (strcmp(arg[3],"z") == 0) axis = 2;
else error->all("Illegal fix wiggle command");
amplitude = atof(arg[4]);
period = atof(arg[5]);
// perform initial allocation of atom-based array
// register with Atom class
original = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
// oscillation parameters
double PI = 4.0 * atan(1.0);
omega = 2.0*PI / period;
time_origin = update->ntimestep;
// original = initial atom coord in wiggled direction
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) original[i] = x[i][axis];
else original[i] = 0.0;
}
/* ---------------------------------------------------------------------- */
FixWiggle::~FixWiggle()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
// delete locally stored array
memory->sfree(original);
}
/* ---------------------------------------------------------------------- */
int FixWiggle::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixWiggle::init()
{
dt = update->dt;
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixWiggle::post_force(int vflag)
{
double arg = omega * (update->ntimestep - time_origin) * dt;
double cosine = cos(arg);
double sine = sin(arg);
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][axis] = original[i] + amplitude - amplitude*cosine;
v[i][axis] = amplitude*omega*sine;
f[i][axis] = 0.0;
}
}
}
/* ---------------------------------------------------------------------- */
void FixWiggle::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixWiggle::memory_usage()
{
double bytes = atom->nmax * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate local atom-based array
------------------------------------------------------------------------- */
void FixWiggle::grow_arrays(int nmax)
{
original = (double *)
memory->srealloc(original,nmax*sizeof(double),"wiggle:original");
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
void FixWiggle::copy_arrays(int i, int j)
{
original[j] = original[i];
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixWiggle::pack_exchange(int i, double *buf)
{
buf[0] = original[i];
return 1;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixWiggle::unpack_exchange(int nlocal, double *buf)
{
original[nlocal] = buf[0];
return 1;
}
/* ---------------------------------------------------------------------- */
void FixWiggle::reset_dt()
{
dt = update->dt;
}

47
src/fix_wiggle.h
View File

@ -1,47 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef FIX_WIGGLE_H
#define FIX_WIGGLE_H
#include "fix.h"
namespace LAMMPS_NS {
class FixWiggle : public Fix {
public:
FixWiggle(class LAMMPS *, int, char **);
~FixWiggle();
int setmask();
void init();
void post_force(int);
void post_force_respa(int, int, int);
double memory_usage();
void grow_arrays(int);
void copy_arrays(int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
void reset_dt();
private:
double dt;
double *original;
double amplitude,period,omega;
int axis,time_origin;
int nlevels_respa;
};
}
#endif

4
src/style.h
View File

@ -174,6 +174,7 @@ DumpStyle(xyz,DumpXYZ)
#include "fix_minimize.h"
#include "fix_msd.h"
#include "fix_momentum.h"
#include "fix_move.h"
#include "fix_nph.h"
#include "fix_npt.h"
#include "fix_npt_sphere.h"
@ -208,7 +209,6 @@ DumpStyle(xyz,DumpXYZ)
#include "fix_wall_lj126.h"
#include "fix_wall_lj93.h"
#include "fix_wall_reflect.h"
#include "fix_wiggle.h"
#endif
#ifdef FixClass
@ -235,6 +235,7 @@ FixStyle(langevin,FixLangevin)
FixStyle(lineforce,FixLineForce)
FixStyle(MINIMIZE,FixMinimize)
FixStyle(momentum,FixMomentum)
FixStyle(move,FixMove)
FixStyle(msd,FixMSD)
FixStyle(nph,FixNPH)
FixStyle(npt,FixNPT)
@ -270,7 +271,6 @@ FixStyle(viscous,FixViscous)
FixStyle(wall/lj126,FixWallLJ126)
FixStyle(wall/lj93,FixWallLJ93)
FixStyle(wall/reflect,FixWallReflect)
FixStyle(wiggle,FixWiggle)
#endif
#ifdef ImproperInclude

3
src/variable.cpp
View File

@ -398,7 +398,8 @@ double Variable::compute_equal(int ivar)
/* ----------------------------------------------------------------------
compute result of atom-style variable evaluation
stride used since result may not be contiguous memory locs
only computed for atoms in igroup, else result is 0.0
answers are placed every stride locations into result
if sumflag, add variable values to existing result
------------------------------------------------------------------------- */