ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

sync with GH

Browse Source 
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15746 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
2016-10-11 18:43:51 +00:00
parent 7a4da54a71
commit e710053de6
14 changed files with 205 additions and 130 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/utils/converters/lammpsdoc/txt2rst.py
View File

@ -104,7 +104,7 @@ class RSTMarkup(Markup):
anchor_pos = href.find('#')
if anchor_pos >= 0:
if anchor_pos >= 0 and not href.startswith('http'):
href = href[anchor_pos+1:]
return ":ref:`%s <%s>`" % (content, href)

5
doc/utils/converters/tests/test_txt2rst.py
View File

@ -424,6 +424,11 @@ class TestSpecialCommands(unittest.TestCase):
"one \n\n"
"a :ref:`link <name>` to above\n\n", s)
def test_external_anchor_link(self):
s = self.txt2rst.convert('some text "containing a\n'
'link"_http://lammps.sandia.gov/movies.html#granregion with an anchor')
self.assertEqual('some text `containing a link <http://lammps.sandia.gov/movies.html#granregion>`_ with an anchor\n\n', s)
def test_define_link_alias(self):
s = self.txt2rst.convert("one :link(alias,value)\n"
"\"test\"_alias\n")

2
lib/colvars/colvarmodule.cpp
View File

@ -1557,5 +1557,5 @@ size_t const colvarmodule::cv_prec = 14;
size_t const colvarmodule::cv_width = 21;
size_t const colvarmodule::en_prec = 14;
size_t const colvarmodule::en_width = 21;
std::string const colvarmodule::line_marker =
const char * const colvarmodule::line_marker = (const char *)
"----------------------------------------------------------------------\n";

2
lib/colvars/colvarmodule.h
View File

@ -366,7 +366,7 @@ public:
/// Number of characters to represent the collective variables energy
static size_t const en_width;
/// Line separator in the log output
static std::string const line_marker;
static const char * const line_marker;
// proxy functions

3
src/CLASS2/improper_class2.cpp
View File

@ -819,7 +819,8 @@ void ImproperClass2::angleangle(int eflag, int vflag)
if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,
fabcd[0],fabcd[2],fabcd[3],
delxAB,delyAB,delzAB,delxBC,delyBC,delzBC,delxBD,delyBD,delzBD);
delxAB,delyAB,delzAB,delxBC,delyBC,delzBC,
delxBD-delxBC,delyBD-delyBC,delzBD-delzBC);
}
}

194
src/GRANULAR/fix_wall_gran.cpp
View File

