Commit 6eddc1a2 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

coding style and whitespace cleanup to match LAMMPS style

parent 9fa4588e
Loading
Loading
Loading
Loading
+39 −27
Original line number Diff line number Diff line
@@ -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;
+266 −288
Original line number Diff line number Diff line
/* -*- 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,12 +35,8 @@
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)
{
@@ -57,8 +57,6 @@ FixWallRegionEES::FixWallRegionEES(LAMMPS *lmp, int narg, char **arg) :
  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]);
@@ -98,7 +96,6 @@ void FixWallRegionEES::init()
  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");
@@ -115,8 +112,7 @@ void FixWallRegionEES::init()
      if (ellipsoid[i] < 0)
        error->one(FLERR,"Fix wall/region/ees requires extended particles");


    // setup coefficients for each style
  // setup coefficients

  coeff1 = ( 2. / 4725. ) * epsilon * pow(sigma,12.0);
  coeff2 = ( 1. / 24. ) * epsilon * pow(sigma,6.0);
@@ -155,26 +151,21 @@ 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
  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 *radius = atom->radius;

    double **tor = atom->torque; //me
  double **tor = atom->torque;

    //avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");//me
    AtomVecEllipsoid::Bonus *bonus = avec->bonus;//me
    int *ellipsoid = atom->ellipsoid;//me
  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
  int *ellipsoid = atom->ellipsoid;

  int *mask = atom->mask;
  int nlocal = atom->nlocal;
@@ -214,8 +205,6 @@ void FixWallRegionEES::post_force(int vflag)
        tooclose[which] = sn;
      }



      n = region->surface(x[i][0],x[i][1],x[i][2],cutoff);

      for (m = 0; m < n; m++) {
@@ -223,17 +212,15 @@ void FixWallRegionEES::post_force(int vflag)
        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;
        } else rinv = 1.0/region->contact[m].r;

                ees(m,i);//me
        ees(m,i);

        ewall[0] += eng;
        fx = fwall * region->contact[m].delx * rinv;
@@ -250,7 +237,6 @@ void FixWallRegionEES::post_force(int vflag)
        tor[i][0] += torque[0];
        tor[i][1] += torque[1];
        tor[i][2] += torque[2];

      }
    }

@@ -302,9 +288,6 @@ double FixWallRegionEES::compute_vector(int n)
  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
@@ -315,42 +298,38 @@ 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 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 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
  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);

@@ -374,15 +353,15 @@ void FixWallRegionEES::ees(int m, int i)
  hps = delta + sigman;
  hms = delta - sigman;

    fwall =  -1*coeff4/hhss2 +
            coeff3 * (  21*delta6 + 63*delta4*sigman2 + 27*delta2*sigman4 + sigman6  ) / hhss8
            ;
  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    ;
  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);
@@ -401,5 +380,4 @@ void FixWallRegionEES::ees(int m, int i)

  for(int j = 0; j<3 ; j++)
    torque[j] = twall * that[j];

}
+4 −4
Original line number Diff line number Diff line
@@ -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);
};

}