diff --git a/src/atom.cpp b/src/atom.cpp index df4db0a84..e46b1a724 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1,2255 +1,2256 @@ /* ---------------------------------------------------------------------- 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 #include #include #include #include #include #include "atom.h" #include "style_atom.h" #include "atom_vec.h" #include "atom_vec_ellipsoid.h" #include "comm.h" #include "neighbor.h" #include "force.h" #include "modify.h" #include "fix.h" #include "compute.h" #include "output.h" #include "thermo.h" #include "update.h" #include "domain.h" #include "group.h" #include "input.h" #include "variable.h" #include "molecule.h" #include "atom_masks.h" #include "math_const.h" #include "memory.h" #include "error.h" #ifdef LMP_USER_INTEL #include "neigh_request.h" #endif using namespace LAMMPS_NS; using namespace MathConst; #define DELTA 1 #define DELTA_MEMSTR 1024 #define EPSILON 1.0e-6 enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files /* ---------------------------------------------------------------------- */ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) { natoms = 0; nlocal = nghost = nmax = 0; ntypes = 0; nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; nbonds = nangles = ndihedrals = nimpropers = 0; firstgroupname = NULL; sortfreq = 1000; nextsort = 0; userbinsize = 0.0; maxbin = maxnext = 0; binhead = NULL; next = permute = NULL; // initialize atom arrays // customize by adding new array tag = NULL; type = mask = NULL; image = NULL; x = v = f = NULL; molecule = NULL; molindex = molatom = NULL; q = NULL; mu = NULL; omega = angmom = torque = NULL; radius = rmass = NULL; ellipsoid = line = tri = body = NULL; vfrac = s0 = NULL; x0 = NULL; spin = NULL; eradius = ervel = erforce = NULL; cs = csforce = vforce = ervelforce = NULL; etag = NULL; rho = drho = e = de = cv = NULL; vest = NULL; // USER-DPD uCond = uMech = uChem = uCG = uCGnew = NULL; duChem = NULL; dpdTheta = NULL; ssaAIR = NULL; // USER-SMD contact_radius = NULL; smd_data_9 = NULL; smd_stress = NULL; eff_plastic_strain = NULL; eff_plastic_strain_rate = NULL; damage = NULL; // molecular info bond_per_atom = extra_bond_per_atom = 0; num_bond = NULL; bond_type = NULL; bond_atom = NULL; angle_per_atom = extra_angle_per_atom = 0; num_angle = NULL; angle_type = NULL; angle_atom1 = angle_atom2 = angle_atom3 = NULL; dihedral_per_atom = extra_dihedral_per_atom = 0; num_dihedral = NULL; dihedral_type = NULL; dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = NULL; improper_per_atom = extra_improper_per_atom = 0; num_improper = NULL; improper_type = NULL; improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = NULL; maxspecial = 1; nspecial = NULL; special = NULL; // user-defined molecules nmolecule = 0; molecules = NULL; // custom atom arrays nivector = ndvector = 0; ivector = NULL; dvector = NULL; iname = dname = NULL; // initialize atom style and array existence flags // customize by adding new flag sphere_flag = peri_flag = electron_flag = 0; wavepacket_flag = sph_flag = 0; molecule_flag = 0; q_flag = mu_flag = 0; omega_flag = torque_flag = angmom_flag = 0; radius_flag = rmass_flag = 0; ellipsoid_flag = line_flag = tri_flag = body_flag = 0; vfrac_flag = 0; spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0; cs_flag = csforce_flag = vforce_flag = etag_flag = 0; rho_flag = e_flag = cv_flag = vest_flag = 0; dpd_flag = 0; // USER-SMD smd_flag = 0; contact_radius_flag = 0; smd_data_9_flag = 0; smd_stress_flag = 0; x0_flag = 0; eff_plastic_strain_flag = 0; eff_plastic_strain_rate_flag = 0; damage_flag = 0; // Peridynamic scale factor pdscale = 1.0; // ntype-length arrays mass = NULL; mass_setflag = NULL; // callback lists & extra restart info nextra_grow = nextra_restart = nextra_border = 0; extra_grow = extra_restart = extra_border = NULL; nextra_grow_max = nextra_restart_max = nextra_border_max = 0; nextra_store = 0; extra = NULL; // default atom ID and mapping values tag_enable = 1; map_style = map_user = 0; map_tag_max = -1; map_maxarray = map_nhash = -1; max_same = 0; sametag = NULL; map_array = NULL; map_bucket = NULL; map_hash = NULL; atom_style = NULL; avec = NULL; avec_map = new AtomVecCreatorMap(); #define ATOM_CLASS #define AtomStyle(key,Class) \ (*avec_map)[#key] = &avec_creator; #include "style_atom.h" #undef AtomStyle #undef ATOM_CLASS } /* ---------------------------------------------------------------------- */ Atom::~Atom() { delete [] atom_style; delete avec; delete avec_map; delete [] firstgroupname; memory->destroy(binhead); memory->destroy(next); memory->destroy(permute); // delete atom arrays // customize by adding new array memory->destroy(tag); memory->destroy(type); memory->destroy(mask); memory->destroy(image); memory->destroy(x); memory->destroy(v); memory->destroy(f); memory->destroy(molecule); memory->destroy(molindex); memory->destroy(molatom); memory->destroy(q); memory->destroy(mu); memory->destroy(omega); memory->destroy(angmom); memory->destroy(torque); memory->destroy(radius); memory->destroy(rmass); memory->destroy(ellipsoid); memory->destroy(line); memory->destroy(tri); memory->destroy(body); memory->destroy(vfrac); memory->destroy(s0); memory->destroy(x0); memory->destroy(spin); memory->destroy(eradius); memory->destroy(ervel); memory->destroy(erforce); memory->destroy(ervelforce); memory->destroy(cs); memory->destroy(csforce); memory->destroy(vforce); memory->destroy(etag); memory->destroy(rho); memory->destroy(drho); memory->destroy(e); memory->destroy(de); memory->destroy(cv); memory->destroy(vest); memory->destroy(contact_radius); memory->destroy(smd_data_9); memory->destroy(smd_stress); memory->destroy(eff_plastic_strain); memory->destroy(eff_plastic_strain_rate); memory->destroy(damage); memory->destroy(dpdTheta); memory->destroy(uCond); memory->destroy(uMech); memory->destroy(uChem); memory->destroy(uCG); memory->destroy(uCGnew); memory->destroy(duChem); memory->destroy(ssaAIR); memory->destroy(nspecial); memory->destroy(special); 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); // delete custom atom arrays for (int i = 0; i < nivector; i++) { delete [] iname[i]; memory->destroy(ivector[i]); } for (int i = 0; i < ndvector; i++) { delete [] dname[i]; memory->destroy(dvector[i]); } memory->sfree(iname); memory->sfree(dname); memory->sfree(ivector); memory->sfree(dvector); // delete user-defined molecules for (int i = 0; i < nmolecule; i++) delete molecules[i]; memory->sfree(molecules); // delete per-type arrays delete [] mass; delete [] mass_setflag; // delete extra arrays memory->destroy(extra_grow); memory->destroy(extra_restart); memory->destroy(extra_border); memory->destroy(extra); // delete mapping data structures map_delete(); } /* ---------------------------------------------------------------------- copy modify settings from old Atom class to current Atom class ------------------------------------------------------------------------- */ void Atom::settings(Atom *old) { tag_enable = old->tag_enable; map_user = old->map_user; map_style = old->map_style; sortfreq = old->sortfreq; userbinsize = old->userbinsize; if (old->firstgroupname) { int n = strlen(old->firstgroupname) + 1; firstgroupname = new char[n]; strcpy(firstgroupname,old->firstgroupname); } } /* ---------------------------------------------------------------------- create an AtomVec style called from lammps.cpp, input script, restart file, replicate ------------------------------------------------------------------------- */ void Atom::create_avec(const char *style, int narg, char **arg, int trysuffix) { delete [] atom_style; if (avec) delete avec; atom_style = NULL; avec = NULL; // unset atom style and array existence flags // may have been set by old avec // customize by adding new flag sphere_flag = peri_flag = electron_flag = 0; wavepacket_flag = sph_flag = 0; molecule_flag = 0; q_flag = mu_flag = 0; omega_flag = torque_flag = angmom_flag = 0; radius_flag = rmass_flag = 0; ellipsoid_flag = line_flag = tri_flag = body_flag = 0; vfrac_flag = 0; spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0; cs_flag = csforce_flag = vforce_flag = etag_flag = 0; rho_flag = e_flag = cv_flag = vest_flag = 0; // create instance of AtomVec // use grow() to initialize atom-based arrays to length 1 // so that x[0][0] can always be referenced even if proc has no atoms int sflag; avec = new_avec(style,trysuffix,sflag); avec->store_args(narg,arg); avec->process_args(narg,arg); avec->grow(1); if (sflag) { char estyle[256]; if (sflag == 1) sprintf(estyle,"%s/%s",style,lmp->suffix); else sprintf(estyle,"%s/%s",style,lmp->suffix2); int n = strlen(estyle) + 1; atom_style = new char[n]; strcpy(atom_style,estyle); } else { int n = strlen(style) + 1; atom_style = new char[n]; strcpy(atom_style,style); } // if molecular system: // atom IDs must be defined // force atom map to be created // map style may be reset by map_init() and its call to map_style_set() molecular = avec->molecular; if (molecular && tag_enable == 0) error->all(FLERR,"Atom IDs must be used for molecular systems"); if (molecular) map_style = 1; } /* ---------------------------------------------------------------------- generate an AtomVec class, first with suffix appended ------------------------------------------------------------------------- */ AtomVec *Atom::new_avec(const char *style, int trysuffix, int &sflag) { if (trysuffix && lmp->suffix_enable) { if (lmp->suffix) { sflag = 1; char estyle[256]; sprintf(estyle,"%s/%s",style,lmp->suffix); if (avec_map->find(estyle) != avec_map->end()) { AtomVecCreator avec_creator = (*avec_map)[estyle]; return avec_creator(lmp); } } if (lmp->suffix2) { sflag = 2; char estyle[256]; sprintf(estyle,"%s/%s",style,lmp->suffix2); if (avec_map->find(estyle) != avec_map->end()) { AtomVecCreator avec_creator = (*avec_map)[estyle]; return avec_creator(lmp); } } } sflag = 0; if (avec_map->find(style) != avec_map->end()) { AtomVecCreator avec_creator = (*avec_map)[style]; return avec_creator(lmp); } error->all(FLERR,"Unknown atom style"); return NULL; } /* ---------------------------------------------------------------------- one instance per AtomVec style in style_atom.h ------------------------------------------------------------------------- */ template AtomVec *Atom::avec_creator(LAMMPS *lmp) { return new T(lmp); } /* ---------------------------------------------------------------------- */ void Atom::init() { // delete extra array since it doesn't persist past first run if (nextra_store) { memory->destroy(extra); extra = NULL; nextra_store = 0; } // check arrays that are atom type in length check_mass(FLERR); // setup of firstgroup if (firstgroupname) { firstgroup = group->find(firstgroupname); if (firstgroup < 0) error->all(FLERR,"Could not find atom_modify first group ID"); } else firstgroup = -1; // init AtomVec avec->init(); } /* ---------------------------------------------------------------------- */ void Atom::setup() { // setup bins for sorting // cannot do this in init() because uses neighbor cutoff if (sortfreq > 0) setup_sort_bins(); } /* ---------------------------------------------------------------------- return ptr to AtomVec class if matches style or to matching hybrid sub-class return NULL if no match ------------------------------------------------------------------------- */ AtomVec *Atom::style_match(const char *style) { if (strcmp(atom_style,style) == 0) return avec; else if (strcmp(atom_style,"hybrid") == 0) { AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) avec; for (int i = 0; i < avec_hybrid->nstyles; i++) if (strcmp(avec_hybrid->keywords[i],style) == 0) return avec_hybrid->styles[i]; } return NULL; } /* ---------------------------------------------------------------------- modify parameters of the atom style some options can only be invoked before simulation box is defined first and sort options cannot be used together ------------------------------------------------------------------------- */ void Atom::modify_params(int narg, char **arg) { if (narg == 0) error->all(FLERR,"Illegal atom_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"id") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); if (domain->box_exist) error->all(FLERR, "Atom_modify id command after simulation box is defined"); if (strcmp(arg[iarg+1],"yes") == 0) tag_enable = 1; else if (strcmp(arg[iarg+1],"no") == 0) tag_enable = 0; else error->all(FLERR,"Illegal atom_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"map") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); if (domain->box_exist) error->all(FLERR, "Atom_modify map command after simulation box is defined"); if (strcmp(arg[iarg+1],"array") == 0) map_user = 1; else if (strcmp(arg[iarg+1],"hash") == 0) map_user = 2; else error->all(FLERR,"Illegal atom_modify command"); map_style = map_user; iarg += 2; } else if (strcmp(arg[iarg],"first") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); if (strcmp(arg[iarg+1],"all") == 0) { delete [] firstgroupname; firstgroupname = NULL; } else { int n = strlen(arg[iarg+1]) + 1; firstgroupname = new char[n]; strcpy(firstgroupname,arg[iarg+1]); sortfreq = 0; } iarg += 2; } else if (strcmp(arg[iarg],"sort") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal atom_modify command"); sortfreq = force->inumeric(FLERR,arg[iarg+1]); userbinsize = force->numeric(FLERR,arg[iarg+2]); if (sortfreq < 0 || userbinsize < 0.0) error->all(FLERR,"Illegal atom_modify command"); if (sortfreq >= 0 && firstgroupname) error->all(FLERR,"Atom_modify sort and first options " "cannot be used together"); iarg += 3; } else error->all(FLERR,"Illegal atom_modify command"); } } /* ---------------------------------------------------------------------- check that atom IDs are valid error if any atom ID < 0 or atom ID = MAXTAGINT if any atom ID > 0, error if any atom ID == 0 if any atom ID > 0, error if tag_enable = 0 if all atom IDs = 0, tag_enable must be 0 if max atom IDs < natoms, must be duplicates OK if max atom IDs > natoms NOTE: not fully checking that atom IDs are unique ------------------------------------------------------------------------- */ void Atom::tag_check() { tagint min = MAXTAGINT; tagint max = 0; for (int i = 0; i < nlocal; i++) { min = MIN(min,tag[i]); max = MAX(max,tag[i]); } tagint minall,maxall; MPI_Allreduce(&min,&minall,1,MPI_LMP_TAGINT,MPI_MIN,world); MPI_Allreduce(&max,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world); if (minall < 0) error->all(FLERR,"One or more Atom IDs is negative"); if (maxall >= MAXTAGINT) error->all(FLERR,"One or more atom IDs is too big"); if (maxall > 0 && minall == 0) error->all(FLERR,"One or more atom IDs is zero"); if (maxall > 0 && tag_enable == 0) error->all(FLERR,"Non-zero atom IDs with atom_modify id = no"); if (maxall == 0 && natoms && tag_enable) error->all(FLERR,"All atom IDs = 0 but atom_modify id = yes"); if (tag_enable && maxall < natoms) error->all(FLERR,"Duplicate atom IDs exist"); } /* ---------------------------------------------------------------------- add unique tags to any atoms with tag = 0 new tags are grouped by proc and start after max current tag called after creating new atoms error if new tags will exceed MAXTAGINT ------------------------------------------------------------------------- */ void Atom::tag_extend() { // maxtag_all = max tag for all atoms tagint maxtag = 0; for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]); tagint maxtag_all; MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world); // DEBUG: useful for generating 64-bit IDs even for small systems // use only when LAMMPS is compiled with BIGBIG //maxtag_all += 1000000000000; // notag = # of atoms I own with no tag (tag = 0) // notag_sum = # of total atoms on procs <= me with no tag bigint notag = 0; for (int i = 0; i < nlocal; i++) if (tag[i] == 0) notag++; bigint notag_total; MPI_Allreduce(¬ag,¬ag_total,1,MPI_LMP_BIGINT,MPI_SUM,world); if (notag_total >= MAXTAGINT) error->all(FLERR,"New atom IDs exceed maximum allowed ID"); bigint notag_sum; MPI_Scan(¬ag,¬ag_sum,1,MPI_LMP_BIGINT,MPI_SUM,world); // itag = 1st new tag that my untagged atoms should use tagint itag = maxtag_all + notag_sum - notag + 1; for (int i = 0; i < nlocal; i++) if (tag[i] == 0) tag[i] = itag++; } /* ---------------------------------------------------------------------- check that atom IDs span range from 1 to Natoms inclusive return 0 if mintag != 1 or maxtag != Natoms return 1 if OK doesn't actually check if all tag values are used ------------------------------------------------------------------------- */ int Atom::tag_consecutive() { tagint idmin = MAXTAGINT; tagint idmax = 0; for (int i = 0; i < nlocal; i++) { idmin = MIN(idmin,tag[i]); idmax = MAX(idmax,tag[i]); } tagint idminall,idmaxall; MPI_Allreduce(&idmin,&idminall,1,MPI_LMP_TAGINT,MPI_MIN,world); MPI_Allreduce(&idmax,&idmaxall,1,MPI_LMP_TAGINT,MPI_MAX,world); if (idminall != 1 || idmaxall != natoms) return 0; return 1; } /* ---------------------------------------------------------------------- count and return words in a single line make copy of line before using strtok so as not to change line trim anything from '#' onward ------------------------------------------------------------------------- */ int Atom::count_words(const char *line) { int n = strlen(line) + 1; char *copy; memory->create(copy,n,"atom:copy"); strcpy(copy,line); char *ptr; if ((ptr = strchr(copy,'#'))) *ptr = '\0'; if (strtok(copy," \t\n\r\f") == NULL) { memory->destroy(copy); return 0; } n = 1; while (strtok(NULL," \t\n\r\f")) n++; memory->destroy(copy); return n; } /* ---------------------------------------------------------------------- count and return words in a single line using provided copy buf make copy of line before using strtok so as not to change line trim anything from '#' onward ------------------------------------------------------------------------- */ int Atom::count_words(const char *line, char *copy) { strcpy(copy,line); char *ptr; if ((ptr = strchr(copy,'#'))) *ptr = '\0'; if (strtok(copy," \t\n\r\f") == NULL) { memory->destroy(copy); return 0; } int n = 1; while (strtok(NULL," \t\n\r\f")) n++; return n; } /* ---------------------------------------------------------------------- deallocate molecular topology arrays done before realloc with (possibly) new 2nd dimension set to correctly initialized per-atom values, e.g. bond_per_atom needs to be called whenever 2nd dimensions are changed and these arrays are already pre-allocated, e.g. due to grow(1) in create_avec() ------------------------------------------------------------------------- */ void Atom::deallocate_topology() { memory->destroy(atom->bond_type); memory->destroy(atom->bond_atom); atom->bond_type = NULL; atom->bond_atom = NULL; memory->destroy(atom->angle_type); memory->destroy(atom->angle_atom1); memory->destroy(atom->angle_atom2); memory->destroy(atom->angle_atom3); atom->angle_type = NULL; atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL; memory->destroy(atom->dihedral_type); memory->destroy(atom->dihedral_atom1); memory->destroy(atom->dihedral_atom2); memory->destroy(atom->dihedral_atom3); memory->destroy(atom->dihedral_atom4); atom->dihedral_type = NULL; atom->dihedral_atom1 = atom->dihedral_atom2 = atom->dihedral_atom3 = atom->dihedral_atom4 = NULL; memory->destroy(atom->improper_type); memory->destroy(atom->improper_atom1); memory->destroy(atom->improper_atom2); memory->destroy(atom->improper_atom3); memory->destroy(atom->improper_atom4); atom->improper_type = NULL; atom->improper_atom1 = atom->improper_atom2 = atom->improper_atom3 = atom->improper_atom4 = NULL; } /* ---------------------------------------------------------------------- unpack N lines from Atom section of data file call style-specific routine to parse line ------------------------------------------------------------------------- */ void Atom::data_atoms(int n, char *buf, tagint id_offset, int type_offset, int shiftflag, double *shift) { int m,xptr,iptr; imageint imagedata; double xdata[3],lamda[3]; double *coord; char *next; next = strchr(buf,'\n'); *next = '\0'; int nwords = count_words(buf); *next = '\n'; if (nwords != avec->size_data_atom && nwords != avec->size_data_atom + 3) error->all(FLERR,"Incorrect atom format in data file"); char **values = new char*[nwords]; // set bounds for my proc // if periodic and I am lo/hi proc, adjust bounds by EPSILON // insures all data atoms will be owned even with round-off int triclinic = domain->triclinic; double epsilon[3]; if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON; else { epsilon[0] = domain->prd[0] * EPSILON; epsilon[1] = domain->prd[1] * EPSILON; epsilon[2] = domain->prd[2] * EPSILON; } double sublo[3],subhi[3]; if (triclinic == 0) { sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0]; sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1]; sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2]; } else { sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0]; sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1]; sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2]; } if (comm->layout != LAYOUT_TILED) { if (domain->xperiodic) { if (comm->myloc[0] == 0) sublo[0] -= epsilon[0]; if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += epsilon[0]; } if (domain->yperiodic) { if (comm->myloc[1] == 0) sublo[1] -= epsilon[1]; if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += epsilon[1]; } if (domain->zperiodic) { if (comm->myloc[2] == 0) sublo[2] -= epsilon[2]; if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += epsilon[2]; } } else { if (domain->xperiodic) { if (comm->mysplit[0][0] == 0.0) sublo[0] -= epsilon[0]; if (comm->mysplit[0][1] == 1.0) subhi[0] += epsilon[0]; } if (domain->yperiodic) { if (comm->mysplit[1][0] == 0.0) sublo[1] -= epsilon[1]; if (comm->mysplit[1][1] == 1.0) subhi[1] += epsilon[1]; } if (domain->zperiodic) { if (comm->mysplit[2][0] == 0.0) sublo[2] -= epsilon[2]; if (comm->mysplit[2][1] == 1.0) subhi[2] += epsilon[2]; } } // xptr = which word in line starts xyz coords // iptr = which word in line starts ix,iy,iz image flags xptr = avec->xcol_data - 1; int imageflag = 0; if (nwords > avec->size_data_atom) imageflag = 1; if (imageflag) iptr = nwords - 3; // loop over lines of atom data // tokenize the line into values // extract xyz coords and image flags // remap atom into simulation box // if atom is in my sub-domain, unpack its values for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); values[0] = strtok(buf," \t\n\r\f"); if (values[0] == NULL) error->all(FLERR,"Incorrect atom format in data file"); for (m = 1; m < nwords; m++) { values[m] = strtok(NULL," \t\n\r\f"); if (values[m] == NULL) error->all(FLERR,"Incorrect atom format in data file"); } if (imageflag) imagedata = ((imageint) (atoi(values[iptr]) + IMGMAX) & IMGMASK) | (((imageint) (atoi(values[iptr+1]) + IMGMAX) & IMGMASK) << IMGBITS) | (((imageint) (atoi(values[iptr+2]) + IMGMAX) & IMGMASK) << IMG2BITS); else imagedata = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; xdata[0] = atof(values[xptr]); xdata[1] = atof(values[xptr+1]); xdata[2] = atof(values[xptr+2]); if (shiftflag) { xdata[0] += shift[0]; xdata[1] += shift[1]; xdata[2] += shift[2]; } domain->remap(xdata,imagedata); if (triclinic) { domain->x2lamda(xdata,lamda); coord = lamda; } else coord = xdata; if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[2] >= sublo[2] && coord[2] < subhi[2]) { avec->data_atom(xdata,imagedata,values); if (id_offset) tag[nlocal-1] += id_offset; if (type_offset) { type[nlocal-1] += type_offset; if (type[nlocal-1] > ntypes) error->one(FLERR,"Invalid atom type in Atoms section of data file"); } } buf = next + 1; } delete [] values; } /* ---------------------------------------------------------------------- unpack N lines from Velocity section of data file check that atom IDs are > 0 and <= map_tag_max call style-specific routine to parse line ------------------------------------------------------------------------- */ void Atom::data_vels(int n, char *buf, tagint id_offset) { int j,m; tagint tagdata; char *next; next = strchr(buf,'\n'); *next = '\0'; int nwords = count_words(buf); *next = '\n'; if (nwords != avec->size_data_vel) error->all(FLERR,"Incorrect velocity format in data file"); char **values = new char*[nwords]; // loop over lines of atom velocities // tokenize the line into values // if I own atom tag, unpack its values for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); values[0] = strtok(buf," \t\n\r\f"); for (j = 1; j < nwords; j++) values[j] = strtok(NULL," \t\n\r\f"); tagdata = ATOTAGINT(values[0]) + id_offset; if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Velocities section of data file"); if ((m = map(tagdata)) >= 0) avec->data_vel(m,&values[1]); buf = next + 1; } delete [] values; } /* ---------------------------------------------------------------------- process N bonds read into buf from data files if count is non-NULL, just count bonds per atom else store them with atoms check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset, int type_offset) { int m,tmp,itype; tagint atom1,atom2; char *next; int newton_bond = force->newton_bond; for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2); if (id_offset) { atom1 += id_offset; atom2 += id_offset; } itype += type_offset; if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max) error->one(FLERR,"Invalid atom ID in Bonds section of data file"); if (itype <= 0 || itype > nbondtypes) error->one(FLERR,"Invalid bond type in Bonds section of data file"); if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { bond_type[m][num_bond[m]] = itype; bond_atom[m][num_bond[m]] = atom2; num_bond[m]++; } } if (newton_bond == 0) { if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { bond_type[m][num_bond[m]] = itype; bond_atom[m][num_bond[m]] = atom1; num_bond[m]++; } } } buf = next + 1; } } /* ---------------------------------------------------------------------- process N angles read into buf from data files if count is non-NULL, just count angles per atom else store them with atoms check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset, int type_offset) { int m,tmp,itype; tagint atom1,atom2,atom3; char *next; int newton_bond = force->newton_bond; for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3); if (id_offset) { atom1 += id_offset; atom2 += id_offset; atom3 += id_offset; } itype += type_offset; if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max) error->one(FLERR,"Invalid atom ID in Angles section of data file"); if (itype <= 0 || itype > nangletypes) error->one(FLERR,"Invalid angle type in Angles section of data file"); if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { 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) { if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { 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 ((m = map(atom3)) >= 0) { if (count) count[m]++; else { 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]++; } } } buf = next + 1; } } /* ---------------------------------------------------------------------- process N dihedrals read into buf from data files if count is non-NULL, just count diihedrals per atom else store them with atoms check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset, int type_offset) { int m,tmp,itype; tagint atom1,atom2,atom3,atom4; char *next; int newton_bond = force->newton_bond; for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3,&atom4); if (id_offset) { atom1 += id_offset; atom2 += id_offset; atom3 += id_offset; atom4 += id_offset; } itype += type_offset; if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max || atom4 <= 0 || atom4 > map_tag_max) error->one(FLERR,"Invalid atom ID in Dihedrals section of data file"); if (itype <= 0 || itype > ndihedraltypes) error->one(FLERR, "Invalid dihedral type in Dihedrals section of data file"); if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { 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; dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } if (newton_bond == 0) { if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { 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; dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } if ((m = map(atom3)) >= 0) { if (count) count[m]++; else { 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; dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } if ((m = map(atom4)) >= 0) { if (count) count[m]++; else { 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; dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } } buf = next + 1; } } /* ---------------------------------------------------------------------- process N impropers read into buf from data files if count is non-NULL, just count impropers per atom else store them with atoms check that atom IDs are > 0 and <= map_tag_max ------------------------------------------------------------------------- */ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset, int type_offset) { int m,tmp,itype; tagint atom1,atom2,atom3,atom4; char *next; int newton_bond = force->newton_bond; for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3,&atom4); if (id_offset) { atom1 += id_offset; atom2 += id_offset; atom3 += id_offset; atom4 += id_offset; } itype += type_offset; if (atom1 <= 0 || atom1 > map_tag_max || atom2 <= 0 || atom2 > map_tag_max || atom3 <= 0 || atom3 > map_tag_max || atom4 <= 0 || atom4 > map_tag_max) error->one(FLERR,"Invalid atom ID in Impropers section of data file"); if (itype <= 0 || itype > nimpropertypes) error->one(FLERR, "Invalid improper type in Impropers section of data file"); if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { 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; improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } if (newton_bond == 0) { if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { 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; improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } if ((m = map(atom3)) >= 0) { if (count) count[m]++; else { 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; improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } if ((m = map(atom4)) >= 0) { if (count) count[m]++; else { 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; improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } } buf = next + 1; } } /* ---------------------------------------------------------------------- unpack N lines from atom-style specific bonus section of data file check that atom IDs are > 0 and <= map_tag_max call style-specific routine to parse line ------------------------------------------------------------------------- */ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset) { int j,m,tagdata; char *next; next = strchr(buf,'\n'); *next = '\0'; int nwords = count_words(buf); *next = '\n'; if (nwords != avec_bonus->size_data_bonus) error->all(FLERR,"Incorrect bonus data format in data file"); char **values = new char*[nwords]; // loop over lines of bonus atom data // tokenize the line into values // if I own atom tag, unpack its values for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); values[0] = strtok(buf," \t\n\r\f"); for (j = 1; j < nwords; j++) values[j] = strtok(NULL," \t\n\r\f"); tagdata = ATOTAGINT(values[0]) + id_offset; if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Bonus section of data file"); // ok to call child's data_atom_bonus() method thru parent avec_bonus, // since data_bonus() was called with child ptr, and method is virtual if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,&values[1]); buf = next + 1; } delete [] values; } /* ---------------------------------------------------------------------- unpack N bodies from Bodies section of data file each body spans multiple lines check that atom IDs are > 0 and <= map_tag_max call style-specific routine to parse line ------------------------------------------------------------------------- */ void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body, tagint id_offset) { int j,m,nvalues,tagdata,ninteger,ndouble; int maxint = 0; int maxdouble = 0; int *ivalues = NULL; double *dvalues = NULL; // loop over lines of body data // if I own atom tag, tokenize lines into ivalues/dvalues, call data_body() // else skip values for (int i = 0; i < n; i++) { if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")) + id_offset; else tagdata = ATOTAGINT(strtok(NULL," \t\n\r\f")) + id_offset; if (tagdata <= 0 || tagdata > map_tag_max) error->one(FLERR,"Invalid atom ID in Bodies section of data file"); ninteger = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); ndouble = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); if ((m = map(tagdata)) >= 0) { if (ninteger > maxint) { delete [] ivalues; maxint = ninteger; ivalues = new int[maxint]; } if (ndouble > maxdouble) { delete [] dvalues; maxdouble = ndouble; dvalues = new double[maxdouble]; } for (j = 0; j < ninteger; j++) ivalues[j] = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); for (j = 0; j < ndouble; j++) dvalues[j] = force->numeric(FLERR,strtok(NULL," \t\n\r\f")); avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues); } else { nvalues = ninteger + ndouble; // number of values to skip for (j = 0; j < nvalues; j++) strtok(NULL," \t\n\r\f"); } } delete [] ivalues; delete [] dvalues; } /* ---------------------------------------------------------------------- init per-atom fix/compute/variable values for newly created atoms called from create_atoms, read_data, read_dump, lib::lammps_create_atoms() fixes, computes, variables may or may not exist when called ------------------------------------------------------------------------- */ void Atom::data_fix_compute_variable(int nprev, int nnew) { for (int m = 0; m < modify->nfix; m++) { Fix *fix = modify->fix[m]; if (fix->create_attribute) for (int i = nprev; i < nnew; i++) fix->set_arrays(i); } for (int m = 0; m < modify->ncompute; m++) { Compute *compute = modify->compute[m]; if (compute->create_attribute) for (int i = nprev; i < nnew; i++) compute->set_arrays(i); } for (int i = nprev; i < nnew; i++) input->variable->set_arrays(i); } /* ---------------------------------------------------------------------- allocate arrays of length ntypes only done after ntypes is set ------------------------------------------------------------------------- */ void Atom::allocate_type_arrays() { if (avec->mass_type) { mass = new double[ntypes+1]; mass_setflag = new int[ntypes+1]; for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0; } } /* ---------------------------------------------------------------------- set a mass and flag it as set called from reading of data file type_offset may be used when reading multiple data files ------------------------------------------------------------------------- */ void Atom::set_mass(const char *file, int line, const char *str, int type_offset) { if (mass == NULL) error->all(file,line,"Cannot set mass for this atom style"); int itype; double mass_one; int n = sscanf(str,"%d %lg",&itype,&mass_one); if (n != 2) error->all(file,line,"Invalid mass line in data file"); itype += type_offset; if (itype < 1 || itype > ntypes) error->all(file,line,"Invalid type for mass set"); mass[itype] = mass_one; mass_setflag[itype] = 1; if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value"); } /* ---------------------------------------------------------------------- set a mass and flag it as set called from EAM pair routine ------------------------------------------------------------------------- */ void Atom::set_mass(const char *file, int line, int itype, double value) { if (mass == NULL) error->all(file,line,"Cannot set mass for this atom style"); if (itype < 1 || itype > ntypes) error->all(file,line,"Invalid type for mass set"); mass[itype] = value; mass_setflag[itype] = 1; if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value"); } /* ---------------------------------------------------------------------- set one or more masses and flag them as set called from reading of input script ------------------------------------------------------------------------- */ void Atom::set_mass(const char *file, int line, int narg, char **arg) { if (mass == NULL) error->all(file,line,"Cannot set mass for this atom style"); int lo,hi; force->bounds(file,line,arg[0],ntypes,lo,hi); if (lo < 1 || hi > ntypes) error->all(file,line,"Invalid type for mass set"); for (int itype = lo; itype <= hi; itype++) { mass[itype] = atof(arg[1]); mass_setflag[itype] = 1; if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value"); } } /* ---------------------------------------------------------------------- set all masses as read in from restart file ------------------------------------------------------------------------- */ void Atom::set_mass(double *values) { for (int itype = 1; itype <= ntypes; itype++) { mass[itype] = values[itype]; mass_setflag[itype] = 1; } } /* ---------------------------------------------------------------------- - check that all masses have been set + check that all per-atom-type masses have been set ------------------------------------------------------------------------- */ void Atom::check_mass(const char *file, int line) { if (mass == NULL) return; + if (rmass_flag) return; for (int itype = 1; itype <= ntypes; itype++) if (mass_setflag[itype] == 0) error->all(file,line,"Not all per-type masses are set"); } /* ---------------------------------------------------------------------- check that radii of all particles of itype are the same return 1 if true, else return 0 also return the radius value for that type ------------------------------------------------------------------------- */ int Atom::radius_consistency(int itype, double &rad) { double value = -1.0; int flag = 0; for (int i = 0; i < nlocal; i++) { if (type[i] != itype) continue; if (value < 0.0) value = radius[i]; else if (value != radius[i]) flag = 1; } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall) return 0; MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world); return 1; } /* ---------------------------------------------------------------------- check that shape of all particles of itype are the same return 1 if true, else return 0 also return the 3 shape params for itype ------------------------------------------------------------------------- */ int Atom::shape_consistency(int itype, double &shapex, double &shapey, double &shapez) { double zero[3] = {0.0, 0.0, 0.0}; double one[3] = {-1.0, -1.0, -1.0}; double *shape; AtomVecEllipsoid *avec_ellipsoid = (AtomVecEllipsoid *) style_match("ellipsoid"); AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus; int flag = 0; for (int i = 0; i < nlocal; i++) { if (type[i] != itype) continue; if (ellipsoid[i] < 0) shape = zero; else shape = bonus[ellipsoid[i]].shape; if (one[0] < 0.0) { one[0] = shape[0]; one[1] = shape[1]; one[2] = shape[2]; } else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2]) flag = 1; } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall) return 0; double oneall[3]; MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world); shapex = oneall[0]; shapey = oneall[1]; shapez = oneall[2]; return 1; } /* ---------------------------------------------------------------------- add a new molecule template = set of molecules ------------------------------------------------------------------------- */ 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 template ID"); // 1st molecule in set stores nset = # of mols, others store nset = 0 // ifile = count of molecules in set // index = argument index where next molecule starts, updated by constructor int ifile = 1; int index = 1; while (1) { molecules = (Molecule **) memory->srealloc(molecules,(nmolecule+1)*sizeof(Molecule *), "atom::molecules"); molecules[nmolecule] = new Molecule(lmp,narg,arg,index); molecules[nmolecule]->nset = 0; molecules[nmolecule-ifile+1]->nset++; nmolecule++; if (molecules[nmolecule-1]->last) break; ifile++; } } /* ---------------------------------------------------------------------- find first molecule in set with template ID return -1 if does not exist ------------------------------------------------------------------------- */ int Atom::find_molecule(char *id) { if(id == NULL) return -1; int imol; for (imol = 0; imol < nmolecule; imol++) if (strcmp(id,molecules[imol]->id) == 0) return imol; return -1; } /* ---------------------------------------------------------------------- add info to current atom ilocal from molecule template onemol and its iatom offset = atom ID preceding IDs of atoms in this molecule called by fixes and commands that add molecules ------------------------------------------------------------------------- */ void Atom::add_molecule_atom(Molecule *onemol, int iatom, int ilocal, tagint offset) { if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom]; if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom]; if (onemol->rmassflag && rmass_flag) 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->bodyflag) { body[ilocal] = 0; // as if a body read from data file onemol->avec_body->data_body(ilocal,onemol->nibody,onemol->ndbody, onemol->ibodyparams,onemol->dbodyparams); onemol->avec_body->set_quat(ilocal,onemol->quat_external); } if (molecular != 1) return; // add bond topology info // for molecular atom styles, but not atom style template if (avec->bonds_allow) { 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 (avec->angles_allow) { 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 (avec->dihedrals_allow) { 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 (avec->impropers_allow) { 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 only owned atoms exist at this point, no ghost atoms ------------------------------------------------------------------------- */ void Atom::first_reorder() { // insure there is one extra atom location at end of arrays for swaps if (nlocal == nmax) avec->grow(0); // loop over owned atoms // nfirst = index of first atom not in firstgroup // when find firstgroup atom out of place, swap it with atom nfirst int bitmask = group->bitmask[firstgroup]; nfirst = 0; while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++; for (int i = 0; i < nlocal; i++) { if (mask[i] & bitmask && i > nfirst) { avec->copy(i,nlocal,0); avec->copy(nfirst,i,0); avec->copy(nlocal,nfirst,0); while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++; } } } /* ---------------------------------------------------------------------- perform spatial sort of atoms within my sub-domain always called between comm->exchange() and comm->borders() don't have to worry about clearing/setting atom->map since done in comm ------------------------------------------------------------------------- */ void Atom::sort() { int i,m,n,ix,iy,iz,ibin,empty; // set next timestep for sorting to take place nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq; // re-setup sort bins if needed if (domain->box_change) setup_sort_bins(); if (nbins == 1) return; // reallocate per-atom vectors if needed if (nlocal > maxnext) { memory->destroy(next); memory->destroy(permute); maxnext = atom->nmax; memory->create(next,maxnext,"atom:next"); memory->create(permute,maxnext,"atom:permute"); } // insure there is one extra atom location at end of arrays for swaps if (nlocal == nmax) avec->grow(0); // bin atoms in reverse order so linked list will be in forward order for (i = 0; i < nbins; i++) binhead[i] = -1; for (i = nlocal-1; i >= 0; i--) { ix = static_cast ((x[i][0]-bboxlo[0])*bininvx); iy = static_cast ((x[i][1]-bboxlo[1])*bininvy); iz = static_cast ((x[i][2]-bboxlo[2])*bininvz); ix = MAX(ix,0); iy = MAX(iy,0); iz = MAX(iz,0); ix = MIN(ix,nbinx-1); iy = MIN(iy,nbiny-1); iz = MIN(iz,nbinz-1); ibin = iz*nbiny*nbinx + iy*nbinx + ix; next[i] = binhead[ibin]; binhead[ibin] = i; } // permute = desired permutation of atoms // permute[I] = J means Ith new atom will be Jth old atom n = 0; for (m = 0; m < nbins; m++) { i = binhead[m]; while (i >= 0) { permute[n++] = i; i = next[i]; } } // current = current permutation, just reuse next vector // current[I] = J means Ith current atom is Jth old atom int *current = next; for (i = 0; i < nlocal; i++) current[i] = i; // reorder local atom list, when done, current = permute // perform "in place" using copy() to extra atom location at end of list // inner while loop processes one cycle of the permutation // copy before inner-loop moves an atom to end of atom list // copy after inner-loop moves atom at end of list back into list // empty = location in atom list that is currently empty for (i = 0; i < nlocal; i++) { if (current[i] == permute[i]) continue; avec->copy(i,nlocal,0); empty = i; while (permute[empty] != i) { avec->copy(permute[empty],empty,0); empty = current[empty] = permute[empty]; } avec->copy(nlocal,empty,0); current[empty] = permute[empty]; } // sanity check that current = permute //int flag = 0; //for (i = 0; i < nlocal; i++) // if (current[i] != permute[i]) flag = 1; //int flagall; //MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); //if (flagall) error->all(FLERR,"Atom sort did not operate correctly"); } /* ---------------------------------------------------------------------- setup bins for spatial sorting of atoms ------------------------------------------------------------------------- */ void Atom::setup_sort_bins() { // binsize: // user setting if explicitly set // default = 1/2 of neighbor cutoff // check if neighbor cutoff = 0.0 double binsize; if (userbinsize > 0.0) binsize = userbinsize; else binsize = 0.5 * neighbor->cutneighmax; if (binsize == 0.0) error->all(FLERR,"Atom sorting has bin size = 0.0"); double bininv = 1.0/binsize; // nbin xyz = local bins // bbox lo/hi = bounding box of my sub-domain if (domain->triclinic) domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi); else { bboxlo[0] = domain->sublo[0]; bboxlo[1] = domain->sublo[1]; bboxlo[2] = domain->sublo[2]; bboxhi[0] = domain->subhi[0]; bboxhi[1] = domain->subhi[1]; bboxhi[2] = domain->subhi[2]; } nbinx = static_cast ((bboxhi[0]-bboxlo[0]) * bininv); nbiny = static_cast ((bboxhi[1]-bboxlo[1]) * bininv); nbinz = static_cast ((bboxhi[2]-bboxlo[2]) * bininv); if (domain->dimension == 2) nbinz = 1; if (nbinx == 0) nbinx = 1; if (nbiny == 0) nbiny = 1; if (nbinz == 0) nbinz = 1; bininvx = nbinx / (bboxhi[0]-bboxlo[0]); bininvy = nbiny / (bboxhi[1]-bboxlo[1]); bininvz = nbinz / (bboxhi[2]-bboxlo[2]); #ifdef LMP_USER_INTEL int intel_neigh = 0; if (neighbor->nrequest) { if (neighbor->requests[0]->intel) intel_neigh = 1; } else if (neighbor->old_nrequest) if (neighbor->old_requests[0]->intel) intel_neigh = 1; if (intel_neigh && userbinsize == 0.0) { if (neighbor->binsizeflag) bininv = 1.0/neighbor->binsize_user; double nx_low = neighbor->bboxlo[0]; double ny_low = neighbor->bboxlo[1]; double nz_low = neighbor->bboxlo[2]; double nxbbox = neighbor->bboxhi[0] - nx_low; double nybbox = neighbor->bboxhi[1] - ny_low; double nzbbox = neighbor->bboxhi[2] - nz_low; int nnbinx = static_cast (nxbbox * bininv); int nnbiny = static_cast (nybbox * bininv); int nnbinz = static_cast (nzbbox * bininv); if (domain->dimension == 2) nnbinz = 1; if (nnbinx == 0) nnbinx = 1; if (nnbiny == 0) nnbiny = 1; if (nnbinz == 0) nnbinz = 1; double binsizex = nxbbox/nnbinx; double binsizey = nybbox/nnbiny; double binsizez = nzbbox/nnbinz; bininvx = 1.0 / binsizex; bininvy = 1.0 / binsizey; bininvz = 1.0 / binsizez; int lxo = (bboxlo[0] - nx_low) * bininvx; int lyo = (bboxlo[1] - ny_low) * bininvy; int lzo = (bboxlo[2] - nz_low) * bininvz; bboxlo[0] = nx_low + static_cast(lxo) / bininvx; bboxlo[1] = ny_low + static_cast(lyo) / bininvy; bboxlo[2] = nz_low + static_cast(lzo) / bininvz; nbinx = static_cast((bboxhi[0] - bboxlo[0]) * bininvx) + 1; nbiny = static_cast((bboxhi[1] - bboxlo[1]) * bininvy) + 1; nbinz = static_cast((bboxhi[2] - bboxlo[2]) * bininvz) + 1; bboxhi[0] = bboxlo[0] + static_cast(nbinx) / bininvx; bboxhi[1] = bboxlo[1] + static_cast(nbiny) / bininvy; bboxhi[2] = bboxlo[2] + static_cast(nbinz) / bininvz; } #endif if (1.0*nbinx*nbiny*nbinz > INT_MAX) error->one(FLERR,"Too many atom sorting bins"); nbins = nbinx*nbiny*nbinz; // reallocate per-bin memory if needed if (nbins > maxbin) { memory->destroy(binhead); maxbin = nbins; memory->create(binhead,maxbin,"atom:binhead"); } } /* ---------------------------------------------------------------------- register a callback to a fix so it can manage atom-based arrays happens when fix is created flag = 0 for grow, 1 for restart, 2 for border comm ------------------------------------------------------------------------- */ void Atom::add_callback(int flag) { int ifix; // find the fix // if find NULL ptr: // it's this one, since it is being replaced and has just been deleted // at this point in re-creation // if don't find NULL ptr: // i is set to nfix = new one currently being added at end of list for (ifix = 0; ifix < modify->nfix; ifix++) if (modify->fix[ifix] == NULL) break; // add callback to lists, reallocating if necessary if (flag == 0) { if (nextra_grow == nextra_grow_max) { nextra_grow_max += DELTA; memory->grow(extra_grow,nextra_grow_max,"atom:extra_grow"); } extra_grow[nextra_grow] = ifix; nextra_grow++; } else if (flag == 1) { if (nextra_restart == nextra_restart_max) { nextra_restart_max += DELTA; memory->grow(extra_restart,nextra_restart_max,"atom:extra_restart"); } extra_restart[nextra_restart] = ifix; nextra_restart++; } else if (flag == 2) { if (nextra_border == nextra_border_max) { nextra_border_max += DELTA; memory->grow(extra_border,nextra_border_max,"atom:extra_border"); } extra_border[nextra_border] = ifix; nextra_border++; } } /* ---------------------------------------------------------------------- unregister a callback to a fix happens when fix is deleted, called by its destructor flag = 0 for grow, 1 for restart ------------------------------------------------------------------------- */ void Atom::delete_callback(const char *id, int flag) { if (id == NULL) return; int ifix; for (ifix = 0; ifix < modify->nfix; ifix++) if (strcmp(id,modify->fix[ifix]->id) == 0) break; // compact the list of callbacks if (flag == 0) { int match; for (match = 0; match < nextra_grow; match++) if (extra_grow[match] == ifix) break; for (int i = match; i < nextra_grow-1; i++) extra_grow[i] = extra_grow[i+1]; nextra_grow--; } else if (flag == 1) { int match; for (match = 0; match < nextra_restart; match++) if (extra_restart[match] == ifix) break; for (int i = match; i < nextra_restart-1; i++) extra_restart[i] = extra_restart[i+1]; nextra_restart--; } else if (flag == 2) { int match; for (match = 0; match < nextra_border; match++) if (extra_border[match] == ifix) break; for (int i = match; i < nextra_border-1; i++) extra_border[i] = extra_border[i+1]; nextra_border--; } } /* ---------------------------------------------------------------------- decrement ptrs in callback lists to fixes beyond the deleted ifix happens after fix is deleted ------------------------------------------------------------------------- */ void Atom::update_callback(int ifix) { for (int i = 0; i < nextra_grow; i++) if (extra_grow[i] > ifix) extra_grow[i]--; for (int i = 0; i < nextra_restart; i++) if (extra_restart[i] > ifix) extra_restart[i]--; for (int i = 0; i < nextra_border; i++) if (extra_border[i] > ifix) extra_border[i]--; } /* ---------------------------------------------------------------------- find custom per-atom vector with name return index if found, and flag = 0/1 for int/double return -1 if not found ------------------------------------------------------------------------- */ int Atom::find_custom(const char *name, int &flag) { if(name == NULL) return -1; for (int i = 0; i < nivector; i++) if (iname[i] && strcmp(iname[i],name) == 0) { flag = 0; return i; } for (int i = 0; i < ndvector; i++) if (dname[i] && strcmp(dname[i],name) == 0) { flag = 1; return i; } return -1; } /* ---------------------------------------------------------------------- add a custom variable with name of type flag = 0/1 for int/double assumes name does not already exist return index in ivector or dvector of its location ------------------------------------------------------------------------- */ int Atom::add_custom(const char *name, int flag) { int index; if (flag == 0) { index = nivector; nivector++; iname = (char **) memory->srealloc(iname,nivector*sizeof(char *), "atom:iname"); int n = strlen(name) + 1; iname[index] = new char[n]; strcpy(iname[index],name); ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *), "atom:ivector"); memory->create(ivector[index],nmax,"atom:ivector"); } else { index = ndvector; ndvector++; dname = (char **) memory->srealloc(dname,ndvector*sizeof(char *), "atom:dname"); int n = strlen(name) + 1; dname[index] = new char[n]; strcpy(dname[index],name); dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *), "atom:dvector"); memory->create(dvector[index],nmax,"atom:dvector"); } return index; } /* ---------------------------------------------------------------------- remove a custom variable of type flag = 0/1 for int/double at index free memory for vector and name and set ptrs to NULL ivector/dvector and iname/dname lists never shrink ------------------------------------------------------------------------- */ void Atom::remove_custom(int flag, int index) { if (flag == 0) { memory->destroy(ivector[index]); ivector[index] = NULL; delete [] iname[index]; iname[index] = NULL; } else { memory->destroy(dvector[index]); dvector[index] = NULL; delete [] dname[index]; dname[index] = NULL; } } /* ---------------------------------------------------------------------- return a pointer to a named internal variable if don't recognize name, return NULL customize by adding names ------------------------------------------------------------------------- */ void *Atom::extract(char *name) { if (strcmp(name,"mass") == 0) return (void *) mass; if (strcmp(name,"id") == 0) return (void *) tag; if (strcmp(name,"type") == 0) return (void *) type; if (strcmp(name,"mask") == 0) return (void *) mask; if (strcmp(name,"image") == 0) return (void *) image; if (strcmp(name,"x") == 0) return (void *) x; if (strcmp(name,"v") == 0) return (void *) v; if (strcmp(name,"f") == 0) return (void *) f; if (strcmp(name,"molecule") == 0) return (void *) molecule; if (strcmp(name,"q") == 0) return (void *) q; if (strcmp(name,"mu") == 0) return (void *) mu; if (strcmp(name,"omega") == 0) return (void *) omega; if (strcmp(name,"angmom") == 0) return (void *) angmom; if (strcmp(name,"torque") == 0) return (void *) torque; if (strcmp(name,"radius") == 0) return (void *) radius; if (strcmp(name,"rmass") == 0) return (void *) rmass; if (strcmp(name,"ellipsoid") == 0) return (void *) ellipsoid; if (strcmp(name,"line") == 0) return (void *) line; if (strcmp(name,"tri") == 0) return (void *) tri; if (strcmp(name,"vfrac") == 0) return (void *) vfrac; if (strcmp(name,"s0") == 0) return (void *) s0; if (strcmp(name,"x0") == 0) return (void *) x0; if (strcmp(name,"spin") == 0) return (void *) spin; if (strcmp(name,"eradius") == 0) return (void *) eradius; if (strcmp(name,"ervel") == 0) return (void *) ervel; if (strcmp(name,"erforce") == 0) return (void *) erforce; if (strcmp(name,"ervelforce") == 0) return (void *) ervelforce; if (strcmp(name,"cs") == 0) return (void *) cs; if (strcmp(name,"csforce") == 0) return (void *) csforce; if (strcmp(name,"vforce") == 0) return (void *) vforce; if (strcmp(name,"etag") == 0) return (void *) etag; if (strcmp(name,"rho") == 0) return (void *) rho; if (strcmp(name,"drho") == 0) return (void *) drho; if (strcmp(name,"e") == 0) return (void *) e; if (strcmp(name,"de") == 0) return (void *) de; if (strcmp(name,"cv") == 0) return (void *) cv; if (strcmp(name,"vest") == 0) return (void *) vest; if (strcmp(name, "contact_radius") == 0) return (void *) contact_radius; if (strcmp(name, "smd_data_9") == 0) return (void *) smd_data_9; if (strcmp(name, "smd_stress") == 0) return (void *) smd_stress; if (strcmp(name, "eff_plastic_strain") == 0) return (void *) eff_plastic_strain; if (strcmp(name, "eff_plastic_strain_rate") == 0) return (void *) eff_plastic_strain_rate; if (strcmp(name, "damage") == 0) return (void *) damage; if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta; return NULL; } /* ---------------------------------------------------------------------- return # of bytes of allocated memory call to avec tallies per-atom vectors add in global to local mapping storage ------------------------------------------------------------------------- */ bigint Atom::memory_usage() { memlength = DELTA_MEMSTR; memory->create(memstr,memlength,"atom:memstr"); memstr[0] = '\0'; bigint bytes = avec->memory_usage(); memory->destroy(memstr); bytes += max_same*sizeof(int); if (map_style == 1) bytes += memory->usage(map_array,map_maxarray); else if (map_style == 2) { bytes += map_nbucket*sizeof(int); bytes += map_nhash*sizeof(HashElem); } if (maxnext) { bytes += memory->usage(next,maxnext); bytes += memory->usage(permute,maxnext); } return bytes; } /* ---------------------------------------------------------------------- accumulate per-atom vec names in memstr, padded by spaces return 1 if padded str is not already in memlist, else 0 ------------------------------------------------------------------------- */ int Atom::memcheck(const char *str) { int n = strlen(str) + 3; char *padded = new char[n]; strcpy(padded," "); strcat(padded,str); strcat(padded," "); if (strstr(memstr,padded)) { delete [] padded; return 0; } if (strlen(memstr) + n >= memlength) { memlength += DELTA_MEMSTR; memory->grow(memstr,memlength,"atom:memstr"); } strcat(memstr,padded); delete [] padded; return 1; } diff --git a/src/molecule.cpp b/src/molecule.cpp index 76b28e3d4..e0e9ec8aa 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -1,1813 +1,1813 @@ /* ---------------------------------------------------------------------- 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 #include #include "molecule.h" #include "atom.h" #include "atom_vec.h" #include "atom_vec_body.h" #include "force.h" #include "comm.h" #include "domain.h" #include "math_extra.h" #include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; #define MAXLINE 256 #define EPSILON 1.0e-7 #define BIG 1.0e20 #define SINERTIA 0.4 // moment of inertia prefactor for sphere /* ---------------------------------------------------------------------- */ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : Pointers(lmp), id(NULL), 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(NULL), angle_atom2(NULL), angle_atom3(NULL), num_dihedral(NULL), dihedral_type(NULL), dihedral_atom1(NULL), dihedral_atom2(NULL), dihedral_atom3(NULL), dihedral_atom4(NULL), num_improper(NULL), improper_type(NULL), improper_atom1(NULL), improper_atom2(NULL), improper_atom3(NULL), improper_atom4(NULL), nspecial(NULL), special(NULL), shake_flag(NULL), shake_atom(NULL), shake_type(NULL), avec_body(NULL), ibodyparams(NULL), dbodyparams(NULL), dx(NULL), dxcom(NULL), dxbody(NULL), quat_external(NULL), fp(NULL), count(NULL) { me = comm->me; if (index >= narg) 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 template ID must be " "alphanumeric or underscore characters"); // parse args until reach unknown arg (next file) toffset = 0; boffset = aoffset = doffset = ioffset = 0; sizescale = 1.0; int ifile = index; int iarg = ifile+1; while (iarg < narg) { if (strcmp(arg[iarg],"offset") == 0) { if (iarg+6 > narg) error->all(FLERR,"Illegal molecule command"); toffset = force->inumeric(FLERR,arg[iarg+1]); boffset = force->inumeric(FLERR,arg[iarg+2]); aoffset = force->inumeric(FLERR,arg[iarg+3]); doffset = force->inumeric(FLERR,arg[iarg+4]); ioffset = force->inumeric(FLERR,arg[iarg+5]); if (toffset < 0 || boffset < 0 || aoffset < 0 || doffset < 0 || ioffset < 0) error->all(FLERR,"Illegal molecule command"); iarg += 6; } else if (strcmp(arg[iarg],"toff") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command"); toffset = force->inumeric(FLERR,arg[iarg+1]); if (toffset < 0) error->all(FLERR,"Illegal molecule command"); iarg += 2; } else if (strcmp(arg[iarg],"boff") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command"); boffset = force->inumeric(FLERR,arg[iarg+1]); if (boffset < 0) error->all(FLERR,"Illegal molecule command"); iarg += 2; } else if (strcmp(arg[iarg],"aoff") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command"); aoffset = force->inumeric(FLERR,arg[iarg+1]); if (aoffset < 0) error->all(FLERR,"Illegal molecule command"); iarg += 2; } else if (strcmp(arg[iarg],"doff") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command"); doffset = force->inumeric(FLERR,arg[iarg+1]); if (doffset < 0) error->all(FLERR,"Illegal molecule command"); iarg += 2; } else if (strcmp(arg[iarg],"ioff") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command"); ioffset = force->inumeric(FLERR,arg[iarg+1]); if (ioffset < 0) error->all(FLERR,"Illegal molecule command"); iarg += 2; } else if (strcmp(arg[iarg],"scale") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command"); sizescale = force->numeric(FLERR,arg[iarg+1]); if (sizescale <= 0.0) error->all(FLERR,"Illegal molecule command"); iarg += 2; } else break; } index = iarg; // last molecule if have scanned all args if (iarg == narg) last = 1; else last = 0; // initialize all fields to empty initialize(); // scan file for sizes of all fields and allocate them if (me == 0) open(arg[ifile]); read(0); if (me == 0) fclose(fp); allocate(); // read file again to populate all fields if (me == 0) open(arg[ifile]); read(1); if (me == 0) fclose(fp); // stats if (me == 0) { if (screen) fprintf(screen,"Read molecule %s:\n" " %d atoms with %d types\n %d bonds with %d types\n" " %d angles with %d types\n %d dihedrals with %d types\n" " %d impropers with %d types\n", id,natoms,ntypes, nbonds,nbondtypes,nangles,nangletypes, ndihedrals,ndihedraltypes,nimpropers,nimpropertypes); if (logfile) fprintf(logfile,"Read molecule %s:\n" " %d atoms with %d types\n %d bonds with %d types\n" " %d angles with %d types\n %d dihedrals with %d types\n" " %d impropers with %d types\n", id,natoms,ntypes, nbonds,nbondtypes,nangles,nangletypes, ndihedrals,ndihedraltypes,nimpropers,nimpropertypes); } } /* ---------------------------------------------------------------------- */ Molecule::~Molecule() { delete [] id; deallocate(); } /* ---------------------------------------------------------------------- compute center = geometric center of molecule also compute: dx = displacement of each atom from center molradius = radius of molecule from center including finite-size particles or body particles ------------------------------------------------------------------------- */ 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]; } molradius = 0.0; for (int i = 0; i < natoms; i++) { double rad = MathExtra::len3(dx[i]); if (radiusflag) rad += radius[i]; molradius = MAX(molradius,rad); } } /* ---------------------------------------------------------------------- compute masstotal = total mass of molecule could have been set by user, otherwise calculate it ------------------------------------------------------------------------- */ void Molecule::compute_mass() { if (massflag) return; massflag = 1; - if (!rmassflag) atom->check_mass(FLERR); + atom->check_mass(FLERR); masstotal = 0.0; for (int i = 0; i < natoms; i++) { if (rmassflag) masstotal += rmass[i]; else masstotal += atom->mass[type[i]]; } } /* ---------------------------------------------------------------------- compute com = center of mass of molecule could have been set by user, otherwise calculate it works for finite size particles assuming no overlap also compute: dxcom = displacement of each atom from COM comatom = which atom (1-Natom) is nearest the COM maxextent = furthest any atom in molecule is from comatom (not COM) ------------------------------------------------------------------------- */ void Molecule::compute_com() { if (!comflag) { comflag = 1; - if (!rmassflag) atom->check_mass(FLERR); + atom->check_mass(FLERR); double onemass; com[0] = com[1] = com[2] = 0.0; for (int i = 0; i < natoms; i++) { if (rmassflag) onemass = rmass[i]; else onemass = atom->mass[type[i]]; com[0] += x[i][0]*onemass; com[1] += x[i][1]*onemass; com[2] += x[i][2]*onemass; } if (masstotal > 0.0) { com[0] /= masstotal; com[1] /= masstotal; com[2] /= masstotal; } } memory->destroy(dxcom); memory->create(dxcom,natoms,3,"molecule:dxcom"); for (int i = 0; i < natoms; i++) { dxcom[i][0] = x[i][0] - com[0]; dxcom[i][1] = x[i][1] - com[1]; dxcom[i][2] = x[i][2] - com[2]; } double rsqmin = BIG; for (int i = 0; i < natoms; i++) { double rsq = MathExtra::lensq3(dxcom[i]); if (rsq < rsqmin) { comatom = i; rsqmin = rsq; } } double rsqmax = 0.0; for (int i = 0; i < natoms; i++) { double dx = x[comatom][0] - x[i][0]; double dy = x[comatom][1] - x[i][1]; double dz = x[comatom][2] - x[i][2]; double rsq = dx*dx + dy*dy + dz*dz; rsqmax = MAX(rsqmax,rsq); } comatom++; maxextent = sqrt(rsqmax); } /* ---------------------------------------------------------------------- compute itensor = 6 moments of inertia of molecule around xyz axes could have been set by user, otherwise calculate it accounts for finite size spheres, assuming no overlap also compute: inertia = 3 principal components of inertia ex,ey,ez = principal axes in space coords quat = quaternion for orientation of molecule dxbody = displacement of each atom from COM in body frame ------------------------------------------------------------------------- */ void Molecule::compute_inertia() { if (!inertiaflag) { inertiaflag = 1; - if (!rmassflag) atom->check_mass(FLERR); + atom->check_mass(FLERR); double onemass,dx,dy,dz; for (int i = 0; i < 6; i++) itensor[i] = 0.0; for (int i = 0; i < natoms; i++) { if (rmassflag) onemass = rmass[i]; else onemass = atom->mass[type[i]]; dx = dxcom[i][0]; dy = dxcom[i][1]; dz = dxcom[i][2]; itensor[0] += onemass * (dy*dy + dz*dz); itensor[1] += onemass * (dx*dx + dz*dz); itensor[2] += onemass * (dx*dx + dy*dy); itensor[3] -= onemass * dy*dz; itensor[4] -= onemass * dx*dz; itensor[5] -= onemass * dx*dy; } if (radiusflag) { for (int i = 0; i < natoms; i++) { if (rmassflag) onemass = rmass[i]; else onemass = atom->mass[type[i]]; itensor[0] += SINERTIA*onemass * radius[i]*radius[i]; itensor[1] += SINERTIA*onemass * radius[i]*radius[i]; itensor[2] += SINERTIA*onemass * radius[i]*radius[i]; } } } // diagonalize inertia tensor for each body via Jacobi rotations // inertia = 3 eigenvalues = principal moments of inertia // evectors and exzy = 3 evectors = principal axes of rigid body double cross[3]; double tensor[3][3],evectors[3][3]; tensor[0][0] = itensor[0]; tensor[1][1] = itensor[1]; tensor[2][2] = itensor[2]; tensor[1][2] = tensor[2][1] = itensor[3]; tensor[0][2] = tensor[2][0] = itensor[4]; tensor[0][1] = tensor[1][0] = itensor[5]; if (MathExtra::jacobi(tensor,inertia,evectors)) error->all(FLERR,"Insufficient Jacobi rotations for rigid molecule"); ex[0] = evectors[0][0]; ex[1] = evectors[1][0]; ex[2] = evectors[2][0]; ey[0] = evectors[0][1]; ey[1] = evectors[1][1]; ey[2] = evectors[2][1]; ez[0] = evectors[0][2]; ez[1] = evectors[1][2]; ez[2] = evectors[2][2]; // if any principal moment < scaled EPSILON, set to 0.0 double max; max = MAX(inertia[0],inertia[1]); max = MAX(max,inertia[2]); if (inertia[0] < EPSILON*max) inertia[0] = 0.0; if (inertia[1] < EPSILON*max) inertia[1] = 0.0; if (inertia[2] < EPSILON*max) inertia[2] = 0.0; // enforce 3 evectors as a right-handed coordinate system // flip 3rd vector if needed MathExtra::cross3(ex,ey,cross); if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez); // create quaternion MathExtra::exyz_to_q(ex,ey,ez,quat); // compute displacements in body frame defined by quat memory->destroy(dxbody); memory->create(dxbody,natoms,3,"molecule:dxbody"); for (int i = 0; i < natoms; i++) MathExtra::transpose_matvec(ex,ey,ez,dxcom[i],dxbody[i]); } /* ---------------------------------------------------------------------- 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 if (strstr(line,"mass")) { massflag = 1; sscanf(line,"%lg",&masstotal); masstotal *= sizescale*sizescale*sizescale; } else if (strstr(line,"com")) { comflag = 1; sscanf(line,"%lg %lg %lg",&com[0],&com[1],&com[2]); com[0] *= sizescale; com[1] *= sizescale; com[2] *= sizescale; if (domain->dimension == 2 && com[2] != 0.0) error->all(FLERR,"Molecule file z center-of-mass must be 0.0 for 2d"); } else if (strstr(line,"inertia")) { inertiaflag = 1; sscanf(line,"%lg %lg %lg %lg %lg %lg", &itensor[0],&itensor[1],&itensor[2], &itensor[3],&itensor[4],&itensor[5]); itensor[0] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[1] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[2] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[3] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[4] *= sizescale*sizescale*sizescale*sizescale*sizescale; itensor[5] *= sizescale*sizescale*sizescale*sizescale*sizescale; } else if (strstr(line,"body")) { bodyflag = 1; avec_body = (AtomVecBody *) atom->style_match("body"); if (!avec_body) error->all(FLERR,"Molecule file requires atom style body"); sscanf(line,"%d %d",&nibody,&ndbody); } else break; } // error checks if (natoms < 1) error->all(FLERR,"No count or invalid atom count in molecule file"); if (nbonds < 0) error->all(FLERR,"Invalid bond count in molecule file"); if (nangles < 0) error->all(FLERR,"Invalid angle count in molecule file"); if (ndihedrals < 0) error->all(FLERR,"Invalid dihedral count in molecule file"); if (nimpropers < 0) error->all(FLERR,"Invalid improper count in molecule file"); // 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 = tag_require = 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 = tag_require = 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 = tag_require = 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 = tag_require = 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 = tag_require = 1; if (flag) special_read(line); else skip_lines(natoms,line); } else if (strcmp(keyword,"Shake Flags") == 0) { shakeflagflag = 1; if (flag) shakeflag_read(line); else skip_lines(natoms,line); } else if (strcmp(keyword,"Shake Atoms") == 0) { shakeatomflag = tag_require = 1; if (shaketypeflag) shakeflag = 1; if (!shakeflagflag) error->all(FLERR,"Molecule file shake flags not before shake atoms"); if (flag) shakeatom_read(line); else skip_lines(natoms,line); } else if (strcmp(keyword,"Shake Bond Types") == 0) { shaketypeflag = 1; if (shakeatomflag) shakeflag = 1; if (!shakeflagflag) error->all(FLERR,"Molecule file shake flags not before shake bonds"); if (flag) shaketype_read(line); else skip_lines(natoms,line); } else if (strcmp(keyword,"Body Integers") == 0) { if (bodyflag == 0 || nibody == 0) error->all(FLERR,"Molecule file has body params " "but no setting for them"); ibodyflag = 1; body(flag,0,line); } else if (strcmp(keyword,"Body Doubles") == 0) { if (bodyflag == 0 || ndbody == 0) error->all(FLERR,"Molecule file has body params " "but no setting for them"); dbodyflag = 1; body(flag,1,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 ((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 ((shakeflagflag || shakeatomflag || shaketypeflag) && !shakeflag) error->all(FLERR,"Molecule file shake info is incomplete"); if (bodyflag && nibody && ibodyflag == 0) error->all(FLERR,"Molecule file has no Body Integers section"); if (bodyflag && ndbody && dbodyflag == 0) error->all(FLERR,"Molecule file has no Body Doubles section"); } // auto-generate special bonds if needed and not in file // set maxspecial on first pass, so allocate() has a size if (bondflag && specialflag == 0) { if (domain->box_exist == 0) error->all(FLERR,"Cannot auto-generate special bonds before " "simulation box is defined"); maxspecial = atom->maxspecial; if (flag) { special_generate(); specialflag = 1; nspecialflag = 1; } } // body particle must have natom = 1 // set radius by having body class compute its own radius if (bodyflag) { radiusflag = 1; if (natoms != 1) error->all(FLERR,"Molecule natoms must be 1 for body particle"); if (sizescale != 1.0) error->all(FLERR,"Molecule sizescale must be 1.0 for body particle"); if (flag) { radius[0] = avec_body->radius_body(nibody,ndbody,ibodyparams,dbodyparams); maxradius = radius[0]; } } } /* ---------------------------------------------------------------------- read coords from file ------------------------------------------------------------------------- */ void Molecule::coords(char *line) { int tmp; for (int i = 0; i < natoms; i++) { readline(line); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 4) error->all(FLERR,"Invalid Coords section in molecule file"); } sscanf(line,"%d %lg %lg %lg",&tmp,&x[i][0],&x[i][1],&x[i][2]); x[i][0] *= sizescale; x[i][1] *= sizescale; x[i][2] *= sizescale; } 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 set ntypes = max of any atom type ------------------------------------------------------------------------- */ void Molecule::types(char *line) { int tmp; for (int i = 0; i < natoms; i++) { readline(line); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 2) error->all(FLERR,"Invalid Types section in molecule file"); } sscanf(line,"%d %d",&tmp,&type[i]); type[i] += toffset; } for (int i = 0; i < natoms; i++) if (type[i] <= 0) error->all(FLERR,"Invalid atom type in molecule file"); for (int i = 0; i < natoms; i++) ntypes = MAX(ntypes,type[i]); } /* ---------------------------------------------------------------------- read charges from file ------------------------------------------------------------------------- */ void Molecule::charges(char *line) { int tmp; for (int i = 0; i < natoms; i++) { readline(line); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 2) error->all(FLERR,"Invalid Charges section in molecule file"); } sscanf(line,"%d %lg",&tmp,&q[i]); } } /* ---------------------------------------------------------------------- read diameters from file and set radii ------------------------------------------------------------------------- */ void Molecule::diameters(char *line) { int tmp; maxradius = 0.0; for (int i = 0; i < natoms; i++) { readline(line); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 2) error->all(FLERR,"Invalid Diameters section in molecule file"); } sscanf(line,"%d %lg",&tmp,&radius[i]); radius[i] *= sizescale; radius[i] *= 0.5; maxradius = MAX(maxradius,radius[i]); } for (int i = 0; i < natoms; i++) if (radius[i] < 0.0) error->all(FLERR,"Invalid atom diameter in molecule file"); } /* ---------------------------------------------------------------------- read masses from file ------------------------------------------------------------------------- */ void Molecule::masses(char *line) { int tmp; for (int i = 0; i < natoms; i++) { readline(line); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 2) error->all(FLERR,"Invalid Masses section in molecule file"); } sscanf(line,"%d %lg",&tmp,&rmass[i]); rmass[i] *= sizescale*sizescale*sizescale; } for (int i = 0; i < natoms; i++) if (rmass[i] <= 0.0) error->all(FLERR,"Invalid atom mass in molecule file"); } /* ---------------------------------------------------------------------- read bonds from file set nbondtypes = max type of any bond 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 tmp,itype; tagint m,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); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 4) error->all(FLERR,"Invalid Bonds section in molecule file"); } sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2); itype += boffset; if (atom1 <= 0 || atom1 > natoms || atom2 <= 0 || atom2 > natoms) error->one(FLERR,"Invalid atom ID in Bonds section of molecule file"); if (itype <= 0) error->one(FLERR,"Invalid bond type in Bonds section of molecule file"); if (flag) { m = atom1-1; nbondtypes = MAX(nbondtypes,itype); 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 tmp,itype; tagint m,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); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 5) error->all(FLERR,"Invalid Angles section in molecule file"); } sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3); itype += aoffset; 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) error->one(FLERR,"Invalid angle type in Angles section of molecule file"); if (flag) { m = atom2-1; nangletypes = MAX(nangletypes,itype); 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 tmp,itype; tagint m,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); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 6) error->all(FLERR,"Invalid Dihedrals section in molecule file"); } sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " ", &tmp,&itype,&atom1,&atom2,&atom3,&atom4); itype += doffset; 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) error->one(FLERR, "Invalid dihedral type in dihedrals section of molecule file"); if (flag) { m = atom2-1; ndihedraltypes = MAX(ndihedraltypes,itype); 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; dihedral_atom4[m][num_dihedral[m]] = atom4; 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; dihedral_atom4[m][num_dihedral[m]] = atom4; 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; dihedral_atom4[m][num_dihedral[m]] = atom4; 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; dihedral_atom4[m][num_dihedral[m]] = atom4; 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 tmp,itype; tagint m,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); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 6) error->all(FLERR,"Invalid Impropers section in molecule file"); } sscanf(line,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " ", &tmp,&itype,&atom1,&atom2,&atom3,&atom4); itype += ioffset; 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) error->one(FLERR, "Invalid improper type in impropers section of molecule file"); if (flag) { m = atom2-1; nimpropertypes = MAX(nimpropertypes,itype); 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; improper_atom4[m][num_improper[m]] = atom4; 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; improper_atom4[m][num_improper[m]] = atom4; 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; improper_atom4[m][num_improper[m]] = atom4; 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; improper_atom4[m][num_improper[m]] = atom4; 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); if (i == 0) { int nwords = atom->count_words(line); if (nwords != 4) error->all(FLERR,"Invalid Special Bond Counts section in " "molecule file"); } 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] = ATOTAGINT(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; } /* ---------------------------------------------------------------------- auto generate special bond info ------------------------------------------------------------------------- */ void Molecule::special_generate() { int newton_bond = force->newton_bond; tagint atom1,atom2; int count[natoms]; for (int i = 0; i < natoms; i++) count[i] = 0; // 1-2 neighbors if (newton_bond) { for (int i = 0; i < natoms; i++) { for (int j = 0; j < num_bond[i]; j++) { atom1 = i; atom2 = bond_atom[i][j]-1; nspecial[i][0]++; nspecial[atom2][0]++; if (count[i] >= maxspecial || count[atom2] >= maxspecial) error->one(FLERR,"Molecule auto special bond generation overflow"); special[i][count[i]++] = atom2 + 1; special[atom2][count[atom2]++] = i + 1; } } } else { for (int i = 0; i < natoms; i++) { nspecial[i][0] = num_bond[i]; for (int j = 0; j < num_bond[i]; j++) { atom1 = i; atom2 = bond_atom[i][j]; if (count[atom1] >= maxspecial) error->one(FLERR,"Molecule auto special bond generation overflow"); special[i][count[atom1]++] = atom2; } } } // 1-3 neighbors with no duplicates for (int i = 0; i < natoms; i++) nspecial[i][1] = nspecial[i][0]; int dedup; for (int i = 0; i < natoms; i++) { for (int m = 0; m < nspecial[i][0]; m++) { for (int j = 0; j < nspecial[special[i][m]-1][0]; j++) { dedup = 0; for (int k =0; k < count[i]; k++) { if (special[special[i][m]-1][j] == special[i][k] || special[special[i][m]-1][j] == i+1) { dedup = 1; } } if (!dedup) { if (count[i] >= maxspecial) error->one(FLERR,"Molecule auto special bond generation overflow"); special[i][count[i]++] = special[special[i][m]-1][j]; nspecial[i][1]++; } } } } // 1-4 neighbors with no duplicates for (int i = 0; i < natoms; i++) nspecial[i][2] = nspecial[i][1]; for (int i = 0; i < natoms; i++) { for (int m = nspecial[i][0]; m < nspecial[i][1]; m++) { for (int j = 0; j < nspecial[special[i][m]-1][0]; j++) { dedup = 0; for (int k =0; k < count[i]; k++) { if (special[special[i][m]-1][j] == special[i][k] || special[special[i][m]-1][j] == i+1) { dedup = 1; } } if (!dedup) { if (count[i] >= maxspecial) error->one(FLERR,"Molecule auto special bond generation overflow"); special[i][count[i]++] = special[special[i][m]-1][j]; nspecial[i][2]++; } } } } } /* ---------------------------------------------------------------------- read SHAKE flags from file ------------------------------------------------------------------------- */ void Molecule::shakeflag_read(char *line) { int tmp; for (int i = 0; i < natoms; i++) { readline(line); sscanf(line,"%d %d",&tmp,&shake_flag[i]); } for (int i = 0; i < natoms; i++) if (shake_flag[i] < 0 || shake_flag[i] > 4) error->all(FLERR,"Invalid shake flag in molecule file"); } /* ---------------------------------------------------------------------- read SHAKE atom info from file ------------------------------------------------------------------------- */ void Molecule::shakeatom_read(char *line) { int tmp; for (int i = 0; i < natoms; i++) { readline(line); if (shake_flag[i] == 1) sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&shake_atom[i][0],&shake_atom[i][1],&shake_atom[i][2]); else if (shake_flag[i] == 2) sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&shake_atom[i][0],&shake_atom[i][1]); else if (shake_flag[i] == 3) sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&shake_atom[i][0],&shake_atom[i][1],&shake_atom[i][2]); else if (shake_flag[i] == 4) sscanf(line,"%d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&shake_atom[i][0],&shake_atom[i][1], &shake_atom[i][2],&shake_atom[i][3]); } for (int i = 0; i < natoms; i++) { int m = shake_flag[i]; if (m == 1) m = 3; for (int j = 0; j < m; j++) if (shake_atom[i][j] <= 0 || shake_atom[i][j] > natoms) error->all(FLERR,"Invalid shake atom in molecule file"); } } /* ---------------------------------------------------------------------- read SHAKE bond type info from file ------------------------------------------------------------------------- */ void Molecule::shaketype_read(char *line) { int tmp; for (int i = 0; i < natoms; i++) { readline(line); if (shake_flag[i] == 1) sscanf(line,"%d %d %d %d",&tmp, &shake_type[i][0],&shake_type[i][1],&shake_type[i][2]); else if (shake_flag[i] == 2) sscanf(line,"%d %d",&tmp,&shake_type[i][0]); else if (shake_flag[i] == 3) sscanf(line,"%d %d %d",&tmp,&shake_type[i][0],&shake_type[i][1]); else if (shake_flag[i] == 4) sscanf(line,"%d %d %d %d",&tmp, &shake_type[i][0],&shake_type[i][1],&shake_type[i][2]); } for (int i = 0; i < natoms; i++) { int m = shake_flag[i]; if (m == 1) m = 3; for (int j = 0; j < m-1; j++) if (shake_type[i][j] <= 0) error->all(FLERR,"Invalid shake bond type in molecule file"); if (shake_flag[i] == 1) if (shake_type[i][2] <= 0) error->all(FLERR,"Invalid shake angle type in molecule file"); } } /* ---------------------------------------------------------------------- read body params from file pflag = 0/1 for integer/double params ------------------------------------------------------------------------- */ void Molecule::body(int flag, int pflag, char *line) { int i,ncount; int nparam = nibody; if (pflag) nparam = ndbody; int nword = 0; while (nword < nparam) { readline(line); ncount = atom->count_words(line); if (ncount == 0) error->one(FLERR,"Too few values in body section of molecule file"); if (nword+ncount > nparam) error->all(FLERR,"Too many values in body section of molecule file"); if (flag) { if (pflag == 0) { ibodyparams[nword++] = force->inumeric(FLERR,strtok(line," \t\n\r\f")); for (i = 1; i < ncount; i++) ibodyparams[nword++] = force->inumeric(FLERR,strtok(NULL," \t\n\r\f")); } else { dbodyparams[nword++] = force->numeric(FLERR,strtok(line," \t\n\r\f")); for (i = 1; i < ncount; i++) dbodyparams[nword++] = force->numeric(FLERR,strtok(NULL," \t\n\r\f")); } } else nword += ncount; } } /* ---------------------------------------------------------------------- error check molecule attributes and topology against system settings flag = 0, just check this molecule flag = 1, check all molecules in set, this is 1st molecule in set ------------------------------------------------------------------------- */ void Molecule::check_attributes(int flag) { int n = 1; if (flag) n = nset; int imol = atom->find_molecule(id); for (int i = imol; i < imol+n; i++) { Molecule *onemol = atom->molecules[imol]; // check per-atom attributes of molecule // warn if not a match int mismatch = 0; if (onemol->qflag && !atom->q_flag) mismatch = 1; if (onemol->radiusflag && !atom->radius_flag) mismatch = 1; if (onemol->rmassflag && !atom->rmass_flag) mismatch = 1; if (mismatch && me == 0) error->warning(FLERR, "Molecule attributes do not match system attributes"); // for all atom styles, check nbondtype,etc mismatch = 0; if (atom->nbondtypes < onemol->nbondtypes) mismatch = 1; if (atom->nangletypes < onemol->nangletypes) mismatch = 1; if (atom->ndihedraltypes < onemol->ndihedraltypes) mismatch = 1; if (atom->nimpropertypes < onemol->nimpropertypes) mismatch = 1; if (mismatch) error->all(FLERR,"Molecule topology type exceeds system topology type"); // for molecular atom styles, check bond_per_atom,etc + maxspecial // do not check for atom style template, since nothing stored per atom if (atom->molecular == 1) { if (atom->avec->bonds_allow && atom->bond_per_atom < onemol->bond_per_atom) mismatch = 1; if (atom->avec->angles_allow && atom->angle_per_atom < onemol->angle_per_atom) mismatch = 1; if (atom->avec->dihedrals_allow && atom->dihedral_per_atom < onemol->dihedral_per_atom) mismatch = 1; if (atom->avec->impropers_allow && atom->improper_per_atom < onemol->improper_per_atom) mismatch = 1; if (atom->maxspecial < onemol->maxspecial) mismatch = 1; if (mismatch) error->all(FLERR,"Molecule topology/atom exceeds system topology/atom"); } // warn if molecule topology defined but no special settings if (onemol->bondflag && !onemol->specialflag) if (me == 0) error->warning(FLERR,"Molecule has bond topology " "but no special bond settings"); } } /* ---------------------------------------------------------------------- init all data structures to empty ------------------------------------------------------------------------- */ void Molecule::initialize() { natoms = 0; nbonds = nangles = ndihedrals = nimpropers = 0; ntypes = 0; nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; nibody = ndbody = 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; shakeflag = shakeflagflag = shakeatomflag = shaketypeflag = 0; bodyflag = ibodyflag = dbodyflag = 0; centerflag = massflag = comflag = inertiaflag = 0; tag_require = 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; shake_flag = NULL; shake_atom = NULL; shake_type = NULL; ibodyparams = NULL; dbodyparams = NULL; dx = NULL; dxcom = NULL; dxbody = NULL; } /* ---------------------------------------------------------------------- allocate all data structures also initialize values for data structures that are always allocated ------------------------------------------------------------------------- */ 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"); // always allocate num_bond,num_angle,etc and special+nspecial // even if not in molecule file, initialize to 0 // this is so methods that use these arrays don't have to check they exist memory->create(num_bond,natoms,"molecule:num_bond"); for (int i = 0; i < natoms; i++) num_bond[i] = 0; memory->create(num_angle,natoms,"molecule:num_angle"); for (int i = 0; i < natoms; i++) num_angle[i] = 0; memory->create(num_dihedral,natoms,"molecule:num_dihedral"); for (int i = 0; i < natoms; i++) num_dihedral[i] = 0; memory->create(num_improper,natoms,"molecule:num_improper"); for (int i = 0; i < natoms; i++) num_improper[i] = 0; memory->create(special,natoms,maxspecial,"molecule:special"); memory->create(nspecial,natoms,3,"molecule:nspecial"); for (int i = 0; i < natoms; i++) nspecial[i][0] = nspecial[i][1] = nspecial[i][2] = 0; if (bondflag) { 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(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(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(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 (shakeflag) { memory->create(shake_flag,natoms,"molecule:shake_flag"); memory->create(shake_atom,natoms,4,"molecule:shake_flag"); memory->create(shake_type,natoms,3,"molecule:shake_flag"); } if (bodyflag) { if (nibody) memory->create(ibodyparams,nibody,"molecule:ibodyparams"); if (ndbody) memory->create(dbodyparams,ndbody,"molecule:dbodyparams"); } } /* ---------------------------------------------------------------------- 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(shake_flag); memory->destroy(shake_atom); memory->destroy(shake_type); memory->destroy(dx); memory->destroy(dxcom); memory->destroy(dxbody); memory->destroy(ibodyparams); memory->destroy(dbodyparams); } /* ---------------------------------------------------------------------- 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"); } } } */