diff --git a/src/COLLOID/Install.csh b/src/COLLOID/Install.csh index c7a6fd25f3..504510afbe 100644 --- a/src/COLLOID/Install.csh +++ b/src/COLLOID/Install.csh @@ -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 diff --git a/src/COLLOID/atom_vec_colloid.cpp b/src/COLLOID/atom_vec_colloid.cpp new file mode 100644 index 0000000000..899c870edf --- /dev/null +++ b/src/COLLOID/atom_vec_colloid.cpp @@ -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 (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (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 (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (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 (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (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 (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; +} diff --git a/src/COLLOID/atom_vec_colloid.h b/src/COLLOID/atom_vec_colloid.h new file mode 100644 index 0000000000..aeff7c9088 --- /dev/null +++ b/src/COLLOID/atom_vec_colloid.h @@ -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 diff --git a/src/COLLOID/fix_wall_colloid.cpp b/src/COLLOID/fix_wall_colloid.cpp new file mode 100644 index 0000000000..64fc48fc1a --- /dev/null +++ b/src/COLLOID/fix_wall_colloid.cpp @@ -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]; +} diff --git a/src/COLLOID/fix_wall_colloid.h b/src/COLLOID/fix_wall_colloid.h new file mode 100644 index 0000000000..75e07cac84 --- /dev/null +++ b/src/COLLOID/fix_wall_colloid.h @@ -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 diff --git a/src/COLLOID/pair_yukawa_colloid.cpp b/src/COLLOID/pair_yukawa_colloid.cpp new file mode 100644 index 0000000000..442636794c --- /dev/null +++ b/src/COLLOID/pair_yukawa_colloid.cpp @@ -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; +} diff --git a/src/COLLOID/pair_yukawa_colloid.h b/src/COLLOID/pair_yukawa_colloid.h new file mode 100644 index 0000000000..eacae35a99 --- /dev/null +++ b/src/COLLOID/pair_yukawa_colloid.h @@ -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 diff --git a/src/COLLOID/style_colloid.h b/src/COLLOID/style_colloid.h index 9655340b35..d5864a0c5f 100644 --- a/src/COLLOID/style_colloid.h +++ b/src/COLLOID/style_colloid.h @@ -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 diff --git a/src/pair_yukawa.h b/src/pair_yukawa.h index a8a84fc83f..2d81275212 100644 --- a/src/pair_yukawa.h +++ b/src/pair_yukawa.h @@ -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; diff --git a/src/style_colloid.h b/src/style_colloid.h index 9655340b35..d5864a0c5f 100644 --- a/src/style_colloid.h +++ b/src/style_colloid.h @@ -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