diff --git a/src/displace_atoms.cpp b/src/displace_atoms.cpp index 77e03a0c5..50ffcd86f 100644 --- a/src/displace_atoms.cpp +++ b/src/displace_atoms.cpp @@ -1,242 +1,242 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "stdlib.h" #include "string.h" #include "displace_atoms.h" #include "atom.h" #include "domain.h" #include "lattice.h" #include "comm.h" #include "group.h" #include "random_park.h" #include "error.h" using namespace LAMMPS_NS; enum{MOVE,RAMP,RANDOM}; #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) /* ---------------------------------------------------------------------- */ DisplaceAtoms::DisplaceAtoms(LAMMPS *lmp) : Pointers(lmp) {} /* ---------------------------------------------------------------------- */ void DisplaceAtoms::command(int narg, char **arg) { int i; if (domain->box_exist == 0) error->all("Displace_atoms command before simulation box is defined"); if (narg < 2) error->all("Illegal displace_atoms command"); // init entire system since comm->exchange is done // comm::init needs neighbor::init needs pair::init needs kspace::init, etc if (comm->me == 0 && screen) fprintf(screen,"System init for displace_atoms ...\n"); lmp->init(); if (comm->me == 0 && screen) fprintf(screen,"Displacing atoms ...\n"); // group and style int igroup = group->find(arg[0]); if (igroup == -1) error->all("Could not find displace_atoms group ID"); int groupbit = group->bitmask[igroup]; int style; if (strcmp(arg[1],"move") == 0) style = MOVE; else if (strcmp(arg[1],"ramp") == 0) style = RAMP; else if (strcmp(arg[1],"random") == 0) style = RANDOM; else error->all("Illegal displace_atoms command"); // set option defaults scaleflag = 1; // read options from end of input line if (style == MOVE) options(narg-5,&arg[5]); else if (style == RAMP) options(narg-8,&arg[8]); else if (style == RANDOM) options(narg-6,&arg[6]); // setup scaling if (scaleflag && domain->lattice == NULL) error->all("Use of displace_atoms with undefined lattice"); double xscale,yscale,zscale; if (scaleflag) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; zscale = domain->lattice->zlattice; } else xscale = yscale = zscale = 1.0; // move atoms by 3-vector if (style == MOVE) { double delx = xscale*atof(arg[2]); double dely = yscale*atof(arg[3]); double delz = zscale*atof(arg[4]); double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { x[i][0] += delx; x[i][1] += dely; x[i][2] += delz; } } } // move atoms in ramped fashion if (style == RAMP) { int d_dim; if (strcmp(arg[2],"x") == 0) d_dim = 0; else if (strcmp(arg[2],"y") == 0) d_dim = 1; else if (strcmp(arg[2],"z") == 0) d_dim = 2; else error->all("Illegal displace_atoms ramp command"); double d_lo,d_hi; if (d_dim == 0) { d_lo = xscale*atof(arg[3]); d_hi = xscale*atof(arg[4]); } else if (d_dim == 1) { d_lo = yscale*atof(arg[3]); d_hi = yscale*atof(arg[4]); } else if (d_dim == 2) { d_lo = zscale*atof(arg[3]); d_hi = zscale*atof(arg[4]); } int coord_dim; if (strcmp(arg[5],"x") == 0) coord_dim = 0; else if (strcmp(arg[5],"y") == 0) coord_dim = 1; else if (strcmp(arg[5],"z") == 0) coord_dim = 2; else error->all("Illegal velocity ramp command"); double coord_lo,coord_hi; if (coord_dim == 0) { coord_lo = xscale*atof(arg[6]); coord_hi = xscale*atof(arg[7]); } else if (coord_dim == 1) { coord_lo = yscale*atof(arg[6]); coord_hi = yscale*atof(arg[7]); } else if (coord_dim == 2) { coord_lo = zscale*atof(arg[6]); coord_hi = zscale*atof(arg[7]); } double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; double fraction,dramp; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { fraction = (x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); dramp = d_lo + fraction*(d_hi - d_lo); x[i][d_dim] += dramp; } } } // move atoms randomly if (style == RANDOM) { RanPark *random = new RanPark(lmp,1); double dx = xscale*atof(arg[2]); double dy = yscale*atof(arg[3]); double dz = zscale*atof(arg[4]); int seed = atoi(arg[5]); double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; int j; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { random->reset(seed,x[i]); x[i][0] += dx * 2.0*(random->uniform()-0.5); x[i][1] += dy * 2.0*(random->uniform()-0.5); x[i][2] += dz * 2.0*(random->uniform()-0.5); } } delete random; } // move atoms back inside simulation box and to new processors // use remap() instead of pbc() in case atoms moved a long distance // use irregular() instead of exchange() in case atoms moved a long distance double **x = atom->x; int *image = atom->image; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->reset_box(); comm->setup(); comm->irregular(); if (domain->triclinic) domain->lamda2x(atom->nlocal); // check if any atoms were lost double natoms; double rlocal = atom->nlocal; MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world); if (natoms != atom->natoms) { char str[128]; - sprintf(str,"Lost atoms via displacement: original %.15g current %.15g", + sprintf(str,"Lost atoms via displace_atoms: original %.15g current %.15g", atom->natoms,natoms); error->all(str); } } /* ---------------------------------------------------------------------- parse optional parameters at end of displace_atoms input line ------------------------------------------------------------------------- */ void DisplaceAtoms::options(int narg, char **arg) { if (narg < 0) error->all("Illegal displace_atoms command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all("Illegal displace_atoms command"); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; else error->all("Illegal displace_atoms command"); iarg += 2; } else error->all("Illegal displace_atoms command"); } } diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index f78243fe0..26ec6ca22 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -1,706 +1,706 @@ /* ---------------------------------------------------------------------- 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: Pieter in 't Veld (SNL) ------------------------------------------------------------------------- */ #include "string.h" #include "stdlib.h" #include "math.h" #include "fix_deform.h" #include "atom.h" #include "update.h" #include "comm.h" #include "domain.h" #include "lattice.h" #include "force.h" #include "modify.h" #include "kspace.h" #include "error.h" using namespace LAMMPS_NS; enum{NONE,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME}; enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE}; // same as domain.cpp, fix_nvt_sllod.cpp, compute_temp_deform.cpp enum{NO_REMAP,X_REMAP,V_REMAP}; /* ---------------------------------------------------------------------- */ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg < 4) error->all("Illegal fix deform command"); box_change = 1; nevery = atoi(arg[3]); if (nevery <= 0) error->all("Illegal fix deform command"); // set defaults set = new Set[6]; set[0].style = set[1].style = set[2].style = set[3].style = set[4].style = set[5].style = NONE; // parse arguments triclinic = domain->triclinic; int index; int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"x") == 0 || strcmp(arg[iarg],"y") == 0 || strcmp(arg[iarg],"z") == 0) { if (strcmp(arg[iarg],"x") == 0) index = 0; else if (strcmp(arg[iarg],"y") == 0) index = 1; else if (strcmp(arg[iarg],"z") == 0) index = 2; if (iarg+2 > narg) error->all("Illegal fix deform command"); if (strcmp(arg[iarg+1],"final") == 0) { if (iarg+4 > narg) error->all("Illegal fix deform command"); set[index].style = FINAL; set[index].flo = atof(arg[iarg+2]); set[index].fhi = atof(arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg+1],"delta") == 0) { if (iarg+4 > narg) error->all("Illegal fix deform command"); set[index].style = DELTA; set[index].dlo = atof(arg[iarg+2]); set[index].dhi = atof(arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg+1],"scale") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = SCALE; set[index].scale = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"vel") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = VEL; set[index].vel = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"erate") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = ERATE; set[index].rate = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"trate") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = TRATE; set[index].rate = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"volume") == 0) { set[index].style = VOLUME; iarg += 2; } else error->all("Illegal fix deform command"); } else if (strcmp(arg[iarg],"xy") == 0 || strcmp(arg[iarg],"xz") == 0 || strcmp(arg[iarg],"yz") == 0) { if (triclinic == 0) error->all("Fix deform tilt factors require triclinic box"); if (strcmp(arg[iarg],"xy") == 0) index = 5; else if (strcmp(arg[iarg],"xz") == 0) index = 4; else if (strcmp(arg[iarg],"yz") == 0) index = 3; if (iarg+2 > narg) error->all("Illegal fix deform command"); if (strcmp(arg[iarg+1],"final") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = FINAL; set[index].ftilt = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"delta") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = DELTA; set[index].dtilt = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"vel") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = VEL; set[index].vel = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"erate") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = ERATE; set[index].rate = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"trate") == 0) { if (iarg+3 > narg) error->all("Illegal fix deform command"); set[index].style = TRATE; set[index].rate = atof(arg[iarg+2]); iarg += 3; } else error->all("Illegal fix deform command"); } else break; } // read options from end of input line options(narg-iarg,&arg[iarg]); // check periodicity if ((set[0].style && domain->xperiodic == 0) || (set[1].style && domain->yperiodic == 0) || (set[2].style && domain->zperiodic == 0)) error->all("Cannot fix deform on a non-periodic boundary"); if (set[3].style && (domain->yperiodic == 0 || domain->zperiodic == 0)) error->all("Cannot fix deform on a non-periodic boundary"); if (set[4].style && (domain->xperiodic == 0 || domain->zperiodic == 0)) error->all("Cannot fix deform on a non-periodic boundary"); if (set[5].style && (domain->xperiodic == 0 || domain->yperiodic == 0)) error->all("Cannot fix deform on a non-periodic boundary"); // apply scaling to FINAL,DELTA,VEL since they have distance/vel units int flag = 0; for (int i = 0; i < 6; i++) if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == VEL) flag = 1; if (flag && scaleflag && domain->lattice == NULL) error->all("Use of fix deform with undefined lattice"); if (flag && scaleflag) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; zscale = domain->lattice->zlattice; } else xscale = yscale = zscale = 1.0; // for 3,4,5 scaling is in 1st dimension, e.g. x for xz double map[6]; map[0] = xscale; map[1] = yscale; map[2] = zscale; map[3] = yscale; map[4] = xscale; map[5] = xscale; for (int i = 0; i < 3; i++) { if (set[i].style == FINAL) { set[i].flo *= map[i]; set[i].fhi *= map[i]; } else if (set[i].style == DELTA) { set[i].dlo *= map[i]; set[i].dhi *= map[i]; } else if (set[i].style == VEL) { set[i].vel *= map[i]; } } for (int i = 3; i < 6; i++) { if (set[i].style == FINAL) set[i].ftilt *= map[i]; else if (set[i].style == DELTA) set[i].dtilt *= map[i]; else if (set[i].style == VEL) set[i].vel *= map[i]; } // set initial/final values for box size and shape for FINAL,DELTA,SCALE // final not possible for VEL,ERATE,TRATE since don't know length of run yet // final = initial if no setting for (int i = 0; i < 3; i++) { set[i].lo_target = set[i].lo_stop = set[i].lo_start = domain->boxlo[i]; set[i].hi_target = set[i].hi_stop = set[i].hi_start = domain->boxhi[i]; if (set[i].style == FINAL) { set[i].lo_stop = set[i].flo; set[i].hi_stop = set[i].fhi; } else if (set[i].style == DELTA) { set[i].lo_stop = set[i].lo_start + set[i].dlo; set[i].hi_stop = set[i].hi_start + set[i].dhi; } else if (set[i].style == SCALE) { set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - 0.5*set[i].scale*(set[i].hi_start-set[i].lo_start); set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + 0.5*set[i].scale*(set[i].hi_start-set[i].lo_start); } } for (int i = 3; i < 6; i++) { if (i == 5) set[i].tilt_start = domain->xy; else if (i == 4) set[i].tilt_start = domain->xz; else if (i == 3) set[i].tilt_start = domain->yz; set[i].tilt_flip = set[i].tilt_target = set[i].tilt_stop = set[i].tilt_start; if (set[i].style == FINAL) { set[i].tilt_stop = set[i].ftilt; } else if (set[i].style == DELTA) { set[i].tilt_stop = set[i].tilt_start + set[i].dtilt; } } // for VOLUME, setup links to other dims // fixed, dynamic1,2, vol_start for (int i = 0; i < 3; i++) { set[i].vol_start = domain->xprd * domain->yprd * domain->zprd; if (set[i].style != VOLUME) continue; int other1 = (i+1) % 3; int other2 = (i+2) % 3; if (set[other1].style == NONE) { if (set[other2].style == NONE || set[other2].style == VOLUME) error->all("Fix deform volume setting is invalid"); set[i].substyle = ONE_FROM_ONE; set[i].fixed = other1; set[i].dynamic1 = other2; } else if (set[other2].style == NONE) { if (set[other1].style == NONE || set[other1].style == VOLUME) error->all("Fix deform volume setting is invalid"); set[i].substyle = ONE_FROM_ONE; set[i].fixed = other2; set[i].dynamic1 = other1; } else if (set[other1].style == VOLUME) { if (set[other2].style == NONE || set[other2].style == VOLUME) error->all("Fix deform volume setting is invalid"); set[i].substyle = TWO_FROM_ONE; set[i].fixed = other1; set[i].dynamic1 = other2; } else if (set[other2].style == VOLUME) { if (set[other1].style == NONE || set[other1].style == VOLUME) error->all("Fix deform volume setting is invalid"); set[i].substyle = TWO_FROM_ONE; set[i].fixed = other2; set[i].dynamic1 = other1; } else { set[i].substyle = ONE_FROM_TWO; set[i].dynamic2 = other1; set[i].dynamic2 = other2; } } // initial settings // reneighboring only forced if flips will occur due to shape changes if (set[3].style || set[4].style || set[5].style) force_reneighbor = 1; next_reneighbor = -1; nrigid = 0; rfix = NULL; flip = 0; } /* ---------------------------------------------------------------------- */ FixDeform::~FixDeform() { delete [] set; delete [] rfix; } /* ---------------------------------------------------------------------- */ int FixDeform::setmask() { int mask = 0; mask |= PRE_EXCHANGE; mask |= END_OF_STEP; return mask; } /* ---------------------------------------------------------------------- */ void FixDeform::init() { // error if more than one fix deform // domain, fix nvt/sllod, compute temp/deform only work on single h_rate int count = 0; for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) count++; if (count > 1) error->all("More than one fix deform"); // Kspace setting if (force->kspace) kspace_flag = 1; else kspace_flag = 0; // elapsed time for entire simulation, including multiple runs if defined // cannot be 0.0 since can't deform a finite distance in 0 time double delt = (update->endstep - update->beginstep) * update->dt; if (delt == 0.0) error->all("Cannot use fix deform for 0 timestep run"); // set final values for box size and shape for VEL,ERATE,TRATE // now possible since length of run is known for (int i = 0; i < 3; i++) { if (set[i].style == VEL) { set[i].lo_stop = set[i].lo_start - 0.5*delt*set[i].vel; set[i].hi_stop = set[i].hi_start + 0.5*delt*set[i].vel; } else if (set[i].style == ERATE) { set[i].lo_stop = set[i].lo_start - 0.5*delt*set[i].rate * (set[i].hi_start-set[i].lo_start); set[i].hi_stop = set[i].hi_start + 0.5*delt*set[i].rate * (set[i].hi_start-set[i].lo_start); if (set[i].hi_stop <= set[i].lo_stop) error->all("Final box dimension due to fix deform is < 0.0"); } else if (set[i].style == TRATE) { set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - 0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt)); set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + 0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt)); } } for (int i = 3; i < 6; i++) { if (set[i].style == VEL) { set[i].tilt_stop = set[i].tilt_start + delt*set[i].vel; } else if (set[i].style == ERATE) { if (i == 3) set[i].tilt_stop = set[i].tilt_start + delt*set[i].rate * (set[2].hi_start-set[2].lo_start); if (i == 4) set[i].tilt_stop = set[i].tilt_start + delt*set[i].rate * (set[2].hi_start-set[2].lo_start); if (i == 5) set[i].tilt_stop = set[i].tilt_start + delt*set[i].rate * (set[1].hi_start-set[1].lo_start); } else if (set[i].style == TRATE) { set[i].tilt_stop = set[i].tilt_start * pow(1.0+set[i].rate,delt); } } // if using tilt TRATE, then initial tilt must be non-zero for (int i = 3; i < 6; i++) if (set[i].style == TRATE) if ((i == 5 && domain->xy == 0.0) || (i == 4 && domain->xz == 0.0) || (i == 3 && domain->yz == 0.0)) error->all("Cannot use fix deform trate on a box with zero tilt"); // if yz changes and will cause box flip, then xy cannot be changing // this is b/c the flips would induce continuous changes in xz // in order to keep the edge vectors of the flipped shape matrix // a linear combination of the edge vectors of the unflipped shape matrix if (set[3].style && set[5].style) if (set[3].tilt_stop < -0.5*(set[1].hi_start-set[1].lo_start) || set[3].tilt_stop > 0.5*(set[1].hi_start-set[1].lo_start) || set[3].tilt_stop < -0.5*(set[1].hi_stop-set[1].lo_stop) || set[3].tilt_stop > 0.5*(set[1].hi_stop-set[1].lo_stop)) error->all("Fix deform is changing yz by too much with changing xy"); // set domain->h_rate values for use by domain and other fixes/computes // initialize all rates to 0.0 // cannot set h_rate now for TRATE,VOLUME styles since not constant h_rate = domain->h_rate; h_ratelo = domain->h_ratelo; for (int i = 0; i < 3; i++) { h_rate[i] = h_ratelo[i] = 0.0; if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == SCALE || set[i].style == VEL || set[i].style == ERATE) { double dlo_dt = (set[i].lo_stop - set[i].lo_start) / delt; double dhi_dt = (set[i].hi_stop - set[i].hi_start) / delt; h_rate[i] = dhi_dt - dlo_dt; h_ratelo[i] = dlo_dt; } } for (int i = 3; i < 6; i++) { h_rate[i] = 0.0; if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == VEL || set[i].style == ERATE) h_rate[i] = (set[i].tilt_stop - set[i].tilt_start) / delt; } // detect if any rigid fixes exist so rigid bodies can be rescaled // rfix[] = indices to each fix rigid delete [] rfix; nrigid = 0; rfix = NULL; for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) nrigid++; if (nrigid) { rfix = new int[nrigid]; nrigid = 0; for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; } } /* ---------------------------------------------------------------------- box flipped on previous step perform irregular comm to migrate atoms to new procs reset box tilts for flipped config and create new box in domain remap to put far-away atoms back into new box perform irregular on atoms in lamda coords to get atoms to new procs force reneighboring on next timestep ------------------------------------------------------------------------- */ void FixDeform::pre_exchange() { if (flip == 0) return; domain->yz = set[3].tilt_target = set[3].tilt_flip; domain->xz = set[4].tilt_target = set[4].tilt_flip; domain->xy = set[5].tilt_target = set[5].tilt_flip; domain->set_global_box(); domain->set_local_box(); double **x = atom->x; int *image = atom->image; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); domain->x2lamda(atom->nlocal); comm->irregular(); domain->lamda2x(atom->nlocal); flip = 0; } /* ---------------------------------------------------------------------- */ void FixDeform::end_of_step() { int i; double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; // set new box size // for TRATE, set target directly based on current time and set h_rate // for others except VOLUME, target is linear value between start and stop for (i = 0; i < 3; i++) { if (set[i].style == NONE) continue; if (set[i].style == TRATE) { double delt = (update->ntimestep - update->beginstep) * update->dt; set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - 0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt)); set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + 0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt)); h_rate[i] = set[i].rate * domain->h[i]; h_ratelo[i] = -0.5*h_rate[i]; } else if (set[i].style != VOLUME) { set[i].lo_target = set[i].lo_start + delta*(set[i].lo_stop - set[i].lo_start); set[i].hi_target = set[i].hi_start + delta*(set[i].hi_stop - set[i].hi_start); } } // set new box size for VOLUME dims that are linked to other dims // also need to set h_rate for these dims for (int i = 0; i < 3; i++) { if (set[i].style != VOLUME) continue; if (set[i].substyle == ONE_FROM_ONE) { set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - 0.5*(set[i].vol_start / (set[set[i].dynamic1].hi_target - set[set[i].dynamic1].lo_target) / (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + 0.5*(set[i].vol_start / (set[set[i].dynamic1].hi_target - set[set[i].dynamic1].lo_target) / (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); } else if (set[i].substyle == ONE_FROM_TWO) { set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - 0.5*(set[i].vol_start / (set[set[i].dynamic1].hi_target - set[set[i].dynamic1].lo_target) / (set[set[i].dynamic2].hi_target - set[set[i].dynamic2].lo_target)); set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + 0.5*(set[i].vol_start / (set[set[i].dynamic1].hi_target - set[set[i].dynamic1].lo_target) / (set[set[i].dynamic2].hi_target - set[set[i].dynamic2].lo_target)); } else if (set[i].substyle == TWO_FROM_ONE) { set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - 0.5*sqrt(set[i].vol_start / (set[set[i].dynamic1].hi_target - set[set[i].dynamic1].lo_target) / (set[set[i].fixed].hi_start - set[set[i].fixed].lo_start) * (set[i].hi_start - set[i].lo_start)); set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + 0.5*sqrt(set[i].vol_start / (set[set[i].dynamic1].hi_target - set[set[i].dynamic1].lo_target) / (set[set[i].fixed].hi_start - set[set[i].fixed].lo_start) * (set[i].hi_start - set[i].lo_start)); } } // for triclinic, set new box shape // for TRATE, set target directly based on current time and set h_rate // for other styles, target is linear value between start and stop values if (triclinic) { double *h = domain->h; for (i = 3; i < 6; i++) { if (set[i].style == NONE) continue; if (set[i].style == TRATE) { double delt = (update->ntimestep - update->beginstep) * update->dt; set[i].tilt_target = set[i].tilt_start * pow(1.0+set[i].rate,delt); h_rate[i] = set[i].rate * domain->h[i]; } else if (set[i].style) { set[i].tilt_target = set[i].tilt_start + delta*(set[i].tilt_stop - set[i].tilt_start); } // tilt_target can be large positive or large negative value // add/subtract box lengths until tilt_target is closest to current value int idenom; if (i == 5) idenom = 0; else if (i == 4) idenom = 0; else if (i == 3) idenom = 1; double denom = set[idenom].hi_target - set[idenom].lo_target; double current = h[i]/h[idenom]; while (set[i].tilt_target/denom - current > 0.0) set[i].tilt_target -= denom; while (set[i].tilt_target/denom - current < 0.0) set[i].tilt_target += denom; if (fabs(set[i].tilt_target/denom - 1.0 - current) < fabs(set[i].tilt_target/denom - current)) set[i].tilt_target -= denom; } } // if any tilt targets exceed bounds, set flip flag and new tilt_flip values // flip will be performed on next timestep before reneighboring // when yz flips and xy is non-zero, xz must also change // this is to keep the edge vectors of the flipped shape matrix // a linear combination of the edge vectors of the unflipped shape matrix if (triclinic) { double xprd = set[0].hi_target - set[0].lo_target; double yprd = set[1].hi_target - set[1].lo_target; if (set[3].tilt_target < -0.5*yprd || set[3].tilt_target > 0.5*yprd || set[4].tilt_target < -0.5*xprd || set[4].tilt_target > 0.5*xprd || set[5].tilt_target < -0.5*xprd || set[5].tilt_target > 0.5*xprd) { flip = 1; next_reneighbor = update->ntimestep + 1; set[3].tilt_flip = set[3].tilt_target; set[4].tilt_flip = set[4].tilt_target; set[5].tilt_flip = set[5].tilt_target; if (set[3].tilt_flip < -0.5*yprd) { set[3].tilt_flip += yprd; set[4].tilt_flip += set[5].tilt_flip; } else if (set[3].tilt_flip >= 0.5*yprd) { set[3].tilt_flip -= yprd; set[4].tilt_flip -= set[5].tilt_flip; } if (set[4].tilt_flip < -0.5*xprd) set[4].tilt_flip += xprd; if (set[4].tilt_flip > 0.5*xprd) set[4].tilt_flip -= xprd; if (set[5].tilt_flip < -0.5*xprd) set[5].tilt_flip += xprd; if (set[5].tilt_flip > 0.5*xprd) set[5].tilt_flip -= xprd; } } // convert atoms to lamda coords if (remapflag == X_REMAP) { double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) domain->x2lamda(x[i],x[i]); } // reset global and local box to new size/shape domain->boxlo[0] = set[0].lo_target; domain->boxlo[1] = set[1].lo_target; domain->boxlo[2] = set[2].lo_target; domain->boxhi[0] = set[0].hi_target; domain->boxhi[1] = set[1].hi_target; domain->boxhi[2] = set[2].hi_target; if (triclinic) { domain->yz = set[3].tilt_target; domain->xz = set[4].tilt_target; domain->xy = set[5].tilt_target; } domain->set_global_box(); domain->set_local_box(); - // convert affine atoms back to box coords + // convert atoms back to box coords if (remapflag == X_REMAP) { double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) domain->lamda2x(x[i],x[i]); } // remap centers-of-mass of rigid bodies // needs to be new call, not dilate // check where else fix dilate is called from /* if (nrigid) for (i = 0; i < nrigid; i++) modify->fix[rfix[i]]->dilate(2,oldlo,oldhi,newlo,newhi); */ // redo KSpace coeffs since box has changed if (kspace_flag) force->kspace->setup(); } /* ---------------------------------------------------------------------- */ void FixDeform::options(int narg, char **arg) { if (narg < 0) error->all("Illegal fix deform command"); remapflag = X_REMAP; scaleflag = 1; int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"remap") == 0) { if (iarg+2 > narg) error->all("Illegal fix deform command"); if (strcmp(arg[iarg+1],"x") == 0) remapflag = X_REMAP; else if (strcmp(arg[iarg+1],"v") == 0) remapflag = V_REMAP; else if (strcmp(arg[iarg+1],"none") == 0) remapflag = NO_REMAP; else error->all("Illegal fix deform command"); iarg += 2; } else if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all("Illegal fix deform command"); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; else error->all("Illegal fix deform command"); iarg += 2; } else error->all("Illegal fix deform command"); } } diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index aded32acb..349615246 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -1,673 +1,675 @@ /* ---------------------------------------------------------------------- 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: James Fischer (High Performance Technologies, Inc) Vincent Natoli (Stone Ridge Technology) David Richie (Stone Ridge Technology) ------------------------------------------------------------------------- */ #include "math.h" #include "string.h" #include "ctype.h" #include "pair_hybrid.h" #include "atom.h" #include "force.h" #include "pair.h" #include "neighbor.h" #include "update.h" #include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define NEIGHEXTRA 10000 #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) /* ---------------------------------------------------------------------- */ PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp) { nstyles = 0; } /* ---------------------------------------------------------------------- */ PairHybrid::~PairHybrid() { if (nstyles) { for (int m = 0; m < nstyles; m++) delete styles[m]; delete [] styles; for (int m = 0; m < nstyles; m++) delete [] keywords[m]; delete [] keywords; } if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_int_array(map); memory->destroy_2d_double_array(cutsq); delete [] nnlist; delete [] maxneigh; for (int m = 0; m < nstyles; m++) memory->sfree(nlist[m]); delete [] nlist; delete [] nnlist_full; delete [] maxneigh_full; for (int m = 0; m < nstyles; m++) memory->sfree(nlist_full[m]); delete [] nlist_full; for (int m = 0; m < nstyles; m++) delete [] firstneigh[m]; delete [] firstneigh; for (int m = 0; m < nstyles; m++) delete [] numneigh[m]; delete [] numneigh; for (int m = 0; m < nstyles; m++) delete [] firstneigh_full[m]; delete [] firstneigh_full; for (int m = 0; m < nstyles; m++) delete [] numneigh_full[m]; delete [] numneigh_full; } } /* ---------------------------------------------------------------------- */ void PairHybrid::compute(int eflag, int vflag) { int i,j,k,m,n,jfull,nneigh; int *neighs,*mapi; double **f_original; // save ptrs to original neighbor lists int **firstneigh_original = neighbor->firstneigh; int *numneigh_original = neighbor->numneigh; int **firstneigh_full_original = neighbor->firstneigh_full; int *numneigh_full_original = neighbor->numneigh_full; // if this is re-neighbor step, create sub-style lists if (neighbor->ago == 0) { int *type = atom->type; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; // realloc per-atom per-style firstneigh/numneigh half/full if necessary if (nlocal > maxlocal) { maxlocal = nlocal; if (neigh_half_every) { for (m = 0; m < nstyles; m++) { delete [] firstneigh[m]; delete [] numneigh[m]; } for (m = 0; m < nstyles; m++) { firstneigh[m] = new int*[maxlocal]; numneigh[m] = new int[maxlocal]; } } if (neigh_full_every) { for (m = 0; m < nstyles; m++) { delete [] firstneigh_full[m]; delete [] numneigh_full[m]; } for (m = 0; m < nstyles; m++) { firstneigh_full[m] = new int*[maxlocal]; numneigh_full[m] = new int[maxlocal]; } } } // nnlist[] = length of each sub-style half list // nnlist_full[] = length of each sub-style full list // count from half and/or full list depending on what sub-styles use for (m = 0; m < nstyles; m++) nnlist[m] = nnlist_full[m] = 0; if (neigh_half_every && neigh_full_every) { for (i = 0; i < nlocal; i++) { mapi = map[type[i]]; neighs = firstneigh_original[i]; nneigh = numneigh_original[i]; for (k = 0; k < nneigh; k++) { j = neighs[k]; if (j >= nall) j %= nall; m = mapi[type[j]]; if (styles[m] && styles[m]->neigh_half_every) nnlist[m]++; } neighs = firstneigh_full_original[i]; nneigh = numneigh_full_original[i]; for (k = 0; k < nneigh; k++) { j = neighs[k]; if (j >= nall) j %= nall; m = mapi[type[j]]; if (styles[m] && styles[m]->neigh_full_every) nnlist_full[m]++; } } } else if (neigh_half_every) { for (i = 0; i < nlocal; i++) { mapi = map[type[i]]; neighs = firstneigh_original[i]; nneigh = numneigh_original[i]; for (k = 0; k < nneigh; k++) { j = neighs[k]; if (j >= nall) j %= nall; nnlist[mapi[type[j]]]++; } } } else if (neigh_full_every) { for (i = 0; i < nlocal; i++) { mapi = map[type[i]]; neighs = firstneigh_full_original[i]; nneigh = numneigh_full_original[i]; for (k = 0; k < nneigh; k++) { j = neighs[k]; if (j >= nall) j %= nall; nnlist_full[mapi[type[j]]]++; } } } // realloc sub-style nlist and nlist_full if necessary if (neigh_half_every) { for (m = 0; m < nstyles; m++) { if (nnlist[m] > maxneigh[m]) { memory->sfree(nlist[m]); maxneigh[m] = nnlist[m] + NEIGHEXTRA; nlist[m] = (int *) memory->smalloc(maxneigh[m]*sizeof(int),"pair_hybrid:nlist"); } nnlist[m] = 0; } } if (neigh_full_every) { for (m = 0; m < nstyles; m++) { if (nnlist_full[m] > maxneigh_full[m]) { memory->sfree(nlist_full[m]); maxneigh_full[m] = nnlist_full[m] + NEIGHEXTRA; nlist_full[m] = (int *) memory->smalloc(maxneigh_full[m]*sizeof(int), "pair_hybrid:nlist_full"); } nnlist_full[m] = 0; } } // load sub-style half/full list with values from original lists // load from half and/or full list depending on what sub-styles use if (neigh_half_every && neigh_full_every) { for (i = 0; i < nlocal; i++) { for (m = 0; m < nstyles; m++) { firstneigh[m][i] = &nlist[m][nnlist[m]]; numneigh[m][i] = nnlist[m]; firstneigh_full[m][i] = &nlist_full[m][nnlist_full[m]]; numneigh_full[m][i] = nnlist_full[m]; } mapi = map[type[i]]; neighs = firstneigh_original[i]; nneigh = numneigh_original[i]; for (k = 0; k < nneigh; k++) { j = jfull = neighs[k]; if (j >= nall) j %= nall; m = mapi[type[j]]; if (styles[m] && styles[m]->neigh_half_every) nlist[m][nnlist[m]++] = jfull; } neighs = firstneigh_full_original[i]; nneigh = numneigh_full_original[i]; for (k = 0; k < nneigh; k++) { j = jfull = neighs[k]; if (j >= nall) j %= nall; m = mapi[type[j]]; if (styles[m] && styles[m]->neigh_full_every) nlist_full[m][nnlist_full[m]++] = jfull; } for (m = 0; m < nstyles; m++) { numneigh[m][i] = nnlist[m] - numneigh[m][i]; numneigh_full[m][i] = nnlist_full[m] - numneigh_full[m][i]; } } } else if (neigh_half_every) { for (i = 0; i < nlocal; i++) { for (m = 0; m < nstyles; m++) { firstneigh[m][i] = &nlist[m][nnlist[m]]; numneigh[m][i] = nnlist[m]; } mapi = map[type[i]]; neighs = firstneigh_original[i]; nneigh = numneigh_original[i]; for (k = 0; k < nneigh; k++) { j = jfull = neighs[k]; if (j >= nall) j %= nall; m = mapi[type[j]]; nlist[m][nnlist[m]++] = jfull; } for (m = 0; m < nstyles; m++) numneigh[m][i] = nnlist[m] - numneigh[m][i]; } } else if (neigh_full_every) { for (i = 0; i < nlocal; i++) { for (m = 0; m < nstyles; m++) { firstneigh_full[m][i] = &nlist_full[m][nnlist_full[m]]; numneigh_full[m][i] = nnlist_full[m]; } mapi = map[type[i]]; neighs = firstneigh_full_original[i]; nneigh = numneigh_full_original[i]; for (k = 0; k < nneigh; k++) { j = jfull = neighs[k]; if (j >= nall) j %= nall; m = mapi[type[j]]; nlist_full[m][nnlist_full[m]++] = jfull; } for (m = 0; m < nstyles; m++) numneigh_full[m][i] = nnlist_full[m] - numneigh_full[m][i]; } } } // call each sub-style's compute function // set neighbor->firstneigh/numneigh to sub-style lists before call // set half or full or both depending on what sub-style uses // for vflag = 1: // sub-style accumulates in its virial[6] // sum sub-style virial[6] to hybrid's virial[6] // for vflag = 2: // set atom->f to update->f_pair so sub-style will sum its f to f_pair // call sub-style compute() with vflag % 2 to prevent sub-style // from calling virial_compute() // reset atom->f to stored f_original // call hybrid virial_compute() which will use update->f_pair // accumulate sub-style energy,virial in hybrid's energy,virial eng_vdwl = eng_coul = 0.0; if (vflag) for (n = 0; n < 6; n++) virial[n] = 0.0; if (vflag == 2) { f_original = atom->f; atom->f = update->f_pair; } for (m = 0; m < nstyles; m++) { if (styles[m] == NULL) continue; if (styles[m]->neigh_half_every) { neighbor->firstneigh = firstneigh[m]; neighbor->numneigh = numneigh[m]; } if (styles[m]->neigh_full_every) { neighbor->firstneigh_full = firstneigh_full[m]; neighbor->numneigh_full = numneigh_full[m]; } styles[m]->compute(eflag,vflag % 2); if (eflag) { eng_vdwl += styles[m]->eng_vdwl; eng_coul += styles[m]->eng_coul; } if (vflag == 1) for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n]; } if (vflag == 2) { atom->f = f_original; virial_compute(); } // restore ptrs to original neighbor lists neighbor->firstneigh = firstneigh_original; neighbor->numneigh = numneigh_original; neighbor->firstneigh_full = firstneigh_full_original; neighbor->numneigh_full = numneigh_full_original; } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairHybrid::allocate() { allocated = 1; int n = atom->ntypes; map = memory->create_2d_int_array(n+1,n+1,"pair:map"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) map[i][j] = -1; setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); nnlist = new int[nstyles]; maxneigh = new int[nstyles]; nlist = new int*[nstyles]; for (int m = 0; m < nstyles; m++) maxneigh[m] = 0; for (int m = 0; m < nstyles; m++) nlist[m] = NULL; nnlist_full = new int[nstyles]; maxneigh_full = new int[nstyles]; nlist_full = new int*[nstyles]; for (int m = 0; m < nstyles; m++) maxneigh_full[m] = 0; for (int m = 0; m < nstyles; m++) nlist_full[m] = NULL; maxlocal = 0; firstneigh = new int**[nstyles]; numneigh = new int*[nstyles]; for (int m = 0; m < nstyles; m++) firstneigh[m] = NULL; for (int m = 0; m < nstyles; m++) numneigh[m] = NULL; firstneigh_full = new int**[nstyles]; numneigh_full = new int*[nstyles]; for (int m = 0; m < nstyles; m++) firstneigh_full[m] = NULL; for (int m = 0; m < nstyles; m++) numneigh_full[m] = NULL; } /* ---------------------------------------------------------------------- create one pair style for each arg in list ------------------------------------------------------------------------- */ void PairHybrid::settings(int narg, char **arg) { int i,m,istyle; if (narg < 1) error->all("Illegal pair_style command"); // delete old lists, since cannot just change settings if (nstyles) { for (m = 0; m < nstyles; m++) delete styles[m]; delete [] styles; for (m = 0; m < nstyles; m++) delete [] keywords[m]; delete [] keywords; } if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_int_array(map); memory->destroy_2d_double_array(cutsq); } allocated = 0; // count sub-styles by skipping numeric args // one exception is 1st arg of style "table", which is non-numeric word nstyles = i = 0; while (i < narg) { if (strcmp(arg[i],"table") == 0) i++; i++; while (i < narg && !isalpha(arg[i][0])) i++; nstyles++; } // allocate list of sub-styles styles = new Pair*[nstyles]; keywords = new char*[nstyles]; // allocate each sub-style and call its settings() with subset of args // define subset of sub-styles by skipping numeric args // one exception is 1st arg of style "table", which is non-numeric word nstyles = i = 0; while (i < narg) { for (m = 0; m < nstyles; m++) if (strcmp(arg[i],keywords[m]) == 0) error->all("Pair style hybrid cannot use same pair style twice"); if (strcmp(arg[i],"hybrid") == 0) error->all("Pair style hybrid cannot have hybrid as an argument"); styles[nstyles] = force->new_pair(arg[i]); keywords[nstyles] = new char[strlen(arg[i])+1]; strcpy(keywords[nstyles],arg[i]); istyle = i; if (strcmp(arg[i],"table") == 0) i++; i++; while (i < narg && !isalpha(arg[i][0])) i++; if (styles[nstyles]) styles[nstyles]->settings(i-istyle-1,&arg[istyle+1]); nstyles++; } // set comm_forward, comm_reverse to max of any sub-style for (m = 0; m < nstyles; m++) { if (styles[m]) comm_forward = MAX(comm_forward,styles[m]->comm_forward); if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse); } // neigh_every = 1 if any sub-style = 1 neigh_half_every = neigh_full_every = 0; for (m = 0; m < nstyles; m++) { if (styles[m] && styles[m]->neigh_half_every) neigh_half_every = 1; if (styles[m] && styles[m]->neigh_full_every) neigh_full_every = 1; } // single_enable = 0 if any sub-style = 0 for (m = 0; m < nstyles; m++) if (styles[m] && styles[m]->single_enable == 0) single_enable = 0; } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairHybrid::coeff(int narg, char **arg) { if (narg < 3) 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); // 3rd arg = pair style name int m; for (m = 0; m < nstyles; m++) if (strcmp(arg[2],keywords[m]) == 0) break; if (m == nstyles) error->all("Pair coeff for hybrid has invalid style"); // move 1st/2nd args to 2nd/3rd args sprintf(arg[2],"%s",arg[1]); sprintf(arg[1],"%s",arg[0]); // invoke sub-style coeff() starting with 1st arg if (styles[m]) styles[m]->coeff(narg-1,&arg[1]); // if pair style only allows one pair coeff call (with * * and type mapping) // then unset any old setflag/map assigned to that style first // in case pair coeff for this sub-style is being called for 2nd time if (styles[m] && styles[m]->one_coeff) for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) if (map[i][j] == m) { map[i][j] = -1; setflag[i][j] = 0; } // set hybrid map & setflag only if substyle set its setflag // if sub-style is NULL for "none", still set setflag int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { if (styles[m] == NULL || styles[m]->setflag[i][j]) { map[i][j] = m; setflag[i][j] = 1; count++; } } } if (count == 0) error->all("Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairHybrid::init_one(int i, int j) { // if i,j is set explicity, call its sub-style // if i,j is not set and i,i sub-style = j,j sub-style // then set map[i][j] to this sub-style and call sub-style for init/mixing // else i,j has not been set by user // check for special case = style none double cut = 0.0; if (setflag[i][j]) { if (styles[map[i][j]]) { cut = styles[map[i][j]]->init_one(i,j); styles[map[i][j]]->cutsq[i][j] = styles[map[i][j]]->cutsq[j][i] = cut*cut; if (tail_flag) { etail_ij = styles[map[i][j]]->etail_ij; ptail_ij = styles[map[i][j]]->ptail_ij; } } } else if (map[i][i] == map[j][j]) { map[i][j] = map[i][i]; if (styles[map[i][j]]) { cut = styles[map[i][j]]->init_one(i,j); styles[map[i][j]]->cutsq[i][j] = styles[map[i][j]]->cutsq[j][i] = cut*cut; if (tail_flag) { etail_ij = styles[map[i][j]]->etail_ij; ptail_ij = styles[map[i][j]]->ptail_ij; } } } else error->one("All pair coeffs are not set"); map[j][i] = map[i][j]; return cut; } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairHybrid::init_style() { for (int m = 0; m < nstyles; m++) if (styles[m]) styles[m]->init_style(); } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairHybrid::write_restart(FILE *fp) { fwrite(&nstyles,sizeof(int),1,fp); // each sub-style writes its settings int n; for (int m = 0; m < nstyles; m++) { n = strlen(keywords[m]) + 1; fwrite(&n,sizeof(int),1,fp); fwrite(keywords[m],sizeof(char),n,fp); if (styles[m]) styles[m]->write_restart_settings(fp); } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairHybrid::read_restart(FILE *fp) { allocate(); int me = comm->me; if (me == 0) fread(&nstyles,sizeof(int),1,fp); MPI_Bcast(&nstyles,1,MPI_INT,0,world); styles = new Pair*[nstyles]; keywords = new char*[nstyles]; // each sub-style is created via new_pair() and reads its settings int n; for (int m = 0; m < nstyles; m++) { if (me == 0) fread(&n,sizeof(int),1,fp); MPI_Bcast(&n,1,MPI_INT,0,world); keywords[m] = new char[n]; if (me == 0) fread(keywords[m],sizeof(char),n,fp); MPI_Bcast(keywords[m],n,MPI_CHAR,0,world); styles[m] = force->new_pair(keywords[m]); if (styles[m]) styles[m]->read_restart_settings(fp); } } /* ---------------------------------------------------------------------- */ void PairHybrid::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, int eflag, One &one) { if (map[itype][jtype] == -1) error->one("Invoked pair single on pair style none"); styles[map[itype][jtype]]-> single(i,j,itype,jtype,rsq,factor_coul,factor_lj,eflag,one); } /* ---------------------------------------------------------------------- */ void PairHybrid::single_embed(int i, int itype, double &phi) { if (map[itype][itype] == -1) error->one("Invoked pair single on pair style none"); styles[map[itype][itype]]->single_embed(i,itype,phi); } /* ---------------------------------------------------------------------- modify parameters of the pair style - simply pass command args to each sub-style of hybrid + call modify_params of PairHybrid + also pass command args to each sub-style of hybrid ------------------------------------------------------------------------- */ void PairHybrid::modify_params(int narg, char **arg) { + Pair::modify_params(narg,arg); for (int m = 0; m < nstyles; m++) if (styles[m]) styles[m]->modify_params(narg,arg); } /* ---------------------------------------------------------------------- memory usage of sub-style firstneigh, numneigh, neighbor list add in memory usage of each sub-style itself ------------------------------------------------------------------------- */ int PairHybrid::memory_usage() { int bytes = nstyles*maxlocal * (sizeof(int *) + sizeof(int)); for (int m = 0; m < nstyles; m++) bytes += maxneigh[m] * sizeof(int); for (int m = 0; m < nstyles; m++) if (styles[m]) bytes += styles[m]->memory_usage(); return bytes; } diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index 7d0c40923..dc342a13b 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -1,654 +1,654 @@ /* ---------------------------------------------------------------------- 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 "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_lj_cut.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "update.h" #include "integrate.h" #include "respa.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)) /* ---------------------------------------------------------------------- */ PairLJCut::PairLJCut(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; } /* ---------------------------------------------------------------------- */ PairLJCut::~PairLJCut() { if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); memory->destroy_2d_double_array(cut); memory->destroy_2d_double_array(epsilon); memory->destroy_2d_double_array(sigma); memory->destroy_2d_double_array(lj1); memory->destroy_2d_double_array(lj2); memory->destroy_2d_double_array(lj3); memory->destroy_2d_double_array(lj4); memory->destroy_2d_double_array(offset); } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute(int eflag, int vflag) { int i,j,k,numneigh,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcelj,fforce,factor_lj,philj; int *neighs; double **f; eng_vdwl = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; if (vflag == 2) f = update->f_pair; else f = atom->f; double **x = atom->x; int *type = atom->type; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; // loop over neighbors of my atoms for (i = 0; i < nlocal; i++) { xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; neighs = neighbor->firstneigh[i]; numneigh = neighbor->numneigh[i]; for (k = 0; k < numneigh; k++) { j = neighs[k]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fforce = factor_lj*forcelj*r2inv; f[i][0] += delx*fforce; f[i][1] += dely*fforce; f[i][2] += delz*fforce; if (newton_pair || j < nlocal) { f[j][0] -= delx*fforce; f[j][1] -= dely*fforce; f[j][2] -= delz*fforce; } if (eflag) { philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; if (newton_pair || j < nlocal) eng_vdwl += factor_lj*philj; else eng_vdwl += 0.5*factor_lj*philj; } if (vflag == 1) { if (newton_pair == 0 && j >= nlocal) fforce *= 0.5; virial[0] += delx*delx*fforce; virial[1] += dely*dely*fforce; virial[2] += delz*delz*fforce; virial[3] += delx*dely*fforce; virial[4] += delx*delz*fforce; virial[5] += dely*delz*fforce; } } } } if (vflag == 2) virial_compute(); } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_inner() { int i,j,k,numneigh,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcelj,fforce,factor_lj; double rsw; int *neighs; double **f; f = atom->f; double **x = atom->x; int *type = atom->type; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; // loop over neighbors of my atoms for (i = 0; i < nlocal; i++) { xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; neighs = neighbor->firstneigh_inner[i]; numneigh = neighbor->numneigh_inner[i]; for (k = 0; k < numneigh; k++) { j = neighs[k]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cut_out_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; jtype = type[j]; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fforce = factor_lj*forcelj*r2inv; if (rsq > cut_out_on_sq) { rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fforce *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw); } f[i][0] += delx*fforce; f[i][1] += dely*fforce; f[i][2] += delz*fforce; if (newton_pair || j < nlocal) { f[j][0] -= delx*fforce; f[j][1] -= dely*fforce; f[j][2] -= delz*fforce; } } } } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_middle() { int i,j,k,numneigh,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcelj,fforce,factor_lj; double rsw; int *neighs; double **f; f = atom->f; double **x = atom->x; int *type = atom->type; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; double cut_out_on = cut_respa[2]; double cut_out_off = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; // loop over neighbors of my atoms for (i = 0; i < nlocal; i++) { xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; neighs = neighbor->firstneigh_middle[i]; numneigh = neighbor->numneigh_middle[i]; for (k = 0; k < numneigh; k++) { j = neighs[k]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; jtype = type[j]; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fforce = factor_lj*forcelj*r2inv; if (rsq < cut_in_on_sq) { rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; fforce *= rsw*rsw*(3.0 - 2.0*rsw); } if (rsq > cut_out_on_sq) { rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fforce *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); } f[i][0] += delx*fforce; f[i][1] += dely*fforce; f[i][2] += delz*fforce; if (newton_pair || j < nlocal) { f[j][0] -= delx*fforce; f[j][1] -= dely*fforce; f[j][2] -= delz*fforce; } } } } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_outer(int eflag, int vflag) { int i,j,k,numneigh,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcelj,fforce,factor_lj,philj; double rsw; int *neighs; double **f; eng_vdwl = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; f = atom->f; double **x = atom->x; int *type = atom->type; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; // loop over neighbors of my atoms for (i = 0; i < nlocal; i++) { xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; neighs = neighbor->firstneigh[i]; numneigh = neighbor->numneigh[i]; for (k = 0; k < numneigh; k++) { j = neighs[k]; if (j < nall) factor_lj = 1.0; else { factor_lj = special_lj[j/nall]; j %= nall; } delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { if (rsq > cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fforce = factor_lj*forcelj*r2inv; if (rsq < cut_in_on_sq) { rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; fforce *= rsw*rsw*(3.0 - 2.0*rsw); } f[i][0] += delx*fforce; f[i][1] += dely*fforce; f[i][2] += delz*fforce; if (newton_pair || j < nlocal) { f[j][0] -= delx*fforce; f[j][1] -= dely*fforce; f[j][2] -= delz*fforce; } } if (eflag) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; if (newton_pair || j < nlocal) eng_vdwl += factor_lj*philj; else eng_vdwl += 0.5*factor_lj*philj; } if (vflag) { if (rsq <= cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fforce = factor_lj*forcelj*r2inv; } else if (rsq < cut_in_on_sq) fforce = factor_lj*forcelj*r2inv; if (newton_pair == 0 && j >= nlocal) fforce *= 0.5; virial[0] += delx*delx*fforce; virial[1] += dely*dely*fforce; virial[2] += delz*delz*fforce; virial[3] += delx*dely*fforce; virial[4] += delx*delz*fforce; virial[5] += dely*delz*fforce; } } } } } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairLJCut::allocate() { allocated = 1; int n = atom->ntypes; setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); epsilon = memory->create_2d_double_array(n+1,n+1,"pair:epsilon"); sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); lj1 = memory->create_2d_double_array(n+1,n+1,"pair:lj1"); lj2 = memory->create_2d_double_array(n+1,n+1,"pair:lj2"); lj3 = memory->create_2d_double_array(n+1,n+1,"pair:lj3"); lj4 = memory->create_2d_double_array(n+1,n+1,"pair:lj4"); offset = memory->create_2d_double_array(n+1,n+1,"pair:offset"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairLJCut::settings(int narg, char **arg) { if (narg != 1) error->all("Illegal pair_style command"); cut_global = atof(arg[0]); // 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 PairLJCut::coeff(int narg, char **arg) { if (narg < 4 || narg > 5) 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 epsilon_one = atof(arg[2]); double sigma_one = atof(arg[3]); double cut_one = cut_global; if (narg == 5) cut_one = atof(arg[4]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; 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 for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairLJCut::init_one(int i, int j) { if (setflag[i][j] == 0) { epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], sigma[i][i],sigma[j][j]); sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); cut[i][j] = mix_distance(cut[i][i],cut[j][j]); } lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - + if (offset_flag) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; lj3[j][i] = lj3[i][j]; lj4[j][i] = lj4[i][j]; offset[j][i] = offset[i][j]; // set & error check interior rRESPA cutoff if (strcmp(update->integrate_style,"respa") == 0) { if (((Respa *) update->integrate)->level_inner >= 0) { cut_respa = ((Respa *) update->integrate)->cutoff; if (cut[i][j] < cut_respa[3]) error->all("Pair cutoff < Respa interior cutoff"); } } else cut_respa = NULL; // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce if (tail_flag) { int *type = atom->type; int nlocal = atom->nlocal; double count[2],all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); double PI = 4.0*atan(1.0); double sig2 = sigma[i][j]*sigma[i][j]; double sig6 = sig2*sig2*sig2; double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; double rc6 = rc3*rc3; double rc9 = rc3*rc6; etail_ij = 8.0*PI*all[0]*all[1]*epsilon[i][j] * sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); ptail_ij = 16.0*PI*all[0]*all[1]*epsilon[i][j] * sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); } return cut[i][j]; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJCut::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(&epsilon[i][j],sizeof(double),1,fp); fwrite(&sigma[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJCut::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(&epsilon[i][j],sizeof(double),1,fp); fread(&sigma[i][j],sizeof(double),1,fp); fread(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); 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 PairLJCut::write_restart_settings(FILE *fp) { fwrite(&cut_global,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJCut::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { fread(&cut_global,sizeof(double),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(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- */ void PairLJCut::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, int eflag, One &one) { double r2inv,r6inv,forcelj,philj; r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); one.fforce = factor_lj*forcelj*r2inv; if (eflag) { philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; one.eng_vdwl = factor_lj*philj; one.eng_coul = 0.0; } } diff --git a/src/style_asphere.h b/src/style_asphere.h index e69de29bb..f67f9033e 100644 --- a/src/style_asphere.h +++ b/src/style_asphere.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + 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 AtomInclude +#include "atom_vec_ellipsoid.h" +#endif + +#ifdef AtomClass +AtomStyle(ellipsoid,AtomVecEllipsoid) +# endif + +#ifdef ComputeInclude +#include "compute_temp_asphere.h" +#endif + +#ifdef ComputeClass +ComputeStyle(temp/asphere,ComputeTempAsphere) +#endif + +#ifdef FixInclude +#include "fix_nve_asphere.h" +#include "fix_nvt_asphere.h" +#include "fix_npt_asphere.h" +#endif + +#ifdef FixClass +FixStyle(nve/asphere,FixNVEASphere) +FixStyle(nvt/asphere,FixNVTASphere) +FixStyle(npt/asphere,FixNPTASphere) +#endif + +#ifdef PairInclude +#include "pair_gayberne.h" +#endif + +#ifdef PairClass +PairStyle(gayberne,PairGayBerne) +#endif diff --git a/src/style_class2.h b/src/style_class2.h index e69de29bb..3ba9b32df 100644 --- a/src/style_class2.h +++ b/src/style_class2.h @@ -0,0 +1,56 @@ +/* ---------------------------------------------------------------------- + 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 AngleInclude +#include "angle_class2.h" +#endif + +#ifdef AngleClass +AngleStyle(class2,AngleClass2) +#endif + +#ifdef BondInclude +#include "bond_class2.h" +#endif + +#ifdef BondClass +BondStyle(class2,BondClass2) +#endif + +#ifdef DihedralInclude +#include "dihedral_class2.h" +#endif + +#ifdef DihedralClass +DihedralStyle(class2,DihedralClass2) +#endif + +#ifdef ImproperInclude +#include "improper_class2.h" +#endif + +#ifdef ImproperClass +ImproperStyle(class2,ImproperClass2) +#endif + +#ifdef PairInclude +#include "pair_lj_class2.h" +#include "pair_lj_class2_coul_cut.h" +#include "pair_lj_class2_coul_long.h" +#endif + +#ifdef PairClass +PairStyle(lj/class2,PairLJClass2) +PairStyle(lj/class2/coul/cut,PairLJClass2CoulCut) +PairStyle(lj/class2/coul/long,PairLJClass2CoulLong) +#endif diff --git a/src/style_colloid.h b/src/style_colloid.h index e69de29bb..393605e64 100644 --- a/src/style_colloid.h +++ b/src/style_colloid.h @@ -0,0 +1,20 @@ +/* ---------------------------------------------------------------------- + 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 PairInclude +#include "pair_colloid.h" +#endif + +#ifdef PairClass +PairStyle(colloid,PairColloid) +#endif diff --git a/src/style_dipole.h b/src/style_dipole.h index e69de29bb..7f989773f 100644 --- a/src/style_dipole.h +++ b/src/style_dipole.h @@ -0,0 +1,44 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 AtomInclude +#include "atom_vec_dipole.h" +#endif + +#ifdef AtomClass +AtomStyle(dipole,AtomVecDipole) +#endif + +#ifdef ComputeInclude +#include "compute_temp_dipole.h" +#endif + +#ifdef ComputeClass +ComputeStyle(temp/dipole,ComputeTempDipole) +#endif + +#ifdef FixInclude +#include "fix_nve_dipole.h" +#endif + +#ifdef FixClass +FixStyle(nve/dipole,FixNVEDipole) +#endif + +#ifdef PairInclude +#include "pair_dipole_cut.h" +#endif + +#ifdef PairClass +PairStyle(dipole/cut,PairDipoleCut) +#endif diff --git a/src/style_dpd.h b/src/style_dpd.h index e69de29bb..8ce617c0c 100644 --- a/src/style_dpd.h +++ b/src/style_dpd.h @@ -0,0 +1,28 @@ +/* ---------------------------------------------------------------------- + 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 AtomInclude +#include "atom_vec_dpd.h" +#endif + +#ifdef AtomClass +AtomStyle(dpd,AtomVecDPD) +#endif + +#ifdef PairInclude +#include "pair_dpd.h" +#endif + +#ifdef PairClass +PairStyle(dpd,PairDPD) +#endif diff --git a/src/style_granular.h b/src/style_granular.h index e69de29bb..7b0f4b6a7 100644 --- a/src/style_granular.h +++ b/src/style_granular.h @@ -0,0 +1,50 @@ +/* ---------------------------------------------------------------------- + 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 AtomInclude +#include "atom_vec_granular.h" +#endif + +#ifdef AtomClass +AtomStyle(granular,AtomVecGranular) +# endif + +#ifdef FixInclude +#include "fix_freeze.h" +#include "fix_gran_diag.h" +#include "fix_nve_gran.h" +#include "fix_pour.h" +#include "fix_shear_history.h" +#include "fix_wall_gran.h" +#endif + +#ifdef FixClass +FixStyle(freeze,FixFreeze) +FixStyle(gran/diag,FixGranDiag) +FixStyle(nve/gran,FixNVEGran) +FixStyle(pour,FixPour) +FixStyle(SHEAR_HISTORY,FixShearHistory) +FixStyle(wall/gran,FixWallGran) +#endif + +#ifdef PairInclude +#include "pair_gran_hertzian.h" +#include "pair_gran_history.h" +#include "pair_gran_no_history.h" +#endif + +#ifdef PairClass +PairStyle(gran/hertzian,PairGranHertzian) +PairStyle(gran/history,PairGranHistory) +PairStyle(gran/no_history,PairGranNoHistory) +#endif diff --git a/src/style_opt.h b/src/style_opt.h index e69de29bb..061816e13 100644 --- a/src/style_opt.h +++ b/src/style_opt.h @@ -0,0 +1,30 @@ +/* ---------------------------------------------------------------------- + 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 PairInclude +#include "pair_eam_opt.h" +#include "pair_eam_alloy_opt.h" +#include "pair_eam_fs_opt.h" +#include "pair_lj_charmm_coul_long_opt.h" +#include "pair_lj_cut_opt.h" +#include "pair_morse_opt.h" +#endif + +#ifdef PairClass +PairStyle(eam/opt,PairEAMOpt) +PairStyle(eam/alloy/opt,PairEAMAlloyOpt) +PairStyle(eam/fs/opt,PairEAMFSOpt) +PairStyle(lj/cut/opt,PairLJCutOpt) +PairStyle(lj/charmm/coul/long/opt,PairLJCharmmCoulLongOpt) +PairStyle(morse/opt,PairMorseOpt) +#endif diff --git a/src/style_xtc.h b/src/style_xtc.h index e69de29bb..7110dda31 100644 --- a/src/style_xtc.h +++ b/src/style_xtc.h @@ -0,0 +1,20 @@ +/* ---------------------------------------------------------------------- + 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 DumpInclude +#include "dump_xtc.h" +#endif + +#ifdef DumpClass +DumpStyle(xtc,DumpXTC) +#endif