lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2011-10-24 17:45:39 +00:00
parent fa80bceb30
commit cb5b707c74
13 changed files with 196 additions and 615 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/ASPHERE/Install.sh
View File

@ -8,6 +8,7 @@ if (test $1 = 1) then
cp fix_nph_asphere.cpp ..
cp fix_npt_asphere.cpp ..
cp fix_nve_asphere.cpp ..
cp fix_nve_asphere_noforce.cpp ..
cp fix_nve_line.cpp ..
cp fix_nve_tri.cpp ..
cp fix_nvt_asphere.cpp ..
@ -22,6 +23,7 @@ if (test $1 = 1) then
cp fix_nph_asphere.h ..
cp fix_npt_asphere.h ..
cp fix_nve_asphere.h ..
cp fix_nve_asphere_noforce.h ..
cp fix_nve_line.h ..
cp fix_nve_tri.h ..
cp fix_nvt_asphere.h ..
@ -38,6 +40,7 @@ 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_asphere_noforce.cpp
rm -f ../fix_nve_line.cpp
rm -f ../fix_nve_tri.cpp
rm -f ../fix_nvt_asphere.cpp
@ -52,6 +55,7 @@ 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_asphere_noforce.h
rm -f ../fix_nve_line.h
rm -f ../fix_nve_tri.h
rm -f ../fix_nvt_asphere.h

2
src/ASPHERE/fix_nve_asphere.cpp
View File

@ -45,7 +45,7 @@ FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) :
void FixNVEAsphere::init()
{
// check that all particles are finite-size
// check that all particles are finite-size ellipsoids
// no point particles allowed, spherical is OK
int *ellipsoid = atom->ellipsoid;

109
src/ASPHERE/fix_nve_asphere_noforce.cpp Normal file
View File

@ -0,0 +1,109 @@
/* ----------------------------------------------------------------------
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_nve_asphere_noforce.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "group.h"
#include "update.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixNVEAsphereNoforce::FixNVEAsphereNoforce(LAMMPS *lmp, int narg, char **arg) :
FixNVENoforce(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal fix nve/asphere/noforce command");
time_integrate = 1;
// error check
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!atom->ellipsoid_flag)
error->all(FLERR,"Fix nve/asphere/noforce requires atom style ellipsoid");
}
/* ---------------------------------------------------------------------- */
void FixNVEAsphereNoforce::init()
{
FixNVENoforce::init();
dtq = 0.5 * dtv;
// check that all particles are finite-size ellipsoids
// no point particles allowed, spherical is OK
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (ellipsoid[i] < 0)
error->one(FLERR,"Fix nve/asphere/noforce requires extended particles");
}
/* ---------------------------------------------------------------------- */
void FixNVEAsphereNoforce::initial_integrate(int vflag)
{
AtomVecEllipsoid::Bonus *bonus;
if (avec) bonus = avec->bonus;
double **x = atom->x;
double **v = atom->v;
double **angmom = atom->angmom;
double *rmass = atom->rmass;
int *ellipsoid = atom->ellipsoid;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double *shape,*quat;
double inertia[3],omega[3];
// update positions and quaternions for all particles
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
// principal moments of inertia
shape = bonus[ellipsoid[i]].shape;
quat = bonus[ellipsoid[i]].quat;
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;
// 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],quat,inertia,omega);
MathExtra::richardson(quat,angmom[i],omega,inertia,dtq);
}
}
}

41
src/ASPHERE/fix_nve_asphere_noforce.h Executable file
View File

