From 619f3e50e225300bfaccdacf36cab4b1e0f0efe1 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 15:21:03 +0000
Subject: [PATCH 01/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7465
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/USER-REAXC/pair_reax_c.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp
index 325e15e3de..a83ee1b5cc 100644
--- a/src/USER-REAXC/pair_reax_c.cpp
+++ b/src/USER-REAXC/pair_reax_c.cpp
@@ -48,7 +48,7 @@
#include "reaxc_traj.h"
#include "reaxc_vector.h"
#include "fix_reaxc_bonds.h"
-#include "fix_reaxc_species.h"
+//#include "fix_reaxc_species.h"
using namespace LAMMPS_NS;
From 9e4fd659e9d6d8ec4e9610543d75f2f5cf7c1306 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 15:22:33 +0000
Subject: [PATCH 02/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7466
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/USER-OMP/pair_born_coul_wolf_omp.cpp | 202 +++++++
src/USER-OMP/pair_born_coul_wolf_omp.h | 48 ++
src/USER-OMP/pair_brownian_omp.cpp | 343 ++++++++++++
src/USER-OMP/pair_brownian_omp.h | 52 ++
src/USER-OMP/pair_brownian_poly_omp.cpp | 330 +++++++++++
src/USER-OMP/pair_brownian_poly_omp.h | 52 ++
src/USER-OMP/pair_coul_wolf_omp.cpp | 183 +++++++
src/USER-OMP/pair_coul_wolf_omp.h | 48 ++
src/USER-OMP/pair_lj_class2_coul_pppm_omp.cpp | 226 ++++++++
src/USER-OMP/pair_lj_class2_coul_pppm_omp.h | 54 ++
src/USER-OMP/pair_lubricate_omp.cpp | 413 ++++++++++++++
src/USER-OMP/pair_lubricate_omp.h | 49 ++
src/USER-OMP/pair_lubricate_poly_omp.cpp | 418 ++++++++++++++
src/USER-OMP/pair_lubricate_poly_omp.h | 49 ++
src/USER-OMP/pppm_cg_omp.cpp | 517 ++++++++++++++++++
src/USER-OMP/pppm_cg_omp.h | 51 ++
16 files changed, 3035 insertions(+)
create mode 100644 src/USER-OMP/pair_born_coul_wolf_omp.cpp
create mode 100644 src/USER-OMP/pair_born_coul_wolf_omp.h
create mode 100644 src/USER-OMP/pair_brownian_omp.cpp
create mode 100644 src/USER-OMP/pair_brownian_omp.h
create mode 100644 src/USER-OMP/pair_brownian_poly_omp.cpp
create mode 100644 src/USER-OMP/pair_brownian_poly_omp.h
create mode 100644 src/USER-OMP/pair_coul_wolf_omp.cpp
create mode 100644 src/USER-OMP/pair_coul_wolf_omp.h
create mode 100644 src/USER-OMP/pair_lj_class2_coul_pppm_omp.cpp
create mode 100644 src/USER-OMP/pair_lj_class2_coul_pppm_omp.h
create mode 100644 src/USER-OMP/pair_lubricate_omp.cpp
create mode 100644 src/USER-OMP/pair_lubricate_omp.h
create mode 100644 src/USER-OMP/pair_lubricate_poly_omp.cpp
create mode 100644 src/USER-OMP/pair_lubricate_poly_omp.h
create mode 100644 src/USER-OMP/pppm_cg_omp.cpp
create mode 100644 src/USER-OMP/pppm_cg_omp.h
diff --git a/src/USER-OMP/pair_born_coul_wolf_omp.cpp b/src/USER-OMP/pair_born_coul_wolf_omp.cpp
new file mode 100644
index 0000000000..fe7e7ae49f
--- /dev/null
+++ b/src/USER-OMP/pair_born_coul_wolf_omp.cpp
@@ -0,0 +1,202 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_born_coul_wolf_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairBornCoulWolfOMP::PairBornCoulWolfOMP(LAMMPS *lmp) :
+ PairBornCoulWolf(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ respa_enable = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBornCoulWolfOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0>(ifrom, ito, thr);
+ }
+
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+/* ---------------------------------------------------------------------- */
+
+template
+void PairBornCoulWolfOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+ double rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj;
+ double prefactor;
+ double r,rexp;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double erfcc,erfcd,v_sh,dvdrr,e_self,e_shift,f_shift,qisq;
+
+ evdwl = ecoul = 0.0;
+
+ const double * const * const x = atom->x;
+ double * const * const f = thr->get_f();
+ const double * const q = atom->q;
+ const int * const type = atom->type;
+ int nlocal = atom->nlocal;
+ double *special_coul = force->special_coul;
+ double *special_lj = force->special_lj;
+ double qqrd2e = force->qqrd2e;
+ double fxtmp,fytmp,fztmp;
+
+ // self and shifted coulombic energy
+
+ e_self = v_sh = 0.0;
+ e_shift = erfc(alf*cut_coul)/cut_coul;
+ f_shift = -(e_shift+ 2.0*alf/MY_PIS * exp(-alf*alf*cut_coul*cut_coul)) /
+ cut_coul;
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = iifrom; ii < iito; ++ii) {
+
+ i = ilist[ii];
+ qtmp = q[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+ fxtmp=fytmp=fztmp=0.0;
+
+ qisq = qtmp*qtmp;
+ e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e;
+ if (EVFLAG) ev_tally_thr(this,i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0,thr);
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ 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]) {
+ r2inv = 1.0/rsq;
+ r = sqrt(rsq);
+
+ if (rsq < cut_coulsq) {
+ prefactor = qqrd2e*qtmp*q[j]/r;
+ erfcc = erfc(alf*r);
+ erfcd = exp(-alf*alf*r*r);
+ v_sh = (erfcc - e_shift*r) * prefactor;
+ dvdrr = (erfcc/rsq + 2.0*alf/MY_PIS * erfcd/r) + f_shift;
+ forcecoul = dvdrr*rsq*prefactor;
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ r6inv = r2inv*r2inv*r2inv;
+ rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]);
+ forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv
+ + born3[itype][jtype]*r2inv*r6inv;
+ } else forceborn = 0.0;
+
+ fpair = (forcecoul + factor_lj*forceborn)*r2inv;
+
+ fxtmp += delx*fpair;
+ fytmp += dely*fpair;
+ fztmp += delz*fpair;
+ if (NEWTON_PAIR || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (EFLAG) {
+ if (rsq < cut_coulsq) {
+ ecoul = v_sh;
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else ecoul = 0.0;
+ if (rsq < cut_ljsq[itype][jtype]) {
+ evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv +
+ d[itype][jtype]*r6inv*r2inv - offset[itype][jtype];
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
+ evdwl,ecoul,fpair,delx,dely,delz,thr);
+ }
+ }
+ f[i][0] += fxtmp;
+ f[i][1] += fytmp;
+ f[i][2] += fztmp;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairBornCoulWolfOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairBornCoulWolf::memory_usage();
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_born_coul_wolf_omp.h b/src/USER-OMP/pair_born_coul_wolf_omp.h
new file mode 100644
index 0000000000..a190ad1a07
--- /dev/null
+++ b/src/USER-OMP/pair_born_coul_wolf_omp.h
@@ -0,0 +1,48 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(born/coul/wolf/omp,PairBornCoulWolfOMP)
+
+#else
+
+#ifndef LMP_PAIR_BORN_COUL_WOLF_OMP_H
+#define LMP_PAIR_BORN_COUL_WOLF_OMP_H
+
+#include "pair_born_coul_wolf.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairBornCoulWolfOMP : public PairBornCoulWolf, public ThrOMP {
+
+ public:
+ PairBornCoulWolfOMP(class LAMMPS *);
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ private:
+ template
+ void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pair_brownian_omp.cpp b/src/USER-OMP/pair_brownian_omp.cpp
new file mode 100644
index 0000000000..42594ae1f8
--- /dev/null
+++ b/src/USER-OMP/pair_brownian_omp.cpp
@@ -0,0 +1,343 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_brownian_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "update.h"
+#include "random_mars.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define EPSILON 1.0e-10
+
+/* ---------------------------------------------------------------------- */
+
+PairBrownianOMP::PairBrownianOMP(LAMMPS *lmp) :
+ PairBrownian(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ respa_enable = 0;
+ random_thr = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairBrownianOMP::~PairBrownianOMP()
+{
+ if (random_thr) {
+ for (int i=1; i < comm->nthreads; ++i)
+ delete random_thr[i];
+
+ delete[] random_thr;
+ random_thr = NULL;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBrownianOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+ if (!random_thr)
+ random_thr = new RanMars*[nthreads];
+
+ // to ensure full compatibility with the serial Brownian style
+ // we use is random number generator instance for thread 0
+ random_thr[0] = random;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ // generate a random number generator instance for
+ // all threads != 0. make sure we use unique seeds.
+ if (random_thr && tid > 0)
+ random_thr[tid] = new RanMars(Pair::lmp, seed + comm->me
+ + comm->nprocs*tid);
+ if (flaglog) {
+ if (evflag) {
+ if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (force->newton_pair) eval<0,1,1>(ifrom, ito, thr);
+ else eval<0,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0>(ifrom, ito, thr);
+ }
+ }
+
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+template
+void PairBrownianOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
+ double rsq,r,h_sep,radi;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ const double * const * const x = atom->x;
+ double * const * const f = thr->get_f();
+ double * const * const torque = thr->get_torque();
+ const double * const radius = atom->radius;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+
+ RanMars &rng = *random_thr[thr->get_tid()];
+
+ double vxmu2f = force->vxmu2f;
+ double randr;
+ double prethermostat;
+ double xl[3],a_sq,a_sh,a_pu,Fbmag;
+ double p1[3],p2[3],p3[3];
+ int overlaps = 0;
+
+ // scale factor for Brownian moments
+
+ prethermostat = sqrt(24.0*force->boltz*t_target/update->dt);
+ prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e);
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = iifrom; ii < iito; ++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];
+
+ // FLD contribution to force and torque due to isotropic terms
+
+ if (flagfld) {
+ f[i][0] += prethermostat*sqrt(R0)*(rng.uniform()-0.5);
+ f[i][1] += prethermostat*sqrt(R0)*(rng.uniform()-0.5);
+ f[i][2] += prethermostat*sqrt(R0)*(rng.uniform()-0.5);
+ if (FLAGLOG) {
+ torque[i][0] += prethermostat*sqrt(RT0)*(rng.uniform()-0.5);
+ torque[i][1] += prethermostat*sqrt(RT0)*(rng.uniform()-0.5);
+ torque[i][2] += prethermostat*sqrt(RT0)*(rng.uniform()-0.5);
+ }
+ }
+
+ 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);
+
+ // scalar resistances a_sq and a_sh
+
+ h_sep = r - 2.0*radi;
+
+ // check for overlaps
+
+ if (h_sep < 0.0) overlaps++;
+
+ // if less than minimum gap, use minimum gap instead
+
+ if (r < cut_inner[itype][jtype])
+ h_sep = cut_inner[itype][jtype] - 2.0*radi;
+
+ // scale h_sep by radi
+
+ h_sep = h_sep/radi;
+
+ // scalar resistances
+
+ if (FLAGLOG) {
+ a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep));
+ a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep));
+ a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep));
+ } else
+ a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep);
+
+ // generate the Pairwise Brownian Force: a_sq
+
+ Fbmag = prethermostat*sqrt(a_sq);
+
+ // generate a random number
+
+ randr = rng.uniform()-0.5;
+
+ // contribution due to Brownian motion
+
+ fx = Fbmag*randr*delx/r;
+ fy = Fbmag*randr*dely/r;
+ fz = Fbmag*randr*delz/r;
+
+ // add terms due to a_sh
+
+ if (FLAGLOG) {
+
+ // generate two orthogonal vectors to the line of centers
+
+ p1[0] = delx/r; p1[1] = dely/r; p1[2] = delz/r;
+ set_3_orthogonal_vectors(p1,p2,p3);
+
+ // magnitude
+
+ Fbmag = prethermostat*sqrt(a_sh);
+
+ // force in each of the two directions
+
+ randr = rng.uniform()-0.5;
+ fx += Fbmag*randr*p2[0];
+ fy += Fbmag*randr*p2[1];
+ fz += Fbmag*randr*p2[2];
+
+ randr = rng.uniform()-0.5;
+ fx += Fbmag*randr*p3[0];
+ fy += Fbmag*randr*p3[1];
+ fz += Fbmag*randr*p3[2];
+ }
+
+ // scale forces to appropriate units
+
+ fx = vxmu2f*fx;
+ fy = vxmu2f*fy;
+ fz = vxmu2f*fz;
+
+ // sum to total force
+
+ f[i][0] -= fx;
+ f[i][1] -= fy;
+ f[i][2] -= fz;
+
+ if (NEWTON_PAIR || j < nlocal) {
+ //randr = rng.uniform()-0.5;
+ //fx = Fbmag*randr*delx/r;
+ //fy = Fbmag*randr*dely/r;
+ //fz = Fbmag*randr*delz/r;
+
+ f[j][0] += fx;
+ f[j][1] += fy;
+ f[j][2] += fz;
+ }
+
+ // torque due to the Brownian Force
+
+ if (FLAGLOG) {
+
+ // location of the point of closest approach on I from its center
+
+ xl[0] = -delx/r*radi;
+ xl[1] = -dely/r*radi;
+ xl[2] = -delz/r*radi;
+
+ // torque = xl_cross_F
+
+ tx = xl[1]*fz - xl[2]*fy;
+ ty = xl[2]*fx - xl[0]*fz;
+ tz = xl[0]*fy - xl[1]*fx;
+
+ // torque is same on both particles
+
+ torque[i][0] -= tx;
+ torque[i][1] -= ty;
+ torque[i][2] -= tz;
+
+ if (NEWTON_PAIR || j < nlocal) {
+ torque[j][0] -= tx;
+ torque[j][1] -= ty;
+ torque[j][2] -= tz;
+ }
+
+ // torque due to a_pu
+
+ Fbmag = prethermostat*sqrt(a_pu);
+
+ // force in each direction
+
+ randr = rng.uniform()-0.5;
+ tx = Fbmag*randr*p2[0];
+ ty = Fbmag*randr*p2[1];
+ tz = Fbmag*randr*p2[2];
+
+ randr = rng.uniform()-0.5;
+ tx += Fbmag*randr*p3[0];
+ ty += Fbmag*randr*p3[1];
+ tz += Fbmag*randr*p3[2];
+
+ // torque has opposite sign on two particles
+
+ torque[i][0] -= tx;
+ torque[i][1] -= ty;
+ torque[i][2] -= tz;
+
+ if (NEWTON_PAIR || j < nlocal) {
+ torque[j][0] += tx;
+ torque[j][1] += ty;
+ torque[j][2] += tz;
+ }
+ }
+
+ if (EVFLAG) ev_tally_xyz(i,j,nlocal,NEWTON_PAIR,
+ 0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairBrownianOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairBrownian::memory_usage();
+ bytes += comm->nthreads * sizeof(RanMars*);
+ bytes += comm->nthreads * sizeof(RanMars);
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_brownian_omp.h b/src/USER-OMP/pair_brownian_omp.h
new file mode 100644
index 0000000000..adf6d3362d
--- /dev/null
+++ b/src/USER-OMP/pair_brownian_omp.h
@@ -0,0 +1,52 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(brownian/omp,PairBrownianOMP)
+
+#else
+
+#ifndef LMP_PAIR_BROWNIAN_OMP_H
+#define LMP_PAIR_BROWNIAN_OMP_H
+
+#include "pair_brownian.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairBrownianOMP : public PairBrownian, public ThrOMP {
+
+ public:
+ PairBrownianOMP(class LAMMPS *);
+ virtual ~PairBrownianOMP();
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ protected:
+ class RanMars **random_thr;
+
+ private:
+ template
+ void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pair_brownian_poly_omp.cpp b/src/USER-OMP/pair_brownian_poly_omp.cpp
new file mode 100644
index 0000000000..9d1dbb22e6
--- /dev/null
+++ b/src/USER-OMP/pair_brownian_poly_omp.cpp
@@ -0,0 +1,330 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_brownian_poly_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "update.h"
+#include "random_mars.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define EPSILON 1.0e-10
+
+/* ---------------------------------------------------------------------- */
+
+PairBrownianPolyOMP::PairBrownianPolyOMP(LAMMPS *lmp) :
+ PairBrownianPoly(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ respa_enable = 0;
+ random_thr = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairBrownianPolyOMP::~PairBrownianPolyOMP()
+{
+ if (random_thr) {
+ for (int i=1; i < comm->nthreads; ++i)
+ delete random_thr[i];
+
+ delete[] random_thr;
+ random_thr = NULL;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBrownianPolyOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+ if (!random_thr)
+ random_thr = new RanMars*[nthreads];
+
+ // to ensure full compatibility with the serial BrownianPoly style
+ // we use is random number generator instance for thread 0
+ random_thr[0] = random;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ // generate a random number generator instance for
+ // all threads != 0. make sure we use unique seeds.
+ if (random_thr && tid > 0)
+ random_thr[tid] = new RanMars(Pair::lmp, seed + comm->me
+ + comm->nprocs*tid);
+ if (flaglog) {
+ if (evflag)
+ eval<1,1>(ifrom, ito, thr);
+ else
+ eval<1,0>(ifrom, ito, thr);
+ } else {
+ if (evflag)
+ eval<0,1>(ifrom, ito, thr);
+ else eval<0,0>(ifrom, ito, thr);
+ }
+
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+template
+void PairBrownianPolyOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
+ double rsq,r,h_sep,beta0,beta1,radi,radj;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ const double * const * const x = atom->x;
+ double * const * const f = thr->get_f();
+ double * const * const torque = thr->get_torque();
+ const double * const radius = atom->radius;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+
+ RanMars &rng = *random_thr[thr->get_tid()];
+
+ double vxmu2f = force->vxmu2f;
+ int overlaps = 0;
+ double randr;
+ double prethermostat;
+ double xl[3],a_sq,a_sh,a_pu,Fbmag;
+ double p1[3],p2[3],p3[3];
+
+ // scale factor for Brownian moments
+
+ prethermostat = sqrt(24.0*force->boltz*t_target/update->dt);
+ prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e);
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = iifrom; ii < iito; ++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];
+
+ // FLD contribution to force and torque due to isotropic terms
+
+ if (flagfld) {
+ f[i][0] += prethermostat*sqrt(R0)*(rng.uniform()-0.5);
+ f[i][1] += prethermostat*sqrt(R0)*(rng.uniform()-0.5);
+ f[i][2] += prethermostat*sqrt(R0)*(rng.uniform()-0.5);
+ if (FLAGLOG) {
+ const double rad3 = radi*radi*radi;
+ torque[i][0] += prethermostat*sqrt(RT0*rad3)*(rng.uniform()-0.5);
+ torque[i][1] += prethermostat*sqrt(RT0*rad3)*(rng.uniform()-0.5);
+ torque[i][2] += prethermostat*sqrt(RT0*rad3)*(rng.uniform()-0.5);
+ }
+ }
+
+ 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];
+ radj = radius[j];
+
+ if (rsq < cutsq[itype][jtype]) {
+ r = sqrt(rsq);
+
+ // scalar resistances a_sq and a_sh
+
+ h_sep = r - radi-radj;
+
+ // check for overlaps
+
+ if (h_sep < 0.0) overlaps++;
+
+ // if less than minimum gap, use minimum gap instead
+
+ if (r < cut_inner[itype][jtype])
+ h_sep = cut_inner[itype][jtype] - radi-radj;
+
+ // scale h_sep by radi
+
+ h_sep = h_sep/radi;
+ beta0 = radj/radi;
+ beta1 = 1.0 + beta0;
+
+ // scalar resistances
+
+ if (FLAGLOG) {
+ a_sq = beta0*beta0/beta1/beta1/h_sep +
+ (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep);
+ a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3) +
+ pow(beta0,4))/21.0/pow(beta1,4)*h_sep*log(1.0/h_sep);
+ a_sq *= 6.0*MY_PI*mu*radi;
+ a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3) *
+ log(1.0/h_sep);
+ a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) +
+ 16.0*pow(beta0,4))/375.0/pow(beta1,4) *
+ h_sep*log(1.0/h_sep);
+ a_sh *= 6.0*MY_PI*mu*radi;
+ a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep);
+ a_pu += (32.0-33.0*beta0+83.0*beta0*beta0+43.0 *
+ pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep);
+ a_pu *= 8.0*MY_PI*mu*pow(radi,3);
+
+ } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep);
+
+ // generate the Pairwise Brownian Force: a_sq
+
+ Fbmag = prethermostat*sqrt(a_sq);
+
+ // generate a random number
+
+ randr = rng.uniform()-0.5;
+
+ // contribution due to Brownian motion
+
+ fx = Fbmag*randr*delx/r;
+ fy = Fbmag*randr*dely/r;
+ fz = Fbmag*randr*delz/r;
+
+ // add terms due to a_sh
+
+ if (FLAGLOG) {
+
+ // generate two orthogonal vectors to the line of centers
+
+ p1[0] = delx/r; p1[1] = dely/r; p1[2] = delz/r;
+ set_3_orthogonal_vectors(p1,p2,p3);
+
+ // magnitude
+
+ Fbmag = prethermostat*sqrt(a_sh);
+
+ // force in each of the two directions
+
+ randr = rng.uniform()-0.5;
+ fx += Fbmag*randr*p2[0];
+ fy += Fbmag*randr*p2[1];
+ fz += Fbmag*randr*p2[2];
+
+ randr = rng.uniform()-0.5;
+ fx += Fbmag*randr*p3[0];
+ fy += Fbmag*randr*p3[1];
+ fz += Fbmag*randr*p3[2];
+ }
+
+ // scale forces to appropriate units
+
+ fx = vxmu2f*fx;
+ fy = vxmu2f*fy;
+ fz = vxmu2f*fz;
+
+ // sum to total force
+
+ f[i][0] -= fx;
+ f[i][1] -= fy;
+ f[i][2] -= fz;
+
+ // torque due to the Brownian Force
+
+ if (FLAGLOG) {
+
+ // location of the point of closest approach on I from its center
+
+ xl[0] = -delx/r*radi;
+ xl[1] = -dely/r*radi;
+ xl[2] = -delz/r*radi;
+
+ // torque = xl_cross_F
+
+ tx = xl[1]*fz - xl[2]*fy;
+ ty = xl[2]*fx - xl[0]*fz;
+ tz = xl[0]*fy - xl[1]*fx;
+
+ // torque is same on both particles
+
+ torque[i][0] -= tx;
+ torque[i][1] -= ty;
+ torque[i][2] -= tz;
+
+ // torque due to a_pu
+
+ Fbmag = prethermostat*sqrt(a_pu);
+
+ // force in each direction
+
+ randr = rng.uniform()-0.5;
+ tx = Fbmag*randr*p2[0];
+ ty = Fbmag*randr*p2[1];
+ tz = Fbmag*randr*p2[2];
+
+ randr = rng.uniform()-0.5;
+ tx += Fbmag*randr*p3[0];
+ ty += Fbmag*randr*p3[1];
+ tz += Fbmag*randr*p3[2];
+
+ // torque has opposite sign on two particles
+
+ torque[i][0] -= tx;
+ torque[i][1] -= ty;
+ torque[i][2] -= tz;
+
+ }
+
+ if (EVFLAG) ev_tally_xyz(i,j,nlocal,/* newton_pair */ 0,
+ 0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairBrownianPolyOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairBrownianPoly::memory_usage();
+ bytes += comm->nthreads * sizeof(RanMars*);
+ bytes += comm->nthreads * sizeof(RanMars);
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_brownian_poly_omp.h b/src/USER-OMP/pair_brownian_poly_omp.h
new file mode 100644
index 0000000000..23443e51e0
--- /dev/null
+++ b/src/USER-OMP/pair_brownian_poly_omp.h
@@ -0,0 +1,52 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(brownian/poly/omp,PairBrownianPolyOMP)
+
+#else
+
+#ifndef LMP_PAIR_BROWNIAN_POLY_OMP_H
+#define LMP_PAIR_BROWNIAN_POLY_OMP_H
+
+#include "pair_brownian_poly.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairBrownianPolyOMP : public PairBrownianPoly, public ThrOMP {
+
+ public:
+ PairBrownianPolyOMP(class LAMMPS *);
+ virtual ~PairBrownianPolyOMP();
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ protected:
+ class RanMars **random_thr;
+
+ private:
+ template
+ void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pair_coul_wolf_omp.cpp b/src/USER-OMP/pair_coul_wolf_omp.cpp
new file mode 100644
index 0000000000..6bf9a616cf
--- /dev/null
+++ b/src/USER-OMP/pair_coul_wolf_omp.cpp
@@ -0,0 +1,183 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_coul_wolf_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairCoulWolfOMP::PairCoulWolfOMP(LAMMPS *lmp) :
+ PairCoulWolf(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ respa_enable = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairCoulWolfOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0>(ifrom, ito, thr);
+ }
+
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+/* ---------------------------------------------------------------------- */
+
+template
+void PairCoulWolfOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
+ double rsq,forcecoul,factor_coul;
+ double prefactor;
+ double r,rexp;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double erfcc,erfcd,v_sh,dvdrr,e_self,e_shift,f_shift,qisq;
+
+ ecoul = 0.0;
+
+ const double * const * const x = atom->x;
+ double * const * const f = thr->get_f();
+ const double * const q = atom->q;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+ const double * const special_coul = force->special_coul;
+ const double qqrd2e = force->qqrd2e;
+ double fxtmp,fytmp,fztmp;
+
+ // self and shifted coulombic energy
+
+ e_self = v_sh = 0.0;
+ e_shift = erfc(alf*cut_coul)/cut_coul;
+ f_shift = -(e_shift+ 2.0*alf/MY_PIS * exp(-alf*alf*cut_coul*cut_coul)) /
+ cut_coul;
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = iifrom; ii < iito; ++ii) {
+
+ i = ilist[ii];
+ qtmp = q[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+ fxtmp=fytmp=fztmp=0.0;
+
+ qisq = qtmp*qtmp;
+ e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e;
+ if (EVFLAG) ev_tally_thr(this,i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0,thr);
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_coul = special_coul[sbmask(j)];
+ 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 < cut_coulsq) {
+ r = sqrt(rsq);
+ prefactor = qqrd2e*qtmp*q[j]/r;
+ erfcc = erfc(alf*r);
+ erfcd = exp(-alf*alf*r*r);
+ v_sh = (erfcc - e_shift*r) * prefactor;
+ dvdrr = (erfcc/rsq + 2.0*alf/MY_PIS * erfcd/r) + f_shift;
+ forcecoul = dvdrr*rsq*prefactor;
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ fpair = forcecoul / rsq;
+
+ fxtmp += delx*fpair;
+ fytmp += dely*fpair;
+ fztmp += delz*fpair;
+ if (NEWTON_PAIR || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (EFLAG) {
+ if (rsq < cut_coulsq) {
+ ecoul = v_sh;
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else ecoul = 0.0;
+ }
+
+ if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
+ 0.0,ecoul,fpair,delx,dely,delz,thr);
+ }
+ }
+ f[i][0] += fxtmp;
+ f[i][1] += fytmp;
+ f[i][2] += fztmp;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairCoulWolfOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairCoulWolf::memory_usage();
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_coul_wolf_omp.h b/src/USER-OMP/pair_coul_wolf_omp.h
new file mode 100644
index 0000000000..0b474601f2
--- /dev/null
+++ b/src/USER-OMP/pair_coul_wolf_omp.h
@@ -0,0 +1,48 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(coul/wolf/omp,PairCoulWolfOMP)
+
+#else
+
+#ifndef LMP_PAIR_COUL_WOLF_OMP_H
+#define LMP_PAIR_COUL_WOLF_OMP_H
+
+#include "pair_coul_wolf.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairCoulWolfOMP : public PairCoulWolf, public ThrOMP {
+
+ public:
+ PairCoulWolfOMP(class LAMMPS *);
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ private:
+ template
+ void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pair_lj_class2_coul_pppm_omp.cpp b/src/USER-OMP/pair_lj_class2_coul_pppm_omp.cpp
new file mode 100644
index 0000000000..56011e17b1
--- /dev/null
+++ b/src/USER-OMP/pair_lj_class2_coul_pppm_omp.cpp
@@ -0,0 +1,226 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_lj_class2_coul_pppm_omp.h"
+#include "pppm_proxy.h"
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "update.h"
+
+#include
+
+using namespace LAMMPS_NS;
+
+#define EWALD_F 1.12837917
+#define EWALD_P 0.3275911
+#define A1 0.254829592
+#define A2 -0.284496736
+#define A3 1.421413741
+#define A4 -1.453152027
+#define A5 1.061405429
+
+/* ---------------------------------------------------------------------- */
+
+PairLJClass2CoulPPPMOMP::PairLJClass2CoulPPPMOMP(LAMMPS *lmp) :
+ PairLJClass2CoulLong(lmp), ThrOMP(lmp, THR_PAIR|THR_PROXY)
+{
+ respa_enable = 0;
+ nproxy=1;
+
+ kspace = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJClass2CoulPPPMOMP::init_style()
+{
+ if (comm->nthreads < 2)
+ error->all(FLERR,"need at least two threads per MPI task for this pair style");
+
+ if (strcmp(force->kspace_style,"pppm/proxy") != 0)
+ error->all(FLERR,"kspace style pppm/proxy is required with this pair style");
+
+ kspace = static_cast(force->kspace);
+
+ PairLJClass2CoulLong::init_style();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJClass2CoulPPPMOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads, nproxy);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ // thread id 0 runs pppm, the rest the pair style
+ if (tid < nproxy) {
+ kspace->compute_proxy(eflag,vflag);
+ } else {
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0>(ifrom, ito, thr);
+ }
+ }
+
+ sync_threads();
+ reduce_thr(this, eflag, vflag, thr, nproxy);
+ } // end of omp parallel region
+}
+
+/* ---------------------------------------------------------------------- */
+
+template
+void PairLJClass2CoulPPPMOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+ double r,rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
+ double grij,expm2,prefactor,t,erfc;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = ecoul = 0.0;
+
+ const double * const * const x = atom->x;
+ double * const * const f = thr->get_f();
+ const double * const q = atom->q;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+ const double * const special_coul = force->special_coul;
+ const double * const special_lj = force->special_lj;
+ const double qqrd2e = force->qqrd2e;
+ double fxtmp,fytmp,fztmp;
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = iifrom; ii < iito; ++ii) {
+
+ i = ilist[ii];
+ qtmp = q[i];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+ fxtmp=fytmp=fztmp=0.0;
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ factor_lj = special_lj[sbmask(j)];
+ factor_coul = special_coul[sbmask(j)];
+ 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]) {
+ r2inv = 1.0/rsq;
+
+ if (rsq < cut_coulsq) {
+ r = sqrt(rsq);
+ grij = g_ewald * r;
+ expm2 = exp(-grij*grij);
+ t = 1.0 / (1.0 + EWALD_P*grij);
+ erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+ prefactor = qqrd2e * qtmp*q[j]/r;
+ forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+ if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+ } else forcecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ rinv = sqrt(r2inv);
+ r3inv = r2inv*rinv;
+ r6inv = r3inv*r3inv;
+ forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
+ forcelj *= factor_lj;
+ } else forcelj = 0.0;
+
+ fpair = (forcecoul + forcelj) * r2inv;
+
+ fxtmp += delx*fpair;
+ fytmp += dely*fpair;
+ fztmp += delz*fpair;
+ if (NEWTON_PAIR || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (EFLAG) {
+ if (rsq < cut_coulsq) {
+ ecoul = prefactor*erfc;
+ if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+ } else ecoul = 0.0;
+
+ if (rsq < cut_ljsq[itype][jtype]) {
+ evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
+ offset[itype][jtype];
+ evdwl *= factor_lj;
+ } else evdwl = 0.0;
+ }
+
+ if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
+ evdwl,ecoul,fpair,delx,dely,delz,thr);
+ }
+ }
+ f[i][0] += fxtmp;
+ f[i][1] += fytmp;
+ f[i][2] += fztmp;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJClass2CoulPPPMOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairLJClass2CoulLong::memory_usage();
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_lj_class2_coul_pppm_omp.h b/src/USER-OMP/pair_lj_class2_coul_pppm_omp.h
new file mode 100644
index 0000000000..3f77bf4f9a
--- /dev/null
+++ b/src/USER-OMP/pair_lj_class2_coul_pppm_omp.h
@@ -0,0 +1,54 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(lj/class2/coul/pppm/omp,PairLJClass2CoulPPPMOMP)
+
+#else
+
+#ifndef LMP_PAIR_LJ_CLASS2_COUL_PPPM_OMP_H
+#define LMP_PAIR_LJ_CLASS2_COUL_PPPM_OMP_H
+
+#include "pair_lj_class2_coul_long.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairLJClass2CoulPPPMOMP : public PairLJClass2CoulLong, public ThrOMP {
+
+ public:
+ PairLJClass2CoulPPPMOMP(class LAMMPS *);
+ virtual ~PairLJClass2CoulPPPMOMP() {};
+
+ virtual void init_style();
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ private:
+ template
+ void eval(int ifrom, int ito, ThrData * const thr);
+
+ class PPPMProxy *kspace;
+ int nproxy;
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pair_lubricate_omp.cpp b/src/USER-OMP/pair_lubricate_omp.cpp
new file mode 100644
index 0000000000..89762acdc5
--- /dev/null
+++ b/src/USER-OMP/pair_lubricate_omp.cpp
@@ -0,0 +1,413 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_lubricate_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "domain.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "update.h"
+#include "random_mars.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairLubricateOMP::PairLubricateOMP(LAMMPS *lmp) :
+ PairLubricate(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ respa_enable = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairLubricateOMP::~PairLubricateOMP()
+{}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLubricateOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ if (flaglog) {
+ if (evflag) {
+ if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (evflag) {
+ if (force->newton_pair) eval<0,1,1>(ifrom, ito, thr);
+ else eval<0,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0>(ifrom, ito, thr);
+ }
+ }
+
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+template
+void PairLubricateOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
+ double rsq,r,h_sep,radi;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
+ double vt1,vt2,vt3,wt1,wt2,wt3,wdotn;
+ double vRS0;
+ double vi[3],vj[3],wi[3],wj[3],xl[3];
+ double a_sq,a_sh,a_pu;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double lamda[3],vstream[3];
+
+ double vxmu2f = force->vxmu2f;
+
+ double * const * const x = atom->x;
+ double * const * const v = atom->v;
+ double * const * const f = thr->get_f();
+ double * const * const omega = atom->omega;
+ double * const * const torque = thr->get_torque();
+ const double * const radius = atom->radius;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+
+ int overlaps = 0;
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // subtract streaming component of velocity, omega, angmom
+ // assume fluid streaming velocity = box deformation rate
+ // vstream = (ux,uy,uz)
+ // ux = h_rate[0]*x + h_rate[5]*y + h_rate[4]*z
+ // uy = h_rate[1]*y + h_rate[3]*z
+ // uz = h_rate[2]*z
+ // omega_new = omega - curl(vstream)/2
+ // angmom_new = angmom - I*curl(vstream)/2
+ // Ef = (grad(vstream) + (grad(vstream))^T) / 2
+
+ if (shearing) {
+ double *h_rate = domain->h_rate;
+ double *h_ratelo = domain->h_ratelo;
+
+ for (ii = iifrom; ii < iito; ii++) {
+ i = ilist[ii];
+ itype = type[i];
+ radi = radius[i];
+
+ domain->x2lamda(x[i],lamda);
+ vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
+ h_rate[4]*lamda[2] + h_ratelo[0];
+ vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
+ vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
+ v[i][0] -= vstream[0];
+ v[i][1] -= vstream[1];
+ v[i][2] -= vstream[2];
+
+ omega[i][0] += 0.5*h_rate[3];
+ omega[i][1] -= 0.5*h_rate[4];
+ omega[i][2] += 0.5*h_rate[5];
+ }
+
+ // set Ef from h_rate in strain units
+
+ Ef[0][0] = h_rate[0]/domain->xprd;
+ Ef[1][1] = h_rate[1]/domain->yprd;
+ Ef[2][2] = h_rate[2]/domain->zprd;
+ Ef[0][1] = Ef[1][0] = 0.5 * h_rate[5]/domain->yprd;
+ Ef[0][2] = Ef[2][0] = 0.5 * h_rate[4]/domain->zprd;
+ Ef[1][2] = Ef[2][1] = 0.5 * h_rate[3]/domain->zprd;
+
+ // copy updated velocity/omega/angmom to the ghost particles
+ // no need to do this if not shearing since comm->ghost_velocity is set
+
+ sync_threads();
+ comm->forward_comm_pair(this);
+ }
+
+ // loop over neighbors of my atoms
+
+ for (ii = iifrom; ii < iito; ++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];
+
+ // angular velocity
+
+ wi[0] = omega[i][0];
+ wi[1] = omega[i][1];
+ wi[2] = omega[i][2];
+
+ // FLD contribution to force and torque due to isotropic terms
+ // FLD contribution to stress from isotropic RS0
+
+ if (flagfld) {
+ f[i][0] -= vxmu2f*R0*v[i][0];
+ f[i][1] -= vxmu2f*R0*v[i][1];
+ f[i][2] -= vxmu2f*R0*v[i][2];
+ torque[i][0] -= vxmu2f*RT0*wi[0];
+ torque[i][1] -= vxmu2f*RT0*wi[1];
+ torque[i][2] -= vxmu2f*RT0*wi[2];
+
+ if (shearing && vflag_either) {
+ vRS0 = -vxmu2f * RS0;
+ v_tally_tensor(i,i,nlocal,NEWTON_PAIR,
+ vRS0*Ef[0][0],vRS0*Ef[1][1],vRS0*Ef[2][2],
+ vRS0*Ef[0][1],vRS0*Ef[0][2],vRS0*Ef[1][2]);
+ }
+ }
+
+ 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);
+
+ // angular momentum = I*omega = 2/5 * M*R^2 * omega
+
+ wj[0] = omega[j][0];
+ wj[1] = omega[j][1];
+ wj[2] = omega[j][2];
+
+ // xl = point of closest approach on particle i from its center
+
+ xl[0] = -delx/r*radi;
+ xl[1] = -dely/r*radi;
+ xl[2] = -delz/r*radi;
+
+ // velocity at the point of closest approach on both particles
+ // v = v + omega_cross_xl - Ef.xl
+
+ // particle i
+
+ vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1])
+ - (Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]);
+
+ vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2])
+ - (Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]);
+
+ vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0])
+ - (Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]);
+
+ // particle j
+
+ vj[0] = v[j][0] - (wj[1]*xl[2] - wj[2]*xl[1])
+ + (Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]);
+
+ vj[1] = v[j][1] - (wj[2]*xl[0] - wj[0]*xl[2])
+ + (Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]);
+
+ vj[2] = v[j][2] - (wj[0]*xl[1] - wj[1]*xl[0])
+ + (Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]);
+
+ // scalar resistances XA and YA
+
+ h_sep = r - 2.0*radi;
+
+ // check for overlaps
+
+ if (h_sep < 0.0) overlaps++;
+
+ // if less than the minimum gap use the minimum gap instead
+
+ if (r < cut_inner[itype][jtype])
+ h_sep = cut_inner[itype][jtype] - 2.0*radi;
+
+ // scale h_sep by radi
+
+ h_sep = h_sep/radi;
+
+ // scalar resistances
+
+ if (FLAGLOG) {
+ a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep));
+ a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep));
+ a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep));
+ } else
+ a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep);
+
+ // relative velocity at the point of closest approach
+ // includes fluid velocity
+
+ vr1 = vi[0] - vj[0];
+ vr2 = vi[1] - vj[1];
+ vr3 = vi[2] - vj[2];
+
+ // normal component (vr.n)n
+
+ vnnr = (vr1*delx + vr2*dely + vr3*delz)/r;
+ vn1 = vnnr*delx/r;
+ vn2 = vnnr*dely/r;
+ vn3 = vnnr*delz/r;
+
+ // tangential component vr - (vr.n)n
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // force due to squeeze type motion
+
+ fx = a_sq*vn1;
+ fy = a_sq*vn2;
+ fz = a_sq*vn3;
+
+ // force due to all shear kind of motions
+
+ if (FLAGLOG) {
+ fx = fx + a_sh*vt1;
+ fy = fy + a_sh*vt2;
+ fz = fz + a_sh*vt3;
+ }
+
+ // scale forces for appropriate units
+
+ fx *= vxmu2f;
+ fy *= vxmu2f;
+ fz *= vxmu2f;
+
+ // add to total force
+
+ 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;
+ }
+
+ // torque due to this force
+
+ if (FLAGLOG) {
+ tx = xl[1]*fz - xl[2]*fy;
+ ty = xl[2]*fx - xl[0]*fz;
+ tz = xl[0]*fy - xl[1]*fx;
+
+ torque[i][0] -= vxmu2f*tx;
+ torque[i][1] -= vxmu2f*ty;
+ torque[i][2] -= vxmu2f*tz;
+
+ if (NEWTON_PAIR || j < nlocal) {
+ torque[j][0] -= vxmu2f*tx;
+ torque[j][1] -= vxmu2f*ty;
+ torque[j][2] -= vxmu2f*tz;
+ }
+
+ // torque due to a_pu
+
+ wdotn = ((wi[0]-wj[0])*delx + (wi[1]-wj[1])*dely +
+ (wi[2]-wj[2])*delz)/r;
+ wt1 = (wi[0]-wj[0]) - wdotn*delx/r;
+ wt2 = (wi[1]-wj[1]) - wdotn*dely/r;
+ wt3 = (wi[2]-wj[2]) - wdotn*delz/r;
+
+ tx = a_pu*wt1;
+ ty = a_pu*wt2;
+ tz = a_pu*wt3;
+
+ torque[i][0] -= vxmu2f*tx;
+ torque[i][1] -= vxmu2f*ty;
+ torque[i][2] -= vxmu2f*tz;
+
+ if (NEWTON_PAIR || j < nlocal) {
+ 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);
+ }
+ }
+ }
+
+ // restore streaming component of velocity, omega, angmom
+
+ if (shearing) {
+ double *h_rate = domain->h_rate;
+ double *h_ratelo = domain->h_ratelo;
+
+ for (ii = iifrom; ii < iito; ii++) {
+ i = ilist[ii];
+ itype = type[i];
+ radi = radius[i];
+
+ domain->x2lamda(x[i],lamda);
+ vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
+ h_rate[4]*lamda[2] + h_ratelo[0];
+ vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
+ vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
+ v[i][0] += vstream[0];
+ v[i][1] += vstream[1];
+ v[i][2] += vstream[2];
+
+ omega[i][0] -= 0.5*h_rate[3];
+ omega[i][1] += 0.5*h_rate[4];
+ omega[i][2] -= 0.5*h_rate[5];
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLubricateOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairLubricate::memory_usage();
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_lubricate_omp.h b/src/USER-OMP/pair_lubricate_omp.h
new file mode 100644
index 0000000000..045421e566
--- /dev/null
+++ b/src/USER-OMP/pair_lubricate_omp.h
@@ -0,0 +1,49 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(lubricate/omp,PairLubricateOMP)
+
+#else
+
+#ifndef LMP_PAIR_LUBRICATE_OMP_H
+#define LMP_PAIR_LUBRICATE_OMP_H
+
+#include "pair_lubricate.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairLubricateOMP : public PairLubricate, public ThrOMP {
+
+ public:
+ PairLubricateOMP(class LAMMPS *);
+ virtual ~PairLubricateOMP();
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ private:
+ template
+ void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pair_lubricate_poly_omp.cpp b/src/USER-OMP/pair_lubricate_poly_omp.cpp
new file mode 100644
index 0000000000..888b01f679
--- /dev/null
+++ b/src/USER-OMP/pair_lubricate_poly_omp.cpp
@@ -0,0 +1,418 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_lubricate_poly_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "domain.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "update.h"
+#include "random_mars.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairLubricatePolyOMP::PairLubricatePolyOMP(LAMMPS *lmp) :
+ PairLubricatePoly(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ respa_enable = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairLubricatePolyOMP::~PairLubricatePolyOMP()
+{}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLubricatePolyOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ if (flaglog) {
+ if (shearing) {
+ if (evflag)
+ eval<1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0>(ifrom, ito, thr);
+ } else {
+ if (evflag)
+ eval<1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (shearing) {
+ if (evflag)
+ eval<0,1,1>(ifrom, ito, thr);
+ else eval<0,1,0>(ifrom, ito, thr);
+ } else {
+ if (evflag)
+ eval<0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0>(ifrom, ito, thr);
+ }
+ }
+
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+template
+void PairLubricatePolyOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,ii,jj,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
+ double rsq,r,h_sep,beta0,beta1,radi,radj;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
+ double vt1,vt2,vt3,wt1,wt2,wt3,wdotn;
+ double vRS0;
+ double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3];
+ double a_sq,a_sh,a_pu;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double lamda[3],vstream[3];
+
+ double vxmu2f = force->vxmu2f;
+
+ double * const * const x = atom->x;
+ double * const * const v = atom->v;
+ double * const * const f = thr->get_f();
+ double * const * const omega = atom->omega;
+ double * const * const torque = thr->get_torque();
+ const double * const radius = atom->radius;
+ const int * const type = atom->type;
+ const int nlocal = atom->nlocal;
+
+ int overlaps = 0;
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // subtract streaming component of velocity, omega, angmom
+ // assume fluid streaming velocity = box deformation rate
+ // vstream = (ux,uy,uz)
+ // ux = h_rate[0]*x + h_rate[5]*y + h_rate[4]*z
+ // uy = h_rate[1]*y + h_rate[3]*z
+ // uz = h_rate[2]*z
+ // omega_new = omega - curl(vstream)/2
+ // angmom_new = angmom - I*curl(vstream)/2
+ // Ef = (grad(vstream) + (grad(vstream))^T) / 2
+
+ if (shearing) {
+ double *h_rate = domain->h_rate;
+ double *h_ratelo = domain->h_ratelo;
+
+ for (ii = iifrom; ii < iito; ii++) {
+ i = ilist[ii];
+ itype = type[i];
+ radi = radius[i];
+ domain->x2lamda(x[i],lamda);
+ vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
+ h_rate[4]*lamda[2] + h_ratelo[0];
+ vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
+ vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
+ v[i][0] -= vstream[0];
+ v[i][1] -= vstream[1];
+ v[i][2] -= vstream[2];
+
+ omega[i][0] += 0.5*h_rate[3];
+ omega[i][1] -= 0.5*h_rate[4];
+ omega[i][2] += 0.5*h_rate[5];
+ }
+
+ // set Ef from h_rate in strain units
+
+ Ef[0][0] = h_rate[0]/domain->xprd;
+ Ef[1][1] = h_rate[1]/domain->yprd;
+ Ef[2][2] = h_rate[2]/domain->zprd;
+ Ef[0][1] = Ef[1][0] = 0.5 * h_rate[5]/domain->yprd;
+ Ef[0][2] = Ef[2][0] = 0.5 * h_rate[4]/domain->zprd;
+ Ef[1][2] = Ef[2][1] = 0.5 * h_rate[3]/domain->zprd;
+
+ // copy updated omega to the ghost particles
+ // no need to do this if not shearing since comm->ghost_velocity is set
+
+ sync_threads();
+ comm->forward_comm_pair(this);
+ }
+
+ // loop over neighbors of my atoms
+
+ for (ii = iifrom; ii < iito; ++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];
+
+ // angular velocity
+
+ wi[0] = omega[i][0];
+ wi[1] = omega[i][1];
+ wi[2] = omega[i][2];
+
+ // FLD contribution to force and torque due to isotropic terms
+ // FLD contribution to stress from isotropic RS0
+
+ if (flagfld) {
+ f[i][0] -= vxmu2f*R0*radi*v[i][0];
+ f[i][1] -= vxmu2f*R0*radi*v[i][1];
+ f[i][2] -= vxmu2f*R0*radi*v[i][2];
+ const double rad3 = radi*radi*radi;
+ torque[i][0] -= vxmu2f*RT0*rad3*wi[0];
+ torque[i][1] -= vxmu2f*RT0*rad3*wi[1];
+ torque[i][2] -= vxmu2f*RT0*rad3*wi[2];
+
+ if (SHEARING && vflag_either) {
+ vRS0 = -vxmu2f * RS0*rad3;
+ v_tally_tensor(i,i,nlocal,/* newton_pair */ 0,
+ vRS0*Ef[0][0],vRS0*Ef[1][1],vRS0*Ef[2][2],
+ vRS0*Ef[0][1],vRS0*Ef[0][2],vRS0*Ef[1][2]);
+ }
+ }
+
+ 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];
+ radj = atom->radius[j];
+
+ if (rsq < cutsq[itype][jtype]) {
+ r = sqrt(rsq);
+
+ // angular momentum = I*omega = 2/5 * M*R^2 * omega
+
+ wj[0] = omega[j][0];
+ wj[1] = omega[j][1];
+ wj[2] = omega[j][2];
+
+ // xl = point of closest approach on particle i from its center
+
+ xl[0] = -delx/r*radi;
+ xl[1] = -dely/r*radi;
+ xl[2] = -delz/r*radi;
+ jl[0] = -delx/r*radj;
+ jl[1] = -dely/r*radj;
+ jl[2] = -delz/r*radj;
+
+ // velocity at the point of closest approach on both particles
+ // v = v + omega_cross_xl - Ef.xl
+
+ // particle i
+
+ vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1])
+ - (Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]);
+
+ vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2])
+ - (Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]);
+
+ vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0])
+ - (Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]);
+
+ // particle j
+
+ vj[0] = v[j][0] - (wj[1]*jl[2] - wj[2]*jl[1])
+ + (Ef[0][0]*jl[0] + Ef[0][1]*jl[1] + Ef[0][2]*jl[2]);
+
+ vj[1] = v[j][1] - (wj[2]*jl[0] - wj[0]*jl[2])
+ + (Ef[1][0]*jl[0] + Ef[1][1]*jl[1] + Ef[1][2]*jl[2]);
+
+ vj[2] = v[j][2] - (wj[0]*jl[1] - wj[1]*jl[0])
+ + (Ef[2][0]*jl[0] + Ef[2][1]*jl[1] + Ef[2][2]*jl[2]);
+
+ // scalar resistances XA and YA
+
+ h_sep = r - radi-radj;
+
+ // check for overlaps
+
+ if (h_sep < 0.0) overlaps++;
+
+ // if less than the minimum gap use the minimum gap instead
+
+ if (r < cut_inner[itype][jtype])
+ h_sep = cut_inner[itype][jtype] - radi-radj;
+
+ // scale h_sep by radi
+
+ h_sep = h_sep/radi;
+ beta0 = radj/radi;
+ beta1 = 1.0 + beta0;
+
+ // scalar resistances
+
+ if (FLAGLOG) {
+ a_sq = beta0*beta0/beta1/beta1/h_sep +
+ (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep);
+ a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0 *
+ pow(beta0,3)+pow(beta0,4))/21.0/pow(beta1,4) *
+ h_sep*log(1.0/h_sep);
+ a_sq *= 6.0*MY_PI*mu*radi;
+ a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3) *
+ log(1.0/h_sep);
+ a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) +
+ 16.0*pow(beta0,4))/375.0/pow(beta1,4) *
+ h_sep*log(1.0/h_sep);
+ a_sh *= 6.0*MY_PI*mu*radi;
+ a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep);
+ a_pu += (32.0-33.0*beta0+83.0*beta0*beta0+43.0 *
+ pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep);
+ a_pu *= 8.0*MY_PI*mu*pow(radi,3);
+ } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep);
+
+ // relative velocity at the point of closest approach
+ // includes fluid velocity
+
+ vr1 = vi[0] - vj[0];
+ vr2 = vi[1] - vj[1];
+ vr3 = vi[2] - vj[2];
+
+ // normal component (vr.n)n
+
+ vnnr = (vr1*delx + vr2*dely + vr3*delz)/r;
+ vn1 = vnnr*delx/r;
+ vn2 = vnnr*dely/r;
+ vn3 = vnnr*delz/r;
+
+ // tangential component vr - (vr.n)n
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // force due to squeeze type motion
+
+ fx = a_sq*vn1;
+ fy = a_sq*vn2;
+ fz = a_sq*vn3;
+
+ // force due to all shear kind of motions
+
+ if (FLAGLOG) {
+ fx = fx + a_sh*vt1;
+ fy = fy + a_sh*vt2;
+ fz = fz + a_sh*vt3;
+ }
+
+ // scale forces for appropriate units
+
+ fx *= vxmu2f;
+ fy *= vxmu2f;
+ fz *= vxmu2f;
+
+ // add to total force
+
+ f[i][0] -= fx;
+ f[i][1] -= fy;
+ f[i][2] -= fz;
+
+ // torque due to this force
+
+ if (FLAGLOG) {
+ tx = xl[1]*fz - xl[2]*fy;
+ ty = xl[2]*fx - xl[0]*fz;
+ tz = xl[0]*fy - xl[1]*fx;
+
+ torque[i][0] -= vxmu2f*tx;
+ torque[i][1] -= vxmu2f*ty;
+ torque[i][2] -= vxmu2f*tz;
+
+ // torque due to a_pu
+
+ wdotn = ((wi[0]-wj[0])*delx + (wi[1]-wj[1])*dely +
+ (wi[2]-wj[2])*delz)/r;
+ wt1 = (wi[0]-wj[0]) - wdotn*delx/r;
+ wt2 = (wi[1]-wj[1]) - wdotn*dely/r;
+ wt3 = (wi[2]-wj[2]) - wdotn*delz/r;
+
+ tx = a_pu*wt1;
+ ty = a_pu*wt2;
+ tz = a_pu*wt3;
+
+ torque[i][0] -= vxmu2f*tx;
+ torque[i][1] -= vxmu2f*ty;
+ torque[i][2] -= vxmu2f*tz;
+
+ }
+
+ if (EVFLAG) ev_tally_xyz(i,nlocal,nlocal, /* newton_pair */ 0,
+ 0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
+ }
+ }
+ }
+
+ // restore streaming component of velocity, omega, angmom
+
+ if (SHEARING) {
+ double *h_rate = domain->h_rate;
+ double *h_ratelo = domain->h_ratelo;
+
+ for (ii = iifrom; ii < iito; ii++) {
+ i = ilist[ii];
+ itype = type[i];
+ radi = radius[i];
+
+ domain->x2lamda(x[i],lamda);
+ vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
+ h_rate[4]*lamda[2] + h_ratelo[0];
+ vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
+ vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2];
+ v[i][0] += vstream[0];
+ v[i][1] += vstream[1];
+ v[i][2] += vstream[2];
+
+ omega[i][0] -= 0.5*h_rate[3];
+ omega[i][1] += 0.5*h_rate[4];
+ omega[i][2] -= 0.5*h_rate[5];
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLubricatePolyOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairLubricatePoly::memory_usage();
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_lubricate_poly_omp.h b/src/USER-OMP/pair_lubricate_poly_omp.h
new file mode 100644
index 0000000000..8e44d94093
--- /dev/null
+++ b/src/USER-OMP/pair_lubricate_poly_omp.h
@@ -0,0 +1,49 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(lubricate/poly/omp,PairLubricateOMP)
+
+#else
+
+#ifndef LMP_PAIR_LUBRICATE_POLY_OMP_H
+#define LMP_PAIR_LUBRICATE_POLY_OMP_H
+
+#include "pair_lubricate_poly.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairLubricatePolyOMP : public PairLubricatePoly, public ThrOMP {
+
+ public:
+ PairLubricatePolyOMP(class LAMMPS *);
+ virtual ~PairLubricatePolyOMP();
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ private:
+ template
+ void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pppm_cg_omp.cpp b/src/USER-OMP/pppm_cg_omp.cpp
new file mode 100644
index 0000000000..d89c7d4358
--- /dev/null
+++ b/src/USER-OMP/pppm_cg_omp.cpp
@@ -0,0 +1,517 @@
+/* ----------------------------------------------------------------------
+ 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "pppm_cg_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "domain.h"
+#include "force.h"
+#include "memory.h"
+#include "math_const.h"
+
+#include
+#include
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define SMALLQ 0.00001
+#if defined(FFT_SINGLE)
+#define ZEROF 0.0f
+#else
+#define ZEROF 0.0
+#endif
+
+#define EPS_HOC 1.0e-7
+
+/* ---------------------------------------------------------------------- */
+
+PPPMCGOMP::PPPMCGOMP(LAMMPS *lmp, int narg, char **arg) :
+ PPPMCG(lmp, narg, arg), ThrOMP(lmp, THR_KSPACE)
+{
+}
+
+/* ----------------------------------------------------------------------
+ allocate memory that depends on # of K-vectors and order
+------------------------------------------------------------------------- */
+
+void PPPMCGOMP::allocate()
+{
+ PPPMCG::allocate();
+
+ const int nthreads = comm->nthreads;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
+ {
+#if defined(_OPENMP)
+ const int tid = omp_get_thread_num();
+#else
+ const int tid = 0;
+#endif
+
+ FFT_SCALAR **rho1d_thr;
+ memory->create2d_offset(rho1d_thr,3,-order/2,order/2,"pppm:rho1d_thr");
+ ThrData *thr = fix->get_thr(tid);
+ thr->init_pppm(static_cast(rho1d_thr));
+ }
+
+ const int nzend = (nzhi_out-nzlo_out+1)*nthreads + nzlo_out -1;
+
+ // reallocate density brick, so it fits our needs
+ memory->destroy3d_offset(density_brick,nzlo_out,nylo_out,nxlo_out);
+ memory->create3d_offset(density_brick,nzlo_out,nzend,nylo_out,nyhi_out,
+ nxlo_out,nxhi_out,"pppm:density_brick");
+}
+
+// NOTE: special version of reduce_data for FFT_SCALAR data type.
+// reduce per thread data into the first part of the data
+// array that is used for the non-threaded parts and reset
+// the temporary storage to 0.0. this routine depends on
+// multi-dimensional arrays like force stored in this order
+// x1,y1,z1,x2,y2,z2,...
+// we need to post a barrier to wait until all threads are done
+// writing to the array.
+static void data_reduce_fft(FFT_SCALAR *dall, int nall, int nthreads, int ndim, int tid)
+{
+#if defined(_OPENMP)
+ // NOOP in non-threaded execution.
+ if (nthreads == 1) return;
+#pragma omp barrier
+ {
+ const int nvals = ndim*nall;
+ const int idelta = nvals/nthreads + 1;
+ const int ifrom = tid*idelta;
+ const int ito = ((ifrom + idelta) > nvals) ? nvals : (ifrom + idelta);
+
+ // this if protects against having more threads than atoms
+ if (ifrom < nall) {
+ for (int m = ifrom; m < ito; ++m) {
+ for (int n = 1; n < nthreads; ++n) {
+ dall[m] += dall[n*nvals + m];
+ dall[n*nvals + m] = 0.0;
+ }
+ }
+ }
+ }
+#else
+ // NOOP in non-threaded execution.
+ return;
+#endif
+}
+
+/* ----------------------------------------------------------------------
+ free memory that depends on # of K-vectors and order
+------------------------------------------------------------------------- */
+
+void PPPMCGOMP::deallocate()
+{
+ PPPMCG::deallocate();
+ for (int i=0; i < comm->nthreads; ++i) {
+ ThrData * thr = fix->get_thr(i);
+ FFT_SCALAR ** rho1d_thr = static_cast(thr->get_rho1d());
+ memory->destroy2d_offset(rho1d_thr,-order/2);
+ }
+}
+
+/* ----------------------------------------------------------------------
+ adjust PPPMCG coeffs, called initially and whenever volume has changed
+------------------------------------------------------------------------- */
+
+void PPPMCGOMP::setup()
+{
+ int i,j,k,n;
+ double *prd;
+
+ // volume-dependent factors
+ // adjust z dimension for 2d slab PPPM
+ // z dimension for 3d PPPM is zprd since slab_volfactor = 1.0
+
+ if (triclinic == 0) prd = domain->prd;
+ else prd = domain->prd_lamda;
+
+ const double xprd = prd[0];
+ const double yprd = prd[1];
+ const double zprd = prd[2];
+ const double zprd_slab = zprd*slab_volfactor;
+ volume = xprd * yprd * zprd_slab;
+
+ delxinv = nx_pppm/xprd;
+ delyinv = ny_pppm/yprd;
+ delzinv = nz_pppm/zprd_slab;
+
+ delvolinv = delxinv*delyinv*delzinv;
+
+ const double unitkx = (2.0*MY_PI/xprd);
+ const double unitky = (2.0*MY_PI/yprd);
+ const double unitkz = (2.0*MY_PI/zprd_slab);
+
+ // fkx,fky,fkz for my FFT grid pts
+
+ double per;
+
+ for (i = nxlo_fft; i <= nxhi_fft; i++) {
+ per = i - nx_pppm*(2*i/nx_pppm);
+ fkx[i] = unitkx*per;
+ }
+
+ for (i = nylo_fft; i <= nyhi_fft; i++) {
+ per = i - ny_pppm*(2*i/ny_pppm);
+ fky[i] = unitky*per;
+ }
+
+ for (i = nzlo_fft; i <= nzhi_fft; i++) {
+ per = i - nz_pppm*(2*i/nz_pppm);
+ fkz[i] = unitkz*per;
+ }
+
+ // virial coefficients
+
+ double sqk,vterm;
+
+ n = 0;
+ for (k = nzlo_fft; k <= nzhi_fft; k++) {
+ for (j = nylo_fft; j <= nyhi_fft; j++) {
+ for (i = nxlo_fft; i <= nxhi_fft; i++) {
+ sqk = fkx[i]*fkx[i] + fky[j]*fky[j] + fkz[k]*fkz[k];
+ if (sqk == 0.0) {
+ vg[n][0] = 0.0;
+ vg[n][1] = 0.0;
+ vg[n][2] = 0.0;
+ vg[n][3] = 0.0;
+ vg[n][4] = 0.0;
+ vg[n][5] = 0.0;
+ } else {
+ vterm = -2.0 * (1.0/sqk + 0.25/(g_ewald*g_ewald));
+ vg[n][0] = 1.0 + vterm*fkx[i]*fkx[i];
+ vg[n][1] = 1.0 + vterm*fky[j]*fky[j];
+ vg[n][2] = 1.0 + vterm*fkz[k]*fkz[k];
+ vg[n][3] = vterm*fkx[i]*fky[j];
+ vg[n][4] = vterm*fkx[i]*fkz[k];
+ vg[n][5] = vterm*fky[j]*fkz[k];
+ }
+ n++;
+ }
+ }
+ }
+
+ // modified (Hockney-Eastwood) Coulomb Green's function
+
+ const int nbx = static_cast ((g_ewald*xprd/(MY_PI*nx_pppm)) *
+ pow(-log(EPS_HOC),0.25));
+ const int nby = static_cast ((g_ewald*yprd/(MY_PI*ny_pppm)) *
+ pow(-log(EPS_HOC),0.25));
+ const int nbz = static_cast ((g_ewald*zprd_slab/(MY_PI*nz_pppm)) *
+ pow(-log(EPS_HOC),0.25));
+
+ const double form = 1.0;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
+ {
+ int tid,nn,nnfrom,nnto,nx,ny,nz,k,l,m;
+ double snx,sny,snz,sqk;
+ double argx,argy,argz,wx,wy,wz,sx,sy,sz,qx,qy,qz;
+ double sum1,dot1,dot2;
+ double numerator,denominator;
+ const double gew2 = -4.0*g_ewald*g_ewald;
+
+ const int nnx = nxhi_fft-nxlo_fft+1;
+ const int nny = nyhi_fft-nylo_fft+1;
+
+ loop_setup_thr(nnfrom, nnto, tid, nfft, comm->nthreads);
+
+ for (m = nzlo_fft; m <= nzhi_fft; m++) {
+
+ const double fkzm = fkz[m];
+ snz = sin(0.5*fkzm*zprd_slab/nz_pppm);
+ snz *= snz;
+
+ for (l = nylo_fft; l <= nyhi_fft; l++) {
+ const double fkyl = fky[l];
+ sny = sin(0.5*fkyl*yprd/ny_pppm);
+ sny *= sny;
+
+ for (k = nxlo_fft; k <= nxhi_fft; k++) {
+
+ /* only compute the part designated to this thread */
+ nn = k-nxlo_fft + nnx*(l-nylo_fft + nny*(m-nzlo_fft));
+ if ((nn < nnfrom) || (nn >=nnto)) continue;
+
+ const double fkxk = fkx[k];
+ snx = sin(0.5*fkxk*xprd/nx_pppm);
+ snx *= snx;
+
+ sqk = fkxk*fkxk + fkyl*fkyl + fkzm*fkzm;
+
+ if (sqk != 0.0) {
+ numerator = form*MY_4PI/sqk;
+ denominator = gf_denom(snx,sny,snz);
+ sum1 = 0.0;
+ for (nx = -nbx; nx <= nbx; nx++) {
+ qx = fkxk + unitkx*nx_pppm*nx;
+ sx = exp(qx*qx/gew2);
+ wx = 1.0;
+ argx = 0.5*qx*xprd/nx_pppm;
+ if (argx != 0.0) wx = pow(sin(argx)/argx,order);
+ wx *=wx;
+
+ for (ny = -nby; ny <= nby; ny++) {
+ qy = fkyl + unitky*ny_pppm*ny;
+ sy = exp(qy*qy/gew2);
+ wy = 1.0;
+ argy = 0.5*qy*yprd/ny_pppm;
+ if (argy != 0.0) wy = pow(sin(argy)/argy,order);
+ wy *= wy;
+
+ for (nz = -nbz; nz <= nbz; nz++) {
+ qz = fkzm + unitkz*nz_pppm*nz;
+ sz = exp(qz*qz/gew2);
+ wz = 1.0;
+ argz = 0.5*qz*zprd_slab/nz_pppm;
+ if (argz != 0.0) wz = pow(sin(argz)/argz,order);
+ wz *= wz;
+
+ dot1 = fkxk*qx + fkyl*qy + fkzm*qz;
+ dot2 = qx*qx+qy*qy+qz*qz;
+ sum1 += (dot1/dot2) * sx*sy*sz * wx*wy*wz;
+ }
+ }
+ }
+ greensfn[nn] = numerator*sum1/denominator;
+ } else greensfn[nn] = 0.0;
+ }
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ run the regular toplevel compute method from plain PPPPMCG
+ which will have individual methods replaced by our threaded
+ versions and then call the obligatory force reduction.
+------------------------------------------------------------------------- */
+
+void PPPMCGOMP::compute(int eflag, int vflag)
+{
+
+ PPPMCG::compute(eflag,vflag);
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+#if defined(_OPENMP)
+ const int tid = omp_get_thread_num();
+#else
+ const int tid = 0;
+#endif
+ ThrData *thr = fix->get_thr(tid);
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+/* ----------------------------------------------------------------------
+ create discretized "density" on section of global grid due to my particles
+ density(x,y,z) = charge "density" at grid points of my 3d brick
+ (nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts)
+ in global grid
+------------------------------------------------------------------------- */
+
+void PPPMCGOMP::make_rho()
+{
+ const double * const q = atom->q;
+ const double * const * const x = atom->x;
+ const int nthreads = comm->nthreads;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+ {
+ // each thread works on a fixed chunk of atoms.
+ const int tid = omp_get_thread_num();
+ const int inum = num_charged;
+ const int idelta = 1 + inum/nthreads;
+ const int ifrom = tid*idelta;
+ const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
+#else
+ const int tid = 0;
+ const int ifrom = 0;
+ const int ito = num_charged;
+#endif
+
+ int i,j,l,m,n,nx,ny,nz,mx,my,mz;
+ FFT_SCALAR dx,dy,dz,x0,y0,z0;
+
+ // set up clear 3d density array
+ const int nzoffs = (nzhi_out-nzlo_out+1)*tid;
+ FFT_SCALAR * const * const * const db = &(density_brick[nzoffs]);
+ memset(&(db[nzlo_out][nylo_out][nxlo_out]),0,ngrid*sizeof(FFT_SCALAR));
+
+ ThrData *thr = fix->get_thr(tid);
+ FFT_SCALAR * const * const r1d = static_cast(thr->get_rho1d());
+
+ // loop over my charges, add their contribution to nearby grid points
+ // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
+ // (dx,dy,dz) = distance to "lower left" grid pt
+ // (mx,my,mz) = global coords of moving stencil pt
+
+ // this if protects against having more threads than charged local atoms
+ if (ifrom < num_charged) {
+ for (j = ifrom; j < ito; j++) {
+ i = is_charged[j];
+
+ nx = part2grid[i][0];
+ ny = part2grid[i][1];
+ nz = part2grid[i][2];
+ dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
+ dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
+ dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
+
+ compute_rho1d_thr(r1d,dx,dy,dz);
+
+ z0 = delvolinv * q[i];
+ for (n = nlower; n <= nupper; n++) {
+ mz = n+nz;
+ y0 = z0*r1d[2][n];
+ for (m = nlower; m <= nupper; m++) {
+ my = m+ny;
+ x0 = y0*r1d[1][m];
+ for (l = nlower; l <= nupper; l++) {
+ mx = l+nx;
+ db[mz][my][mx] += x0*r1d[0][l];
+ }
+ }
+ }
+ }
+ }
+#if defined(_OPENMP)
+ // reduce 3d density array
+ if (nthreads > 1) {
+ data_reduce_fft(&(density_brick[nzlo_out][nylo_out][nxlo_out]),ngrid,nthreads,1,tid);
+ }
+ }
+#endif
+}
+
+/* ----------------------------------------------------------------------
+ interpolate from grid to get electric field & force on my particles
+------------------------------------------------------------------------- */
+
+void PPPMCGOMP::fieldforce()
+{
+ // loop over my charges, interpolate electric field from nearby grid points
+ // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
+ // (dx,dy,dz) = distance to "lower left" grid pt
+ // (mx,my,mz) = global coords of moving stencil pt
+ // ek = 3 components of E-field on particle
+
+ const double * const q = atom->q;
+ const double * const * const x = atom->x;
+ const int nthreads = comm->nthreads;
+ const double qqrd2e = force->qqrd2e;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+ {
+ // each thread works on a fixed chunk of atoms.
+ const int tid = omp_get_thread_num();
+ const int inum = num_charged;
+ const int idelta = 1 + inum/nthreads;
+ const int ifrom = tid*idelta;
+ const int ito = ((ifrom + idelta) > inum) ? inum : ifrom + idelta;
+#else
+ const int ifrom = 0;
+ const int ito = num_charged;
+ const int tid = 0;
+#endif
+ ThrData *thr = fix->get_thr(tid);
+ double * const * const f = thr->get_f();
+ FFT_SCALAR * const * const r1d = static_cast(thr->get_rho1d());
+
+ int i,j,l,m,n,nx,ny,nz,mx,my,mz;
+ FFT_SCALAR dx,dy,dz,x0,y0,z0;
+ FFT_SCALAR ekx,eky,ekz;
+
+ // this if protects against having more threads than charged local atoms
+ if (ifrom < num_charged) {
+ for (j = ifrom; j < ito; j++) {
+ i = is_charged[j];
+
+ nx = part2grid[i][0];
+ ny = part2grid[i][1];
+ nz = part2grid[i][2];
+ dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
+ dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
+ dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
+
+ compute_rho1d_thr(r1d,dx,dy,dz);
+
+ ekx = eky = ekz = ZEROF;
+ for (n = nlower; n <= nupper; n++) {
+ mz = n+nz;
+ z0 = r1d[2][n];
+ for (m = nlower; m <= nupper; m++) {
+ my = m+ny;
+ y0 = z0*r1d[1][m];
+ for (l = nlower; l <= nupper; l++) {
+ mx = l+nx;
+ x0 = y0*r1d[0][l];
+ ekx -= x0*vdx_brick[mz][my][mx];
+ eky -= x0*vdy_brick[mz][my][mx];
+ ekz -= x0*vdz_brick[mz][my][mx];
+ }
+ }
+ }
+
+ // convert E-field to force
+ const double qfactor = qqrd2e*scale*q[i];
+ f[i][0] += qfactor*ekx;
+ f[i][1] += qfactor*eky;
+ f[i][2] += qfactor*ekz;
+ }
+ }
+#if defined(_OPENMP)
+ }
+#endif
+}
+
+/* ----------------------------------------------------------------------
+ charge assignment into rho1d
+ dx,dy,dz = distance of particle from "lower left" grid point
+------------------------------------------------------------------------- */
+void PPPMCGOMP::compute_rho1d_thr(FFT_SCALAR * const * const r1d, const FFT_SCALAR &dx,
+ const FFT_SCALAR &dy, const FFT_SCALAR &dz)
+{
+ int k,l;
+ FFT_SCALAR r1,r2,r3;
+
+ for (k = (1-order)/2; k <= order/2; k++) {
+ r1 = r2 = r3 = ZEROF;
+
+ for (l = order-1; l >= 0; l--) {
+ r1 = rho_coeff[l][k] + r1*dx;
+ r2 = rho_coeff[l][k] + r2*dy;
+ r3 = rho_coeff[l][k] + r3*dz;
+ }
+ r1d[0][k] = r1;
+ r1d[1][k] = r2;
+ r1d[2][k] = r3;
+ }
+}
+
diff --git a/src/USER-OMP/pppm_cg_omp.h b/src/USER-OMP/pppm_cg_omp.h
new file mode 100644
index 0000000000..4c32630271
--- /dev/null
+++ b/src/USER-OMP/pppm_cg_omp.h
@@ -0,0 +1,51 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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 KSPACE_CLASS
+
+KSpaceStyle(pppm/cg/omp,PPPMCGOMP)
+
+#else
+
+#ifndef LMP_PPPM_CG_OMP_H
+#define LMP_PPPM_CG_OMP_H
+
+#include "pppm_cg.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+ class PPPMCGOMP : public PPPMCG, public ThrOMP {
+ public:
+ PPPMCGOMP(class LAMMPS *, int, char **);
+ virtual void compute(int, int);
+ virtual void setup();
+ virtual ~PPPMCGOMP () {};
+
+ protected:
+ virtual void allocate();
+ virtual void deallocate();
+ virtual void fieldforce();
+ virtual void make_rho();
+
+ void compute_rho1d_thr(FFT_SCALAR * const * const, const FFT_SCALAR &,
+ const FFT_SCALAR &, const FFT_SCALAR &);
+// void compute_rho_coeff();
+// void slabcorr(int);
+
+};
+
+}
+
+#endif
+#endif
From a9f695113b108248daecbcffc709782a951d7202 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 15:25:06 +0000
Subject: [PATCH 03/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7467
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/USER-REAXC/pair_reax_c.cpp | 2 +-
src/USER-REAXC/pair_reax_c.h | 2 +-
2 files changed, 2 insertions(+), 2 deletions(-)
diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp
index a83ee1b5cc..f67100d789 100644
--- a/src/USER-REAXC/pair_reax_c.cpp
+++ b/src/USER-REAXC/pair_reax_c.cpp
@@ -763,7 +763,7 @@ void PairReaxC::read_reax_forces()
/* ---------------------------------------------------------------------- */
-void *PairReaxC::extract(char *str, int &dim)
+void *PairReaxC::extract(const char *str, int &dim)
{
dim = 1;
if (strcmp(str,"chi") == 0 && chi) {
diff --git a/src/USER-REAXC/pair_reax_c.h b/src/USER-REAXC/pair_reax_c.h
index 8fd14f6053..cfd9e35937 100644
--- a/src/USER-REAXC/pair_reax_c.h
+++ b/src/USER-REAXC/pair_reax_c.h
@@ -44,7 +44,7 @@ class PairReaxC : public Pair {
void coeff(int, char **);
void init_style();
double init_one(int, int);
- void *extract(char *, int &);
+ void *extract(const char *, int &);
int fixbond_flag, fixspecies_flag;
control_params *control;
From 46d332f7119864129db7a788103f6f9f7d17604d Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 15:28:00 +0000
Subject: [PATCH 04/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7468
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/version.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/version.h b/src/version.h
index cb03811359..abdad986fa 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "7 Jan 2012"
+#define LAMMPS_VERSION "8 Jan 2012"
From 75a49712c772be25fc11f9a80b1bb1d7bac95754 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 16:19:08 +0000
Subject: [PATCH 05/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7470
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/neigh_derive.cpp | 8 +++---
src/neigh_respa.cpp | 60 ++++++++++++++++++++++++--------------------
2 files changed, 37 insertions(+), 31 deletions(-)
diff --git a/src/neigh_derive.cpp b/src/neigh_derive.cpp
index d8e3c3edfb..7c83088b82 100644
--- a/src/neigh_derive.cpp
+++ b/src/neigh_derive.cpp
@@ -470,15 +470,15 @@ void Neighbor::skip_from_respa(NeighList *list)
firstneigh_inner[i] = neighptr_inner;
numneigh_inner[i] = n_inner;
npnt_inner += n_inner;
- if (npnt_inner >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_inner > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
if (respamiddle) {
firstneigh_middle[i] = neighptr_middle;
numneigh_middle[i] = n_middle;
npnt_middle += n_middle;
- if (npnt_middle >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_middle > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
}
diff --git a/src/neigh_respa.cpp b/src/neigh_respa.cpp
index 0af9eeec7a..c71d3c583a 100644
--- a/src/neigh_respa.cpp
+++ b/src/neigh_respa.cpp
@@ -32,8 +32,6 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
- int which = 0;
-
// loop over each atom, storing neighbors
int **special = atom->special;
@@ -82,6 +80,8 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
int npage_middle = 0;
int npnt_middle = 0;
+ int which = 0;
+
for (i = 0; i < nlocal; i++) {
if (pgsize - npnt < oneatom) {
@@ -159,16 +159,16 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
firstneigh_inner[i] = neighptr_inner;
numneigh_inner[i] = n_inner;
npnt_inner += n_inner;
- if (npnt_inner >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_inner > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
if (respamiddle) {
ilist_middle[inum] = i;
firstneigh_middle[i] = neighptr_middle;
numneigh_middle[i] = n_middle;
npnt_middle += n_middle;
- if (npnt_middle >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_middle > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
inum++;
@@ -189,7 +189,7 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
void Neighbor::respa_nsq_newton(NeighList *list)
{
- int i,j,n,itype,jtype,itag,jtag,which,n_inner,n_middle,bitmask;
+ int i,j,n,itype,jtype,itag,jtag,n_inner,n_middle,bitmask;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
@@ -241,6 +241,8 @@ void Neighbor::respa_nsq_newton(NeighList *list)
int npage_middle = 0;
int npnt_middle = 0;
+ int which = 0;
+
for (i = 0; i < nlocal; i++) {
if (pgsize - npnt < oneatom) {
@@ -336,16 +338,16 @@ void Neighbor::respa_nsq_newton(NeighList *list)
firstneigh_inner[i] = neighptr_inner;
numneigh_inner[i] = n_inner;
npnt_inner += n_inner;
- if (npnt_inner >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_inner > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
if (respamiddle) {
ilist_middle[inum] = i;
firstneigh_middle[i] = neighptr_middle;
numneigh_middle[i] = n_middle;
npnt_middle += n_middle;
- if (npnt_middle >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_middle > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
inum++;
@@ -366,7 +368,7 @@ void Neighbor::respa_nsq_newton(NeighList *list)
void Neighbor::respa_bin_no_newton(NeighList *list)
{
- int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle;
+ int i,j,k,n,itype,jtype,ibin,n_inner,n_middle;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
@@ -420,6 +422,8 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
int npage_middle = 0;
int npnt_middle = 0;
+ int which = 0;
+
for (i = 0; i < nlocal; i++) {
if (pgsize - npnt < oneatom) {
@@ -506,16 +510,16 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
firstneigh_inner[i] = neighptr_inner;
numneigh_inner[i] = n_inner;
npnt_inner += n_inner;
- if (npnt_inner >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_inner > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
if (respamiddle) {
ilist_middle[inum] = i;
firstneigh_middle[i] = neighptr_middle;
numneigh_middle[i] = n_middle;
npnt_middle += n_middle;
- if (npnt_middle >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_middle > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
inum++;
@@ -539,8 +543,6 @@ void Neighbor::respa_bin_newton(NeighList *list)
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
- int which = 0;
-
// bin local & ghost atoms
bin_atoms();
@@ -591,6 +593,8 @@ void Neighbor::respa_bin_newton(NeighList *list)
int npage_middle = 0;
int npnt_middle = 0;
+ int which = 0;
+
for (i = 0; i < nlocal; i++) {
if (pgsize - npnt < oneatom) {
@@ -713,16 +717,16 @@ void Neighbor::respa_bin_newton(NeighList *list)
firstneigh_inner[i] = neighptr_inner;
numneigh_inner[i] = n_inner;
npnt_inner += n_inner;
- if (npnt_inner >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_inner > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
if (respamiddle) {
ilist_middle[inum] = i;
firstneigh_middle[i] = neighptr_middle;
numneigh_middle[i] = n_middle;
npnt_middle += n_middle;
- if (npnt_middle >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_middle > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
inum++;
@@ -742,7 +746,7 @@ void Neighbor::respa_bin_newton(NeighList *list)
void Neighbor::respa_bin_newton_tri(NeighList *list)
{
- int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle;
+ int i,j,k,n,itype,jtype,ibin,n_inner,n_middle;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
@@ -796,6 +800,8 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
int npage_middle = 0;
int npnt_middle = 0;
+ int which = 0;
+
for (i = 0; i < nlocal; i++) {
if (pgsize - npnt < oneatom) {
@@ -890,16 +896,16 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
firstneigh_inner[i] = neighptr_inner;
numneigh_inner[i] = n_inner;
npnt_inner += n_inner;
- if (npnt_inner >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_inner > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
if (respamiddle) {
ilist_middle[inum] = i;
firstneigh_middle[i] = neighptr_middle;
numneigh_middle[i] = n_middle;
npnt_middle += n_middle;
- if (npnt_middle >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ if (n_middle > oneatom)
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
inum++;
From 8c1481250b572a54812884c4b829bd3f0499482a Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 17:04:58 +0000
Subject: [PATCH 06/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7471
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/image.cpp | 72 ++++++++++++++++++++++++++++++---------------------
1 file changed, 42 insertions(+), 30 deletions(-)
diff --git a/src/image.cpp b/src/image.cpp
index 897ee147bc..9fa87eae93 100644
--- a/src/image.cpp
+++ b/src/image.cpp
@@ -37,6 +37,7 @@ using namespace MathConst;
#define NCOLORS 140
#define NELEMENTS 109
#define BIG 1.0e20
+#define EPSILON 1.0e-6
enum{NUMERIC,MINVALUE,MAXVALUE};
enum{CONTINUOUS,DISCRETE,SEQUENTIAL};
@@ -164,6 +165,43 @@ void Image::view_params(double boxxlo, double boxxhi, double boxylo,
camDir[1] = sin(theta)*sin(phi);
camDir[2] = cos(theta);
+ // normalize up vector
+
+ if (up[0] == 0.0 && up[1] == 0.0 && up[2] == 0.0)
+ error->all(FLERR,"Invalid image up vector");
+ MathExtra::norm3(up);
+
+ // adjust camDir by epsilon if camDir and up are parallel
+ // do this by tweaking view direction, not up direction
+ // try to insure continuous images as changing view passes thru up
+ // sufficient to handle common cases where theta = 0 or 180 is degenerate?
+
+ double dot = MathExtra::dot3(up,camDir);
+ if (fabs(dot) > 1.0-EPSILON) {
+ if (theta == 0.0) {
+ camDir[0] = sin(EPSILON)*cos(phi);
+ camDir[1] = sin(EPSILON)*sin(phi);
+ camDir[2] = cos(EPSILON);
+ } else if (theta == MY_PI) {
+ camDir[0] = sin(theta-EPSILON)*cos(phi);
+ camDir[1] = sin(theta-EPSILON)*sin(phi);
+ camDir[2] = cos(theta-EPSILON);
+ } else {
+ camDir[0] = sin(theta+EPSILON)*cos(phi);
+ camDir[1] = sin(theta+EPSILON)*sin(phi);
+ camDir[2] = cos(theta+EPSILON);
+ }
+ }
+
+ // camUp = camDir x (Up x camDir)
+
+ MathExtra::cross3(up,camDir,camRight);
+ MathExtra::norm3(camRight);
+ MathExtra::cross3(camDir,camRight,camUp);
+ if (camUp[0] == 0.0 && camUp[1] == 0.0 && camUp[2] == 0.0)
+ error->all(FLERR,"Invalid image up vector");
+ MathExtra::norm3(camUp);
+
// zdist = camera distance = function of zoom & bounding box
// camPos = camera position = function of camDir and zdist
@@ -182,32 +220,6 @@ void Image::view_params(double boxxlo, double boxxhi, double boxylo,
camPos[1] = camDir[1] * zdist;
camPos[2] = camDir[2] * zdist;
- // normalize up vector
-
- if (up[0] == 0.0 && up[1] == 0.0 && up[2] == 0.0)
- error->all(FLERR,"Invalid image up vector");
- MathExtra::norm3(up);
-
- // rotate up if camDir and up point in same direction
- // should also check for opposite direction
- // NOTE: this is not yet right for preserving view as rotate thru up
-
- if (camDir[0] == up[0] && camDir[1] == up[1] && camDir[2] == up[2]) {
- double tmp = up[0];
- up[0] = up[1];
- up[1] = up[2];
- up[2] = tmp;
- }
-
- // camUp = camDir x (Up x camDir)
-
- MathExtra::cross3(up,camDir,camRight);
- MathExtra::norm3(camRight);
- MathExtra::cross3(camDir,camRight,camUp);
- if (camUp[0] == 0.0 && camUp[1] == 0.0 && camUp[2] == 0.0)
- error->all(FLERR,"Invalid image up vector");
- MathExtra::norm3(camUp);
-
// light directions in terms of -camDir = z
keyLightDir[0] = cos(keyLightTheta) * sin(keyLightPhi);
@@ -932,13 +944,13 @@ void Image::compute_SSAO()
double sy = surfaceBuffer[index * 2 + 1];
double sin_t = -sqrt(sx*sx + sy*sy);
- double theta = random->uniform() * SSAOJitter;
+ double mytheta = random->uniform() * SSAOJitter;
double ao = 0.0;
for (s = 0; s < SSAOSamples; s ++) {
- double hx = cos(theta);
- double hy = sin(theta);
- theta += delTheta;
+ double hx = cos(mytheta);
+ double hy = sin(mytheta);
+ mytheta += delTheta;
// multiply by z cross surface tangent
// so that dot (aka cos) works here
From fc5482c2401c9181110a9421b754aa6b75d219ab Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 17:10:34 +0000
Subject: [PATCH 07/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7472
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/version.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/version.h b/src/version.h
index abdad986fa..78003d8d4f 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "8 Jan 2012"
+#define LAMMPS_VERSION "9 Jan 2012"
From 0fba6fd65816dab8097322c88ad9f4b7e87c76cc Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 17:49:15 +0000
Subject: [PATCH 08/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7476
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
doc/Section_commands.html | 12 ++--
doc/Section_commands.txt | 1 +
doc/pair_coeff.html | 1 +
doc/pair_coeff.txt | 1 +
doc/pair_lj_cut_smooth.html | 111 +++++++++++++++++++++++++++++++++
doc/pair_lj_smooth.html | 3 +-
doc/pair_lj_smooth.txt | 3 +-
doc/pair_lj_smooth_linear.html | 86 +++++++++++++++++++++++++
doc/pair_lj_smooth_linear.txt | 81 ++++++++++++++++++++++++
doc/pair_style.html | 1 +
doc/pair_style.txt | 1 +
11 files changed, 293 insertions(+), 8 deletions(-)
create mode 100644 doc/pair_lj_cut_smooth.html
create mode 100644 doc/pair_lj_smooth_linear.html
create mode 100644 doc/pair_lj_smooth_linear.txt
diff --git a/doc/Section_commands.html b/doc/Section_commands.html
index 49dd4a7fdc..e2c38ffaea 100644
--- a/doc/Section_commands.html
+++ b/doc/Section_commands.html
@@ -426,12 +426,12 @@ potentials. Click on the style itself for a full description:
| lj/charmm/coul/charmm/implicit | lj/charmm/coul/long | lj/class2 | lj/class2/coul/cut |
| lj/class2/coul/long | lj/cut | lj/cut/coul/cut | lj/cut/coul/debye |
| lj/cut/coul/long | lj/cut/coul/long/tip4p | lj/expand | lj/gromacs |
-| lj/gromacs/coul/gromacs | lj/smooth | lj96/cut | lubricate |
-| lubricate/poly | lubricateU | lubricateU/poly | meam |
-| morse | peri/lps | peri/pmb | reax |
-| rebo | resquared | soft | sw |
-| table | tersoff | tersoff/zbl | tri/lj |
-| yukawa | yukawa/colloid
+ |
| lj/gromacs/coul/gromacs | lj/smooth | lj/smooth/linear | lj96/cut |
+| lubricate | lubricate/poly | lubricateU | lubricateU/poly |
+| meam | morse | peri/lps | peri/pmb |
+| reax | rebo | resquared | soft |
+| sw | table | tersoff | tersoff/zbl |
+| tri/lj | yukawa | yukawa/colloid
|
These are pair styles contributed by users, which can be used if
diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt
index cd619d11fe..a976a5df28 100644
--- a/doc/Section_commands.txt
+++ b/doc/Section_commands.txt
@@ -664,6 +664,7 @@ potentials. Click on the style itself for a full description:
"lj/gromacs"_pair_gromacs.html,
"lj/gromacs/coul/gromacs"_pair_gromacs.html,
"lj/smooth"_pair_lj_smooth.html,
+"lj/smooth/linear"_pair_lj_smooth_linear.html,
"lj96/cut"_pair_lj96.html,
"lubricate"_pair_lubricate.html,
"lubricate/poly"_pair_lubricate.html,
diff --git a/doc/pair_coeff.html b/doc/pair_coeff.html
index 0ac1a684b6..54420db500 100644
--- a/doc/pair_coeff.html
+++ b/doc/pair_coeff.html
@@ -134,6 +134,7 @@ the pair_style command, and coefficients specified by the associated
pair_style lj/gromacs - GROMACS-style Lennard-Jones potential
pair_style lj/gromacs/coul/gromacs - GROMACS-style LJ and Coulombic potential
pair_style lj/smooth - smoothed Lennard-Jones potential
+pair_style lj/smooth/linear - linear smoothed Lennard-Jones potential
pair_style lj96/cut - Lennard-Jones 9/6 potential
pair_style lubricate - hydrodynamic lubrication forces
pair_style lubricate/poly - hydrodynamic lubrication forces with polydispersity
diff --git a/doc/pair_coeff.txt b/doc/pair_coeff.txt
index 1b66fe2a33..4902a8e214 100644
--- a/doc/pair_coeff.txt
+++ b/doc/pair_coeff.txt
@@ -131,6 +131,7 @@ the pair_style command, and coefficients specified by the associated
"pair_style lj/gromacs"_pair_gromacs.html - GROMACS-style Lennard-Jones potential
"pair_style lj/gromacs/coul/gromacs"_pair_gromacs.html - GROMACS-style LJ and Coulombic potential
"pair_style lj/smooth"_pair_lj_smooth.html - smoothed Lennard-Jones potential
+"pair_style lj/smooth/linear"_pair_lj_smooth_linear.html - linear smoothed Lennard-Jones potential
"pair_style lj96/cut"_pair_lj96.html - Lennard-Jones 9/6 potential
"pair_style lubricate"_pair_lubricate.html - hydrodynamic lubrication forces
"pair_style lubricate/poly"_pair_lubricate.html - hydrodynamic lubrication forces with polydispersity
diff --git a/doc/pair_lj_cut_smooth.html b/doc/pair_lj_cut_smooth.html
new file mode 100644
index 0000000000..1430313a80
--- /dev/null
+++ b/doc/pair_lj_cut_smooth.html
@@ -0,0 +1,111 @@
+
+LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands
+
+
+
+
+
+
+
+
+
+pair_style lj/cut/smooth command
+
+pair_style lj/cut/smooth/cuda command
+
+pair_style lj/cut/smooth/omp command
+
+Syntax:
+
+pair_style lj/cut/smooth Rc
+
+- Rc = cutoff for lj/cut/smooth interactions (distance units)
+
+Examples:
+
+pair_style lj/cut/smooth 5.456108274435118
+pair_coeff * * 0.7242785984051078 2.598146797350056
+pair_coeff 1 1 20.0 1.3 9.0
+
+Description:
+
+Style lj/cut/smooth computes a LJ interaction that combines the standard
+12/6 Lennard-Jones function and subtracts a linear term that includes a
+cutoff distance Rc.
+
+
+
+At the cutoff Rc, the energy and force (its 1st derivative) will be 0.0.
+
+The following coefficients must be defined for each pair of atoms
+types via the pair_coeff command as in the examples
+above, or in the data file or restart files read by the
+read_data or read_restart
+commands, or by mixing as described below:
+
+- epsilon (energy units)
+
- sigma (distance units)
+
- cutoff (distance units)
+
+If not specified, the global value for Rc is used.
+
+
+
+Styles with a cuda, gpu, omp, or opt suffix are functionally
+the same as the corresponding style without the suffix. They have
+been optimized to run faster, depending on your available hardware, as
+discussed in Section_accelerate of the
+manual. The accelerated styles take the same arguments and should
+produce the same results, except for round-off and precision issues.
+
+These accelerated styles are part of the USER-CUDA, GPU, USER-OMP and OPT
+packages, respectively. They are only enabled if LAMMPS was built with
+those packages. See the Making LAMMPS
+section for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the -suffix command-line
+switch when you invoke LAMMPS, or you can
+use the suffix command in your input script.
+
+See Section_accelerate of the manual for
+more instructions on how to use the accelerated styles effectively.
+
+
+
+Mixing, shift, table, tail correction, restart, rRESPA info:
+
+For atom type pairs I,J and I != J, the epsilon and sigma coefficients and
+cutoff distance can be mixed. The default mix value is geometric.
+See the "pair_modify" command for details.
+
+This pair style does not support the pair_modify shift
+option for the energy of the pair interaction.
+
+The pair_modify table option is not relevant
+for this pair style.
+
+This pair style does not support the pair_modify
+tail option for adding long-range tail corrections to energy and
+pressure, since the energy of the pair interaction is smoothed to 0.0
+at the cutoff.
+
+This pair style writes its information to binary restart
+files, so pair_style and pair_coeff commands do not need
+to be specified in an input script that reads a restart file.
+
+This pair style can only be used via the pair keyword of the
+run_style respa command. It does not support the
+inner, middle, outer keywords.
+
+
+
+Restrictions: none
+
+Related commands:
+
+pair_coeff
+
+Default: none
+
+
diff --git a/doc/pair_lj_smooth.html b/doc/pair_lj_smooth.html
index 57af043e75..9cefda53eb 100644
--- a/doc/pair_lj_smooth.html
+++ b/doc/pair_lj_smooth.html
@@ -121,7 +121,8 @@ to be specified in an input script that reads a restart file.
Related commands:
-pair_coeff
+
pair_coeff, pair
+lj/smooth/linear
Default: none
diff --git a/doc/pair_lj_smooth.txt b/doc/pair_lj_smooth.txt
index a2d37902c9..65dd323a3c 100644
--- a/doc/pair_lj_smooth.txt
+++ b/doc/pair_lj_smooth.txt
@@ -116,6 +116,7 @@ This pair style can only be used via the {pair} keyword of the
[Related commands:]
-"pair_coeff"_pair_coeff.html
+"pair_coeff"_pair_coeff.html, "pair
+lj/smooth/linear"_pair_lj_smooth_linear.html
[Default:] none
diff --git a/doc/pair_lj_smooth_linear.html b/doc/pair_lj_smooth_linear.html
new file mode 100644
index 0000000000..5a43b896d9
--- /dev/null
+++ b/doc/pair_lj_smooth_linear.html
@@ -0,0 +1,86 @@
+
+LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands
+
+
+
+
+
+
+
+
+
+pair_style lj/smooth/linear command
+
+Syntax:
+
+pair_style lj/smooth/linear Rc
+
+- Rc = cutoff for lj/smooth/linear interactions (distance units)
+
+Examples:
+
+pair_style lj/smooth/linear 5.456108274435118
+pair_coeff * * 0.7242785984051078 2.598146797350056
+pair_coeff 1 1 20.0 1.3 9.0
+
+Description:
+
+Style lj/smooth/linear computes a LJ interaction that combines the
+standard 12/6 Lennard-Jones function and subtracts a linear term that
+includes the cutoff distance Rc, as in this formula:
+
+
+
+At the cutoff Rc, the energy and force (its 1st derivative) will be 0.0.
+
+The following coefficients must be defined for each pair of atoms
+types via the pair_coeff command as in the examples
+above, or in the data file or restart files read by the
+read_data or read_restart
+commands, or by mixing as described below:
+
+- epsilon (energy units)
+
- sigma (distance units)
+
- cutoff (distance units)
+
+The last coefficient is optional. If not specified, the global value
+for Rc is used.
+
+
+
+Mixing, shift, table, tail correction, restart, rRESPA info:
+
+For atom type pairs I,J and I != J, the epsilon and sigma coefficients
+and cutoff distance can be mixed. The default mix value is geometric.
+See the "pair_modify" command for details.
+
+This pair style does not support the pair_modify
+shift option for the energy of the pair interaction.
+
+The pair_modify table option is not relevant for
+this pair style.
+
+This pair style does not support the pair_modify
+tail option for adding long-range tail corrections to energy and
+pressure, since the energy of the pair interaction is smoothed to 0.0
+at the cutoff.
+
+This pair style writes its information to binary restart
+files, so pair_style and pair_coeff commands do not need
+to be specified in an input script that reads a restart file.
+
+This pair style can only be used via the pair keyword of the
+run_style respa command. It does not support the
+inner, middle, outer keywords.
+
+
+
+Restrictions: none
+
+Related commands:
+
+pair_coeff, pair lj/smooth
+
+Default: none
+
+
diff --git a/doc/pair_lj_smooth_linear.txt b/doc/pair_lj_smooth_linear.txt
new file mode 100644
index 0000000000..81899d8a5a
--- /dev/null
+++ b/doc/pair_lj_smooth_linear.txt
@@ -0,0 +1,81 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+pair_style lj/smooth/linear command :h3
+
+[Syntax:]
+
+pair_style lj/smooth/linear Rc :pre
+
+Rc = cutoff for lj/smooth/linear interactions (distance units) :ul
+
+[Examples:]
+
+pair_style lj/smooth/linear 5.456108274435118
+pair_coeff * * 0.7242785984051078 2.598146797350056
+pair_coeff 1 1 20.0 1.3 9.0 :pre
+
+[Description:]
+
+Style {lj/smooth/linear} computes a LJ interaction that combines the
+standard 12/6 Lennard-Jones function and subtracts a linear term that
+includes the cutoff distance Rc, as in this formula:
+
+:c,image(Eqs/pair_lj_cut_smooth.jpg)
+
+At the cutoff Rc, the energy and force (its 1st derivative) will be 0.0.
+
+The following coefficients must be defined for each pair of atoms
+types via the "pair_coeff"_pair_coeff.html command as in the examples
+above, or in the data file or restart files read by the
+"read_data"_read_data.html or "read_restart"_read_restart.html
+commands, or by mixing as described below:
+
+epsilon (energy units)
+sigma (distance units)
+cutoff (distance units) :ul
+
+The last coefficient is optional. If not specified, the global value
+for Rc is used.
+
+:line
+
+[Mixing, shift, table, tail correction, restart, rRESPA info]:
+
+For atom type pairs I,J and I != J, the epsilon and sigma coefficients
+and cutoff distance can be mixed. The default mix value is geometric.
+See the "pair_modify" command for details.
+
+This pair style does not support the "pair_modify"_pair_modify.html
+shift option for the energy of the pair interaction.
+
+The "pair_modify"_pair_modify.html table option is not relevant for
+this pair style.
+
+This pair style does not support the "pair_modify"_pair_modify.html
+tail option for adding long-range tail corrections to energy and
+pressure, since the energy of the pair interaction is smoothed to 0.0
+at the cutoff.
+
+This pair style writes its information to "binary restart
+files"_restart.html, so pair_style and pair_coeff commands do not need
+to be specified in an input script that reads a restart file.
+
+This pair style can only be used via the {pair} keyword of the
+"run_style respa"_run_style.html command. It does not support the
+{inner}, {middle}, {outer} keywords.
+
+:line
+
+[Restrictions:] none
+
+[Related commands:]
+
+"pair_coeff"_pair_coeff.html, "pair lj/smooth"_pair_lj_smooth.html
+
+[Default:] none
diff --git a/doc/pair_style.html b/doc/pair_style.html
index c146e10b2f..a5ff6192c0 100644
--- a/doc/pair_style.html
+++ b/doc/pair_style.html
@@ -136,6 +136,7 @@ the pair_style command, and coefficients specified by the associated
pair_style lj/gromacs - GROMACS-style Lennard-Jones potential
pair_style lj/gromacs/coul/gromacs - GROMACS-style LJ and Coulombic potential
pair_style lj/smooth - smoothed Lennard-Jones potential
+pair_style lj/smooth/linear - linear smoothed Lennard-Jones potential
pair_style lj96/cut - Lennard-Jones 9/6 potential
pair_style lubricate - hydrodynamic lubrication forces
pair_style lubricate/poly - hydrodynamic lubrication forces with polydispersity
diff --git a/doc/pair_style.txt b/doc/pair_style.txt
index 6535e69e2d..1050b13d3c 100644
--- a/doc/pair_style.txt
+++ b/doc/pair_style.txt
@@ -133,6 +133,7 @@ the pair_style command, and coefficients specified by the associated
"pair_style lj/gromacs"_pair_gromacs.html - GROMACS-style Lennard-Jones potential
"pair_style lj/gromacs/coul/gromacs"_pair_gromacs.html - GROMACS-style LJ and Coulombic potential
"pair_style lj/smooth"_pair_lj_smooth.html - smoothed Lennard-Jones potential
+"pair_style lj/smooth/linear"_pair_lj_smooth_linear.html - linear smoothed Lennard-Jones potential
"pair_style lj96/cut"_pair_lj96.html - Lennard-Jones 9/6 potential
"pair_style lubricate"_pair_lubricate.html - hydrodynamic lubrication forces
"pair_style lubricate/poly"_pair_lubricate.html - hydrodynamic lubrication forces with polydispersity
From 411a21d6b0c98ca7798c997decbe2276bd996cb4 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 17:49:36 +0000
Subject: [PATCH 09/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7477
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
doc/Eqs/pair_lj_smooth_linear.jpg | Bin 0 -> 10372 bytes
doc/Eqs/pair_lj_smooth_linear.tex | 13 +++++++++++++
2 files changed, 13 insertions(+)
create mode 100644 doc/Eqs/pair_lj_smooth_linear.jpg
create mode 100644 doc/Eqs/pair_lj_smooth_linear.tex
diff --git a/doc/Eqs/pair_lj_smooth_linear.jpg b/doc/Eqs/pair_lj_smooth_linear.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..b0626abae1199aec35611c9974f7d63ece5b768a
GIT binary patch
literal 10372
zcmcI~cT`i|wr?;`UyWYP0B$QQDk%a82nYbruP?yW
z1VA2e6G%u%2)ub++`M^{h?ty&_!?AXWTfQORCn%BQ&CgX&@t1~+`UIjP0hf@a1R7x
zWo4zIe*k&F0%2xhWx4tVAR__f-7qI0paxvE0GI&;0N^!=Ya;)`^{|A%8w4cRy)SM9
z2ncTw+#tMhlbG}-AsOKM?HfSCo76;i*obNPG$VaPhC86L_lrHXu57
z7$UEgKE#gx{DtG8$TJPo501{6Wmn?>lD~+^04_V9{eQu$^X3@is@S_r0yy)c7plwz
z{cmGkYm`_RFcr}iy$(F~0D!la4AqjTw9Y|M2~zDMF0=b|iDQ_zez-7ITG;twupxQ<
z=RryBjRy~Sd%hj)J5v_e>2783j%7F%eJU*(cyeQjn4;1rk)`tw^)=|bN?~^qXdNb|
zd&0}2D^ycyC$2nA4%&GtU&VM5*_(KCO*ZiRxWoIA7t)zqHa=6~6T!Q-BHT?87^1G9
zkK(iZYTyIE4%?YJSw{F3EB>%`-wz*4^^;oc6~Lz0x5Ws4W3I>C;aE~C
zzv1oRHf}IAqF*>X^Y5FUbJ%kq_)KS$7MY@ajf}58n^!)j-
zQ^;Qt$t$A8Dak9R=V3c@deQNf2ALL_C1Vk`0i62=mu1xU6D%WE&MiBWMF6{Jw?~A!
zX*xyb_kX-TmEV21CflY9rNwkYpKl%|y81gvFUwM-8WE@6f(o^EiBVExwq3HSlbROR
zr#dZJFfdLV%FKuH8^GK{k^X*st@-li0k=CfwXDv9F}{gMAR^~roeu+TS6SMJDMIZs
zzpreoI?fdd5*0HBngI)(mJ6Cs{=^NFQ3QXdTG6Qy#j>YqMMVB
zSfcuq@|9$`opYrP+ZizmD%H1;XNv+73_9OX3NQ)XedfHrb=$Uqa>0Sokis6fty!R%
z2KJ&TvT#*}@P@%)U^2XttIFtX{Sn|{?~O(gcj7-Tj>i$}vj6wz-t9PmG@a*`e&}*8
zLsvv+%JIvGv9N9KQOnj>BeNKr!S4#YS}d;dF-Uu_BepQ0$qNJFR$PnB6hnq%v#;;V
zScT-@#>H0vg|6(bVFpJr{2^<1g=ADHv=ROD7Z|Ya+tY?ArmDeJ8amt5f$4d`p>tNc
zw!Q{IpQ$IJBY|CCO8k!QNGaVv(qv_TTq>JAx#W<2;V`#k&fLd>p3iZ!Nt)d>iSGK$
zY6M{&;Dl-ZNO~%MKHkfH=E%)HA!L#RQM>J~NHPr`WBDlei&**cY;5|2(_I#5Q%5;(
ztND5FSA)4;{*Cb|2(6p>lTw84Tif)&B?yU$*ET>;R8vv^Lv4&H#%<8f3H
zM#{{dm5g|9`CaA?8_Z1=`)9_@i?>UnyLSQqfiWJH_osuv6Sh+9gz4WH9nY)#!Pime
z(z?T(Wq=@~`grz1xfjri7+ikRj2k>kWH3Vk4%Zy1zLUivO&=PZKL21FB*3Keb;IRW
zHrcK)U66`b#|d0Y(4b9kGWGQmz&kd_qrhe;C;J1*a)F76f}py{drxM;rRYr($YW@r
zNd6lOJ3lj96_c&I<>P(boZ+cbjJiv2qOp^|6C+YxpUrt|#_?ii>f;gCg-;K8A&(PK
zWg4A4XqLjX5mJW3l(;)`IIKqPN|g68inVpeKuiTSYxJc5Thr+HGQB8TC25&vZCNyq
zW>ZO1?LTwVQ&}dc>_Y=K_Y-9;enTA^d7+aDfQ`QmIP&tlD;Dt0kNzO=<4UL3sQnew1yO|02K?RKx)@URwR^c1EzmF)v
zLowLg_*K#sfRvSeW^7Z2%lJ>bX-;0|o)^MOCjSn}W3pdddnY0RczGet^Y)z9r(>y&
z^(<+RH;TA
zyLeHQK_hM3_yv`3d-d75u=ke>`j&b7I_LVtcgJ6aXD5bSq;NI|c-O<2mg)1L>5B(v
zL7NgimY30PpSMC>d~S`t+uHQ`F2VhXp2K>f>Rh3=Pa|Q!_~|MmoVv1OM);o09ZxRa
z*1BQ1o@a=hu!tn0w|j>!zgL*&`DK|14!%sZ;jv$t9l0qF#~lr%>;NCNooIz7c?QjX
zU%$=AzqDIn2S@u?^*;=H$g+ev3JSh+B#`qGKn${#P?}>%81<4c3ML5ttpuKg95{G9
zarPi~CYfJjPp@UqT%6X`^!+=(i}c?;I41rWKe-UQq_8V7xp{hlTwe0L`ORW;3_G_}
ztwH0;E9fmHxY|J1$$=+W6_AS9AcZ18{%%I+=rv$+QT%}c>pQoci^$p3
zrCOQQwS~!jWE#PUJtaYTION%2Hb_;n$|%zawK8$*o_Va0M`%tPgjwbgKsO*rwq^ZE$t)Uk{LN+>
zct{!IWg!yjtL<{|c364wc5SpHo`{7)JcSZFJ_ZdOeGc~IeCm5=Io^4VP<48WsV3uk
zH92mfzaBff@1WY2%id(P;)5=Eo3VHuKrI<98E`x$*~nEMTfClTRZgApjc?dcy)llxQ-noNZ`G1}4SFC*lIzAbK~sJbDjkcR|V-e9qPB82PW*5NTAzmT&z>@9}l>4;za%!W*re`v$X%7DRmdC
zJ-lH`;bM1mus%Jq125&e0{klKo{~G$Pr3w)_x)Wq{
zHkISmlh9wbLEP@aD{2$~WD8=L=a+1F8RS1Szcyc_G3sg%@)yEoc1HGGQI{h|SAdXT
z3EN+m&pZ$#<|tvWCR5JA1o^yKs10Rxr0Cwx-U>q6lWd
z%$^_nzG9<}-$5w0G)uZMX2n+hUDv7u!SskkU)9?i1|mT|!g{UeKR3E1EUIP7va6B2y5jl#t6mP_}X_9+X&;dKpCAn6gYo4&0z8_iv+1l2)
z&=NT^Cg#`WR(xqz1jhjS9}gZpj*E~q@G+gm7!oTgwnlrYrm{W~-~ID!mWyY)c#mN+
zHTto9^CiqMb{lWWR>T;*9!0DQY}ls(4syM`;E)u{*?ITo3XtWUA0*Ih`xE}9dN$!P
zq5y@g+d|6wv_v->$iiVRVhjFw+bu(q$IT&gbgR2u2p2506P{8}KT-|MVo_~3?DHBC
z=&zT8W*78dQH1-7BKm)^ugoo-vfeIfVWE1oWWf}IW>dS
zmX|ku5W7_48_PS4IKO`!$NyOT`x%jTS1tOe1qwn^lw*y*#9HA)fb1&N$yiR;HMtE~
zW`f?}=$C#Bmj#yzSnQNqy$!PrdlK~OIGbC14=kOIJRlgipF(}>vSgWi(dtUVp4!$r
ztNmpx`(%a#z)^{Vck^bA-jZjb0-WGQy55s$0g
z(D`DQ+~2+Q4(r$2b}D2?EW(C1^OGi#aOS^HrIoPlAz$&0GRMxiz-bD~szMW<7mmRp
z3m!-PytR!H93AZg!mqV}{$D29f3`sX<>4r3?#M@)iK#Ia!DxI6$c8^S?FucH0l$S8
z)@{TWINn=AlNu4fByRk%HoNWk+^6EasD0;bV&$}QO#39y$$PV*k^h}|<@P!phde6D
zyKybj{t5rZ?LEfKd{LwdP-&)io%qpHUAE>*HCd$9srOO+R%OQgc4MTOXeCn#jN7N+
z3|H#PFU(#5Ob;(?*^2MBVP3vX;0{DUD}Pk{YsCL^&D+k>Ho6LN)W!tpErS5W)dqqH
zje6o#MLb=`8uxmT?BRw-B+7YXxg@*5F@dZnf%HHeDt~mt_XBE?@vTu6PD6}~C;sc2
zqps6?abD`;x^h+%=j;dnBb)x)!ykI6q#bMR^Ui{{v(F{@XvvuQ{8gJKd6$(epKy`P
zjLDQ_WP$BNs?2%@QLI0ML^Vtko2Z%9W#HLqf&?4)BpuMC_2KG7Cgk!eq1{i^7KRXG
z&c`JTYaigT3qu@c)Y4DEVOpYbhd)LSRSdE$TD
zB3@C;0LzZHZ%EkF%TReeQqywKz*h8Nf|Ga1Ns+=b`Swe*+uh7TBC$$+Gah1MOS|_xhgMOiTNH*&=nFhO0Z`T65i9b{nmeQEJ4)7
zXd)FH2%p9&?>~EHVByAO{w1P6bP6j~5=WUMlyzf9>iNxX!b|W{{aZR*B^=Qj@nJR>K(o;oq6$j$aJKU;B2+;e$Q
zg{}N1#__#fI*+El2E^s9N?aDsAoX0|YmbX<4*wHA{8EK8$*anXn**84RULhz50dUL
zH8tLJ&D0s1Ppft3x$^?r-wefUj*39iTi7!NTT(k2R%68`4z@CUJ^Eh%wWj|CI;vTN
z2f+__U@hkbwT#diJ{{4LyKFA%zrr4+g}Wr5Fa0!Fx=br#*knasYC0fwu-`hH-^g(;
zhn_G_l9#koEu79pFF4zcgl<;&8f@3()_$S6m16)4gp_*Y-lb&G$40K}pXIU$2@BiO%9Qd`n>|vv
z9V9{B!%npVlERaY{h8hwFE|dBD2>0A`aSwbKF$fZE9+4cY2Q2Cy
zE={S!OW?^HpOn#~Y*8>S|75eMKB7y6BV5}27v;juy@;OKQdi4tREI`|MXmf|*dzHQ
z99q^T_~-a+qiQ+om8NCz%MD(s9{a-X49?cLRzb>s1v-kJKftZEZ&%`gX@=toDSXWFAMbg!8W_G9>y8%lHbQCWtOnX
zmydOSd10)N#COQ(Db}lB(uTFk8Z493^FpN~Zkvc_
zDGM!12jkTlyx7~gI8i^q7<5eCizD%m7kS<0dT?qRX$4vGM0*D5X};-#WGO@@2i#ke
z;LA4`ft7i1sSxpA8xKEG4O;yuD7n@6{zlOPnvy^5$
zhSypI#$4zWz1k!TXJWLpSFEUj?wj;{kh}trcTRMTRd~{}d>E4SxM=VS~P;-aPAPr@f3|0-W054K(*%d^z;>phbg3(6^9L4(fQWC+$d+43m4_Sif>A4
zExXTwrukKu2fcH9Dbu9Yfti+HcLHkQlk4A4sQRW$%QOBV_LthHvm?L1pDqKr$x*Dk
z_l=1IyxOHZgQV54kP78|oi3*s0}5xw(q{eTV8Qx&2nbY8PvPPBeB*A-*>D2EGPV!bp>y-6#W=c_)
z*PW~0M|K}_a-GW9uT@G?|FUp*uFF5+`UEqBR&W>$A&_EvtKh+pe&xdaJu{Iu)@m!Q
zNek|mjB%CmW!8s}B^STD9h`YDS0Ziv!eF(hKJ)oU;ztMKH4jh=)z+gdVlh2<@d9Kn
zcXjRhOysA{D?k({#kN{ggb3tR{P+2{beb~B9|aMgYClTkiw`CpN?M^5#e#JYbro6p
zD9$fo^e9`NoG@E^RW3^-ZSeOdp;0=4wCF29m)*Cp0=FrdjTAm`<5yJd9C_`2e>tYo
zciK_rb^oD$iohb&RjUpvlwKdt}=otw;>!d+M@3iogld#RUQBlleo
z!0TAkBko6cT_3fnhR%!W6<9BI4dtquh{DdWe4l>TH<-P%lxCa{tP-J&?M9XJwJRrUtZrXX6`>Z{%y_
z;oov3J>jV}%~_?1a7sZta0XDEParDei{hqGzS7g|#PdtyKXsw4ZfTUU1{%pZpcTBQ
zcqA)7KJeguny>Ex_v3G56dEF?B3g`ZnIpe~3xr}7fzx+ASBz!Y%EsOl25e=$gkYrU
z95@>sfTy@Mch8fHG3VW_J*~uo(o@IM&wb7B`m5w2{frFQ^;sSxM-@NuQA)&xu#_fj
zWBTqD;HRFzd9LM!oQv^0H-Ea9i~^F05qfE)d+E8@oFdLMkWe$Pkkst
zBPx*w^Pw6(+#R2*p3q{CIzpO~1$uh54fJ>kV2SM!Anrs~EAwj3nH~#7%I?cqRvfVM
z4-;e;Yy}GRkR@ED&e!RffOCC+DTqgLVMSUvRz-`}0Vj_L1Q51dp}*af&D`jZe;ao>
zZ->-~2)P7d);0*KVzK6Lw(Uqmmb4pOI|M{ePj1y!N?5Zkik|a9G%MF8W%pXNcXCjncYeGfv5TYVB0My=Gp=u&hz=|(
zl_p3}pB~!2NSsSO>wh-6VoqQ{M00JYorCvOTsZn7bkw6N0sOH0c>khoj_KDB&tmC`emQ0C
zkj=i#ng10af8#6Q4+de5!@{t+(H@^9ad
zboY+hnR&BJN38D<&EI;&_T9C;7#4PxHV10GjCtM52EVWM%j~M#d=pV(LnNlT?3I%K
zsVDm@Q*a%bu#7l)?zSB4j?_Mq)NWEQv4|orDPBos_MJ0Zz1-d2o@nl30s7f*NOh^|
zARJEnmF}V<=wVi^gIs=B0B!Wv&iRdcuCiDYckl^CCv(G3dGroUX&4?$r*Ydqm6J9o
zt5gb^6gJZ#N>WhvO5$=1MgnV*K+!%$QR+o
zFV-)WEHUcnNwB@*r}+DGm+^K5b_5}(Q@|A&O}WM6(#(JnK{!Z3RRJ-X>0=(VyPx?F
z^etIs5}6=QYB#fyk$=P;nMVLiRjeii`XT9J!L&^u@TMj4$6RAsht+5;vDin~X(EP_gLnCgAUIY{D_VjX5~S+6--eSg
z3|I^!x09T03Wl@>o@fWbSiwot99_I-$%FxiEP*>th^hP8S4ZF@Ov-ddk0X)CFnsr)g`{30pT#8lf`uM6J
zqj(rwT|Uce9~HP_hSmY?tx#_@8rwfpP%3J&zXEh55`0Fy#b{-L0~io$v>#Tp9CM5!
zi^x}p?9qMVYOhD7hB3+3Lu>Nt$!>Wgp8`ZPLj5Bm)2xIFKz_iV0l6PMP)&RJzCsZ+
zFO`Go;|5(9TSu1U#4Z@7%xvldEgX@7xBe`h(?{%SGwH^tm-&==_>NsRO4_3bhejAi
zw+CjtnH8@9v}|_BIi8)+W8Bm7xE+`K-tKeNCQ6fFP<^n}Ca-iu)`knNoR&9cn@EoL
zZXz{&O7R;jl>BhD)2>h_Er>K82w&AB%J{O=KJa{62@*$8YtA^6&Y~saOF={@recvI
z-?vT^`$XV}A+J#suA+dT7Xpy*9&Vk*Jzlgb9oxv_31&~NLv)e_w;J5XC=kN-6Ytut
zjUXyz4)W9=JhlsRX2s;GH9>FYW9qkB>nfjG^SC=1dD75peV&mL*Pa7-Y6gbSlt@~l
z27k@u{!9`6?4#*+Z&@O}p`j2GQ~I(Uk_&mJiO`%qqFrQEu9jR>!Q~Idmkrx@z#}Ub
zt@_ZyrLs7g{RDz@EnkkgDrpr=FH85JGv0ThjVk%lRG;xC>W*E$7#o?BR!x#FQl$cKgeZlR5j(OcO#%jkj5
z!N!$f$%@8|zI>IrcP2pskO6*Siub+t0K3Xuog6^hLV;W398*rOeu-a8O|~do3HSw{
zwwJ}b8YiN_XVD@spEI~_c4+BO==0xTe{N&@G2n=tR^1zSR5PU(SY8Q1H%rh>2}<+*
zE){ZqtMni8^g-l86_*;Mxjm$N@O=lr(%qs@J&y9p%#qQFtM^x=BE|Q``$r*{okjVTfA=CwjpD*ondsX4J!Yspv4OLjAL6MG${
z?XL+qreRu_xA)3Svm)v17n}2jgp$EOk2JPtN6&ftFHYs1?F=)V@9(eISB;eHD&
zSls~yo%qSXS2(1IlLaZ_wX9
zdOwLN3M-rIz3e!#y8?Wyy#j1|>@MS5&(&X8A@21kRdWA+6=J|JbCN}kVI1DB+d#uK
zJm3RL`^0&m9$gr{gT~sf&l1l{-0PR3=rU9Blhokv
zQ_GFbxu)Qe;SFXYg0Z%W*8UL3R~IVippXxn7irTb<2$8w&7`0Py(8(IouqHf>r8WP
z!UmWfL0IDy<^$8){wKt0OfM*Qbw?{c(ueA27Q-iz+4&j&Wa0l*@3@g^6wd1rz-o!oIG^>2cHn=gpXAd%wJPby?hE9*;{+o;c
KAKo=yjr|v7Zw)d4
literal 0
HcmV?d00001
diff --git a/doc/Eqs/pair_lj_smooth_linear.tex b/doc/Eqs/pair_lj_smooth_linear.tex
new file mode 100644
index 0000000000..b8980f8d25
--- /dev/null
+++ b/doc/Eqs/pair_lj_smooth_linear.tex
@@ -0,0 +1,13 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+
+\begin{eqnarray*}
+\phi\left(r\right) & = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
+ \left(\frac{\sigma}{r}\right)^6 \right] \\
+E\left(r\right) & = & \phi\left(r\right) - \phi\left(R_c\right) - \left(r - R_c\right) \left.\frac{d\phi}{d r} \right|_{r=R_c} \qquad r < R_c
+\end{eqnarray*}
+
+\end{document}
+
+
From 4c61ce2fb45ec3e4229bd7252e5386e594c2e1d1 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 17:55:39 +0000
Subject: [PATCH 10/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7478
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/pair_lj_smooth_linear.cpp | 353 ++++++++++++++++++++++++++++++++++
src/pair_lj_smooth_linear.h | 54 ++++++
2 files changed, 407 insertions(+)
create mode 100644 src/pair_lj_smooth_linear.cpp
create mode 100644 src/pair_lj_smooth_linear.h
diff --git a/src/pair_lj_smooth_linear.cpp b/src/pair_lj_smooth_linear.cpp
new file mode 100644
index 0000000000..f752b44aa0
--- /dev/null
+++ b/src/pair_lj_smooth_linear.cpp
@@ -0,0 +1,353 @@
+/* ----------------------------------------------------------------------
+ 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: Jonathan Zimmerman (Sandia)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "pair_lj_smooth_linear.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+PairLJSmoothLinear::PairLJSmoothLinear(LAMMPS *lmp) : Pair(lmp) {}
+
+/* ---------------------------------------------------------------------- */
+
+PairLJSmoothLinear::~PairLJSmoothLinear()
+{
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+ memory->destroy(cut);
+ memory->destroy(epsilon);
+ memory->destroy(sigma);
+ memory->destroy(ljcut);
+ memory->destroy(dljcut);
+ memory->destroy(lj1);
+ memory->destroy(lj2);
+ memory->destroy(lj3);
+ memory->destroy(lj4);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairLJSmoothLinear::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+ double rsq,r2inv,r6inv,forcelj,factor_lj;
+ double r,rinv;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ int nall = nlocal + atom->nghost;
+ double *special_lj = force->special_lj;
+ int newton_pair = force->newton_pair;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+
+ if (j < nall) factor_lj = 1.0;
+ else {
+ factor_lj = special_lj[j/nall];
+ j %= nall;
+ }
+
+ 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]) {
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ rinv = sqrt(r2inv);
+ forcelj = r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
+ forcelj = rinv*forcelj - dljcut[itype][jtype];
+
+ fpair = factor_lj*forcelj*rinv;
+
+ f[i][0] += delx*fpair;
+ f[i][1] += dely*fpair;
+ f[i][2] += delz*fpair;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= delx*fpair;
+ f[j][1] -= dely*fpair;
+ f[j][2] -= delz*fpair;
+ }
+
+ if (eflag) {
+ r = sqrt(rsq);
+ evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+ evdwl = evdwl - ljcut[itype][jtype]
+ + (r-cut[itype][jtype])*dljcut[itype][jtype];
+ }
+
+ if (evflag) ev_tally(i,j,nlocal,newton_pair,
+ evdwl,0.0,fpair,delx,dely,delz);
+ }
+ }
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairLJSmoothLinear::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(epsilon,n+1,n+1,"pair:epsilon");
+ memory->create(sigma,n+1,n+1,"pair:sigma");
+ memory->create(ljcut,n+1,n+1,"pair:ljcut");
+ memory->create(dljcut,n+1,n+1,"pair:dljcut");
+ memory->create(lj1,n+1,n+1,"pair:lj1");
+ memory->create(lj2,n+1,n+1,"pair:lj2");
+ memory->create(lj3,n+1,n+1,"pair:lj3");
+ memory->create(lj4,n+1,n+1,"pair:lj4");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairLJSmoothLinear::settings(int narg, char **arg)
+{
+ if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+
+ cut_global = atof(arg[0]);
+
+ // 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[i][j] = cut_global;
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairLJSmoothLinear::coeff(int narg, char **arg)
+{
+ if (narg != 4 && narg != 5)
+ 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 epsilon_one = atof(arg[2]);
+ double sigma_one = atof(arg[3]);
+
+ double cut_one = cut_global;
+ if (narg == 5) {
+ cut_one = atof(arg[4]);
+ }
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ epsilon[i][j] = epsilon_one;
+ sigma[i][j] = sigma_one;
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairLJSmoothLinear::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) {
+ epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
+ sigma[i][i],sigma[j][j]);
+ sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
+ cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
+ }
+
+ lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+ lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+ lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+ double cut6inv = pow(cut[i][j],-6.0);
+ double cutinv = 1.0/cut[i][j];
+ ljcut[i][j] = cut6inv*(lj3[i][j]*cut6inv-lj4[i][j]);
+ dljcut[i][j] = cutinv*cut6inv*(lj1[i][j]*cut6inv-lj2[i][j]);
+
+ lj1[j][i] = lj1[i][j];
+ lj2[j][i] = lj2[i][j];
+ lj3[j][i] = lj3[i][j];
+ lj4[j][i] = lj4[i][j];
+ ljcut[j][i] = ljcut[i][j];
+ dljcut[j][i] = dljcut[i][j];
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJSmoothLinear::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(&epsilon[i][j],sizeof(double),1,fp);
+ fwrite(&sigma[i][j],sizeof(double),1,fp);
+ fwrite(&cut[i][j],sizeof(double),1,fp);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJSmoothLinear::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(&epsilon[i][j],sizeof(double),1,fp);
+ fread(&sigma[i][j],sizeof(double),1,fp);
+ fread(&cut[i][j],sizeof(double),1,fp);
+ }
+ MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairLJSmoothLinear::write_restart_settings(FILE *fp)
+{
+ fwrite(&cut_global,sizeof(double),1,fp);
+ fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairLJSmoothLinear::read_restart_settings(FILE *fp)
+{
+ int me = comm->me;
+ if (me == 0) {
+ fread(&cut_global,sizeof(double),1,fp);
+ fread(&mix_flag,sizeof(int),1,fp);
+ }
+ MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairLJSmoothLinear::single(int i, int j, int itype, int jtype,
+ double rsq,
+ double factor_coul, double factor_lj,
+ double &fforce)
+{
+ double r2inv,r6inv,forcelj,philj,r,rinv;
+
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ rinv = sqrt(r2inv);
+ r = sqrt(rsq);
+ forcelj = r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
+ forcelj = rinv*forcelj - dljcut[itype][jtype];
+ fforce = factor_lj*forcelj*rinv;
+
+ philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+ philj = philj - ljcut[itype][jtype]
+ + (r-cut[itype][jtype])*dljcut[itype][jtype];
+
+ return factor_lj*philj;
+}
diff --git a/src/pair_lj_smooth_linear.h b/src/pair_lj_smooth_linear.h
new file mode 100644
index 0000000000..4f6a5a04fd
--- /dev/null
+++ b/src/pair_lj_smooth_linear.h
@@ -0,0 +1,54 @@
+/* ----------------------------------------------------------------------
+ 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(lj/smooth/linear,PairLJSmoothLinear)
+
+#else
+
+#ifndef PAIR_LJ_SMOOTH_LINEAR_H
+#define PAIR_LJ_SMOOTH_LINEAR_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairLJSmoothLinear : public Pair {
+ public:
+ PairLJSmoothLinear(class LAMMPS *);
+ ~PairLJSmoothLinear();
+ void compute(int, int);
+ void settings(int, char **);
+ void coeff(int, char **);
+ double init_one(int, int);
+ void write_restart(FILE *);
+ void read_restart(FILE *);
+ void write_restart_settings(FILE *);
+ void read_restart_settings(FILE *);
+ double single(int, int, int, int, double, double, double, double &);
+
+ private:
+ double cut_global;
+ double **cut;
+ double **epsilon,**sigma;
+ double **ljcut,**dljcut;
+ double **lj1,**lj2,**lj3,**lj4;
+
+ void allocate();
+};
+
+}
+
+#endif
+#endif
From 8de7c3e24181c2275004ebd461073066c689b3a4 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 17:58:18 +0000
Subject: [PATCH 11/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7480
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/version.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/src/version.h b/src/version.h
index 78003d8d4f..f995c78690 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "9 Jan 2012"
+#define LAMMPS_VERSION "10 Jan 2012"
From f9d24c40127ecef70977c37adb584ecc564a8fe6 Mon Sep 17 00:00:00 2001
From: sjplimp
Date: Mon, 9 Jan 2012 23:18:10 +0000
Subject: [PATCH 12/23] git-svn-id:
svn://svn.icms.temple.edu/lammps-ro/trunk@7487
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
doc/Manual.pdf | Bin 4605553 -> 6305519 bytes
doc/Section_packages.html | 12 ++++++------
doc/Section_packages.txt | 12 ++++++------
3 files changed, 12 insertions(+), 12 deletions(-)
diff --git a/doc/Manual.pdf b/doc/Manual.pdf
index fba5069280233fce2927d42cd410ae6257c2609e..7ff854c4b184860a7164cf98029b77989e386d21 100644
GIT binary patch
literal 6305519
zcmY(ob9mfq_dOiDVZ#Q~*c02fjfNB3PQ%7F8mqCLGw%O>Ho}TY{pLeeL2X#_?91HjhMf`~u^F3~T^V14nm&vaN}etAPUu@R#Z2NDC0OvH~dm!8!tz
zK#m{>7mzVM6;Rj#WZ-0OYa?Rd1fmq-Vq#!?e=~40vN16-vd}QFGcZsARm`2NK$J3q
zva$-w0M&P@BS6-`#@WD%3MlC8WM=y(EhjrJAkfOd+S<;M-qFCu*xZ2L)Yb(+2T=ZJ
z0+9PBih&hC#=y|l!NAGZ!5s9S)cL0l!27bXxucVVxuLTY$Qa;kV+?WtIGKR};&Q40agYtj;XS5;v!Ru_5kSV=
z2xQ|3qNn2H1KAkAD~rJNpTa8cb|9drjghUfxs54M#N5OL?jX
z9v}x>psfuE=;Ue(bTV@Qfq*8q&JI8ma~BZM(cBH_2y%Jnfc`iPv@w4N8rfRe+Pp7Y
zoBthIfgBxyAp3XUfTn+32RXcF2Xg$QFwnpWCkq0vb9P7+HeetKVJx)VBI#v5~XG
zAA{Y2?(eD@I@nr*Yz*HWbOagwBZ#s2yZ?^ne;oPu5d+#f7@NF{`H!}y&gNEDAZy#d
zSGTe?H8(P_vaxjnzVE-==m0W(FUt25jDgk$Mh@=+8-w08e=mV|lm9~AZFF@0$6<%R
zcK$arFnTX#ptbX#DlvH1$0R{h27by6jID!s;
zq!bhd3jbe9qJPPvf61c%Ll*syCHWU4`4=PkKN!jXFrrRoK)Ju+LykKF&s@b2~>GUzWE^grZx>;A1Q^S>DLznFiYk$-eH
z{|^Itm%-+5INQH)+rMzz|2fk3U%-DXV{;et_XKVKis1a0==_)H{69qJ{}}FnH~%zq
z8=(9DE$jZ*D4G9jlrkV2QztV<03#DSC-Cot3Bbw3$-?ySj2Xbn!U7aCxB63VVpi{c
zMg;WURzd%612~yE85jY7UVMD-jT>ZOjo=2J(dw9BY+cK%(XhBntM0tmV7>hBdvz9+
zuAAc1`P90eX?r%xU(2Na{rV(iL0%+-z_gsFIW9rrOQ_h-h4BfInen;QV85jA7ywHH
zqk}C{BU2jLJTt>RhzG=)?wB}gZ4T}9XCMty$&b*y_@GR>ecvZ|&NLj3H|W3*8E-8#
z>8)T;I-_7in(0-0C
z^{lPH>+4%t89t+Y`{7p~iYYMDz0|$7HjXLqS>!Q`))xbm{rJ=cA&h{^T+M3#tUc8U
z@n)V>oSBz8)M=4pYGe$dkiyNq>^L=<`pTO6YLx5DR8Y@qU-&d4No{CJYixQg&nCyf
z@CJpb&5?$~xypI
z?BSVWvUB-refN##%i7Paw~z8~{7;e#X-nW7uccv~r%(KeUb_PYiRzY>
zzs*XAy-i_moZ~Em=Yzi?)ztfKzumD{Raf0T_=Y~{eYYTg8zuj`^%gGKy|{560>6-^
znf~0q@a9A)V7f_av1wySYrAP>{0dZkA=2raTj-r!>Yfu>8z0)r@`U~x?DyIk@M^8X
z*2qv(!)ovFn`Yk^I&A^(Z67lSqnfpbxuJ>`l9KZcv(7K?veVBOox#E0>6zQ%E$(Q2
z;q6Ph*0bznZ|`a^T>nhS%e2DRII@($o1cEjVPDesCE(yh1#S=!Ah3&l9?e@RW4Q$A
zTDZ(J9Nfo93=L})KQKx$ZGeJkb7v4a1(db2?g!BDa+Z!u7$I&B^W%759k%V3Dt*8i
zZZQL!k@u7LtNcgV?zr$y-XCJ2+I*QwCsTvmek97a2RhFLerSh
zys!avWKjiwluph#X-}{rgb>C~#Lx)#mVG?!RPxPBnc_NvzeD~8;~T%|EuW=xF~H6G
zC{@~$Sh@j4S$aGAe1?Yl3EK-3{p;Az{v0Bx%&W579M3?9&dfD1l8en%LNEukFxU
zUi)q5ds5>xbIY5@QU_YeJoriTgzk?%1=E(aybfC|rjdTh<$!irxngbpizGjCPG
zx#0@HjF5`9Xd8Ea+AfTjm?@Mng3f}bb5D&mKWmTt+$ZDK<|YaG&=XwmBSG`zj&=(L
zg1B8`q$xM#@Kwo6H{4p&O#*aa%4(UUC*>LDN~;Txp-ZD#JJMK$!5=rNP}7X%09jOd
z=S-u8Np3B@G)a2hIH@myPZL!hqF=w7b@&(QdQGCjW~IZB2R$eZQ}RjwAVRDY
zA{&c5#5%o}ZeH_~ldHHLw*&MStZ>5SFqA47EkF$nxxPiG41VzHlXQ+sUaVya>cM2Q
zt1|Peo_Jj^JJZ?izPGA0(Ylo%PkK%P|EBrMs;-6^8grybVP;S@qp~a-XGEp&!2IlN
zTs%Q3hC0!nK-Z_t{N$134*JE+5ane@F}o0*848*@Nrp`?WWL$n7bp0P$vKToe67%*
zAMCN583lCC{bi4`bQN{wooeTW_}iYA#(W~Y+NPulXX*gSQ?SNPb_CcLV)Kz;#WMW-
zW4iXDM+J6B?L6BS14Qw7A5bqpl4L2Ie+fS>_~E=Ul4eKM?vJ)nRIE=b9mRwhL|yZV
z1P)#EGPtBD02c)N@hDH88P1;@ImjS7>DpV|9XXX%`lY^7;f`e)eL$1gA&_|W4m)uN
zmldgR8O3$k^!W!=j&h=%N0vLk@LpF4{{gQuk|Tq?E7f;mJcvDwfdf
zQMc6^Ny76S%h+*pj-{7?HCr*@ddoqk4QKgmy^iFqLQ61@Lq2fV+NEu4tuGn@ZjVf|
z@7(jK5<{V(>?1)j_YQ$u#XXvI#+iC@zis+?y&9|C5_?MXcAI*Ox#0|IJ(#jfFxdY<@M$9{c5X@6jp*te5Hho26lJoB3*(y;7db>q>-zikkiJhM
z)y3LpQOi25{gb#=EL{o>ODlwd_fR>YN37yFLs(Zpv&|Sa_2Sy(&fZ2gN;N=(6Kv6z
zX^
z`d{BC&PHf^+{%y6Zc3t%#_szaOBxkKvK*H4RWfNE>c8vLRf~{E?8$QI8|5<3yJuP_
zdm{RXI>qduL}sxl;){?hBY&e=WJb1GsRi#y0Ck5ga4N8NF61M&xB(F|1$E9oO9vkZ
zvJp6Q1`g{NCrj$Bc_Yim&ll$$h&zVl?E&|^ksA$iOUDL)H@~*6w;y78s
z9EAav;f<^wQ{*1j$RzZRHyKGzUcFXIDl2`lT<?7)-HLC&dPJ#9`m;^sC+22Ll|D
z&h5P+u?PGEb9!ZJHGNYra;~}q?MJoTlUs1U*5wVSX|^mY(g!{9Ui+(t7s3lE`8T!6
zVp~lO7(NIst8oGL*A<(b)_x{KoeY;@-PrUOlH@g>Pnbx^eLP!XHr(y5V5r9uKLCer
z-!JxKkpW-
z=372_j*dn&vX?Jxa<{)n;TvF(9=DJ75fu2QiH0#mct-iHnu_jv6C91iA=}Jv@15G5
zACtVT{qF}YVH-B$X$KPWrzxW`)>9lmt9$8*8Xm;WU*ELvI7czoE)*T9gOHYeh6D^*
z2Nejyf>MTF?$62?ka+?9)8Ir9f9!uUrz0c@P#Dy8>kZB$w@w
zcxU8ixQ~qu(4s0e?O1PE0#Bs{Qr&kCMlH%$??S!CK7@{LDWy+MxRKd^?E+U%qEus<
z56P_hGD44aiiopsi@fV+d${56gH(pG|<%);9Mzhc-^K
z^s#<1EpN_!Z5h1(1W619(Rf+dxzG)ta2qnI^{c0vWGJm{E}i#6JleKzv0_F$`8?jC
zme>uiW#GeooW?CKwuDG2#I2nMmyig
zsx8VrhAXf!fmUO-@MoI=H6pss7`5?u{wavRA}n+=9esSNhd763gWMw00T7$^kxw;_
zUT|BvcDfYvp2ni4Kr&yNir;(Os>OsGJ>__dO61NfoMyH}mc_%7V!mv!-=dQ_w&K_$
z1Lub-BZ^~4EiS)AY|g0=QmD_q4@-Q%Oak)_^pw%8vD0fGL&XS3vuzl9yMJAmi0uw4
z6D2YbK`}vYXT{%QmA-m%nFMmj3p?2rOsOAjko`AV>AHD>Qj{3jh-vaOmI>|o{^VH6
zibYgKv+h-Wttok39*T~P%+!96=A7W5OEe+Vsc
zYRSIouszJedro)`<7Pp;kXTu*yYGHNS+mqhZAofz9od?h7J7w8TI;goM@kyzyi&ML
z7FK?^k(;(e$R9N(jHTaDpb}?T^jmaeS=T%8PbPBTJ9&92Y#}d%jwL__jSbQ56xz1?
zoP}2RqxJcGu{$(FDI6!by}g>$JgIly7k+3r<+=k!XH`;`Om0R|9kxZ4hV$Kj
zSQ~BJLrN(Q8m`NS*B-U6*77=#UY3)~KcXjlj6i)}vfiSDRX1Wkj~(1Qa6d|fA)*a>
zR%W2KjVAztLM7IurW{FnwKLX_;-z*eJ9i4~Ef4yCey9%g`$4X+LSp351!wQ)M;ZvK
z{H60-d$fe6w6@lYfil+{k}QtGd~(AKy=XafiCzFBQ~9WBx}&>bNd+vNPFjftudNs&
zQIN0{OgM#}ye+`9i7s^E=y(+1hC7ip4NJS=Um%fNZU6cTIfd
zurfA36dVU`^b9O%YXgLMdB>N3Fv)e8eto5gF4E?NZS1E7boHM<-WNUj6Lhnr$V7tY
zjwzk~RGD`6o6Dd)T^Vf^>N_7q7ftuaK@&KEZ6bkb;9TX{lFIiC`}x_m!7vf1F<1
z&{eM`LFX3sErsT0`0M&iR$^=Fm9bdIMl0S0%?MndF}yOTC$nH=fSQBWu|$?v47apYBodv)&)
zw*rCla)&~yvN5niT&Buzx_!kkE)DL;F4-D7P`mtGo2%|B`)l(cwl8b-iOE&1{Z&}8
zm+-$k%EM2a&v5k#3M#gVf`b^Ga3{2jYl=Q!i~yia7s9BD=W~|=inRL9g}z@Ip*FzR
z<91=96ipvysKtuBK&zr&_4FhvF$dWjqC3MVwh;>nFd^xu!!N($NLzm-
z>T7gUv5nhgXn{7HS(1X1K2Il`319ZWc&SmD%`?E?dhQmYid7)M@J+%+?~OkSgV6JW
z=4E0DxQg$^b9^uXs$nh%B%;dlP(|%?d?2~UxBDLc(-H_>La53N8H^2eO*AJ}@7_ap*83TdE#Eg6nu85<5i
znrks?yAAX8o0i>O%pO_E=Xa1h*cCK5dPow%QnZ`Mn+nUPM2F?6@2H0glm&2vAYg8!
z_;m|Qv9!eOj#c6+l%1#>2`*vEAeLoe@Q>5EB!BJS>nn%=c&DSIpIe3GLhY&QRoaSay>=6~W<1Hc*~We1AxBwmvQ=NdwQ?e+woBCF`mUlu
z*YuQ@mRVTKXXEfkv+l`eA`(Vg4J9y=!sg!i_cEw4N%e3IDRij!MBx97P_2~g+XK3W
z7sH%KAWfoc2`+7^O$feiUMvkwMRIWLet3RLaW`K>(e21XXrq1h#M3Nf`*BM{T^DGM
z`1Q~Z0sLeZ&r~3xp9aUYL!>DsklU*VBb-dPxIUV{pn!D59!&g6?EHf{QAoDu(n{kT
z`?r?e%ts3!X8;RSp%zJ%fnuiMWaXJrDEDC!R_NR)e6#eFmZUG+TUXBwC!+O6YcF#1
zWF)ICJn_-4HXBUQmjUTwU6wdeL9>Y46y|kOZZ~f2U%Kn~jhA_@;~6LwBIAas^GL}2
zM_PpsP+eVMYnc{1XE9MlShbapIsq@B6~D>#h}lR?EEaYxzSCVxlqZrTw_yX3&V(Ju
zO;glvYGx8yG`l0#Tz$^Rea%pHW5B~8dD3tS=Sr06Sgn%k`7r82zEWc
zi!XlTN>IyS2P|1}&QDx3L*NAr7u2k)I(HRD9BUBvUrtRa17__fZ+Cu3-p|y+ZX5nTZ@8TNoi*23Rt>GQbBzj)DPj&+ESPkmgva
zUU05h`FW_Qpylgr^r^LiAD_sRo;NKE?1Zja7fFZxtf*RcKRB7Kp?80}%k|m%R=^J_
zU7NGpl}sdJ7-^T%6!=-;WSUissX}zNmVKrB84SUIEVO_MxeDJKPQLAs)o>m^&N}U6
zSX)T?SGh*=w~feMnN}a+RvGb17M7RWtdR4vpxJ|Aww})e=f^|Lhq!*3-*Q<@GBo^4v(-1O
zwD>Z0&T$huJVnhTXX%Jk)F*2{1{NFF?K$7D#^^!74Hfbcd6hdB}fg!NBKdNb3I4xB=$wqYMAYh(y2%X=Td+n
zR&8bNavMLNWm&Wj)of=Ft(tvTZpBp)&Z_7;Kt4Ta@cz8N@69mkJH`E+mDF^`TPhJU
zUBo=ypOR|MzB-}rdU)*qw@Z86I;-n+L@@;r}heCD$vhoJp_
z3rXZtW+f=HluiERr=w70bU@juKszm##Kl&?XiKKwuZ630ZDZj%6+2ifi~Yu$?;Eo)
zhzvSS1^XhWbju>GvJbtfDU}6qD4EG)5RiFsS~-;0;BH~#+3+lr-dM_EzV+H0WfBJL
zR1&i<##!Ro^~J_=ZJ~ZK7G(`_)pDb=65pOoe~OSXdOe=1QG^8@oX^N)A|FuDAZe(g
zeZ=9nSS2F(DOm;A2-yx7k%4>P0f#^`H$+w>;1~syx>cKF@GC@Bk!;JCBfTPf9GOet
z;7p>Xl5`zzwhk$wl-JxnULhj58u}SQwz%TEOT}KZ8{IF4BWvhS-$$`Af^9uF8{^G?
zwxTr3*Hw>ZA*Zz?@W+7N_;QqK4LoL%9#Dm}Ys@dt4jUPUmixw&8?EkX!;sa?f>7+Y@M{QV#O-Yz=!+o6Ae}bd|
zu8)CHs9`D;$)I82qK@?^)D{1JDR^2%2
z54fYDa?vvB8q}gJTrld8YrsXngYm4?1J73c6VDq9c|ymrb=bJMdclqAnJjRtYCSBC
z$=a2wsI)t!Ffa&r!L4NuPaPhSBz7z9Qg%lRZ!BGoUw%wwvVw#+ClE&a%=`?|xNgMq
z-l(r#Kvr?zjY>w;HD>qkZf|D~+ZTJ?Pt+pyyEc|*$Ok`CwIcwJaaz~Ai8{mlV?GiQ
zIX;rACVSAA1Bj~vAHqjestq5%aJYIUW_GiDVq{Wpl99qoCKai#ES7{QPP@HvrmptL
zoR!(&-QPHPwK>uL@l>l}{{qKpKY`c({jdd=T-6IhV@aiQ8mSwG{WmB6B1Su`at$Fz
z^Aa8tCB~-vq>b&%GnRJ%z0nilULK$Qv76Cg#-PwPaHBT8+!cQ*ElHYzycaY)ygr_I
zB8Zypr|_c(&xD@pqrbw1EItwt!xFT2XY}3gwn-(WyTaLiyaY43Nqgf*6!N|VI7Kgk
zTKv#JtjKfBKIm5Oa5r-f7~qed?gH#*(x0l?4?=N5wDqXp6rg*KFEgk^mcRe-rNv|4
z?-t5WyVhhzm32%i2@gb?I_^`yGmB66{YPmAv#w5-ws#ee$?Fl$sM$aD?LK7Fzv6n6u75w_C>equPtrXD6L;d>^*d7`e&RAw8I{y=Ho7+YL
z+6!tKwSx2!=sv1TLbNt!tfFSkM*$)obonK|YXmyj8)^1j63`&|Wwr0?@Rr_c
z9Ba0{YY^3AoqpXQG|9E1ZR6DE{*f7bs~sPbvK+`jBWB1~DMw}N3EBY#_xRXk>+Q&`
zC;k{LpvNf@XT~Dpt6C5&owVXUu)Um(Yo@8m+LP0lp0CciMJ
zh7UVQ>9iQ(pw5khl{29cu+lusqMO{)cmqU-PV(qr=Er%G+e-o7XJfR2i10WV8dQN#
zzgR1{_Vyk6LrZi&j(~4<3c@TcR*`&9c%Ty1NPe1bN9SS!RgbZR&NQ_pIj&Hvzl92?UM!Q{@!Zg~2uIh?)O81LNFw
zqt^@$RNEB%kj=T_B#5nc_H8J6alky)y|1>1)nIQcmxZHSCmvxN{O?d
zH{&{CIL(JX38=b{yf{F7C2G0~c~lOvq3&XQTM>*8wBjFGz9ilhZ<(!fbAEb-`i
z6_}%*_(a8pWV*u|=ylgRRly~Of7>xuB;Bp9I?YznB
zvUf;k-1LK^(?I4?80w%Wvd`lSpNlk0MVfRz)8>Jk9$Qz!&+j5H`XD#P;|mXOyT>Is
znb#4&`|N|XPF27+VhjzS$PVXMply29Qqq-rAVydA$qTdHnU>iilDGarzW}25#l1-Y
zkqBHb7@JtnV)%Rz%qVnDRWdp0tC$Q0hD{zN1cf~>0&KuV*z7r+QXgs2Y4^=$xCz_R
zcD_!Es8OwVOZaq}e6^i$VoFE6M$>egxzi^4cWl`n@imk8zp0ihx3MEVBHL-G;bR5-
zHPk6uT_PRhWDbv}mY0Wlea|&Mdn0#TP6zvEs9E^5MG4ayz7gJ$5TFDB{ds0>o<17-
z@H46EHGA96L@U`g<1@B&f}0{J6Tw(U*^&m!01R#NW}m*KXhl({k}^oG^w>jo-U%IO{y>c525vjT4c}S_xdKm`Q@4yIrDRdZ~#pBbf!G^M{$2}dO}BT#UWQ8BA@3@DD}KA
z7~7g=uAO_zFNq@@le=HG@FP0AJV)4>B5^)hRUtrW`@*?g{iX=g7Wq(Ao+u1Xn+q$*
z&^|eo)T;=(7l$KQ$jAEx1IZu+M5tme%YPrIUo(ap)A>0uqzMHu@7XDh0h^Ka(&G_!
z3r=g;u_e^*X>c(=wLhXJfkfUYi6^Vb~&p2@`dL|4u3EQ^+$v;_SW_c!fF*N>%jS
zmIndCKsI`oGmv)&5CBRC^L+TNpThv*PW(g7$oF>omk7uUV=ap(Rs#
z?hP-szMWjE@|8S;&dRg<&RY3Kqzj86!efg|%*iG!!Qj@}y%qMVS#D%#hIdzF+{G(R
zq%Jv|C_26lZhqX-Y#;$$+#Mf3)F3F~p}$^3YVE-gI`q}tIT8&pppWxbfD)OWzfzmR
z;Vfs`b1{)^*D%^qVAN7`D#z8|8F0
z=Qvk#Kt>n}Ql2}4er~-*<3~vB&2qm*(JxJ`Mcgn;SBi990>9PSN=#$`_@Nj@Jp>ER
z^CsDb5U!%^)N;vA9MoxN08`ILfc;rAM2VtFTv7p|i2Jv_+{Z6`nRPjw=P3OtC&IAy
zOOcrjSOUh^oU(+?7>@WD4sZ^zc$KoSsgIH210@7OpJ`;^YC7igXQh2VcClSIx*t5ep!LK>rv)By;Z)Y^
zZHRnC-CH(`JW0>R2CR>9B7Uq1!zs|LzZ{fuUUowSR^-Vm${l@y{!losi
z)X}0(f+IYDkn5xbYE*~gT=@1o55{l5#hm8eT(qEmOm;W3yT?S<%J2HAww3ZbPczR5
zz45GZ4wY_RlnFih!K3pZgX{G^CL1*O-x*|kX0Tgc?{J^8->ZrbC&vU9RhjnyZ7NQM
zgm?h4-;>Rbd$k)uiBpl#p9RadJ6T%xH+g#iZKcG1dZ|xEzB{mcq@TZQEJ&yWuYN#P
zzrfY)4#D!*lF4Wr5i;Gu;(9Z}f=jQi38EsEZ4>er@ajGZ>qLL$Ym@(6x61XV-3aBn
zl*<%UgHJm?c?&Jyh;M!X;csvVQUT;-dh)K!?lr
z@T6iXXMhq#g*sJ-=_>u(Hagb-%J2!hr);e^5w~Ajy+UIQ561hOTT0D#l_%1R9%fpX
zTc2XHmCbCp!i&FNTLxk~9=%6Yr*K{MvS<$8^0qrMo-1+`CUE9uM`JXn{pMafw&wX_
zuaKQA{cf|=%795;;tX_M=ZnD#lQAJ6-by;|BW}e@w(w-@DV?JfA!f3eJb#Qhj8Ah%
zg&&33G(Ku#Q@3=FYAuYOIHST+=Py4rO@hoH_*-D
z`z51-z&ts-zp
zTPM{p?1W*khEFdO^OXnt*=VG5MPi#B(3MJX=Ypc}gE42@IOy?wQW1kn6A+Ak!rDkZ
zi60z@ZH40>bzvNtIW)H{1p6TDZcg>!I&M@ih(C~&*{oOuXEy7j#+KirWX7)v>2v*U
zIezrIrlfeGm1u$U*R*oJwT5^Nsm-QFVfEobN9c!bjX-QZiH~Gdq(*kONzg<)8gg)Z
z2DSO<5O00*X(^7N~Eat)up$CTQlN3Fix
zC`_MM`fZ$7VYG5jL9fs0MrAK1R9G(YELXiYcovYk^}*zNGL3Kjp`6Jq)m?^NYE&y6#>eSmvEOoL?cr%%&`pawkK
z{7L;{oWR^eQS10m3YIy@li((5iPf!!hN;{7DZ929*9e|jVrSm?G2Yobw{8#55ZQwC
zy@sSQb>+UQCgiCr0XYk}pDt|>+qRX_@L9@&h7lCe+k)X+nm)0c{E}tc{b@RTJ$awC
z$fzb*FMfL06Nu$nZ>wYXqFHG()s-DR@K8-8J8w)iFwA^8(dDQ@J8?gibKe5`)n?N6
z^Kim}R>>$n7LN1>Cn6-{*m-6Gn#66RaWvuavM@;!*+hhraX57b1EY`oA
zM@++*%80vNNxTq)0i%B71nnyfSlIwj`18@`Jt~d`6I@i&Y8;|5L^H-iF$?X`
ze1n+mr*>l86Fhz@4rbroL%zUh`9VK%+t*5P$b?s^+keM+_C2dRAjEVg^VD%CtXsmJ
zFzLr*m!DaRv~W}kfA=^~3f(4J%&@|Vuaae!6getcjWE3LW0iTL|HGigy}-U=(5o5o
zD;}>L!NYAXI>YnF^XPwX*e})@sCH&3!W^g8xaOFihZrR~OIpkG&
zm-m5J3#-asvCZ;B0t`2_*6T1^3y
zVVeQW6IsR|2AU$ztJDWV#02mxo(PQ4>IdEOTPd@mf>T-D6sF|23l|N2@;kDtyw1sG
z?g|yxcOWOOY{XGFsC{OhD;+nr_L07P1xmMoUdOUQY(0w=)&!WJ--ynaPPXGmaU=Ip7*515j=;%{C$SihduS^e<`Xqh
z%q5T1azUqU;VZ{86(hZa5*+gET4uCiOghT?T~XCa$5PVB;xLGXmAf?BIvm5_h=*Z7
zC4vU;!Bx>JRhG2j9{Fk|TCi{ft{yPGa$HlAec7)y*tLG{KCA%W7#w4V9gcq4GS{9w
zU8d^K!b4YH$H#ee0EOy%L7?IWEHH+$W
zkexndS}1~A(D=QnhVJm^3?GY&R?C&A=H(AI^7!a*7#fR!#FPkpzz-nZasFxH6avySQm-Mt!yh&gmcDpw26tL7b9>fl*8t75Q;bfc+-x-^qNT$
zZIj>!bouq3!V+clPp3PvAJ4wlwa3Xkkrd>b7bNrLQ2FHknst)?seSK~pv=)!KBA%n
z*%NWb9PjEv>sG}o*ON8+f^-;AkY80v*YTXT@Kd0Wn|?C=x)5`#k_R`F~SOw6u%7#gA)=$q;I6Efk(HiD#+y80K;!QEotr`*YHrhx&dc7wQ_=oxvqJ+4_v8dVP4#cSBk-EG$9
zbhnanNUWkQaTmGDY#d$j08Q*OHE)B7-Qov>{thLf)Hq1!Dg)d)jf5}~hM?Mo?HhQ^
zRPH1jqaD4twLWD*FkuB|KE;tcqnO;XoVw+o$1#aP6*DxXDrhRqK|
z38~Z$KJ)~l%A0d
zPx%lQ_VGGP7QR2?_jw(+$Y9IJU8t6({O6h2Ry#GDW`lu8W$3`d9-dq;x?e~36=3hPubG>lP`J#wqx-JkV2b-m6Y1-w!UKG?D7^FVHcTol#=SXHci_D7@<(7nV}eARX)2d
z2AC}k!}q1qS%NX#aGQw~{w1RdiXPzxLOM{6)*{}|HGM2Y?)sXEn)`;lkg-3Z?WLMk
zI7=|85ahz7M(yhxX3-HJZGH_}L?SiW9D2tRgxVs+M|fSc0>Q}6(cm|cRrsqh2A*lL
z`za#xgT!Zpo65i6sQi|@PQGf#5Fh7qsoWM78G!dWOlQ?-X3}b`&2df9-1L*;z#v=I
zH1NsWI@Go+_N(s2uFL=OWBFvCC2gPQ(ED2j61zwC7qrupY1MD6*hE&wo@oa}+)t2n-dk{`)TOI?26
zNXvw^=5FXs_PbA_$Fq}y{|Kckkk)1ec6M*LRkvy
zHRiM0j)Wy^B=f*VI`~-Y-xw6hcV1J?vmthr5Jz+q>JzMT!qtmVpp9QtebjmO4VgU(
zg}q>k-Qkn|#__MDmSeZnP%ubWup#Vl=E#mPPGY%$jELf3h+Pr0dv1mpKN(oJ1my2q
zPBCeITV47yHO0E}Ck)JrXAx}JSf9l|t>#n4i(Oe#C?ajLUszU@
zy657$33pEpTrFvu)YI?}z^S~!zpnVXNf1swkJ~lg-0E;!DNavs-$;m+>@ojnvA<(9
zNAsoI>%8)R-d_e@JlBp2%<+}JSskE|j$<JGar*XtQoEx4$ib)A+ft7G+V<
zkbZVb>pP7%_~|TyOGt2S@WmCHi3{@frX@t(0MO(
zbpXN>feiFZ4WHcG`H)}}9}fFXeTO?y&t^BKrG-iD1X9aM{Z*4sKPMMVw
z#x8iQoLnq3dy`3|0Z>6ywpFLOC|!)@X79kOOuY55aHOj;TcC-Z)7Rx3}
z)xISNNMVycozE>Ri;LY*=0;*y(G036k)i&OzTK`7`OLOtuH8!QeV9ypB`+MY|6<^E
z!cI1;F2ltynxA)s5d1V;MsVQBuCGD`TmHDBVnl%L=k;_pOp&8JhSO>^e1Oc;(Xj9w
z7yBE>yG~(zu7EjzAFL0+9~1xm&1CcImWbKMP);x+*4b>itFau4>>jJ{6PLcof3wjqE#x0NP3i7&C8-m-jfiMRUAuZIUr$s>m
z-$vf@f|GsI~1LJUYeL8AjP7Dhcp&p&*-=R9q>SvrfJ
z(4kxfRAHSTG^StheFefC<=gMCYH3=?x4EbW%jg@Y9
zlDIAYdmCcg$Ig%li!vViOWbkstw;}z5j3{X=*zB5tSV0D
zWkai9=_x(rShD6Mleo*vQLUPO9AO!^^n4~H58;y^G&j;y4BwA+s~aMfak~SF+6{-DJxbS-q6D
zPEnRehYOX$4|R4g6w%Xv;APEC?K?@ntB~BjplUgMz)Wz(e;;s77mqb*-*4PLCKin^
z)kB8e)u_4K@qhV>Z>y0Sj59?)@|Y0O(|&PuZ~8%`1||5!gxxYsWQyFNtvw&-es#xXg~Q(O+I3@mlYsIsrtHX1
z-dZb2PQg>1nOR|R>`b9VU;T^ys_g{BD#D_)HL4BC*hXH5mRognPspHdxh~EFvnE
zSrvMO5=rtV22ij$kHB}WZ06?heR@-Px&2m;38}9(fA9GFUpi6`jtC@|gK7S)?^QUwhMCZjmVWH<+5&u;=u1Vs4diwTfYf*W&kjQk?yF$*RZ|Q!^^U0SsuBnuCzy(j4wxDh+%_|EsM@JUisCu#>4?pHSZjs~%f3j-F2*BbcpO^GN-
zK9*5=*gkU5D~~)vj?7)xb?A`kchPIUyb@pua}TQYFsjid@rLn!fIq=lXBS$P+Am+T
zMO74(M%VB!?u|9-?ab+R|J|hz;GWZZ~AO
zQ<_j1gvrDblg}vj@t*#}c0{VO=PI3G{Jv7ME690Gt&$d7>7oHhu}8?%^=iI+*rCLaM;ere+!YDZ53X=H04MX-(Gfi({8N8r-o(YXNq|{~afL<|Z
zHlbjk>JPil5!j!p4p@OSbO-275Ewl)eqnL^wjjSfvW`OC{iUx^4|8JimT%td^37igOr+>Xk6N#%-MB`caqX-T<%v_
z*e${YkcDR*e_8PpWMmI}>IV?lgF<;946#+FkBY>h12oPD4qo7*J^fGq7chsEoQ935si}uBHnPr+>RrA%=@7)H8`kz`IdN6IX+i)NkLg~(4NJ3
zwRLBDaty7WQEnw{Ezb3Zk5VF6!n+rzY|O2rEm-DQxhQz>6A0g!zA4@JmC$I^3FLFI
zKq=j9LOafAuTO^B^n|hqQBhbjK~&!5DAHQk*b6x2-GC&m&j8;UOzP8cv?}koA-x;r
z@{eNw7XUjz#J>;Zu~xzSAVGDTyeG?}kF+h;-NZM&Ff&ZHrGsnbr7r#U#{pkfvPQuU
zSZ!^&f4*>ns?c5kxY(Vy-Sw~%L!X2)5~rlz=vsB>MxWE}dEkb!oix9dYgRdy`cvKZ
z{gn5b+b9lP!-YQRoHKB1H;?FAPzv3>c{6CLStPDD)ckMDh2dPk_mEMYB(dV6cy}9h
zjwdl~BZT)u_uEOa>G$-UWYH!WGK0y8Q=OE`E{OF26BYyS8rQp#qsro=?px7|c1Y{3
zlzZ33cJdWVyn#DsoD6~nMD?6%+H)BniYOjq6r*e$V)jv=8jh*9RBB#IdLuW;*bemb
z42s9ElNOHR0zN}e$e4iCsM}f5{JvNr%24f?n!Eb276{<500
z>8Cn6#1B@4hWL!0vZrJGUno*Ta120g36^Y-an%5;c6MED1C`sLWR(rYr_zv16~ygF
zP1}C_S&_~ER0;GcMP8YzBeIIY(F~wd#www+Q#%|xUZ+U|biVb-<8phUnjl!2!Y6*q
zI6ku7W*OP5Yl#LGXj}LkxVtlP?)?jo*?ukl+bb&|#a8fR<4uMNR-=e==?Jdlv&6HH
z5xgus*Uj&Z^LzC{;)?PA_&;ixfaQ8+%n9l#PSLN-{rgSWk{ljTPph-MpS=DdhMeDJNqN8*X4&3*7XcZvKa8r5bzvzXloNHsj);|F1n0Q
zaoX(3LPKm1^N%IvXb|u(AOlaZ9e+}*k@KUN8azV
zSRp4jNcW&}U@`^O#PdixG4kCR=@T`;}q?nL2;Lt`&J1In;
zh0_0X$dp8r_-aH7d6|=V75d3ST03I)9$y%h3c5*<{g)WTP=mA>}VTxK;CGIlS&Ick?
zDIL+oVR&(_2IclCO2t#z9+~ShAz!Yu6gspFgNc42&LnFm8t|
z%9F4bn{;{05xAeskKnJ0h`#kRSz&8;DVCQ=4ORC;4dz_evZQ@iN*hOl^GkEb1K3*(4UjetEXm*9JIU-MX9eR)Em6eABr6GC_-pN!3ETXu8Jt
zXRj>(s;Bb&!)7tl$lsN@4G5XV9~(+Xk2t}
zH+#D0?RaEQrMi>9;IL$P!&u%MMX5d2Y;nR|SYIDz`(V7D+h>4~)>X>Fdact3H#{%j
zK&|2yTXOqLWI)aRyQhPw@1#hKJ5Ziyd*}FOY
zrvt`P*ANjM^%_ws6=qn*{g?U35K}ATPzsbNk$5l~QtGY1cG&{sz2b!1P2L|2u1=
zQ1eX1@Nl9w%b^AG{pp(pVe5n0m{bAz+xTyOTzb`k+kL8qO)nMqt%1q@&0gb>OtIU~
zw>p#HQMe-?9=MG`IVn`fx^FbZ2@!K9@FCH3sIuE!j);-sw+z?nq_x?%d>oDmXeY-i+w2
zGCy@=Z#A2{U-;iH(?(ih-O3Nq@}IzF@$E|f)~s+z8L#=LMzTlW9R${mlUas&CDk)>
zIlj-LR+g7Nxz%UvIzhhR1P@!lNJ&JtvhFb`phKtY|v1RZQ
zZ(Lw912}W7Wt=k)4dy#KQb~`NF_5G(zY^O!Wr5TYM(?yo=yB{;MEWVHzAHit;?mBC
z@f&|YvzL!Otk5;#tEEPKyXC)-jhL_RhdY*qCnfdt
zTRcl`G4wj#X0{Tl(T~)$AVpC9l{SO0pwPU}%pt|?^G7JNuO
z^kejD@yvs1Id|4gQwQ6JLU?*clj2m88;m<^Sy-<8K#s-_Cv|o?fzM{#4d8gjf5%jm
z#^!}{Gv(*VCw~peWd}-))Ztxay{GwM0WXT0bT@V^RgBQCU3&n61>1`Qxi`a+v$!oPAK$3E9Y8Diq
zy=Fqsme}-psj_I0ZB>`t^{Z+y^p#vkSmb&iA%vuu1H6i
z+Qqy2=d|BD+k+U4&iGcp$^RGzFq@3^W=aG}kPC{h3RM7R%ZpBgef8LmV3($3q{Oo#
zg-KVevvlLeGzy4aRgVoL$Eb=K5mlc7Wa*u95YXh1Rsz}#n}$9sDD8(|PV)~0U;}?^
z0(7|flNf%8TH?lffTI4b9cK8>*XK33XFJ&2TTouEuy$?3I3dGwtn8h%pN*D&6puI@
zm5+U7R0m)%`PD`k_Q-;6Ti#I?V87@(a|%X6&Pmy5c^+1@J_PTTgpdXfBmv7mNBm&C
zDd*eznRuSjH!d>W9$s>p!;xz1yPn;4MnM_R62!H9uoR2U?w{pFqmtQQgbgpn{ov+U
zBFYi~zn{%Iys#u|TaeFObEIAh?$)B=|SEzG>mG`mqf=CK!2Ml8BRny|W
z-_O#~@yCzWOK+d&P{+?=75@qSf3(7*gz@2$A-`&iv7CBqL1{Z
zhY%t=aa2rY7Ck9?$!M*7=q%p1b=9ONI3_^vyxY5kwl5A{-@m)DuqR3oB1)^#NeN*;
zq&*1jF1$RC%)+u&n`j(;N=H_=Zhhhx)h{4%U6A7UM1$^eL3_o9y1UMEB;Q;jM5q|y
zBjTdL^`REikD8+EMCZJSSUQ2D&v#PqJ4=?HW-zlSF;(JDp2H`{)Rq+j#-Xu+P`L~8
zw?{W09@?rA#hB|nPrrZJ`aX{2DSK>ENtvSSE
zIt(McPK_m{Nw#T5)u`#W12hBm?(EJ&X5Uf-yOPGGTtW!8zBJRLZ=SBw7*KQlTtNFG
zgfHer^~J7kKA82)dP}(@u3t*6Exj-YyaQ3w6JSg@SCM!Q)GQdSl9<2r6{+?lvasbV6Fooyq7M+9@Uv<%00H0hR`XtP^(R&N3HJ
znpyKhP58yIF`$49pn(v0fDsK=|MA{I2la+;V5jL#!iloM9|b+X|6^r*QID=gjf0}j
zM!JhOl-&Xy7#}{TIJqzJ)qSM&b@#*#QlwF>#lV$eWZ*_OkK7l!p@g*fffP+pksUO3
zR5%+n{)YR6-08H3cfdig5TG&U7uI)irxCjf*Q$H$?fe5Zva$N;+cRp6b&^>YLdZCb
zO{sPB1-$c%K}h;j(t!hY47s`!yIL10psxEFAOa9NVLC0mLgZqwNn+wec5a7U1dCJ@
zTIWg`|qitn4
z$PQd}k9sGO5ZX&HCrs@S&D@n#Kv32a~NA&Ok!+-fh05a1q7)QT05(_
zeJ*E+)I?Ra
ze;V^nm+>wrarx%X7`PWMrY<9m4Z}%ybC|7rWNXPv6i7dXtqTQ@9H!TZdJlR$aepM
zdk~v^BQB+syp|=Lb0kHos@dlOAM
zj?yJ_wnDYYW=h(-D{b9P&+l9HmJ>LRd@Kj2bT~(@Fz*AU(4sK@M9hS*583p}E{q^6
zdhZ$jz>)pr;lx$9zfhomo6X~NZ@Ga_7W2Dt_1(D%zT&p)pO*+)sxwx#Gn#z%mL$E{
zn@Fb;=PsrOTeBwR6X6(CqB{9p`2kUyW8}C8H@KPZ^ebUMp3LI3un;Z;?WJX2E}*Kx
z(kkRJvQ!p%->U0_|BKDtm;o@OmL(ody7GcQ+M7EjM^0>OXWoIDlQRVLAlEly#uxv4
zf`ck35bu+34UcM3izpuBMo`SX?igZk!|^&FudauTnIw`*Dx
zhu@1=rZ<*sheFlyls&d1q!J8iDXzi*ojCd)!1eAhu-`VQW>OyhCLqp7V`V8ZcPP-@fCYPv(o
zkv%8N;9-5X-2(C}QIcgRKWbX|1xXr~!8s1mI|eitQK>LG0B~bJ+~txdEFJK7^+pld
znrP|hwK<6NX(Zv-<<9sePeEAW1xQ~OUVp-q5$t{aTU)%Wj=8bNLW|9!p&qe#t-i9AYDHf
zFhdx|QQW1D*h>3_x9C(}A_{h2^~}ow@oOaES=V{dn7vz??xb7nM(p?5TQ;|*zuR#d
z?N1kCr!w^d*v?kL9VV#^BFa7i2(GM{f_a|Ai>-fdimRk06(lh|fG_U+56L3A1go-=
zMthFlvnS^jEagZ=HpI5$JfTzTaUHj=O~BzbII7VsV7140t=Ch5XLQ)&+ZNp>C>qOXKO&Ybl
zuZQW)3X}3lsjZ1EO)+*#7Pe;1awKKb(1;A>JDrq{>;!n@oYD|gqOe&s29RSW>d&aZ
z=4;opxBATHXuNUO>YQWfM2g8sG4Wndsknfk^|i49WL~|`g-e3>!Hx1jMYTFws(LIU
zwsR=xLdR!DdE~g=a`_Duzd0ofb7ZqAaNS(}6}D%&K0IV7hLbH6$9elb>MzpVZ_m{!
zpZ~xAXa4B9Io>j34+-W%ETBrv(*T-#>Jhqgr=U0gM7T__^Ps;nZ6p25^Hp(&&bjcY4#bYuYHe(|~;cgz|s{~Hr
z`ueRLbj+~Ebzh@mpj-!kGS&q*?0xb*{4MkM*jNE=Y_rMc!;uVr
z&c<#xx2{+xX?g*dRKjhpV6@OzEn!uAFRFY&P!Tdvf(V!rH-~Y_jMuOjh#Ini5hOXK
z+2$^EpAh@crAaN|i?u(+=Pf?DXUGu8?U#Ek^dod+fkKE>#dvSlW3Swt(a0N?!xUI^
zr>2m;zAEZ}8kR8smuP!9Uyi~OD352Y`TRanY8Xmb=Jy9Tc9e{XjxLHBBS-E7S16`
zS+})_`!}g)am`aY{UR=53LPSqGjG70hmR>;fsFPMqH+g_hp+XzTmO;*X&xm_CFrYFoU^zh@IAM9XzKSFQYQka4+7roKkE(i|r&g2obomou7Wr
zahmY|3{Gk;YBb!Hpd~tJkx#9a#K$F*R>2Xf>n=7y
z+NIeaanh!&E>87p&f4~9uI8JW#e=A6J%}UBd=sB%r=&O>l7k!baaC>Aqh6_>;n{&7!%FBS*M&3O`4B(`qW1=$cS5*=H
zzH0gibNHaYyFch(MicyOM&0-f6Kx$&)^*EI`W7CKJu
zN7OFHSTX$K41LXPZ#KLV;-^0X1k2XS2$=bLKZysWSBs9jf8Rm}B4qAUR)9rG?7e6B
z-YfgH0w+Y6DCBcu29%2`hK>MEL4AR4ncq`s7J_+8tKxMR99Z}a_E*Hsi=RV@83nueOhG2GmU2)6a#^Q_|+}?BPTLiMcWam=iBAkiUmSq1SYqI)4V
zz9EnQI2tmMfmlD>FtI@#H)n_3%tSanq&ro7CU4Bw<7dthAB{Lh4(=2rUW_0KAq^Ws
zk#Ab$v(HWHDPD(fZ|Bx}?wje&Me70{+Y-Hl>xSB?e|da$?^><$0G<9Add84t!h|OY
zc;b&;0+SRcI++Nnap4*$ES{e~lA8rwX;4(TN1k;w#q|1JhB7An)&d`94m$ryx^j`l
zpS*yoT-LRPZ+=r(z<;y)H^vp_X8#9-wRp;Zv00F!=^&9xJRC9}NL2{#6@TJVzVee}
zDuW7k+4xiyYWAv(X&2?KccZH^&YXVK?`4o)Y|ua!0?7nFYfz-@#y5E@5aE3+or`bh
z49Ab`V)N`j;zwPW0_DLKhKFZ{&cr25jO5D?vF)+XM;ijR$!G3o)hs-nd-DY1|02^ZpLL6W|;njr34IW;(QqPddhOLtv^N3Y)PUaTy<)J6mveW^jHz
z=E7t_Dw
zp27GXKgROBDo4zAIX_2zHG-v8-eL|;7X>fx1KCaCI&7cxIvNvuE~>wpxq|*fMLF0J
zYaM-5069U7Z|R=Ca)B}D?bFd|;{6vSS1Dr7^5k4+*eih-W4L5pn*D0SZH98pM=qb(
zM#zK>$Fn7vH>!8CoB^z+s)dh-Yc^?o*^puHm{Z(dE@O#hIyfQGgAub$+wceSD@e-`^qNjGA
zssatu{59y&ylKg+Q?a@Np>hY1u(on5Tv)WnptAd?Yu3`t
zYaXY{)QZ1V;=9#V)4elX8-K~jEf!@Y+ruVFgBoe&=KAl-@?Du_%3#Uj+!=f{5e2p{
z-V9`j&IWjJ2OH?JO||TuzgF|XSw=P|0`syJfZ_yATQVqQa7~dls>*)U`C?vwp9jdN
zHMXq$Rb^H)@GZQIg?p2?c+qFNGw^?Z*>&f_
z>ksrRrgqQn82Ug1K3WYdb@lTOt?~s6ai^GZK(zK(FEt&3LD!R-`IZd}MjdvP2C(_h
z2h)5amuu^tS%L#>?tmbzKU%X#!l0b+EHFQ_tISepZUVF6=R6o7x`wJT_yV3EwQ*Gq
zpb@B=Cd3R9z1Jy#cznz8C&Qo0ADe738Z(R_`9R;Du7#9eDnZ%S4vr*BWailu?ri;_
z&zNz983!j7t@J$n!E|6fwT3bP{_D`5&5`WaE-l>XxCW#AA6o+{>vsXIyDFkfWv5Cm#n
zK`9M&n9@cs7~#zo3P6D{(>SokVa`vRVfi&t(CWmS9NEQ>?&cyy4cKLgX_eWOybMzm8n>TlNb9~+65zPihmEd
zBWQeHjAz8Rn*dJ<7aqL-G(9Tv;RzZeB2nHS##i1)7(#d_2hRJ?AUP8Yr1yEP3EWIJ
zhWT)sqLD=Nqy+(x0{@Tl@P>J}EZyJ8pFIdy*cyn5PR;#k=#gkxct=Lm_Z0m~e5#h+uv@xZE-x)U~d<
zmr&c-^)lm~zLTzlB~Z#JqJfGcw$O=C;G~Rf+vkU$biWSA$Nc#f`?)`4UDsz!WRk&Ynfdl7>Up#%IQIupr
z3OH{M;d``zu)Ng|UFgZ95JLg#$}u00pK6@0VBatqkXJ}^rn{7j6Vmg70={#qtRat!
zoNn!4eus7&
z<(WxqiwRcyFXFB!&GhC@ICrK&@Q1uT=SpQ|v?VL<8PQSFn&45Iz}NQ
zXv~ljMtK+JR&&p$f0W1r{z#ZIV}+*nPF3jgE5@W{8S|;x^=E8&yOjDY>|s~A8JRSi
zkuRk?#^S1f#M&jH{l`XFf0I2fq+P`0g*n-mM>il));G}6y;^G8JjfLB@qEEpgc5%2
zsc){w3Hprc5b#C9BcR-Q5N`Z6n0gH2h#&$hY;sUQzJa|`;Z%ov
zG&q`%A#Q20U2^~L+4xkoeQN~T)oAK<;
z6c$l*ON>0Ja%?2ghxX`{*YKXnQ$|rP=R4NweNV2P(cAi)X-UC2(eTTjq;$Q<#+$8D%xc0Ze?F4qEBCRxk+;{6q(WVrb5QjJiC*BQG~E;DWM-
z;%IQ0Qyk%V{toG-0)WOBH9osKqlj)7jM>mVf
zx~kI+T5NrpFP@6BJ1oI9gXHue#F*$}=>Rj{qbLCTD&3WX9PXUB$Q>|ia^0Ym{EZXo
zsqIqUsxAnr$T&X~J9HIE5^D>;6Jwsg;171kyd``ud>mFT9gq^*{V(1DRP6&~e;c<@
zLX+s~k`O1+Sb>P35za9Z3c!l}||!#}9ON!GrTOb{(85dvmU7@kE1L(IdI@
z&Y9xMWUz!&$^2Ml$sl&mm~LzqfUk1B-oLZU9P|g=WVIBmMKbzqmTH}5F~&IB+d00E
zB{}R{dYrZAI*JjKSVajEq%S~!6SkEnZC2C0LLd5WZy_9k>Txv3(zk`juC#Lui*+*T
z8&cx9_hg;=G=)`tI0JI+C!%G=jnjwL)dKL_r*`$y?e1fr{CiDbP_x_69sZ^sH;;ek
zcpo$2z^;e!NcrK?Q6v7~>GXZTJ%5*@1c#RfOFuN)5WM`mPF%MU-htdkRARK*Z22)v
zT2ncJ72pgp_o;GZb@dH-HTrHJ(#N~%)xpG?ouXL7+n~9!n2|F>cKbA84H`%T3%Tze
z+W&g?;6T!nY^AnO0eO9cEBQ7AUv%bucK;w#fzG+0H27ye22`pvSAW&qn80`aD+zy$
z2PSK-QzAJ>J+BKDQ{>ot**C*f^YUZU&UP3Xhwcf<$BeXj6-q4g7aleyxDa=DVqOB#
z+NRTm_f?W?sKIss;msvp`{A1wb^`{05HfA`!EwG}?09Ce@bmR4hr=uYoAw0M&zZ4`
zK%nH^*0z0rgdi?T|D?Q?sm=7_Xa5-&GOnr`Qqn58)^_XdW0zTL#J==5c&PFDhKP;j
z#7h}0;3N%-Py57Qky186F_f-5digeI?DTl;A*R_Xu<36SY{oEVa9Bo!9#Y?s!|>4-
zLg_5Hc!E&pSoA0mTYoj<~uX5aKIVx
zV|!vfEDLsb&@;xDTSPg?Vzp|bCwVXAl^rZyD+M1jx{aUyC29+77#{ea-d7d#;v@Hi
zP5-9MVpLrA_A}1Ag|Rq~YAu|NDXSvItq+%~@2_!t;@A#=`}ezoH4$_+1`4Jugo;NI
zSU4sraos9Z^2(siz+R%8FVQ7?6jlok18@gPJMWz0;(gvt_ynWGgX9QQBD{M%wJEXC
zb{+{EB=Pc-q?PImNE{*`gU2NxAT$(jx2!?cnF&4Lh=%Qa&3Bu>C}XVPoR&$6;1+~G
z@jkSic*?B9{%+T$YWCTRbjkT+Af;Q29#X_4997Tkoczw*c8fLH2h8hkZP?`dwg5gs
zVxK|2N`@?s$?ZGcI=~}4DFC4`07HuycHE5uSv5>KnlAUj({-y87nj{7GpVQ^$d?5f
zSIqT@_-7LG8DJ;B@d!VY}?2B
zA-8eS4&8s#BchLPSFDAPBiUZnlXVaOd*{g|ZTxVqabpB$Ef2iN=Wb?cd!#zoz{`Y7
zGAQ2Gq=I<~ari=<2Fcp-d~F+371(G*i9-85kZ|KqjdMb!-sw{e`O3aVB{p&{mHA{$
zEGmkobUS5yF%|pd=*xdHUci9%(dkNuGVe)p3jc64j;~PXI0Genc=m5LqBlXk-Fg?q
zme3TPa;dcRx5Qm=jrksoLB6o)&OtHEUJv2&jb&pSrGFGuTdfq(=Ja6w1vxyXkFQ6?
zaeuvetq6&9xajPB>XA{qI4;AuY8)RIM~@&(exK#ImO8$sYg+9RqjMwG9e?J1NZ-G=
zjWRv@{kNQ|LA!Offd`>f&s&+QsnnummjNH-yTvPgrZgIkVlpQ85TlqwU!lJZ&zb66
z4*i`%G%AgPxv}sQ-A5FN;I!^Qg)WRpQrLE&hHlJ-6Ct!CFb3F*)P2UPg9BtZ|G4FO
zY(t9&XS~sVGgpPN)=7$01e)6U{q2uxYmaji4BjmjN6udDyEsSwYHH6l7tpH;L@^)%>^9nj3)%oL;@I8(b{e~9->J1uVtA2}3#hqLq$TOhRNq@nrme#a^b)Od7
zM%{1Lj`%U4@CI~^RG&lo4=&a)kVO+7&7@Q&JGzb#i_lv&<4A?)>
zvR@DB1%zj{>mX}+Z)-@om8NSqdO)7>7*+=Urp@QJwiLR}fAP3WrjI{j!~>Ipi0^VM
zw_Q0npE*xuhvW?-FmCN-W4}8*$^Y(5MCke%^!Y%?KCXW(K8j#Y^-!jsaOF0RdXO4W(8P@x9K_p
zQN{{NUtCuB;N*(mGI{jg{yMZLx>?uy-xfq@e}+jwFLk0DqPY%lZG;DF*}gyhCc+P{
zYF-+$65pPx#%Y1J(W7!Pbhs=6+pCz!ntM@rb0y6vF=W1FAYSs>z%<F46X@5aPyYaZmhh)2zXlJoii9ej
zn2C4O#t@$sn6!S7CXE+1^0)oDD@i~LO^7+|5oR#Mbm0Os=Q;+=g?_@0oXU&92AuPx
z3`LAxRVRz%#1t7ii5w%A6y>Si{|bj*>jev24g(Q)FpB&Y4xoTk`%-zvr#rUg7!tDd
zd(n3He$5w%qw|c->qp^{W_Gj
zgu7Iyv^M77GJ(Tp010dOLtCp`4k)!w+l}TCQDb6ijx5GG|0b2#en{HF5k(h%9#u0Q
zjCfw2*Rp>2N^vLf<$p2&P|$}sFs>6y`A8r7RyUYzp|x%49PUwI0qqvB=__W@Vg67w
zpa&K&Gb7N#?1}Wlw<&YW>Ww2Jj7O3T1W-Zt+feoOGv1d(Q@&qx+;jTGp}1Im4EwYA
zjpph>vK0t5*w)z`4RjMM4D(JMW`pkYz!|-3Va;nD;gs76=ZxOJi^kReFVO*c?CzDE
zE3nnPRrf+;BqDG9^Gd$N|rCGOb
z0X4Sp2K{~OV;Xnzh5KQb+r81!Zh$1LBRy0Toin*#U10D%V}Y3dH6BDbq9YUJ{*7&o
zPMO4+QesuuYQ+&Cf+-XQ`;3_C_RzrU+RDxZq~
z5ZcFC?YnjrZ!TO9Rl`SI+@+sKn8Z-yLlHfNGby#y35z%0KB?8*%QV;RMr~Re*fDs^bQ0y(eVSyJXAE2yj|ax4g^WylbtdEzgi+bk^VY*
zm{uW{ZeaG5{rgK4z|5(Adj%rI5`tV~H^#Mo{=seG&`%Ol(PsfW8wd4bbJzNdwbmi1
z*_x!dN@jwHx}!Kiks>>Ek+Ko0O=cKu7Z@wdaE%Pp<$p=GPnZ<89-6CTxQwqfZkyDcL)7l(3*&z
zmH|JCq(-Jj_=SI{-L%p29Q7&0bC5J}>-aTNhWnMpDUalj3xT{20hZ}~daEkGa$mA3
z`o*w8%j>SkHNK_DnEDJBhDM{OAPRSVXL|jif
z?wQbx{(~Guq+U$k(%D2tAyXs5{~#Q1qhB9ShR`FrfTivbv&-=KnP=auR4{)KMKNO>
zsL%9Pp{w?Q(;kTq@s-$_q^%E(uI5Sti4W$W*mXn`0p#bo%gn9+(=#=(T9Tins$MYLPf3W7f<&syw3rOx
znr6(QNc7-QTJDm_3ye#-+%@|Yi_(d%M&6ml&XU-s_YQ?@@D0-g6%e|5M5;viFB8|Y
zV%FP(Ss>{hRJ42c+t#C5iy^0o;MJ|}uk7T*&A*Kk;Z_nB&hw0t^2uF0{?JE~Tt5*6
zC-?As7a&xL+iyr&bzH#v2#5Z9=417D?$mk>21@v=7|q3*>Qom-TH=gxLV<>l
zth*-z>Ya`G8W8=W2P3K|s3K{vA9N*wpga|YOzm6e%
z8?$T)*=g21Nz}C-AXM?QLq!88Arr?UU_a5Ry1yu(TWD=E2@~GnmZf40i)3K5NJ&M)
z*=&7i^381K!B1T-MOtNaL`ISdUo|qX4xSP>;Wa^DNpV6clp$rN=;3VRS2}BZ{mGei
z1+Eu+cOUkmOBztp(0n5tpkGMk6=+I>`JBqG+^#&rq%^Cr*sLR{=%n(he{5_s=RSgH
ziqGVeiq9TgDZ$jzIgx#EG(6RQ9V}c4)r{o~xg7j#B$iJP=5i~Y$HQWev{m-A@Q=D)34LG`HS_Emc7*imnEF{j*8UwmSSI$A5g#80Xr}W(-tIrG|M+W(Ec}1zu??&=qVCot#wK+B_Pdn&+$N
z0_rwVgaEb2Q*PZ77yv?;Z0X_EE2E#W7+%uVVZSWR#t*ruFc@mABC!u&&Z$7tmhC&A
zf%M2Y1DO%-@!8C1W}|#zyW9*(cmEM&k#P{yU|&+0H^EV|qJzrU#C4F9sYGy4Q8w;&RWpI^wTl
z_h?t78JyW73ga==O-Y3!sR;!R4*K6c2}H1wV5P?;pQD-*wm~`^(8Yp-2B2hD@6Q31
z6^~`L;!x4-um5~p<*cATI|n3h*Zn#CKa(hL)ZDf{hpATqgmcw-n$9DOhH~#Ejy=@j
zIYWUM3RQJbS7U19z&Xw7_4@&6Bdf{6hZn;)Cti9)l-;8<{|3()hw#%U3?aSjs3dv(
z_sIQFF^$n(4^Dak9xV`munh}A94lRC^@@KAIKin_ng}JxjlL41=5Z-ay=qpu)wKXq
zsQtrUO*IAD(y)FVt0MFCp54*pyLVbv0sb%^VxlE;(qrhPJX_QN)97A}WaUD1G2A8z#w>nf%&b8vYI1u$--V
z(bj_3=$SGZ=CO~v6}mX7W>v6r#s=v2QUgq(^;eMO(nD?h^7&n3
z^9AYcJ4udTN-ehQJB4@WhnAJjumbZ|`aXb~lm^l|Q5Hb%=lTqlSN{G(Jq5J6omGw-
z4Rq78xS(vSNZN@Brt6e!MyF0Fqw`dG1(AN%eyQszhqi$VF5sA`Wlr`@$cUZcP=RsN
zDuRYQ^@EA0vlWbz(r7z@|1N@a?li-s?*z3A=jj6ue8mAI;Z&y|T4!%Sf$fBJ_Ez7$
zSY99n)U3KfcBUHRBAkC-`JOq0b6L7IV&sjGw7b+mPvV?C7hX4tzMMC;b>{M__0r$F
zQ<+pyqDNnq{^B?RCD2@$OhjwAQC0|1zQb^XvR+3_R984!SC&@ihHpui*w}5YX4CAG
zH!7Ox=v6qYGI3TUzRq5LUWT0`z1MM&g9|cIjClorw$pILl(%6i7{6h2>0Zzf1*w>D
zKmP+pxM3cSVwg>{sus!qOCR$DbuQ4BUc*Me@#q5z`+oTO71{l$GQaeHDxQEF4G
zpOU{5`6?gXDKCKObK-a36+hcpRS-nGS~Vq|&|}~0zl!iQ%|z=iJg5qBQxd0X;n`C;
z6uZqeqXaBF34sE^%OnJ@Fmr~Tt`Y~t@2d&*g0dmq
zw{@~ReC29%>5#-D;MHS2%!9$Qp2AukEX?@=8JMm`C-q?zC$)9ZGkciE4wN-mO0XP|
zHcVmY8DGe4i(>n2p?L+U?H(-i5mUS1aMKWZ$K5FmS$;nNSwN=0cCxN`*hg@ws_@)P
zey#>}3k9gEOzf-EMKacN_v|GisC-iLRw4r)OQ>#jv`V`p36O}UK}jClwclf8SWS{l
z4fl3%TsDc#c}%6|2Sbc}o!u3G`fTeO++#*&qrJ&grzofjHf~0^JL=c>*u;g3Ez0`D
zZ0HuBHqOG3K`xOys$gh#BOb?$%1Q(yua|C$wHg6y^xvpAEi&67kY^$bBV!$4-A@NE
zpUvWTs^tw@xXl7#K`#lk_*(nv-v1`l*1{a<)~G}bNrA|?bi?dw)jwBEa148QjP5%$
zORx(+iEZ@7o{co%WjtdPd46U2`DZ9)aHxjtxuh*~0pzMKake}>>;TBrI(4#VGF7$0
z_GNTuQkT4S{aj|-(zcsoOJE9?qDQsAtEeev1mArET