forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6796 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
9eb5aebed9
commit
3a9ba7dd28
|
|
@ -0,0 +1,15 @@
|
|||
# Install/unInstall package files in LAMMPS
|
||||
|
||||
if (test $1 = 1) then
|
||||
|
||||
cp pair_dsmc.cpp ..
|
||||
|
||||
cp pair_dsmc.h ..
|
||||
|
||||
elif (test $1 = 0) then
|
||||
|
||||
rm -f ../pair_dsmc.cpp
|
||||
|
||||
rm -f ../pair_dsmc.h
|
||||
|
||||
fi
|
||||
|
|
@ -0,0 +1,392 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "mpi.h"
|
||||
#include "string.h"
|
||||
#include "stdlib.h"
|
||||
#include "fix_bond_break.h"
|
||||
#include "update.h"
|
||||
#include "respa.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "comm.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "random_mars.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
|
||||
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 6) error->all("Illegal fix bond/break command");
|
||||
|
||||
MPI_Comm_rank(world,&me);
|
||||
|
||||
nevery = atoi(arg[3]);
|
||||
if (nevery <= 0) error->all("Illegal fix bond/break command");
|
||||
|
||||
force_reneighbor = 1;
|
||||
next_reneighbor = -1;
|
||||
vector_flag = 1;
|
||||
size_vector = 2;
|
||||
global_freq = 1;
|
||||
extvector = 0;
|
||||
|
||||
btype = atoi(arg[4]);
|
||||
double cutoff = atof(arg[5]);
|
||||
|
||||
if (btype < 1 || btype > atom->nbondtypes)
|
||||
error->all("Invalid bond type in fix bond/break command");
|
||||
if (cutoff < 0.0) error->all("Illegal fix bond/break command");
|
||||
|
||||
cutsq = cutoff*cutoff;
|
||||
|
||||
// optional keywords
|
||||
|
||||
fraction = 1.0;
|
||||
int seed = 12345;
|
||||
|
||||
int iarg = 6;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"prob") == 0) {
|
||||
if (iarg+3 > narg) error->all("Illegal fix bond/break command");
|
||||
fraction = atof(arg[iarg+1]);
|
||||
seed = atoi(arg[iarg+2]);
|
||||
if (fraction < 0.0 || fraction > 1.0)
|
||||
error->all("Illegal fix bond/break command");
|
||||
if (seed <= 0) error->all("Illegal fix bond/break command");
|
||||
iarg += 3;
|
||||
} else error->all("Illegal fix bond/break command");
|
||||
}
|
||||
|
||||
// error check
|
||||
|
||||
if (atom->molecular == 0)
|
||||
error->all("Cannot use fix bond/break with non-molecular systems");
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
|
||||
random = new RanMars(lmp,seed + me);
|
||||
|
||||
// set comm sizes needed by this fix
|
||||
|
||||
comm_forward = 2;
|
||||
comm_reverse = 2;
|
||||
|
||||
// allocate arrays local to this fix
|
||||
|
||||
nmax = 0;
|
||||
partner = NULL;
|
||||
distsq = NULL;
|
||||
|
||||
// zero out stats
|
||||
|
||||
breakcount = 0;
|
||||
breakcounttotal = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBondBreak::~FixBondBreak()
|
||||
{
|
||||
delete random;
|
||||
|
||||
// delete locally stored arrays
|
||||
|
||||
memory->destroy(partner);
|
||||
memory->destroy(distsq);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBondBreak::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_INTEGRATE;
|
||||
mask |= POST_INTEGRATE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondBreak::init()
|
||||
{
|
||||
// require special bonds = 0,1,1
|
||||
|
||||
int flag = 0;
|
||||
if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 ||
|
||||
force->special_lj[3] != 1.0) flag = 1;
|
||||
if (force->special_coul[1] != 0.0 || force->special_coul[2] != 1.0 ||
|
||||
force->special_coul[3] != 1.0) flag = 1;
|
||||
if (flag) error->all("Fix bond/break requires special_bonds = 0,1,1");
|
||||
|
||||
// warn if angles, dihedrals, impropers are being used
|
||||
|
||||
if (force->angle || force->dihedral || force->improper) {
|
||||
if (me == 0)
|
||||
error->warning("Broken bonds will not alter angles, "
|
||||
"dihedrals, or impropers");
|
||||
}
|
||||
|
||||
if (strstr(update->integrate_style,"respa"))
|
||||
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondBreak::post_integrate()
|
||||
{
|
||||
int i,j,k,m,n,i1,i2,n1,n3,type;
|
||||
double delx,dely,delz,rsq,min,max;
|
||||
int *slist;
|
||||
|
||||
if (update->ntimestep % nevery) return;
|
||||
|
||||
// need updated ghost atom positions
|
||||
|
||||
comm->forward_comm();
|
||||
|
||||
// resize bond partner list and initialize it
|
||||
// probability array overlays distsq array
|
||||
// needs to be atom->nmax in length
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(partner);
|
||||
memory->destroy(distsq);
|
||||
nmax = atom->nmax;
|
||||
memory->create(partner,nmax,"bond/break:partner");
|
||||
memory->create(distsq,nmax,"bond/break:distsq");
|
||||
probability = distsq;
|
||||
}
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
int nall = atom->nlocal + atom->nghost;
|
||||
|
||||
for (i = 0; i < nall; i++) {
|
||||
partner[i] = 0;
|
||||
distsq[i] = 0.0;
|
||||
}
|
||||
|
||||
// loop over bond list
|
||||
// setup possible partner list of bonds to break
|
||||
|
||||
double **x = atom->x;
|
||||
int *tag = atom->tag;
|
||||
int *mask = atom->mask;
|
||||
int **bondlist = neighbor->bondlist;
|
||||
int nbondlist = neighbor->nbondlist;
|
||||
|
||||
for (n = 0; n < nbondlist; n++) {
|
||||
i1 = bondlist[n][0];
|
||||
i2 = bondlist[n][1];
|
||||
type = bondlist[n][2];
|
||||
if (!(mask[i1] & groupbit)) continue;
|
||||
if (!(mask[i2] & groupbit)) continue;
|
||||
if (type != btype) continue;
|
||||
|
||||
delx = x[i1][0] - x[i2][0];
|
||||
dely = x[i1][1] - x[i2][1];
|
||||
delz = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
if (rsq <= cutsq) continue;
|
||||
|
||||
if (rsq > distsq[i1]) {
|
||||
partner[i1] = tag[i2];
|
||||
distsq[i1] = rsq;
|
||||
}
|
||||
if (rsq > distsq[i2]) {
|
||||
partner[i2] = tag[i1];
|
||||
distsq[i2] = rsq;
|
||||
}
|
||||
}
|
||||
|
||||
// reverse comm of partner info
|
||||
|
||||
if (force->newton_bond) comm->reverse_comm_fix(this);
|
||||
|
||||
// each atom now knows its winning partner
|
||||
// for prob check, generate random value for each atom with a bond partner
|
||||
// forward comm of partner and random value, so ghosts have it
|
||||
|
||||
if (fraction < 1.0) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (partner[i]) probability[i] = random->uniform();
|
||||
}
|
||||
|
||||
comm->forward_comm_fix(this);
|
||||
|
||||
// break bonds
|
||||
// if both atoms list each other as winning bond partner
|
||||
// and probability constraint is satisfied
|
||||
|
||||
int **bond_type = atom->bond_type;
|
||||
int **bond_atom = atom->bond_atom;
|
||||
int *num_bond = atom->num_bond;
|
||||
int **nspecial = atom->nspecial;
|
||||
int **special = atom->special;
|
||||
|
||||
int nbreak = 0;
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (partner[i] == 0) continue;
|
||||
j = atom->map(partner[i]);
|
||||
if (partner[j] != tag[i]) continue;
|
||||
|
||||
// apply probability constraint
|
||||
// MIN,MAX insures values are added in same order on different procs
|
||||
|
||||
if (fraction < 1.0) {
|
||||
min = MIN(probability[i],probability[j]);
|
||||
max = MAX(probability[i],probability[j]);
|
||||
if (0.5*(min+max) >= fraction) continue;
|
||||
}
|
||||
|
||||
// delete bond from atom I if I stores it
|
||||
// atom J will also do this
|
||||
|
||||
for (m = 0; m < num_bond[i]; m++) {
|
||||
if (bond_atom[i][m] == partner[i]) {
|
||||
for (k = m; k < num_bond[i]-1; k++) {
|
||||
bond_atom[i][k] = bond_atom[i][k+1];
|
||||
bond_type[i][k] = bond_type[i][k+1];
|
||||
}
|
||||
num_bond[i]--;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// remove J from special bond list for atom I
|
||||
// atom J will also do this
|
||||
|
||||
slist = special[i];
|
||||
n1 = nspecial[i][0];
|
||||
n3 = nspecial[i][2];
|
||||
for (m = 0; m < n1; m++)
|
||||
if (slist[m] == partner[i]) break;
|
||||
for (; m < n3-1; m++) slist[m] = slist[m+1];
|
||||
nspecial[i][0]--;
|
||||
nspecial[i][1]--;
|
||||
nspecial[i][2]--;
|
||||
|
||||
// count the broken bond once
|
||||
|
||||
if (tag[i] < tag[j]) nbreak++;
|
||||
}
|
||||
|
||||
// tally stats
|
||||
|
||||
MPI_Allreduce(&nbreak,&breakcount,1,MPI_INT,MPI_SUM,world);
|
||||
breakcounttotal += breakcount;
|
||||
atom->nbonds -= breakcount;
|
||||
|
||||
// trigger reneighboring if any bonds were formed
|
||||
|
||||
if (breakcount) next_reneighbor = update->ntimestep;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondBreak::post_integrate_respa(int ilevel, int iloop)
|
||||
{
|
||||
if (ilevel == nlevels_respa-1) post_integrate();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBondBreak::pack_comm(int n, int *list, double *buf,
|
||||
int pbc_flag, int *pbc)
|
||||
{
|
||||
int i,j,m;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
buf[m++] = partner[j];
|
||||
buf[m++] = probability[j];
|
||||
}
|
||||
return 2;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondBreak::unpack_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i,m,last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
partner[i] = static_cast<int> (buf[m++]);
|
||||
probability[i] = buf[m++];
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBondBreak::pack_reverse_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i,m,last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
buf[m++] = partner[i];
|
||||
buf[m++] = distsq[i];
|
||||
}
|
||||
return 2;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondBreak::unpack_reverse_comm(int n, int *list, double *buf)
|
||||
{
|
||||
int i,j,m;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
if (buf[m+1] > distsq[j]) {
|
||||
partner[j] = static_cast<int> (buf[m++]);
|
||||
distsq[j] = buf[m++];
|
||||
} else m += 2;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixBondBreak::compute_vector(int n)
|
||||
{
|
||||
if (n == 1) return (double) breakcount;
|
||||
return (double) breakcounttotal;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixBondBreak::memory_usage()
|
||||
{
|
||||
int nmax = atom->nmax;
|
||||
double bytes = nmax * sizeof(int);
|
||||
bytes += nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,60 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 FIX_CLASS
|
||||
|
||||
FixStyle(bond/break,FixBondBreak)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_BOND_BREAK_H
|
||||
#define LMP_FIX_BOND_BREAK_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixBondBreak : public Fix {
|
||||
public:
|
||||
FixBondBreak(class LAMMPS *, int, char **);
|
||||
~FixBondBreak();
|
||||
int setmask();
|
||||
void init();
|
||||
void post_integrate();
|
||||
void post_integrate_respa(int,int);
|
||||
|
||||
int pack_comm(int, int *, double *, int, int *);
|
||||
void unpack_comm(int, int, double *);
|
||||
int pack_reverse_comm(int, int, double *);
|
||||
void unpack_reverse_comm(int, int *, double *);
|
||||
double compute_vector(int);
|
||||
double memory_usage();
|
||||
|
||||
private:
|
||||
int me;
|
||||
int btype,seed;
|
||||
double cutsq,fraction;
|
||||
|
||||
int breakcount,breakcounttotal;
|
||||
int nmax;
|
||||
int *partner;
|
||||
double *distsq,*probability;
|
||||
|
||||
class RanMars *random;
|
||||
int nlevels_respa;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,617 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 "mpi.h"
|
||||
#include "string.h"
|
||||
#include "stdlib.h"
|
||||
#include "fix_bond_create.h"
|
||||
#include "update.h"
|
||||
#include "respa.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "comm.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "random_mars.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define BIG 1.0e20
|
||||
|
||||
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
|
||||
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBondCreate::FixBondCreate(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 8) error->all("Illegal fix bond/create command");
|
||||
|
||||
MPI_Comm_rank(world,&me);
|
||||
|
||||
nevery = atoi(arg[3]);
|
||||
if (nevery <= 0) error->all("Illegal fix bond/create command");
|
||||
|
||||
force_reneighbor = 1;
|
||||
next_reneighbor = -1;
|
||||
vector_flag = 1;
|
||||
size_vector = 2;
|
||||
global_freq = 1;
|
||||
extvector = 0;
|
||||
|
||||
iatomtype = atoi(arg[4]);
|
||||
jatomtype = atoi(arg[5]);
|
||||
double cutoff = atof(arg[6]);
|
||||
btype = atoi(arg[7]);
|
||||
|
||||
if (iatomtype < 1 || iatomtype > atom->ntypes ||
|
||||
jatomtype < 1 || jatomtype > atom->ntypes)
|
||||
error->all("Invalid atom type in fix bond/create command");
|
||||
if (cutoff < 0.0) error->all("Illegal fix bond/create command");
|
||||
if (btype < 1 || btype > atom->nbondtypes)
|
||||
error->all("Invalid bond type in fix bond/create command");
|
||||
|
||||
cutsq = cutoff*cutoff;
|
||||
|
||||
// optional keywords
|
||||
|
||||
imaxbond = 0;
|
||||
inewtype = iatomtype;
|
||||
jmaxbond = 0;
|
||||
jnewtype = jatomtype;
|
||||
fraction = 1.0;
|
||||
int seed = 12345;
|
||||
|
||||
int iarg = 8;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"iparam") == 0) {
|
||||
if (iarg+3 > narg) error->all("Illegal fix bond/create command");
|
||||
imaxbond = atoi(arg[iarg+1]);
|
||||
inewtype = atoi(arg[iarg+2]);
|
||||
if (imaxbond < 0) error->all("Illegal fix bond/create command");
|
||||
if (inewtype < 1 || inewtype > atom->ntypes)
|
||||
error->all("Invalid atom type in fix bond/create command");
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg],"jparam") == 0) {
|
||||
if (iarg+3 > narg) error->all("Illegal fix bond/create command");
|
||||
jmaxbond = atoi(arg[iarg+1]);
|
||||
jnewtype = atoi(arg[iarg+2]);
|
||||
if (jmaxbond < 0) error->all("Illegal fix bond/create command");
|
||||
if (jnewtype < 1 || jnewtype > atom->ntypes)
|
||||
error->all("Invalid atom type in fix bond/create command");
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg],"prob") == 0) {
|
||||
if (iarg+3 > narg) error->all("Illegal fix bond/create command");
|
||||
fraction = atof(arg[iarg+1]);
|
||||
seed = atoi(arg[iarg+2]);
|
||||
if (fraction < 0.0 || fraction > 1.0)
|
||||
error->all("Illegal fix bond/create command");
|
||||
if (seed <= 0) error->all("Illegal fix bond/create command");
|
||||
iarg += 3;
|
||||
} else error->all("Illegal fix bond/create command");
|
||||
}
|
||||
|
||||
// error check
|
||||
|
||||
if (atom->molecular == 0)
|
||||
error->all("Cannot use fix bond/create with non-molecular systems");
|
||||
if (iatomtype == jatomtype &&
|
||||
((imaxbond != jmaxbond) || (inewtype != jnewtype)))
|
||||
error->all("Inconsistent iparam/jparam values in fix bond/create command");
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
|
||||
random = new RanMars(lmp,seed + me);
|
||||
|
||||
// perform initial allocation of atom-based arrays
|
||||
// register with Atom class
|
||||
// bondcount values will be initialized in setup()
|
||||
|
||||
bondcount = NULL;
|
||||
grow_arrays(atom->nmax);
|
||||
atom->add_callback(0);
|
||||
countflag = 0;
|
||||
|
||||
// set comm sizes needed by this fix
|
||||
|
||||
comm_forward = 2;
|
||||
comm_reverse = 2;
|
||||
|
||||
// allocate arrays local to this fix
|
||||
|
||||
nmax = 0;
|
||||
partner = NULL;
|
||||
distsq = NULL;
|
||||
|
||||
// zero out stats
|
||||
|
||||
createcount = 0;
|
||||
createcounttotal = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBondCreate::~FixBondCreate()
|
||||
{
|
||||
// unregister callbacks to this fix from Atom class
|
||||
|
||||
atom->delete_callback(id,0);
|
||||
|
||||
delete random;
|
||||
|
||||
// delete locally stored arrays
|
||||
|
||||
memory->destroy(bondcount);
|
||||
memory->destroy(partner);
|
||||
memory->destroy(distsq);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBondCreate::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_INTEGRATE;
|
||||
mask |= POST_INTEGRATE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::init()
|
||||
{
|
||||
// check cutoff for iatomtype,jatomtype
|
||||
|
||||
if (force->pair == NULL || cutsq > force->pair->cutsq[iatomtype][jatomtype])
|
||||
error->all("Fix bond/create cutoff is longer than pairwise cutoff");
|
||||
|
||||
// require special bonds = 0,1,1
|
||||
|
||||
if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 ||
|
||||
force->special_lj[3] != 1.0)
|
||||
error->all("Fix bond/create requires special_bonds lj = 0,1,1");
|
||||
|
||||
if (atom->q_flag)
|
||||
if (force->special_coul[1] != 0.0 || force->special_coul[2] != 1.0 ||
|
||||
force->special_coul[3] != 1.0)
|
||||
error->all("Fix bond/create requires special_bonds coul = 0,1,1");
|
||||
|
||||
// warn if angles, dihedrals, impropers are being used
|
||||
|
||||
if (force->angle || force->dihedral || force->improper) {
|
||||
if (me == 0)
|
||||
error->warning("Created bonds will not create angles, "
|
||||
"dihedrals, or impropers");
|
||||
}
|
||||
|
||||
// need a half neighbor list, built when ever re-neighboring occurs
|
||||
|
||||
int irequest = neighbor->request((void *) this);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->fix = 1;
|
||||
|
||||
if (strstr(update->integrate_style,"respa"))
|
||||
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::init_list(int id, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::setup(int vflag)
|
||||
{
|
||||
int i,j,m;
|
||||
|
||||
// compute initial bondcount if this is first run
|
||||
// can't do this earlier, like in constructor or init, b/c need ghost info
|
||||
|
||||
if (countflag) return;
|
||||
countflag = 1;
|
||||
|
||||
// count bonds stored with each bond I own
|
||||
// if newton bond is not set, just increment count on atom I
|
||||
// if newton bond is set, also increment count on atom J even if ghost
|
||||
// bondcount is long enough to tally ghost atom counts
|
||||
|
||||
int *num_bond = atom->num_bond;
|
||||
int **bond_type = atom->bond_type;
|
||||
int **bond_atom = atom->bond_atom;
|
||||
int nlocal = atom->nlocal;
|
||||
int nghost = atom->nghost;
|
||||
int nall = nlocal + nghost;
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
for (i = 0; i < nall; i++) bondcount[i] = 0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
for (j = 0; j < num_bond[i]; j++) {
|
||||
if (bond_type[i][j] == btype) {
|
||||
bondcount[i]++;
|
||||
if (newton_bond) {
|
||||
m = atom->map(bond_atom[i][j]);
|
||||
if (m < 0)
|
||||
error->one("Could not count initial bonds in fix bond/create");
|
||||
bondcount[m]++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// if newton_bond is set, need to sum bondcount
|
||||
|
||||
commflag = 0;
|
||||
if (newton_bond) comm->reverse_comm_fix(this);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::post_integrate()
|
||||
{
|
||||
int i,j,m,ii,jj,inum,jnum,itype,jtype,n1,n3,possible;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,min,max;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh,*slist;
|
||||
|
||||
if (update->ntimestep % nevery) return;
|
||||
|
||||
// need updated ghost atom positions
|
||||
|
||||
comm->forward_comm();
|
||||
|
||||
// forward comm of bondcount, so ghosts have it
|
||||
|
||||
commflag = 0;
|
||||
comm->forward_comm_fix(this);
|
||||
|
||||
// resize bond partner list and initialize it
|
||||
// probability array overlays distsq array
|
||||
// needs to be atom->nmax in length
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(partner);
|
||||
memory->destroy(distsq);
|
||||
nmax = atom->nmax;
|
||||
memory->create(partner,nmax,"bond/create:partner");
|
||||
memory->create(distsq,nmax,"bond/create:distsq");
|
||||
probability = distsq;
|
||||
}
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
int nall = atom->nlocal + atom->nghost;
|
||||
|
||||
for (i = 0; i < nall; i++) {
|
||||
partner[i] = 0;
|
||||
distsq[i] = BIG;
|
||||
}
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
// each atom sets one closest eligible partner atom ID to bond with
|
||||
|
||||
double **x = atom->x;
|
||||
int *tag = atom->tag;
|
||||
int *mask = atom->mask;
|
||||
int *type = atom->type;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
itype = type[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
if (!(mask[j] & groupbit)) continue;
|
||||
jtype = type[j];
|
||||
|
||||
possible = 0;
|
||||
if (itype == iatomtype && jtype == jatomtype) {
|
||||
if ((imaxbond == 0 || bondcount[i] < imaxbond) &&
|
||||
(jmaxbond == 0 || bondcount[j] < jmaxbond))
|
||||
possible = 1;
|
||||
} else if (itype == jatomtype && jtype == iatomtype) {
|
||||
if ((jmaxbond == 0 || bondcount[i] < jmaxbond) &&
|
||||
(imaxbond == 0 || bondcount[j] < imaxbond))
|
||||
possible = 1;
|
||||
}
|
||||
if (!possible) continue;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
if (rsq >= cutsq) continue;
|
||||
|
||||
if (rsq < distsq[i]) {
|
||||
partner[i] = tag[j];
|
||||
distsq[i] = rsq;
|
||||
}
|
||||
if (rsq < distsq[j]) {
|
||||
partner[j] = tag[i];
|
||||
distsq[j] = rsq;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// reverse comm of distsq and partner
|
||||
// not needed if newton_pair off since I,J pair was seen by both procs
|
||||
|
||||
commflag = 1;
|
||||
if (force->newton_pair) comm->reverse_comm_fix(this);
|
||||
|
||||
// each atom now knows its winning partner
|
||||
// for prob check, generate random value for each atom with a bond partner
|
||||
// forward comm of partner and random value, so ghosts have it
|
||||
|
||||
if (fraction < 1.0) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (partner[i]) probability[i] = random->uniform();
|
||||
}
|
||||
|
||||
commflag = 1;
|
||||
comm->forward_comm_fix(this);
|
||||
|
||||
// create bonds for atoms I own
|
||||
// if other atom is owned by another proc, it should create same bond
|
||||
// if both atoms list each other as winning bond partner
|
||||
// and probability constraint is satisfied
|
||||
|
||||
int **bond_type = atom->bond_type;
|
||||
int **bond_atom = atom->bond_atom;
|
||||
int *num_bond = atom->num_bond;
|
||||
int **nspecial = atom->nspecial;
|
||||
int **special = atom->special;
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
int ncreate = 0;
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (partner[i] == 0) continue;
|
||||
j = atom->map(partner[i]);
|
||||
if (partner[j] != tag[i]) continue;
|
||||
|
||||
// apply probability constraint
|
||||
// MIN,MAX insures values are added in same order on different procs
|
||||
|
||||
if (fraction < 1.0) {
|
||||
min = MIN(probability[i],probability[j]);
|
||||
max = MAX(probability[i],probability[j]);
|
||||
if (0.5*(min+max) >= fraction) continue;
|
||||
}
|
||||
|
||||
// if newton_bond is set, only store with I or J
|
||||
// if not newton_bond, store bond with both I and J
|
||||
|
||||
if (!newton_bond || tag[i] < tag[j]) {
|
||||
if (num_bond[i] == atom->bond_per_atom)
|
||||
error->one("New bond exceeded bonds per atom in fix bond/create");
|
||||
bond_type[i][num_bond[i]] = btype;
|
||||
bond_atom[i][num_bond[i]] = tag[j];
|
||||
num_bond[i]++;
|
||||
}
|
||||
|
||||
// add a 1-2 neighbor to special bond list for atom I
|
||||
// atom J will also do this
|
||||
|
||||
slist = special[i];
|
||||
n1 = nspecial[i][0];
|
||||
n3 = nspecial[i][2];
|
||||
if (n3 == atom->maxspecial)
|
||||
error->one("New bond exceeded special list size in fix bond/create");
|
||||
for (m = n3; m > n1; m--) slist[m+1] = slist[m];
|
||||
slist[n1] = tag[j];
|
||||
nspecial[i][0]++;
|
||||
nspecial[i][1]++;
|
||||
nspecial[i][2]++;
|
||||
|
||||
// increment bondcount, convert atom to new type if limit reached
|
||||
|
||||
bondcount[i]++;
|
||||
if (type[i] == iatomtype) {
|
||||
if (bondcount[i] == imaxbond) type[i] = inewtype;
|
||||
} else {
|
||||
if (bondcount[i] == jmaxbond) type[i] = jnewtype;
|
||||
}
|
||||
|
||||
// count the created bond once
|
||||
|
||||
if (tag[i] < tag[j]) ncreate++;
|
||||
}
|
||||
|
||||
// tally stats
|
||||
|
||||
MPI_Allreduce(&ncreate,&createcount,1,MPI_INT,MPI_SUM,world);
|
||||
createcounttotal += createcount;
|
||||
atom->nbonds += createcount;
|
||||
|
||||
// trigger reneighboring if any bonds were formed
|
||||
|
||||
if (createcount) next_reneighbor = update->ntimestep;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::post_integrate_respa(int ilevel, int iloop)
|
||||
{
|
||||
if (ilevel == nlevels_respa-1) post_integrate();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBondCreate::pack_comm(int n, int *list, double *buf,
|
||||
int pbc_flag, int *pbc)
|
||||
{
|
||||
int i,j,m;
|
||||
|
||||
m = 0;
|
||||
|
||||
if (commflag == 0) {
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
buf[m++] = bondcount[j];
|
||||
}
|
||||
return 1;
|
||||
|
||||
} else {
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
buf[m++] = partner[j];
|
||||
buf[m++] = probability[j];
|
||||
}
|
||||
return 2;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::unpack_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i,m,last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
|
||||
if (commflag == 0) {
|
||||
for (i = first; i < last; i++)
|
||||
bondcount[i] = static_cast<int> (buf[m++]);
|
||||
|
||||
} else {
|
||||
for (i = first; i < last; i++) {
|
||||
partner[i] = static_cast<int> (buf[m++]);
|
||||
probability[i] = buf[m++];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBondCreate::pack_reverse_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i,m,last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
|
||||
if (commflag == 0) {
|
||||
for (i = first; i < last; i++)
|
||||
buf[m++] = bondcount[i];
|
||||
return 1;
|
||||
|
||||
} else {
|
||||
for (i = first; i < last; i++) {
|
||||
buf[m++] = distsq[i];
|
||||
buf[m++] = partner[i];
|
||||
}
|
||||
return 2;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::unpack_reverse_comm(int n, int *list, double *buf)
|
||||
{
|
||||
int i,j,m;
|
||||
|
||||
m = 0;
|
||||
|
||||
if (commflag == 0) {
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
bondcount[j] += static_cast<int> (buf[m++]);
|
||||
}
|
||||
|
||||
} else {
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
if (buf[m] < distsq[j]) {
|
||||
distsq[j] = buf[m++];
|
||||
partner[j] = static_cast<int> (buf[m++]);
|
||||
} else m += 2;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::grow_arrays(int nmax)
|
||||
{
|
||||
memory->grow(bondcount,nmax,"bond/create:bondcount");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
copy values within local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixBondCreate::copy_arrays(int i, int j)
|
||||
{
|
||||
bondcount[j] = bondcount[i];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
pack values in local atom-based arrays for exchange with another proc
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int FixBondCreate::pack_exchange(int i, double *buf)
|
||||
{
|
||||
buf[0] = bondcount[i];
|
||||
return 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
unpack values in local atom-based arrays from exchange with another proc
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int FixBondCreate::unpack_exchange(int nlocal, double *buf)
|
||||
{
|
||||
bondcount[nlocal] = static_cast<int> (buf[0]);
|
||||
return 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixBondCreate::compute_vector(int n)
|
||||
{
|
||||
if (n == 1) return (double) createcount;
|
||||
return (double) createcounttotal;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixBondCreate::memory_usage()
|
||||
{
|
||||
int nmax = atom->nmax;
|
||||
double bytes = nmax*2 * sizeof(int);
|
||||
bytes += nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,74 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 FIX_CLASS
|
||||
|
||||
FixStyle(bond/create,FixBondCreate)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_BOND_CREATE_H
|
||||
#define LMP_FIX_BOND_CREATE_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixBondCreate : public Fix {
|
||||
public:
|
||||
FixBondCreate(class LAMMPS *, int, char **);
|
||||
~FixBondCreate();
|
||||
int setmask();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void setup(int);
|
||||
void post_integrate();
|
||||
void post_integrate_respa(int, int);
|
||||
|
||||
int pack_comm(int, int *, double *, int, int *);
|
||||
void unpack_comm(int, int, double *);
|
||||
int pack_reverse_comm(int, int, double *);
|
||||
void unpack_reverse_comm(int, int *, double *);
|
||||
void grow_arrays(int);
|
||||
void copy_arrays(int, int);
|
||||
int pack_exchange(int, double *);
|
||||
int unpack_exchange(int, double *);
|
||||
double compute_vector(int);
|
||||
double memory_usage();
|
||||
|
||||
private:
|
||||
int me;
|
||||
int iatomtype,jatomtype;
|
||||
int btype,seed;
|
||||
int imaxbond,jmaxbond;
|
||||
int inewtype,jnewtype;
|
||||
double cutsq,fraction;
|
||||
|
||||
int createcount,createcounttotal; // bond formation stats
|
||||
|
||||
int nmax;
|
||||
int *bondcount; // count of created bonds this atom is part of
|
||||
int *partner; // ID of preferred atom for this atom to bond to
|
||||
double *distsq; // distance to preferred bond partner
|
||||
double *probability; // random # to use in decision to form bond
|
||||
|
||||
class RanMars *random;
|
||||
class NeighList *list;
|
||||
int countflag,commflag;
|
||||
int nlevels_respa;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,689 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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_bond_swap.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "bond.h"
|
||||
#include "angle.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "group.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "modify.h"
|
||||
#include "compute.h"
|
||||
#include "random_mars.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
#include "update.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg != 6) error->all("Illegal fix bond/swap command");
|
||||
|
||||
vector_flag = 1;
|
||||
size_vector = 2;
|
||||
global_freq = 1;
|
||||
extvector = 0;
|
||||
|
||||
fraction = atof(arg[3]);
|
||||
double cutoff = atof(arg[4]);
|
||||
cutsq = cutoff*cutoff;
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
|
||||
int seed = atoi(arg[5]);
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
// create a new compute temp style
|
||||
// id = fix-ID + temp, compute group = fix group
|
||||
|
||||
int n = strlen(id) + 6;
|
||||
id_temp = new char[n];
|
||||
strcpy(id_temp,id);
|
||||
strcat(id_temp,"_temp");
|
||||
|
||||
char **newarg = new char*[3];
|
||||
newarg[0] = id_temp;
|
||||
newarg[1] = "all";
|
||||
newarg[2] = "temp";
|
||||
modify->add_compute(3,newarg);
|
||||
delete [] newarg;
|
||||
tflag = 1;
|
||||
|
||||
// initialize atom list
|
||||
|
||||
nmax = 0;
|
||||
alist = NULL;
|
||||
|
||||
naccept = foursome = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBondSwap::~FixBondSwap()
|
||||
{
|
||||
delete random;
|
||||
|
||||
// delete temperature if fix created it
|
||||
|
||||
if (tflag) modify->delete_compute(id_temp);
|
||||
delete [] id_temp;
|
||||
|
||||
memory->destroy(alist);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBondSwap::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= PRE_NEIGHBOR;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondSwap::init()
|
||||
{
|
||||
// require an atom style with molecule IDs
|
||||
|
||||
if (atom->molecule == NULL)
|
||||
error->all("Must use atom style with molecule IDs with fix bond/swap");
|
||||
|
||||
int icompute = modify->find_compute(id_temp);
|
||||
if (icompute < 0)
|
||||
error->all("Temperature ID for fix bond/swap does not exist");
|
||||
temperature = modify->compute[icompute];
|
||||
|
||||
// pair and bonds must be defined
|
||||
// no dihedral or improper potentials allowed
|
||||
// special bonds must be 0 1 1
|
||||
|
||||
if (force->pair == NULL || force->bond == NULL)
|
||||
error->all("Fix bond/swap requires pair and bond styles");
|
||||
|
||||
if (force->pair->single_enable == 0)
|
||||
error->all("Pair style does not support fix bond/swap");
|
||||
|
||||
if (force->angle == NULL && atom->nangles > 0 && comm->me == 0)
|
||||
error->warning("Fix bond/swap will ignore defined angles");
|
||||
|
||||
if (force->dihedral || force->improper)
|
||||
error->all("Fix bond/swap cannot use dihedral or improper styles");
|
||||
|
||||
if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 ||
|
||||
force->special_lj[3] != 1.0)
|
||||
error->all("Fix bond/swap requires special_bonds = 0,1,1");
|
||||
|
||||
// need a half neighbor list, built when ever re-neighboring occurs
|
||||
|
||||
int irequest = neighbor->request((void *) this);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->fix = 1;
|
||||
|
||||
// zero out stats
|
||||
|
||||
naccept = foursome = 0;
|
||||
angleflag = 0;
|
||||
if (force->angle) angleflag = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondSwap::init_list(int id, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBondSwap::pre_neighbor()
|
||||
{
|
||||
int i,j,ii,jj,m,inum,jnum;
|
||||
int inext,iprev,ilast,jnext,jprev,jlast,ibond,iangle,jbond,jangle;
|
||||
int itag,inexttag,iprevtag,ilasttag,jtag,jnexttag,jprevtag,jlasttag;
|
||||
int ibondtype,jbondtype,iangletype,inextangletype,jangletype,jnextangletype;
|
||||
int i1,i2,i3,j1,j2,j3,tmp;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double delta,factor;
|
||||
|
||||
// compute current temp for Boltzmann factor test
|
||||
|
||||
double t_current = temperature->compute_scalar();
|
||||
|
||||
// local ptrs to atom arrays
|
||||
|
||||
int *tag = atom->tag;
|
||||
int *mask = atom->mask;
|
||||
int *molecule = atom->molecule;
|
||||
int *num_bond = atom->num_bond;
|
||||
int **bond_atom = atom->bond_atom;
|
||||
int **bond_type = atom->bond_type;
|
||||
int *num_angle = atom->num_angle;
|
||||
int **angle_atom1 = atom->angle_atom1;
|
||||
int **angle_atom2 = atom->angle_atom2;
|
||||
int **angle_atom3 = atom->angle_atom3;
|
||||
int **angle_type = atom->angle_type;
|
||||
int **nspecial = atom->nspecial;
|
||||
int **special = atom->special;
|
||||
int newton_bond = force->newton_bond;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
type = atom->type;
|
||||
x = atom->x;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// randomize list of my owned atoms that are in fix group
|
||||
// grow atom list if necessary
|
||||
|
||||
if (nlocal > nmax) {
|
||||
memory->destroy(alist);
|
||||
nmax = atom->nmax;
|
||||
memory->create(alist,nmax,"bondswap:alist");
|
||||
}
|
||||
|
||||
int neligible = 0;
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
if (mask[i] & groupbit)
|
||||
alist[neligible++] = i;
|
||||
}
|
||||
|
||||
for (i = 0; i < neligible; i++) {
|
||||
j = static_cast<int> (random->uniform() * neligible);
|
||||
tmp = alist[i];
|
||||
alist[i] = alist[j];
|
||||
alist[j] = tmp;
|
||||
}
|
||||
|
||||
// examine ntest of my eligible atoms for potential swaps
|
||||
// atom i is randomly selected via atom list
|
||||
// look at all j neighbors of atom i
|
||||
// atom j must be on-processor (j < nlocal)
|
||||
// atom j must be in fix group
|
||||
// i and j must be same distance from chain end (mol[i] = mol[j])
|
||||
// NOTE: must use extra parens in if test on mask[j] & groupbit
|
||||
|
||||
int ntest = static_cast<int> (fraction * neligible);
|
||||
int accept = 0;
|
||||
|
||||
for (int itest = 0; itest < ntest; itest++) {
|
||||
i = alist[itest];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
if (j >= nlocal) continue;
|
||||
if ((mask[j] & groupbit) == 0) continue;
|
||||
if (molecule[i] != molecule[j]) continue;
|
||||
|
||||
// look at all bond partners of atoms i and j
|
||||
// use num_bond for this, not special list, so also find bondtypes
|
||||
// inext,jnext = bonded atoms
|
||||
// inext,jnext must be on-processor (inext,jnext < nlocal)
|
||||
// inext,jnext must be same dist from chain end (mol[inext] = mol[jnext])
|
||||
// since swaps may occur between two ends of a single chain, insure
|
||||
// the 4 atoms are unique (no duplicates): inext != jnext, inext != j
|
||||
// all 4 old and new bonds must have length < cutoff
|
||||
|
||||
for (ibond = 0; ibond < num_bond[i]; ibond++) {
|
||||
inext = atom->map(bond_atom[i][ibond]);
|
||||
if (inext >= nlocal || inext < 0) continue;
|
||||
ibondtype = bond_type[i][ibond];
|
||||
|
||||
for (jbond = 0; jbond < num_bond[j]; jbond++) {
|
||||
jnext = atom->map(bond_atom[j][jbond]);
|
||||
if (jnext >= nlocal || jnext < 0) continue;
|
||||
jbondtype = bond_type[j][jbond];
|
||||
|
||||
if (molecule[inext] != molecule[jnext]) continue;
|
||||
if (inext == jnext || inext == j) continue;
|
||||
if (dist_rsq(i,inext) >= cutsq) continue;
|
||||
if (dist_rsq(j,jnext) >= cutsq) continue;
|
||||
if (dist_rsq(i,jnext) >= cutsq) continue;
|
||||
if (dist_rsq(j,inext) >= cutsq) continue;
|
||||
|
||||
// if angles are enabled:
|
||||
// find other atoms i,inext,j,jnext are in angles with
|
||||
// and angletypes: i/j angletype, i/j nextangletype
|
||||
// use num_angle for this, not special list, so also find angletypes
|
||||
// 4 atoms consecutively along 1st chain: iprev,i,inext,ilast
|
||||
// 4 atoms consecutively along 2nd chain: jprev,j,jnext,jlast
|
||||
// prev or last atom can be non-existent at end of chain
|
||||
// set prev/last = -1 in this case
|
||||
// if newton bond = 0, then angles are stored by all 4 atoms
|
||||
// so require that iprev,ilast,jprev,jlast be owned by this proc
|
||||
// so all copies of angles can be updated if a swap takes place
|
||||
|
||||
if (angleflag) {
|
||||
itag = tag[i];
|
||||
inexttag = tag[inext];
|
||||
jtag = tag[j];
|
||||
jnexttag = tag[jnext];
|
||||
|
||||
iprev = -1;
|
||||
for (iangle = 0; iangle < num_angle[i]; iangle++) {
|
||||
i1 = angle_atom1[i][iangle];
|
||||
i2 = angle_atom2[i][iangle];
|
||||
i3 = angle_atom3[i][iangle];
|
||||
if (i2 == itag && i3 == inexttag) iprev = atom->map(i1);
|
||||
else if (i1 == inexttag && i2 == itag) iprev = atom->map(i3);
|
||||
if (iprev >= 0) {
|
||||
iangletype = angle_type[i][iangle];
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (!newton_bond && iprev >= nlocal) continue;
|
||||
|
||||
ilast = -1;
|
||||
for (iangle = 0; iangle < num_angle[inext]; iangle++) {
|
||||
i1 = angle_atom1[inext][iangle];
|
||||
i2 = angle_atom2[inext][iangle];
|
||||
i3 = angle_atom3[inext][iangle];
|
||||
if (i1 == itag && i2 == inexttag) ilast = atom->map(i3);
|
||||
else if (i2 == inexttag && i3 == itag) ilast = atom->map(i1);
|
||||
if (ilast >= 0) {
|
||||
inextangletype = angle_type[inext][iangle];
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (!newton_bond && ilast >= nlocal) continue;
|
||||
|
||||
jprev = -1;
|
||||
for (jangle = 0; jangle < num_angle[j]; jangle++) {
|
||||
j1 = angle_atom1[j][jangle];
|
||||
j2 = angle_atom2[j][jangle];
|
||||
j3 = angle_atom3[j][jangle];
|
||||
if (j2 == jtag && j3 == jnexttag) jprev = atom->map(j1);
|
||||
else if (j1 == jnexttag && j2 == jtag) jprev = atom->map(j3);
|
||||
if (jprev >= 0) {
|
||||
jangletype = angle_type[j][jangle];
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (!newton_bond && jprev >= nlocal) continue;
|
||||
|
||||
jlast = -1;
|
||||
for (jangle = 0; jangle < num_angle[jnext]; jangle++) {
|
||||
j1 = angle_atom1[jnext][jangle];
|
||||
j2 = angle_atom2[jnext][jangle];
|
||||
j3 = angle_atom3[jnext][jangle];
|
||||
if (j1 == jtag && j2 == jnexttag) jlast = atom->map(j3);
|
||||
else if (j2 == jnexttag && j3 == jtag) jlast = atom->map(j1);
|
||||
if (jlast >= 0) {
|
||||
jnextangletype = angle_type[jnext][jangle];
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (!newton_bond && jlast >= nlocal) continue;
|
||||
}
|
||||
|
||||
// valid foursome found between 2 chains:
|
||||
// chains = iprev-i-inext-ilast and jprev-j-jnext-jlast
|
||||
// prev or last values are -1 if do not exist due to end of chain
|
||||
// OK to call angle_eng with -1 atom, since just return 0.0
|
||||
// current energy of foursome =
|
||||
// E_nb(i,j) + E_nb(i,jnext) + E_nb(inext,j) + E_nb(inext,jnext) +
|
||||
// E_bond(i,inext) + E_bond(j,jnext) +
|
||||
// E_angle(iprev,i,inext) + E_angle(i,inext,ilast) +
|
||||
// E_angle(jprev,j,jnext) + E_angle(j,jnext,jlast)
|
||||
// new energy of foursome with swapped bonds =
|
||||
// E_nb(i,j) + E_nb(i,inext) + E_nb(j,jnext) + E_nb(inext,jnext) +
|
||||
// E_bond(i,jnext) + E_bond(j,inext) +
|
||||
// E_angle(iprev,i,jnext) + E_angle(i,jnext,jlast) +
|
||||
// E_angle(jprev,j,inext) + E_angle(j,inext,ilast)
|
||||
// energy delta = add/subtract differing terms between 2 formulas
|
||||
|
||||
foursome++;
|
||||
|
||||
delta = pair_eng(i,inext) + pair_eng(j,jnext) -
|
||||
pair_eng(i,jnext) - pair_eng(inext,j);
|
||||
delta += bond_eng(ibondtype,i,jnext) + bond_eng(jbondtype,j,inext) -
|
||||
bond_eng(ibondtype,i,inext) - bond_eng(jbondtype,j,jnext);
|
||||
if (angleflag)
|
||||
delta += angle_eng(iangletype,iprev,i,jnext) +
|
||||
angle_eng(jnextangletype,i,jnext,jlast) +
|
||||
angle_eng(jangletype,jprev,j,inext) +
|
||||
angle_eng(inextangletype,j,inext,ilast) -
|
||||
angle_eng(iangletype,iprev,i,inext) -
|
||||
angle_eng(inextangletype,i,inext,ilast) -
|
||||
angle_eng(jangletype,jprev,j,jnext) -
|
||||
angle_eng(jnextangletype,j,jnext,jlast);
|
||||
|
||||
// if delta <= 0, accept swap
|
||||
// if delta > 0, compute Boltzmann factor with current temperature
|
||||
// only accept if greater than random value
|
||||
// whether accept or not, exit test loop
|
||||
|
||||
if (delta < 0.0) accept = 1;
|
||||
else {
|
||||
factor = exp(-delta/force->boltz/t_current);
|
||||
if (random->uniform() < factor) accept = 1;
|
||||
}
|
||||
goto done;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
done:
|
||||
if (!accept) return;
|
||||
naccept++;
|
||||
|
||||
// change bond partners of affected atoms
|
||||
// on atom i: bond i-inext changes to i-jnext
|
||||
// on atom j: bond j-jnext changes to j-inext
|
||||
// on atom inext: bond inext-i changes to inext-j
|
||||
// on atom jnext: bond jnext-j changes to jnext-i
|
||||
|
||||
for (ibond = 0; ibond < num_bond[i]; ibond++)
|
||||
if (bond_atom[i][ibond] == tag[inext]) bond_atom[i][ibond] = tag[jnext];
|
||||
for (jbond = 0; jbond < num_bond[j]; jbond++)
|
||||
if (bond_atom[j][jbond] == tag[jnext]) bond_atom[j][jbond] = tag[inext];
|
||||
for (ibond = 0; ibond < num_bond[inext]; ibond++)
|
||||
if (bond_atom[inext][ibond] == tag[i]) bond_atom[inext][ibond] = tag[j];
|
||||
for (jbond = 0; jbond < num_bond[jnext]; jbond++)
|
||||
if (bond_atom[jnext][jbond] == tag[j]) bond_atom[jnext][jbond] = tag[i];
|
||||
|
||||
// set global tags of 4 atoms in bonds
|
||||
|
||||
itag = tag[i];
|
||||
inexttag = tag[inext];
|
||||
|
||||
jtag = tag[j];
|
||||
jnexttag = tag[jnext];
|
||||
|
||||
// change 1st special neighbors of affected atoms: i,j,inext,jnext
|
||||
// don't need to change 2nd/3rd special neighbors for any atom
|
||||
// since special bonds = 0 1 1 means they are never used
|
||||
|
||||
for (m = 0; m < nspecial[i][0]; m++)
|
||||
if (special[i][m] == inexttag) special[i][m] = jnexttag;
|
||||
for (m = 0; m < nspecial[j][0]; m++)
|
||||
if (special[j][m] == jnexttag) special[j][m] = inexttag;
|
||||
for (m = 0; m < nspecial[inext][0]; m++)
|
||||
if (special[inext][m] == itag) special[inext][m] = jtag;
|
||||
for (m = 0; m < nspecial[jnext][0]; m++)
|
||||
if (special[jnext][m] == jtag) special[jnext][m] = itag;
|
||||
|
||||
// done if no angles
|
||||
|
||||
if (!angleflag) return;
|
||||
|
||||
// set global tags of 4 additional atoms in angles, 0 if no angle
|
||||
|
||||
if (iprev >= 0) iprevtag = tag[iprev];
|
||||
else iprevtag = 0;
|
||||
if (ilast >= 0) ilasttag = tag[ilast];
|
||||
else ilasttag = 0;
|
||||
|
||||
if (jprev >= 0) jprevtag = tag[jprev];
|
||||
else jprevtag = 0;
|
||||
if (jlast >= 0) jlasttag = tag[jlast];
|
||||
else jlasttag = 0;
|
||||
|
||||
// change angle partners of affected atoms
|
||||
// must check if each angle is stored as a-b-c or c-b-a
|
||||
// on atom i:
|
||||
// angle iprev-i-inext changes to iprev-i-jnext
|
||||
// angle i-inext-ilast changes to i-jnext-jlast
|
||||
// on atom j:
|
||||
// angle jprev-j-jnext changes to jprev-j-inext
|
||||
// angle j-jnext-jlast changes to j-inext-ilast
|
||||
// on atom inext:
|
||||
// angle iprev-i-inext changes to jprev-j-inext
|
||||
// angle i-inext-ilast changes to j-inext-ilast
|
||||
// on atom jnext:
|
||||
// angle jprev-j-jnext changes to iprev-i-jnext
|
||||
// angle j-jnext-jlast changes to i-jnext-jlast
|
||||
|
||||
for (iangle = 0; iangle < num_angle[i]; iangle++) {
|
||||
i1 = angle_atom1[i][iangle];
|
||||
i2 = angle_atom2[i][iangle];
|
||||
i3 = angle_atom3[i][iangle];
|
||||
|
||||
if (i1 == iprevtag && i2 == itag && i3 == inexttag)
|
||||
angle_atom3[i][iangle] = jnexttag;
|
||||
else if (i1 == inexttag && i2 == itag && i3 == iprevtag)
|
||||
angle_atom1[i][iangle] = jnexttag;
|
||||
else if (i1 == itag && i2 == inexttag && i3 == ilasttag) {
|
||||
angle_atom2[i][iangle] = jnexttag;
|
||||
angle_atom3[i][iangle] = jlasttag;
|
||||
} else if (i1 == ilasttag && i2 == inexttag && i3 == itag) {
|
||||
angle_atom1[i][iangle] = jlasttag;
|
||||
angle_atom2[i][iangle] = jnexttag;
|
||||
}
|
||||
}
|
||||
|
||||
for (jangle = 0; jangle < num_angle[j]; jangle++) {
|
||||
j1 = angle_atom1[j][jangle];
|
||||
j2 = angle_atom2[j][jangle];
|
||||
j3 = angle_atom3[j][jangle];
|
||||
|
||||
if (j1 == jprevtag && j2 == jtag && j3 == jnexttag)
|
||||
angle_atom3[j][jangle] = inexttag;
|
||||
else if (j1 == jnexttag && j2 == jtag && j3 == jprevtag)
|
||||
angle_atom1[j][jangle] = inexttag;
|
||||
else if (j1 == jtag && j2 == jnexttag && j3 == jlasttag) {
|
||||
angle_atom2[j][jangle] = inexttag;
|
||||
angle_atom3[j][jangle] = ilasttag;
|
||||
} else if (j1 == jlasttag && j2 == jnexttag && j3 == jtag) {
|
||||
angle_atom1[j][jangle] = ilasttag;
|
||||
angle_atom2[j][jangle] = inexttag;
|
||||
}
|
||||
}
|
||||
|
||||
for (iangle = 0; iangle < num_angle[inext]; iangle++) {
|
||||
i1 = angle_atom1[inext][iangle];
|
||||
i2 = angle_atom2[inext][iangle];
|
||||
i3 = angle_atom3[inext][iangle];
|
||||
|
||||
if (i1 == iprevtag && i2 == itag && i3 == inexttag) {
|
||||
angle_atom1[inext][iangle] = jprevtag;
|
||||
angle_atom2[inext][iangle] = jtag;
|
||||
} else if (i1 == inexttag && i2 == itag && i3 == iprevtag) {
|
||||
angle_atom2[inext][iangle] = jtag;
|
||||
angle_atom3[inext][iangle] = jprevtag;
|
||||
} else if (i1 == itag && i2 == inexttag && i3 == ilasttag)
|
||||
angle_atom1[inext][iangle] = jtag;
|
||||
else if (i1 == ilasttag && i2 == inexttag && i3 == itag)
|
||||
angle_atom3[inext][iangle] = jtag;
|
||||
}
|
||||
|
||||
for (jangle = 0; jangle < num_angle[jnext]; jangle++) {
|
||||
j1 = angle_atom1[jnext][jangle];
|
||||
j2 = angle_atom2[jnext][jangle];
|
||||
j3 = angle_atom3[jnext][jangle];
|
||||
|
||||
if (j1 == jprevtag && j2 == jtag && j3 == jnexttag) {
|
||||
angle_atom1[jnext][jangle] = iprevtag;
|
||||
angle_atom2[jnext][jangle] = itag;
|
||||
} else if (j1 == jnexttag && j2 == jtag && j3 == jprevtag) {
|
||||
angle_atom2[jnext][jangle] = itag;
|
||||
angle_atom3[jnext][jangle] = iprevtag;
|
||||
} else if (j1 == jtag && j2 == jnexttag && j3 == jlasttag)
|
||||
angle_atom1[jnext][jangle] = itag;
|
||||
else if (j1 == jlasttag && j2 == jnexttag && j3 == jtag)
|
||||
angle_atom3[jnext][jangle] = itag;
|
||||
}
|
||||
|
||||
// done if newton bond set
|
||||
|
||||
if (newton_bond) return;
|
||||
|
||||
// change angles stored by iprev,ilast,jprev,jlast
|
||||
// on atom iprev: angle iprev-i-inext changes to iprev-i-jnext
|
||||
// on atom jprev: angle jprev-j-jnext changes to jprev-j-inext
|
||||
// on atom ilast: angle i-inext-ilast changes to j-inext-ilast
|
||||
// on atom jlast: angle j-jnext-jlast changes to i-jnext-jlast
|
||||
|
||||
for (iangle = 0; iangle < num_angle[iprev]; iangle++) {
|
||||
i1 = angle_atom1[iprev][iangle];
|
||||
i2 = angle_atom2[iprev][iangle];
|
||||
i3 = angle_atom3[iprev][iangle];
|
||||
|
||||
if (i1 == iprevtag && i2 == itag && i3 == inexttag)
|
||||
angle_atom3[iprev][iangle] = jnexttag;
|
||||
else if (i1 == inexttag && i2 == itag && i3 == iprevtag)
|
||||
angle_atom1[iprev][iangle] = jnexttag;
|
||||
}
|
||||
|
||||
for (jangle = 0; jangle < num_angle[jprev]; jangle++) {
|
||||
j1 = angle_atom1[jprev][jangle];
|
||||
j2 = angle_atom2[jprev][jangle];
|
||||
j3 = angle_atom3[jprev][jangle];
|
||||
|
||||
if (j1 == jprevtag && j2 == jtag && j3 == jnexttag)
|
||||
angle_atom3[jprev][jangle] = inexttag;
|
||||
else if (j1 == jnexttag && j2 == jtag && j3 == jprevtag)
|
||||
angle_atom1[jprev][jangle] = inexttag;
|
||||
}
|
||||
|
||||
for (iangle = 0; iangle < num_angle[ilast]; iangle++) {
|
||||
i1 = angle_atom1[ilast][iangle];
|
||||
i2 = angle_atom2[ilast][iangle];
|
||||
i3 = angle_atom3[ilast][iangle];
|
||||
|
||||
if (i1 == itag && i2 == inexttag && i3 == ilasttag)
|
||||
angle_atom1[ilast][iangle] = jtag;
|
||||
else if (i1 == ilasttag && i2 == inexttag && i3 == itag)
|
||||
angle_atom3[ilast][iangle] = jtag;
|
||||
}
|
||||
|
||||
for (jangle = 0; jangle < num_angle[jlast]; jangle++) {
|
||||
j1 = angle_atom1[jlast][jangle];
|
||||
j2 = angle_atom2[jlast][jangle];
|
||||
j3 = angle_atom3[jlast][jangle];
|
||||
|
||||
if (j1 == jtag && j2 == jnexttag && j3 == jlasttag)
|
||||
angle_atom1[jlast][jangle] = itag;
|
||||
else if (j1 == jlasttag && j2 == jnexttag && j3 == jtag)
|
||||
angle_atom3[jlast][jangle] = itag;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBondSwap::modify_param(int narg, char **arg)
|
||||
{
|
||||
if (strcmp(arg[0],"temp") == 0) {
|
||||
if (narg < 2) error->all("Illegal fix_modify command");
|
||||
if (tflag) {
|
||||
modify->delete_compute(id_temp);
|
||||
tflag = 0;
|
||||
}
|
||||
delete [] id_temp;
|
||||
int n = strlen(arg[1]) + 1;
|
||||
id_temp = new char[n];
|
||||
strcpy(id_temp,arg[1]);
|
||||
|
||||
int icompute = modify->find_compute(id_temp);
|
||||
if (icompute < 0) error->all("Could not find fix_modify temperature ID");
|
||||
temperature = modify->compute[icompute];
|
||||
|
||||
if (temperature->tempflag == 0)
|
||||
error->all("Fix_modify temperature ID does not compute temperature");
|
||||
if (temperature->igroup != igroup && comm->me == 0)
|
||||
error->warning("Group for fix_modify temp != fix group");
|
||||
return 2;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute squared distance between atoms I,J
|
||||
must use minimum_image since J was found thru atom->map()
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixBondSwap::dist_rsq(int i, int j)
|
||||
{
|
||||
double delx = x[i][0] - x[j][0];
|
||||
double dely = x[i][1] - x[j][1];
|
||||
double delz = x[i][2] - x[j][2];
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
return (delx*delx + dely*dely + delz*delz);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return pairwise interaction energy between atoms I,J
|
||||
will always be full non-bond interaction, so factors = 1 in single() call
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixBondSwap::pair_eng(int i, int j)
|
||||
{
|
||||
double tmp;
|
||||
double rsq = dist_rsq(i,j);
|
||||
return force->pair->single(i,j,type[i],type[j],rsq,1.0,1.0,tmp);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixBondSwap::bond_eng(int btype, int i, int j)
|
||||
{
|
||||
double rsq = dist_rsq(i,j);
|
||||
return force->bond->single(btype,rsq,i,j);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixBondSwap::angle_eng(int atype, int i, int j, int k)
|
||||
{
|
||||
// test for non-existent angle at end of chain
|
||||
|
||||
if (i == -1 || k == -1) return 0.0;
|
||||
return force->angle->single(atype,i,j,k);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return bond swapping stats
|
||||
n = 1 is # of swaps
|
||||
n = 2 is # of attempted swaps
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixBondSwap::compute_vector(int n)
|
||||
{
|
||||
double one,all;
|
||||
if (n == 0) one = naccept;
|
||||
else one = foursome;
|
||||
MPI_Allreduce(&one,&all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
return all;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of alist
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixBondSwap::memory_usage()
|
||||
{
|
||||
double bytes = nmax * sizeof(int);
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -0,0 +1,63 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 FIX_CLASS
|
||||
|
||||
FixStyle(bond/swap,FixBondSwap)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_BONDSWAP_H
|
||||
#define LMP_FIX_BONDSWAP_H
|
||||
|
||||
#include "fix.h"
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixBondSwap : public Fix {
|
||||
public:
|
||||
FixBondSwap(class LAMMPS *, int, char **);
|
||||
~FixBondSwap();
|
||||
int setmask();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void pre_neighbor();
|
||||
int modify_param(int, char **);
|
||||
double compute_vector(int);
|
||||
double memory_usage();
|
||||
|
||||
private:
|
||||
double fraction,cutsq;
|
||||
int nmax,tflag;
|
||||
int *alist;
|
||||
int naccept,foursome;
|
||||
int angleflag;
|
||||
char *id_temp;
|
||||
int *type;
|
||||
double **x;
|
||||
|
||||
class NeighList *list;
|
||||
class Compute *temperature;
|
||||
class RanMars *random;
|
||||
|
||||
double dist_rsq(int, int);
|
||||
double pair_eng(int, int);
|
||||
double bond_eng(int, int, int);
|
||||
double angle_eng(int, int, int, int);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,562 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 author: Paul Crozier (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "fix_gcmc.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "comm.h"
|
||||
#include "group.h"
|
||||
#include "domain.h"
|
||||
#include "random_park.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 11) error->all("Illegal fix GCMC command");
|
||||
|
||||
vector_flag = 1;
|
||||
size_vector = 6;
|
||||
global_freq = 1;
|
||||
extvector = 0;
|
||||
restart_global = 1;
|
||||
time_depend = 1;
|
||||
|
||||
// required args
|
||||
|
||||
nevery = atoi(arg[3]);
|
||||
nexchanges = atoi(arg[4]);
|
||||
nmcmoves = atoi(arg[5]);
|
||||
ntype = atoi(arg[6]);
|
||||
seed = atoi(arg[7]);
|
||||
reservoir_temperature = atof(arg[8]);
|
||||
chemical_potential = atof(arg[9]);
|
||||
displace = atof(arg[10]);
|
||||
|
||||
if (ntype <= 0 || ntype > atom->ntypes)
|
||||
error->all("Invalid atom type in fix GCMC command");
|
||||
if (nexchanges < 0) error->all("Illegal fix GCMC command");
|
||||
if (nmcmoves < 0) error->all("Illegal fix GCMC command");
|
||||
if (seed <= 0) error->all("Illegal fix GCMC command");
|
||||
if (reservoir_temperature < 0.0) error->all("Illegal fix GCMC command");
|
||||
if (displace < 0.0) error->all("Illegal fix GCMC command");
|
||||
|
||||
// compute beta, lambda, sigma, and the zz factor
|
||||
|
||||
beta = 1.0/(force->boltz*reservoir_temperature);
|
||||
double PI = 4.0*atan(1.0);
|
||||
double gas_mass = atom->mass[ntype];
|
||||
double lambda = sqrt(force->hplanck*force->hplanck/
|
||||
(2.0*PI*gas_mass*force->mvv2e*
|
||||
force->boltz*reservoir_temperature));
|
||||
sigma = sqrt(force->boltz*reservoir_temperature/gas_mass/force->mvv2e);
|
||||
zz = exp(beta*chemical_potential)/(pow(lambda,3));
|
||||
|
||||
// set defaults
|
||||
|
||||
molflag = 0;
|
||||
|
||||
// read options from end of input line
|
||||
|
||||
options(narg-11,&arg[11]);
|
||||
|
||||
// random number generator, same for all procs
|
||||
|
||||
random_equal = new RanPark(lmp,seed);
|
||||
|
||||
// random number generator, not the same for all procs
|
||||
|
||||
random_unequal = new RanPark(lmp,seed);
|
||||
|
||||
// compute the number of MC cycles that occur nevery timesteps
|
||||
|
||||
ncycles = nexchanges + nmcmoves;
|
||||
|
||||
// set up reneighboring
|
||||
|
||||
force_reneighbor = 1;
|
||||
next_reneighbor = update->ntimestep + 1;
|
||||
|
||||
// zero out counters
|
||||
|
||||
nmove_attempts = 0.0;
|
||||
nmove_successes = 0.0;
|
||||
ndel_attempts = 0.0;
|
||||
ndel_successes = 0.0;
|
||||
ninsert_attempts = 0.0;
|
||||
ninsert_successes = 0.0;
|
||||
|
||||
nmax = 0;
|
||||
local_gas_list = NULL;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixGCMC::~FixGCMC()
|
||||
{
|
||||
delete random_equal;
|
||||
delete random_unequal;
|
||||
memory->sfree(local_gas_list);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixGCMC::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= PRE_EXCHANGE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixGCMC::init()
|
||||
{
|
||||
// check that no deletable atoms are in atom->firstgroup
|
||||
// deleting such an atom would not leave firstgroup atoms first
|
||||
|
||||
int *type = atom->type;
|
||||
|
||||
if (atom->firstgroup >= 0) {
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int firstgroupbit = group->bitmask[atom->firstgroup];
|
||||
|
||||
int flag = 0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if ((type[i] == ntype) && (mask[i] && firstgroupbit)) flag = 1;
|
||||
|
||||
int flagall;
|
||||
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
|
||||
|
||||
if (flagall)
|
||||
error->all("Cannot do GCMC on atoms in atom_modify first group");
|
||||
}
|
||||
|
||||
// if molflag not set, warn if any deletable atom has a mol ID
|
||||
|
||||
if (molflag == 0 && atom->molecule_flag) {
|
||||
int *molecule = atom->molecule;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int flag = 0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (type[i] == ntype)
|
||||
if (molecule[i]) flag = 1;
|
||||
int flagall;
|
||||
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
|
||||
if (flagall && comm->me == 0)
|
||||
error->warning("Fix GCMC may delete atom with non-zero molecule ID");
|
||||
}
|
||||
|
||||
if (molflag && atom->molecule_flag == 0)
|
||||
error->all("Fix GCMC molecule command requires atom attribute molecule");
|
||||
|
||||
if (molflag != 0) error->all("Fix GCMC molecule feature does not yet work");
|
||||
|
||||
if (force->pair->single_enable == 0)
|
||||
error->all("Fix GCMC incompatible with given pair_style");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
attempt particle insertions and deletions
|
||||
done before exchange, borders, reneighbor
|
||||
so that ghost atoms and neighbor lists will be correct
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGCMC::pre_exchange()
|
||||
{
|
||||
// just return if should not be called on this timestep
|
||||
|
||||
if (next_reneighbor != update->ntimestep) return;
|
||||
|
||||
if (domain->triclinic == 0) {
|
||||
xlo = domain->boxlo[0];
|
||||
xhi = domain->boxhi[0];
|
||||
ylo = domain->boxlo[1];
|
||||
yhi = domain->boxhi[1];
|
||||
zlo = domain->boxlo[2];
|
||||
zhi = domain->boxhi[2];
|
||||
sublo = domain->sublo;
|
||||
subhi = domain->subhi;
|
||||
} else {
|
||||
xlo = domain->boxlo_bound[0];
|
||||
xhi = domain->boxhi_bound[0];
|
||||
ylo = domain->boxlo_bound[1];
|
||||
yhi = domain->boxhi_bound[1];
|
||||
zlo = domain->boxlo_bound[2];
|
||||
zhi = domain->boxhi_bound[2];
|
||||
sublo = domain->sublo_lamda;
|
||||
subhi = domain->subhi_lamda;
|
||||
}
|
||||
|
||||
volume = domain->xprd * domain->yprd * domain->zprd;
|
||||
|
||||
// grow local_gas_list array if necessary
|
||||
|
||||
if (atom->nlocal > nmax) {
|
||||
memory->sfree(local_gas_list);
|
||||
nmax = atom->nmax;
|
||||
local_gas_list = (int *) memory->smalloc(nmax*sizeof(int),
|
||||
"GCMC:local_gas_list");
|
||||
}
|
||||
|
||||
int *type = atom->type;
|
||||
ngas_local = 0;
|
||||
for (int i = 0; i < atom->nlocal; i++)
|
||||
if (type[i] == ntype)
|
||||
local_gas_list[ngas_local++] = i;
|
||||
|
||||
MPI_Allreduce(&ngas_local,&ngas,1,MPI_INT,MPI_SUM,world);
|
||||
MPI_Scan(&ngas_local,&ngas_before,1,MPI_INT,MPI_SUM,world);
|
||||
ngas_before -= ngas_local;
|
||||
|
||||
// perform ncycles MC cycles
|
||||
|
||||
for (int i = 0; i < ncycles; i++) {
|
||||
int random_int_fraction =
|
||||
static_cast<int>(random_equal->uniform()*ncycles) + 1;
|
||||
if (random_int_fraction <= nmcmoves) {
|
||||
attempt_move();
|
||||
} else {
|
||||
if (random_equal->uniform() < 0.5) attempt_deletion();
|
||||
else attempt_insertion();
|
||||
}
|
||||
}
|
||||
|
||||
next_reneighbor = update->ntimestep + nevery;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
choose particle randomly across all procs and attempt displacement
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGCMC::attempt_move()
|
||||
{
|
||||
int i,iwhichglobal,iwhichlocal;
|
||||
double rx,ry,rz;
|
||||
double coord[3];
|
||||
double **x = atom->x;
|
||||
|
||||
nmove_attempts += 1.0;
|
||||
|
||||
int success = 0;
|
||||
iwhichglobal = static_cast<int> (ngas*random_equal->uniform());
|
||||
if ((iwhichglobal >= ngas_before) &&
|
||||
(iwhichglobal < ngas_before + ngas_local)) {
|
||||
iwhichlocal = iwhichglobal - ngas_before;
|
||||
i = local_gas_list[iwhichlocal];
|
||||
double energy_before = energy(i,x[i]);
|
||||
double rsq = 1.1;
|
||||
while (rsq > 1.0) {
|
||||
rx = 2*random_unequal->uniform() - 1.0;
|
||||
ry = 2*random_unequal->uniform() - 1.0;
|
||||
if (domain->dimension == 3) rz = 2*random_unequal->uniform() - 1.0;
|
||||
else rz = 0.0;
|
||||
rsq = rx*rx + ry*ry + rz*rz;
|
||||
}
|
||||
coord[0] = x[i][0] + displace*rx;
|
||||
coord[1] = x[i][1] + displace*ry;
|
||||
coord[2] = x[i][2] + displace*rz;
|
||||
double energy_after = energy(i,coord);
|
||||
if (random_unequal->uniform() < exp(-beta*(energy_after - energy_before))) {
|
||||
x[i][0] = coord[0];
|
||||
x[i][1] = coord[1];
|
||||
x[i][2] = coord[2];
|
||||
success = 1;
|
||||
}
|
||||
}
|
||||
|
||||
int success_all = 0;
|
||||
MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
|
||||
|
||||
if (success_all) {
|
||||
nmove_successes += 1.0;
|
||||
comm->borders();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
attempt particle deletion
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGCMC::attempt_deletion()
|
||||
{
|
||||
ndel_attempts += 1.0;
|
||||
|
||||
if (ngas == 0) return;
|
||||
|
||||
int i,iwhichglobal,iwhichlocal;
|
||||
AtomVec *avec = atom->avec;
|
||||
|
||||
// choose particle randomly across all procs and delete it
|
||||
// keep ngas, ngas_local, ngas_before, and local_gas_list current
|
||||
// after each deletion
|
||||
|
||||
int success = 0;
|
||||
iwhichglobal = static_cast<int> (ngas*random_equal->uniform());
|
||||
if ((iwhichglobal >= ngas_before) &&
|
||||
(iwhichglobal < ngas_before + ngas_local)) {
|
||||
iwhichlocal = iwhichglobal - ngas_before;
|
||||
i = local_gas_list[iwhichlocal];
|
||||
double deletion_energy = energy(i,atom->x[i]);
|
||||
if (random_unequal->uniform() <
|
||||
ngas*exp(beta*deletion_energy)/(zz*volume)) {
|
||||
avec->copy(atom->nlocal-1,i,1);
|
||||
atom->nlocal--;
|
||||
local_gas_list[iwhichlocal] = local_gas_list[ngas_local-1];
|
||||
ngas_local--;
|
||||
success = 1;
|
||||
}
|
||||
}
|
||||
|
||||
int success_all = 0;
|
||||
MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
|
||||
|
||||
if (success_all) {
|
||||
ngas--;
|
||||
ndel_successes += 1.0;
|
||||
atom->natoms--;
|
||||
if (iwhichglobal < ngas_before) ngas_before--;
|
||||
comm->borders();
|
||||
if (atom->tag_enable) {
|
||||
atom->tag_extend();
|
||||
if (atom->map_style) {
|
||||
atom->map_init();
|
||||
atom->map_set();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
attempt particle insertion
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGCMC::attempt_insertion()
|
||||
{
|
||||
int flag,success;
|
||||
double coord[3],lamda[3];
|
||||
double *newcoord;
|
||||
|
||||
ninsert_attempts += 1.0;
|
||||
|
||||
// choose random position for new atom within box
|
||||
|
||||
coord[0] = xlo + random_equal->uniform() * (xhi-xlo);
|
||||
coord[1] = ylo + random_equal->uniform() * (yhi-ylo);
|
||||
coord[2] = zlo + random_equal->uniform() * (zhi-zlo);
|
||||
|
||||
// check if new atom is in my sub-box or above it if I'm highest proc
|
||||
// if so, add to my list via create_atom()
|
||||
// initialize info about the atoms
|
||||
// set group mask to "all" plus fix group
|
||||
|
||||
if (domain->triclinic) {
|
||||
domain->x2lamda(coord,lamda);
|
||||
newcoord = lamda;
|
||||
} else newcoord = coord;
|
||||
|
||||
flag = 0;
|
||||
if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
|
||||
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] &&
|
||||
newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1;
|
||||
|
||||
success = 0;
|
||||
if (flag) {
|
||||
int nall = atom->nlocal + atom->nghost;
|
||||
double insertion_energy = energy(nall,coord);
|
||||
if (random_unequal->uniform() <
|
||||
zz*volume*exp(-beta*insertion_energy)/(ngas+1)) {
|
||||
atom->avec->create_atom(ntype,coord);
|
||||
int m = atom->nlocal - 1;
|
||||
atom->type[m] = ntype;
|
||||
atom->mask[m] = 1 | groupbit;
|
||||
atom->v[m][0] = random_unequal->gaussian()*sigma;
|
||||
atom->v[m][1] = random_unequal->gaussian()*sigma;
|
||||
atom->v[m][2] = random_unequal->gaussian()*sigma;
|
||||
|
||||
int nfix = modify->nfix;
|
||||
Fix **fix = modify->fix;
|
||||
for (int j = 0; j < nfix; j++)
|
||||
if (fix[j]->create_attribute) fix[j]->set_arrays(m);
|
||||
|
||||
if (atom->nlocal > nmax) {
|
||||
nmax = atom->nmax;
|
||||
local_gas_list = (int *)
|
||||
memory->srealloc(local_gas_list,nmax*sizeof(int),
|
||||
"GCMC:local_gas_list");
|
||||
}
|
||||
|
||||
local_gas_list[ngas_local] = atom->nlocal;
|
||||
ngas_local++;
|
||||
success = 1;
|
||||
}
|
||||
}
|
||||
|
||||
int success_all = 0;
|
||||
MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
|
||||
|
||||
if (success_all) {
|
||||
ngas++;
|
||||
ninsert_successes += 1.0;
|
||||
MPI_Scan(&ngas_local,&ngas_before,1,MPI_INT,MPI_SUM,world);
|
||||
ngas_before -= ngas_local;
|
||||
comm->borders();
|
||||
if (atom->tag_enable) {
|
||||
atom->natoms++;
|
||||
atom->tag_extend();
|
||||
if (atom->map_style) {
|
||||
atom->map_init();
|
||||
atom->map_set();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute particle's interaction energy with the rest of the system
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixGCMC::energy(int i, double *coord)
|
||||
{
|
||||
double delx,dely,delz,rsq;
|
||||
|
||||
double **x = atom->x;
|
||||
int *type = atom->type;
|
||||
int nall = atom->nlocal + atom->nghost;
|
||||
pair = force->pair;
|
||||
cutsq = force->pair->cutsq;
|
||||
|
||||
double fpair = 0.0;
|
||||
double factor_coul = 1.0;
|
||||
double factor_lj = 1.0;
|
||||
|
||||
double total_energy = 0.0;
|
||||
for (int j = 0; j < nall; j++) {
|
||||
|
||||
if (i == j) continue;
|
||||
|
||||
delx = coord[0] - x[j][0];
|
||||
dely = coord[1] - x[j][1];
|
||||
delz = coord[2] - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
int jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[ntype][jtype])
|
||||
total_energy +=
|
||||
pair->single(i,j,ntype,jtype,rsq,factor_coul,factor_lj,fpair);
|
||||
}
|
||||
|
||||
return total_energy;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
parse optional parameters at end of input line
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGCMC::options(int narg, char **arg)
|
||||
{
|
||||
if (narg < 0) error->all("Illegal fix GCMC command");
|
||||
|
||||
int iarg = 0;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"molecule") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal fix GCMC command");
|
||||
if (strcmp(arg[iarg+1],"no") == 0) molflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) molflag = 1;
|
||||
else error->all("Illegal fix evaporate command");
|
||||
iarg += 2;
|
||||
} else error->all("Illegal fix GCMC command");
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return acceptance ratios
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixGCMC::compute_vector(int n)
|
||||
{
|
||||
if (n == 0) return nmove_attempts;
|
||||
if (n == 1) return nmove_successes;
|
||||
if (n == 2) return ndel_attempts;
|
||||
if (n == 3) return ndel_successes;
|
||||
if (n == 4) return ninsert_attempts;
|
||||
if (n == 5) return ninsert_successes;
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixGCMC::memory_usage()
|
||||
{
|
||||
double bytes = nmax * sizeof(int);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
pack entire state of Fix into one write
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGCMC::write_restart(FILE *fp)
|
||||
{
|
||||
int n = 0;
|
||||
double list[4];
|
||||
list[n++] = random_equal->state();
|
||||
list[n++] = random_unequal->state();
|
||||
list[n++] = next_reneighbor;
|
||||
|
||||
if (comm->me == 0) {
|
||||
int size = n * sizeof(double);
|
||||
fwrite(&size,sizeof(int),1,fp);
|
||||
fwrite(list,sizeof(double),n,fp);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
use state info from restart file to restart the Fix
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGCMC::restart(char *buf)
|
||||
{
|
||||
int n = 0;
|
||||
double *list = (double *) buf;
|
||||
|
||||
seed = static_cast<int> (list[n++]);
|
||||
random_equal->reset(seed);
|
||||
|
||||
seed = static_cast<int> (list[n++]);
|
||||
random_unequal->reset(seed);
|
||||
|
||||
next_reneighbor = static_cast<int> (list[n++]);
|
||||
}
|
||||
|
|
@ -0,0 +1,78 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 FIX_CLASS
|
||||
|
||||
FixStyle(gcmc,FixGCMC)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_GCMC_H
|
||||
#define LMP_FIX_GCMC_H
|
||||
|
||||
#include "stdio.h"
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixGCMC : public Fix {
|
||||
public:
|
||||
FixGCMC(class LAMMPS *, int, char **);
|
||||
~FixGCMC();
|
||||
int setmask();
|
||||
void init();
|
||||
void pre_exchange();
|
||||
void attempt_move();
|
||||
void attempt_deletion();
|
||||
void attempt_insertion();
|
||||
double energy(int, double *);
|
||||
double compute_vector(int);
|
||||
double memory_usage();
|
||||
void write_restart(FILE *);
|
||||
void restart(char *);
|
||||
|
||||
private:
|
||||
int ntype,nevery,seed;
|
||||
int ncycles,nexchanges,nmcmoves;
|
||||
int ngas; // # of gas molecules (or atoms) on all procs
|
||||
int ngas_local; // # of gas molecules (or atoms) on this proc
|
||||
int ngas_before; // # of gas molecules (or atoms) on procs < this proc
|
||||
int molflag; // 0 = atomic, 1 = molecular system
|
||||
double nmove_attempts;
|
||||
double nmove_successes;
|
||||
double ndel_attempts;
|
||||
double ndel_successes;
|
||||
double ninsert_attempts;
|
||||
double ninsert_successes;
|
||||
|
||||
int nmax;
|
||||
double reservoir_temperature;
|
||||
double chemical_potential;
|
||||
double displace;
|
||||
double beta,zz,sigma,volume;
|
||||
double xlo,xhi,ylo,yhi,zlo,zhi;
|
||||
double *sublo,*subhi;
|
||||
int *local_gas_list;
|
||||
double **cutsq;
|
||||
class Pair *pair;
|
||||
|
||||
class RanPark *random_equal;
|
||||
class RanPark *random_unequal;
|
||||
|
||||
void options(int, char **);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,525 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Paul Crozier (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "pair_dsmc.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "domain.h"
|
||||
#include "update.h"
|
||||
#include "random_mars.h"
|
||||
#include "limits.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define MIN(a,b) ((a) < (b) ? (a) : (b))
|
||||
#define MAX(a,b) ((a) > (b) ? (a) : (b))
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairDSMC::PairDSMC(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
single_enable = 0;
|
||||
|
||||
total_number_of_collisions = 0;
|
||||
max_particles = max_particle_list = 0;
|
||||
next_particle = NULL;
|
||||
random = NULL;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairDSMC::~PairDSMC()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
memory->destroy(sigma);
|
||||
memory->destroy(cut);
|
||||
memory->destroy(V_sigma_max);
|
||||
memory->destroy(particle_list);
|
||||
memory->destroy(first);
|
||||
memory->destroy(number);
|
||||
}
|
||||
|
||||
delete [] next_particle;
|
||||
delete random;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::compute(int eflag, int vflag)
|
||||
{
|
||||
double **x = atom->x;
|
||||
double *mass = atom->mass;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 1; i <= atom->ntypes; ++i)
|
||||
for (int j = 0; j < total_ncells; ++j) {
|
||||
first[i][j] = -1;
|
||||
number[i][j] = 0;
|
||||
}
|
||||
|
||||
if (atom->nmax > max_particles) {
|
||||
delete [] next_particle;
|
||||
max_particles = atom->nmax;
|
||||
next_particle = new int[max_particles];
|
||||
}
|
||||
|
||||
// find each particle's cell and sort by type
|
||||
// assume a constant volume and shape simulation domain
|
||||
// skip particle if outside processor domain
|
||||
|
||||
for (int i = 0; i < nlocal; ++i) {
|
||||
int xcell = static_cast<int>((x[i][0] - domain->boxlo[0])/cellx);
|
||||
int ycell = static_cast<int>((x[i][1] - domain->boxlo[1])/celly);
|
||||
int zcell = static_cast<int>((x[i][2] - domain->boxlo[2])/cellz);
|
||||
|
||||
if ((xcell < 0) or (xcell > ncellsx-1) or
|
||||
(ycell < 0) or (ycell > ncellsy-1) or
|
||||
(zcell < 0) or (zcell > ncellsz-1)) continue;
|
||||
|
||||
int icell = xcell + ycell*ncellsx + zcell*ncellsx*ncellsy;
|
||||
itype = type[i];
|
||||
next_particle[i] = first[itype][icell];
|
||||
first[itype][icell] = i;
|
||||
number[itype][icell]++;
|
||||
}
|
||||
|
||||
for (int icell = 0; icell < total_ncells; ++icell) {
|
||||
|
||||
for (itype = 1; itype <= atom->ntypes; ++itype) {
|
||||
number_of_A = number[itype][icell];
|
||||
if (number_of_A > max_particle_list) {
|
||||
max_particle_list = number_of_A;
|
||||
memory->grow(particle_list,atom->ntypes+1,max_particle_list,
|
||||
"pair:particle_list");
|
||||
}
|
||||
|
||||
int m = first[itype][icell];
|
||||
for (int k = 0; k < number_of_A; k++) {
|
||||
particle_list[itype][k] = m;
|
||||
m = next_particle[m];
|
||||
}
|
||||
}
|
||||
|
||||
for (itype = 1; itype <= atom->ntypes; ++itype) {
|
||||
imass = mass[itype];
|
||||
number_of_A = number[itype][icell];
|
||||
|
||||
for (jtype = itype; jtype <= atom->ntypes; ++jtype) {
|
||||
jmass = mass[jtype];
|
||||
number_of_B = number[jtype][icell];
|
||||
|
||||
reduced_mass = imass*jmass/(imass + jmass);
|
||||
total_mass = imass + jmass;
|
||||
jmass_tmass = jmass/total_mass;
|
||||
imass_tmass = imass/total_mass;
|
||||
|
||||
// if necessary, recompute V_sigma_max values
|
||||
|
||||
if (recompute_vsigmamax_stride &&
|
||||
(update->ntimestep % recompute_vsigmamax_stride == 0))
|
||||
recompute_V_sigma_max(icell);
|
||||
|
||||
// # of collisions to perform for itype-jtype pairs
|
||||
|
||||
double &Vs_max = V_sigma_max[itype][jtype];
|
||||
double num_of_collisions_double = number_of_A * number_of_B *
|
||||
weighting * Vs_max * update->dt / vol;
|
||||
|
||||
if ((itype == jtype) and number_of_B)
|
||||
num_of_collisions_double *=
|
||||
0.5 * double(number_of_B - 1) / double(number_of_B);
|
||||
|
||||
int num_of_collisions =
|
||||
convert_double_to_equivalent_int(num_of_collisions_double);
|
||||
|
||||
if (num_of_collisions > number_of_A)
|
||||
error->warning("Pair dsmc: num_of_collisions > number_of_A",0);
|
||||
if (num_of_collisions > number_of_B)
|
||||
error->warning("Pair dsmc: num_of_collisions > number_of_B",0);
|
||||
|
||||
// perform collisions on pairs of particles in icell
|
||||
|
||||
for (int k = 0; k < num_of_collisions; k++) {
|
||||
if ((number_of_A < 1) or (number_of_B < 1)) break;
|
||||
int ith_A = static_cast<int>(random->uniform()*number_of_A);
|
||||
int jth_B = static_cast<int>(random->uniform()*number_of_B);
|
||||
int i = particle_list[itype][ith_A];
|
||||
int j = particle_list[jtype][jth_B];
|
||||
if (i == j) {
|
||||
k--;
|
||||
continue;
|
||||
}
|
||||
double probability = V_sigma(i,j)/Vs_max;
|
||||
if (probability > random->uniform()) scatter_random(i,j,icell);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::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(sigma,n+1,n+1,"pair:sigma");
|
||||
memory->create(V_sigma_max,n+1,n+1,"pair:V_sigma_max");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 6) error->all("Illegal pair_style command");
|
||||
|
||||
cut_global = 0.0;
|
||||
max_cell_size = force->numeric(arg[0]);
|
||||
seed = force->inumeric(arg[1]);
|
||||
weighting = force->numeric(arg[2]);
|
||||
T_ref = force->numeric(arg[3]);
|
||||
recompute_vsigmamax_stride = force->inumeric(arg[4]);
|
||||
vsigmamax_samples = force->inumeric(arg[5]);
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
|
||||
if (max_cell_size <= 0.0) error->all("Illegal pair_style command");
|
||||
if (seed <= 0) error->all("Illegal pair_style command");
|
||||
if (random) delete random;
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
kT_ref = force->boltz*T_ref;
|
||||
|
||||
// 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[i][j] = cut_global;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 3 || narg > 4) error->all("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 sigma_one = force->numeric(arg[2]);
|
||||
|
||||
double cut_one = cut_global;
|
||||
if (narg == 4) cut_one = force->numeric(arg[3]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
sigma[i][j] = sigma_one;
|
||||
cut[i][j] = cut_one;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all("Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::init_style()
|
||||
{
|
||||
ncellsx = ncellsy = ncellsz = 1;
|
||||
while (((domain->boxhi[0] - domain->boxlo[0])/ncellsx) > max_cell_size)
|
||||
ncellsx++;
|
||||
while (((domain->boxhi[1] - domain->boxlo[1])/ncellsy) > max_cell_size)
|
||||
ncellsy++;
|
||||
while (((domain->boxhi[2] - domain->boxlo[2])/ncellsz) > max_cell_size)
|
||||
ncellsz++;
|
||||
|
||||
cellx = (domain->boxhi[0] - domain->boxlo[0])/ncellsx;
|
||||
celly = (domain->boxhi[1] - domain->boxlo[1])/ncellsy;
|
||||
cellz = (domain->boxhi[2] - domain->boxlo[2])/ncellsz;
|
||||
|
||||
if (comm->me == 0) {
|
||||
if (screen) fprintf(screen,"DSMC cell size = %g x %g x %g\n",
|
||||
cellx,celly,cellz);
|
||||
if (logfile) fprintf(logfile,"DSMC cell size = %g x %g x %g\n",
|
||||
cellx,celly,cellz);
|
||||
}
|
||||
|
||||
total_ncells = ncellsx*ncellsy*ncellsz;
|
||||
vol = cellx*celly*cellz;
|
||||
|
||||
memory->create(particle_list,atom->ntypes+1,0,"pair:particle_list");
|
||||
memory->create(first,atom->ntypes+1,total_ncells,"pair:first");
|
||||
memory->create(number,atom->ntypes+1,total_ncells,"pair:number");
|
||||
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
for (int j = 1; j <= atom->ntypes; j++)
|
||||
V_sigma_max[i][j] = 0.0;
|
||||
|
||||
two_pi = 8.0*atan(1.0);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairDSMC::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) cut[i][j] = 0.0;
|
||||
return cut[i][j];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::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(&sigma[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::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(&sigma[i][j],sizeof(double),1,fp);
|
||||
fread(&cut[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_global,sizeof(double),1,fp);
|
||||
fwrite(&max_cell_size,sizeof(double),1,fp);
|
||||
fwrite(&seed,sizeof(int),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairDSMC::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_global,sizeof(double),1,fp);
|
||||
fread(&max_cell_size,sizeof(double),1,fp);
|
||||
fread(&seed,sizeof(int),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&max_cell_size,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);
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
// same seed that pair_style command initially specified
|
||||
|
||||
if (random) delete random;
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
}
|
||||
|
||||
/*-------------------------------------------------------------------------
|
||||
rezero and recompute the V_sigma_max values this timestep for use during
|
||||
the next nrezero timesteps
|
||||
-------------------------------------------------------------------------*/
|
||||
|
||||
void PairDSMC::recompute_V_sigma_max(int icell)
|
||||
{
|
||||
int i,j,k;
|
||||
double Vsigma_max = 0;
|
||||
|
||||
if (number_of_A && number_of_B) {
|
||||
for (k = 0; k < vsigmamax_samples; k++) {
|
||||
i = particle_list[itype]
|
||||
[static_cast<int>(random->uniform()*number_of_A)];
|
||||
j = particle_list[jtype]
|
||||
[static_cast<int>(random->uniform()*number_of_B)];
|
||||
if (i == j) continue;
|
||||
Vsigma_max = MAX(Vsigma_max,V_sigma(i,j));
|
||||
}
|
||||
}
|
||||
V_sigma_max[itype][jtype] = Vsigma_max;
|
||||
}
|
||||
|
||||
/*-------------------------------------------------------------------------
|
||||
VHS model
|
||||
compute the velocity vector difference between i and j and multiply by
|
||||
their combined collision cross section, sigma, for neutral-neutral
|
||||
collisions using the Variable Hard Sphere model
|
||||
-------------------------------------------------------------------------*/
|
||||
|
||||
double PairDSMC::V_sigma(int i, int j)
|
||||
{
|
||||
double relative_velocity_sq,relative_velocity,pair_sigma;
|
||||
double delv[3];
|
||||
double *vi = atom->v[i];
|
||||
double *vj = atom->v[j];
|
||||
|
||||
subtract3d(vi,vj,delv);
|
||||
relative_velocity_sq = dot3d(delv,delv);
|
||||
relative_velocity = sqrt(relative_velocity_sq);
|
||||
|
||||
// from Bird eq 4.63, and omega=0.67
|
||||
// (omega - 0.5) = 0.17
|
||||
// 1/GAMMA(2.5 - omega) = 1.06418029298371
|
||||
|
||||
if (relative_velocity_sq != 0.0)
|
||||
pair_sigma = sigma[itype][jtype]*
|
||||
pow(kT_ref/(0.5*reduced_mass*relative_velocity_sq),0.17) *
|
||||
1.06418029298371;
|
||||
else
|
||||
pair_sigma = 0.0;
|
||||
|
||||
return relative_velocity*pair_sigma;
|
||||
}
|
||||
|
||||
/*-------------------------------------------------------------------------
|
||||
generate new velocities for collided particles
|
||||
-------------------------------------------------------------------------*/
|
||||
|
||||
void PairDSMC::scatter_random(int i, int j, int icell)
|
||||
{
|
||||
double mag_delv,cos_phi,cos_squared,r,theta;
|
||||
double delv[3],vcm[3];
|
||||
double *vi = atom->v[i];
|
||||
double *vj = atom->v[j];
|
||||
|
||||
subtract3d(vi,vj,delv);
|
||||
if (itype == jtype) mag_delv = sqrt(dot3d(delv,delv))*0.5;
|
||||
else mag_delv = sqrt(dot3d(delv,delv));
|
||||
|
||||
cos_phi = 1.0 - (2.0*random->uniform());
|
||||
cos_squared = MIN(1.0,cos_phi*cos_phi);
|
||||
r = sqrt(1.0 - cos_squared);
|
||||
delv[0] = cos_phi*mag_delv;
|
||||
theta = two_pi*random->uniform();
|
||||
delv[1] = r*mag_delv*cos(theta);
|
||||
delv[2] = r*mag_delv*sin(theta);
|
||||
|
||||
if (itype == jtype) {
|
||||
vcm[0] = (vi[0]+vj[0])*0.5;
|
||||
vcm[1] = (vi[1]+vj[1])*0.5;
|
||||
vcm[2] = (vi[2]+vj[2])*0.5;
|
||||
vi[0] = vcm[0] + delv[0];
|
||||
vi[1] = vcm[1] + delv[1];
|
||||
vi[2] = vcm[2] + delv[2];
|
||||
vj[0] = vcm[0] - delv[0];
|
||||
vj[1] = vcm[1] - delv[1];
|
||||
vj[2] = vcm[2] - delv[2];
|
||||
} else {
|
||||
vcm[0] = vi[0]*imass_tmass + vj[0]*jmass_tmass;
|
||||
vcm[1] = vi[1]*imass_tmass + vj[1]*jmass_tmass;
|
||||
vcm[2] = vi[2]*imass_tmass + vj[2]*jmass_tmass;
|
||||
vi[0] = vcm[0] + delv[0]*jmass_tmass;
|
||||
vi[1] = vcm[1] + delv[1]*jmass_tmass;
|
||||
vi[2] = vcm[2] + delv[2]*jmass_tmass;
|
||||
vj[0] = vcm[0] - delv[0]*imass_tmass;
|
||||
vj[1] = vcm[1] - delv[1]*imass_tmass;
|
||||
vj[2] = vcm[2] - delv[2]*imass_tmass;
|
||||
}
|
||||
|
||||
total_number_of_collisions++;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
This method converts the double supplied by the calling function into
|
||||
an int, which is returned. By adding a random number between 0 and 1
|
||||
to the double before converting it to an int, we ensure that,
|
||||
statistically, we round down with probability identical to the
|
||||
remainder and up the rest of the time. So even though we're using an
|
||||
integer, we're statistically matching the exact expression represented
|
||||
by the double.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int PairDSMC::convert_double_to_equivalent_int(double input_double)
|
||||
{
|
||||
if (input_double > INT_MAX)
|
||||
error->all("Tried to convert a double to int, but input_double > INT_MAX");
|
||||
|
||||
int output_int = static_cast<int>(input_double + random->uniform());
|
||||
return output_int;
|
||||
}
|
||||
|
|
@ -0,0 +1,109 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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(dsmc,PairDSMC)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_DSMC_H
|
||||
#define LMP_PAIR_DSMC_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairDSMC : public Pair {
|
||||
public:
|
||||
PairDSMC(class LAMMPS *);
|
||||
virtual ~PairDSMC();
|
||||
virtual void compute(int, int);
|
||||
virtual void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
virtual void write_restart_settings(FILE *);
|
||||
virtual void read_restart_settings(FILE *);
|
||||
|
||||
private:
|
||||
double cut_global;
|
||||
double **cut;
|
||||
double **sigma;
|
||||
|
||||
double cellx;
|
||||
double celly;
|
||||
double cellz;
|
||||
int ncellsx;
|
||||
int ncellsy;
|
||||
int ncellsz;
|
||||
int total_ncells;
|
||||
int total_number_of_collisions;
|
||||
int recompute_vsigmamax_stride;
|
||||
int vsigmamax_samples;
|
||||
double T_ref;
|
||||
double kT_ref;
|
||||
double two_pi;
|
||||
double max_cell_size;
|
||||
|
||||
int seed;
|
||||
int number_of_A;
|
||||
int number_of_B;
|
||||
int max_particle_list;
|
||||
|
||||
class RanMars *random;
|
||||
|
||||
int **particle_list;
|
||||
int **first;
|
||||
int **number;
|
||||
|
||||
double **V_sigma_max;
|
||||
|
||||
int max_particles;
|
||||
int *next_particle;
|
||||
|
||||
int itype;
|
||||
int jtype;
|
||||
|
||||
double imass;
|
||||
double jmass;
|
||||
double total_mass;
|
||||
double reduced_mass;
|
||||
double imass_tmass;
|
||||
double jmass_tmass;
|
||||
double vol;
|
||||
double weighting;
|
||||
|
||||
void allocate();
|
||||
void recompute_V_sigma_max(int);
|
||||
double V_sigma(int, int);
|
||||
void scatter_random(int, int, int);
|
||||
int convert_double_to_equivalent_int(double);
|
||||
|
||||
inline void subtract3d(const double *v1, const double *v2, double *v3) {
|
||||
v3[0] = v2[0] - v1[0];
|
||||
v3[1] = v2[1] - v1[1];
|
||||
v3[2] = v2[2] - v1[2];
|
||||
}
|
||||
|
||||
inline double dot3d(const double *v1, const double *v2) {
|
||||
return( v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] );
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
Loading…
Reference in New Issue