@ -37,7 +37,7 @@ using namespace MathConst;
// XYZ PLANE need to be 0,1,2
enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION};
enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION};
enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,BONDED_HISTORY};
enum{NONE,CONSTANT,EQUAL};
@ -78,21 +78,21 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
kn = force->numeric(FLERR,arg[4]);
if (strcmp(arg[5],"NULL") == 0) kt = kn * 2.0/7.0;
else kt = force->numeric(FLERR,arg[5]);
gamman = force->numeric(FLERR,arg[6]);
if (strcmp(arg[7],"NULL") == 0) gammat = 0.5 * gamman;
else gammat = force->numeric(FLERR,arg[7]);
xmu = force->numeric(FLERR,arg[8]);
int dampflag = force->inumeric(FLERR,arg[9]);
if (dampflag == 0) gammat = 0.0;
if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 ||
xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1)
error->all(FLERR,"Illegal fix wall/gran command");
// convert Kn and Kt from pressure units to force/distance^2 if Hertzian
if (pairstyle == HERTZ_HISTORY) {
kn /= force->nktv2p;
kt /= force->nktv2p;
@ -119,7 +119,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
}
// wallstyle args
idregion = NULL;
if (strcmp(arg[iarg],"xplane") == 0) {
@ -160,12 +160,12 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
strcpy(idregion,arg[iarg+1]);
iarg += 2;
}
// optional args
wiggle = 0;
wshear = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"wiggle") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix wall/gran command");
@ -188,7 +188,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
iarg += 3;
} else error->all(FLERR,"Illegal fix wall/gran command");
}
if (wallstyle == XPLANE && domain->xperiodic)
error->all(FLERR,"Cannot use wall in periodic dimension");
if (wallstyle == YPLANE && domain->yperiodic)
@ -197,7 +197,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Cannot use wall in periodic dimension");
if (wallstyle == ZCYLINDER && (domain->xperiodic || domain->yperiodic))
error->all(FLERR,"Cannot use wall in periodic dimension");
if (wiggle && wshear)
error->all(FLERR,"Cannot wiggle and shear fix wall/gran");
if (wiggle && wallstyle == ZCYLINDER && axis != 2)
@ -212,15 +212,15 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Cannot wiggle or shear with fix wall/gran/region");
// setup oscillations
if (wiggle) omega = 2.0*MY_PI / period;
// perform initial allocation of atom-based arrays
// register with Atom class
if (pairstyle == BONDED_HISTORY) sheardim = 7;
else sheardim = 3;
shearone = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
@ -348,7 +348,7 @@ void FixWallGran::post_force(int vflag)
}
vwall[axis] = amplitude*omega*sin(arg);
} else if (wshear) vwall[axis] = vshear;
// loop over all my atoms
// rsq = distance from wall
// dx,dy,dz = signed distance from wall
@ -358,7 +358,7 @@ void FixWallGran::post_force(int vflag)
// compute force and torque on atom if close enough to wall
// via wall potential matched to pair potential
// set shear if pair potential stores history
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
@ -368,14 +368,14 @@ void FixWallGran::post_force(int vflag)
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
rwall = 0.0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dx = dy = dz = 0.0;
if (wallstyle == XPLANE) {
del1 = x[i][0] - wlo;
del2 = whi - x[i][0];
@ -401,7 +401,7 @@ void FixWallGran::post_force(int vflag)
dx = -delr/delxy * x[i][0];
dy = -delr/delxy * x[i][1];
// rwall = -2r_c if inside cylinder, 2r_c outside
rwall = 2*(1-2*(delxy < cylradius))*cylradius;
rwall = (delxy < cylradius) ? -2*cylradius : 2*cylradius;
if (wshear && axis != 2) {
vwall[0] += vshear * x[i][1]/delxy;
vwall[1] += -vshear * x[i][0]/delxy;
@ -409,12 +409,12 @@ void FixWallGran::post_force(int vflag)
}
}
}
rsq = dx*dx + dy*dy + dz*dz;
if (rsq > radius[i]*radius[i]) {
if (history)
for (j = 0; j < sheardim; j++)
for (j = 0; j < sheardim; j++)
shearone[i][j] = 0.0;
} else {
@ -467,66 +467,66 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz,
rsqinv = 1.0/rsq;
// relative translational velocity
vr1 = v[0] - vwall[0];
vr2 = v[1] - vwall[1];
vr3 = v[2] - vwall[2];
// normal component
vnnr = vr1*dx + vr2*dy + vr3*dz;
vn1 = dx*vnnr * rsqinv;
vn2 = dy*vnnr * rsqinv;
vn3 = dz*vnnr * rsqinv;
// tangential component
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// relative rotational velocity
wr1 = radius*omega[0] * rinv;
wr2 = radius*omega[1] * rinv;
wr3 = radius*omega[2] * rinv;
// normal forces = Hookian contact + normal velocity damping
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radius-r)*rinv - damp;
// relative velocities
vtr1 = vt1 - (dz*wr2-dy*wr3);
vtr2 = vt2 - (dx*wr3-dz*wr1);
vtr3 = vt3 - (dy*wr1-dx*wr2);
vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
vrel = sqrt(vrel);
// force normalization
fn = xmu * fabs(ccel*r);
fs = meff*gammat*vrel;
if (vrel != 0.0) ft = MIN(fn,fs) / vrel;
else ft = 0.0;
// tangential force due to tangential velocity damping
fs1 = -ft*vtr1;
fs2 = -ft*vtr2;
fs3 = -ft*vtr3;
// forces & torques
fx = dx*ccel + fs1;
fy = dy*ccel + fs2;
fz = dz*ccel + fs3;
f[0] += fx;
f[1] += fy;
f[2] += fz;
tor1 = rinv * (dy*fs3 - dz*fs2);
tor2 = rinv * (dz*fs1 - dx*fs3);
tor3 = rinv * (dx*fs2 - dy*fs1);
@ -546,60 +546,60 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
double wr1,wr2,wr3,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;
double shrmag,rsht,rinv,rsqinv;
r = sqrt(rsq);
rinv = 1.0/r;
rsqinv = 1.0/rsq;
// relative translational velocity
vr1 = v[0] - vwall[0];
vr2 = v[1] - vwall[1];
vr3 = v[2] - vwall[2];
// normal component
vnnr = vr1*dx + vr2*dy + vr3*dz;
vn1 = dx*vnnr * rsqinv;
vn2 = dy*vnnr * rsqinv;
vn3 = dz*vnnr * rsqinv;
// tangential component
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// relative rotational velocity
wr1 = radius*omega[0] * rinv;
wr2 = radius*omega[1] * rinv;
wr3 = radius*omega[2] * rinv;
// normal forces = Hookian contact + normal velocity damping
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radius-r)*rinv - damp;
// relative velocities
vtr1 = vt1 - (dz*wr2-dy*wr3);
vtr2 = vt2 - (dx*wr3-dz*wr1);
vtr3 = vt3 - (dy*wr1-dx*wr2);
vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
vrel = sqrt(vrel);
// shear history effects
if (shearupdate) {
shear[0] += vtr1*dt;
shear[1] += vtr2*dt;
shear[2] += vtr3*dt;
}
shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]);
// rotate shear displacements
rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz;
rsht = rsht*rsqinv;
if (shearupdate) {
@ -607,18 +607,18 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
shear[1] -= rsht*dy;
shear[2] -= rsht*dz;
}
// tangential forces = shear + tangential velocity damping
fs1 = - (kt*shear[0] + meff*gammat*vtr1);
fs2 = - (kt*shear[1] + meff*gammat*vtr2);
fs3 = - (kt*shear[2] + meff*gammat*vtr3);
// rescale frictional displacements and forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
fn = xmu * fabs(ccel*r);
if (fs > fn) {
if (shrmag != 0.0) {
shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) -
@ -632,17 +632,17 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
fs3 *= fn/fs;
} else fs1 = fs2 = fs3 = 0.0;
}
// forces & torques
fx = dx*ccel + fs1;
fy = dy*ccel + fs2;
fz = dz*ccel + fs3;
f[0] += fx;
f[1] += fy;
f[2] += fz;
tor1 = rinv * (dy*fs3 - dz*fs2);
tor2 = rinv * (dz*fs1 - dx*fs3);
tor3 = rinv * (dx*fs2 - dy*fs1);
@ -662,56 +662,56 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
double wr1,wr2,wr3,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;
double shrmag,rsht,polyhertz,rinv,rsqinv;
r = sqrt(rsq);
rinv = 1.0/r;
rsqinv = 1.0/rsq;
// relative translational velocity
vr1 = v[0] - vwall[0];
vr2 = v[1] - vwall[1];
vr3 = v[2] - vwall[2];
// normal component
vnnr = vr1*dx + vr2*dy + vr3*dz;
vn1 = dx*vnnr / rsq;
vn2 = dy*vnnr / rsq;
vn3 = dz*vnnr / rsq;
// tangential component
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// relative rotational velocity
wr1 = radius*omega[0] * rinv;
wr2 = radius*omega[1] * rinv;
wr3 = radius*omega[2] * rinv;
// normal forces = Hertzian contact + normal velocity damping
// rwall = 0 is flat wall case
// rwall positive or negative is curved wall
// will break (as it should) if rwall is negative and
// will break (as it should) if rwall is negative and
// its absolute value < radius of particle
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radius-r)*rinv - damp;
if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius);
else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius));
ccel *= polyhertz;
// relative velocities
vtr1 = vt1 - (dz*wr2-dy*wr3);
vtr2 = vt2 - (dx*wr3-dz*wr1);
vtr3 = vt3 - (dy*wr1-dx*wr2);
vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
vrel = sqrt(vrel);
// shear history effects
if (shearupdate) {
@ -720,9 +720,9 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
shear[2] += vtr3*dt;
}
shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]);
// rotate shear displacements
rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz;
rsht = rsht*rsqinv;
if (shearupdate) {
@ -730,18 +730,18 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
shear[1] -= rsht*dy;
shear[2] -= rsht*dz;
}
// tangential forces = shear + tangential velocity damping
fs1 = -polyhertz * (kt*shear[0] + meff*gammat*vtr1);
fs2 = -polyhertz * (kt*shear[1] + meff*gammat*vtr2);
fs3 = -polyhertz * (kt*shear[2] + meff*gammat*vtr3);
// rescale frictional displacements and forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
fn = xmu * fabs(ccel*r);
if (fs > fn) {
if (shrmag != 0.0) {
shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) -
@ -757,7 +757,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
}
// forces & torques
fx = dx*ccel + fs1;
fy = dy*ccel + fs2;
fz = dz*ccel + fs3;
@ -765,7 +765,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
f[0] += fx;
f[1] += fy;
f[2] += fz;
tor1 = rinv * (dy*fs3 - dz*fs2);
tor2 = rinv * (dz*fs1 - dx*fs3);
tor3 = rinv * (dx*fs2 - dy*fs1);
@ -864,8 +864,8 @@ void FixWallGran::bonded_history(double rsq, double dx, double dy, double dz,
// use Tsuji et al form
polyhertz = 1.2728- 4.2783*0.9 + 11.087*0.9*0.9 - 22.348*0.9*0.9*0.9 +
27.467*0.9*0.9*0.9*0.9 - 18.022*0.9*0.9*0.9*0.9*0.9 +
polyhertz = 1.2728- 4.2783*0.9 + 11.087*0.9*0.9 - 22.348*0.9*0.9*0.9 +
27.467*0.9*0.9*0.9*0.9 - 18.022*0.9*0.9*0.9*0.9*0.9 +
4.8218*0.9*0.9*0.9*0.9*0.9*0.9;
gammatsuji = 0.2*sqrt(meff*kn);
@ -969,7 +969,7 @@ void FixWallGran::bonded_history(double rsq, double dx, double dy, double dz,
shear[6] += magtwist*dt;
ktwist = 0.5*kt*(a0*aovera0)*(a0*aovera0);
magtortwist = -ktwist*shear[6] -
magtortwist = -ktwist*shear[6] -
0.5*polyhertz*gammatsuji*(a0*aovera0)*(a0*aovera0)*magtwist;
twistcrit=TWOTHIRDS*a0*aovera0*Fcrit;
@ -1113,11 +1113,11 @@ void FixWallGran::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++;
for (int i = 0; i < sheardim; i++)
shearone[nlocal][i] = extra[nlocal][m++];
}

18
src/MOLECULE/improper_umbrella.cpp
View File

@ -231,9 +231,25 @@ void ImproperUmbrella::compute(int eflag, int vflag)
f[i4][2] += f4[2]*a;
}
if (evflag)
if (evflag) {
// get correct 4-body geometry for virial tally
vb1x = x[i1][0] - x[i2][0];
vb1y = x[i1][1] - x[i2][1];
vb1z = x[i1][2] - x[i2][2];
vb2x = x[i3][0] - x[i2][0];
vb2y = x[i3][1] - x[i2][1];
vb2z = x[i3][2] - x[i2][2];
vb3x = x[i4][0] - x[i3][0];
vb3y = x[i4][1] - x[i3][1];
vb3z = x[i4][2] - x[i3][2];
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
}
}
}

