diff --git a/src/USER-MISC/fix_wall_ees.cpp b/src/USER-MISC/fix_wall_ees.cpp
index 8ccadf274ac9cbd07bfd592ca7bb9eb18db5132a..9b65d1e30f37b74dd226d19956ee1d7e57864a0d 100644
--- a/src/USER-MISC/fix_wall_ees.cpp
+++ b/src/USER-MISC/fix_wall_ees.cpp
@@ -11,13 +11,25 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
-#include "math.h"
+/* ----------------------------------------------------------------------
+ Contributing author: Abdoreza Ershadinia, a.ershadinia at gmail.com
+------------------------------------------------------------------------- */
+
+#include <math.h>
#include "math_extra.h"
#include "fix_wall_ees.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
+#include "domain.h"
+#include "region.h"
+#include "force.h"
+#include "lattice.h"
+#include "update.h"
+#include "output.h"
+#include "respa.h"
#include "error.h"
+#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@@ -96,13 +108,13 @@ void FixWallEES::wall_particle(int m, int which, double coord)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
- if (side < 0)
- delta = x[i][dim] - coord;
+ if (side < 0)
+ delta = x[i][dim] - coord;
else
- delta = coord - x[i][dim];
+ delta = coord - x[i][dim];
if (delta >= cutoff[m])
- continue;
+ continue;
double A[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
double tempvec[3]= {0,0,0};
@@ -168,37 +180,37 @@ void FixWallEES::wall_particle(int m, int which, double coord)
hms = delta - sigman;
fwall = side*(
- -1*coeff4[m]/hhss2 +
- coeff3[m] * ( 21*delta6 + 63*delta4*sigman2 + 27*delta2*sigman4 + sigman6 ) / hhss8
- );
- f[i][dim] -= fwall;
+ -1*coeff4[m]/hhss2 +
+ coeff3[m] * ( 21*delta6 + 63*delta4*sigman2 + 27*delta2*sigman4 + sigman6 ) / hhss8
+ );
+ f[i][dim] -= fwall;
- ewall[0] += -1*coeff2[m] * ( 4*delta/sigman2/hhss + 2*log(hms/hps)/sigman3 ) +
- coeff1[m] * ( 35*delta5 + 70*delta3*sigman2 + 15*delta*sigman4 ) / hhss7;
+ ewall[0] += -1*coeff2[m] * ( 4*delta/sigman2/hhss + 2*log(hms/hps)/sigman3 ) +
+ coeff1[m] * ( 35*delta5 + 70*delta3*sigman2 + 15*delta*sigman4 ) / hhss7;
- ewall[m+1] += fwall;
+ ewall[m+1] += fwall;
- twall = coeff6[m] * ( 6.*delta3/sigman4/hhss2 - 10.*delta/sigman2/hhss2 + 3.*log(hms/hps)/sigman5 ) +
- coeff5[m] * ( 21.*delta5 + 30.*delta3*sigman2 + 5.*delta*sigman4 ) / hhss8 ;
+ twall = coeff6[m] * ( 6.*delta3/sigman4/hhss2 - 10.*delta/sigman2/hhss2 + 3.*log(hms/hps)/sigman5 ) +
+ coeff5[m] * ( 21.*delta5 + 30.*delta3*sigman2 + 5.*delta*sigman4 ) / hhss8 ;
- MathExtra::matvec(Lx,nhat,tempvec);
- MathExtra::transpose_matvec(A,tempvec,tempvec2);
- for(int k = 0; k<3; k++) tempvec2[k] *= shape[k];
- that[0] = MathExtra::dot3(SAn,tempvec2);
+ MathExtra::matvec(Lx,nhat,tempvec);
+ MathExtra::transpose_matvec(A,tempvec,tempvec2);
+ for(int k = 0; k<3; k++) tempvec2[k] *= shape[k];
+ that[0] = MathExtra::dot3(SAn,tempvec2);
- MathExtra::matvec(Ly,nhat,tempvec);
- MathExtra::transpose_matvec(A,tempvec,tempvec2);
- for(int k = 0; k<3; k++) tempvec2[k] *= shape[k];
- that[1] = MathExtra::dot3(SAn,tempvec2);
+ MathExtra::matvec(Ly,nhat,tempvec);
+ MathExtra::transpose_matvec(A,tempvec,tempvec2);
+ for(int k = 0; k<3; k++) tempvec2[k] *= shape[k];
+ that[1] = MathExtra::dot3(SAn,tempvec2);
- MathExtra::matvec(Lz,nhat,tempvec);
- MathExtra::transpose_matvec(A,tempvec,tempvec2);
- for(int k = 0; k < 3; k++) tempvec2[k] *= shape[k];
- that[2] = MathExtra::dot3(SAn,tempvec2);
+ MathExtra::matvec(Lz,nhat,tempvec);
+ MathExtra::transpose_matvec(A,tempvec,tempvec2);
+ for(int k = 0; k < 3; k++) tempvec2[k] *= shape[k];
+ that[2] = MathExtra::dot3(SAn,tempvec2);
- for(int j = 0; j<3 ; j++)
+ for(int j = 0; j<3 ; j++)
tor[i][j] += twall * that[j];
}
diff --git a/src/USER-MISC/fix_wall_region_ees.cpp b/src/USER-MISC/fix_wall_region_ees.cpp
index 20a5f74c1508969c8ecf85837b43f6496b8cce7e..e070e2a2216e81c63e897fc99bafc771ecadf1fb 100644
--- a/src/USER-MISC/fix_wall_region_ees.cpp
+++ b/src/USER-MISC/fix_wall_region_ees.cpp
@@ -1,4 +1,4 @@
-/* -*- c++ -*- ----------------------------------------------------------
+/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -10,10 +10,14 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
-#include <iostream>
-#include "math.h"
-#include "stdlib.h"
-#include "string.h"
+
+/* ----------------------------------------------------------------------
+ Contributing author: Abdoreza Ershadinia, a.ershadinia at gmail.com
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
#include "fix_wall_region_ees.h"
#include "atom.h"
#include "atom_vec.h"
@@ -31,117 +35,109 @@
using namespace LAMMPS_NS;
using namespace FixConst;
-enum{LJ93,LJ126,COLLOID,HARMONIC,EES};//me
-
/* ---------------------------------------------------------------------- */
-/// USAGE:
-/// fix ID group-ID wall/region/ees region-ID epsilon sigma cutoff
-///
+
FixWallRegionEES::FixWallRegionEES(LAMMPS *lmp, int narg, char **arg) :
- Fix(lmp, narg, arg)
+ Fix(lmp, narg, arg)
{
- if (narg != 7) error->all(FLERR,"Illegal fix wall/region/ees command");
- scalar_flag = 1;
- vector_flag = 1;
- size_vector = 3;
- global_freq = 1;
- extscalar = 1;
- extvector = 1;
-
- // parse args
-
- iregion = domain->find_region(arg[3]);
- if (iregion == -1)
- error->all(FLERR,"Region ID for fix wall/region/ees does not exist");
- int n = strlen(arg[3]) + 1;
- idregion = new char[n];
- strcpy(idregion,arg[3]);
-
-
-
- epsilon = force->numeric(FLERR,arg[4]);
- sigma = force->numeric(FLERR,arg[5]);
- cutoff = force->numeric(FLERR,arg[6]);
-
- if (cutoff <= 0.0) error->all(FLERR,"Fix wall/region/ees cutoff <= 0.0");
-
- eflag = 0;
- ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
+ if (narg != 7) error->all(FLERR,"Illegal fix wall/region/ees command");
+ scalar_flag = 1;
+ vector_flag = 1;
+ size_vector = 3;
+ global_freq = 1;
+ extscalar = 1;
+ extvector = 1;
+
+ // parse args
+
+ iregion = domain->find_region(arg[3]);
+ if (iregion == -1)
+ error->all(FLERR,"Region ID for fix wall/region/ees does not exist");
+ int n = strlen(arg[3]) + 1;
+ idregion = new char[n];
+ strcpy(idregion,arg[3]);
+
+ epsilon = force->numeric(FLERR,arg[4]);
+ sigma = force->numeric(FLERR,arg[5]);
+ cutoff = force->numeric(FLERR,arg[6]);
+
+ if (cutoff <= 0.0) error->all(FLERR,"Fix wall/region/ees cutoff <= 0.0");
+
+ eflag = 0;
+ ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
}
/* ---------------------------------------------------------------------- */
FixWallRegionEES::~FixWallRegionEES()
{
- delete [] idregion;
+ delete [] idregion;
}
/* ---------------------------------------------------------------------- */
int FixWallRegionEES::setmask()
{
- int mask = 0;
- mask |= POST_FORCE;
- mask |= THERMO_ENERGY;
- mask |= POST_FORCE_RESPA;
- mask |= MIN_POST_FORCE;
- return mask;
+ int mask = 0;
+ mask |= POST_FORCE;
+ mask |= THERMO_ENERGY;
+ mask |= POST_FORCE_RESPA;
+ mask |= MIN_POST_FORCE;
+ return mask;
}
/* ---------------------------------------------------------------------- */
void FixWallRegionEES::init()
{
- // set index and check validity of region
-
- iregion = domain->find_region(idregion);
- if (iregion == -1)
- error->all(FLERR,"Region ID for fix wall/region/ees does not exist");
+ // set index and check validity of region
+ iregion = domain->find_region(idregion);
+ if (iregion == -1)
+ error->all(FLERR,"Region ID for fix wall/region/ees does not exist");
- avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- if (!avec)
- error->all(FLERR,"Fix wall/region/ees requires atom style ellipsoid");
+ avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec)
+ error->all(FLERR,"Fix wall/region/ees requires atom style ellipsoid");
- // check that all particles are finite-size ellipsoids
- // no point particles allowed, spherical is OK
+ // check that all particles are finite-size ellipsoids
+ // no point particles allowed, spherical is OK
- int *ellipsoid = atom->ellipsoid;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
+ int *ellipsoid = atom->ellipsoid;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
- for (int i = 0; i < nlocal; i++)
- if (mask[i] & groupbit)
- if (ellipsoid[i] < 0)
- error->one(FLERR,"Fix wall/region/ees requires extended particles");
+ for (int i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ if (ellipsoid[i] < 0)
+ error->one(FLERR,"Fix wall/region/ees requires extended particles");
+ // setup coefficients
- // setup coefficients for each style
+ coeff1 = ( 2. / 4725. ) * epsilon * pow(sigma,12.0);
+ coeff2 = ( 1. / 24. ) * epsilon * pow(sigma,6.0);
+ coeff3 = ( 2. / 315. ) * epsilon * pow(sigma,12.0);
+ coeff4 = ( 1. / 3. ) * epsilon * pow(sigma,6.0);
+ coeff5 = ( 4. / 315. ) * epsilon * pow(sigma,12.0);
+ coeff6 = ( 1. / 12. ) * epsilon * pow(sigma,6.0);
+ offset = 0;
- coeff1 = ( 2. / 4725. ) * epsilon * pow(sigma,12.0);
- coeff2 = ( 1. / 24. ) * epsilon * pow(sigma,6.0);
- coeff3 = ( 2. / 315. ) * epsilon * pow(sigma,12.0);
- coeff4 = ( 1. / 3. ) * epsilon * pow(sigma,6.0);
- coeff5 = ( 4. / 315. ) * epsilon * pow(sigma,12.0);
- coeff6 = ( 1. / 12. ) * epsilon * pow(sigma,6.0);
- offset = 0;
-
- if (strstr(update->integrate_style,"respa"))
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ if (strstr(update->integrate_style,"respa"))
+ nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixWallRegionEES::setup(int vflag)
{
- if (strstr(update->integrate_style,"verlet"))
- post_force(vflag);
- else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
- }
+ if (strstr(update->integrate_style,"verlet"))
+ post_force(vflag);
+ else {
+ ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+ post_force_respa(vflag,nlevels_respa-1,0);
+ ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -155,121 +151,111 @@ void FixWallRegionEES::min_setup(int vflag)
void FixWallRegionEES::post_force(int vflag)
{
- //me
- //sth is needed here, but I dont know what
- //that is calculation of sn
-
- int i,m,n;
- double rinv,fx,fy,fz,tooclose[3];//me
- double sn;//me
-
- eflag = 0;
- ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
-
- double **x = atom->x;
- double **f = atom->f;
- double *radius = atom->radius;
-
- double **tor = atom->torque; //me
-
- //avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");//me
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;//me
- int *ellipsoid = atom->ellipsoid;//me
-
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
-
- Region *region = domain->regions[iregion];
- region->prematch();
-
- int onflag = 0;
-
- // region->match() insures particle is in region or on surface, else error
- // if returned contact dist r = 0, is on surface, also an error
- // in COLLOID case, r <= radius is an error
-
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) {
- if (!region->match(x[i][0],x[i][1],x[i][2])) {
- onflag = 1;
- continue;
- }
-
- double A[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
- double tempvec[3]= {0,0,0};
- double sn2 = 0.0;
- double nhat[3] = {0,0,0};
- double* shape = bonus[ellipsoid[i]].shape;;
- MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A);
-
- for(int which = 0 ; which < 3; which ++){//me
- nhat[which]=1;
- nhat[(which+1)%3] = 0 ;
- nhat[(which+2)%3] = 0 ;
- sn2 = 0 ;
- MathExtra::transpose_matvec(A,nhat,tempvec);
- for(int k = 0; k<3; k++) tempvec[k] *= shape[k];
- for(int k = 0; k<3 ; k++) sn2 += tempvec[k]*tempvec[k];
- sn = sqrt(sn2);
- tooclose[which] = sn;
- }
-
-
-
- n = region->surface(x[i][0],x[i][1],x[i][2],cutoff);
-
- for (m = 0; m < n; m++) {
-
- if(region->contact[m].delx != 0 && region->contact[m].r <= tooclose[0] ){
- onflag = 1;
- continue;
-
- }else if (region->contact[m].dely != 0 && region->contact[m].r <= tooclose[1]){
- onflag = 1;
- continue;
- }else if (region->contact[m].delz !=0 && region->contact[m].r <= tooclose[2]){
- onflag = 1;
- continue;
- }
- else rinv = 1.0/region->contact[m].r;
-
- ees(m,i);//me
-
- ewall[0] += eng;
- fx = fwall * region->contact[m].delx * rinv;
- fy = fwall * region->contact[m].dely * rinv;
- fz = fwall * region->contact[m].delz * rinv;
- f[i][0] += fx;
- f[i][1] += fy;
- f[i][2] += fz;
-
- ewall[1] -= fx;
- ewall[2] -= fy;
- ewall[3] -= fz;
-
- tor[i][0] += torque[0];
- tor[i][1] += torque[1];
- tor[i][2] += torque[2];
-
- }
- }
-
- if (onflag) error->one(FLERR,"Particle on or inside surface of region "
- "used in fix wall/region/ees");
+ //sth is needed here, but I dont know what
+ //that is calculation of sn
+
+ int i,m,n;
+ double rinv,fx,fy,fz,sn,tooclose[3];
+
+ eflag = 0;
+ ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double **tor = atom->torque;
+
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
+
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ Region *region = domain->regions[iregion];
+ region->prematch();
+
+ int onflag = 0;
+
+ // region->match() insures particle is in region or on surface, else error
+ // if returned contact dist r = 0, is on surface, also an error
+ // in COLLOID case, r <= radius is an error
+
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ if (!region->match(x[i][0],x[i][1],x[i][2])) {
+ onflag = 1;
+ continue;
+ }
+
+ double A[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
+ double tempvec[3]= {0,0,0};
+ double sn2 = 0.0;
+ double nhat[3] = {0,0,0};
+ double* shape = bonus[ellipsoid[i]].shape;;
+ MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A);
+
+ for(int which = 0 ; which < 3; which ++){//me
+ nhat[which]=1;
+ nhat[(which+1)%3] = 0 ;
+ nhat[(which+2)%3] = 0 ;
+ sn2 = 0 ;
+ MathExtra::transpose_matvec(A,nhat,tempvec);
+ for(int k = 0; k<3; k++) tempvec[k] *= shape[k];
+ for(int k = 0; k<3 ; k++) sn2 += tempvec[k]*tempvec[k];
+ sn = sqrt(sn2);
+ tooclose[which] = sn;
+ }
+
+ n = region->surface(x[i][0],x[i][1],x[i][2],cutoff);
+
+ for (m = 0; m < n; m++) {
+
+ if (region->contact[m].delx != 0 && region->contact[m].r <= tooclose[0]){
+ onflag = 1;
+ continue;
+ } else if (region->contact[m].dely != 0 && region->contact[m].r <= tooclose[1]){
+ onflag = 1;
+ continue;
+ } else if (region->contact[m].delz !=0 && region->contact[m].r <= tooclose[2]){
+ onflag = 1;
+ continue;
+ } else rinv = 1.0/region->contact[m].r;
+
+ ees(m,i);
+
+ ewall[0] += eng;
+ fx = fwall * region->contact[m].delx * rinv;
+ fy = fwall * region->contact[m].dely * rinv;
+ fz = fwall * region->contact[m].delz * rinv;
+ f[i][0] += fx;
+ f[i][1] += fy;
+ f[i][2] += fz;
+
+ ewall[1] -= fx;
+ ewall[2] -= fy;
+ ewall[3] -= fz;
+
+ tor[i][0] += torque[0];
+ tor[i][1] += torque[1];
+ tor[i][2] += torque[2];
+ }
+ }
+
+ if (onflag) error->one(FLERR,"Particle on or inside surface of region "
+ "used in fix wall/region/ees");
}
/* ---------------------------------------------------------------------- */
void FixWallRegionEES::post_force_respa(int vflag, int ilevel, int iloop)
{
- if (ilevel == nlevels_respa-1) post_force(vflag);
+ if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixWallRegionEES::min_post_force(int vflag)
{
- post_force(vflag);
+ post_force(vflag);
}
/* ----------------------------------------------------------------------
@@ -278,13 +264,13 @@ void FixWallRegionEES::min_post_force(int vflag)
double FixWallRegionEES::compute_scalar()
{
- // only sum across procs one time
+ // only sum across procs one time
- if (eflag == 0) {
- MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
- eflag = 1;
- }
- return ewall_all[0];
+ if (eflag == 0) {
+ MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
+ eflag = 1;
+ }
+ return ewall_all[0];
}
/* ----------------------------------------------------------------------
@@ -293,18 +279,15 @@ double FixWallRegionEES::compute_scalar()
double FixWallRegionEES::compute_vector(int n)
{
- // only sum across procs one time
+ // only sum across procs one time
- if (eflag == 0) {
- MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
- eflag = 1;
- }
- return ewall_all[n+1];
+ if (eflag == 0) {
+ MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
+ eflag = 1;
+ }
+ return ewall_all[n+1];
}
-
-
-//me
/* ----------------------------------------------------------------------
EES interaction for ellipsoid particle with wall
compute eng and fwall and twall = magnitude of wall force and torque
@@ -312,94 +295,89 @@ double FixWallRegionEES::compute_vector(int n)
void FixWallRegionEES::ees(int m, int i)
{
- Region *region = domain->regions[iregion];
- region->prematch();
-
- double delta = 0.0, delta2 = 0.0, delta3 = 0.0, delta4 = 0.0, delta5 = 0.0, delta6 = 0.0;
- double sigman = 0.0, sigman2 = 0.0 , sigman3 = 0.0, sigman4 = 0.0, sigman5 = 0.0, sigman6 = 0.0;
- double hhss = 0.0, hhss2 = 0.0, hhss4 = 0.0, hhss7 = 0.0, hhss8 = 0.0; //h^2 - s_n^2
- double hps = 0.0; //h+s_n
- double hms = 0.0; //h-s_n
- double twall = 0.0;
- double tempvec[3]={0,0,0};
- double tempvec2[3]= {0,0,0};
-
- double SAn[3] = {0,0,0};
- double that[3] = {0,0,0};
-
- double Lx[3][3] = {{0,0,0},{0,0,-1},{0,1,0}};
- double Ly[3][3] = {{0,0,1},{0,0,0},{-1,0,0}};
- double Lz[3][3] = {{0,-1,0},{1,0,0},{0,0,0}};
-
- double A[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
- double nhat[3] = {0,0,0};
-
- nhat[0] = region->contact[m].delx / region->contact[m].r;
- nhat[1] = region->contact[m].dely / region->contact[m].r;
- nhat[2] = region->contact[m].delz / region->contact[m].r;
-
- //avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;
- int *ellipsoid = atom->ellipsoid;//me
-
- double* shape = bonus[ellipsoid[i]].shape;;
- MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A);
-
- MathExtra::transpose_matvec(A,nhat,tempvec);
- for(int k = 0; k<3; k++) tempvec[k] *= shape[k];
- for(int k = 0; k<3 ; k++) sigman2 += tempvec[k]*tempvec[k];
- for(int k = 0; k<3; k++) SAn[k] = tempvec[k];
-
-
- sigman = sqrt(sigman2);
- delta = fabs(region->contact[m].r);
-
- sigman3 = sigman2 * sigman;
- sigman4 = sigman2 * sigman2;
- sigman5 = sigman4 * sigman;
- sigman6 = sigman3 * sigman3;
-
- delta2 = delta * delta;
- delta3 = delta2 * delta;
- delta4 = delta2 * delta2;
- delta5 = delta3 * delta2;
- delta6 = delta3 * delta3;
-
- hhss = delta2 - sigman2;
- hhss2 = hhss * hhss;
- hhss4 = hhss2 * hhss2;
- hhss8 = hhss4 * hhss4;
- hhss7 = hhss4 * hhss2 * hhss;
-
- hps = delta + sigman;
- hms = delta - sigman;
-
- fwall = -1*coeff4/hhss2 +
- coeff3 * ( 21*delta6 + 63*delta4*sigman2 + 27*delta2*sigman4 + sigman6 ) / hhss8
- ;
-
- eng = -1*coeff2 * ( 4*delta/sigman2/hhss + 2*log(hms/hps)/sigman3 ) +
- coeff1 * ( 35*delta5 + 70*delta3*sigman2 + 15*delta*sigman4 ) / hhss7;
-
- twall = coeff6 * ( 6.*delta3/sigman4/hhss2 - 10.*delta/sigman2/hhss2 + 3.*log(hms/hps)/sigman5 ) +
- coeff5 * ( 21.*delta5 + 30.*delta3*sigman2 + 5.*delta*sigman4 ) / hhss8 ;
-
- MathExtra::matvec(Lx,nhat,tempvec);
- MathExtra::transpose_matvec(A,tempvec,tempvec2);
- for(int k = 0; k<3; k++) tempvec2[k] *= shape[k];
- that[0] = MathExtra::dot3(SAn,tempvec2);
-
- MathExtra::matvec(Ly,nhat,tempvec);
- MathExtra::transpose_matvec(A,tempvec,tempvec2);
- for(int k = 0; k<3; k++) tempvec2[k] *= shape[k];
- that[1] = MathExtra::dot3(SAn,tempvec2);
-
- MathExtra::matvec(Lz,nhat,tempvec);
- MathExtra::transpose_matvec(A,tempvec,tempvec2);
- for(int k = 0; k < 3; k++) tempvec2[k] *= shape[k];
- that[2] = MathExtra::dot3(SAn,tempvec2);
-
- for(int j = 0; j<3 ; j++)
- torque[j] = twall * that[j];
-
+ Region *region = domain->regions[iregion];
+ region->prematch();
+
+ double delta, delta2, delta3, delta4, delta5, delta6;
+ double sigman, sigman2 , sigman3, sigman4, sigman5, sigman6;
+ double hhss, hhss2, hhss4, hhss7, hhss8; //h^2 - s_n^2
+ double hps; //h+s_n
+ double hms; //h-s_n
+ double twall;
+
+ double A[3][3], nhat[3], SAn[3], that[3];
+
+ double tempvec[3]= {0,0,0};
+ double tempvec2[3]= {0,0,0};
+
+ double Lx[3][3] = {{0,0,0},{0,0,-1},{0,1,0}};
+ double Ly[3][3] = {{0,0,1},{0,0,0},{-1,0,0}};
+ double Lz[3][3] = {{0,-1,0},{1,0,0},{0,0,0}};
+
+ nhat[0] = region->contact[m].delx / region->contact[m].r;
+ nhat[1] = region->contact[m].dely / region->contact[m].r;
+ nhat[2] = region->contact[m].delz / region->contact[m].r;
+
+ AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ int *ellipsoid = atom->ellipsoid;
+
+ double* shape = bonus[ellipsoid[i]].shape;;
+ MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A);
+
+ sigman2 = 0.0;
+ MathExtra::transpose_matvec(A,nhat,tempvec);
+ for(int k = 0; k<3; k++) tempvec[k] *= shape[k];
+ for(int k = 0; k<3; k++) sigman2 += tempvec[k]*tempvec[k];
+ for(int k = 0; k<3; k++) SAn[k] = tempvec[k];
+
+ sigman = sqrt(sigman2);
+ delta = fabs(region->contact[m].r);
+
+ sigman3 = sigman2 * sigman;
+ sigman4 = sigman2 * sigman2;
+ sigman5 = sigman4 * sigman;
+ sigman6 = sigman3 * sigman3;
+
+ delta2 = delta * delta;
+ delta3 = delta2 * delta;
+ delta4 = delta2 * delta2;
+ delta5 = delta3 * delta2;
+ delta6 = delta3 * delta3;
+
+ hhss = delta2 - sigman2;
+ hhss2 = hhss * hhss;
+ hhss4 = hhss2 * hhss2;
+ hhss8 = hhss4 * hhss4;
+ hhss7 = hhss4 * hhss2 * hhss;
+
+ hps = delta + sigman;
+ hms = delta - sigman;
+
+ fwall = -1*coeff4/hhss2 + coeff3
+ * (21*delta6 + 63*delta4*sigman2 + 27*delta2*sigman4 + sigman6) / hhss8;
+
+ eng = -1*coeff2 * (4*delta/sigman2/hhss + 2*log(hms/hps)/sigman3) +
+ coeff1 * (35*delta5 + 70*delta3*sigman2 + 15*delta*sigman4) / hhss7;
+
+ twall = coeff6 * (6*delta3/sigman4/hhss2 - 10*delta/sigman2/hhss2
+ + 3*log(hms/hps)/sigman5)
+ + coeff5 * (21.*delta5 + 30.*delta3*sigman2 + 5.*delta*sigman4) / hhss8;
+
+ MathExtra::matvec(Lx,nhat,tempvec);
+ MathExtra::transpose_matvec(A,tempvec,tempvec2);
+ for(int k = 0; k<3; k++) tempvec2[k] *= shape[k];
+ that[0] = MathExtra::dot3(SAn,tempvec2);
+
+ MathExtra::matvec(Ly,nhat,tempvec);
+ MathExtra::transpose_matvec(A,tempvec,tempvec2);
+ for(int k = 0; k<3; k++) tempvec2[k] *= shape[k];
+ that[1] = MathExtra::dot3(SAn,tempvec2);
+
+ MathExtra::matvec(Lz,nhat,tempvec);
+ MathExtra::transpose_matvec(A,tempvec,tempvec2);
+ for(int k = 0; k < 3; k++) tempvec2[k] *= shape[k];
+ that[2] = MathExtra::dot3(SAn,tempvec2);
+
+ for(int j = 0; j<3 ; j++)
+ torque[j] = twall * that[j];
}
diff --git a/src/USER-MISC/fix_wall_region_ees.h b/src/USER-MISC/fix_wall_region_ees.h
index 59679a0b41ce6dfd4af3dbcc8a6069cb0139f812..9da2cccce624ed68f054aacccb401089721eda62 100644
--- a/src/USER-MISC/fix_wall_region_ees.h
+++ b/src/USER-MISC/fix_wall_region_ees.h
@@ -40,7 +40,7 @@ class FixWallRegionEES : public Fix {
double compute_vector(int);
private:
- class AtomVecEllipsoid *avec;//me
+ class AtomVecEllipsoid *avec;
int iregion;
double epsilon,sigma,cutoff;
@@ -50,11 +50,11 @@ class FixWallRegionEES : public Fix {
char *idregion;
double coeff1,coeff2,coeff3,coeff4,offset;
- double coeff5, coeff6;//me
+ double coeff5, coeff6;
double eng,fwall;
- double torque[3];//me
+ double torque[3];
- void ees(int, int);//me
+ void ees(int, int);
};
}