diff --git a/src/ASPHERE/Install.sh b/src/ASPHERE/Install.sh index 77a640197c..aca41e02b9 100644 --- a/src/ASPHERE/Install.sh +++ b/src/ASPHERE/Install.sh @@ -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 diff --git a/src/ASPHERE/fix_nve_asphere.cpp b/src/ASPHERE/fix_nve_asphere.cpp index 17ae8950e1..ac3abea6bc 100755 --- a/src/ASPHERE/fix_nve_asphere.cpp +++ b/src/ASPHERE/fix_nve_asphere.cpp @@ -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; diff --git a/src/ASPHERE/fix_nve_asphere_noforce.cpp b/src/ASPHERE/fix_nve_asphere_noforce.cpp new file mode 100644 index 0000000000..42793723b5 --- /dev/null +++ b/src/ASPHERE/fix_nve_asphere_noforce.cpp @@ -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); + } + } +} diff --git a/src/ASPHERE/fix_nve_asphere_noforce.h b/src/ASPHERE/fix_nve_asphere_noforce.h new file mode 100755 index 0000000000..c1708ff7bb --- /dev/null +++ b/src/ASPHERE/fix_nve_asphere_noforce.h @@ -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 diff --git a/src/COLLOID/Install.sh b/src/COLLOID/Install.sh index 4286a636f8..d5bac42feb 100644 --- a/src/COLLOID/Install.sh +++ b/src/COLLOID/Install.sh @@ -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 diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp deleted file mode 100644 index f3c9a0dc06..0000000000 --- a/src/COLLOID/pair_lubricate.cpp +++ /dev/null @@ -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 *) μ - return NULL; -} diff --git a/src/COLLOID/pair_lubricate.h b/src/COLLOID/pair_lubricate.h deleted file mode 100644 index 30c9f71baf..0000000000 --- a/src/COLLOID/pair_lubricate.h +++ /dev/null @@ -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 diff --git a/src/domain.h b/src/domain.h index a6fee60963..c07808ca02 100644 --- a/src/domain.h +++ b/src/domain.h @@ -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 **); diff --git a/src/fix_nve_noforce.h b/src/fix_nve_noforce.h index d56ecd96a1..90989869b8 100644 --- a/src/fix_nve_noforce.h +++ b/src/fix_nve_noforce.h @@ -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; }; diff --git a/src/fix_nve_sphere.cpp b/src/fix_nve_sphere.cpp index c81452994c..a00fa0bdf1 100644 --- a/src/fix_nve_sphere.cpp +++ b/src/fix_nve_sphere.cpp @@ -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(); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_nve_sphere.h b/src/fix_nve_sphere.h index 28fc47c591..c3148c2c54 100644 --- a/src/fix_nve_sphere.h +++ b/src/fix_nve_sphere.h @@ -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(); diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index 3cc68017cc..2be0651690 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -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) diff --git a/src/fix_wall.h b/src/fix_wall.h index 8c0c39d024..62ee69bc54 100644 --- a/src/fix_wall.h +++ b/src/fix_wall.h @@ -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; }; }