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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3342 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2009-11-06 21:02:28 +00:00
parent 21b5ea51a7
commit cff9030e66
10 changed files with 1244 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/COLLOID/Install.csh
View File

@ -4,21 +4,33 @@ if ($1 == 1) then
cp style_colloid.h ..
cp atom_vec_colloid.cpp ..
cp fix_wall_colloid.cpp ..
cp pair_colloid.cpp ..
cp pair_lubricate.cpp ..
cp pair_yukawa_colloid.cpp ..
cp atom_vec_colloid.h ..
cp fix_wall_colloid.h ..
cp pair_colloid.h ..
cp pair_lubricate.h ..
cp pair_yukawa_colloid.h ..
else if ($1 == 0) then
rm ../style_colloid.h
touch ../style_colloid.h
rm ../atom_vec_colloid.cpp
rm ../fix_wall_colloid.cpp
rm ../pair_colloid.cpp
rm ../pair_lubricate.cpp
rm ../pair_yukawa_colloid.cpp
rm ../atom_vec_colloid.h
rm ../fix_wall_colloid.h
rm ../pair_colloid.h
rm ../pair_lubricate.h
rm ../pair_yukawa_colloid.h
endif

644
src/COLLOID/atom_vec_colloid.cpp Normal file
View File

@ -0,0 +1,644 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "atom_vec_colloid.h"
#include "atom.h"
#include "force.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define DELTA 10000
/* ---------------------------------------------------------------------- */
AtomVecColloid::AtomVecColloid(LAMMPS *lmp, int narg, char **arg) :
AtomVec(lmp, narg, arg)
{
mass_type = 1;
shape_type = 1;
comm_x_only = comm_f_only = 0;
ghost_velocity = 1;
size_comm = 9;
size_reverse = 6;
size_border = 12;
size_data_atom = 5;
size_data_vel = 7;
xcol_data = 3;
atom->omega_flag = atom->torque_flag = 1;
}
/* ----------------------------------------------------------------------
grow atom arrays
n = 0 grows arrays by DELTA
n > 0 allocates arrays to size n
------------------------------------------------------------------------- */
void AtomVecColloid::grow(int n)
{
if (n == 0) nmax += DELTA;
else nmax = n;
atom->nmax = nmax;
tag = atom->tag = (int *)
memory->srealloc(atom->tag,nmax*sizeof(int),"atom:tag");
type = atom->type = (int *)
memory->srealloc(atom->type,nmax*sizeof(int),"atom:type");
mask = atom->mask = (int *)
memory->srealloc(atom->mask,nmax*sizeof(int),"atom:mask");
image = atom->image = (int *)
memory->srealloc(atom->image,nmax*sizeof(int),"atom:image");
x = atom->x = memory->grow_2d_double_array(atom->x,nmax,3,"atom:x");
v = atom->v = memory->grow_2d_double_array(atom->v,nmax,3,"atom:v");
f = atom->f = memory->grow_2d_double_array(atom->f,nmax,3,"atom:f");
omega = atom->omega =
memory->grow_2d_double_array(atom->omega,nmax,3,"atom:omega");
torque = atom->torque =
memory->grow_2d_double_array(atom->torque,nmax,3,"atom:torque");
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
}
/* ---------------------------------------------------------------------- */
void AtomVecColloid::copy(int i, int j)
{
tag[j] = tag[i];
type[j] = type[i];
mask[j] = mask[i];
image[j] = image[i];
x[j][0] = x[i][0];
x[j][1] = x[i][1];
x[j][2] = x[i][2];
v[j][0] = v[i][0];
v[j][1] = v[i][1];
v[j][2] = v[i][2];
omega[j][0] = omega[i][0];
omega[j][1] = omega[i][1];
omega[j][2] = omega[i][2];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j);
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::pack_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
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];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
dz = pbc[2]*domain->zprd;
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
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;
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::pack_comm_one(int i, double *buf)
{
buf[0] = v[i][0];
buf[1] = v[i][1];
buf[2] = v[i][2];
buf[3] = omega[i][0];
buf[4] = omega[i][1];
buf[5] = omega[i][2];
return 6;
}
/* ---------------------------------------------------------------------- */
void AtomVecColloid::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
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++];
}
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::unpack_comm_one(int i, double *buf)
{
v[i][0] = buf[0];
v[i][1] = buf[1];
v[i][2] = buf[2];
omega[i][0] = buf[3];
omega[i][1] = buf[4];
omega[i][2] = buf[5];
return 6;
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::pack_reverse(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = f[i][0];
buf[m++] = f[i][1];
buf[m++] = f[i][2];
buf[m++] = torque[i][0];
buf[m++] = torque[i][1];
buf[m++] = torque[i][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::pack_reverse_one(int i, double *buf)
{
buf[0] = torque[i][0];
buf[1] = torque[i][1];
buf[2] = torque[i][2];
return 3;
}
/* ---------------------------------------------------------------------- */
void AtomVecColloid::unpack_reverse(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
f[j][0] += buf[m++];
f[j][1] += buf[m++];
f[j][2] += buf[m++];
torque[j][0] += buf[m++];
torque[j][1] += buf[m++];
torque[j][2] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::unpack_reverse_one(int i, double *buf)
{
torque[i][0] += buf[0];
torque[i][1] += buf[1];
torque[i][2] += buf[2];
return 3;
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::pack_border(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
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];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0];
dy = pbc[1];
dz = pbc[2];
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
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;
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::pack_border_one(int i, double *buf)
{
buf[0] = v[i][0];
buf[1] = v[i][1];
buf[2] = v[i][2];
buf[3] = omega[i][0];
buf[4] = omega[i][1];
buf[5] = omega[i][2];
return 6;
}
/* ---------------------------------------------------------------------- */
void AtomVecColloid::unpack_border(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (i == nmax) grow(0);
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
tag[i] = static_cast<int> (buf[m++]);
type[i] = static_cast<int> (buf[m++]);
mask[i] = static_cast<int> (buf[m++]);
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++];
}
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::unpack_border_one(int i, double *buf)
{
v[i][0] = buf[0];
v[i][1] = buf[1];
v[i][2] = buf[2];
omega[i][0] = buf[3];
omega[i][1] = buf[4];
omega[i][2] = buf[5];
return 6;
}
/* ----------------------------------------------------------------------
pack data for atom I for sending to another proc
xyz must be 1st 3 values, so comm::exchange() can test on them
------------------------------------------------------------------------- */
int AtomVecColloid::pack_exchange(int i, double *buf)
{
int m = 1;
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
buf[m++] = v[i][0];
buf[m++] = v[i][1];
buf[m++] = v[i][2];
buf[m++] = tag[i];
buf[m++] = type[i];
buf[m++] = mask[i];
buf[m++] = image[i];
buf[m++] = omega[i][0];
buf[m++] = omega[i][1];
buf[m++] = omega[i][2];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
buf[0] = m;
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecColloid::unpack_exchange(double *buf)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
int m = 1;
x[nlocal][0] = buf[m++];
x[nlocal][1] = buf[m++];
x[nlocal][2] = buf[m++];
v[nlocal][0] = buf[m++];
v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++];
tag[nlocal] = static_cast<int> (buf[m++]);
type[nlocal] = static_cast<int> (buf[m++]);
mask[nlocal] = static_cast<int> (buf[m++]);
image[nlocal] = static_cast<int> (buf[m++]);
omega[nlocal][0] = buf[m++];
omega[nlocal][1] = buf[m++];
omega[nlocal][2] = buf[m++];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->
unpack_exchange(nlocal,&buf[m]);
atom->nlocal++;
return m;
}
/* ----------------------------------------------------------------------
size of restart data for all atoms owned by this proc
include extra data stored by fixes
------------------------------------------------------------------------- */
int AtomVecColloid::size_restart()
{
int i;
int nlocal = atom->nlocal;
int n = 14 * nlocal;
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
for (i = 0; i < nlocal; i++)
n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
return n;
}
/* ----------------------------------------------------------------------
pack atom I's data for restart file including extra quantities
xyz must be 1st 3 values, so that read_restart can test on them
molecular types may be negative, but write as positive
------------------------------------------------------------------------- */
int AtomVecColloid::pack_restart(int i, double *buf)
{
int m = 1;
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
buf[m++] = tag[i];
buf[m++] = type[i];
buf[m++] = mask[i];
buf[m++] = image[i];
buf[m++] = 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];
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
buf[0] = m;
return m;
}
/* ----------------------------------------------------------------------
unpack data for one atom from restart file including extra quantities
------------------------------------------------------------------------- */
int AtomVecColloid::unpack_restart(double *buf)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) {
grow(0);
if (atom->nextra_store)
atom->extra = memory->grow_2d_double_array(atom->extra,nmax,
atom->nextra_store,
"atom:extra");
}
int m = 1;
x[nlocal][0] = buf[m++];
x[nlocal][1] = buf[m++];
x[nlocal][2] = buf[m++];
tag[nlocal] = static_cast<int> (buf[m++]);
type[nlocal] = static_cast<int> (buf[m++]);
mask[nlocal] = static_cast<int> (buf[m++]);
image[nlocal] = static_cast<int> (buf[m++]);
v[nlocal][0] = buf[m++];
v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++];
omega[nlocal][0] = buf[m++];
omega[nlocal][1] = buf[m++];
omega[nlocal][2] = buf[m++];
double **extra = atom->extra;
if (atom->nextra_store) {
int size = static_cast<int> (buf[0]) - m;
for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
}
atom->nlocal++;
return m;
}
/* ----------------------------------------------------------------------
create one atom of itype at coord
set other values to defaults
------------------------------------------------------------------------- */
void AtomVecColloid::create_atom(int itype, double *coord)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = 0;
type[nlocal] = itype;
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = (512 << 20) | (512 << 10) | 512;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
omega[nlocal][0] = 0.0;
omega[nlocal][1] = 0.0;
omega[nlocal][2] = 0.0;
atom->nlocal++;
}
/* ----------------------------------------------------------------------
unpack one line from Atoms section of data file
initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecColloid::data_atom(double *coord, int imagetmp, char **values)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = atoi(values[0]);
if (tag[nlocal] <= 0)
error->one("Invalid atom ID in Atoms section of data file");
type[nlocal] = atoi(values[1]);
if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
error->one("Invalid atom type in Atoms section of data file");
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
image[nlocal] = imagetmp;
mask[nlocal] = 1;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
omega[nlocal][0] = 0.0;
omega[nlocal][1] = 0.0;
omega[nlocal][2] = 0.0;
atom->nlocal++;
}
/* ----------------------------------------------------------------------
unpack hybrid quantities from one line in Atoms section of data file
initialize other atom quantities for this sub-style
------------------------------------------------------------------------- */
int AtomVecColloid::data_atom_hybrid(int nlocal, char **values)
{
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
omega[nlocal][0] = 0.0;
omega[nlocal][1] = 0.0;
omega[nlocal][2] = 0.0;
return 0;
}
/* ----------------------------------------------------------------------
unpack one line from Velocities section of data file
------------------------------------------------------------------------- */
void AtomVecColloid::data_vel(int m, char **values)
{
v[m][0] = atof(values[0]);
v[m][1] = atof(values[1]);
v[m][2] = atof(values[2]);
omega[m][0] = atof(values[3]);
omega[m][1] = atof(values[4]);
omega[m][2] = atof(values[5]);
}
/* ----------------------------------------------------------------------
unpack hybrid quantities from one line in Velocities section of data file
------------------------------------------------------------------------- */
int AtomVecColloid::data_vel_hybrid(int m, char **values)
{
omega[m][0] = atof(values[0]);
omega[m][1] = atof(values[1]);
omega[m][2] = atof(values[2]);
return 3;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */
double AtomVecColloid::memory_usage()
{
double bytes = 0.0;
if (atom->memcheck("tag")) bytes += nmax * sizeof(int);
if (atom->memcheck("type")) bytes += nmax * sizeof(int);
if (atom->memcheck("mask")) bytes += nmax * sizeof(int);
if (atom->memcheck("image")) bytes += nmax * sizeof(int);
if (atom->memcheck("x")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("v")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("omega")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("torque")) bytes += nmax*3 * sizeof(double);
return bytes;
}

59
src/COLLOID/atom_vec_colloid.h Normal file
View File

@ -0,0 +1,59 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifndef ATOM_VEC_COLLOID_H
#define ATOM_VEC_COLLOID_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecColloid : public AtomVec {
public:
AtomVecColloid(class LAMMPS *, int, char **);
virtual ~AtomVecColloid() {}
void grow(int);
void copy(int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_one(int, double *);
void unpack_comm(int, int, double *);
int unpack_comm_one(int, double *);
int pack_reverse(int, int, double *);
int pack_reverse_one(int, double *);
void unpack_reverse(int, int *, double *);
int unpack_reverse_one(int, double *);
int pack_border(int, int *, double *, int, int *);
int pack_border_one(int, double *);
void unpack_border(int, int, double *);
int unpack_border_one(int, double *);
int pack_exchange(int, double *);
int unpack_exchange(double *);
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double *);
void data_atom(double *, int, char **);
int data_atom_hybrid(int, char **);
void data_vel(int, char **);
int data_vel_hybrid(int, char **);
double memory_usage();
private:
int *tag,*type,*mask,*image;
double **x,**v,**f;
double **omega,**torque;
};
}
#endif