4
src/Makefile
View File

@ -55,8 +55,8 @@ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
user-quip user-reaxc user-smd user-smtbq user-sph user-tally \
user-vtk
PACKLIB = compress gpu kim kokkos meam poems python reax voronoi \
user-atc user-awpmd user-colvars user-h5md user-molfile \
PACKLIB = compress gpu kim kokkos meam mpiio poems python reax voronoi \
user-atc user-awpmd user-colvars user-h5md user-lb user-molfile \
user-qmmm user-quip user-vtk
PACKALL = $(PACKAGE) $(PACKUSER)

2
src/USER-COLVARS/colvarproxy_lammps.cpp
View File

@ -193,7 +193,7 @@ double colvarproxy_lammps::compute()
previous_step = _lmp->update->ntimestep;
if (cvm::debug()) {
cvm::log(cvm::line_marker+
cvm::log(std::string(cvm::line_marker)+
"colvarproxy_lammps, step no. "+cvm::to_str(colvars->it)+"\n"+
"Updating internal data.\n");
}

31
src/USER-MISC/fix_flow_gauss.cpp
View File

@ -27,6 +27,7 @@
#include "domain.h"
#include "error.h"
#include "citeme.h"
#include "respa.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -63,6 +64,8 @@ FixFlowGauss::FixFlowGauss(LAMMPS *lmp, int narg, char **arg) :
extvector = 1;
size_vector = 3;
global_freq = 1; //data available every timestep
respa_level_support = 1;
//default respa level=outermost level is set in init()
dimension = domain->dimension;
@ -109,9 +112,23 @@ int FixFlowGauss::setmask()
int mask = 0;
mask |= POST_FORCE;
mask |= THERMO_ENERGY;
mask |= POST_FORCE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixFlowGauss::init()
{
//if respa level specified by fix_modify, then override default (outermost)
//if specified level too high, set to max level
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0)
ilevel_respa = MIN(respa_level,ilevel_respa);
}
}
/* ----------------------------------------------------------------------
setup is called after the initial evaluation of forces before a run, so we
must remove the total force here too
@ -127,7 +144,13 @@ void FixFlowGauss::setup(int vflag)
if (mTot <= 0.0)
error->all(FLERR,"Invalid group mass in fix flow/gauss");
post_force(vflag);
if (strstr(update->integrate_style,"respa")) {
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
else
post_force(vflag);
}
/* ----------------------------------------------------------------------
@ -148,7 +171,6 @@ void FixFlowGauss::post_force(int vflag)
int ii,jj;
//find the total force on all atoms
//initialize to zero
double f_thisProc[3];
for (ii=0; ii<3; ii++)
@ -201,6 +223,11 @@ void FixFlowGauss::post_force(int vflag)
}
void FixFlowGauss::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ----------------------------------------------------------------------
negative of work done by this fix
This is only computed if requested, either with fix_modify energy yes, or with the energy keyword. Otherwise returns 0.

7
src/USER-MISC/fix_flow_gauss.h
View File

@ -30,10 +30,12 @@ FixStyle(flow/gauss,FixFlowGauss)
public:
FixFlowGauss(class LAMMPS *, int, char **);
int setmask();
void init();
void setup(int);
void post_force(int);
void post_force_respa(int, int, int);
double compute_scalar();
double compute_vector(int n);
void post_force(int);
void setup(int);
protected:
int dimension;
@ -44,6 +46,7 @@ FixStyle(flow/gauss,FixFlowGauss)
double pe_tot; //total added energy
double dt; //timestep
bool workflag; //if calculate work done by fix
int ilevel_respa;
};

42
src/USER-OMP/fix_wall_gran_omp.cpp
View File

@ -25,8 +25,8 @@
using namespace LAMMPS_NS;
using namespace FixConst;
enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER}; // XYZ PLANE need to be 0,1,2
enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY};
enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; // XYZ PLANE need to be 0,1,2
enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,BONDED_HISTORY};
#define BIG 1.0e20
@ -107,6 +107,7 @@ void FixWallGranOMP::post_force(int vflag)
if (mask[i] & groupbit) {
double dx,dy,dz,del1,del2,delxy,delr,rsq;
double rwall = 0.0;
dx = dy = dz = 0.0;
@ -128,13 +129,16 @@ void FixWallGranOMP::post_force(int vflag)
} else if (wallstyle == ZCYLINDER) {
delxy = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1]);
delr = cylradius - delxy;
if (delr > radius[i]) dz = cylradius;
else {
if (delr > radius[i]) {
dz = cylradius;
rwall = 0.0;
} else {
dx = -delr/delxy * x[i][0];
dy = -delr/delxy * x[i][1];
rwall = (delxy < cylradius) ? -2*cylradius : 2*cylradius;
if (wshear && axis != 2) {
vwall[0] = vshear * x[i][1]/delxy;
vwall[1] = -vshear * x[i][0]/delxy;
vwall[0] += vshear * x[i][1]/delxy;
vwall[1] += -vshear * x[i][0]/delxy;
vwall[2] = 0.0;
}
}
@ -143,11 +147,10 @@ void FixWallGranOMP::post_force(int vflag)
rsq = dx*dx + dy*dy + dz*dz;
if (rsq > radius[i]*radius[i]) {
if (pairstyle != HOOKE) {
shear[i][0] = 0.0;
shear[i][1] = 0.0;
shear[i][2] = 0.0;
}
if (history)
for (int j = 0; j < sheardim; j++)
shearone[i][j] = 0.0;
} else {
// meff = effective mass of sphere
@ -156,17 +159,20 @@ void FixWallGranOMP::post_force(int vflag)
double meff = rmass[i];
if (fix_rigid && mass_rigid[i] > 0.0) meff = mass_rigid[i];
// inovke sphere/wall interaction
// invoke sphere/wall interaction
if (pairstyle == HOOKE)
hooke(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],meff);
hooke(rsq,dx,dy,dz,vwall,v[i],f[i],
omega[i],torque[i],radius[i],meff);
else if (pairstyle == HOOKE_HISTORY)
hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],meff,shear[i]);
hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],
omega[i],torque[i],radius[i],meff,shearone[i]);
else if (pairstyle == HERTZ_HISTORY)
hertz_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],meff,shear[i]);
hertz_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i],
omega[i],torque[i],radius[i],meff,shearone[i]);
else if (pairstyle == BONDED_HISTORY)
bonded_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i],
omega[i],torque[i],radius[i],meff,shearone[i]);
}
}
}

5
src/USER-OMP/improper_class2_omp.cpp
View File

@ -696,7 +696,8 @@ void ImproperClass2OMP::angleangle_thr(int nfrom, int nto, ThrData * const thr)
if (EVFLAG)
ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,
fabcd[0],fabcd[2],fabcd[3],delxAB,delyAB,delzAB,
delxBC,delyBC,delzBC,delxBD,delyBD,delzBD,thr);
fabcd[0],fabcd[2],fabcd[3],
delxAB,delyAB,delzAB,delxBC,delyBC,delzBC,
delxBD-delxBC,delyBD-delyBC,delzBD-delzBC,thr);
}
}

18
src/USER-OMP/improper_umbrella_omp.cpp
View File

@ -254,8 +254,24 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr)
f[i4].z += f4[2]*a;
}
if (EVFLAG)
if (EVFLAG) {
// get correct 4-body geometry for virial tally
vb1x = x[i1].x - x[i2].x;
vb1y = x[i1].y - x[i2].y;
vb1z = x[i1].z - x[i2].z;
vb2x = x[i3].x - x[i2].x;
vb2y = x[i3].y - x[i2].y;
vb2z = x[i3].z - x[i2].z;
vb3x = x[i4].x - x[i3].x;
vb3y = x[i4].y - x[i3].y;
vb3z = x[i4].z - x[i3].z;
ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f3,f4,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,thr);
}
}
}