diff --git a/src/FLD/pair_brownian.cpp b/src/FLD/pair_brownian.cpp
deleted file mode 100755
index c1cc523d3262c1885a1dcb485a907b8f36094f5c..0000000000000000000000000000000000000000
--- a/src/FLD/pair_brownian.cpp
+++ /dev/null
@@ -1,733 +0,0 @@
-/* ----------------------------------------------------------------------
- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
- http://lammps.sandia.gov, Sandia National Laboratories
- Steve Plimpton, sjplimp@sandia.gov
-
- Copyright (2003) Sandia Corporation. Under the terms of Contract
- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
- certain rights in this software. This software is distributed under
- the GNU General Public License.
-
- See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing authors: Amit Kumar and Michael Bybee (UIUC)
-------------------------------------------------------------------------- */
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include "pair_brownian.h"
-#include "atom.h"
-#include "atom_vec.h"
-#include "comm.h"
-#include "force.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
-#include "domain.h"
-#include "update.h"
-#include "modify.h"
-#include "fix.h"
-#include "fix_deform.h"
-#include "fix_wall.h"
-#include "input.h"
-#include "variable.h"
-#include "random_mars.h"
-#include "math_const.h"
-#include "math_special.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-using namespace MathConst;
-using namespace MathSpecial;
-
-// same as fix_wall.cpp
-
-enum{EDGE,CONSTANT,VARIABLE};
-
-/* ---------------------------------------------------------------------- */
-
-PairBrownian::PairBrownian(LAMMPS *lmp) : Pair(lmp)
-{
- single_enable = 0;
- random = NULL;
-}
-
-/* ---------------------------------------------------------------------- */
-
-PairBrownian::~PairBrownian()
-{
- if (allocated) {
- memory->destroy(setflag);
- memory->destroy(cutsq);
-
- memory->destroy(cut);
- memory->destroy(cut_inner);
- }
- delete random;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairBrownian::compute(int eflag, int vflag)
-{
- int i,j,ii,jj,inum,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;
-
- if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = vflag_fdotr = 0;
-
- double **x = atom->x;
- double **f = atom->f;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int newton_pair = force->newton_pair;
-
- double vxmu2f = force->vxmu2f;
- double randr;
- double prethermostat;
- double xl[3],a_sq,a_sh,a_pu,Fbmag;
- double p1[3],p2[3],p3[3];
- int overlaps = 0;
-
- // This section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (int j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu*cube(rad);
- //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*cube(rad)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
-
- // scale factor for Brownian moments
-
- prethermostat = sqrt(24.0*force->boltz*t_target/update->dt);
- prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e);
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- // FLD contribution to force and torque due to isotropic terms
-
- if (flagfld) {
- f[i][0] += prethermostat*sqrt(R0)*(random->uniform()-0.5);
- f[i][1] += prethermostat*sqrt(R0)*(random->uniform()-0.5);
- f[i][2] += prethermostat*sqrt(R0)*(random->uniform()-0.5);
- if (flaglog) {
- torque[i][0] += prethermostat*sqrt(RT0)*(random->uniform()-0.5);
- torque[i][1] += prethermostat*sqrt(RT0)*(random->uniform()-0.5);
- torque[i][2] += prethermostat*sqrt(RT0)*(random->uniform()-0.5);
- }
- }
-
- if (!flagHI) continue;
-
- 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*cube(radi)*(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 = random->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 = random->uniform()-0.5;
- fx += Fbmag*randr*p2[0];
- fy += Fbmag*randr*p2[1];
- fz += Fbmag*randr*p2[2];
-
- randr = random->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 = random->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 = random->uniform()-0.5;
- tx = Fbmag*randr*p2[0];
- ty = Fbmag*randr*p2[1];
- tz = Fbmag*randr*p2[2];
-
- randr = random->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);
- }
- }
- }
-
- int print_overlaps = 0;
- if (print_overlaps && overlaps)
- printf("Number of overlaps=%d\n",overlaps);
-
- if (vflag_fdotr) virial_fdotr_compute();
-}
-
-/* ----------------------------------------------------------------------
- allocate all arrays
-------------------------------------------------------------------------- */
-
-void PairBrownian::allocate()
-{
- allocated = 1;
- int n = atom->ntypes;
-
- memory->create(setflag,n+1,n+1,"pair:setflag");
- for (int i = 1; i <= n; i++)
- for (int j = i; j <= n; j++)
- setflag[i][j] = 0;
-
- memory->create(cutsq,n+1,n+1,"pair:cutsq");
-
- memory->create(cut,n+1,n+1,"pair:cut");
- memory->create(cut_inner,n+1,n+1,"pair:cut_inner");
-}
-
-/* ----------------------------------------------------------------------
- global settings
-------------------------------------------------------------------------- */
-
-void PairBrownian::settings(int narg, char **arg)
-{
- if (narg != 7 && narg != 9) error->all(FLERR,"Illegal pair_style command");
-
- mu = force->numeric(FLERR,arg[0]);
- flaglog = force->inumeric(FLERR,arg[1]);
- flagfld = force->inumeric(FLERR,arg[2]);
- cut_inner_global = force->numeric(FLERR,arg[3]);
- cut_global = force->numeric(FLERR,arg[4]);
- t_target = force->numeric(FLERR,arg[5]);
- seed = force->inumeric(FLERR,arg[6]);
-
- flagHI = flagVF = 1;
- if (narg == 9) {
- flagHI = force->inumeric(FLERR,arg[7]);
- flagVF = force->inumeric(FLERR,arg[8]);
- }
-
- if (flaglog == 1 && flagHI == 0) {
- error->warning(FLERR,"Cannot include log terms without 1/r terms; "
- "setting flagHI to 1");
- flagHI = 1;
- }
-
- // initialize Marsaglia RNG with processor-unique seed
-
- delete random;
- random = new RanMars(lmp,seed + comm->me);
-
- // reset cutoffs that have been explicitly set
-
- if (allocated) {
- for (int i = 1; i <= atom->ntypes; i++)
- for (int j = i+1; j <= atom->ntypes; j++)
- if (setflag[i][j]) {
- cut_inner[i][j] = cut_inner_global;
- cut[i][j] = cut_global;
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- set coeffs for one or more type pairs
-------------------------------------------------------------------------- */
-
-void PairBrownian::coeff(int narg, char **arg)
-{
- if (narg != 2 && narg != 4)
- error->all(FLERR,"Incorrect args for pair coefficients");
-
- if (!allocated) allocate();
-
- int ilo,ihi,jlo,jhi;
- force->bounds(arg[0],atom->ntypes,ilo,ihi);
- force->bounds(arg[1],atom->ntypes,jlo,jhi);
-
- double cut_inner_one = cut_inner_global;
- double cut_one = cut_global;
-
- if (narg == 4) {
- cut_inner_one = force->numeric(FLERR,arg[2]);
- cut_one = force->numeric(FLERR,arg[3]);
- }
-
- int count = 0;
- for (int i = ilo; i <= ihi; i++)
- for (int j = MAX(jlo,i); j <= jhi; j++) {
- cut_inner[i][j] = cut_inner_one;
- cut[i][j] = cut_one;
- setflag[i][j] = 1;
- count++;
- }
-
- if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
-}
-
-/* ----------------------------------------------------------------------
- init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairBrownian::init_style()
-{
- if (!atom->sphere_flag)
- error->all(FLERR,"Pair brownian requires atom style sphere");
-
- // if newton off, forces between atoms ij will be double computed
- // using different random numbers
-
- if (force->newton_pair == 0 && comm->me == 0)
- error->warning(FLERR,
- "Pair brownian needs newton pair on for "
- "momentum conservation");
-
- neighbor->request(this,instance_me);
-
- // insure all particles are finite-size
- // for pair hybrid, should limit test to types using the pair style
-
- double *radius = atom->radius;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (radius[i] == 0.0)
- error->one(FLERR,"Pair brownian requires extended particles");
-
- // require monodisperse system with same radii for all types
-
- double radtype;
- for (int i = 1; i <= atom->ntypes; i++) {
- if (!atom->radius_consistency(i,radtype))
- error->all(FLERR,"Pair brownian requires monodisperse particles");
- if (i > 1 && radtype != rad)
- error->all(FLERR,"Pair brownian requires monodisperse particles");
- rad = radtype;
- }
-
- // set the isotropic constants that depend on the volume fraction
- // vol_T = total volume
- // check for fix deform, if exists it must use "remap v"
- // If box will change volume, set appropriate flag so that volume
- // and v.f. corrections are re-calculated at every step.
- //
- // If available volume is different from box volume
- // due to walls, set volume appropriately; if walls will
- // move, set appropriate flag so that volume and v.f. corrections
- // are re-calculated at every step.
-
- flagdeform = flagwall = 0;
- for (int i = 0; i < modify->nfix; i++){
- if (strcmp(modify->fix[i]->style,"deform") == 0)
- flagdeform = 1;
- else if (strstr(modify->fix[i]->style,"wall") != NULL) {
- if (flagwall)
- error->all(FLERR,
- "Cannot use multiple fix wall commands with pair brownian");
- flagwall = 1; // Walls exist
- wallfix = (FixWall *) modify->fix[i];
- if (wallfix->xflag) flagwall = 2; // Moving walls exist
- }
- }
-
- // set the isotropic constants depending on the volume fraction
- // vol_T = total volumeshearing = flagdeform = flagwall = 0;
-
- double vol_T, wallcoord;
- if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
- else {
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallfix->xindex[m] = input->variable->find(wallfix->xstr[m]);
- // Since fix->wall->init happens after pair->init_style
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
-
- else wallcoord = wallfix->coord0[m];
-
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
- (wallhi[2] - walllo[2]);
- }
-
- // vol_P = volume of particles, assuming mono-dispersity
- // vol_f = volume fraction
-
- vol_P = atom->natoms*(4.0/3.0)*MY_PI*cube(rad);
-
- double vol_f = vol_P/vol_T;
-
- // set isotropic constants
- if (!flagVF) vol_f = 0;
-
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu*cube(rad); // not actually needed
- } else {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*cube(rad)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- }
-}
-
-/* ----------------------------------------------------------------------
- init for one type pair i,j and corresponding j,i
-------------------------------------------------------------------------- */
-
-double PairBrownian::init_one(int i, int j)
-{
- if (setflag[i][j] == 0) {
- cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
- cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
- }
-
- cut_inner[j][i] = cut_inner[i][j];
-
- return cut[i][j];
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes to restart file
-------------------------------------------------------------------------- */
-
-void PairBrownian::write_restart(FILE *fp)
-{
- write_restart_settings(fp);
-
- int i,j;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i; j <= atom->ntypes; j++) {
- fwrite(&setflag[i][j],sizeof(int),1,fp);
- if (setflag[i][j]) {
- fwrite(&cut_inner[i][j],sizeof(double),1,fp);
- fwrite(&cut[i][j],sizeof(double),1,fp);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- proc 0 reads from restart file, bcasts
-------------------------------------------------------------------------- */
-
-void PairBrownian::read_restart(FILE *fp)
-{
- read_restart_settings(fp);
- allocate();
-
- int i,j;
- int me = comm->me;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i; j <= atom->ntypes; j++) {
- if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
- MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
- if (setflag[i][j]) {
- if (me == 0) {
- fread(&cut_inner[i][j],sizeof(double),1,fp);
- fread(&cut[i][j],sizeof(double),1,fp);
- }
- MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes to restart file
-------------------------------------------------------------------------- */
-
-void PairBrownian::write_restart_settings(FILE *fp)
-{
- fwrite(&mu,sizeof(double),1,fp);
- fwrite(&flaglog,sizeof(int),1,fp);
- fwrite(&flagfld,sizeof(int),1,fp);
- fwrite(&cut_inner_global,sizeof(double),1,fp);
- fwrite(&cut_global,sizeof(double),1,fp);
- fwrite(&t_target,sizeof(double),1,fp);
- fwrite(&seed,sizeof(int),1,fp);
- fwrite(&offset_flag,sizeof(int),1,fp);
- fwrite(&mix_flag,sizeof(int),1,fp);
- fwrite(&flagHI,sizeof(int),1,fp);
- fwrite(&flagVF,sizeof(int),1,fp);
-}
-
-/* ----------------------------------------------------------------------
- proc 0 reads from restart file, bcasts
-------------------------------------------------------------------------- */
-
-void PairBrownian::read_restart_settings(FILE *fp)
-{
- int me = comm->me;
- if (me == 0) {
- fread(&mu,sizeof(double),1,fp);
- fread(&flaglog,sizeof(int),1,fp);
- fread(&flagfld,sizeof(int),1,fp);
- fread(&cut_inner_global,sizeof(double),1,fp);
- fread(&cut_global,sizeof(double),1,fp);
- fread(&t_target, sizeof(double),1,fp);
- fread(&seed, sizeof(int),1,fp);
- fread(&offset_flag,sizeof(int),1,fp);
- fread(&mix_flag,sizeof(int),1,fp);
- fread(&flagHI,sizeof(int),1,fp);
- fread(&flagVF,sizeof(int),1,fp);
- }
- MPI_Bcast(&mu,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&flaglog,1,MPI_INT,0,world);
- MPI_Bcast(&flagfld,1,MPI_INT,0,world);
- MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&t_target,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&seed,1,MPI_INT,0,world);
- MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
- MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
- MPI_Bcast(&flagHI,1,MPI_INT,0,world);
- MPI_Bcast(&flagVF,1,MPI_INT,0,world);
-
- // additional setup based on restart parameters
-
- delete random;
- random = new RanMars(lmp,seed + comm->me);
-}
-
-/* ----------------------------------------------------------------------*/
-
-void PairBrownian::set_3_orthogonal_vectors(double p1[3],
- double p2[3], double p3[3])
-{
- double norm;
- int ix,iy,iz;
-
- // find the index of maximum magnitude and store it in iz
-
- if (fabs(p1[0]) > fabs(p1[1])) {
- iz=0;
- ix=1;
- iy=2;
- } else {
- iz=1;
- ix=2;
- iy=0;
- }
-
- if (iz==0) {
- if (fabs(p1[0]) < fabs(p1[2])) {
- iz = 2;
- ix = 0;
- iy = 1;
- }
- } else {
- if (fabs(p1[1]) < fabs(p1[2])) {
- iz = 2;
- ix = 0;
- iy = 1;
- }
- }
-
- // set p2 arbitrarily such that it's orthogonal to p1
-
- p2[ix]=1.0;
- p2[iy]=1.0;
- p2[iz] = -(p1[ix]*p2[ix] + p1[iy]*p2[iy])/p1[iz];
-
- // normalize p2
-
- norm = sqrt(p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2]);
-
- p2[0] = p2[0]/norm;
- p2[1] = p2[1]/norm;
- p2[2] = p2[2]/norm;
-
- // Set p3 by taking the cross product p3=p2xp1
-
- p3[0] = p1[1]*p2[2] - p1[2]*p2[1];
- p3[1] = p1[2]*p2[0] - p1[0]*p2[2];
- p3[2] = p1[0]*p2[1] - p1[1]*p2[0];
-}
diff --git a/src/FLD/pair_brownian.h b/src/FLD/pair_brownian.h
deleted file mode 100644
index 874507852268f49ab24b2f0b5b1aabcacd89c613..0000000000000000000000000000000000000000
--- a/src/FLD/pair_brownian.h
+++ /dev/null
@@ -1,102 +0,0 @@
-/* -*- 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 PAIR_CLASS
-
-PairStyle(brownian,PairBrownian)
-
-#else
-
-#ifndef LMP_PAIR_BROWNIAN_H
-#define LMP_PAIR_BROWNIAN_H
-
-#include "pair.h"
-
-namespace LAMMPS_NS {
-
-class PairBrownian : public Pair {
- public:
- PairBrownian(class LAMMPS *);
- virtual ~PairBrownian();
- virtual void compute(int, int);
- void settings(int, char **);
- void coeff(int, char **);
- virtual double init_one(int, int);
- virtual void init_style();
- void write_restart(FILE *);
- void read_restart(FILE *);
- void write_restart_settings(FILE *);
- void read_restart_settings(FILE *);
-
- protected:
- double cut_inner_global,cut_global;
- double t_target,mu;
- int flaglog,flagfld;
- int flagHI, flagVF;
- int flagdeform, flagwall;
- double vol_P;
- double rad;
- class FixWall *wallfix;
-
- int seed;
- double **cut_inner,**cut;
- double R0,RT0;
-
- class RanMars *random;
-
- void set_3_orthogonal_vectors(double*,double*,double*);
- void allocate();
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Illegal ... command
-
-Self-explanatory. Check the input script syntax and compare to the
-documentation for the command. You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-W: Cannot include log terms without 1/r terms; setting flagHI to 1
-
-Self-explanatory.
-
-E: Incorrect args for pair coefficients
-
-Self-explanatory. Check the input script or data file.
-
-E: Pair brownian requires atom style sphere
-
-Self-explanatory.
-
-W: Pair brownian needs newton pair on for momentum conservation
-
-Self-explanatory.
-
-E: Pair brownian requires extended particles
-
-One of the particles has radius 0.0.
-
-E: Pair brownian requires monodisperse particles
-
-All particles must be the same finite size.
-
-E: Cannot use multiple fix wall commands with pair brownian
-
-Self-explanatory.
-
-*/
diff --git a/src/FLD/pair_brownian_poly.cpp b/src/FLD/pair_brownian_poly.cpp
deleted file mode 100644
index 84c56109d13cbb6f788cd65135a8ca9a268ca853..0000000000000000000000000000000000000000
--- a/src/FLD/pair_brownian_poly.cpp
+++ /dev/null
@@ -1,441 +0,0 @@
-/* ----------------------------------------------------------------------
- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
- http://lammps.sandia.gov, Sandia National Laboratories
- Steve Plimpton, sjplimp@sandia.gov
-
- Copyright (2003) Sandia Corporation. Under the terms of Contract
- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
- certain rights in this software. This software is distributed under
- the GNU General Public License.
-
- See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing authors: Amit Kumar and Michael Bybee (UIUC)
- Dave Heine (Corning), polydispersity
-------------------------------------------------------------------------- */
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include "pair_brownian_poly.h"
-#include "atom.h"
-#include "atom_vec.h"
-#include "comm.h"
-#include "force.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
-#include "domain.h"
-#include "update.h"
-#include "modify.h"
-#include "fix.h"
-#include "fix_deform.h"
-#include "fix_wall.h"
-#include "input.h"
-#include "variable.h"
-#include "random_mars.h"
-#include "math_const.h"
-#include "math_special.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-using namespace MathConst;
-using namespace MathSpecial;
-
-// same as fix_wall.cpp
-
-enum{EDGE,CONSTANT,VARIABLE};
-
-/* ---------------------------------------------------------------------- */
-
-PairBrownianPoly::PairBrownianPoly(LAMMPS *lmp) : PairBrownian(lmp)
-{
- no_virial_fdotr_compute = 1;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairBrownianPoly::compute(int eflag, int vflag)
-{
- int i,j,ii,jj,inum,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;
-
- if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = vflag_fdotr = 0;
-
- double **x = atom->x;
- double **f = atom->f;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
-
- 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];
-
- // this section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu*cube(rad);
- //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*cube(rad)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
- // scale factor for Brownian moments
-
- prethermostat = sqrt(24.0*force->boltz*t_target/update->dt);
- prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e);
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- // FLD contribution to force and torque due to isotropic terms
-
- if (flagfld) {
- f[i][0] += prethermostat*sqrt(R0*radi)*(random->uniform()-0.5);
- f[i][1] += prethermostat*sqrt(R0*radi)*(random->uniform()-0.5);
- f[i][2] += prethermostat*sqrt(R0*radi)*(random->uniform()-0.5);
- if (flaglog) {
- const double radi3 = radi*radi*radi;
- torque[i][0] += prethermostat*sqrt(RT0*radi3) *
- (random->uniform()-0.5);
- torque[i][1] += prethermostat*sqrt(RT0*radi3) *
- (random->uniform()-0.5);
- torque[i][2] += prethermostat*sqrt(RT0*radi3) *
- (random->uniform()-0.5);
- }
- }
-
- if (!flagHI) continue;
-
- 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/cube(beta1)*log(1.0/h_sep);
- a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*cube(beta0) +
- powint(beta0,4))/21.0/powint(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/cube(beta1) *
- log(1.0/h_sep);
- a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*cube(beta0) +
- 16.0*powint(beta0,4))/375.0/powint(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 *
- cube(beta0))/250.0/cube(beta1)*h_sep*log(1.0/h_sep);
- a_pu *= 8.0*MY_PI*mu*cube(radi);
-
- } 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 = random->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 = random->uniform()-0.5;
- fx += Fbmag*randr*p2[0];
- fy += Fbmag*randr*p2[1];
- fz += Fbmag*randr*p2[2];
-
- randr = random->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 = random->uniform()-0.5;
- tx = Fbmag*randr*p2[0];
- ty = Fbmag*randr*p2[1];
- tz = Fbmag*randr*p2[2];
-
- randr = random->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;
-
- }
-
- // set j = nlocal so that only I gets tallied
-
- if (evflag) ev_tally_xyz(i,nlocal,nlocal,0,
- 0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
- }
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairBrownianPoly::init_style()
-{
- if (force->newton_pair == 1)
- error->all(FLERR,"Pair brownian/poly requires newton pair off");
- if (!atom->sphere_flag)
- error->all(FLERR,"Pair brownian/poly requires atom style sphere");
-
- // insure all particles are finite-size
- // for pair hybrid, should limit test to types using the pair style
-
- double *radius = atom->radius;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (radius[i] == 0.0)
- error->one(FLERR,"Pair brownian/poly requires extended particles");
-
- int irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->full = 1;
-
- // set the isotropic constants that depend on the volume fraction
- // vol_T = total volume
- // check for fix deform, if exists it must use "remap v"
- // If box will change volume, set appropriate flag so that volume
- // and v.f. corrections are re-calculated at every step.
- //
- // If available volume is different from box volume
- // due to walls, set volume appropriately; if walls will
- // move, set appropriate flag so that volume and v.f. corrections
- // are re-calculated at every step.
-
- flagdeform = flagwall = 0;
- for (int i = 0; i < modify->nfix; i++){
- if (strcmp(modify->fix[i]->style,"deform") == 0)
- flagdeform = 1;
- else if (strstr(modify->fix[i]->style,"wall") != NULL) {
- if (flagwall)
- error->all(FLERR,
- "Cannot use multiple fix wall commands with pair brownian");
- flagwall = 1; // Walls exist
- wallfix = (FixWall *) modify->fix[i];
- if (wallfix->xflag) flagwall = 2; // Moving walls exist
- }
- }
-
- // set the isotropic constants that depend on the volume fraction
- // vol_T = total volume
-
- double vol_T, wallcoord;
- if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
- else {
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallfix->xindex[m] = input->variable->find(wallfix->xstr[m]);
- // Since fix->wall->init happens after pair->init_style
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
-
- else wallcoord = wallfix->coord0[m];
-
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
- (wallhi[2] - walllo[2]);
- }
-
- // vol_P = volume of particles, assuming mono-dispersity
- // vol_f = volume fraction
-
- double volP = 0.0;
- for (int i = 0; i < nlocal; i++)
- volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3.0);
- MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world);
-
- double vol_f = vol_P/vol_T;
-
- if (!flagVF) vol_f = 0;
- // set isotropic constants
-
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu;
- } else {
- R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- }
-}
-
-/* ----------------------------------------------------------------------
- init for one type pair i,j and corresponding j,i
-------------------------------------------------------------------------- */
-
-double PairBrownianPoly::init_one(int i, int j)
-{
- if (setflag[i][j] == 0) {
- cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
- cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
- }
-
- cut_inner[j][i] = cut_inner[i][j];
- return cut[i][j];
-}
diff --git a/src/FLD/pair_brownian_poly.h b/src/FLD/pair_brownian_poly.h
deleted file mode 100644
index 7808583876f0ede91d0c0c5023ae8e26c04ffe71..0000000000000000000000000000000000000000
--- a/src/FLD/pair_brownian_poly.h
+++ /dev/null
@@ -1,59 +0,0 @@
-/* -*- 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 PAIR_CLASS
-
-PairStyle(brownian/poly,PairBrownianPoly)
-
-#else
-
-#ifndef LMP_PAIR_BROWNIAN_POLY_H
-#define LMP_PAIR_BROWNIAN_POLY_H
-
-#include "pair_brownian.h"
-
-namespace LAMMPS_NS {
-
-class PairBrownianPoly : public PairBrownian {
- public:
- PairBrownianPoly(class LAMMPS *);
- ~PairBrownianPoly() {}
- void compute(int, int);
- double init_one(int, int);
- void init_style();
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Pair brownian/poly requires newton pair off
-
-Self-explanatory.
-
-E: Pair brownian/poly requires atom style sphere
-
-Self-explanatory.
-
-E: Pair brownian/poly requires extended particles
-
-One of the particles has radius 0.0.
-
-E: Cannot use multiple fix wall commands with pair brownian
-
-Self-explanatory.
-
-*/
diff --git a/src/FLD/pair_lubricate.cpp b/src/FLD/pair_lubricate.cpp
deleted file mode 100755
index ea398c340a9467711ad95fb19cab7573b5eae20e..0000000000000000000000000000000000000000
--- a/src/FLD/pair_lubricate.cpp
+++ /dev/null
@@ -1,820 +0,0 @@
-/* ----------------------------------------------------------------------
- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
- http://lammps.sandia.gov, Sandia National Laboratories
- Steve Plimpton, sjplimp@sandia.gov
-
- Copyright (2003) Sandia Corporation. Under the terms of Contract
- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
- certain rights in this software. This software is distributed under
- the GNU General Public License.
-
- See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing authors: Randy Schunk (SNL)
- Amit Kumar and Michael Bybee (UIUC)
-------------------------------------------------------------------------- */
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include "pair_lubricate.h"
-#include "atom.h"
-#include "atom_vec.h"
-#include "comm.h"
-#include "force.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
-#include "domain.h"
-#include "modify.h"
-#include "fix.h"
-#include "fix_deform.h"
-#include "fix_wall.h"
-#include "input.h"
-#include "variable.h"
-#include "random_mars.h"
-#include "math_const.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-using namespace MathConst;
-
-// same as fix_deform.cpp
-
-enum{NO_REMAP,X_REMAP,V_REMAP};
-
-// same as fix_wall.cpp
-
-enum{EDGE,CONSTANT,VARIABLE};
-
-/* ---------------------------------------------------------------------- */
-
-PairLubricate::PairLubricate(LAMMPS *lmp) : Pair(lmp)
-{
- single_enable = 0;
-
- // set comm size needed by this Pair
-
- comm_forward = 6;
-}
-
-/* ---------------------------------------------------------------------- */
-
-PairLubricate::~PairLubricate()
-{
- if (allocated) {
- memory->destroy(setflag);
- memory->destroy(cutsq);
-
- memory->destroy(cut);
- memory->destroy(cut_inner);
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairLubricate::compute(int eflag, int vflag)
-{
- int i,j,ii,jj,inum,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;
-
- if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = vflag_fdotr = 0;
-
- double **x = atom->x;
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int newton_pair = force->newton_pair;
-
- inum = list->inum;
- 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 = 0; ii < inum; 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
-
- comm->forward_comm_pair(this);
- }
-
- // This section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (int j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0);
- RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*
- (1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*
- (1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
-
- // end of R0 adjustment code
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- // 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]);
- }
- }
-
- if (!flagHI) continue;
-
- 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;
-
- // 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.0)*(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 = 0; ii < inum; 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];
- }
- }
-
- if (vflag_fdotr) virial_fdotr_compute();
-}
-
-/* ----------------------------------------------------------------------
- allocate all arrays
-------------------------------------------------------------------------- */
-
-void PairLubricate::allocate()
-{
- allocated = 1;
- int n = atom->ntypes;
-
- memory->create(setflag,n+1,n+1,"pair:setflag");
- for (int i = 1; i <= n; i++)
- for (int j = i; j <= n; j++)
- setflag[i][j] = 0;
-
- memory->create(cutsq,n+1,n+1,"pair:cutsq");
-
- memory->create(cut,n+1,n+1,"pair:cut");
- memory->create(cut_inner,n+1,n+1,"pair:cut_inner");
-}
-
-/* ----------------------------------------------------------------------
- global settings
-------------------------------------------------------------------------- */
-
-void PairLubricate::settings(int narg, char **arg)
-{
- if (narg != 5 && narg != 7) error->all(FLERR,"Illegal pair_style command");
-
- mu = force->numeric(FLERR,arg[0]);
- flaglog = force->inumeric(FLERR,arg[1]);
- flagfld = force->inumeric(FLERR,arg[2]);
- cut_inner_global = force->numeric(FLERR,arg[3]);
- cut_global = force->numeric(FLERR,arg[4]);
-
- flagHI = flagVF = 1;
- if (narg == 7) {
- flagHI = force->inumeric(FLERR,arg[5]);
- flagVF = force->inumeric(FLERR,arg[6]);
- }
-
- if (flaglog == 1 && flagHI == 0) {
- error->warning(FLERR,"Cannot include log terms without 1/r terms; "
- "setting flagHI to 1");
- flagHI = 1;
- }
-
- // reset cutoffs that have been explicitly set
-
- if (allocated) {
- for (int i = 1; i <= atom->ntypes; i++)
- for (int j = i+1; j <= atom->ntypes; j++)
- if (setflag[i][j]) {
- cut_inner[i][j] = cut_inner_global;
- cut[i][j] = cut_global;
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- set coeffs for one or more type pairs
-------------------------------------------------------------------------- */
-
-void PairLubricate::coeff(int narg, char **arg)
-{
- if (narg != 2 && narg != 4)
- error->all(FLERR,"Incorrect args for pair coefficients");
-
- if (!allocated) allocate();
-
- int ilo,ihi,jlo,jhi;
- force->bounds(arg[0],atom->ntypes,ilo,ihi);
- force->bounds(arg[1],atom->ntypes,jlo,jhi);
-
- double cut_inner_one = cut_inner_global;
- double cut_one = cut_global;
- if (narg == 4) {
- cut_inner_one = force->numeric(FLERR,arg[2]);
- cut_one = force->numeric(FLERR,arg[3]);
- }
-
- int count = 0;
- for (int i = ilo; i <= ihi; i++) {
- for (int j = MAX(jlo,i); j <= jhi; j++) {
- cut_inner[i][j] = cut_inner_one;
- cut[i][j] = cut_one;
- setflag[i][j] = 1;
- count++;
- }
- }
-
- if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
-}
-
-/* ----------------------------------------------------------------------
- init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairLubricate::init_style()
-{
- if (!atom->sphere_flag)
- error->all(FLERR,"Pair lubricate requires atom style sphere");
- if (comm->ghost_velocity == 0)
- error->all(FLERR,"Pair lubricate requires ghost atoms store velocity");
-
- neighbor->request(this,instance_me);
-
- // require that atom radii are identical within each type
- // require monodisperse system with same radii for all types
-
- double radtype;
- for (int i = 1; i <= atom->ntypes; i++) {
- if (!atom->radius_consistency(i,radtype))
- error->all(FLERR,"Pair lubricate requires monodisperse particles");
- if (i > 1 && radtype != rad)
- error->all(FLERR,"Pair lubricate requires monodisperse particles");
- rad = radtype;
- }
-
- // check for fix deform, if exists it must use "remap v"
- // If box will change volume, set appropriate flag so that volume
- // and v.f. corrections are re-calculated at every step.
- //
- // If available volume is different from box volume
- // due to walls, set volume appropriately; if walls will
- // move, set appropriate flag so that volume and v.f. corrections
- // are re-calculated at every step.
-
- shearing = flagdeform = flagwall = 0;
- for (int i = 0; i < modify->nfix; i++){
- if (strcmp(modify->fix[i]->style,"deform") == 0) {
- shearing = flagdeform = 1;
- if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP)
- error->all(FLERR,"Using pair lubricate with inconsistent "
- "fix deform remap option");
- }
- if (strstr(modify->fix[i]->style,"wall") != NULL) {
- if (flagwall)
- error->all(FLERR,
- "Cannot use multiple fix wall commands with pair lubricate");
- flagwall = 1; // Walls exist
- wallfix = (FixWall *) modify->fix[i];
- if (wallfix->xflag) flagwall = 2; // Moving walls exist
- }
- }
-
- // set the isotropic constants that depend on the volume fraction
- // vol_T = total volume
-
- double vol_T;
- double wallcoord;
- if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
- else {
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
-
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallfix->xindex[m] = input->variable->find(wallfix->xstr[m]);
- //Since fix->wall->init happens after pair->init_style
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
-
- else wallcoord = wallfix->coord0[m];
-
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
- (wallhi[2] - walllo[2]);
- }
-
- // vol_P = volume of particles, assuming monodispersity
- // vol_f = volume fraction
-
- vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3.0);
- double vol_f = vol_P/vol_T;
-
- if (!flagVF) vol_f = 0;
-
- // set isotropic constants for FLD
-
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0);
- RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
-
-
- // set Ef = 0 since used whether shearing or not
-
- Ef[0][0] = Ef[0][1] = Ef[0][2] = 0.0;
- Ef[1][0] = Ef[1][1] = Ef[1][2] = 0.0;
- Ef[2][0] = Ef[2][1] = Ef[2][2] = 0.0;
-}
-
-/* ----------------------------------------------------------------------
- init for one type pair i,j and corresponding j,i
-------------------------------------------------------------------------- */
-
-double PairLubricate::init_one(int i, int j)
-{
- if (setflag[i][j] == 0) {
- cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
- cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
- }
-
- cut_inner[j][i] = cut_inner[i][j];
-
- return cut[i][j];
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes to restart file
-------------------------------------------------------------------------- */
-
-void PairLubricate::write_restart(FILE *fp)
-{
- write_restart_settings(fp);
-
- int i,j;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i; j <= atom->ntypes; j++) {
- fwrite(&setflag[i][j],sizeof(int),1,fp);
- if (setflag[i][j]) {
- fwrite(&cut_inner[i][j],sizeof(double),1,fp);
- fwrite(&cut[i][j],sizeof(double),1,fp);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- proc 0 reads from restart file, bcasts
-------------------------------------------------------------------------- */
-
-void PairLubricate::read_restart(FILE *fp)
-{
- read_restart_settings(fp);
- allocate();
-
- int i,j;
- int me = comm->me;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i; j <= atom->ntypes; j++) {
- if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
- MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
- if (setflag[i][j]) {
- if (me == 0) {
- fread(&cut_inner[i][j],sizeof(double),1,fp);
- fread(&cut[i][j],sizeof(double),1,fp);
- }
- MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes to restart file
-------------------------------------------------------------------------- */
-
-void PairLubricate::write_restart_settings(FILE *fp)
-{
- fwrite(&mu,sizeof(double),1,fp);
- fwrite(&flaglog,sizeof(int),1,fp);
- fwrite(&flagfld,sizeof(int),1,fp);
- fwrite(&cut_inner_global,sizeof(double),1,fp);
- fwrite(&cut_global,sizeof(double),1,fp);
- fwrite(&offset_flag,sizeof(int),1,fp);
- fwrite(&mix_flag,sizeof(int),1,fp);
- fwrite(&flagHI,sizeof(int),1,fp);
- fwrite(&flagVF,sizeof(int),1,fp);
-}
-
-/* ----------------------------------------------------------------------
- proc 0 reads from restart file, bcasts
-------------------------------------------------------------------------- */
-
-void PairLubricate::read_restart_settings(FILE *fp)
-{
- int me = comm->me;
- if (me == 0) {
- fread(&mu,sizeof(double),1,fp);
- fread(&flaglog,sizeof(int),1,fp);
- fread(&flagfld,sizeof(int),1,fp);
- fread(&cut_inner_global,sizeof(double),1,fp);
- fread(&cut_global,sizeof(double),1,fp);
- fread(&offset_flag,sizeof(int),1,fp);
- fread(&mix_flag,sizeof(int),1,fp);
- fread(&flagHI,sizeof(int),1,fp);
- fread(&flagVF,sizeof(int),1,fp);
- }
- MPI_Bcast(&mu,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&flaglog,1,MPI_INT,0,world);
- MPI_Bcast(&flagfld,1,MPI_INT,0,world);
- MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
- MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
- MPI_Bcast(&flagHI,1,MPI_INT,0,world);
- MPI_Bcast(&flagVF,1,MPI_INT,0,world);
-}
-
-/* ---------------------------------------------------------------------- */
-
-int PairLubricate::pack_forward_comm(int n, int *list, double *buf,
- int pbc_flag, int *pbc)
-{
- int i,j,m;
-
- double **v = atom->v;
- double **omega = atom->omega;
-
- m = 0;
- for (i = 0; i < n; i++) {
- j = list[i];
- buf[m++] = v[j][0];
- buf[m++] = v[j][1];
- buf[m++] = v[j][2];
- buf[m++] = omega[j][0];
- buf[m++] = omega[j][1];
- buf[m++] = omega[j][2];
- }
-
- return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairLubricate::unpack_forward_comm(int n, int first, double *buf)
-{
- int i,m,last;
-
- double **v = atom->v;
- double **omega = atom->omega;
-
- m = 0;
- last = first + n;
- for (i = first; i < last; i++) {
- v[i][0] = buf[m++];
- v[i][1] = buf[m++];
- v[i][2] = buf[m++];
- omega[i][0] = buf[m++];
- omega[i][1] = buf[m++];
- omega[i][2] = buf[m++];
- }
-}
-
-/* ----------------------------------------------------------------------
- check if name is recognized, return integer index for that name
- if name not recognized, return -1
- if type pair setting, return -2 if no type pairs are set
-------------------------------------------------------------------------- */
-
-int PairLubricate::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi)
-{
- if (strcmp(name,"mu") == 0) return 0;
- return -1;
-}
-
-/* ----------------------------------------------------------------------
- adapt parameter indexed by which
- change all pair variables affected by the reset parameter
- if type pair setting, set I-J and J-I coeffs
-------------------------------------------------------------------------- */
-
-void PairLubricate::adapt(int which, int ilo, int ihi, int jlo, int jhi,
- double value)
-{
- mu = value;
-}
diff --git a/src/FLD/pair_lubricate.h b/src/FLD/pair_lubricate.h
deleted file mode 100644
index c2c4c6c417a941b081571dac9722004603af9bc4..0000000000000000000000000000000000000000
--- a/src/FLD/pair_lubricate.h
+++ /dev/null
@@ -1,103 +0,0 @@
-/* -*- 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 PAIR_CLASS
-
-PairStyle(lubricate,PairLubricate)
-
-#else
-
-#ifndef LMP_PAIR_LUBRICATE_H
-#define LMP_PAIR_LUBRICATE_H
-
-#include "pair.h"
-
-namespace LAMMPS_NS {
-
-class PairLubricate : public Pair {
- public:
- PairLubricate(class LAMMPS *);
- virtual ~PairLubricate();
- virtual void compute(int, int);
- void settings(int, char **);
- void coeff(int, char **);
- double init_one(int, int);
- virtual void init_style();
- void write_restart(FILE *);
- void read_restart(FILE *);
- void write_restart_settings(FILE *);
- void read_restart_settings(FILE *);
- int pre_adapt(char *, int, int, int, int);
- void adapt(int, int, int, int, int, double);
-
- int pack_forward_comm(int, int *, double *, int, int *);
- void unpack_forward_comm(int, int, double *);
-
- protected:
- double mu,cut_inner_global,cut_global;
- double rad;
- int flaglog,flagfld,shearing;
- int flagdeform, flagwall;
- double vol_P;
- class FixWall *wallfix;
- int flagVF, flagHI;
-
- double Ef[3][3];
- double R0,RT0,RS0;
- double **cut_inner,**cut;
-
- void allocate();
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Illegal ... command
-
-Self-explanatory. Check the input script syntax and compare to the
-documentation for the command. You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-W: Cannot include log terms without 1/r terms; setting flagHI to 1
-
-Self-explanatory.
-
-E: Incorrect args for pair coefficients
-
-Self-explanatory. Check the input script or data file.
-
-E: Pair lubricate requires atom style sphere
-
-Self-explanatory.
-
-E: Pair lubricate requires ghost atoms store velocity
-
-Use the comm_modify vel yes command to enable this.
-
-E: Pair lubricate requires monodisperse particles
-
-All particles must be the same finite size.
-
-E: Using pair lubricate with inconsistent fix deform remap option
-
-Must use remap v option with fix deform with this pair style.
-
-E: Cannot use multiple fix wall commands with pair lubricate
-
-Self-explanatory.
-
-*/
diff --git a/src/FLD/pair_lubricateU.cpp b/src/FLD/pair_lubricateU.cpp
deleted file mode 100644
index 214bc3c2f3ed68e71fb87d2d234087582d6165f0..0000000000000000000000000000000000000000
--- a/src/FLD/pair_lubricateU.cpp
+++ /dev/null
@@ -1,2061 +0,0 @@
-/* ----------------------------------------------------------------------
- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
- http://lammps.sandia.gov, Sandia National Laboratories
- Steve Plimpton, sjplimp@sandia.gov
-
- Copyright (2003) Sandia Corporation. Under the terms of Contract
- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
- certain rights in this software. This software is distributed under
- the GNU General Public License.
-
- See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing authors: Amit Kumar and Michael Bybee (UIUC)
-------------------------------------------------------------------------- */
-
-#include <mpi.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include "pair_lubricateU.h"
-#include "atom.h"
-#include "atom_vec.h"
-#include "comm.h"
-#include "force.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
-#include "domain.h"
-#include "update.h"
-#include "math_const.h"
-#include "modify.h"
-#include "fix.h"
-#include "fix_deform.h"
-#include "fix_wall.h"
-#include "input.h"
-#include "variable.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-using namespace MathConst;
-
-#define TOL 1E-4 // tolerance for conjugate gradient
-
-// same as fix_wall.cpp
-
-enum{EDGE,CONSTANT,VARIABLE};
-
-/* ---------------------------------------------------------------------- */
-
-PairLubricateU::PairLubricateU(LAMMPS *lmp) : Pair(lmp)
-{
- single_enable = 0;
-
- // pair lubricateU cannot compute virial as F dot r
- // due to how drag forces are applied to atoms
- // correct method is how per-atom virial does it
-
- no_virial_fdotr_compute = 1;
-
- nmax = 0;
- fl = Tl = xl = NULL;
-
- cgmax = 0;
- bcg = xcg = rcg = rcg1 = pcg = RU = NULL;
-
- // set comm size needed by this Pair
-
- comm_forward = 6;
-}
-
-/* ---------------------------------------------------------------------- */
-
-PairLubricateU::~PairLubricateU()
-{
- memory->destroy(fl);
- memory->destroy(Tl);
- memory->destroy(xl);
-
- memory->destroy(bcg);
- memory->destroy(xcg);
- memory->destroy(rcg);
- memory->destroy(rcg1);
- memory->destroy(pcg);
- memory->destroy(RU);
-
- if (allocated) {
- memory->destroy(setflag);
- memory->destroy(cutsq);
-
- memory->destroy(cut);
- memory->destroy(cut_inner);
- }
-}
-
-/* ----------------------------------------------------------------------
- It first has to solve for the velocity of the particles such that
- the net force on the particles is zero. NOTE: it has to be the last
- type of pair interaction specified in the input file. Also, it
- assumes that no other types of interactions, like k-space, is
- present. As already mentioned, the net force on the particles after
- this pair interaction would be identically zero.
- ---------------------------------------------------------------------- */
-
-void PairLubricateU::compute(int eflag, int vflag)
-{
- int i,j;
-
- double **x = atom->x;
- double **f = atom->f;
- double **torque = atom->torque;
- int nlocal = atom->nlocal;
- int nghost = atom->nghost;
- int nall = nlocal + nghost;
-
- if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = vflag_fdotr = 0;
-
- // skip compute() if called from integrate::setup()
- // this is b/c do not want compute() to update velocities twice on a restart
- // when restarting, call compute on step N (last step of prev run),
- // again on step N (setup of restart run),
- // then on step N+1 (first step of restart)
- // so this is one extra time which leads to bad dynamics
-
- if (update->setupflag) return;
-
- // grow per-atom arrays if necessary
- // need to be atom->nmax in length
-
- if (atom->nmax > nmax) {
- memory->destroy(fl);
- memory->destroy(Tl);
- memory->destroy(xl);
- nmax = atom->nmax;
- memory->create(fl,nmax,3,"pair:fl");
- memory->create(Tl,nmax,3,"pair:Tl");
- memory->create(xl,nmax,3,"pair:xl");
- }
-
- // Added to implement midpoint integration scheme
- // Save force, torque found so far. Also save the positions
-
- for (i=0;i<nlocal+nghost;i++) {
- for (j=0;j<3;j++) {
- fl[i][j] = f[i][j];
- Tl[i][j] = torque[i][j];
- xl[i][j] = x[i][j];
- }
- }
-
- // Stage one of Midpoint method
- // Solve for velocities based on intial positions
-
- stage_one();
-
- // find positions at half the timestep and store in xl
-
- intermediates(nall,xl);
-
- // store back the saved forces and torques in original arrays
-
- for(i=0;i<nlocal+nghost;i++) {
- for(j=0;j<3;j++) {
- f[i][j] = fl[i][j];
- torque[i][j] = Tl[i][j];
- }
- }
-
- // stage two: this will give the final velocities
-
- stage_two(xl);
-}
-
-/* ------------------------------------------------------------------------
- Stage one of midpoint method
-------------------------------------------------------------------------- */
-
-void PairLubricateU::stage_one()
-{
- int i,j,ii,inum;
-
- double **x = atom->x;
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
-
- int newton_pair = force->newton_pair;
- int *ilist;
-
- inum = list->inum;
- ilist = list->ilist;
-
- if (6*inum > cgmax) {
- memory->destroy(bcg);
- memory->destroy(xcg);
- memory->destroy(rcg);
- memory->destroy(rcg1);
- memory->destroy(pcg);
- memory->destroy(RU);
- cgmax = 6*inum;
- memory->create(bcg,cgmax,"pair:bcg");
- memory->create(xcg,cgmax,"pair:bcg");
- memory->create(rcg,cgmax,"pair:bcg");
- memory->create(rcg1,cgmax,"pair:bcg");
- memory->create(pcg,cgmax,"pair:bcg");
- memory->create(RU,cgmax,"pair:bcg");
- }
-
- double alpha,beta;
- double normi,error,normig;
- double send[2],recv[2],rcg_dot_rcg;
-
- // First compute R_FE*E
-
- compute_RE();
-
- // Reverse communication of forces and torques to
- // accumulate the net force on each of the particles
-
- if (newton_pair) comm->reverse_comm();
-
- // CONJUGATE GRADIENT
- // Find the right hand side= -ve of all forces/torques
- // b = 6*Npart in overall size
-
- for(ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- for (j = 0; j < 3; j++) {
- bcg[6*ii+j] = -f[i][j];
- bcg[6*ii+j+3] = -torque[i][j];
- }
- }
-
- // Start solving the equation : F^H = -F^P -F^B - F^H_{Ef}
- // Store initial guess for velocity and angular-velocities/angular momentum
- // NOTE velocities and angular velocities are assumed relative to the fluid
-
- for (ii=0;ii<inum;ii++)
- for (j=0;j<3;j++) {
- xcg[6*ii+j] = 0.0;
- xcg[6*ii+j+3] = 0.0;
- }
-
- // Copy initial guess to the global arrays to be acted upon by R_{FU}
- // and returned by f and torque arrays
-
- copy_vec_uo(inum,xcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- // Find initial residual
-
- compute_RU();
-
- // reverse communication of forces and torques
-
- if (newton_pair) comm->reverse_comm();
-
- copy_uo_vec(inum,f,torque,RU);
-
- for (i=0;i<6*inum;i++)
- rcg[i] = bcg[i] - RU[i];
-
- // Set initial conjugate direction
-
- for (i=0;i<6*inum;i++)
- pcg[i] = rcg[i];
-
- // Find initial norm of the residual or norm of the RHS (either is fine)
-
- normi = dot_vec_vec(6*inum,bcg,bcg);
-
- MPI_Allreduce(&normi,&normig,1,MPI_DOUBLE,MPI_SUM,world);
-
- // Loop until convergence
-
- do {
- // find R*p
-
- copy_vec_uo(inum,pcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- compute_RU();
-
- // reverse communication of forces and torques
-
- if (newton_pair) comm->reverse_comm();
-
-
- copy_uo_vec(inum,f,torque,RU);
-
- // Find alpha
-
- send[0] = dot_vec_vec(6*inum,rcg,rcg);
- send[1] = dot_vec_vec(6*inum,RU,pcg);
-
- MPI_Allreduce(send,recv,2,MPI_DOUBLE,MPI_SUM,world);
-
- alpha = recv[0]/recv[1];
- rcg_dot_rcg = recv[0];
-
- // Find new x
-
- for (i=0;i<6*inum;i++)
- xcg[i] = xcg[i] + alpha*pcg[i];
-
- // find new residual
-
- for (i=0;i<6*inum;i++)
- rcg1[i] = rcg[i] - alpha*RU[i];
-
- // find beta
-
- send[0] = dot_vec_vec(6*inum,rcg1,rcg1);
-
- MPI_Allreduce(send,recv,1,MPI_DOUBLE,MPI_SUM,world);
-
- beta = recv[0]/rcg_dot_rcg;
-
- // Find new conjugate direction
-
- for (i=0;i<6*inum;i++)
- pcg[i] = rcg1[i] + beta*pcg[i];
-
- for (i=0;i<6*inum;i++)
- rcg[i] = rcg1[i];
-
- // Find relative error
-
- error = sqrt(recv[0]/normig);
-
- } while (error > TOL);
-
- // update the final converged velocities in respective arrays
-
- copy_vec_uo(inum,xcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- // Find actual particle's velocities from relative velocities
- // Only non-zero component of fluid's vel : vx=gdot*y and wz=-gdot/2
-
- for (ii=0;ii<inum;ii++) {
- i = ilist[ii];
- v[i][0] = v[i][0] + gdot*x[i][1];
- omega[i][2] = omega[i][2] - gdot/2.0;
- }
-}
-
-/*---------------------------------------------------------------
- Finds the position of the particles at half the time step
-----------------------------------------------------------------*/
-
-void PairLubricateU::intermediates(int nall, double **xl)
-{
- int i;
- double **x = atom->x;
- double **v = atom->v;
- double dtv = update->dt;
-
- for (i=0;i<nall;i++) {
- xl[i][0] = x[i][0] + 0.5*dtv*v[i][0];
- xl[i][1] = x[i][1] + 0.5*dtv*v[i][1];
- xl[i][2] = x[i][2] + 0.5*dtv*v[i][2];
- }
-}
-
-/* ------------------------------------------------------------------------
- Stage one of midpoint method
-------------------------------------------------------------------------- */
-
-void PairLubricateU::stage_two(double **x)
-{
- int i,j,ii,inum;
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
-
- int newton_pair = force->newton_pair;
- int *ilist;
-
- inum = list->inum;
- ilist = list->ilist;
-
- double alpha,beta;
- double normi,error,normig;
- double send[2],recv[2],rcg_dot_rcg;
-
- // First compute R_FE*E
-
- compute_RE(x);
-
- // Reverse communication of forces and torques to
- // accumulate the net force on each of the particles
-
- if (newton_pair) comm->reverse_comm();
-
- // CONJUGATE GRADIENT
- // Find the right hand side= -ve of all forces/torques
- // b = 6*Npart in overall size
-
- for(ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- for (j = 0; j < 3; j++) {
- bcg[6*ii+j] = -f[i][j];
- bcg[6*ii+j+3] = -torque[i][j];
- }
- }
-
- // Start solving the equation : F^H = -F^P -F^B - F^H_{Ef}
- // Store initial guess for velocity and angular-velocities/angular momentum
- // NOTE velocities and angular velocities are assumed relative to the fluid
-
- for (ii=0;ii<inum;ii++)
- for (j=0;j<3;j++) {
- xcg[6*ii+j] = 0.0;
- xcg[6*ii+j+3] = 0.0;
- }
-
- // Copy initial guess to the global arrays to be acted upon by R_{FU}
- // and returned by f and torque arrays
-
- copy_vec_uo(inum,xcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- // Find initial residual
-
- compute_RU(x);
-
- // reverse communication of forces and torques
-
- if (newton_pair) comm->reverse_comm();
-
- copy_uo_vec(inum,f,torque,RU);
-
- for (i=0;i<6*inum;i++)
- rcg[i] = bcg[i] - RU[i];
-
- // Set initial conjugate direction
-
- for (i=0;i<6*inum;i++)
- pcg[i] = rcg[i];
-
- // Find initial norm of the residual or norm of the RHS (either is fine)
-
- normi = dot_vec_vec(6*inum,bcg,bcg);
-
- MPI_Allreduce(&normi,&normig,1,MPI_DOUBLE,MPI_SUM,world);
-
- // Loop until convergence
-
- do {
- // find R*p
-
- copy_vec_uo(inum,pcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- compute_RU(x);
-
- // reverse communication of forces and torques
-
- if (newton_pair) comm->reverse_comm();
-
- copy_uo_vec(inum,f,torque,RU);
-
- // Find alpha
-
- send[0] = dot_vec_vec(6*inum,rcg,rcg);
- send[1] = dot_vec_vec(6*inum,RU,pcg);
-
- MPI_Allreduce(send,recv,2,MPI_DOUBLE,MPI_SUM,world);
-
- alpha = recv[0]/recv[1];
- rcg_dot_rcg = recv[0];
-
- // Find new x
-
- for (i=0;i<6*inum;i++)
- xcg[i] = xcg[i] + alpha*pcg[i];
-
- // find new residual
-
- for (i=0;i<6*inum;i++)
- rcg1[i] = rcg[i] - alpha*RU[i];
-
- // find beta
-
- send[0] = dot_vec_vec(6*inum,rcg1,rcg1);
-
- MPI_Allreduce(send,recv,1,MPI_DOUBLE,MPI_SUM,world);
-
- beta = recv[0]/rcg_dot_rcg;
-
- // Find new conjugate direction
-
- for (i=0;i<6*inum;i++)
- pcg[i] = rcg1[i] + beta*pcg[i];
-
- for (i=0;i<6*inum;i++)
- rcg[i] = rcg1[i];
-
- // Find relative error
-
- error = sqrt(recv[0]/normig);
-
- } while (error > TOL);
-
-
- // update the final converged velocities in respective arrays
-
- copy_vec_uo(inum,xcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- // Compute the viscosity/pressure
-
- if (evflag) compute_Fh(x);
-
- // Find actual particle's velocities from relative velocities
- // Only non-zero component of fluid's vel : vx=gdot*y and wz=-gdot/2
-
- for (ii=0;ii<inum;ii++) {
- i = ilist[ii];
- v[i][0] = v[i][0] + gdot*x[i][1];
- omega[i][2] = omega[i][2] - gdot/2.0;
- }
-}
-
-/* ------------------------------------------------------------------------
- This function computes the final hydrodynamic force once the
- velocities have converged.
- ------------------------------------------------------------------------- */
-
-void PairLubricateU::compute_Fh(double **x)
-{
- int i,j,ii,jj,inum,jnum,itype,jtype;
- double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
- double rsq,r,h_sep;
- double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
- double vt1,vt2,vt3;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int nghost = atom->nghost;
- int newton_pair = force->newton_pair;
-
- double radi;
-
- double vxmu2f = force->vxmu2f;
- double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- // This section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (int j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- // R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- // RT0 = 8*MY_PI*mu*pow(rad,3);
- RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*
- (1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- // R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- // RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*
- (1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
-
-
- // end of R0 adjustment code
-
- // Set force to zero which is the final value after this pair interaction
- for (i=0;i<nlocal+nghost;i++)
- for (j=0;j<3;j++) {
- f[i][j] = 0.0;
- torque[i][j] = 0.0;
- }
-
- // reverse communication of forces and torques
-
- if (newton_pair) comm->reverse_comm(); // not really needed
-
- // Find additional contribution from the stresslets
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- // Find the contribution to stress from isotropic RS0
- // Set psuedo force to obtain the required contribution
- // need to set delx and fy only
-
- fx = 0.0; delx = radi;
- fy = vxmu2f*RS0*gdot/2.0/radi; dely = 0.0;
- fz = 0.0; delz = 0.0;
- if (evflag)
- ev_tally_xyz(i,i,nlocal,newton_pair,0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
-
- // Find angular velocity
-
- wi[0] = omega[i][0];
- wi[1] = omega[i][1];
- wi[2] = omega[i][2];
-
- if (!flagHI) continue;
-
- 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);
-
- // Use omega directly if it exists, else angmom
- // 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];
-
- // loc of the point of closest approach on particle i from its cente
-
- 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
-
- // particle i
-
- vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]);
- vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]);
- vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]);
-
- // particle j
-
- vj[0] = v[j][0] - (wj[1]*xl[2] - wj[2]*xl[1]);
- vj[1] = v[j][1] - (wj[2]*xl[0] - wj[0]*xl[2]);
- vj[2] = v[j][2] - (wj[0]*xl[1] - wj[1]*xl[0]);
-
-
- // Relative velocity at the point of closest approach
- // include contribution from Einf of the fluid
-
- vr1 = vi[0] - vj[0] -
- 2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]);
- vr2 = vi[1] - vj[1] -
- 2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]);
- vr3 = vi[2] - vj[2] -
- 2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[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;
-
- // Find the scalar resistances a_sq, a_sh and a_pu
-
- h_sep = r - 2.0*radi;
-
- // 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));
- } else
- a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep);
-
- // Find force due to squeeze type motion
-
- fx = a_sq*vn1;
- fy = a_sq*vn2;
- fz = a_sq*vn3;
-
- // Find 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 to obtain in appropriate units
-
- fx = vxmu2f*fx;
- fy = vxmu2f*fy;
- fz = vxmu2f*fz;
-
- if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
- 0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
- }
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- computes R_FU * U
----------------------------------------------------------------------- */
-
-void PairLubricateU::compute_RU()
-{
- int i,j,ii,jj,inum,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,wdotn,wt1,wt2,wt3;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- double **x = atom->x;
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int nghost = atom->nghost;
- int newton_pair = force->newton_pair;
-
- double vxmu2f = force->vxmu2f;
- double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- // This section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (int j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0);
- // RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- // RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
-
- // end of R0 adjustment code
-
- // Initialize f to zero
- for (i=0;i<nlocal+nghost;i++)
- for (j=0;j<3;j++) {
- f[i][j] = 0.0;
- torque[i][j] = 0.0;
- }
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- // Find angular velocity
-
- wi[0] = omega[i][0];
- wi[1] = omega[i][1];
- wi[2] = omega[i][2];
-
- // Contribution due to the isotropic terms
-
- 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 (!flagHI) continue;
-
- 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);
-
- // Use omega directly if it exists, else angmom
- // 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];
-
- // loc of the 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
-
- // particle i
-
- vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]);
- vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]);
- vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]);
-
- // particle j
-
- vj[0] = v[j][0] - (wj[1]*xl[2] - wj[2]*xl[1]);
- vj[1] = v[j][1] - (wj[2]*xl[0] - wj[0]*xl[2]);
- vj[2] = v[j][2] - (wj[0]*xl[1] - wj[1]*xl[0]);
-
- // Find the scalar resistances a_sq and a_sh
-
- h_sep = r - 2.0*radi;
-
- // 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.0)*(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
-
- 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;
-
- // Find force due to squeeze type motion
-
- fx = a_sq*vn1;
- fy = a_sq*vn2;
- fz = a_sq*vn3;
-
- // Find 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 to obtain in appropriate units
-
- fx = vxmu2f*fx;
- fy = vxmu2f*fy;
- fz = vxmu2f*fz;
-
- // Add to the total forc
-
- 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;
- }
-
- // Find 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;
-
- // Why a scale factor ?
-
- 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;
-
- // add to total
-
- 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;
- }
- }
- }
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- computes R_FU * U
----------------------------------------------------------------------- */
-
-void PairLubricateU::compute_RU(double **x)
-{
- int i,j,ii,jj,inum,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,wdotn,wt1,wt2,wt3;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int nghost = atom->nghost;
- int newton_pair = force->newton_pair;
-
- double vxmu2f = force->vxmu2f;
- double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- // This section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (int j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0);
- // RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- // RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
-
- // end of R0 adjustment code
-
- // Initialize f to zero
- for (i=0;i<nlocal+nghost;i++)
- for (j=0;j<3;j++) {
- f[i][j] = 0.0;
- torque[i][j] = 0.0;
- }
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- // Find angular velocity
-
- wi[0] = omega[i][0];
- wi[1] = omega[i][1];
- wi[2] = omega[i][2];
-
- // Contribution due to the isotropic terms
-
- 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 (!flagHI) continue;
-
- 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);
-
- // Use omega directly if it exists, else angmom
- // 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];
-
- // loc of the 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
-
- // particle i
-
- vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]);
- vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]);
- vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]);
-
- // particle j
-
- vj[0] = v[j][0] - (wj[1]*xl[2] - wj[2]*xl[1]);
- vj[1] = v[j][1] - (wj[2]*xl[0] - wj[0]*xl[2]);
- vj[2] = v[j][2] - (wj[0]*xl[1] - wj[1]*xl[0]);
-
- // Find the scalar resistances a_sq and a_sh
-
- h_sep = r - 2.0*radi;
-
- // 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.0)*(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
-
- 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;
-
- // Find force due to squeeze type motion
-
- fx = a_sq*vn1;
- fy = a_sq*vn2;
- fz = a_sq*vn3;
-
- // Find 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 to obtain in appropriate units
-
- fx = vxmu2f*fx;
- fy = vxmu2f*fy;
- fz = vxmu2f*fz;
-
- // Add to the 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;
- }
-
- // Find 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;
-
- // Why a scale factor ?
-
- 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;
-
- // add to total
-
- 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;
- }
- }
- }
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- This computes R_{FE}*E , where E is the rate of strain of tensor which is
- known apriori, as it depends only on the known fluid velocity.
- So, this part of the hydrodynamic interaction can be pre computed and
- transferred to the RHS
- ---------------------------------------------------------------------- */
-
-void PairLubricateU::compute_RE()
-{
- int i,j,ii,jj,inum,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;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- double **x = atom->x;
- double **f = atom->f;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int newton_pair = force->newton_pair;
-
- double vxmu2f = force->vxmu2f;
- double xl[3],a_sq,a_sh;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- if (!flagHI) return;
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- // No contribution from isotropic terms due to E
-
- 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);
-
- // loc of the 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;
-
- // Find the scalar resistances a_sq and a_sh
-
- h_sep = r - 2.0*radi;
-
- // 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 resistance for Squeeze type motions
-
- if (flaglog)
- a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1/h_sep));
- else
- a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep);
-
- // Scalar resistance for Shear type motions
-
- if (flaglog) {
- a_sh = 6*MY_PI*mu*radi*(1.0/6.0*log(1/h_sep));
- }
-
- // Relative velocity at the point of closest approach due to Ef only
-
- vr1 = -2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]);
- vr2 = -2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]);
- vr3 = -2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[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;
-
- // Find force due to squeeze type motion
-
- fx = a_sq*vn1;
- fy = a_sq*vn2;
- fz = a_sq*vn3;
-
- // Find 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 to obtain in appropriate units
-
- fx = vxmu2f*fx;
- fy = vxmu2f*fy;
- fz = vxmu2f*fz;
-
- // Add to the total forc
-
- 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;
- }
-
- // Find 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;
-
- // Why a scale factor ?
-
- 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;
- }
- }
- }
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- This computes R_{FE}*E , where E is the rate of strain of tensor which is
- known apriori, as it depends only on the known fluid velocity.
- So, this part of the hydrodynamic interaction can be pre computed and
- transferred to the RHS
- ---------------------------------------------------------------------- */
-
-void PairLubricateU::compute_RE(double **x)
-{
- int i,j,ii,jj,inum,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;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- double **f = atom->f;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int newton_pair = force->newton_pair;
-
- double vxmu2f = force->vxmu2f;
- double xl[3],a_sq,a_sh;
-
- if (!flagHI) return;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- // No contribution from isotropic terms due to E
-
- 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);
-
- // loc of the 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;
-
- // Find the scalar resistances a_sq and a_sh
-
- h_sep = r - 2.0*radi;
-
- // 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 resistance for Squeeze type motions
-
- if (flaglog)
- a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1/h_sep));
- else
- a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep);
-
- // Scalar resistance for Shear type motions
-
- if (flaglog) {
- a_sh = 6*MY_PI*mu*radi*(1.0/6.0*log(1/h_sep));
- }
-
- // Relative velocity at the point of closest approach due to Ef only
-
- vr1 = -2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]);
- vr2 = -2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]);
- vr3 = -2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[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;
-
- // Find force due to squeeze type motion
-
- fx = a_sq*vn1;
- fy = a_sq*vn2;
- fz = a_sq*vn3;
-
- // Find 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 to obtain in appropriate units
-
- fx = vxmu2f*fx;
- fy = vxmu2f*fy;
- fz = vxmu2f*fz;
-
- // Add to the total forc
-
- 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;
- }
-
- // Find 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;
-
- // Why a scale factor ?
-
- 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;
- }
- }
- }
- }
- }
-}
-
-
-/* ----------------------------------------------------------------------
- allocate all arrays
-------------------------------------------------------------------------- */
-
-void PairLubricateU::allocate()
-{
- allocated = 1;
- int n = atom->ntypes;
-
- setflag = 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;
-
- cutsq = memory->create(cutsq,n+1,n+1,"pair:cutsq");
-
- memory->create(cut,n+1,n+1,"pair:cut");
- memory->create(cut_inner,n+1,n+1,"pair:cut_inner");
-}
-
-/*-----------------------------------------------------------------------
- global settings
-------------------------------------------------------------------------- */
-
-void PairLubricateU::settings(int narg, char **arg)
-{
- if (narg != 5 && narg != 7) error->all(FLERR,"Illegal pair_style command");
-
- mu = force->numeric(FLERR,arg[0]);
- flaglog = force->inumeric(FLERR,arg[1]);
- cut_inner_global = force->numeric(FLERR,arg[2]);
- cut_global = force->numeric(FLERR,arg[3]);
- gdot = force->numeric(FLERR,arg[4]);
-
- flagHI = flagVF = 1;
- if (narg == 7) {
- flagHI = force->inumeric(FLERR,arg[5]);
- flagVF = force->inumeric(FLERR,arg[6]);
- }
-
- if (flaglog == 1 && flagHI == 0) {
- error->warning(FLERR,"Cannot include log terms without 1/r terms; "
- "setting flagHI to 1.");
- flagHI = 1;
- }
-
- // reset cutoffs that have been explicitly set
-
- if (allocated) {
- int i,j;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i+1; j <= atom->ntypes; j++)
- if (setflag[i][j]) {
- cut_inner[i][j] = cut_inner_global;
- cut[i][j] = cut_global;
- }
- }
-
- // store the rate of strain tensor
-
- Ef[0][0] = 0.0;
- Ef[0][1] = 0.5*gdot;
- Ef[0][2] = 0.0;
- Ef[1][0] = 0.5*gdot;
- Ef[1][1] = 0.0;
- Ef[1][2] = 0.0;
- Ef[2][0] = 0.0;
- Ef[2][1] = 0.0;
- Ef[2][2] = 0.0;
-}
-
-/*-----------------------------------------------------------------------
- set coeffs for one or more type pairs
-------------------------------------------------------------------------- */
-
-void PairLubricateU::coeff(int narg, char **arg)
-{
- if (narg != 2 && narg != 4)
- error->all(FLERR,"Incorrect args for pair coefficients");
-
- if (!allocated) allocate();
-
- int ilo,ihi,jlo,jhi;
- force->bounds(arg[0],atom->ntypes,ilo,ihi);
- force->bounds(arg[1],atom->ntypes,jlo,jhi);
-
- double cut_inner_one = cut_inner_global;
- double cut_one = cut_global;
- if (narg == 4) {
- cut_inner_one = force->numeric(FLERR,arg[2]);
- cut_one = force->numeric(FLERR,arg[3]);
- }
-
- int count = 0;
- for (int i = ilo; i <= ihi; i++) {
- for (int j = MAX(jlo,i); j <= jhi; j++) {
- cut_inner[i][j] = cut_inner_one;
- cut[i][j] = cut_one;
- setflag[i][j] = 1;
- count++;
- }
- }
-
- if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
-}
-
-/* ----------------------------------------------------------------------
- init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairLubricateU::init_style()
-{
- if (!atom->sphere_flag)
- error->all(FLERR,"Pair lubricateU requires atom style sphere");
- if (comm->ghost_velocity == 0)
- error->all(FLERR,"Pair lubricateU requires ghost atoms store velocity");
-
- neighbor->request(this,instance_me);
-
- // require that atom radii are identical within each type
- // require monodisperse system with same radii for all types
-
- double radtype;
- for (int i = 1; i <= atom->ntypes; i++) {
- if (!atom->radius_consistency(i,radtype))
- error->all(FLERR,"Pair lubricateU requires monodisperse particles");
- if (i > 1 && radtype != rad)
- error->all(FLERR,"Pair lubricateU requires monodisperse particles");
- }
-
- // check for fix deform, if exists it must use "remap v"
- // If box will change volume, set appropriate flag so that volume
- // and v.f. corrections are re-calculated at every step.
- //
- // If available volume is different from box volume
- // due to walls, set volume appropriately; if walls will
- // move, set appropriate flag so that volume and v.f. corrections
- // are re-calculated at every step.
-
- flagdeform = flagwall = 0;
- for (int i = 0; i < modify->nfix; i++){
- if (strcmp(modify->fix[i]->style,"deform") == 0)
- flagdeform = 1;
- else if (strstr(modify->fix[i]->style,"wall") != NULL) {
- if (flagwall)
- error->all(FLERR,
- "Cannot use multiple fix wall commands with "
- "pair lubricateU");
- flagwall = 1; // Walls exist
- wallfix = (FixWall *) modify->fix[i];
- if (wallfix->xflag) flagwall = 2; // Moving walls exist
- }
- }
-
- // set the isotropic constants depending on the volume fraction
- // vol_T = total volumeshearing = flagdeform = flagwall = 0;
- double vol_T, wallcoord;
- if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
- else {
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallfix->xindex[m] = input->variable->find(wallfix->xstr[m]);
- //Since fix->wall->init happens after pair->init_style
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
-
- else wallcoord = wallfix->coord0[m];
-
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
- (wallhi[2] - walllo[2]);
- }
-
-
- // assuming monodisperse spheres, vol_P = volume of the particles
-
- double tmp = 0.0;
- if (atom->radius) tmp = atom->radius[0];
- MPI_Allreduce(&tmp,&rad,1,MPI_DOUBLE,MPI_MAX,world);
-
- vol_P = atom->natoms * (4.0/3.0)*MY_PI*pow(rad,3.0);
-
- // vol_f = volume fraction
-
- double vol_f = vol_P/vol_T;
-
- if (!flagVF) vol_f = 0;
-
- // set the isotropic constant
-
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0); // not actually needed
- RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*pow(rad,3.0)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3.0)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
-}
-
-/* ----------------------------------------------------------------------
- init for one type pair i,j and corresponding j,i
-------------------------------------------------------------------------- */
-
-double PairLubricateU::init_one(int i, int j)
-{
- if (setflag[i][j] == 0) {
- cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
- cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
- }
-
- cut_inner[j][i] = cut_inner[i][j];
-
- return cut[i][j];
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes to restart file
-------------------------------------------------------------------------- */
-
-void PairLubricateU::write_restart(FILE *fp)
-{
- write_restart_settings(fp);
-
- int i,j;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i; j <= atom->ntypes; j++) {
- fwrite(&setflag[i][j],sizeof(int),1,fp);
- if (setflag[i][j]) {
- fwrite(&cut_inner[i][j],sizeof(double),1,fp);
- fwrite(&cut[i][j],sizeof(double),1,fp);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- proc 0 reads from restart file, bcasts
-------------------------------------------------------------------------- */
-
-void PairLubricateU::read_restart(FILE *fp)
-{
- read_restart_settings(fp);
- allocate();
-
- int i,j;
- int me = comm->me;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i; j <= atom->ntypes; j++) {
- if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
- MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
- if (setflag[i][j]) {
- if (me == 0) {
- fread(&cut_inner[i][j],sizeof(double),1,fp);
- fread(&cut[i][j],sizeof(double),1,fp);
- }
- MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- proc 0 writes to restart file
-------------------------------------------------------------------------- */
-
-void PairLubricateU::write_restart_settings(FILE *fp)
-{
- fwrite(&mu,sizeof(double),1,fp);
- fwrite(&flaglog,sizeof(int),1,fp);
- fwrite(&cut_inner_global,sizeof(double),1,fp);
- fwrite(&cut_global,sizeof(double),1,fp);
- fwrite(&offset_flag,sizeof(int),1,fp);
- fwrite(&mix_flag,sizeof(int),1,fp);
- fwrite(&flagHI,sizeof(int),1,fp);
- fwrite(&flagVF,sizeof(int),1,fp);
-}
-
-/* ----------------------------------------------------------------------
- proc 0 reads from restart file, bcasts
-------------------------------------------------------------------------- */
-
-void PairLubricateU::read_restart_settings(FILE *fp)
-{
- int me = comm->me;
- if (me == 0) {
- fread(&mu,sizeof(double),1,fp);
- fread(&flaglog,sizeof(int),1,fp);
- fread(&cut_inner_global,sizeof(double),1,fp);
- fread(&cut_global,sizeof(double),1,fp);
- fread(&offset_flag,sizeof(int),1,fp);
- fread(&mix_flag,sizeof(int),1,fp);
- fread(&flagHI,sizeof(int),1,fp);
- fread(&flagVF,sizeof(int),1,fp);
- }
- MPI_Bcast(&mu,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&flaglog,1,MPI_INT,0,world);
- MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
- MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
- MPI_Bcast(&flagHI,1,MPI_INT,0,world);
- MPI_Bcast(&flagVF,1,MPI_INT,0,world);
-}
-
-/*---------------------------------------------------------------------------*/
-
-void PairLubricateU::copy_vec_uo(int inum, double *xcg,
- double **v, double **omega)
-{
- int i,j,ii;
- int *ilist = list->ilist;
-
- for (ii=0;ii<inum;ii++) {
- i = ilist[ii];
-
- for (j=0;j<3;j++) {
- v[i][j] = xcg[6*ii+j];
- omega[i][j] = xcg[6*ii+j+3];
- }
- }
-}
-
-/*---------------------------------------------------------------------------*/
-
-void PairLubricateU::copy_uo_vec(int inum, double **f, double **torque,
- double *RU)
-{
- int i,j,ii;
- int *ilist;
-
- ilist = list->ilist;
-
- for (ii=0;ii<inum;ii++) {
- i = ilist[ii];
- for (j=0;j<3;j++) {
- RU[6*ii+j] = f[i][j];
- RU[6*ii+j+3] = torque[i][j];
- }
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int PairLubricateU::pack_forward_comm(int n, int *list, double *buf,
- int pbc_flag, int *pbc)
-{
- int i,j,m;
-
- double **v = atom->v;
- double **omega = atom->omega;
-
- m = 0;
- for (i = 0; i < n; i++) {
- j = list[i];
- buf[m++] = v[j][0];
- buf[m++] = v[j][1];
- buf[m++] = v[j][2];
- buf[m++] = omega[j][0];
- buf[m++] = omega[j][1];
- buf[m++] = omega[j][2];
- }
- return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairLubricateU::unpack_forward_comm(int n, int first, double *buf)
-{
- int i,m,last;
-
- double **v = atom->v;
- double **omega = atom->omega;
-
- m = 0;
- last = first + n;
- for (i = first; i < last; i++) {
- v[i][0] = buf[m++];
- v[i][1] = buf[m++];
- v[i][2] = buf[m++];
- omega[i][0] = buf[m++];
- omega[i][1] = buf[m++];
- omega[i][2] = buf[m++];
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairLubricateU::dot_vec_vec(int N, double *x, double *y)
-{
- double dotp=0.0;
- for (int i = 0; i < N; i++) dotp += x[i]*y[i];
- return dotp;
-}
diff --git a/src/FLD/pair_lubricateU.h b/src/FLD/pair_lubricateU.h
deleted file mode 100644
index 4968c665a61d54290c6cd0854dde3a82a428dd41..0000000000000000000000000000000000000000
--- a/src/FLD/pair_lubricateU.h
+++ /dev/null
@@ -1,114 +0,0 @@
-/* -*- 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 PAIR_CLASS
-
-PairStyle(lubricateU,PairLubricateU)
-
-#else
-
-#ifndef LMP_PAIR_LUBRICATEU_H
-#define LMP_PAIR_LUBRICATEU_H
-
-#include "pair.h"
-
-namespace LAMMPS_NS {
-
-class PairLubricateU : public Pair {
- public:
- PairLubricateU(class LAMMPS *);
- virtual ~PairLubricateU();
- virtual void compute(int, int);
- virtual void settings(int, char **);
- void coeff(int, char **);
- double init_one(int, int);
- virtual void init_style();
- void write_restart(FILE *);
- void read_restart(FILE *);
- void write_restart_settings(FILE *);
- void read_restart_settings(FILE *);
- int pack_forward_comm(int, int *, double *, int, int *);
- void unpack_forward_comm(int, int, double *);
-
- protected:
- double cut_inner_global,cut_global;
- double mu;
- double rad;
- int flaglog;
- int flagdeform, flagwall;
- int flagVF, flagHI;
- double vol_P;
- class FixWall *wallfix;
-
- double gdot,Ef[3][3];
- double **cut_inner,**cut;
- void allocate();
- double R0,RT0,RS0;
-
- int nmax;
- double **fl,**Tl,**xl;
-
- int cgmax;
- double *bcg,*xcg,*rcg,*rcg1,*pcg,*RU;
-
- void compute_RE();
- virtual void compute_RE(double **);
- void compute_RU();
- virtual void compute_RU(double **);
- virtual void compute_Fh(double **);
- void stage_one();
- void intermediates(int, double **);
- void stage_two(double **);
- void copy_vec_uo(int, double *, double **, double **);
- void copy_uo_vec(int, double **, double **, double *);
- double dot_vec_vec(int , double *, double *);
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Illegal ... command
-
-Self-explanatory. Check the input script syntax and compare to the
-documentation for the command. You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-W: Cannot include log terms without 1/r terms; setting flagHI to 1.
-
-Self-explanatory.
-
-E: Incorrect args for pair coefficients
-
-Self-explanatory. Check the input script or data file.
-
-E: Pair lubricateU requires atom style sphere
-
-Self-explanatory.
-
-E: Pair lubricateU requires ghost atoms store velocity
-
-Use the comm_modify vel yes command to enable this.
-
-E: Pair lubricateU requires monodisperse particles
-
-All particles must be the same finite size.
-
-E: Cannot use multiple fix wall commands with pair lubricateU
-
-Self-explanatory.
-
-*/
diff --git a/src/FLD/pair_lubricateU_poly.cpp b/src/FLD/pair_lubricateU_poly.cpp
deleted file mode 100644
index 6a032987d97680ac4e07c0a9d54d7b51e09bcbfd..0000000000000000000000000000000000000000
--- a/src/FLD/pair_lubricateU_poly.cpp
+++ /dev/null
@@ -1,1234 +0,0 @@
-/* ----------------------------------------------------------------------
- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
- http://lammps.sandia.gov, Sandia National Laboratories
- Steve Plimpton, sjplimp@sandia.gov
-
- Copyright (2003) Sandia Corporation. Under the terms of Contract
- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
- certain rights in this software. This software is distributed under
- the GNU General Public License.
-
- See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing authors: Amit Kumar and Michael Bybee (UIUC)
- Pieter in 't Veld (BASF), code restructuring
- Dave Heine (Corning), polydispersity
-------------------------------------------------------------------------- */
-
-#include <mpi.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include "pair_lubricateU_poly.h"
-#include "atom.h"
-#include "atom_vec.h"
-#include "comm.h"
-#include "force.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
-#include "domain.h"
-#include "update.h"
-#include "modify.h"
-#include "fix.h"
-#include "fix_deform.h"
-#include "fix_wall.h"
-#include "input.h"
-#include "variable.h"
-#include "math_const.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-using namespace MathConst;
-
-#define TOL 1E-3 // tolerance for conjugate gradient
-
-// same as fix_wall.cpp
-
-enum{EDGE,CONSTANT,VARIABLE};
-
-
-/* ---------------------------------------------------------------------- */
-
-PairLubricateUPoly::PairLubricateUPoly(LAMMPS *lmp) :
- PairLubricateU(lmp) {}
-
-/* ----------------------------------------------------------------------
- It first has to solve for the velocity of the particles such that
- the net force on the particles is zero. NOTE: it has to be the last
- type of pair interaction specified in the input file. Also, it
- assumes that no other types of interactions, like k-space, is
- present. As already mentioned, the net force on the particles after
- this pair interaction would be identically zero.
- ---------------------------------------------------------------------- */
-
-void PairLubricateUPoly::compute(int eflag, int vflag)
-{
- int i,j;
-
- int nlocal = atom->nlocal;
- int nghost = atom->nghost;
- int nall = nlocal + nghost;
-
- double **x = atom->x;
- double **f = atom->f;
- double **torque = atom->torque;
-
- if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = vflag_fdotr = 0;
-
- // grow per-atom arrays if necessary
- // need to be atom->nmax in length
-
- if (atom->nmax > nmax) {
- memory->destroy(fl);
- memory->destroy(Tl);
- memory->destroy(xl);
- nmax = atom->nmax;
- memory->create(fl,nmax,3,"pair:fl");
- memory->create(Tl,nmax,3,"pair:Tl");
- memory->create(xl,nmax,3,"pair:xl");
- }
-
- if (6*list->inum > cgmax) {
- memory->sfree(bcg);
- memory->sfree(xcg);
- memory->sfree(rcg);
- memory->sfree(rcg1);
- memory->sfree(pcg);
- memory->sfree(RU);
- cgmax = 6*list->inum;
- memory->create(bcg,cgmax,"pair:bcg");
- memory->create(xcg,cgmax,"pair:bcg");
- memory->create(rcg,cgmax,"pair:bcg");
- memory->create(rcg1,cgmax,"pair:bcg");
- memory->create(pcg,cgmax,"pair:bcg");
- memory->create(RU,cgmax,"pair:bcg");
- }
-
- // Added to implement midpoint integration scheme
- // Save force, torque found so far. Also save the positions
-
- for (i=0;i<nlocal+nghost;i++) {
- for (j=0;j<3;j++) {
- fl[i][j] = f[i][j];
- Tl[i][j] = torque[i][j];
- xl[i][j] = x[i][j];
- }
- }
-
- // Stage one of Midpoint method
- // Solve for velocities based on intial positions
-
- iterate(atom->x,1);
-
- // Find positions at half the timestep and store in xl
-
- intermediates(nall,xl);
-
- // Store back the saved forces and torques in original arrays
-
- for(i=0;i<nlocal+nghost;i++) {
- for(j=0;j<3;j++) {
- f[i][j] = fl[i][j];
- torque[i][j] = Tl[i][j];
- }
- }
-
- // Stage two: This will give the final velocities
-
- iterate(xl,2);
-}
-
-/* ------------------------------------------------------------------------
- Stage one of midpoint method
-------------------------------------------------------------------------- */
-
-void PairLubricateUPoly::iterate(double **x, int stage)
-{
- int i,j,ii;
-
- int inum = list->inum;
- int *ilist = list->ilist;
- int newton_pair = force->newton_pair;
-
- double alpha,beta;
- double normi,error,normig;
- double send[2],recv[2],rcg_dot_rcg;
-
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
-
- // First compute R_FE*E
-
- compute_RE(x);
-
- // Reverse communication of forces and torques to
- // accumulate the net force on each of the particles
-
- if (newton_pair) comm->reverse_comm();
-
- // CONJUGATE GRADIENT
- // Find the right hand side= -ve of all forces/torques
- // b = 6*Npart in overall size
-
- for(ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- for (j = 0; j < 3; j++) {
- bcg[6*ii+j] = -f[i][j];
- bcg[6*ii+j+3] = -torque[i][j];
- }
- }
-
- // Start solving the equation : F^H = -F^P -F^B - F^H_{Ef}
- // Store initial guess for velocity and angular-velocities/angular momentum
- // NOTE velocities and angular velocities are assumed relative to the fluid
-
- for (ii=0;ii<inum;ii++)
- for (j=0;j<3;j++) {
- xcg[6*ii+j] = 0.0;
- xcg[6*ii+j+3] = 0.0;
- }
-
- // Copy initial guess to the global arrays to be acted upon by R_{FU}
- // and returned by f and torque arrays
-
- copy_vec_uo(inum,xcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- // Find initial residual
-
- compute_RU(x);
-
- // reverse communication of forces and torques
-
- if (newton_pair) comm->reverse_comm();
-
- copy_uo_vec(inum,f,torque,RU);
-
- for (i=0;i<6*inum;i++)
- rcg[i] = bcg[i] - RU[i];
-
- // Set initial conjugate direction
-
- for (i=0;i<6*inum;i++)
- pcg[i] = rcg[i];
-
- // Find initial norm of the residual or norm of the RHS (either is fine)
-
- normi = dot_vec_vec(6*inum,bcg,bcg);
-
- MPI_Allreduce(&normi,&normig,1,MPI_DOUBLE,MPI_SUM,world);
-
- // Loop until convergence
-
- do {
- // find R*p
-
- copy_vec_uo(inum,pcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- compute_RU(x);
-
- // reverse communication of forces and torques
-
- if (newton_pair) comm->reverse_comm();
-
- copy_uo_vec(inum,f,torque,RU);
-
- // Find alpha
-
- send[0] = dot_vec_vec(6*inum,rcg,rcg);
- send[1] = dot_vec_vec(6*inum,RU,pcg);
-
- MPI_Allreduce(send,recv,2,MPI_DOUBLE,MPI_SUM,world);
-
- alpha = recv[0]/recv[1];
- rcg_dot_rcg = recv[0];
-
- // Find new x
-
- for (i=0;i<6*inum;i++)
- xcg[i] = xcg[i] + alpha*pcg[i];
-
- // find new residual
-
- for (i=0;i<6*inum;i++)
- rcg1[i] = rcg[i] - alpha*RU[i];
-
- // find beta
-
- send[0] = dot_vec_vec(6*inum,rcg1,rcg1);
-
- MPI_Allreduce(send,recv,1,MPI_DOUBLE,MPI_SUM,world);
-
- beta = recv[0]/rcg_dot_rcg;
-
- // Find new conjugate direction
-
- for (i=0;i<6*inum;i++)
- pcg[i] = rcg1[i] + beta*pcg[i];
-
- for (i=0;i<6*inum;i++)
- rcg[i] = rcg1[i];
-
- // Find relative error
-
- error = sqrt(recv[0]/normig);
-
- } while (error > TOL);
-
-
- // update the final converged velocities in respective arrays
-
- copy_vec_uo(inum,xcg,v,omega);
-
- // set velocities for ghost particles
-
- comm->forward_comm_pair(this);
-
- // compute the viscosity/pressure
-
- if (evflag && stage == 2) compute_Fh(x);
-
- // find actual particle's velocities from relative velocities
- // only non-zero component of fluid's vel : vx=gdot*y and wz=-gdot/2
-
- for (ii=0;ii<inum;ii++) {
- i = ilist[ii];
- v[i][0] = v[i][0] + gdot*x[i][1];
- omega[i][2] = omega[i][2] - gdot/2.0;
- }
-}
-
-/* ------------------------------------------------------------------------
- This function computes the final hydrodynamic force once the
- velocities have converged.
- ------------------------------------------------------------------------- */
-
-void PairLubricateUPoly::compute_Fh(double **x)
-{
- int i,j,ii,jj,inum,jnum,itype,jtype;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int nghost = atom->nghost;
- int newton_pair = force->newton_pair;
-
- double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
- double rsq,r,h_sep,radi,radj;
- double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
- double vt1,vt2,vt3;
- double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3],pre[2];
-
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
- double *radius = atom->radius;
-
- double beta[2][5];
- double vxmu2f = force->vxmu2f;
- double a_sq = 0.0;
- double a_sh = 0.0;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- beta[0][0] = beta[1][0] = beta[1][4] = 0.0;
-
- // This section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (int j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- //R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f);
- //RT0 = 8*MY_PI*mu;
- RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- //R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- //RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
-
- // end of R0 adjustment code
- // Set force to zero which is the final value after this pair interaction
-
- for (i=0;i<nlocal+nghost;i++)
- for (j=0;j<3;j++) {
- f[i][j] = 0.0;
- torque[i][j] = 0.0;
- }
-
- // reverse communication of forces and torques
-
- if (newton_pair) comm->reverse_comm(); // not really needed
-
- // Find additional contribution from the stresslets
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
- pre[1] = 8.0*(pre[0] = MY_PI*mu*radi)*radi*radi; // BROKEN?? Should be "+"??
- pre[0] *= 6.0;
-
- // Find the contribution to stress from isotropic RS0
- // Set psuedo force to obtain the required contribution
- // need to set delx and fy only
-
- fx = 0.0; delx = radi;
- fy = vxmu2f*RS0*pow(radi,3.0)*gdot/2.0/radi; dely = 0.0;
- fz = 0.0; delz = 0.0;
- if (evflag)
- ev_tally_xyz(i,i,nlocal,newton_pair,0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
-
- // Find angular velocity
-
- wi[0] = omega[i][0];
- wi[1] = omega[i][1];
- wi[2] = omega[i][2];
-
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- 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);
-
- // Use omega directly if it exists, else angmom
- // 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];
-
- // loc of the point of closest approach on particle i from its cente
- // POC for j is in opposite direction as for i
-
- 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;
-
- h_sep = r - radi-radj;
-
- // velocity at the point of closest approach on both particles
- // v = v + omega_cross_xl
-
- // particle i
-
- vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]);
- vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]);
- vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]);
-
- // particle j
-
- vj[0] = v[j][0] + (wj[1]*jl[2] - wj[2]*jl[1]);
- vj[1] = v[j][1] + (wj[2]*jl[0] - wj[0]*jl[2]);
- vj[2] = v[j][2] + (wj[0]*jl[1] - wj[1]*jl[0]);
-
-
- // Relative velocity at the point of closest approach
- // include contribution from Einf of the fluid
-
- vr1 = vi[0] - vj[0] -
- 2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]);
- vr2 = vi[1] - vj[1] -
- 2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]);
- vr3 = vi[2] - vj[2] -
- 2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[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;
-
- // Find the scalar resistances a_sq, a_sh and a_pu
-
- // 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;
- beta[0][1] = radj/radi;
- beta[1][1] = 1.0 + beta[0][1];
-
- /*beta0 = radj/radi;
- beta1 = 1.0 + beta0;*/
-
- // Scalar resistances
-
- if (flaglog) {
- beta[0][2] = beta[0][1]*beta[0][1];
- beta[0][3] = beta[0][2]*beta[0][1];
- beta[0][4] = beta[0][3]*beta[0][1];
- beta[1][2] = beta[1][1]*beta[1][1];
- beta[1][3] = beta[1][2]*beta[1][1];
- double log_h_sep_beta13 = log(1.0/h_sep)/beta[1][3];
- double h_sep_beta11 = h_sep/beta[1][1];
-
- a_sq = pre[0]*(beta[0][2]/beta[1][2]/h_sep
- +((0.2+1.4*beta[0][1]+0.2*beta[0][2])
- +(1.0+18.0*(beta[0][1]+beta[0][3])-29.0*beta[0][2]
- +beta[0][4])*h_sep_beta11/21.0)*log_h_sep_beta13);
-
- a_sh = pre[0]*((8.0*(beta[0][1]+beta[0][3])+4.0*beta[0][2])/15.0
- +(64.0-180.0*(beta[0][1]+beta[0][3])+232.0*beta[0][2]
- +64.0*beta[0][4])*h_sep_beta11/375.0)*log_h_sep_beta13;
-
- /*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;*/
- } else {
- //a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep);
- a_sq = pre[0]*(beta[0][1]*beta[0][1]/(beta[1][1]*beta[1][1]*h_sep));
- }
-
- // Find force due to squeeze type motion
-
- fx = a_sq*vn1;
- fy = a_sq*vn2;
- fz = a_sq*vn3;
-
- // Find 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 to obtain in appropriate units
-
- fx = vxmu2f*fx;
- fy = vxmu2f*fy;
- fz = vxmu2f*fz;
-
- // set j = nlocal so that only I gets tallied
-
- if (evflag) ev_tally_xyz(i,nlocal,nlocal,0,
- 0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
- }
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- computes R_FU * U
----------------------------------------------------------------------- */
-
-void PairLubricateUPoly::compute_RU(double **x)
-{
- int i,j,ii,jj,inum,jnum,itype,jtype;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int nghost = atom->nghost;
-
- double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
- double rsq,r,radi,radj,h_sep;
- double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
- double vt1,vt2,vt3,wdotn,wt1,wt2,wt3;
- double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3],pre[2];
-
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
- double *radius = atom->radius;
-
- double beta[2][5];
- double vxmu2f = force->vxmu2f;
- double a_sq = 0.0;
- double a_sh = 0.0;
- double a_pu = 0.0;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- beta[0][0] = beta[1][0] = beta[1][4] = 0.0;
-
- // This section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu;
- // RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- // RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
-
- // end of R0 adjustment code
-
- // Initialize f to zero
-
- for (i=0;i<nlocal+nghost;i++)
- for (j=0;j<3;j++) {
- f[i][j] = 0.0;
- torque[i][j] = 0.0;
- }
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
- pre[1] = 8.0*(pre[0] = MY_PI*mu*radi)*radi*radi;
- pre[0] *= 6.0;
-
- // Find angular velocity
-
- wi[0] = omega[i][0];
- wi[1] = omega[i][1];
- wi[2] = omega[i][2];
-
- // Contribution due to the isotropic terms
-
- 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 radi3 = radi*radi*radi;
- torque[i][0] += -vxmu2f*RT0*radi3*wi[0];
- torque[i][1] += -vxmu2f*RT0*radi3*wi[1];
- torque[i][2] += -vxmu2f*RT0*radi3*wi[2];
-
- if (!flagHI) continue;
-
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- 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);
-
- wj[0] = omega[j][0];
- wj[1] = omega[j][1];
- wj[2] = omega[j][2];
-
- // loc of the 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
-
- // particle i
-
- vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]);
- vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]);
- vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]);
-
- // particle j
-
- vj[0] = v[j][0] + (wj[1]*jl[2] - wj[2]*jl[1]);
- vj[1] = v[j][1] + (wj[2]*jl[0] - wj[0]*jl[2]);
- vj[2] = v[j][2] + (wj[0]*jl[1] - wj[1]*jl[0]);
-
- // Find the scalar resistances a_sq and a_sh
-
- h_sep = r - radi-radj;
-
- // 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;
- beta[0][1] = radj/radi;
- beta[1][1] = 1.0 + beta[0][1];
-
- // Scalar resistances
-
- if (flaglog) {
- beta[0][2] = beta[0][1]*beta[0][1];
- beta[0][3] = beta[0][2]*beta[0][1];
- beta[0][4] = beta[0][3]*beta[0][1];
- beta[1][2] = beta[1][1]*beta[1][1];
- beta[1][3] = beta[1][2]*beta[1][1];
- double log_h_sep_beta13 = log(1.0/h_sep)/beta[1][3];
- double h_sep_beta11 = h_sep/beta[1][1];
-
- a_sq = pre[0]*(beta[0][2]/beta[1][2]/h_sep
- +((0.2+1.4*beta[0][1]+0.2*beta[0][2])
- +(1.0+18.0*(beta[0][1]+beta[0][3])-29.0*beta[0][2]
- +beta[0][4])*h_sep_beta11/21.0)*log_h_sep_beta13);
-
- a_sh = pre[0]*((8.0*(beta[0][1]+beta[0][3])+4.0*beta[0][2])/15.0
- +(64.0-180.0*(beta[0][1]+beta[0][3])+232.0*beta[0][2]
- +64.0*beta[0][4])*h_sep_beta11/375.0)*log_h_sep_beta13;
-
- a_pu = pre[1]*((0.4*beta[0][1]+0.1*beta[0][2])*beta[1][1]
- +(0.128-0.132*beta[0][1]+0.332*beta[0][2]
- +0.172*beta[0][3])*h_sep)*log_h_sep_beta13;
-
- /*//a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1/h_sep));
- 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 = pre[0]*(beta[0][1]*beta[0][1]/(beta[1][1]*beta[1][1]*h_sep));
-
- // Relative velocity at the point of closest approach
-
- 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;
-
- // Find force due to squeeze type motion
-
- fx = a_sq*vn1;
- fy = a_sq*vn2;
- fz = a_sq*vn3;
-
- // Find 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 to obtain in appropriate units
-
- fx = vxmu2f*fx;
- fy = vxmu2f*fy;
- fz = vxmu2f*fz;
-
- // Add to the total forc
-
- f[i][0] -= fx;
- f[i][1] -= fy;
- f[i][2] -= fz;
-
- // Find 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;
-
- // Why a scale factor ?
-
- 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;
-
- // add to total
-
- torque[i][0] -= vxmu2f*tx;
- torque[i][1] -= vxmu2f*ty;
- torque[i][2] -= vxmu2f*tz;
- }
- }
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- This computes R_{FE}*E , where E is the rate of strain of tensor which is
- known apriori, as it depends only on the known fluid velocity.
- So, this part of the hydrodynamic interaction can be pre computed and
- transferred to the RHS
- ---------------------------------------------------------------------- */
-
-void PairLubricateUPoly::compute_RE(double **x)
-{
- int i,j,ii,jj,inum,jnum,itype,jtype;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
- int *type = atom->type;
-
- double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
- double rsq,r,h_sep,radi,radj;
- //double beta0,beta1,lhsep;
- double vr1,vr2,vr3,vnnr,vn1,vn2,vn3;
- double vt1,vt2,vt3;
- double xl[3],pre[2];
-
- double **f = atom->f;
- double **torque = atom->torque;
- double *radius = atom->radius;
-
- double beta[2][5];
- double vxmu2f = force->vxmu2f;
- double a_sq = 0.0;
- double a_sh = 0.0;
-
- if (!flagHI) return;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- beta[0][0] = beta[1][0] = beta[1][4] = 0.0;
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- radi = radius[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
- pre[1] = 8.0*(pre[0] = MY_PI*mu*radi)*radi*radi;
- pre[0] *= 6.0;
-
- // No contribution from isotropic terms due to E
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- 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);
-
- // loc of the 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;
-
- // Find the scalar resistances a_sq and a_sh
-
- h_sep = r - radi-radj;
-
- // 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;
- beta[0][1] = radj/radi;
- beta[1][1] = 1.0 + beta[0][1];
-
- /*beta0 = radj/radi;
- beta1 = 1.0 + beta0;
- lhsep = log(1.0/h_sep);*/
-
- // Scalar resistance for Squeeze type motions
-
-
- if (flaglog) {
- beta[0][2] = beta[0][1]*beta[0][1];
- beta[0][3] = beta[0][2]*beta[0][1];
- beta[0][4] = beta[0][3]*beta[0][1];
- beta[1][2] = beta[1][1]*beta[1][1];
- beta[1][3] = beta[1][2]*beta[1][1];
- double log_h_sep_beta13 = log(1.0/h_sep)/beta[1][3];
- double h_sep_beta11 = h_sep/beta[1][1];
-
- a_sq = pre[0]*(beta[0][2]/beta[1][2]/h_sep
- +((0.2+1.4*beta[0][1]+0.2*beta[0][2])
- +(1.0+18.0*(beta[0][1]+beta[0][3])-29.0*beta[0][2]
- +beta[0][4])*h_sep_beta11/21.0)*log_h_sep_beta13);
-
- a_sh = pre[0]*((8.0*(beta[0][1]+beta[0][3])+4.0*beta[0][2])/15.0
- +(64.0-180.0*(beta[0][1]+beta[0][3])+232.0*beta[0][2]
- +64.0*beta[0][4])*h_sep_beta11/375.0)*log_h_sep_beta13;
- } else
- a_sq = pre[0]*(beta[0][1]*beta[0][1]/(beta[1][1]*beta[1][1]*h_sep));
-
- // Relative velocity at the point of closest approach due to Ef only
-
- vr1 = -2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]);
- vr2 = -2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]);
- vr3 = -2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[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;
-
- // Find force due to squeeze type motion
-
- fx = a_sq*vn1;
- fy = a_sq*vn2;
- fz = a_sq*vn3;
-
- // Find 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 to obtain in appropriate units
-
- fx = vxmu2f*fx;
- fy = vxmu2f*fy;
- fz = vxmu2f*fz;
-
- // Add to the total forc
-
- f[i][0] -= fx;
- f[i][1] -= fy;
- f[i][2] -= fz;
-
- // Find 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;
-
- // Why a scale factor ?
-
- torque[i][0] -= vxmu2f*tx;
- torque[i][1] -= vxmu2f*ty;
- torque[i][2] -= vxmu2f*tz;
-
- }
- }
- }
- }
-}
-
-/*-----------------------------------------------------------------------
- global settings
-------------------------------------------------------------------------- */
-
-void PairLubricateUPoly::settings(int narg, char **arg)
-{
- if (narg < 5 || narg > 7) error->all(FLERR,"Illegal pair_style command");
-
- mu = force->numeric(FLERR,arg[0]);
- flaglog = force->inumeric(FLERR,arg[1]);
- cut_inner_global = force->numeric(FLERR,arg[2]);
- cut_global = force->numeric(FLERR,arg[3]);
- gdot = force->numeric(FLERR,arg[4]);
-
- flagHI = flagVF = 1;
- if (narg >= 6) flagHI = force->inumeric(FLERR,arg[5]);
- if (narg == 7) flagVF = force->inumeric(FLERR,arg[6]);
-
- if (flaglog == 1 && flagHI == 0) {
- error->warning(FLERR,"Cannot include log terms without 1/r terms; "
- "setting flagHI to 1");
- flagHI = 1;
- }
-
- // reset cutoffs that have been explicitly set
-
- if (allocated) {
- int i,j;
- for (i = 1; i <= atom->ntypes; i++)
- for (j = i+1; j <= atom->ntypes; j++)
- if (setflag[i][j]) {
- cut_inner[i][j] = cut_inner_global;
- cut[i][j] = cut_global;
- }
- }
-
- // Store the rate of strain tensor
-
- Ef[0][0] = 0.0;
- Ef[0][1] = 0.5*gdot;
- Ef[0][2] = 0.0;
- Ef[1][0] = 0.5*gdot;
- Ef[1][1] = 0.0;
- Ef[1][2] = 0.0;
- Ef[2][0] = 0.0;
- Ef[2][1] = 0.0;
- Ef[2][2] = 0.0;
-}
-
-/* ----------------------------------------------------------------------
- init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairLubricateUPoly::init_style()
-{
- if (force->newton_pair == 1)
- error->all(FLERR,"Pair lubricateU/poly requires newton pair off");
- if (comm->ghost_velocity == 0)
- error->all(FLERR,
- "Pair lubricateU/poly requires ghost atoms store velocity");
- if (!atom->sphere_flag)
- error->all(FLERR,"Pair lubricate/poly requires atom style sphere");
-
- // insure all particles are finite-size
- // for pair hybrid, should limit test to types using the pair style
-
- double *radius = atom->radius;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (radius[i] == 0.0)
- error->one(FLERR,"Pair lubricate/poly requires extended particles");
-
- // Set the isotropic constants depending on the volume fraction
-
- // Find the total volume
- // check for fix deform, if exists it must use "remap v"
- // If box will change volume, set appropriate flag so that volume
- // and v.f. corrections are re-calculated at every step.
- //
- // If available volume is different from box volume
- // due to walls, set volume appropriately; if walls will
- // move, set appropriate flag so that volume and v.f. corrections
- // are re-calculated at every step.
-
- flagdeform = flagwall = 0;
- for (int i = 0; i < modify->nfix; i++){
- if (strcmp(modify->fix[i]->style,"deform") == 0)
- flagdeform = 1;
- else if (strstr(modify->fix[i]->style,"wall") != NULL){
- if (flagwall)
- error->all(FLERR,
- "Cannot use multiple fix wall commands with "
- "pair lubricateU");
- flagwall = 1; // Walls exist
- wallfix = (FixWall *) modify->fix[i];
- if (wallfix->xflag) flagwall = 2; // Moving walls exist
- }
- }
-
- // set the isotropic constants depending on the volume fraction
- // vol_T = total volumeshearing = flagdeform = flagwall = 0;
-
- double vol_T, wallcoord;
- if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
- else {
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallfix->xindex[m] = input->variable->find(wallfix->xstr[m]);
- //Since fix->wall->init happens after pair->init_style
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
-
- else wallcoord = wallfix->coord0[m];
-
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
- (wallhi[2] - walllo[2]);
- }
-
- // Assuming monodisperse spheres, find the volume of the particles
-
- double volP = 0.0;
- for (int i = 0; i < nlocal; i++)
- volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3.0);
- MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world);
-
- double vol_f = vol_P/vol_T;
-
- //DRH volume fraction needs to be defined manually
- // if excluded volume regions are present
- // vol_f=0.5;
-
- if (!flagVF) vol_f = 0;
-
- if (!comm->me) {
- if(logfile)
- fprintf(logfile, "lubricateU: vol_f = %g, vol_p = %g, vol_T = %g\n",
- vol_f,vol_P,vol_T);
- if (screen)
- fprintf(screen, "lubricateU: vol_f = %g, vol_p = %g, vol_T = %g\n",
- vol_f,vol_P,vol_T);
- }
-
- // Set the isotropic constant
-
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu; // Not needed actually
- RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
-
- int irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->full = 1;
-}
diff --git a/src/FLD/pair_lubricateU_poly.h b/src/FLD/pair_lubricateU_poly.h
deleted file mode 100644
index 0f22809cf212c56677d39dc3e6efbfcf23f39157..0000000000000000000000000000000000000000
--- a/src/FLD/pair_lubricateU_poly.h
+++ /dev/null
@@ -1,83 +0,0 @@
-/* -*- 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 PAIR_CLASS
-
-PairStyle(lubricateU/poly,PairLubricateUPoly)
-
-#else
-
-#ifndef LMP_PAIR_LUBRICATEU_POLY_H
-#define LMP_PAIR_LUBRICATEU_POLY_H
-
-#include "pair_lubricateU.h"
-
-namespace LAMMPS_NS {
-
-class PairLubricateUPoly : public PairLubricateU {
- public:
- PairLubricateUPoly(class LAMMPS *);
- ~PairLubricateUPoly() {}
- void compute(int, int);
- void settings(int, char **);
- void init_style();
-
- private:
- double vol_P;
- int flagdeform, flagwall, flagVF, flagHI;
- class FixWall *wallfix;
-
- void iterate(double **, int);
- void compute_RE(double **);
- void compute_RU(double **);
- void compute_Fh(double **);
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Illegal ... command
-
-Self-explanatory. Check the input script syntax and compare to the
-documentation for the command. You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-W: Cannot include log terms without 1/r terms; setting flagHI to 1
-
-Self-explanatory.
-
-E: Pair lubricateU/poly requires newton pair off
-
-Self-explanatory.
-
-E: Pair lubricateU/poly requires ghost atoms store velocity
-
-Use the comm_modify vel yes command to enable this.
-
-E: Pair lubricate/poly requires atom style sphere
-
-Self-explanatory.
-
-E: Pair lubricate/poly requires extended particles
-
-One of the particles has radius 0.0.
-
-E: Cannot use multiple fix wall commands with pair lubricateU
-
-Self-explanatory.
-
-*/
diff --git a/src/FLD/pair_lubricate_poly.cpp b/src/FLD/pair_lubricate_poly.cpp
deleted file mode 100644
index 3701e2d78984e8b1e67fdab3b51c6e503f010b82..0000000000000000000000000000000000000000
--- a/src/FLD/pair_lubricate_poly.cpp
+++ /dev/null
@@ -1,555 +0,0 @@
-/* ----------------------------------------------------------------------
- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
- http://lammps.sandia.gov, Sandia National Laboratories
- Steve Plimpton, sjplimp@sandia.gov
-
- Copyright (2003) Sandia Corporation. Under the terms of Contract
- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
- certain rights in this software. This software is distributed under
- the GNU General Public License.
-
- See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing authors: Randy Schunk (SNL)
- Amit Kumar and Michael Bybee (UIUC)
- Dave Heine (Corning), polydispersity
-------------------------------------------------------------------------- */
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include "pair_lubricate_poly.h"
-#include "atom.h"
-#include "atom_vec.h"
-#include "comm.h"
-#include "force.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
-#include "domain.h"
-#include "modify.h"
-#include "fix.h"
-#include "fix_deform.h"
-#include "memory.h"
-#include "random_mars.h"
-#include "fix_wall.h"
-#include "input.h"
-#include "variable.h"
-#include "math_const.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-using namespace MathConst;
-
-// same as fix_deform.cpp
-
-enum{NO_REMAP,X_REMAP,V_REMAP};
-
-
-// same as fix_wall.cpp
-
-enum{EDGE,CONSTANT,VARIABLE};
-
-/* ---------------------------------------------------------------------- */
-
-PairLubricatePoly::PairLubricatePoly(LAMMPS *lmp) : PairLubricate(lmp)
-{
- no_virial_fdotr_compute = 1;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairLubricatePoly::compute(int eflag, int vflag)
-{
- int i,j,ii,jj,inum,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;
-
- if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = vflag_fdotr = 0;
-
- double **x = atom->x;
- double **v = atom->v;
- double **f = atom->f;
- double **omega = atom->omega;
- double **torque = atom->torque;
- double *radius = atom->radius;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int newton_pair = force->newton_pair;
-
- inum = list->inum;
- 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 = 0; ii < inum; 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
-
- comm->forward_comm_pair(this);
- }
-
- // This section of code adjusts R0/RT0/RS0 if necessary due to changes
- // in the volume fraction as a result of fix deform or moving walls
-
- double dims[3], wallcoord;
- if (flagVF) // Flag for volume fraction corrections
- if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
- if (flagdeform && !flagwall)
- for (j = 0; j < 3; j++)
- dims[j] = domain->prd[j];
- else if (flagwall == 2 || (flagdeform && flagwall == 1)){
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- for (int j = 0; j < 3; j++)
- dims[j] = wallhi[j] - walllo[j];
- }
- double vol_T = dims[0]*dims[1]*dims[2];
- double vol_f = vol_P/vol_T;
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu;
- RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
- }
-
-
- // end of R0 adjustment code
-
- 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];
- radi = radius[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 radi3 = radi*radi*radi;
- torque[i][0] -= vxmu2f*RT0*radi3*wi[0];
- torque[i][1] -= vxmu2f*RT0*radi3*wi[1];
- torque[i][2] -= vxmu2f*RT0*radi3*wi[2];
-
- if (shearing && vflag_either) {
- vRS0 = -vxmu2f * RS0*radi3;
- 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]);
- }
- }
-
- if (!flagHI) continue;
-
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- 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;
-
- // 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.0)*log(1.0/h_sep);
- a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0 *
- pow(beta0,3.0)+pow(beta0,4.0))/21.0/pow(beta1,4.0) *
- 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.0) *
- log(1.0/h_sep);
- a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3.0) +
- 16.0*pow(beta0,4.0))/375.0/pow(beta1,4.0) *
- 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.0))/250.0/pow(beta1,3.0)*h_sep*log(1.0/h_sep);
- a_pu *= 8.0*MY_PI*mu*pow(radi,3.0);
- } 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;
-
- }
-
- // set j = nlocal so that only I gets tallied
-
- if (evflag) ev_tally_xyz(i,nlocal,nlocal,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 = 0; ii < inum; ii++) {
- i = ilist[ii];
- itype = type[i];
- radi = atom->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];
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairLubricatePoly::init_style()
-{
- if (force->newton_pair == 1)
- error->all(FLERR,"Pair lubricate/poly requires newton pair off");
- if (comm->ghost_velocity == 0)
- error->all(FLERR,
- "Pair lubricate/poly requires ghost atoms store velocity");
- if (!atom->sphere_flag)
- error->all(FLERR,"Pair lubricate/poly requires atom style sphere");
-
- // ensure all particles are finite-size
- // for pair hybrid, should limit test to types using the pair style
-
- double *radius = atom->radius;
- int nlocal = atom->nlocal;
-
- for (int i = 0; i < nlocal; i++)
- if (radius[i] == 0.0)
- error->one(FLERR,"Pair lubricate/poly requires extended particles");
-
- int irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->full = 1;
-
- // set the isotropic constants that depend on the volume fraction
- // vol_T = total volume
-
- // check for fix deform, if exists it must use "remap v"
- // If box will change volume, set appropriate flag so that volume
- // and v.f. corrections are re-calculated at every step.
-
- // if available volume is different from box volume
- // due to walls, set volume appropriately; if walls will
- // move, set appropriate flag so that volume and v.f. corrections
- // are re-calculated at every step.
-
- shearing = flagdeform = flagwall = 0;
- for (int i = 0; i < modify->nfix; i++){
- if (strcmp(modify->fix[i]->style,"deform") == 0) {
- shearing = flagdeform = 1;
- if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP)
- error->all(FLERR,"Using pair lubricate with inconsistent "
- "fix deform remap option");
- }
- if (strstr(modify->fix[i]->style,"wall") != NULL) {
- if (flagwall)
- error->all(FLERR,
- "Cannot use multiple fix wall commands with "
- "pair lubricate/poly");
- flagwall = 1; // Walls exist
- wallfix = (FixWall *) modify->fix[i];
- if (wallfix->xflag) flagwall = 2; // Moving walls exist
- }
-
- if (strstr(modify->fix[i]->style,"wall") != NULL){
- flagwall = 1; // Walls exist
- if (((FixWall *) modify->fix[i])->xflag ) {
- flagwall = 2; // Moving walls exist
- wallfix = (FixWall *) modify->fix[i];
- }
- }
- }
-
- double vol_T;
- double wallcoord;
- if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
- else {
- double wallhi[3], walllo[3];
- for (int j = 0; j < 3; j++){
- wallhi[j] = domain->prd[j];
- walllo[j] = 0;
- }
- for (int m = 0; m < wallfix->nwall; m++){
- int dim = wallfix->wallwhich[m] / 2;
- int side = wallfix->wallwhich[m] % 2;
- if (wallfix->xstyle[m] == VARIABLE){
- wallfix->xindex[m] = input->variable->find(wallfix->xstr[m]);
- //Since fix->wall->init happens after pair->init_style
- wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
- }
- else wallcoord = wallfix->coord0[m];
-
- if (side == 0) walllo[dim] = wallcoord;
- else wallhi[dim] = wallcoord;
- }
- vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
- (wallhi[2] - walllo[2]);
- }
-
- double volP = 0.0;
- for (int i = 0; i < nlocal; i++)
- volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3.0);
- MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world);
-
- double vol_f = vol_P/vol_T;
-
- if (!flagVF) vol_f = 0;
-
- // set isotropic constants
-
- if (flaglog == 0) {
- R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f);
- RT0 = 8*MY_PI*mu;
- RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
- } else {
- R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
- RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
- RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
- }
-
- // check for fix deform, if exists it must use "remap v"
-
- shearing = 0;
- for (int i = 0; i < modify->nfix; i++)
- if (strcmp(modify->fix[i]->style,"deform") == 0) {
- shearing = 1;
- if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP)
- error->all(FLERR,"Using pair lubricate/poly with inconsistent "
- "fix deform remap option");
- }
-
- // set Ef = 0 since used whether shearing or not
-
- Ef[0][0] = Ef[0][1] = Ef[0][2] = 0.0;
- Ef[1][0] = Ef[1][1] = Ef[1][2] = 0.0;
- Ef[2][0] = Ef[2][1] = Ef[2][2] = 0.0;
-}
diff --git a/src/FLD/pair_lubricate_poly.h b/src/FLD/pair_lubricate_poly.h
deleted file mode 100644
index 5e58b24254358d5685302ca68047de1b083f67c2..0000000000000000000000000000000000000000
--- a/src/FLD/pair_lubricate_poly.h
+++ /dev/null
@@ -1,70 +0,0 @@
-/* -*- 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 PAIR_CLASS
-
-PairStyle(lubricate/poly,PairLubricatePoly)
-
-#else
-
-#ifndef LMP_PAIR_LUBRICATE_POLY_H
-#define LMP_PAIR_LUBRICATE_POLY_H
-
-#include "pair_lubricate.h"
-
-namespace LAMMPS_NS {
-
-class PairLubricatePoly : public PairLubricate {
- public:
- PairLubricatePoly(class LAMMPS *);
- ~PairLubricatePoly() {}
- void compute(int, int);
- void init_style();
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Pair lubricate/poly requires newton pair off
-
-Self-explanatory.
-
-E: Pair lubricate/poly requires ghost atoms store velocity
-
-Use the comm_modify vel yes command to enable this.
-
-E: Pair lubricate/poly requires atom style sphere
-
-Self-explanatory.
-
-E: Pair lubricate/poly requires extended particles
-
-One of the particles has radius 0.0.
-
-E: Using pair lubricate with inconsistent fix deform remap option
-
-Must use remap v option with fix deform with this pair style.
-
-E: Cannot use multiple fix wall commands with pair lubricate/poly
-
-Self-explanatory.
-
-E: Using pair lubricate/poly with inconsistent fix deform remap option
-
-If fix deform is used, the remap v option is required.
-
-*/
diff --git a/src/Makefile b/src/Makefile
index aaa93bdbc8b79b5f69b6d79bf5eada3e76540139..34838d7fd5c4bde7697fa2d74d019b2f967f1aae 100755
--- a/src/Makefile
+++ b/src/Makefile
@@ -41,7 +41,7 @@ endif
# Package variables
-PACKAGE = asphere body class2 colloid compress coreshell dipole fld gpu \
+PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \
granular kim \
kokkos kspace manybody mc meam misc molecule mpiio opt peri poems \
python qeq reax replica rigid shock snap srd voronoi xtc