diff --git a/src/atom.cpp b/src/atom.cpp index 67ea7daeb8..53300484d0 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -31,14 +31,18 @@ #include "update.h" #include "domain.h" #include "group.h" +#include "molecule.h" #include "accelerator_cuda.h" #include "atom_masks.h" +#include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; +using namespace MathConst; #define DELTA 1 +#define DELTA_MOLECULE 1 #define DELTA_MEMSTR 1024 #define EPSILON 1.0e-6 #define CUDA_CHUNK 3000 @@ -106,6 +110,11 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) improper_type = improper_atom1 = improper_atom2 = NULL; improper_atom3 = improper_atom4 = NULL; + // user-defined molecules + + nmolecule = maxmol = 0; + molecules = NULL; + // custom atom arrays nivector = ndvector = 0; @@ -230,6 +239,11 @@ Atom::~Atom() memory->destroy(improper_atom3); memory->destroy(improper_atom4); + // delete user-defined molecules + + for (int i = 0; i < nmolecule; i++) delete molecules[i]; + memory->sfree(molecules); + // delete custom atom arrays for (int i = 0; i < nivector; i++) { @@ -1152,6 +1166,99 @@ int Atom::shape_consistency(int itype, return 1; } +/* ---------------------------------------------------------------------- + add a new molecule template +------------------------------------------------------------------------- */ + +void Atom::add_molecule(int narg, char **arg) +{ + if (narg < 1) error->all(FLERR,"Illegal molecule command"); + if (find_molecule(arg[0]) >= 0) error->all(FLERR,"Reuse of molecule ID"); + + // extend molecule list if necessary + + if (nmolecule == maxmol) { + maxmol += DELTA_MOLECULE; + molecules = (Molecule **) + memory->srealloc(molecules,maxmol*sizeof(Molecule *),"atom::molecules"); + } + + molecules[nmolecule++] = new Molecule(lmp,narg,arg); +} + +/* ---------------------------------------------------------------------- + find which molecule has molecule ID + return -1 if does not exist +------------------------------------------------------------------------- */ + +int Atom::find_molecule(char *id) +{ + int imol; + for (imol = 0; imol < nmolecule; imol++) + if (strcmp(id,molecules[imol]->id) == 0) return imol; + return -1; +} + +/* ---------------------------------------------------------------------- + add info for iatom from molecule template onemol to current atom ilocal +------------------------------------------------------------------------- */ + +void Atom::add_molecule_atom(Molecule *onemol, int iatom, + int ilocal, int offset) +{ + if (onemol->qflag) q[ilocal] = onemol->q[iatom]; + if (onemol->radiusflag) radius[ilocal] = onemol->radius[iatom]; + if (onemol->rmassflag) rmass[ilocal] = onemol->rmass[iatom]; + else if (rmass_flag) + rmass[ilocal] = 4.0*MY_PI/3.0 * + radius[ilocal]*radius[ilocal]*radius[ilocal]; + + if (onemol->bondflag) { + num_bond[ilocal] = onemol->num_bond[iatom]; + for (int i = 0; i < num_bond[ilocal]; i++) { + bond_type[ilocal][i] = onemol->bond_type[iatom][i]; + bond_atom[ilocal][i] = onemol->bond_atom[iatom][i] + offset; + } + } + if (onemol->angleflag) { + num_angle[ilocal] = onemol->num_angle[iatom]; + for (int i = 0; i < num_angle[ilocal]; i++) { + angle_type[ilocal][i] = onemol->angle_type[iatom][i]; + angle_atom1[ilocal][i] = onemol->angle_atom1[iatom][i] + offset; + angle_atom2[ilocal][i] = onemol->angle_atom2[iatom][i] + offset; + angle_atom3[ilocal][i] = onemol->angle_atom3[iatom][i] + offset; + } + } + if (onemol->dihedralflag) { + num_dihedral[ilocal] = onemol->num_dihedral[iatom]; + for (int i = 0; i < num_dihedral[ilocal]; i++) { + dihedral_type[ilocal][i] = onemol->dihedral_type[iatom][i]; + dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset; + dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset; + dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset; + dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset; + } + } + if (onemol->improperflag) { + num_improper[ilocal] = onemol->num_improper[iatom]; + for (int i = 0; i < num_improper[ilocal]; i++) { + improper_type[ilocal][i] = onemol->improper_type[iatom][i]; + improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset; + improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset; + improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset; + improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset; + } + } + + if (onemol->specialflag) { + nspecial[ilocal][0] = onemol->nspecial[iatom][0]; + nspecial[ilocal][1] = onemol->nspecial[iatom][1]; + int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2]; + for (int i = 0; i < n; i++) + special[ilocal][i] = onemol->special[iatom][i] + offset; + } +} + /* ---------------------------------------------------------------------- reorder owned atoms so those in firstgroup appear first called by comm->exchange() if atom_modify first group is set diff --git a/src/atom.h b/src/atom.h index 6b915f4eed..80dadd0c6c 100644 --- a/src/atom.h +++ b/src/atom.h @@ -109,6 +109,11 @@ class Atom : protected Pointers { int cs_flag,csforce_flag,vforce_flag,ervelforce_flag,etag_flag; int rho_flag,e_flag,cv_flag,vest_flag; + // molecules + + int nmolecule,maxmol; + class Molecule **molecules; + // extra peratom info in restart file destined for fix & diag double **extra; @@ -178,6 +183,10 @@ class Atom : protected Pointers { int radius_consistency(int, double &); int shape_consistency(int, double &, double &, double &); + void add_molecule(int, char **); + int find_molecule(char *); + void add_molecule_atom(class Molecule *, int, int, int); + void first_reorder(); void sort(); diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp index 9ee902e0c2..a9342333a5 100644 --- a/src/fix_deposit.cpp +++ b/src/fix_deposit.cpp @@ -17,6 +17,7 @@ #include "fix_deposit.h" #include "atom.h" #include "atom_vec.h" +#include "molecule.h" #include "force.h" #include "update.h" #include "modify.h" @@ -26,11 +27,16 @@ #include "lattice.h" #include "region.h" #include "random_park.h" +#include "math_extra.h" +#include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; +using namespace MathConst; + +enum{ATOM,MOLECULE}; /* ---------------------------------------------------------------------- */ @@ -55,6 +61,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : iregion = -1; idregion = NULL; + mode = ATOM; idnext = 0; globalflag = localflag = 0; lo = hi = deltasq = 0.0; @@ -96,6 +103,24 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Deposition region extends outside simulation box"); } + // error check and further setup for mode = MOLECULE + + if (mode == MOLECULE) { + if (onemol->xflag == 0) + error->all(FLERR,"Fix deposit molecule must have coordinates"); + if (onemol->typeflag == 0) + error->all(FLERR,"Fix deposit molecule must have atom types"); + + onemol->compute_center(); + } + + // setup of coords and imagesflags array + + if (mode == ATOM) natom = 1; + else natom = onemol->natoms; + memory->create(coords,natom,3,"deposit:coords"); + memory->create(imageflags,natom,"deposit:imageflags"); + // setup scaling double xscale,yscale,zscale; @@ -129,16 +154,9 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : ty *= yscale; tz *= zscale; - // maxtag_all = current max tag for all atoms + // find current max atom and molecule IDs if necessary - if (idnext) { - int *tag = atom->tag; - int nlocal = atom->nlocal; - - int maxtag = 0; - for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]); - MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_INT,MPI_MAX,world); - } + if (idnext) find_maxid(); // random number generator, same for all procs @@ -158,6 +176,8 @@ FixDeposit::~FixDeposit() { delete random; delete [] idregion; + memory->destroy(coords); + memory->destroy(imageflags); } /* ---------------------------------------------------------------------- */ @@ -186,15 +206,19 @@ void FixDeposit::init() void FixDeposit::pre_exchange() { - int i,j; + int i,j,m; int flag,flagall; double coord[3],lamda[3],delx,dely,delz,rsq; - double *newcoord; + double alpha,beta,gamma; + double r[3],rotmat[3][3],quat[4]; + double *newcoord,*center; // just return if should not be called on this timestep if (next_reneighbor != update->ntimestep) return; + int dimension = domain->dimension; + // compute current offset = bottom of insertion volume double offset = 0.0; @@ -209,6 +233,10 @@ void FixDeposit::pre_exchange() subhi = domain->subhi_lamda; } + // find current max atom and molecule IDs if necessary + + if (!idnext) find_maxid(); + // attempt an insertion until successful int nfix = modify->nfix; @@ -219,7 +247,7 @@ void FixDeposit::pre_exchange() while (attempt < maxattempt) { attempt++; - // choose random position for new atom within region + // choose random position for new particle within region coord[0] = xlo + random->uniform() * (xhi-xlo); coord[1] = ylo + random->uniform() * (yhi-ylo); @@ -232,18 +260,19 @@ void FixDeposit::pre_exchange() // adjust vertical coord by offset - if (domain->dimension == 2) coord[1] += offset; + if (dimension == 2) coord[1] += offset; else coord[2] += offset; // if global, reset vertical coord to be lo-hi above highest atom // if local, reset vertical coord to be lo-hi above highest "nearby" atom // local computation computes lateral distance between 2 particles w/ PBC + // when done, have final coord of atom or center pt of molecule if (globalflag || localflag) { int dim; double max,maxall,delx,dely,delz,rsq; - if (domain->dimension == 2) { + if (dimension == 2) { dim = 1; max = domain->boxlo[1]; } else { @@ -259,7 +288,7 @@ void FixDeposit::pre_exchange() dely = coord[1] - x[i][1]; delz = 0.0; domain->minimum_image(delx,dely,delz); - if (domain->dimension == 2) rsq = delx*delx; + if (dimension == 2) rsq = delx*delx; else rsq = delx*delx + dely*dely; if (rsq > deltasq) continue; } @@ -267,38 +296,76 @@ void FixDeposit::pre_exchange() } MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world); - if (domain->dimension == 2) + if (dimension == 2) coord[1] = maxall + lo + random->uniform()*(hi-lo); else coord[2] = maxall + lo + random->uniform()*(hi-lo); } - // now have final coord - // if distance to any atom is less than near, try again + // coords = coords of all atoms + // for molecule, perform random rotation around center pt + // apply PBC so final coords are inside box + // also store image flag modified due to PBC + + if (mode == ATOM) { + coords[0][0] = coord[0]; + coords[0][1] = coord[1]; + coords[0][2] = coord[2]; + } else { + if (dimension == 3) { + r[0] = random->uniform() - 0.5; + r[1] = random->uniform() - 0.5; + r[2] = random->uniform() - 0.5; + } else { + r[0] = r[1] = 0.0; + r[2] = 1.0; + } + double theta = random->uniform() * MY_2PI; + MathExtra::norm3(r); + MathExtra::axisangle_to_quat(r,theta,quat); + MathExtra::quat_to_mat(quat,rotmat); + center = onemol->center; + for (i = 0; i < natom; i++) { + MathExtra::matvec(rotmat,onemol->dx[i],coords[i]); + coords[i][0] += center[0] + coord[0]; + coords[i][1] += center[1] + coord[1]; + coords[i][2] += center[2] + coord[2]; + + imageflags[i] = ((tagint) IMGMAX << IMG2BITS) | + ((tagint) IMGMAX << IMGBITS) | IMGMAX; + domain->remap(coords[i],imageflags[i]); + } + } + + // if distance to any inserted atom is less than near, try again double **x = atom->x; int nlocal = atom->nlocal; flag = 0; - for (i = 0; i < nlocal; i++) { - delx = coord[0] - x[i][0]; - dely = coord[1] - x[i][1]; - delz = coord[2] - x[i][2]; - domain->minimum_image(delx,dely,delz); - rsq = delx*delx + dely*dely + delz*delz; - if (rsq < nearsq) flag = 1; + for (m = 0; m < natom; m++) { + for (i = 0; i < nlocal; i++) { + delx = coords[m][0] - x[i][0]; + dely = coords[m][1] - x[i][1]; + delz = coords[m][2] - x[i][2]; + domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + if (rsq < nearsq) flag = 1; + } } MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world); if (flagall) continue; - // insertion will proceed - // choose random velocity for new atom + // proceed with insertion + // choose random velocity for new particle + // used for every atom in molecule double vxtmp = vxlo + random->uniform() * (vxhi-vxlo); double vytmp = vylo + random->uniform() * (vyhi-vylo); double vztmp = vzlo + random->uniform() * (vzhi-vzlo); - // if we have a sputter target change velocity vector accordingly + // if target specified, change velocity vector accordingly + if (targetflag) { double vel = sqrt(vxtmp*vxtmp + vytmp*vytmp + vztmp*vztmp); delx = tx - coord[0]; @@ -313,67 +380,84 @@ void FixDeposit::pre_exchange() } } - // check if new atom is in my sub-box or above it if I'm highest proc - // if so, add to my list via create_atom() - // initialize info about the atoms + // check if new atoms are in my sub-box or above it if I am highest proc + // if so, add atom to my list via create_atom() + // initialize additional info about the atoms // set group mask to "all" plus fix group - if (domain->triclinic) { - domain->x2lamda(coord,lamda); - newcoord = lamda; - } else newcoord = coord; + for (m = 0; m < natom; m++) { + if (domain->triclinic) { + domain->x2lamda(coords[m],lamda); + newcoord = lamda; + } else newcoord = coords[m]; + + flag = 0; + if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && + newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] && + newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1; + else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] && + comm->myloc[2] == comm->procgrid[2]-1 && + newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && + newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; + else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] && + comm->myloc[1] == comm->procgrid[1]-1 && + newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; - flag = 0; - if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && - newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] && - newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1; - else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] && - comm->myloc[2] == comm->procgrid[2]-1 && - newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && - newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; - else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] && - comm->myloc[1] == comm->procgrid[1]-1 && - newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; - - if (flag) { - atom->avec->create_atom(ntype,coord); - int m = atom->nlocal - 1; - atom->type[m] = ntype; - atom->mask[m] = 1 | groupbit; - atom->v[m][0] = vxtmp; - atom->v[m][1] = vytmp; - atom->v[m][2] = vztmp; - for (j = 0; j < nfix; j++) - if (fix[j]->create_attribute) fix[j]->set_arrays(m); + if (flag) { + if (mode == ATOM) atom->avec->create_atom(ntype,coords[m]); + else atom->avec->create_atom(onemol->type[m],coords[m]); + int n = atom->nlocal - 1; + atom->tag[n] = maxtag_all + m+1; + if (atom->molecule_flag) atom->molecule[n] = maxmol_all; + atom->mask[n] = 1 | groupbit; + atom->image[n] = imageflags[m]; + atom->v[n][0] = vxtmp; + atom->v[n][1] = vytmp; + atom->v[n][2] = vztmp; + if (mode == MOLECULE) atom->add_molecule_atom(onemol,m,n,maxtag_all); + for (j = 0; j < nfix; j++) + if (fix[j]->create_attribute) fix[j]->set_arrays(n); + } } - MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world); + + // old code: unsuccessful if no proc performed insertion of an atom + // don't think that check is necessary + // if get this far, should always be succesful + // would be hard to undo partial insertion for a molecule + // better to check how many atoms could be inserted (w/out inserting) + // then sum to insure all are inserted, before doing actual insertion + // MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world); + + success = 1; break; } - // warn if not successful b/c too many attempts or no proc owned particle + // warn if not successful b/c too many attempts if (!success && comm->me == 0) error->warning(FLERR,"Particle deposition was unsuccessful",0); - // reset global natoms - // if idnext, set new atom ID to incremented maxtag_all - // else set new atom ID to value beyond all current atoms + // reset global natoms,nbonds,etc + // increment maxtag_all and maxmol_all if necessary // if global map exists, reset it now instead of waiting for comm // since adding an atom messes up ghosts if (success) { - atom->natoms += 1; - if (atom->tag_enable) { - if (idnext) { - maxtag_all++; - if (atom->nlocal && atom->tag[atom->nlocal-1] == 0) - atom->tag[atom->nlocal-1] = maxtag_all; - } else atom->tag_extend(); - if (atom->map_style) { - atom->nghost = 0; - atom->map_init(); - atom->map_set(); - } + atom->natoms += natom; + if (mode == MOLECULE) { + atom->nbonds += onemol->nbonds; + atom->nangles += onemol->nangles; + atom->ndihedrals += onemol->ndihedrals; + atom->nimpropers += onemol->nimpropers; + } + if (idnext) { + maxtag_all += natom; + if (mode == MOLECULE) maxmol_all++; + } + if (atom->map_style) { + atom->nghost = 0; + atom->map_init(); + atom->map_set(); } } @@ -385,6 +469,28 @@ void FixDeposit::pre_exchange() else next_reneighbor = 0; } +/* ---------------------------------------------------------------------- + maxtag_all = current max atom ID for all atoms + maxmol_all = current max molecule ID for all atoms +------------------------------------------------------------------------- */ + +void FixDeposit::find_maxid() +{ + int *tag = atom->tag; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + + int max = 0; + for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); + MPI_Allreduce(&max,&maxtag_all,1,MPI_INT,MPI_MAX,world); + + if (mode == MOLECULE) { + max = 0; + for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]); + MPI_Allreduce(&max,&maxmol_all,1,MPI_INT,MPI_MAX,world); + } +} + /* ---------------------------------------------------------------------- parse optional parameters at end of input line ------------------------------------------------------------------------- */ @@ -404,6 +510,14 @@ void FixDeposit::options(int narg, char **arg) idregion = new char[n]; strcpy(idregion,arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"mol") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); + int imol = atom->find_molecule(arg[iarg+1]); + if (imol == -1) + error->all(FLERR,"Molecule ID for fix deposit does not exist"); + mode = MOLECULE; + onemol = atom->molecules[imol]; + iarg += 2; } else if (strcmp(arg[iarg],"id") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); if (strcmp(arg[iarg+1],"max") == 0) idnext = 0; @@ -423,11 +537,13 @@ void FixDeposit::options(int narg, char **arg) globalflag = 0; lo = force->numeric(FLERR,arg[iarg+1]); hi = force->numeric(FLERR,arg[iarg+2]); - deltasq = force->numeric(FLERR,arg[iarg+3])*force->numeric(FLERR,arg[iarg+3]); + deltasq = force->numeric(FLERR,arg[iarg+3]) * + force->numeric(FLERR,arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"near") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - nearsq = force->numeric(FLERR,arg[iarg+1])*force->numeric(FLERR,arg[iarg+1]); + nearsq = force->numeric(FLERR,arg[iarg+1]) * + force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"attempt") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); diff --git a/src/fix_deposit.h b/src/fix_deposit.h index 8f1acf7c9d..2d93effd03 100644 --- a/src/fix_deposit.h +++ b/src/fix_deposit.h @@ -40,15 +40,24 @@ class FixDeposit : public Fix { private: int ninsert,nfreq,seed; int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag,targetflag; - char *idregion; + int mode; double lo,hi,deltasq,nearsq,rate; double vxlo,vxhi,vylo,vyhi,vzlo,vzhi; double xlo,xhi,ylo,yhi,zlo,zhi; double tx,ty,tz; + char *idregion; + class Molecule *onemol; + + int natom; + double **coords; + int *imageflags; + int nfirst,ninserted; - int idnext,maxtag_all; + int idnext; + int maxtag_all,maxmol_all; class RanPark *random; + void find_maxid(); void options(int, char **); }; diff --git a/src/input.cpp b/src/input.cpp index 09cba19c59..5638cf089f 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -564,6 +564,7 @@ int Input::execute_command() else if (!strcmp(command,"mass")) mass(); else if (!strcmp(command,"min_modify")) min_modify(); else if (!strcmp(command,"min_style")) min_style(); + else if (!strcmp(command,"molecule")) molecule(); else if (!strcmp(command,"neigh_modify")) neigh_modify(); else if (!strcmp(command,"neighbor")) neighbor_command(); else if (!strcmp(command,"newton")) newton(); @@ -1311,6 +1312,13 @@ void Input::min_style() /* ---------------------------------------------------------------------- */ +void Input::molecule() +{ + atom->add_molecule(narg,arg); +} + +/* ---------------------------------------------------------------------- */ + void Input::neigh_modify() { neighbor->modify_params(narg,arg); diff --git a/src/input.h b/src/input.h index d0d3620058..f882310b19 100644 --- a/src/input.h +++ b/src/input.h @@ -103,6 +103,7 @@ class Input : protected Pointers { void mass(); void min_modify(); void min_style(); + void molecule(); void neigh_modify(); void neighbor_command(); void newton(); diff --git a/src/molecule.cpp b/src/molecule.cpp new file mode 100644 index 0000000000..69e810845d --- /dev/null +++ b/src/molecule.cpp @@ -0,0 +1,1060 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "string.h" +#include "molecule.h" +#include "atom.h" +#include "atom_vec.h" +#include "force.h" +#include "comm.h" +#include "domain.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define MAXLINE 256 + +/* ---------------------------------------------------------------------- */ + +Molecule::Molecule(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) +{ + me = comm->me; + + if (narg != 2) error->all(FLERR,"Illegal molecule command"); + + int n = strlen(arg[0]) + 1; + id = new char[n]; + strcpy(id,arg[0]); + + for (int i = 0; i < n-1; i++) + if (!isalnum(id[i]) && id[i] != '_') + error->all(FLERR, + "Molecule ID must be alphanumeric or underscore characters"); + + char *file = arg[1]; + + // initialize all fields to empty + + initialize(); + + // scan file for sizes of all fields and allocate them + + if (me == 0) open(file); + read(0); + if (me == 0) fclose(fp); + allocate(); + + // read file again to populate all fields + + if (me == 0) open(file); + read(1); + if (me == 0) fclose(fp); +} + +/* ---------------------------------------------------------------------- */ + +Molecule::~Molecule() +{ + delete [] id; + deallocate(); +} + +/* ---------------------------------------------------------------------- + compute center = geometric center of molecule + also compute dx = displacement of each atom from origin +------------------------------------------------------------------------- */ + +void Molecule::compute_center() +{ + if (centerflag) return; + centerflag = 1; + + center[0] = center[1] = center[2] = 0.0; + for (int i = 0; i < natoms; i++) { + center[0] += x[i][0]; + center[1] += x[i][1]; + center[2] += x[i][2]; + } + center[0] /= natoms; + center[1] /= natoms; + center[2] /= natoms; + + memory->destroy(dx); + memory->create(dx,natoms,3,"molecule:dx"); + + for (int i = 0; i < natoms; i++) { + dx[i][0] = x[i][0] - center[0]; + dx[i][1] = x[i][1] - center[1]; + dx[i][2] = x[i][2] - center[2]; + } +} + +/* ---------------------------------------------------------------------- + compute xcm = center or mass of molecule + account for finite size particles + also compute dxcom = displacement of each atom from COM +------------------------------------------------------------------------- */ + +void Molecule::compute_xcm() +{ + comflag = 1; + + // NOTE this is not yet right + + xcm[0] = xcm[1] = xcm[2] = 0.0; + for (int i = 0; i < natoms; i++) { + xcm[0] += x[i][0]; + xcm[1] += x[i][1]; + xcm[2] += x[i][2]; + } + xcm[0] /= natoms; + xcm[1] /= natoms; + xcm[2] /= natoms; + + memory->destroy(dxcom); + memory->create(dxcom,natoms,3,"molecule:dxcom"); + + for (int i = 0; i < natoms; i++) { + dxcom[i][0] = x[i][0] - center[0]; + dxcom[i][1] = x[i][1] - center[1]; + dxcom[i][2] = x[i][2] - center[2]; + } +} + +/* ---------------------------------------------------------------------- + read molecule info from file + flag = 0, just scan for sizes of fields + flag = 1, read and store fields +------------------------------------------------------------------------- */ + +void Molecule::read(int flag) +{ + char line[MAXLINE],keyword[MAXLINE]; + char *eof,*ptr; + + // skip 1st line of file + + if (me == 0) { + eof = fgets(line,MAXLINE,fp); + if (eof == NULL) error->one(FLERR,"Unexpected end of molecule file"); + } + + // read header lines + // skip blank lines or lines that start with "#" + // stop when read an unrecognized line + + while (1) { + + readline(line); + + // trim anything from '#' onward + // if line is blank, continue + + if (ptr = strchr(line,'#')) *ptr = '\0'; + if (strspn(line," \t\n\r") == strlen(line)) continue; + + // search line for header keywords and set corresponding variable + + if (strstr(line,"atoms")) sscanf(line,"%d",&natoms); + else if (strstr(line,"bonds")) sscanf(line,"%d",&nbonds); + else if (strstr(line,"angles")) sscanf(line,"%d",&nangles); + else if (strstr(line,"dihedrals")) sscanf(line,"%d",&ndihedrals); + else if (strstr(line,"impropers")) sscanf(line,"%d",&nimpropers); + else break; + } + + // error check + + if (flag == 0) { + if (natoms == 0) error->all(FLERR,"No atom count in molecule file"); + + if (nbonds && !atom->avec->bonds_allow) + error->all(FLERR,"Bonds in molecule file not supported by atom style"); + if (nangles && !atom->avec->angles_allow) + error->all(FLERR,"Angles in molecule file not supported by atom style"); + if (ndihedrals && !atom->avec->dihedrals_allow) + error->all(FLERR, + "Dihedrals in molecule file not supported by atom style"); + if (nimpropers && !atom->avec->impropers_allow) + error->all(FLERR, + "Impropers in molecule file not supported by atom style"); + } + + // count = vector for tallying bonds,angles,etc per atom + + if (flag == 0) memory->create(count,natoms,"molecule:count"); + else count = NULL; + + // grab keyword and skip next line + + parse_keyword(0,line,keyword); + readline(line); + + // loop over sections of molecule file + + while (strlen(keyword)) { + if (strcmp(keyword,"Coords") == 0) { + xflag = 1; + if (flag) coords(line); + else skip_lines(natoms,line); + } else if (strcmp(keyword,"Types") == 0) { + typeflag = 1; + if (flag) types(line); + else skip_lines(natoms,line); + } else if (strcmp(keyword,"Charges") == 0) { + qflag = 1; + if (flag) charges(line); + else skip_lines(natoms,line); + } else if (strcmp(keyword,"Diameters") == 0) { + radiusflag = 1; + if (flag) diameters(line); + else skip_lines(natoms,line); + } else if (strcmp(keyword,"Masses") == 0) { + rmassflag = 1; + if (flag) masses(line); + else skip_lines(natoms,line); + + } else if (strcmp(keyword,"Bonds") == 0) { + if (nbonds == 0) + error->all(FLERR,"Molecule file has bonds but no nbonds setting"); + bondflag = 1; + bonds(flag,line); + } else if (strcmp(keyword,"Angles") == 0) { + if (nangles == 0) + error->all(FLERR,"Molecule file has angles but no nangles setting"); + angleflag = 1; + angles(flag,line); + } else if (strcmp(keyword,"Dihedrals") == 0) { + if (ndihedrals == 0) error->all(FLERR,"Molecule file has dihedrals " + "but no ndihedrals setting"); + dihedralflag = 1; + dihedrals(flag,line); + } else if (strcmp(keyword,"Impropers") == 0) { + if (nimpropers == 0) error->all(FLERR,"Molecule file has impropers " + "but no nimpropers setting"); + improperflag = 1; + impropers(flag,line); + + } else if (strcmp(keyword,"Special Bond Counts") == 0) { + nspecialflag = 1; + nspecial_read(flag,line); + } else if (strcmp(keyword,"Special Bonds") == 0) { + specialflag = 1; + if (flag) special_read(line); + else skip_lines(natoms,line); + + } else error->one(FLERR,"Unknown section in molecule file"); + + parse_keyword(1,line,keyword); + } + + // clean up + + memory->destroy(count); + + // error check + + if (flag == 0) { + if (qflag && !atom->q_flag) + error->all(FLERR,"Molecule file has undefined atom property"); + if (radiusflag && !atom->radius_flag) + error->all(FLERR,"Molecule file has undefined atom property"); + if (rmassflag && !atom->rmass_flag) + error->all(FLERR,"Molecule file has undefined atom property"); + + if (bondflag && !atom->avec->bonds_allow) + error->all(FLERR,"Invalid molecule file section: Bonds"); + if (angleflag && !atom->avec->angles_allow) + error->all(FLERR,"Invalid molecule file section: Angles"); + if (dihedralflag && !atom->avec->dihedrals_allow) + error->all(FLERR,"Invalid molecule file section: Dihedrals"); + if (improperflag && !atom->avec->impropers_allow) + error->all(FLERR,"Invalid molecule file section: Impropers"); + + if (bond_per_atom > atom->bond_per_atom || + angle_per_atom > atom->angle_per_atom || + dihedral_per_atom > atom->dihedral_per_atom || + improper_per_atom > atom->improper_per_atom) + error->all(FLERR,"Molecule file bond/angle/etc counts " + "per atom are too large"); + + if ((nspecialflag && !specialflag) || (!nspecialflag && specialflag)) + error->all(FLERR,"Molecule file needs both Special Bond sections"); + if (specialflag && !bondflag) + error->all(FLERR,"Molecule file has special flags but no bonds"); + if (maxspecial > atom->maxspecial) + error->all(FLERR,"Molecule file special bond counts are too large"); + } +} + +/* ---------------------------------------------------------------------- + read coords from file +------------------------------------------------------------------------- */ + +void Molecule::coords(char *line) +{ + int tmp; + for (int i = 0; i < natoms; i++) { + readline(line); + sscanf(line,"%d %lg %lg %lg",&tmp,&x[i][0],&x[i][1],&x[i][2]); + } + + if (domain->dimension == 2) { + for (int i = 0; i < natoms; i++) + if (x[i][2] != 0.0) + error->all(FLERR,"Molecule file z coord must be 0.0 for 2d"); + } +} + +/* ---------------------------------------------------------------------- + read types from file +------------------------------------------------------------------------- */ + +void Molecule::types(char *line) +{ + int tmp; + for (int i = 0; i < natoms; i++) { + readline(line); + sscanf(line,"%d %d",&tmp,&type[i]); + } +} + +/* ---------------------------------------------------------------------- + read charges from file +------------------------------------------------------------------------- */ + +void Molecule::charges(char *line) +{ + int tmp; + for (int i = 0; i < natoms; i++) { + readline(line); + sscanf(line,"%d %lg",&tmp,&q[i]); + } +} + +/* ---------------------------------------------------------------------- + read diameters from file and set radii +------------------------------------------------------------------------- */ + +void Molecule::diameters(char *line) +{ + int tmp; + for (int i = 0; i < natoms; i++) { + readline(line); + sscanf(line,"%d %lg",&tmp,&radius[i]); + radius[i] *= 0.5; + } +} + +/* ---------------------------------------------------------------------- + read masses from file +------------------------------------------------------------------------- */ + +void Molecule::masses(char *line) +{ + int tmp; + for (int i = 0; i < natoms; i++) { + readline(line); + sscanf(line,"%d %lg",&tmp,&rmass[i]); + } +} + +/* ---------------------------------------------------------------------- + read bonds from file + store each with both atoms if newton_bond = 0 + if flag = 0, just count bonds/atom + if flag = 1, store them with atoms +------------------------------------------------------------------------- */ + +void Molecule::bonds(int flag, char *line) +{ + int m,tmp,itype,atom1,atom2; + int newton_bond = force->newton_bond; + + if (flag == 0) + for (int i = 0; i < natoms; i++) count[i] = 0; + else + for (int i = 0; i < natoms; i++) num_bond[i] = 0; + + for (int i = 0; i < nbonds; i++) { + readline(line); + sscanf(line,"%d %d %d %d",&tmp,&itype,&atom1,&atom2); + + if (atom1 <= 0 || atom1 > natoms || + atom2 <= 0 || atom2 > natoms) + error->one(FLERR,"Invalid atom ID in Bonds section of molecule file"); + if (itype <= 0 || itype > atom->nbondtypes) + error->one(FLERR,"Invalid bond type in Bonds section of molecule file"); + + if (flag) { + m = atom1-1; + bond_type[m][num_bond[m]] = itype; + bond_atom[m][num_bond[m]] = atom2; + num_bond[m]++; + if (newton_bond == 0) { + m = atom2-1; + bond_type[m][num_bond[m]] = itype; + bond_atom[m][num_bond[m]] = atom1; + num_bond[m]++; + } + } else { + count[atom1-1]++; + if (newton_bond == 0) count[atom2-1]++; + } + } + + // bond_per_atom = max of count vector + + if (flag == 0) { + bond_per_atom = 0; + for (int i = 0; i < natoms; i++) + bond_per_atom = MAX(bond_per_atom,count[i]); + } +} + +/* ---------------------------------------------------------------------- + read angles from file + store each with all 3 atoms if newton_bond = 0 + if flag = 0, just count angles/atom + if flag = 1, store them with atoms +------------------------------------------------------------------------- */ + +void Molecule::angles(int flag, char *line) +{ + int m,tmp,itype,atom1,atom2,atom3; + int newton_bond = force->newton_bond; + + if (flag == 0) + for (int i = 0; i < natoms; i++) count[i] = 0; + else + for (int i = 0; i < natoms; i++) num_angle[i] = 0; + + for (int i = 0; i < nangles; i++) { + readline(line); + sscanf(line,"%d %d %d %d %d",&tmp,&itype,&atom1,&atom2,&atom3); + + if (atom1 <= 0 || atom1 > natoms || + atom2 <= 0 || atom2 > natoms || + atom3 <= 0 || atom3 > natoms) + error->one(FLERR,"Invalid atom ID in Angles section of molecule file"); + if (itype <= 0 || itype > atom->nangletypes) + error->one(FLERR,"Invalid angle type in Angles section of molecule file"); + + if (flag) { + m = atom2-1; + angle_type[m][num_angle[m]] = itype; + angle_atom1[m][num_angle[m]] = atom1; + angle_atom2[m][num_angle[m]] = atom2; + angle_atom3[m][num_angle[m]] = atom3; + num_angle[m]++; + if (newton_bond == 0) { + m = atom1-1; + angle_type[m][num_angle[m]] = itype; + angle_atom1[m][num_angle[m]] = atom1; + angle_atom2[m][num_angle[m]] = atom2; + angle_atom3[m][num_angle[m]] = atom3; + num_angle[m]++; + m = atom3-1; + angle_type[m][num_angle[m]] = itype; + angle_atom1[m][num_angle[m]] = atom1; + angle_atom2[m][num_angle[m]] = atom2; + angle_atom3[m][num_angle[m]] = atom3; + num_angle[m]++; + } + } else { + count[atom2-1]++; + if (newton_bond == 0) { + count[atom1-1]++; + count[atom3-1]++; + } + } + } + + // angle_per_atom = max of count vector + + if (flag == 0) { + angle_per_atom = 0; + for (int i = 0; i < natoms; i++) + angle_per_atom = MAX(angle_per_atom,count[i]); + } +} + +/* ---------------------------------------------------------------------- + read dihedrals from file + store each with all 4 atoms if newton_bond = 0 + if flag = 0, just count dihedrals/atom + if flag = 1, store them with atoms +------------------------------------------------------------------------- */ + +void Molecule::dihedrals(int flag, char *line) +{ + int m,tmp,itype,atom1,atom2,atom3,atom4; + int newton_bond = force->newton_bond; + + if (flag == 0) + for (int i = 0; i < natoms; i++) count[i] = 0; + else + for (int i = 0; i < natoms; i++) num_dihedral[i] = 0; + + for (int i = 0; i < ndihedrals; i++) { + readline(line); + sscanf(line,"%d %d %d %d %d %d",&tmp,&itype,&atom1,&atom2,&atom3,&atom4); + + if (atom1 <= 0 || atom1 > natoms || + atom2 <= 0 || atom2 > natoms || + atom3 <= 0 || atom3 > natoms || + atom4 <= 0 || atom4 > natoms) + error->one(FLERR, + "Invalid atom ID in dihedrals section of molecule file"); + if (itype <= 0 || itype > atom->ndihedraltypes) + error->one(FLERR, + "Invalid dihedral type in dihedrals section of molecule file"); + + if (flag) { + m = atom2-1; + dihedral_type[m][num_dihedral[m]] = itype; + dihedral_atom1[m][num_dihedral[m]] = atom1; + dihedral_atom2[m][num_dihedral[m]] = atom2; + dihedral_atom3[m][num_dihedral[m]] = atom3; + num_dihedral[m]++; + if (newton_bond == 0) { + m = atom1-1; + dihedral_type[m][num_dihedral[m]] = itype; + dihedral_atom1[m][num_dihedral[m]] = atom1; + dihedral_atom2[m][num_dihedral[m]] = atom2; + dihedral_atom3[m][num_dihedral[m]] = atom3; + num_dihedral[m]++; + m = atom3-1; + dihedral_type[m][num_dihedral[m]] = itype; + dihedral_atom1[m][num_dihedral[m]] = atom1; + dihedral_atom2[m][num_dihedral[m]] = atom2; + dihedral_atom3[m][num_dihedral[m]] = atom3; + num_dihedral[m]++; + m = atom4-1; + dihedral_type[m][num_dihedral[m]] = itype; + dihedral_atom1[m][num_dihedral[m]] = atom1; + dihedral_atom2[m][num_dihedral[m]] = atom2; + dihedral_atom3[m][num_dihedral[m]] = atom3; + num_dihedral[m]++; + } + } else { + count[atom2-1]++; + if (newton_bond == 0) { + count[atom1-1]++; + count[atom3-1]++; + count[atom4-1]++; + } + } + } + + // dihedral_per_atom = max of count vector + + if (flag == 0) { + dihedral_per_atom = 0; + for (int i = 0; i < natoms; i++) + dihedral_per_atom = MAX(dihedral_per_atom,count[i]); + } +} + +/* ---------------------------------------------------------------------- + read impropers from file + store each with all 4 atoms if newton_bond = 0 + if flag = 0, just count impropers/atom + if flag = 1, store them with atoms +------------------------------------------------------------------------- */ + +void Molecule::impropers(int flag, char *line) +{ + int m,tmp,itype,atom1,atom2,atom3,atom4; + int newton_bond = force->newton_bond; + + if (flag == 0) + for (int i = 0; i < natoms; i++) count[i] = 0; + else + for (int i = 0; i < natoms; i++) num_improper[i] = 0; + + for (int i = 0; i < nimpropers; i++) { + readline(line); + sscanf(line,"%d %d %d %d %d %d",&tmp,&itype,&atom1,&atom2,&atom3,&atom4); + + if (atom1 <= 0 || atom1 > natoms || + atom2 <= 0 || atom2 > natoms || + atom3 <= 0 || atom3 > natoms || + atom4 <= 0 || atom4 > natoms) + error->one(FLERR, + "Invalid atom ID in impropers section of molecule file"); + if (itype <= 0 || itype > atom->nimpropertypes) + error->one(FLERR, + "Invalid improper type in impropers section of molecule file"); + + if (flag) { + m = atom2-1; + improper_type[m][num_improper[m]] = itype; + improper_atom1[m][num_improper[m]] = atom1; + improper_atom2[m][num_improper[m]] = atom2; + improper_atom3[m][num_improper[m]] = atom3; + num_improper[m]++; + if (newton_bond == 0) { + m = atom1-1; + improper_type[m][num_improper[m]] = itype; + improper_atom1[m][num_improper[m]] = atom1; + improper_atom2[m][num_improper[m]] = atom2; + improper_atom3[m][num_improper[m]] = atom3; + num_improper[m]++; + m = atom3-1; + improper_type[m][num_improper[m]] = itype; + improper_atom1[m][num_improper[m]] = atom1; + improper_atom2[m][num_improper[m]] = atom2; + improper_atom3[m][num_improper[m]] = atom3; + num_improper[m]++; + m = atom4-1; + improper_type[m][num_improper[m]] = itype; + improper_atom1[m][num_improper[m]] = atom1; + improper_atom2[m][num_improper[m]] = atom2; + improper_atom3[m][num_improper[m]] = atom3; + num_improper[m]++; + } + } else { + count[atom2-1]++; + if (newton_bond == 0) { + count[atom1-1]++; + count[atom3-1]++; + count[atom4-1]++; + } + } + } + + // improper_per_atom = max of count vector + + if (flag == 0) { + improper_per_atom = 0; + for (int i = 0; i < natoms; i++) + improper_per_atom = MAX(improper_per_atom,count[i]); + } +} + +/* ---------------------------------------------------------------------- + read 3 special bonds counts from file + if flag = 0, just tally maxspecial + if flag = 1, store them with atoms +------------------------------------------------------------------------- */ + +void Molecule::nspecial_read(int flag, char *line) +{ + int tmp,c1,c2,c3; + + if (flag == 0) maxspecial = 0; + + for (int i = 0; i < natoms; i++) { + readline(line); + sscanf(line,"%d %d %d %d",&tmp,&c1,&c2,&c3); + + if (flag) { + nspecial[i][0] = c1; + nspecial[i][1] = c1+c2; + nspecial[i][2] = c1+c2+c3; + } else maxspecial = MAX(maxspecial,c1+c2+c3); + } +} + +/* ---------------------------------------------------------------------- + read special bond indices from file +------------------------------------------------------------------------- */ + +void Molecule::special_read(char *line) +{ + int m,nwords; + char **words = new char*[maxspecial+1]; + + for (int i = 0; i < natoms; i++) { + readline(line); + nwords = parse(line,words,maxspecial+1); + if (nwords != nspecial[i][2]+1) + error->all(FLERR,"Molecule file special list " + "does not match special count"); + + for (m = 1; m < nwords; m++) { + special[i][m-1] = atoi(words[m]); + if (special[i][m-1] <= 0 || special[i][m-1] > natoms || + special[i][m-1] == i+1) + error->all(FLERR,"Invalid special atom index in molecule file"); + } + } + + delete [] words; +} + +/* ---------------------------------------------------------------------- + init all data structures to empty +------------------------------------------------------------------------- */ + +void Molecule::initialize() +{ + natoms = 0; + nbonds = nangles = ndihedrals = nimpropers = 0; + bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0; + maxspecial = 0; + + xflag = typeflag = qflag = radiusflag = rmassflag = 0; + bondflag = angleflag = dihedralflag = improperflag = 0; + nspecialflag = specialflag = 0; + + centerflag = comflag = inertiaflag = 0; + + x = NULL; + type = NULL; + q = NULL; + radius = NULL; + rmass = NULL; + + num_bond = NULL; + bond_type = NULL; + bond_atom = NULL; + + num_angle = NULL; + angle_type = NULL; + angle_atom1 = angle_atom2 = angle_atom3 = NULL; + + num_dihedral = NULL; + dihedral_type = NULL; + dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = NULL; + + num_improper = NULL; + improper_type = NULL; + improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = NULL; + + nspecial = NULL; + special = NULL; + + dx = NULL; + dxcom = NULL; +} + +/* ---------------------------------------------------------------------- + allocate all data structures +------------------------------------------------------------------------- */ + +void Molecule::allocate() +{ + if (xflag) memory->create(x,natoms,3,"molecule:x"); + if (typeflag) memory->create(type,natoms,"molecule:type"); + if (qflag) memory->create(q,natoms,"molecule:q"); + if (radiusflag) memory->create(radius,natoms,"molecule:radius"); + if (rmassflag) memory->create(rmass,natoms,"molecule:rmass"); + + if (bondflag) { + memory->create(num_bond,natoms,"molecule:num_bond"); + memory->create(bond_type,natoms,bond_per_atom, + "molecule:bond_type"); + memory->create(bond_atom,natoms,bond_per_atom, + "molecule:bond_atom"); + } + + if (angleflag) { + memory->create(num_angle,natoms,"molecule:num_angle"); + memory->create(angle_type,natoms,angle_per_atom, + "molecule:angle_type"); + memory->create(angle_atom1,natoms,angle_per_atom, + "molecule:angle_atom1"); + memory->create(angle_atom2,natoms,angle_per_atom, + "molecule:angle_atom2"); + memory->create(angle_atom3,natoms,angle_per_atom, + "molecule:angle_atom3"); + } + + if (dihedralflag) { + memory->create(num_dihedral,natoms,"molecule:num_dihedral"); + memory->create(dihedral_type,natoms,dihedral_per_atom, + "molecule:dihedral_type"); + memory->create(dihedral_atom1,natoms,dihedral_per_atom, + "molecule:dihedral_atom1"); + memory->create(dihedral_atom2,natoms,dihedral_per_atom, + "molecule:dihedral_atom2"); + memory->create(dihedral_atom3,natoms,dihedral_per_atom, + "molecule:dihedral_atom3"); + memory->create(dihedral_atom4,natoms,dihedral_per_atom, + "molecule:dihedral_atom4"); + } + + if (improperflag) { + memory->create(num_improper,natoms,"molecule:num_improper"); + memory->create(improper_type,natoms,improper_per_atom, + "molecule:improper_type"); + memory->create(improper_atom1,natoms,improper_per_atom, + "molecule:improper_atom1"); + memory->create(improper_atom2,natoms,improper_per_atom, + "molecule:improper_atom2"); + memory->create(improper_atom3,natoms,improper_per_atom, + "molecule:improper_atom3"); + memory->create(improper_atom4,natoms,improper_per_atom, + "molecule:improper_atom4"); + } + + if (nspecialflag) + memory->create(nspecial,natoms,3,"molecule:nspecial"); + if (specialflag) + memory->create(special,natoms,maxspecial,"molecule:special"); +} + +/* ---------------------------------------------------------------------- + deallocate all data structures +------------------------------------------------------------------------- */ + +void Molecule::deallocate() +{ + memory->destroy(x); + memory->destroy(type); + memory->destroy(q); + memory->destroy(radius); + memory->destroy(rmass); + + memory->destroy(num_bond); + memory->destroy(bond_type); + memory->destroy(bond_atom); + + memory->destroy(num_angle); + memory->destroy(angle_type); + memory->destroy(angle_atom1); + memory->destroy(angle_atom2); + memory->destroy(angle_atom3); + + memory->destroy(num_dihedral); + memory->destroy(dihedral_type); + memory->destroy(dihedral_atom1); + memory->destroy(dihedral_atom2); + memory->destroy(dihedral_atom3); + memory->destroy(dihedral_atom4); + + memory->destroy(num_improper); + memory->destroy(improper_type); + memory->destroy(improper_atom1); + memory->destroy(improper_atom2); + memory->destroy(improper_atom3); + memory->destroy(improper_atom4); + + memory->destroy(nspecial); + memory->destroy(special); + + memory->destroy(dx); + memory->destroy(dxcom); +} + +/* ---------------------------------------------------------------------- + open molecule file +------------------------------------------------------------------------- */ + +void Molecule::open(char *file) +{ + fp = fopen(file,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open molecule file %s",file); + error->one(FLERR,str); + } +} + +/* ---------------------------------------------------------------------- + read and bcast a line +------------------------------------------------------------------------- */ + +void Molecule::readline(char *line) +{ + int n; + if (me == 0) { + if (fgets(line,MAXLINE,fp) == NULL) n = 0; + else n = strlen(line) + 1; + } + MPI_Bcast(&n,1,MPI_INT,0,world); + if (n == 0) error->all(FLERR,"Unexpected end of molecule file"); + MPI_Bcast(line,n,MPI_CHAR,0,world); +} + +/* ---------------------------------------------------------------------- + extract keyword from line + flag = 0, read and bcast line + flag = 1, line has already been read +------------------------------------------------------------------------- */ + +void Molecule::parse_keyword(int flag, char *line, char *keyword) +{ + if (flag) { + + // read upto non-blank line plus 1 following line + // eof is set to 1 if any read hits end-of-file + + int eof = 0; + if (me == 0) { + if (fgets(line,MAXLINE,fp) == NULL) eof = 1; + while (eof == 0 && strspn(line," \t\n\r") == strlen(line)) { + if (fgets(line,MAXLINE,fp) == NULL) eof = 1; + } + if (fgets(keyword,MAXLINE,fp) == NULL) eof = 1; + } + + // if eof, set keyword empty and return + + MPI_Bcast(&eof,1,MPI_INT,0,world); + if (eof) { + keyword[0] = '\0'; + return; + } + + // bcast keyword line to all procs + + int n; + if (me == 0) n = strlen(line) + 1; + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + } + + // copy non-whitespace portion of line into keyword + + int start = strspn(line," \t\n\r"); + int stop = strlen(line) - 1; + while (line[stop] == ' ' || line[stop] == '\t' + || line[stop] == '\n' || line[stop] == '\r') stop--; + line[stop+1] = '\0'; + strcpy(keyword,&line[start]); +} + +/* ---------------------------------------------------------------------- + skip N lines of file +------------------------------------------------------------------------- */ + +void Molecule::skip_lines(int n, char *line) +{ + for (int i = 0; i < n; i++) readline(line); +} + +/* ---------------------------------------------------------------------- + parse line into words separated by whitespace + return # of words + max = max pointers storable in words +------------------------------------------------------------------------- */ + +int Molecule::parse(char *line, char **words, int max) +{ + char *ptr; + + int nwords = 0; + words[nwords++] = strtok(line," \t\n\r\f"); + + while (ptr = strtok(NULL," \t\n\r\f")) { + if (nwords < max) words[nwords] = ptr; + nwords++; + } + + return nwords; +} + +/* ---------------------------------------------------------------------- + proc 0 prints molecule params +------------------------------------------------------------------------- */ + +/* + +void Molecule::print() +{ + printf("MOLECULE %s\n",id); + printf(" %d natoms\n",natoms); + if (nbonds) printf(" %d nbonds\n",nbonds); + if (nangles) printf(" %d nangles\n",nangles); + if (ndihedrals) printf(" %d ndihedrals\n",ndihedrals); + if (nimpropers) printf(" %d nimpropers\n",nimpropers); + + if (xflag) { + printf( "Coords:\n"); + for (int i = 0; i < natoms; i++) + printf(" %d %g %g %g\n",i+1,x[i][0],x[i][1],x[i][2]); + } + if (typeflag) { + printf( "Types:\n"); + for (int i = 0; i < natoms; i++) + printf(" %d %d\n",i+1,type[i]); + } + if (qflag) { + printf( "Charges:\n"); + for (int i = 0; i < natoms; i++) + printf(" %d %g\n",i+1,q[i]); + } + if (radiusflag) { + printf( "Radii:\n"); + for (int i = 0; i < natoms; i++) + printf(" %d %g\n",i+1,radius[i]); + } + if (rmassflag) { + printf( "Masses:\n"); + for (int i = 0; i < natoms; i++) + printf(" %d %g\n",i+1,rmass[i]); + } + + if (bondflag) { + printf( "Bonds:\n"); + for (int i = 0; i < natoms; i++) { + printf(" %d %d\n",i+1,num_bond[i]); + for (int j = 0; j < num_bond[i]; j++) + printf(" %d %d %d %d\n",j+1,bond_type[i][j],i+1,bond_atom[i][j]); + } + } + if (angleflag) { + printf( "Angles:\n"); + for (int i = 0; i < natoms; i++) { + printf(" %d %d\n",i+1,num_angle[i]); + for (int j = 0; j < num_angle[i]; j++) + printf(" %d %d %d %d %d\n", + j+1,angle_type[i][j], + angle_atom1[i][j],angle_atom2[i][j],angle_atom3[i][j]); + } + } + if (dihedralflag) { + printf( "Dihedrals:\n"); + for (int i = 0; i < natoms; i++) { + printf(" %d %d\n",i+1,num_dihedral[i]); + for (int j = 0; j < num_dihedral[i]; j++) + printf(" %d %d %d %d %d %d\n", + j+1,dihedral_type[i][j], + dihedral_atom1[i][j],dihedral_atom2[i][j], + dihedral_atom3[i][j],dihedral_atom4[i][j]); + } + } + if (improperflag) { + printf( "Impropers:\n"); + for (int i = 0; i < natoms; i++) { + printf(" %d %d\n",i+1,num_improper[i]); + for (int j = 0; j < num_improper[i]; j++) + printf(" %d %d %d %d %d %d\n", + j+1,improper_type[i][j], + improper_atom1[i][j],improper_atom2[i][j], + improper_atom3[i][j],improper_atom4[i][j]); + } + } + + if (specialflag) { + printf( "Special neighs:\n"); + for (int i = 0; i < natoms; i++) { + printf(" %d %d %d %d\n",i+1, + nspecial[i][0],nspecial[i][1]-nspecial[i][0], + nspecial[i][2]-nspecial[i][1]); + printf(" "); + for (int j = 0; j < nspecial[i][2]; j++) + printf(" %d",special[i][j]); + printf("\n"); + } + } +} + +*/ diff --git a/src/molecule.h b/src/molecule.h new file mode 100644 index 0000000000..1132a1c207 --- /dev/null +++ b/src/molecule.h @@ -0,0 +1,118 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_ONE_MOLECULE_H +#define LMP_ONE_MOLEUCULE_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class Molecule : protected Pointers { + public: + char *id; + + // number of atoms,bonds,etc in molecule + + int natoms; + int nbonds,nangles,ndihedrals,nimpropers; + + // max bond,angle,etc per atom + + int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom; + int maxspecial; + + // 1 if attribute defined in file, 0 if not + + int xflag,typeflag,qflag,radiusflag,rmassflag; + int bondflag,angleflag,dihedralflag,improperflag; + int nspecialflag,specialflag; + + // 1 if attribute defined or computed, 0 if not + + int centerflag,comflag,inertiaflag; + + // attributes + + double **x; // displacement of each atom from origin + int *type; // type of each atom + double *q; // charge on each atom + double *radius; // radius of each atom + double *rmass; // mass of each atom + + int *num_bond; // bonds, angles, dihedrals, impropers for each atom + int **bond_type; + int **bond_atom; + + int *num_angle; + int **angle_type; + int **angle_atom1,**angle_atom2,**angle_atom3; + + int *num_dihedral; + int **dihedral_type; + int **dihedral_atom1,**dihedral_atom2,**dihedral_atom3,**dihedral_atom4; + + int *num_improper; + int **improper_type; + int **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4; + + int **nspecial; + int **special; + + double center[3]; // geometric center of molecule + double xcm[3]; // center of mass of molecule + double inertia[6]; // moments of inertia of molecule + + double **dx; // displacement of each atom relative to center + double **dxcom; // displacement of each atom relative to COM + + Molecule(class LAMMPS *, int, char **); + ~Molecule(); + void compute_center(); + void compute_xcm(); + void compute_inertia(); + + private: + int me; + FILE *fp; + int *count; + + void read(int); + void coords(char *); + void types(char *); + void charges(char *); + void diameters(char *); + void masses(char *); + void bonds(int, char *); + void angles(int, char *); + void dihedrals(int, char *); + void impropers(int, char *); + void nspecial_read(int, char *); + void special_read(char *); + + void initialize(); + void allocate(); + void deallocate(); + + void open(char *); + void readline(char *); + void parse_keyword(int, char *, char *); + void skip_lines(int, char *); + int parse(char *, char **, int); + + // void print(); +}; + +} + +#endif