223
src/COLLOID/fix_wall_colloid.cpp Normal file
View File

@ -0,0 +1,223 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_wall_colloid.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "output.h"
#include "respa.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixWallColloid::FixWallColloid(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg != 9) error->all("Illegal fix wall/colloid command");
scalar_flag = 1;
vector_flag = 1;
size_vector = 3;
scalar_vector_freq = 1;
extscalar = 1;
extvector = 1;
if (strcmp(arg[3],"xlo") == 0) {
dim = 0;
side = -1;
} else if (strcmp(arg[3],"xhi") == 0) {
dim = 0;
side = 1;
} else if (strcmp(arg[3],"ylo") == 0) {
dim = 1;
side = -1;
} else if (strcmp(arg[3],"yhi") == 0) {
dim = 1;
side = 1;
} else if (strcmp(arg[3],"zlo") == 0) {
dim = 2;
side = -1;
} else if (strcmp(arg[3],"zhi") == 0) {
dim = 2;
side = 1;
} else error->all("Illegal fix wall/colloid command");
coord = atof(arg[4]);
epsilon = atof(arg[5]);
sigma = atof(arg[6]);
diam = atof(arg[7]);
cutoff = atof(arg[8]);
coeff1 = -576.0/315.0 * epsilon * pow(sigma,6.0);
coeff2 = -288.0/3.0 * 0.125*diam*diam*diam* epsilon;
coeff3 = 144.0 * epsilon * pow(sigma,6.0)/7560.0;
coeff4 = 144.0 * epsilon/6.0;
double rinv = 1.0/cutoff;
double r2inv = rinv*rinv;
double r4inv = r2inv*r2inv;
offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
if (dim == 0 && domain->xperiodic)
error->all("Cannot use wall in periodic dimension");
if (dim == 1 && domain->yperiodic)
error->all("Cannot use wall in periodic dimension");
if (dim == 2 && domain->zperiodic)
error->all("Cannot use wall in periodic dimension");
wall_flag = 0;
wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
}
/* ---------------------------------------------------------------------- */
int FixWallColloid::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= THERMO_ENERGY;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixWallColloid::init()
{
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixWallColloid::setup(int vflag)
{
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ---------------------------------------------------------------------- */
void FixWallColloid::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixWallColloid::post_force(int vflag)
{
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double delta,delta2,rinv,r2inv,r4inv,r8inv,fwall;
double r2,rinv2,r2inv2,r4inv2,r6inv2;
double r3,rinv3,r2inv3,r4inv3,r6inv3;
double rad,rad2,rad4,rad8;
wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
wall_flag = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (side == -1) delta = x[i][dim] - coord;
else delta = coord - x[i][dim];
if (delta <= 0.0) continue;
if (delta > cutoff) continue;
rad = 0.5*diam;
rad2 = rad*rad;
rad4 = rad2*rad2;
rad8 = rad4*rad4;
delta2 = rad2 - delta*delta;
rinv = 1.0/delta2;
r2inv = rinv*rinv;
r4inv = r2inv*r2inv;
r8inv = r4inv*r4inv;
fwall = (coeff1*(rad8*rad + 27.0*rad4*rad2*rad*pow(delta,2.0)
+ 63.0*rad4*rad*pow(delta,4.0)
+ 21.0*rad2*rad*pow(delta,6.0))*r8inv
- coeff2*r2inv) * side;
f[i][dim] -= fwall;
r2 = 0.5*diam - delta;
rinv2 = 1.0/r2;
r2inv2 = rinv2*rinv2;
r4inv2 = r2inv2*r2inv2;
r6inv2 = r4inv2*r2inv2;
r3 = delta+0.5*diam;
rinv3 = 1.0/r3;
r2inv3 = rinv3*rinv3;
r4inv3 = r2inv3*r2inv3;
r6inv3 = r4inv3*r2inv3;
wall[0] += coeff3*((-3.5*diam+delta)*r4inv2*r2inv2*rinv2
+ (3.5*diam+delta)*r4inv3*r2inv3*rinv3)
- coeff4*((-diam*delta+r2*r3*(log(-r2)-log(r3)))*(-rinv2)*rinv3) - offset;
wall[dim+1] += fwall;
}
}
/* ---------------------------------------------------------------------- */
void FixWallColloid::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixWallColloid::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
energy of wall interaction
------------------------------------------------------------------------- */
double FixWallColloid::compute_scalar()
{
// only sum across procs one time
if (wall_flag == 0) {
MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
wall_flag = 1;
}
return wall_all[0];
}
/* ----------------------------------------------------------------------
components of force on wall
------------------------------------------------------------------------- */
double FixWallColloid::compute_vector(int n)
{
// only sum across procs one time
if (wall_flag == 0) {
MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
wall_flag = 1;
}
return wall_all[n+1];
}

46
src/COLLOID/fix_wall_colloid.h Normal file
View File

@ -0,0 +1,46 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifndef FIX_WALL_COLLOID_H
#define FIX_WALL_COLLOID_H
#include "fix.h"
namespace LAMMPS_NS {
class FixWallColloid : public Fix {
public:
FixWallColloid(class LAMMPS *, int, char **);
~FixWallColloid() {}
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
double compute_scalar();
double compute_vector(int);
private:
int dim,side,thermo_flag,eflag_enable;
double coord,epsilon,sigma,diam,cutoff;
double coeff1,coeff2,coeff3,coeff4,offset;
double wall[4],wall_all[4];
int wall_flag;
int nlevels_respa;
};
}
#endif

186
src/COLLOID/pair_yukawa_colloid.cpp Normal file
View File

@ -0,0 +1,186 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "pair_yukawa_colloid.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairYukawaColloid::PairYukawaColloid(LAMMPS *lmp) : PairYukawa(lmp) {}
/* ---------------------------------------------------------------------- */
void PairYukawaColloid::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair,radi,radj;
double rsq,r2inv,r,rinv,screening,forceyukawa,factor_coul;
int *ilist,*jlist,*numneigh,**firstneigh;
ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
double *special_coul = force->special_coul;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
radi = atom->shape[itype][0];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (j < nall) factor_coul = 1.0;
else {
factor_coul = special_coul[j/nall];
j %= nall;
}
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
radj = atom->shape[jtype][0];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r = sqrt(rsq);
rinv = 1.0/r;
screening = exp(-kappa*(r-(radi+radj)));
forceyukawa = a[itype][jtype] * screening;
fpair = factor_coul*forceyukawa * rinv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) {
ecoul = a[itype][jtype]/kappa * screening - offset[itype][jtype];
ecoul *= factor_coul;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
0.0,ecoul,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_compute();
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairYukawaColloid::init_style()
{
if (!atom->shape)
error->all("Pair yukawa/colloid requires atom attribute shape");
// insure all particle shapes are spherical
for (int i = 1; i <= atom->ntypes; i++)
if ((atom->shape[i][0] != atom->shape[i][1]) ||
(atom->shape[i][0] != atom->shape[i][2]) ||
(atom->shape[i][1] != atom->shape[i][2]))
error->all("Pair yukawa/colloid requires spherical particles");
if (force->newton_pair != 1)
error->all("Pair lubricate requires newton_pair on");
int irequest = neighbor->request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairYukawaColloid::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
a[i][j] = mix_energy(a[i][i],a[j][j],1.0,1.0);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
}
if (offset_flag) {
double radi = atom->shape[i][0];
double radj = atom->shape[j][0];
double screening = exp(-kappa * (cut[i][j] - (radi+radj)));
offset[i][j] = a[i][j]/kappa * screening;
} else offset[i][j] = 0.0;
a[j][i] = a[i][j];
offset[j][i] = offset[i][j];
return cut[i][j];
}
/* ---------------------------------------------------------------------- */
double PairYukawaColloid::single(int i, int j, int itype, int jtype,
double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
double r2inv,r,rinv,screening,forceyukawa,phi,radi,radj;
int *type = atom->type;
radi = atom->shape[itype][0];
radj = atom->shape[jtype][0];
r2inv = 1.0/rsq;
r = sqrt(rsq);
rinv = 1.0/r;
screening = exp(-kappa*(r-(radi+radj)));
forceyukawa = a[itype][jtype] * screening;
fforce = factor_coul*forceyukawa * rinv;
phi = a[itype][jtype]/kappa * screening - offset[itype][jtype];
return factor_coul*phi;
}