@ -0,0 +1,41 @@
/* ----------------------------------------------------------------------
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/asphere/noforce,FixNVEAsphereNoforce)
#else
#ifndef LMP_FIX_NVE_ASPHERE_NOFORCE_H
#define LMP_FIX_NVE_ASPHERE_NOFORCE_H
#include "fix_nve_noforce.h"
namespace LAMMPS_NS {
class FixNVEAsphereNoforce : public FixNVENoforce {
public:
FixNVEAsphereNoforce(class LAMMPS *, int, char **);
void initial_integrate(int);
void init();
private:
double dtq;
class AtomVecEllipsoid *avec;
};
}
#endif
#endif

4
src/COLLOID/Install.sh
View File

@ -4,24 +4,20 @@ if (test $1 = 1) then
cp fix_wall_colloid.cpp ..
cp pair_colloid.cpp ..
cp pair_lubricate.cpp ..
cp pair_yukawa_colloid.cpp ..
cp fix_wall_colloid.h ..
cp pair_colloid.h ..
cp pair_lubricate.h ..
cp pair_yukawa_colloid.h ..
elif (test $1 = 0) then
rm -f ../fix_wall_colloid.cpp
rm -f ../pair_colloid.cpp
rm -f ../pair_lubricate.cpp
rm -f ../pair_yukawa_colloid.cpp
rm -f ../fix_wall_colloid.h
rm -f ../pair_colloid.h
rm -f ../pair_lubricate.h
rm -f ../pair_yukawa_colloid.h
fi

523
src/COLLOID/pair_lubricate.cpp
View File

@ -1,523 +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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Randy Schunk (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_lubricate.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "random_mars.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairLubricate::PairLubricate(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
random = NULL;
}
/* ---------------------------------------------------------------------- */
PairLubricate::~PairLubricate()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(cut_inner);
}
delete random;
}
/* ---------------------------------------------------------------------- */
void PairLubricate::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz,tx,ty,tz;
double rsq,r,h_sep,radi,tfmag;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
double vt1,vt2,vt3,w1,w2,w3,v_shear1,v_shear2,v_shear3;
double omega_t_1,omega_t_2,omega_t_3;
double n_cross_omega_t_1,n_cross_omega_t_2,n_cross_omega_t_3;
double wr1,wr2,wr3,wnnr,wn1,wn2,wn3;
double P_dot_wrel_1,P_dot_wrel_2,P_dot_wrel_3;
double a_squeeze,a_shear,a_pump,a_twist;
int *ilist,*jlist,*numneigh,**firstneigh;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **omega = atom->omega;
double **torque = atom->torque;
double *radius = atom->radius;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double vxmu2f = force->vxmu2f;
double prethermostat = sqrt(2.0 * force->boltz * t_target / update->dt);
prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
a_squeeze = a_shear = a_pump = a_twist = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
radi = radius[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
// relative translational velocity
vr1 = v[i][0] - v[j][0];
vr2 = v[i][1] - v[j][1];
vr3 = v[i][2] - v[j][2];
// normal component N.(v1-v2) = nn.(v1-v2)
vnnr = vr1*delx + vr2*dely + vr3*delz;
vnnr /= r;
vn1 = delx*vnnr / r;
vn2 = dely*vnnr / r;
vn3 = delz*vnnr / r;
// tangential component -P.(v1-v2)
// P = (I - nn) where n is vector between centers
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// additive rotational velocity = omega_1 + omega_2
w1 = omega[i][0] + omega[j][0];
w2 = omega[i][1] + omega[j][1];
w3 = omega[i][2] + omega[j][2];
// relative velocities n X P . (v1-v2) = n X (I-nn) . (v1-v2)
v_shear1 = (dely*vt3 - delz*vt2) / r;
v_shear2 = -(delx*vt3 - delz*vt1) / r;
v_shear3 = (delx*vt2 - dely*vt1) / r;
// relative rotation rate P.(omega1 + omega2)
omega_t_1 = w1 - delx*(delx*w1) / rsq;
omega_t_2 = w2 - dely*(dely*w2) / rsq;
omega_t_3 = w3 - delz*(delz*w3) / rsq;
// n X omega_t
n_cross_omega_t_1 = (dely*omega_t_3 - delz*omega_t_2) / r;
n_cross_omega_t_2 = -(delx*omega_t_3 - delz*omega_t_1) / r;
n_cross_omega_t_3 = (delx*omega_t_2 - dely*omega_t_1) / r;
// N.(w1-w2) and P.(w1-w2)
wr1 = omega[i][0] - omega[j][0];
wr2 = omega[i][1] - omega[j][1];
wr3 = omega[i][2] - omega[j][2];
wnnr = wr1*delx + wr2*dely + wr3*delz;
wn1 = delx*wnnr / rsq;
wn2 = dely*wnnr / rsq;
wn3 = delz*wnnr / rsq;
P_dot_wrel_1 = wr1 - delx*(delx*wr1)/rsq;
P_dot_wrel_2 = wr2 - dely*(dely*wr2)/rsq;
P_dot_wrel_3 = wr3 - delz*(delz*wr3)/rsq;
// compute components of pair-hydro
h_sep = r - 2.0*radi;
if (flag1)
a_squeeze = (3.0*MY_PI*mu*2.0*radi/2.0) * (2.0*radi/4.0/h_sep);
if (flag2)
a_shear = (MY_PI*mu*2.*radi/2.0) *
log(2.0*radi/2.0/h_sep)*(2.0*radi+h_sep)*(2.0*radi+h_sep)/4.0;
if (flag3)
a_pump = (MY_PI*mu*pow(2.0*radi,4)/8.0) *
((3.0/20.0) * log(2.0*radi/2.0/h_sep) +
(63.0/250.0) * (h_sep/2.0/radi) * log(2.0*radi/2.0/h_sep));
if (flag4)
a_twist = (MY_PI*mu*pow(2.0*radi,4)/4.0) *
(h_sep/2.0/radi) * log(2.0/(2.0*h_sep));
if (h_sep >= cut_inner[itype][jtype]) {
fx = -a_squeeze*vn1 - a_shear*(2.0/r)*(2.0/r)*vt1 +
(2.0/r)*a_shear*n_cross_omega_t_1;
fy = -a_squeeze*vn2 - a_shear*(2.0/r)*(2.0/r)*vt2 +
(2.0/r)*a_shear*n_cross_omega_t_2;
fz = -a_squeeze*vn3 - a_shear*(2.0/r)*(2.0/r)*vt3 +
(2.0/r)*a_shear*n_cross_omega_t_3;
fx *= vxmu2f;
fy *= vxmu2f;
fz *= vxmu2f;
// add in thermostat force
tfmag = prethermostat*sqrt(a_squeeze)*(random->uniform()-0.5);
fx -= tfmag * delx/r;
fy -= tfmag * dely/r;
fz -= tfmag * delz/r;
tx = -(2.0/r)*a_shear*v_shear1 - a_shear*omega_t_1 -
a_pump*P_dot_wrel_1 - a_twist*wn1;
ty = -(2.0/r)*a_shear*v_shear2 - a_shear*omega_t_2 -
a_pump*P_dot_wrel_2 - a_twist*wn2;
tz = -(2.0/r)*a_shear*v_shear3 - a_shear*omega_t_3 -
a_pump*P_dot_wrel_3 - a_twist*wn3;
torque[i][0] += vxmu2f * tx;
torque[i][1] += vxmu2f * ty;
torque[i][2] += vxmu2f * tz;
} else {
a_squeeze = (3.0*MY_PI*mu*2.0*radi/2.0) *
(2.0*radi/4.0/cut_inner[itype][jtype]);
fpair = -a_squeeze*vnnr;
fpair *= vxmu2f;
// add in thermostat force
fpair -= prethermostat*sqrt(a_squeeze)*(random->uniform()-0.5);
fx = fpair * delx/r;
fy = fpair * dely/r;
fz = fpair * delz/r;
}
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
if (h_sep >= cut_inner[itype][jtype]) {
tx = -(2.0/r)*a_shear*v_shear1 - a_shear*omega_t_1 +
a_pump*P_dot_wrel_1 + a_twist*wn1;
ty = -(2.0/r)*a_shear*v_shear2 - a_shear*omega_t_2 +
a_pump*P_dot_wrel_2 + a_twist*wn2;
tz = -(2.0/r)*a_shear*v_shear3 - a_shear*omega_t_3 +
a_pump*P_dot_wrel_3 + a_twist*wn3;
torque[j][0] += vxmu2f * tx;
torque[j][1] += vxmu2f * ty;
torque[j][2] += vxmu2f * tz;
}
}
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
0.0,0.0,fx,fy,fz,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLubricate::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(cut_inner,n+1,n+1,"pair:cut_inner");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLubricate::settings(int narg, char **arg)
{
if (narg != 9) error->all(FLERR,"Illegal pair_style command");
mu = force->numeric(arg[0]);
flag1 = force->inumeric(arg[1]);
flag2 = force->inumeric(arg[2]);
flag3 = force->inumeric(arg[3]);
flag4 = force->inumeric(arg[4]);
cut_inner_global = force->numeric(arg[5]);
cut_global = force->numeric(arg[6]);
t_target = force->numeric(arg[7]);
seed = force->inumeric(arg[8]);
// initialize Marsaglia RNG with processor-unique seed
if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
delete random;
random = new RanMars(lmp,seed + comm->me);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) {
cut_inner[i][j] = cut_inner_global;
cut[i][j] = cut_global;
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLubricate::coeff(int narg, char **arg)
{
if (narg != 2 && narg != 4)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double cut_inner_one = cut_inner_global;
double cut_one = cut_global;
if (narg == 4) {
cut_inner_one = force->numeric(arg[2]);
cut_one = force->numeric(arg[3]);
}
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut_inner[i][j] = cut_inner_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLubricate::init_style()
{
if (!atom->sphere_flag)
error->all(FLERR,"Pair lubricate requires atom style sphere");
if (comm->ghost_velocity == 0)
error->all(FLERR,"Pair lubricate requires ghost atoms store velocity");
neighbor->request(this);
// require that atom radii are identical within each type
// require monodisperse system with same radii for all types
double rad,radtype;
for (int i = 1; i <= atom->ntypes; i++) {
if (!atom->radius_consistency(i,radtype))
error->all(FLERR,"Pair lubricate requires monodisperse particles");
if (i > 1 && radtype != rad)
error->all(FLERR,"Pair lubricate requires monodisperse particles");
rad = radtype;
}
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLubricate::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
}
cut_inner[j][i] = cut_inner[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLubricate::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&cut_inner[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLubricate::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&cut_inner[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLubricate::write_restart_settings(FILE *fp)
{
fwrite(&mu,sizeof(double),1,fp);
fwrite(&flag1,sizeof(int),1,fp);
fwrite(&flag2,sizeof(int),1,fp);
fwrite(&flag3,sizeof(int),1,fp);
fwrite(&flag4,sizeof(int),1,fp);
fwrite(&cut_inner_global,sizeof(double),1,fp);
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&t_target,sizeof(double),1,fp);
fwrite(&seed,sizeof(int),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLubricate::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&mu,sizeof(double),1,fp);
fread(&flag1,sizeof(int),1,fp);
fread(&flag2,sizeof(int),1,fp);
fread(&flag3,sizeof(int),1,fp);
fread(&flag4,sizeof(int),1,fp);
fread(&cut_inner_global,sizeof(double),1,fp);
fread(&cut_global,sizeof(double),1,fp);
fread(&t_target, sizeof(double),1,fp);
fread(&seed, sizeof(int),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&mu,1,MPI_DOUBLE,0,world);
MPI_Bcast(&flag1,1,MPI_INT,0,world);
MPI_Bcast(&flag2,1,MPI_INT,0,world);
MPI_Bcast(&flag3,1,MPI_INT,0,world);
MPI_Bcast(&flag4,1,MPI_INT,0,world);
MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&t_target,1,MPI_DOUBLE,0,world);
MPI_Bcast(&seed,1,MPI_INT,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
// additional setup based on restart parameters
delete random;
random = new RanMars(lmp,seed + comm->me);
}
/* ---------------------------------------------------------------------- */
void *PairLubricate::extract(char *str, int &dim)
{
dim = 0;
if (strcmp(str,"mu") == 0) return (void *) &mu;
return NULL;
}

57
src/COLLOID/pair_lubricate.h
View File

@ -1,57 +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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(lubricate,PairLubricate)
#else
#ifndef LMP_PAIR_LUBRICATE_H
#define LMP_PAIR_LUBRICATE_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLubricate : public Pair {
public:
PairLubricate(class LAMMPS *);
virtual ~PairLubricate();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void *extract(char *, int &);
protected:
double cut_inner_global,cut_global;
double t_target,mu;
int flag1,flag2,flag3,flag4;
int seed;
double **cut_inner,**cut;
class RanMars *random;
void allocate();
};
}
#endif
#endif

19
src/domain.h
View File

@ -93,15 +93,16 @@ class Domain : protected Pointers {
virtual void set_local_box();
virtual void reset_box();
virtual void pbc();
void remap(double *, int &);
void remap(double *);
void remap_near(double *, double *);
void unmap(double *, int);
void unmap(double *, int, double *);
int minimum_image_check(double, double, double);
void minimum_image(double &, double &, double &);
void minimum_image(double *);
void closest_image(double *, double *, double *);
virtual void remap(double *, int &);
virtual void remap(double *);
virtual void remap_near(double *, double *);
virtual void unmap(double *, int);
virtual void unmap(double *, int, double *);
virtual int minimum_image_check(double, double, double);
virtual void minimum_image(double &, double &, double &);
virtual void minimum_image(double *);
virtual void closest_image(double *, double *, double *);
void set_lattice(int, char **);
void add_region(int, char **);
void delete_region(int, char **);

6
src/fix_nve_noforce.h
View File

@ -28,12 +28,12 @@ class FixNVENoforce : public Fix {
public:
FixNVENoforce(class LAMMPS *, int, char **);
int setmask();
void init();
void initial_integrate(int);
virtual void init();
virtual void initial_integrate(int);
void initial_integrate_respa(int, int, int);
void reset_dt();
private:
protected:
double dtv;
double *step_respa;
};

18
src/fix_nve_sphere.cpp
View File

@ -61,21 +61,11 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) :
/* ---------------------------------------------------------------------- */
int FixNVESphere::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNVESphere::init()
{
// check that all particles are finite-size
FixNVE::init();
// check that all particles are finite-size spheres
// no point particles allowed
double *radius = atom->radius;
@ -86,8 +76,6 @@ void FixNVESphere::init()
if (mask[i] & groupbit)
if (radius[i] == 0.0)
error->one(FLERR,"Fix nve/sphere requires extended particles");
FixNVE::init();
}
/* ---------------------------------------------------------------------- */

1
src/fix_nve_sphere.h
View File

@ -28,7 +28,6 @@ class FixNVESphere : public FixNVE {
public:
FixNVESphere(class LAMMPS *, int, char **);
virtual ~FixNVESphere() {}
int setmask();
void init();
virtual void initial_integrate(int);
virtual void final_integrate();

25
src/fix_wall.cpp
View File

@ -45,6 +45,7 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
nwall = 0;
int scaleflag = 1;
fldflag = 0;
int iarg = 3;
while (iarg < narg) {
@ -94,6 +95,12 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
else error->all(FLERR,"Illegal fix wall command");
iarg += 2;
} else if (strcmp(arg[iarg],"fld") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall command");
if (strcmp(arg[iarg+1],"no") == 0) fldflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) fldflag = 1;
else error->all(FLERR,"Illegal fix wall command");
iarg += 2;
} else error->all(FLERR,"Illegal fix wall command");
}
@ -169,7 +176,12 @@ FixWall::~FixWall()
int FixWall::setmask()
{
int mask = 0;
mask |= POST_FORCE;
// FLD implicit needs to invoke wall forces before pair style
if (fldflag) mask != PRE_FORCE;
else mask |= POST_FORCE;
mask |= THERMO_ENERGY;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
@ -204,7 +216,7 @@ void FixWall::init()
void FixWall::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
if (!fldflag) post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
@ -219,6 +231,15 @@ void FixWall::min_setup(int vflag)
post_force(vflag);
}
/* ----------------------------------------------------------------------
only called if fldflag set, in place of post_force
------------------------------------------------------------------------- */
void FixWall::pre_force(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixWall::post_force(int vflag)

2
src/fix_wall.h
View File

@ -26,6 +26,7 @@ class FixWall : public Fix {
virtual void init();
void setup(int);
void min_setup(int);
void pre_force(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
@ -45,6 +46,7 @@ class FixWall : public Fix {
double ewall[7],ewall_all[7];
int nlevels_respa;
double dt;
int fldflag;
};
}