Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80324742
fix_lb_rigid_pc_sphere.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Aug 31, 04:04
Size
53 KB
Mime Type
text/x-c
Expires
Mon, Sep 2, 04:04 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
20292721
Attached To
rLAMMPS lammps
fix_lb_rigid_pc_sphere.cpp
View Options
/* ----------------------------------------------------------------------
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: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
Based on fix_rigid (version from 2008).
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "fix_lb_rigid_pc_sphere.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "update.h"
#include "respa.h"
#include "modify.h"
#include "group.h"
#include "comm.h"
#include "force.h"
#include "output.h"
#include "memory.h"
#include "error.h"
#include "fix_lb_fluid.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* -------------------------------------------------------------------------- */
FixLbRigidPCSphere::FixLbRigidPCSphere(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
int i, ibody;
scalar_flag = 1;
extscalar = 0;
time_integrate = 1;
rigid_flag = 1;
create_attribute = 1;
virial_flag = 1;
// perform initial allocation of atom-based arrays
// register with Atom class
body = NULL;
up = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
// by default assume all of the particles interact with the fluid.
inner_nodes = 0;
// parse command-line args
// set nbody and body[i] for each atom
if (narg < 4) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
// single rigid body
// nbody = 1
// all atoms in fix group are part of body
int iarg;
if (strcmp(arg[3],"single") == 0) {
iarg = 4;
nbody = 1;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
body[i] = -1;
if (mask[i] & groupbit) body[i] = 0;
}
// each molecule in fix group is a rigid body
// maxmol = largest molecule #
// ncount = # of atoms in each molecule (have to sum across procs)
// nbody = # of non-zero ncount values
// use nall as incremented ptr to set body[] values for each atom
} else if (strcmp(arg[3],"molecule") == 0) {
iarg = 4;
if (atom->molecular == 0)
error->all(FLERR,"Must use a molecular atom style with "
"fix lb/rigid/pc/sphere molecule");
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
tagint maxmol_tag = -1;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) maxmol_tag = MAX(maxmol_tag,molecule[i]);
tagint itmp;
MPI_Allreduce(&maxmol_tag,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
if (itmp+1 > MAXSMALLINT)
error->all(FLERR,"Too many molecules for fix lb/rigid/pc/sphere");
int maxmol = (int) itmp;
int *ncount;
memory->create(ncount,maxmol+1,"rigid:ncount");
for (i = 0; i <= maxmol; i++) ncount[i] = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) ncount[molecule[i]]++;
int *nall;
memory->create(nall,maxmol+1,"rigid:ncount");
MPI_Allreduce(ncount,nall,maxmol+1,MPI_LMP_TAGINT,MPI_SUM,world);
nbody = 0;
for (i = 0; i <= maxmol; i++)
if (nall[i]) nall[i] = nbody++;
else nall[i] = -1;
for (i = 0; i < nlocal; i++) {
body[i] = -1;
if (mask[i] & groupbit) body[i] = nall[molecule[i]];
}
memory->destroy(ncount);
memory->destroy(nall);
// each listed group is a rigid body
// check if all listed groups exist
// an atom must belong to fix group and listed group to be in rigid body
// error if atom belongs to more than 1 rigid body
} else if (strcmp(arg[3],"group") == 0) {
if (narg < 5) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
nbody = atoi(arg[4]);
if (nbody <= 0) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
if (narg < 5+nbody)
error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
iarg = 5 + nbody;
int *igroups = new int[nbody];
for (ibody = 0; ibody < nbody; ibody++) {
igroups[ibody] = group->find(arg[5+ibody]);
if (igroups[ibody] == -1)
error->all(FLERR,"Could not find fix lb/rigid/pc/sphere group ID");
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
int flag = 0;
for (i = 0; i < nlocal; i++) {
body[i] = -1;
if (mask[i] & groupbit)
for (ibody = 0; ibody < nbody; ibody++)
if (mask[i] & group->bitmask[igroups[ibody]]) {
if (body[i] >= 0) flag = 1;
body[i] = ibody;
}
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall)
error->all(FLERR,"One or more atoms belong to multiple rigid bodies");
delete [] igroups;
}else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
// error check on nbody
if (nbody == 0) error->all(FLERR,"No rigid bodies defined");
// create all nbody-length arrays
memory->create(nrigid,nbody,"lb/rigid/pc/sphere:nrigid");
memory->create(nrigid_shell,nbody,"lb/rigid/pc/sphere:nrigid_shell");
memory->create(masstotal,nbody,"lb/rigid/pc/sphere:masstotal");
memory->create(masstotal_shell,nbody,"lb/rigid/pc/sphere:masstotal_shell");
memory->create(sphereradius,nbody,"lb/rigid/pc/sphere:sphereradius");
memory->create(xcm,nbody,3,"lb/rigid/pc/sphere:xcm");
memory->create(xcm_old,nbody,3,"lb/rigid/pc/sphere:xcm_old");
memory->create(vcm,nbody,3,"lb/rigid/pc/sphere:vcm");
memory->create(ucm,nbody,3,"lb/rigid/pc/sphere:ucm");
memory->create(ucm_old,nbody,3,"lb/rigid/pc/sphere:ucm_old");
memory->create(fcm,nbody,3,"lb/rigid/pc/sphere:fcm");
memory->create(fcm_old,nbody,3,"lb/rigid/pc/sphere:fcm_old");
memory->create(fcm_fluid,nbody,3,"lb/rigid/pc/sphere:fcm_fluid");
memory->create(omega,nbody,3,"lb/rigid/pc/sphere:omega");
memory->create(torque,nbody,3,"lb/rigid/pc/sphere:torque");
memory->create(torque_old,nbody,3,"lb/rigid/pc/sphere:torque_old");
memory->create(torque_fluid,nbody,3,"lb/rigid/pc/sphere:torque_fluid");
memory->create(torque_fluid_old,nbody,3,"lb/rigid/pc/sphere:torque_fluid_old");
memory->create(rotate,nbody,3,"lb/rigid/pc/sphere:rotate");
memory->create(imagebody,nbody,"lb/rigid/pc/sphere:imagebody");
memory->create(fflag,nbody,3,"lb/rigid/pc/sphere:fflag");
memory->create(tflag,nbody,3,"lb/rigid/pc/sphere:tflag");
memory->create(sum,nbody,6,"lb/rigid/pc/sphere:sum");
memory->create(all,nbody,6,"lb/rigid/pc/sphere:all");
memory->create(remapflag,nbody,4,"lb/rigid/pc/sphere:remapflag");
Gamma_MD = new double[nbody];
// initialize force/torque flags to default = 1.0
array_flag = 1;
size_array_rows = nbody;
size_array_cols = 15;
global_freq = 1;
extarray = 0;
for (i = 0; i < nbody; i++) {
fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0;
tflag[i][0] = tflag[i][1] = tflag[i][2] = 1.0;
}
// parse optional args that set fflag and tflag
while (iarg < narg) {
if (strcmp(arg[iarg],"force") == 0) {
if (iarg+5 > narg) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
int mlo,mhi;
force->bounds(arg[iarg+1],nbody,mlo,mhi);
double xflag,yflag,zflag;
if (strcmp(arg[iarg+2],"off") == 0) xflag = 0.0;
else if (strcmp(arg[iarg+2],"on") == 0) xflag = 1.0;
else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
if (strcmp(arg[iarg+2],"off") == 0) yflag = 0.0;
else if (strcmp(arg[iarg+3],"on") == 0) yflag = 1.0;
else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
if (strcmp(arg[iarg+4],"off") == 0) zflag = 0.0;
else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
int count = 0;
for (int m = mlo; m <= mhi; m++) {
fflag[m-1][0] = xflag;
fflag[m-1][1] = yflag;
fflag[m-1][2] = zflag;
count++;
}
if (count == 0) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
iarg += 5;
} else if (strcmp(arg[iarg],"torque") == 0) {
if (iarg+5 > narg) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
int mlo,mhi;
force->bounds(arg[iarg+1],nbody,mlo,mhi);
double xflag,yflag,zflag;
if (strcmp(arg[iarg+2],"off") == 0) xflag = 0.0;
else if (strcmp(arg[iarg+2],"on") == 0) xflag = 1.0;
else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
if (strcmp(arg[iarg+3],"off") == 0) yflag = 0.0;
else if (strcmp(arg[iarg+3],"on") == 0) yflag = 1.0;
else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
if (strcmp(arg[iarg+4],"off") == 0) zflag = 0.0;
else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
int count = 0;
for (int m = mlo; m <= mhi; m++) {
tflag[m-1][0] = xflag;
tflag[m-1][1] = yflag;
tflag[m-1][2] = zflag;
count++;
}
if (count == 0) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
iarg += 5;
// specify if certain particles are inside the rigid spherical body,
// and therefore should not
} else if(strcmp(arg[iarg],"innerNodes")==0){
inner_nodes = 1;
igroupinner = group->find(arg[iarg+1]);
if(igroupinner == -1)
error->all(FLERR,"Could not find fix lb/rigid/pc/sphere innerNodes group ID");
iarg += 2;
} else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
}
// initialize vector output quantities in case accessed before run
for (i = 0; i < nbody; i++) {
xcm[i][0] = xcm[i][1] = xcm[i][2] = 0.0;
xcm_old[i][0] = xcm_old[i][1] = xcm_old[i][2] = 0.0;
vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0;
ucm[i][0] = ucm[i][1] = ucm[i][2] = 0.0;
ucm_old[i][0] = ucm_old[i][1] = ucm_old[i][2] = 0.0;
fcm[i][0] = fcm[i][1] = fcm[i][2] = 0.0;
fcm_old[i][0] = fcm_old[i][1] = fcm_old[i][2] = 0.0;
fcm_fluid[i][0] = fcm_fluid[i][1] = fcm_fluid[i][2] = 0.0;
torque[i][0] = torque[i][1] = torque[i][2] = 0.0;
torque_old[i][0] = torque_old[i][1] = torque_old[i][2] = 0.0;
torque_fluid[i][0] = torque_fluid[i][1] = torque_fluid[i][2] = 0.0;
torque_fluid_old[i][0] = torque_fluid_old[i][1] = torque_fluid_old[i][2] = 0.0;
}
// nrigid[n] = # of atoms in Nth rigid body
// error if one or zero atoms
int *ncount = new int[nbody];
for (ibody = 0; ibody < nbody; ibody++) ncount[ibody] = 0;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (body[i] >= 0) ncount[body[i]]++;
MPI_Allreduce(ncount,nrigid,nbody,MPI_INT,MPI_SUM,world);
//count the number of atoms in the shell.
if(inner_nodes == 1){
int *mask = atom->mask;
for(ibody=0; ibody<nbody; ibody++) ncount[ibody] = 0;
for(i=0; i<nlocal; i++){
if(!(mask[i] & group->bitmask[igroupinner])){
if(body[i] >= 0) ncount[body[i]]++;
}
}
MPI_Allreduce(ncount,nrigid_shell,nbody,MPI_INT,MPI_SUM,world);
}else {
for(ibody=0; ibody < nbody; ibody++) nrigid_shell[ibody]=nrigid[ibody];
}
delete [] ncount;
for (ibody = 0; ibody < nbody; ibody++)
if (nrigid[ibody] <= 1) error->all(FLERR,"One or zero atoms in rigid body");
// set image flags for each rigid body to default values
// will be reset during init() based on xcm and then by pre_neighbor()
// set here, so image value will persist from run to run
for (ibody = 0; ibody < nbody; ibody++)
imagebody[ibody] = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
// print statistics
int nsum = 0;
for (ibody = 0; ibody < nbody; ibody++) nsum += nrigid[ibody];
if (comm->me == 0) {
if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum);
if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum);
}
int groupbit_lb_fluid = 0;
for(int ifix=0; ifix<modify->nfix; ifix++)
if(strcmp(modify->fix[ifix]->style,"lb/fluid")==0){
fix_lb_fluid = (FixLbFluid *)modify->fix[ifix];
groupbit_lb_fluid = group->bitmask[modify->fix[ifix]->igroup];
}
if(groupbit_lb_fluid == 0)
error->all(FLERR,"the lb/fluid fix must also be used if using the lb/rigid/pc/sphere fix");
int *mask = atom->mask;
if(inner_nodes == 1){
for(int j=0; j<nlocal; j++){
if((mask[j] & groupbit) && !(mask[j] & group->bitmask[igroupinner]) && !(mask[j] & groupbit_lb_fluid))
error->one(FLERR,"use the innerNodes keyword in the lb/rigid/pc/sphere fix for atoms which do not interact with the lb/fluid");
// If inner nodes are present, which should not interact with the fluid, make
// sure these are not used by the lb/fluid fix to apply a force to the fluid.
if((mask[j] & groupbit) && (mask[j] & groupbit_lb_fluid) && (mask[j] & group->bitmask[igroupinner]))
error->one(FLERR,"the inner nodes specified in lb/rigid/pc/sphere should not be included in the lb/fluid fix");
}
}else{
for(int j=0; j<nlocal; j++){
if((mask[j] & groupbit) && !(mask[j] & groupbit_lb_fluid))
error->one(FLERR,"use the innerNodes keyword in the lb/rigid/pc/sphere fix for atoms which do not interact with the lb/fluid");
}
}
}
/* ---------------------------------------------------------------------- */
FixLbRigidPCSphere::~FixLbRigidPCSphere()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
// delete locally stored arrays
memory->destroy(body);
memory->destroy(up);
// delete nbody-length arrays
memory->destroy(nrigid);
memory->destroy(nrigid_shell);
memory->destroy(masstotal);
memory->destroy(masstotal_shell);
memory->destroy(sphereradius);
memory->destroy(xcm);
memory->destroy(xcm_old);
memory->destroy(vcm);
memory->destroy(ucm);
memory->destroy(ucm_old);
memory->destroy(fcm);
memory->destroy(fcm_old);
memory->destroy(fcm_fluid);
memory->destroy(omega);
memory->destroy(torque);
memory->destroy(torque_old);
memory->destroy(torque_fluid);
memory->destroy(torque_fluid_old);
memory->destroy(rotate);
memory->destroy(imagebody);
memory->destroy(fflag);
memory->destroy(tflag);
memory->destroy(sum);
memory->destroy(all);
memory->destroy(remapflag);
delete [] Gamma_MD;
}
/* ---------------------------------------------------------------------- */
int FixLbRigidPCSphere::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= PRE_NEIGHBOR;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixLbRigidPCSphere::init()
{
int i,itype,ibody;
// warn if more than one rigid fix
int count = 0;
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"lb/rigid/pc/sphere") == 0) count++;
if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one fix lb/rigid/pc/sphere");
// timestep info
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
int *type = atom->type;
int nlocal = atom->nlocal;
double **x = atom->x;
imageint *image = atom->image;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *periodicity = domain->periodicity;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double **mu = atom->mu;
double *radius = atom->radius;
int *ellipsoid = atom->ellipsoid;
int extended = 0;
int *mask = atom->mask;
// Warn if any extended particles are included.
if (atom->radius_flag || atom->ellipsoid_flag || atom->mu_flag) {
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
if (radius && radius[i] > 0.0) flag = 1;
if (ellipsoid && ellipsoid[i] >= 0) flag = 1;
if (mu && mu[i][3] > 0.0) flag = 1;
}
MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world);
}
if(extended)
error->warning(FLERR,"Fix lb/rigid/pc/sphere assumes point particles");
// compute masstotal & center-of-mass of each rigid body
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
int xbox,ybox,zbox;
double massone,xunwrap,yunwrap,zunwrap;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
if ((xbox && !periodicity[0]) || (ybox && !periodicity[1]) ||
(zbox && !periodicity[2]))
error->one(FLERR,"Fix lb/rigid/pc/sphere atom has non-zero image flag "
"in a non-periodic dimension");
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
sum[ibody][0] += xunwrap * massone;
sum[ibody][1] += yunwrap * massone;
sum[ibody][2] += zunwrap * massone;
sum[ibody][3] += massone;
if(inner_nodes == 1){
if(!(mask[i] & group->bitmask[igroupinner])){
sum[ibody][4] += massone;
}
}else{
sum[ibody][4] += massone;
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for (ibody = 0; ibody < nbody; ibody++) {
masstotal[ibody] = all[ibody][3];
masstotal_shell[ibody] = all[ibody][4];
xcm[ibody][0] = all[ibody][0]/masstotal[ibody];
xcm[ibody][1] = all[ibody][1]/masstotal[ibody];
xcm[ibody][2] = all[ibody][2]/masstotal[ibody];
}
// Calculate the radius of the rigid body, and assign the value for gamma:
double dx,dy,dz;
double *Gamma = fix_lb_fluid->Gamma;
double dm_lb = fix_lb_fluid->dm_lb;
double dt_lb = fix_lb_fluid->dt_lb;
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
for (i=0; i<nlocal; i++){
if(body[i] < 0) continue;
if(inner_nodes == 1){
if(!(mask[i] & group->bitmask[igroupinner])){
ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
dx = xunwrap - xcm[ibody][0];
dy = yunwrap - xcm[ibody][1];
dz = zunwrap - xcm[ibody][2];
sum[ibody][0] += dx*dx + dy*dy + dz*dz;
sum[ibody][1] += Gamma[type[i]];
}
}else{
ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
dx = xunwrap - xcm[ibody][0];
dy = yunwrap - xcm[ibody][1];
dz = zunwrap - xcm[ibody][2];
sum[ibody][0] += dx*dx + dy*dy + dz*dz;
sum[ibody][1] += Gamma[type[i]];
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for(ibody=0; ibody < nbody; ibody++){
sphereradius[ibody] = sqrt(all[ibody][0]/nrigid_shell[ibody]);
Gamma_MD[ibody] = all[ibody][1]*dm_lb/dt_lb/nrigid_shell[ibody];
}
// Check that all atoms in the rigid body have the same value of gamma.
double eps = 1.0e-7;
for (i=0; i<nlocal; i++){
if(body[i] < 0) continue;
if(inner_nodes == 1){
if(!(mask[i] & group->bitmask[igroupinner])){
ibody = body[i];
if(Gamma_MD[ibody]*dt_lb/dm_lb - Gamma[type[i]] > eps)
error->one(FLERR,"All atoms in a rigid body must have the same gamma value");
}
}else{
ibody = body[i];
if(Gamma_MD[ibody]*dt_lb/dm_lb - Gamma[type[i]] > eps)
error->one(FLERR,"All atoms in a rigid body must have the same gamma value");
}
}
// remap the xcm of each body back into simulation box if needed
// only really necessary the 1st time a run is performed
pre_neighbor();
// temperature scale factor
double ndof = 0.0;
for (ibody = 0; ibody < nbody; ibody++) {
ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2];
ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2];
}
if (ndof > 0.0) tfactor = force->mvv2e / (ndof * force->boltz);
else tfactor = 0.0;
}
/* ---------------------------------------------------------------------- */
void FixLbRigidPCSphere::setup(int vflag)
{
int i,n,ibody;
double massone,radone;
// vcm = velocity of center-of-mass of each rigid body
// fcm = force on center-of-mass of each rigid body
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
imageint *image = atom->image;
int *periodicity = domain->periodicity;
double unwrap[3];
double dx,dy,dz;
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
sum[ibody][0] += v[i][0] * massone;
sum[ibody][1] += v[i][1] * massone;
sum[ibody][2] += v[i][2] * massone;
sum[ibody][3] += f[i][0];
sum[ibody][4] += f[i][1];
sum[ibody][5] += f[i][2];
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for (ibody = 0; ibody < nbody; ibody++) {
vcm[ibody][0] = all[ibody][0]/masstotal[ibody];
vcm[ibody][1] = all[ibody][1]/masstotal[ibody];
vcm[ibody][2] = all[ibody][2]/masstotal[ibody];
fcm[ibody][0] = all[ibody][3];
fcm[ibody][1] = all[ibody][4];
fcm[ibody][2] = all[ibody][5];
}
// omega = angular velocity of each rigid body
// Calculated as the average of the angular velocities of the
// individual atoms comprising the rigid body.
// torque = torque on each rigid body
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
sum[ibody][0] += (dy * (v[i][2]-vcm[ibody][2]) - dz * (v[i][1]-vcm[ibody][1]))/(dx*dx+dy*dy+dz*dz);
sum[ibody][1] += (dz * (v[i][0]-vcm[ibody][0]) - dx * (v[i][2]-vcm[ibody][2]))/(dx*dx+dy*dy+dz*dz);
sum[ibody][2] += (dx * (v[i][1]-vcm[ibody][1]) - dy * (v[i][0]-vcm[ibody][0]))/(dx*dx+dy*dy+dz*dz);
sum[ibody][3] += dy * f[i][2] - dz * f[i][1];
sum[ibody][4] += dz * f[i][0] - dx * f[i][2];
sum[ibody][5] += dx * f[i][1] - dy * f[i][0];
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for (ibody = 0; ibody < nbody; ibody++) {
omega[ibody][0] = all[ibody][0]/nrigid[ibody];
omega[ibody][1] = all[ibody][1]/nrigid[ibody];
omega[ibody][2] = all[ibody][2]/nrigid[ibody];
torque[ibody][0] = all[ibody][3];
torque[ibody][1] = all[ibody][4];
torque[ibody][2] = all[ibody][5];
}
// virial setup before call to set_v
if (vflag) v_setup(vflag);
else evflag = 0;
// Set the velocities
set_v();
if (vflag_global)
for (n = 0; n < 6; n++) virial[n] *= 2.0;
if (vflag_atom) {
for (i = 0; i < nlocal; i++)
for (n = 0; n < 6; n++)
vatom[i][n] *= 2.0;
}
}
/* ---------------------------------------------------------------------- */
void FixLbRigidPCSphere::initial_integrate(int vflag)
{
double dtfm;
int i,n,ibody;
double massone,radone;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
imageint *image = atom->image;
double unwrap[3];
double dx,dy,dz;
int *mask = atom->mask;
// compute the fluid velocity at the initial particle positions
compute_up();
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
// Store the fluid velocity at the center of mass
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
if(inner_nodes == 1){
if(!(mask[i] & group->bitmask[igroupinner])){
sum[ibody][0] += up[i][0]*massone;
sum[ibody][1] += up[i][1]*massone;
sum[ibody][2] += up[i][2]*massone;
}
}else{
sum[ibody][0] += up[i][0]*massone;
sum[ibody][1] += up[i][1]*massone;
sum[ibody][2] += up[i][2]*massone;
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for (ibody = 0; ibody < nbody; ibody++) {
ucm[ibody][0] = all[ibody][0]/masstotal_shell[ibody];
ucm[ibody][1] = all[ibody][1]/masstotal_shell[ibody];
ucm[ibody][2] = all[ibody][2]/masstotal_shell[ibody];
}
//Store the total torque due to the fluid.
for (ibody = 0; ibody < nbody; ibody++)
for(i = 0; i < 6; i++) sum[ibody][i] = 0.0;
for(i = 0; i<nlocal; i++){
if(body[i] < 0) continue;
ibody = body[i];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
if(inner_nodes == 1){
if(!(mask[i] & group->bitmask[igroupinner])){
sum[ibody][0] += Gamma_MD[ibody]*(dy * ((up[i][2]-vcm[ibody][2])) -
dz * ((up[i][1]-vcm[ibody][1])));
sum[ibody][1] += Gamma_MD[ibody]*(dz * ((up[i][0]-vcm[ibody][0])) -
dx * ((up[i][2]-vcm[ibody][2])));
sum[ibody][2] += Gamma_MD[ibody]*(dx * ((up[i][1]-vcm[ibody][1])) -
dy * ((up[i][0]-vcm[ibody][0])));
sum[ibody][3] += -Gamma_MD[ibody]*(v[i][0]-up[i][0]);
sum[ibody][4] += -Gamma_MD[ibody]*(v[i][1]-up[i][1]);
sum[ibody][5] += -Gamma_MD[ibody]*(v[i][2]-up[i][2]);
}
}else{
sum[ibody][0] += Gamma_MD[ibody]*(dy * ((up[i][2]-vcm[ibody][2])) -
dz * ((up[i][1]-vcm[ibody][1])));
sum[ibody][1] += Gamma_MD[ibody]*(dz * ((up[i][0]-vcm[ibody][0])) -
dx * ((up[i][2]-vcm[ibody][2])));
sum[ibody][2] += Gamma_MD[ibody]*(dx * ((up[i][1]-vcm[ibody][1])) -
dy * ((up[i][0]-vcm[ibody][0])));
sum[ibody][3] += -Gamma_MD[ibody]*(v[i][0]-up[i][0]);
sum[ibody][4] += -Gamma_MD[ibody]*(v[i][1]-up[i][1]);
sum[ibody][5] += -Gamma_MD[ibody]*(v[i][2]-up[i][2]);
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for (ibody = 0; ibody < nbody; ibody++) {
torque_fluid[ibody][0] = all[ibody][0];
torque_fluid[ibody][1] = all[ibody][1];
torque_fluid[ibody][2] = all[ibody][2];
fcm_fluid[ibody][0] = all[ibody][3];
fcm_fluid[ibody][1] = all[ibody][4];
fcm_fluid[ibody][2] = all[ibody][5];
}
for (int ibody = 0; ibody < nbody; ibody++) {
fcm_old[ibody][0] = fcm[ibody][0];
fcm_old[ibody][1] = fcm[ibody][1];
fcm_old[ibody][2] = fcm[ibody][2];
torque_old[ibody][0] = torque[ibody][0];
torque_old[ibody][1] = torque[ibody][1];
torque_old[ibody][2] = torque[ibody][2];
torque_fluid_old[ibody][0] = torque_fluid[ibody][0];
torque_fluid_old[ibody][1] = torque_fluid[ibody][1];
torque_fluid_old[ibody][2] = torque_fluid[ibody][2];
ucm_old[ibody][0] = ucm[ibody][0];
ucm_old[ibody][1] = ucm[ibody][1];
ucm_old[ibody][2] = ucm[ibody][2];
xcm_old[ibody][0] = xcm[ibody][0];
xcm_old[ibody][1] = xcm[ibody][1];
xcm_old[ibody][2] = xcm[ibody][2];
// update xcm by full step
dtfm = dtf / masstotal[ibody];
xcm[ibody][0] += dtv * vcm[ibody][0]+(fcm[ibody][0]+fcm_fluid[ibody][0]/force->ftm2v)*fflag[ibody][0]*dtfm*dtv;
xcm[ibody][1] += dtv * vcm[ibody][1]+(fcm[ibody][1]+fcm_fluid[ibody][1]/force->ftm2v)*fflag[ibody][1]*dtfm*dtv;
xcm[ibody][2] += dtv * vcm[ibody][2]+(fcm[ibody][2]+fcm_fluid[ibody][2]/force->ftm2v)*fflag[ibody][2]*dtfm*dtv;
rotate[ibody][0] = omega[ibody][0]*dtv + tflag[ibody][0]*(torque[ibody][0]*force->ftm2v+torque_fluid[ibody][0])*
dtv*dtv*5.0/(4.0*masstotal[ibody]*sphereradius[ibody]*sphereradius[ibody]);
rotate[ibody][1] = omega[ibody][1]*dtv + tflag[ibody][1]*(torque[ibody][1]*force->ftm2v+torque_fluid[ibody][1])*
dtv*dtv*5.0/(4.0*masstotal[ibody]*sphereradius[ibody]*sphereradius[ibody]);
rotate[ibody][2] = omega[ibody][2]*dtv + tflag[ibody][2]*(torque[ibody][2]*force->ftm2v+torque_fluid[ibody][2])*
dtv*dtv*5.0/(4.0*masstotal[ibody]*sphereradius[ibody]*sphereradius[ibody]);
// Approximate vcm
expminusdttimesgamma = exp(-Gamma_MD[ibody]*dtv*nrigid_shell[ibody]/masstotal[ibody]);
force_factor = force->ftm2v/Gamma_MD[ibody]/nrigid_shell[ibody];
if(fflag[ibody][0]==1){
vcm[ibody][0] = expminusdttimesgamma*(vcm[ibody][0] - ucm[ibody][0] - fcm[ibody][0]*force_factor)
+ ucm[ibody][0] + fcm[ibody][0]*force_factor;
}
if(fflag[ibody][1]==1){
vcm[ibody][1] = expminusdttimesgamma*(vcm[ibody][1] - ucm[ibody][1] - fcm[ibody][1]*force_factor) +
ucm[ibody][1] + fcm[ibody][1]*force_factor;
}
if(fflag[ibody][2]==1){
vcm[ibody][2] = expminusdttimesgamma*(vcm[ibody][2] - ucm[ibody][2] - fcm[ibody][2]*force_factor) +
ucm[ibody][2] + fcm[ibody][2]*force_factor;
}
// Approximate angmom
torque_factor = 5.0*Gamma_MD[ibody]*nrigid_shell[ibody]/(3.0*masstotal[ibody]);
expminusdttimesgamma = exp(-dtv*torque_factor);
if(tflag[ibody][0]==1){
omega[ibody][0] = expminusdttimesgamma*(omega[ibody][0] - (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
(force->ftm2v*torque[ibody][0] + torque_fluid[ibody][0])) +
(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
(force->ftm2v*torque[ibody][0] + torque_fluid[ibody][0]);
}
if(tflag[ibody][1]==1){
omega[ibody][1] = expminusdttimesgamma*(omega[ibody][1] - (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
(force->ftm2v*torque[ibody][1] + torque_fluid[ibody][1])) +
(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
(force->ftm2v*torque[ibody][1] + torque_fluid[ibody][1]);
}
if(tflag[ibody][2]==1){
omega[ibody][2] = expminusdttimesgamma*(omega[ibody][2] - (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
(force->ftm2v*torque[ibody][2] + torque_fluid[ibody][2])) +
(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
(force->ftm2v*torque[ibody][2] + torque_fluid[ibody][2]);
}
}
// virial setup before call to set_xv
if (vflag) v_setup(vflag);
else evflag = 0;
set_xv();
}
/* ---------------------------------------------------------------------- */
void FixLbRigidPCSphere::final_integrate()
{
int i,ibody;
double xy,xz,yz;
// sum over atoms to get force and torque on rigid body
double massone;
imageint *image = atom->image;
double **x = atom->x;
double **f = atom->f;
double **v = atom->v;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
double unwrap[3];
double dx,dy,dz;
int *mask = atom->mask;
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
sum[ibody][0] += f[i][0];
sum[ibody][1] += f[i][1];
sum[ibody][2] += f[i][2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2];
sum[ibody][3] += dy*f[i][2] - dz*f[i][1];
sum[ibody][4] += dz*f[i][0] - dx*f[i][2];
sum[ibody][5] += dx*f[i][1] - dy*f[i][0];
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
//Compute the correction to the velocity and angular momentum due to the non-fluid forces:
for (ibody = 0; ibody < nbody; ibody++) {
fcm[ibody][0] = all[ibody][0];
fcm[ibody][1] = all[ibody][1];
fcm[ibody][2] = all[ibody][2];
torque[ibody][0] = all[ibody][3];
torque[ibody][1] = all[ibody][4];
torque[ibody][2] = all[ibody][5];
expminusdttimesgamma = exp(-dtv*Gamma_MD[ibody]*nrigid_shell[ibody]/masstotal[ibody]);
DMDcoeff= (dtv - (masstotal[ibody]/nrigid_shell[ibody])*(1.0-expminusdttimesgamma)/Gamma_MD[ibody]);
vcm[ibody][0] += fflag[ibody][0]*DMDcoeff*force->ftm2v*(fcm[ibody][0]-fcm_old[ibody][0])/Gamma_MD[ibody]/dtv/nrigid_shell[ibody];
vcm[ibody][1] += fflag[ibody][1]*DMDcoeff*force->ftm2v*(fcm[ibody][1]-fcm_old[ibody][1])/Gamma_MD[ibody]/dtv/nrigid_shell[ibody];
vcm[ibody][2] += fflag[ibody][2]*DMDcoeff*force->ftm2v*(fcm[ibody][2]-fcm_old[ibody][2])/Gamma_MD[ibody]/dtv/nrigid_shell[ibody];
torque_factor = 5.0*Gamma_MD[ibody]*nrigid_shell[ibody]/(3.0*masstotal[ibody]);
expminusdttimesgamma = exp(-dtv*torque_factor);
DMDcoeff = (dtv - (1.0-expminusdttimesgamma)/torque_factor);
omega[ibody][0] += tflag[ibody][0]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*DMDcoeff*
force->ftm2v*(torque[ibody][0] - torque_old[ibody][0])/dtv;
omega[ibody][1] += tflag[ibody][1]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*DMDcoeff*
force->ftm2v*(torque[ibody][1] - torque_old[ibody][1])/dtv;
omega[ibody][2] += tflag[ibody][2]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*DMDcoeff*
force->ftm2v*(torque[ibody][2] - torque_old[ibody][2])/dtv;
}
//Next, calculate the correction to the velocity and angular momentum due to the fluid forces:
//Calculate the fluid velocity at the new particle locations.
compute_up();
// store fluid quantities for the total body
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
// Store the fluid velocity at the center of mass, and the total force
// due to the fluid.
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[ibody][0];
dy = unwrap[1] - xcm[ibody][1];
dz = unwrap[2] - xcm[ibody][2];
if(inner_nodes == 1){
if(!(mask[i] & group->bitmask[igroupinner])){
sum[ibody][0] += up[i][0]*massone;
sum[ibody][1] += up[i][1]*massone;
sum[ibody][2] += up[i][2]*massone;
sum[ibody][3] += Gamma_MD[ibody]*(dy * ((up[i][2]-vcm[ibody][2])) -
dz * ((up[i][1]-vcm[ibody][1])));
sum[ibody][4] += Gamma_MD[ibody]*(dz * ((up[i][0]-vcm[ibody][0])) -
dx * ((up[i][2]-vcm[ibody][2])));
sum[ibody][5] += Gamma_MD[ibody]*(dx * ((up[i][1]-vcm[ibody][1])) -
dy * ((up[i][0]-vcm[ibody][0])));
}
}else{
sum[ibody][0] += up[i][0]*massone;
sum[ibody][1] += up[i][1]*massone;
sum[ibody][2] += up[i][2]*massone;
sum[ibody][3] += Gamma_MD[ibody]*(dy * ((up[i][2]-vcm[ibody][2])) -
dz * ((up[i][1]-vcm[ibody][1])));
sum[ibody][4] += Gamma_MD[ibody]*(dz * ((up[i][0]-vcm[ibody][0])) -
dx * ((up[i][2]-vcm[ibody][2])));
sum[ibody][5] += Gamma_MD[ibody]*(dx * ((up[i][1]-vcm[ibody][1])) -
dy * ((up[i][0]-vcm[ibody][0])));
}
}
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
for (ibody = 0; ibody < nbody; ibody++) {
ucm[ibody][0] = all[ibody][0]/masstotal_shell[ibody];
ucm[ibody][1] = all[ibody][1]/masstotal_shell[ibody];
ucm[ibody][2] = all[ibody][2]/masstotal_shell[ibody];
torque_fluid[ibody][0] = all[ibody][3];
torque_fluid[ibody][1] = all[ibody][4];
torque_fluid[ibody][2] = all[ibody][5];
}
for (ibody = 0; ibody < nbody; ibody++) {
expminusdttimesgamma = exp(-dtv*Gamma_MD[ibody]*nrigid_shell[ibody]/masstotal[ibody]);
DMDcoeff= (dtv - (masstotal[ibody]/nrigid_shell[ibody])*(1.0-expminusdttimesgamma)/Gamma_MD[ibody]);
vcm[ibody][0] += DMDcoeff*fflag[ibody][0]*(ucm[ibody][0]-ucm_old[ibody][0])/dtv;
vcm[ibody][1] += DMDcoeff*fflag[ibody][1]*(ucm[ibody][1]-ucm_old[ibody][1])/dtv;
vcm[ibody][2] += DMDcoeff*fflag[ibody][2]*(ucm[ibody][2]-ucm_old[ibody][2])/dtv;
torque_factor = 5.0*Gamma_MD[ibody]*nrigid_shell[ibody]/(3.0*masstotal[ibody]);
expminusdttimesgamma = exp(-dtv*torque_factor);
DMDcoeff = (dtv - (1.0-expminusdttimesgamma)/torque_factor);
omega[ibody][0] += tflag[ibody][0]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
DMDcoeff*(torque_fluid[ibody][0] - torque_fluid_old[ibody][0])/dtv;
omega[ibody][1] += tflag[ibody][1]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
DMDcoeff*(torque_fluid[ibody][1] - torque_fluid_old[ibody][1])/dtv;
omega[ibody][2] += tflag[ibody][2]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*
DMDcoeff*(torque_fluid[ibody][2] - torque_fluid_old[ibody][2])/dtv;
}
set_v();
}
/* ----------------------------------------------------------------------
set space-frame velocity of each atom in a rigid body
set omega and angmom of extended particles
v = Vcm + (W cross (x - Xcm))
------------------------------------------------------------------------- */
void FixLbRigidPCSphere::set_v()
{
int ibody;
int xbox,ybox,zbox;
double dx,dy,dz;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double vr[6];
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
imageint *image = atom->image;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double xunwrap,yunwrap,zunwrap;
// set v of each atom
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
dx = xunwrap - xcm[ibody][0];
dy = yunwrap - xcm[ibody][1];
dz = zunwrap - xcm[ibody][2];
// save old velocities for virial.
if(evflag){
v0 = v[i][0];
v1 = v[i][1];
v2 = v[i][2];
}
v[i][0] = (omega[ibody][1]*dz - omega[ibody][2]*dy) + vcm[ibody][0];
v[i][1] = (omega[ibody][2]*dx - omega[ibody][0]*dz) + vcm[ibody][1];
v[i][2] = (omega[ibody][0]*dy - omega[ibody][1]*dx) + vcm[ibody][2];
// virial = unwrapped coords dotted into body constraint force
// body constraint force = implied force due to v change minus f external
// assume f does not include forces internal to body
// 1/2 factor b/c initial_integrate contributes other half
// assume per-atom contribution is due to constraint force on that atom
if (evflag) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
fc0 = massone*(v[i][0] - v0)/dtf - f[i][0] + Gamma_MD[ibody]*(v0-up[i][0]);
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1] + Gamma_MD[ibody]*(v1-up[i][1]);
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2] + Gamma_MD[ibody]*(v2-up[i][2]);
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
x0 = x[i][0] + xbox*xprd;
x1 = x[i][1] + ybox*yprd;
x2 = x[i][2] + zbox*zprd;
vr[0] = 0.5*x0*fc0;
vr[1] = 0.5*x1*fc1;
vr[2] = 0.5*x2*fc2;
vr[3] = 0.5*x0*fc1;
vr[4] = 0.5*x0*fc2;
vr[5] = 0.5*x1*fc2;
v_tally(1,&i,1.0,vr);
}
}
}
/* ----------------------------------------------------------------------
set space-frame coords and velocity of each atom in each rigid body
set orientation and rotation of extended particles
x = Q displace + Xcm, mapped back to periodic box
v = Vcm + (W cross (x - Xcm))
------------------------------------------------------------------------- */
void FixLbRigidPCSphere::set_xv()
{
int ibody;
int xbox,ybox,zbox;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double vr[6];
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
imageint *image = atom->image;
int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double xunwrap,yunwrap,zunwrap;
double dx,dy,dz;
// set x and v of each atom
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xunwrap = x[i][0] + xbox*xprd;
yunwrap = x[i][1] + ybox*yprd;
zunwrap = x[i][2] + zbox*zprd;
dx = xunwrap - xcm_old[ibody][0];
dy = yunwrap - xcm_old[ibody][1];
dz = zunwrap - xcm_old[ibody][2];
// save old positions and velocities for virial
if(evflag){
x0 = xunwrap;
x1 = yunwrap;
x2 = zunwrap;
v0 = v[i][0];
v1 = v[i][1];
v2 = v[i][2];
}
// x = displacement from center-of-mass, based on body orientation
// v = vcm + omega around center-of-mass
x[i][0] = dx;
x[i][1] = dy;
x[i][2] = dz;
//Perform the rotations:
dx = x[i][0];
dy = x[i][1];
dz = x[i][2];
x[i][0] = cos(rotate[ibody][2])*dx - sin(rotate[ibody][2])*dy;
x[i][1] = sin(rotate[ibody][2])*dx + cos(rotate[ibody][2])*dy;
dx = x[i][0];
dy = x[i][1];
dz = x[i][2];
x[i][0] = cos(rotate[ibody][1])*dx + sin(rotate[ibody][1])*dz;
x[i][2] = -sin(rotate[ibody][1])*dx + cos(rotate[ibody][1])*dz;
dx = x[i][0];
dy = x[i][1];
dz = x[i][2];
x[i][1] = cos(rotate[ibody][0])*dy - sin(rotate[ibody][0])*dz;
x[i][2] = sin(rotate[ibody][0])*dy + cos(rotate[ibody][0])*dz;
v[i][0] = (omega[ibody][1]*x[i][2] - omega[ibody][2]*x[i][1]) + vcm[ibody][0];
v[i][1] = (omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2]) + vcm[ibody][1];
v[i][2] = (omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0]) + vcm[ibody][2];
// add center of mass to displacement
// map back into periodic box via xbox,ybox,zbox
// for triclinic, add in box tilt factors as well
x[i][0] += xcm[ibody][0] - xbox*xprd;
x[i][1] += xcm[ibody][1] - ybox*yprd;
x[i][2] += xcm[ibody][2] - zbox*zprd;
// virial = unwrapped coords dotted into body constraint force
// body constraint force = implied force due to v change minus f external
// assume f does not include forces internal to body
// 1/2 factor b/c final_integrate contributes other half
// assume per-atom contribution is due to constraint force on that atom
if (evflag) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
fc0 = massone*(v[i][0] - v0)/dtf - f[i][0] + Gamma_MD[ibody]*(v0-up[i][0]);
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1] + Gamma_MD[ibody]*(v1-up[i][1]);
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2] + Gamma_MD[ibody]*(v2-up[i][2]);
vr[0] = 0.5*x0*fc0;
vr[1] = 0.5*x1*fc1;
vr[2] = 0.5*x2*fc2;
vr[3] = 0.5*x0*fc1;
vr[4] = 0.5*x0*fc2;
vr[5] = 0.5*x1*fc2;
v_tally(1,&i,1.0,vr);
}
}
}
/* ----------------------------------------------------------------------
remap xcm of each rigid body back into periodic simulation box
done during pre_neighbor so will be after call to pbc()
and after fix_deform::pre_exchange() may have flipped box
use domain->remap() in case xcm is far away from box
due to 1st definition of rigid body or due to box flip
if don't do this, then atoms of a body which drifts far away
from a triclinic box will be remapped back into box
with huge displacements when the box tilt changes via set_x()
------------------------------------------------------------------------- */
void FixLbRigidPCSphere::pre_neighbor()
{
imageint original,oldimage,newimage;
for (int ibody = 0; ibody < nbody; ibody++) {
original = imagebody[ibody];
domain->remap(xcm[ibody],imagebody[ibody]);
if (original == imagebody[ibody]) remapflag[ibody][3] = 0;
else {
oldimage = original & IMGMASK;
newimage = imagebody[ibody] & IMGMASK;
remapflag[ibody][0] = newimage - oldimage;
oldimage = (original >> IMGBITS) & IMGMASK;
newimage = (imagebody[ibody] >> IMGBITS) & IMGMASK;
remapflag[ibody][1] = newimage - oldimage;
oldimage = original >> IMG2BITS;
newimage = imagebody[ibody] >> IMG2BITS;
remapflag[ibody][2] = newimage - oldimage;
remapflag[ibody][3] = 1;
}
}
// adjust image flags of any atom in a rigid body whose xcm was remapped
imageint *atomimage = atom->image;
int nlocal = atom->nlocal;
int ibody;
imageint idim,otherdims;
for (int i = 0; i < nlocal; i++) {
if (body[i] == -1) continue;
if (remapflag[body[i]][3] == 0) continue;
ibody = body[i];
if (remapflag[ibody][0]) {
idim = atomimage[i] & IMGMASK;
otherdims = atomimage[i] ^ idim;
idim -= remapflag[ibody][0];
idim &= IMGMASK;
atomimage[i] = otherdims | idim;
}
if (remapflag[ibody][1]) {
idim = (atomimage[i] >> IMGBITS) & IMGMASK;
otherdims = atomimage[i] ^ (idim << IMGBITS);
idim -= remapflag[ibody][1];
idim &= IMGMASK;
atomimage[i] = otherdims | (idim << IMGBITS);
}
if (remapflag[ibody][2]) {
idim = atomimage[i] >> IMG2BITS;
otherdims = atomimage[i] ^ (idim << IMG2BITS);
idim -= remapflag[ibody][2];
idim &= IMGMASK;
atomimage[i] = otherdims | (idim << IMG2BITS);
}
}
}
/* ----------------------------------------------------------------------
count # of degrees-of-freedom removed by fix_rigid for atoms in igroup
------------------------------------------------------------------------- */
int FixLbRigidPCSphere::dof(int igroup)
{
int groupbit = group->bitmask[igroup];
// ncount = # of atoms in each rigid body that are also in group
int *mask = atom->mask;
int nlocal = atom->nlocal;
int *ncount = new int[nbody];
for (int ibody = 0; ibody < nbody; ibody++) ncount[ibody] = 0;
for (int i = 0; i < nlocal; i++)
if (body[i] >= 0 && mask[i] & groupbit) ncount[body[i]]++;
int *nall = new int[nbody];
MPI_Allreduce(ncount,nall,nbody,MPI_INT,MPI_SUM,world);
// remove 3N - 6 dof for each rigid body if more than 2 atoms in igroup
// remove 3N - 5 dof for each diatomic rigid body in igroup
int n = 0;
for (int ibody = 0; ibody < nbody; ibody++) {
if (nall[ibody] > 2) n += 3*nall[ibody] - 6;
else if (nall[ibody] == 2) n++;
}
delete [] ncount;
delete [] nall;
return n;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixLbRigidPCSphere::memory_usage()
{
int nmax = atom->nmax;
double bytes = nmax * sizeof(int);
bytes += nmax*3 * sizeof(double);
bytes += maxvatom*6 * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate local atom-based arrays
------------------------------------------------------------------------- */
void FixLbRigidPCSphere::grow_arrays(int nmax)
{
memory->grow(body,nmax,"rigid:body");
memory->grow(up,nmax,3,"rigid:up");
}
/* ----------------------------------------------------------------------
copy values within local atom-based arrays
------------------------------------------------------------------------- */
void FixLbRigidPCSphere::copy_arrays(int i, int j, int delflag)
{
body[j] = body[i];
up[j][0] = up[i][0];
up[j][1] = up[i][1];
up[j][2] = up[i][2];
}
/* ----------------------------------------------------------------------
initialize one atom's array values, called when atom is created
------------------------------------------------------------------------- */
void FixLbRigidPCSphere::set_arrays(int i)
{
body[i] = -1;
up[i][0] = 0.0;
up[i][1] = 0.0;
up[i][2] = 0.0;
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for exchange with another proc
------------------------------------------------------------------------- */
int FixLbRigidPCSphere::pack_exchange(int i, double *buf)
{
buf[0] = body[i];
buf[1] = up[i][0];
buf[2] = up[i][1];
buf[3] = up[i][2];
return 4;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based arrays from exchange with another proc
------------------------------------------------------------------------- */
int FixLbRigidPCSphere::unpack_exchange(int nlocal, double *buf)
{
body[nlocal] = static_cast<int> (buf[0]);
up[nlocal][0] = buf[1];
up[nlocal][1] = buf[2];
up[nlocal][2] = buf[3];
return 4;
}
/* ---------------------------------------------------------------------- */
void FixLbRigidPCSphere::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
}
/* ----------------------------------------------------------------------
return temperature of collection of rigid bodies
non-active DOF are removed by fflag/tflag and in tfactor
------------------------------------------------------------------------- */
double FixLbRigidPCSphere::compute_scalar()
{
double inertia;
double t = 0.0;
for (int i = 0; i < nbody; i++) {
t += masstotal[i] * (fflag[i][0]*vcm[i][0]*vcm[i][0] +
fflag[i][1]*vcm[i][1]*vcm[i][1] + \
fflag[i][2]*vcm[i][2]*vcm[i][2]);
// wbody = angular velocity in body frame
inertia = 2.0*masstotal[i]*sphereradius[i]*sphereradius[i]/5.0;
t += tflag[i][0]*inertia*omega[i][0]*omega[i][0] +
tflag[i][1]*inertia*omega[i][1]*omega[i][1] +
tflag[i][2]*inertia*omega[i][2]*omega[i][2];
}
t *= tfactor;
return t;
}
/* ----------------------------------------------------------------------
return attributes of a rigid body
15 values per body
xcm = 0,1,2; vcm = 3,4,5; fcm = 6,7,8; torque = 9,10,11; image = 12,13,14
------------------------------------------------------------------------- */
double FixLbRigidPCSphere::compute_array(int i, int j)
{
if (j < 3) return xcm[i][j];
if (j < 6) return vcm[i][j-3];
if (j < 9) return (fcm[i][j-6]+fcm_fluid[i][j-6]);
if (j < 12) return (torque[i][j-9]+torque_fluid[i][j-9]);
if (j == 12) return (imagebody[i] & IMGMASK) - IMGMAX;
if (j == 13) return (imagebody[i] >> IMGBITS & IMGMASK) - IMGMAX;
return (imagebody[i] >> IMG2BITS) - IMGMAX;
}
/* ---------------------------------------------------------------------- */
void FixLbRigidPCSphere::compute_up(void)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
double **x = atom->x;
int i,j,k,p,m;
int ix,iy,iz;
int ixp,iyp,izp;
double dx1,dy1,dz1;
int isten,ii,jj,kk;
double r,rsq,weightx,weighty,weightz;
double ****u_lb = fix_lb_fluid->u_lb;
int subNbx = fix_lb_fluid->subNbx;
int subNby = fix_lb_fluid->subNby;
int subNbz = fix_lb_fluid->subNbz;
double dx_lb = fix_lb_fluid->dx_lb;
double dt_lb = fix_lb_fluid->dt_lb;
int trilinear_stencil = fix_lb_fluid->trilinear_stencil;
double FfP[64];
for(i=0; i<nlocal; i++){
if(mask[i] & groupbit){
//Calculate nearest leftmost grid point.
//Since array indices from 1 to subNb-2 correspond to the
// local subprocessor domain (not indices from 0), use the
// ceiling value.
ix = (int)ceil((x[i][0]-domain->sublo[0])/dx_lb);
iy = (int)ceil((x[i][1]-domain->sublo[1])/dx_lb);
iz = (int)ceil((x[i][2]-domain->sublo[2])/dx_lb);
//Calculate distances to the nearest points.
dx1 = x[i][0] - (domain->sublo[0] + (ix-1)*dx_lb);
dy1 = x[i][1] - (domain->sublo[1] + (iy-1)*dx_lb);
dz1 = x[i][2] - (domain->sublo[2] + (iz-1)*dx_lb);
// Need to convert these to lattice units:
dx1 = dx1/dx_lb;
dy1 = dy1/dx_lb;
dz1 = dz1/dx_lb;
up[i][0]=0.0; up[i][1]=0.0; up[i][2]=0.0;
if(trilinear_stencil==0){
isten=0;
for(ii=-1; ii<3; ii++){
rsq=(-dx1+ii)*(-dx1+ii);
if(rsq>=4)
weightx=0.0;
else{
r=sqrt(rsq);
if(rsq>1){
weightx=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
} else{
weightx=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
}
}
for(jj=-1; jj<3; jj++){
rsq=(-dy1+jj)*(-dy1+jj);
if(rsq>=4)
weighty=0.0;
else{
r=sqrt(rsq);
if(rsq>1){
weighty=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
} else{
weighty=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
}
}
for(kk=-1; kk<3; kk++){
rsq=(-dz1+kk)*(-dz1+kk);
if(rsq>=4)
weightz=0.0;
else{
r=sqrt(rsq);
if(rsq>1){
weightz=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
} else{
weightz=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
}
}
ixp = ix+ii;
iyp = iy+jj;
izp = iz+kk;
if(ixp==-1) ixp=subNbx+2;
if(iyp==-1) iyp=subNby+2;
if(izp==-1) izp=subNbz+2;
FfP[isten] = weightx*weighty*weightz;
// interpolated velocity based on delta function.
for(k=0; k<3; k++){
up[i][k] += u_lb[ixp][iyp][izp][k]*FfP[isten];
}
}
}
}
}else{
FfP[0] = (1.-dx1)*(1.-dy1)*(1.-dz1);
FfP[1] = (1.-dx1)*(1.-dy1)*dz1;
FfP[2] = (1.-dx1)*dy1*(1.-dz1);
FfP[3] = (1.-dx1)*dy1*dz1;
FfP[4] = dx1*(1.-dy1)*(1.-dz1);
FfP[5] = dx1*(1.-dy1)*dz1;
FfP[6] = dx1*dy1*(1.-dz1);
FfP[7] = dx1*dy1*dz1;
ixp = (ix+1);
iyp = (iy+1);
izp = (iz+1);
for (k=0; k<3; k++) { // tri-linearly interpolated velocity at node
up[i][k] = u_lb[ix][iy][iz][k]*FfP[0]
+ u_lb[ix][iy][izp][k]*FfP[1]
+ u_lb[ix][iyp][iz][k]*FfP[2]
+ u_lb[ix][iyp][izp][k]*FfP[3]
+ u_lb[ixp][iy][iz][k]*FfP[4]
+ u_lb[ixp][iy][izp][k]*FfP[5]
+ u_lb[ixp][iyp][iz][k]*FfP[6]
+ u_lb[ixp][iyp][izp][k]*FfP[7];
}
}
for(k=0; k<3; k++)
up[i][k] = up[i][k]*dx_lb/dt_lb;
}
}
}
Event Timeline
Log In to Comment