33
src/COLLOID/pair_yukawa_colloid.h Normal file
View File

@ -0,0 +1,33 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifndef PAIR_YUKAWA_COLLOID_H
#define PAIR_YUKAWA_COLLOID_H
#include "pair_yukawa.h"
namespace LAMMPS_NS {
class PairYukawaColloid : public PairYukawa {
public:
PairYukawaColloid(class LAMMPS *);
~PairYukawaColloid() {}
void compute(int, int);
void init_style();
double init_one(int, int);
double single(int, int, int, int, double, double, double, double &);
};
}
#endif

18
src/COLLOID/style_colloid.h
View File

@ -11,12 +11,30 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef AtomInclude
#include "atom_vec_colloid.h"
#endif
#ifdef AtomClass
AtomStyle(colloid,AtomVecColloid)
#endif
#ifdef FixInclude
#include "fix_wall_colloid.h"
#endif
#ifdef FixClass
FixStyle(wall/colloid,FixWallColloid)
#endif
#ifdef PairInclude
#include "pair_colloid.h"
#include "pair_lubricate.h"
#include "pair_yukawa_colloid.h"
#endif
#ifdef PairClass
PairStyle(colloid,PairColloid)
PairStyle(lubricate,PairLubricate)
PairStyle(yukawa/colloid,PairYukawaColloid)
#endif

10
src/pair_yukawa.h
View File

@ -21,18 +21,18 @@ namespace LAMMPS_NS {
class PairYukawa : public Pair {
public:
PairYukawa(class LAMMPS *);
~PairYukawa();
void compute(int, int);
virtual ~PairYukawa();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
virtual double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
virtual double single(int, int, int, int, double, double, double, double &);
private:
protected:
double cut_global;
double kappa;
double **cut,**a,**offset;

18
src/style_colloid.h
View File

@ -11,12 +11,30 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef AtomInclude
#include "atom_vec_colloid.h"
#endif
#ifdef AtomClass
AtomStyle(colloid,AtomVecColloid)
#endif
#ifdef FixInclude
#include "fix_wall_colloid.h"
#endif
#ifdef FixClass
FixStyle(wall/colloid,FixWallColloid)
#endif
#ifdef PairInclude
#include "pair_colloid.h"
#include "pair_lubricate.h"
#include "pair_yukawa_colloid.h"
#endif
#ifdef PairClass
PairStyle(colloid,PairColloid)
PairStyle(lubricate,PairLubricate)
PairStyle(yukawa/colloid,PairYukawaColloid)
#endif