mirror of https://github.com/lammps/lammps.git
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7145 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
6cebfc12b5
commit
83b55b844c
|
|
@ -8,6 +8,8 @@ if (test $1 = 1) then
|
|||
cp fix_nph_asphere.cpp ..
|
||||
cp fix_npt_asphere.cpp ..
|
||||
cp fix_nve_asphere.cpp ..
|
||||
cp fix_nve_line.cpp ..
|
||||
cp fix_nve_tri.cpp ..
|
||||
cp fix_nvt_asphere.cpp ..
|
||||
cp pair_gayberne.cpp ..
|
||||
cp pair_resquared.cpp ..
|
||||
|
|
@ -18,6 +20,8 @@ if (test $1 = 1) then
|
|||
cp fix_nph_asphere.h ..
|
||||
cp fix_npt_asphere.h ..
|
||||
cp fix_nve_asphere.h ..
|
||||
cp fix_nve_line.h ..
|
||||
cp fix_nve_tri.h ..
|
||||
cp fix_nvt_asphere.h ..
|
||||
cp pair_gayberne.h ..
|
||||
cp pair_resquared.h ..
|
||||
|
|
@ -30,6 +34,8 @@ elif (test $1 = 0) then
|
|||
rm -f ../fix_nph_asphere.cpp
|
||||
rm -f ../fix_npt_asphere.cpp
|
||||
rm -f ../fix_nve_asphere.cpp
|
||||
rm -f ../fix_nve_line.cpp
|
||||
rm -f ../fix_nve_tri.cpp
|
||||
rm -f ../fix_nvt_asphere.cpp
|
||||
rm -f ../pair_gayberne.cpp
|
||||
rm -f ../pair_resquared.cpp
|
||||
|
|
@ -40,6 +46,8 @@ elif (test $1 = 0) then
|
|||
rm -f ../fix_nph_asphere.h
|
||||
rm -f ../fix_npt_asphere.h
|
||||
rm -f ../fix_nve_asphere.h
|
||||
rm -f ../fix_nve_line.h
|
||||
rm -f ../fix_nve_tri.h
|
||||
rm -f ../fix_nvt_asphere.h
|
||||
rm -f ../pair_gayberne.h
|
||||
rm -f ../pair_resquared.h
|
||||
|
|
|
|||
|
|
@ -16,6 +16,8 @@
|
|||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec_ellipsoid.h"
|
||||
#include "atom_vec_line.h"
|
||||
#include "atom_vec_tri.h"
|
||||
#include "update.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
|
|
@ -36,9 +38,12 @@ ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
// error check
|
||||
|
||||
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
|
||||
if (!avec)
|
||||
error->all(FLERR,"Compute erotate/asphere requires atom style ellipsoid");
|
||||
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
|
||||
avec_line = (AtomVecLine *) atom->style_match("line");
|
||||
avec_tri = (AtomVecTri *) atom->style_match("tri");
|
||||
if (!avec_ellipsoid && !avec_line && !avec_tri)
|
||||
error->all(FLERR,"Compute erotate/asphere requires "
|
||||
"atom style ellipsoid or line or tri");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -49,13 +54,21 @@ void ComputeERotateAsphere::init()
|
|||
// no point particles allowed, spherical is OK
|
||||
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *line = atom->line;
|
||||
int *tri = atom->tri;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
int flag;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
if (ellipsoid[i] < 0)
|
||||
if (mask[i] & groupbit) {
|
||||
flag = 0;
|
||||
if (ellipsoid && ellipsoid[i] >= 0) flag = 1;
|
||||
if (line && line[i] >= 0) flag = 1;
|
||||
if (tri && tri[i] >= 0) flag = 1;
|
||||
if (!flag)
|
||||
error->one(FLERR,"Compute erotate/asphere requires extended particles");
|
||||
}
|
||||
|
||||
pfactor = 0.5 * force->mvv2e;
|
||||
}
|
||||
|
|
@ -66,8 +79,16 @@ double ComputeERotateAsphere::compute_scalar()
|
|||
{
|
||||
invoked_scalar = update->ntimestep;
|
||||
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
AtomVecEllipsoid::Bonus *ebonus;
|
||||
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
|
||||
AtomVecLine::Bonus *lbonus;
|
||||
if (avec_line) lbonus = avec_line->bonus;
|
||||
AtomVecTri::Bonus *tbonus;
|
||||
if (avec_tri) tbonus = avec_tri->bonus;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *line = atom->line;
|
||||
int *tri = atom->tri;
|
||||
double **omega = atom->omega;
|
||||
double **angmom = atom->angmom;
|
||||
double *rmass = atom->rmass;
|
||||
int *mask = atom->mask;
|
||||
|
|
@ -76,6 +97,7 @@ double ComputeERotateAsphere::compute_scalar()
|
|||
// sum rotational energy for each particle
|
||||
// no point particles since divide by inertia
|
||||
|
||||
double length;
|
||||
double *shape,*quat;
|
||||
double wbody[3],inertia[3];
|
||||
double rot[3][3];
|
||||
|
|
@ -83,28 +105,54 @@ double ComputeERotateAsphere::compute_scalar()
|
|||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (ellipsoid && ellipsoid[i] >= 0) {
|
||||
shape = ebonus[ellipsoid[i]].shape;
|
||||
quat = ebonus[ellipsoid[i]].quat;
|
||||
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
// principal moments of inertia
|
||||
|
||||
inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
|
||||
MathExtra::quat_to_mat(quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
|
||||
erotate += inertia[0]*wbody[0]*wbody[0] +
|
||||
inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
|
||||
|
||||
// principal moments of inertia
|
||||
} else if (line && line[i] >= 0) {
|
||||
length = lbonus[line[i]].length;
|
||||
|
||||
inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
|
||||
erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
|
||||
omega[i][2]*omega[i][2]) * length*length*rmass[i] / 12.0;
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
} else if (tri && tri[i] >= 0) {
|
||||
|
||||
MathExtra::quat_to_mat(quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
|
||||
erotate += inertia[0]*wbody[0]*wbody[0] +
|
||||
inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
|
||||
// principal moments of inertia
|
||||
|
||||
inertia[0] = tbonus[tri[i]].inertia[0];
|
||||
inertia[1] = tbonus[tri[i]].inertia[1];
|
||||
inertia[2] = tbonus[tri[i]].inertia[2];
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
|
||||
MathExtra::quat_to_mat(tbonus[tri[i]].quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
|
||||
erotate += inertia[0]*wbody[0]*wbody[0] +
|
||||
inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
scalar *= pfactor;
|
||||
return scalar;
|
||||
|
|
|
|||
|
|
@ -32,7 +32,9 @@ class ComputeERotateAsphere : public Compute {
|
|||
|
||||
private:
|
||||
double pfactor;
|
||||
class AtomVecEllipsoid *avec;
|
||||
class AtomVecEllipsoid *avec_ellipsoid;
|
||||
class AtomVecLine *avec_line;
|
||||
class AtomVecTri *avec_tri;
|
||||
};
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,162 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "stdio.h"
|
||||
#include "string.h"
|
||||
#include "fix_nve_line.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec_line.h"
|
||||
#include "domain.h"
|
||||
#include "math_const.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
#define INERTIA (1.0/12.0) // moment of inertia prefactor for line segment
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixNVELine::FixNVELine(LAMMPS *lmp, int narg, char **arg) :
|
||||
FixNVE(lmp, narg, arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Illegal fix nve/line command");
|
||||
|
||||
time_integrate = 1;
|
||||
|
||||
MINUSPI = -MY_PI;
|
||||
TWOPI = 2.0*MY_PI;
|
||||
|
||||
// error checks
|
||||
|
||||
avec = (AtomVecLine *) atom->style_match("line");
|
||||
if (!avec) error->all(FLERR,"Fix nve/line requires atom style line");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixNVELine::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= FINAL_INTEGRATE;
|
||||
mask |= INITIAL_INTEGRATE_RESPA;
|
||||
mask |= FINAL_INTEGRATE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNVELine::init()
|
||||
{
|
||||
int i,itype;
|
||||
|
||||
if (domain->dimension != 2)
|
||||
error->all(FLERR,"Fix nve/line can only be used for 2d simulations");
|
||||
|
||||
// check that all particles are line segments
|
||||
// no point particles allowed
|
||||
|
||||
int *line = atom->line;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (line[i] < 0) error->one(FLERR,"Fix nve/line requires line particles");
|
||||
}
|
||||
|
||||
FixNVE::init();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNVELine::initial_integrate(int vflag)
|
||||
{
|
||||
double dtfm,dtirotate,delx,dely,length,theta;
|
||||
|
||||
AtomVecLine::Bonus *bonus = avec->bonus;
|
||||
int *line = atom->line;
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double **omega = atom->omega;
|
||||
double **torque = atom->torque;
|
||||
double *rmass = atom->rmass;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
// set timestep here since dt may have changed or come via rRESPA
|
||||
|
||||
double dtfrotate = dtf / INERTIA;
|
||||
|
||||
// update v,x,omega,theta for all particles
|
||||
// d_omega/dt = torque / inertia
|
||||
// bound theta by -PI to PI
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / rmass[i];
|
||||
v[i][0] += dtfm * f[i][0];
|
||||
v[i][1] += dtfm * f[i][1];
|
||||
x[i][0] += dtv * v[i][0];
|
||||
x[i][1] += dtv * v[i][1];
|
||||
|
||||
length = bonus[line[i]].length;
|
||||
theta = bonus[line[i]].theta;
|
||||
dtirotate = dtfrotate / (length*length*rmass[i]);
|
||||
omega[i][2] += dtirotate * torque[i][2];
|
||||
theta += dtv * omega[i][2];
|
||||
while (theta <= MINUSPI) theta += TWOPI;
|
||||
while (theta > MY_PI) theta -= TWOPI;
|
||||
bonus[line[i]].theta = theta;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNVELine::final_integrate()
|
||||
{
|
||||
double dtfm,dtirotate,length;
|
||||
|
||||
AtomVecLine::Bonus *bonus = avec->bonus;
|
||||
int *line = atom->line;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double **omega = atom->omega;
|
||||
double **torque = atom->torque;
|
||||
double *rmass = atom->rmass;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
// set timestep here since dt may have changed or come via rRESPA
|
||||
|
||||
double dtfrotate = dtf / INERTIA;
|
||||
|
||||
// update v,omega for all particles
|
||||
// d_omega/dt = torque / inertia
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / rmass[i];
|
||||
v[i][0] += dtfm * f[i][0];
|
||||
v[i][1] += dtfm * f[i][1];
|
||||
|
||||
length = bonus[line[i]].length;
|
||||
dtirotate = dtfrotate / (length*length*rmass[i]);
|
||||
omega[i][2] += dtirotate * torque[i][2];
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(nve/line,FixNVELine)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_NVE_LINE_H
|
||||
#define LMP_FIX_NVE_LINE_H
|
||||
|
||||
#include "fix_nve.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixNVELine : public FixNVE {
|
||||
public:
|
||||
FixNVELine(class LAMMPS *, int, char **);
|
||||
~FixNVELine() {}
|
||||
int setmask();
|
||||
void init();
|
||||
void initial_integrate(int);
|
||||
void final_integrate();
|
||||
|
||||
private:
|
||||
double MINUSPI,TWOPI;
|
||||
class AtomVecLine *avec;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,156 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "stdio.h"
|
||||
#include "string.h"
|
||||
#include "fix_nve_tri.h"
|
||||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec_tri.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixNVETri::FixNVETri(LAMMPS *lmp, int narg, char **arg) :
|
||||
FixNVE(lmp, narg, arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Illegal fix nve/tri command");
|
||||
|
||||
time_integrate = 1;
|
||||
|
||||
// error checks
|
||||
|
||||
avec = (AtomVecTri *) atom->style_match("tri");
|
||||
if (!avec) error->all(FLERR,"Fix nve/tri requires atom style tri");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixNVETri::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= FINAL_INTEGRATE;
|
||||
mask |= INITIAL_INTEGRATE_RESPA;
|
||||
mask |= FINAL_INTEGRATE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNVETri::init()
|
||||
{
|
||||
int i,itype;
|
||||
|
||||
if (domain->dimension != 3)
|
||||
error->all(FLERR,"Fix nve/line can only be used for 3d simulations");
|
||||
|
||||
// check that all particles are triangles
|
||||
// no point particles allowed
|
||||
|
||||
int *tri = atom->tri;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (tri[i] < 0) error->one(FLERR,"Fix nve/tri requires tri particles");
|
||||
}
|
||||
|
||||
FixNVE::init();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNVETri::initial_integrate(int vflag)
|
||||
{
|
||||
double dtfm;
|
||||
double omega[3];
|
||||
|
||||
AtomVecTri::Bonus *bonus = avec->bonus;
|
||||
int *tri = atom->tri;
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double **angmom = atom->angmom;
|
||||
double **torque = atom->torque;
|
||||
double *rmass = atom->rmass;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
// set timestep here since dt may have changed or come via rRESPA
|
||||
|
||||
dtq = 0.5 * dtv;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / rmass[i];
|
||||
v[i][0] += dtfm * f[i][0];
|
||||
v[i][1] += dtfm * f[i][1];
|
||||
v[i][2] += dtfm * f[i][2];
|
||||
x[i][0] += dtv * v[i][0];
|
||||
x[i][1] += dtv * v[i][1];
|
||||
x[i][2] += dtv * v[i][2];
|
||||
|
||||
// update angular momentum by 1/2 step
|
||||
|
||||
angmom[i][0] += dtf * torque[i][0];
|
||||
angmom[i][1] += dtf * torque[i][1];
|
||||
angmom[i][2] += dtf * torque[i][2];
|
||||
|
||||
// compute omega at 1/2 step from angmom at 1/2 step and current q
|
||||
// update quaternion a full step via Richardson iteration
|
||||
// returns new normalized quaternion
|
||||
|
||||
MathExtra::mq_to_omega(angmom[i],bonus[tri[i]].quat,
|
||||
bonus[tri[i]].inertia,omega);
|
||||
MathExtra::richardson(bonus[tri[i]].quat,angmom[i],omega,
|
||||
bonus[tri[i]].inertia,dtq);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNVETri::final_integrate()
|
||||
{
|
||||
double dtfm;
|
||||
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double **angmom = atom->angmom;
|
||||
double **torque = atom->torque;
|
||||
double *rmass = atom->rmass;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
// update v,omega for all particles
|
||||
// d_omega/dt = torque / inertia
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf / rmass[i];
|
||||
v[i][0] += dtfm * f[i][0];
|
||||
v[i][1] += dtfm * f[i][1];
|
||||
v[i][2] += dtfm * f[i][2];
|
||||
|
||||
angmom[i][0] += dtf * torque[i][0];
|
||||
angmom[i][1] += dtf * torque[i][1];
|
||||
angmom[i][2] += dtf * torque[i][2];
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(nve/tri,FixNVETri)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_NVE_TRI_H
|
||||
#define LMP_FIX_NVE_TRI_H
|
||||
|
||||
#include "fix_nve.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixNVETri : public FixNVE {
|
||||
public:
|
||||
FixNVETri(class LAMMPS *, int, char **);
|
||||
~FixNVETri() {}
|
||||
int setmask();
|
||||
void init();
|
||||
void initial_integrate(int);
|
||||
void final_integrate();
|
||||
|
||||
private:
|
||||
double dtq;
|
||||
class AtomVecTri *avec;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
Loading…
Reference in New Issue