diff --git a/src/pair.h b/src/pair.h index 40f78934e..8eb6854c0 100644 --- a/src/pair.h +++ b/src/pair.h @@ -1,296 +1,297 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifndef LMP_PAIR_H #define LMP_PAIR_H #include "pointers.h" namespace LAMMPS_NS { class Pair : protected Pointers { friend class AngleSDK; friend class AngleSDKOMP; friend class BondQuartic; friend class BondQuarticOMP; friend class DihedralCharmm; friend class DihedralCharmmOMP; friend class FixGPU; friend class FixOMP; friend class ThrOMP; public: double eng_vdwl,eng_coul; // accumulated energies double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial double cutforce; // max cutoff for all atom pairs double **cutsq; // cutoff sq for each atom pair int **setflag; // 0/1 = whether each i,j has been set int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int comm_reverse_off; // size of reverse comm even if newton off int single_enable; // 1 if single() routine exists int restartinfo; // 1 if pair style writes restart info int respa_enable; // 1 if inner/middle/outer rRESPA routines int one_coeff; // 1 if allows only one coeff * * call int no_virial_fdotr_compute; // 1 if does not invoke virial_fdotr_compute() int writedata; // 1 if writes coeffs to data file int ghostneigh; // 1 if pair style needs neighbors of ghosts double **cutghost; // cutoff for each ghost pair int ewaldflag; // 1 if compatible with Ewald solver int pppmflag; // 1 if compatible with PPPM solver int msmflag; // 1 if compatible with MSM solver int dispersionflag; // 1 if compatible with LJ/dispersion solver int tip4pflag; // 1 if compatible with TIP4P solver int tail_flag; // pair_modify flag for LJ tail correction double etail,ptail; // energy/pressure tail corrections double etail_ij,ptail_ij; int evflag; // energy,virial settings int eflag_either,eflag_global,eflag_atom; int vflag_either,vflag_global,vflag_atom; int ncoultablebits; // size of Coulomb table, accessed by KSpace int ndisptablebits; // size of dispersion table double tabinnersq; double tabinnerdispsq; double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable; double *etable,*detable,*ptable,*dptable,*vtable,*dvtable; double *rdisptable, *drdisptable, *fdisptable, *dfdisptable; double *edisptable, *dedisptable; int ncoulshiftbits,ncoulmask; int ndispshiftbits, ndispmask; int nextra; // # of extra quantities pair style calculates double *pvector; // vector of extra pair quantities int single_extra; // number of extra single values calculated double *svector; // vector of extra single quantities class NeighList *list; // standard neighbor list used by most pairs class NeighList *listhalf; // half list used by some pairs class NeighList *listfull; // full list used by some pairs class NeighList *listgranhistory; // granular history list used by some pairs class NeighList *listinner; // rRESPA lists used by some pairs class NeighList *listmiddle; class NeighList *listouter; unsigned int datamask; unsigned int datamask_ext; int compute_flag; // 0 if skip compute() Pair(class LAMMPS *); virtual ~Pair(); // top-level Pair methods void init(); void reinit(); double mix_energy(double, double, double, double); double mix_distance(double, double); void write_file(int, char **); void init_bitmap(double, double, int, int &, int &, int &, int &); virtual void modify_params(int, char **); void compute_dummy(int, int); // need to be public, so can be called by pair_style reaxc void v_tally(int, double *, double *); void ev_tally(int, int, int, int, double, double, double, double, double, double); void ev_tally3(int, int, int, double, double, double *, double *, double *, double *); void v_tally3(int, int, int, double *, double *, double *, double *); void v_tally4(int, int, int, int, double *, double *, double *, double *, double *, double *); void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double, double); // general child-class methods virtual void compute(int, int) = 0; virtual void compute_inner() {} virtual void compute_middle() {} virtual void compute_outer(int, int) {} virtual double single(int, int, int, int, double, double, double, double& fforce) { fforce = 0.0; return 0.0; } virtual void settings(int, char **) = 0; virtual void coeff(int, char **) = 0; virtual void init_style(); virtual void init_list(int, class NeighList *); virtual double init_one(int, int) {return 0.0;} virtual void init_tables(double, double *); virtual void init_tables_disp(double); virtual void free_tables(); virtual void free_disp_tables(); virtual void write_restart(FILE *) {} virtual void read_restart(FILE *) {} virtual void write_restart_settings(FILE *) {} virtual void read_restart_settings(FILE *) {} virtual void write_data(FILE *) {} + virtual void write_data_all(FILE *) {} virtual int pack_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} virtual double memory_usage(); // specific child-class methods for certain Pair styles virtual void *extract(const char *, int &) {return NULL;} virtual void swap_eam(double *, double **) {} virtual void reset_dt() {} virtual void min_xf_pointers(int, double **, double **) {} virtual void min_xf_get(int) {} virtual void min_x_set(int) {} virtual unsigned int data_mask() {return datamask;} virtual unsigned int data_mask_ext() {return datamask_ext;} protected: enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // mixing options int allocated; // 0/1 = whether arrays are allocated int suffix_flag; // suffix compatibility flag // pair_modify settings int offset_flag,mix_flag; // flags for offset and mixing double tabinner; // inner cutoff for Coulomb table double tabinner_disp; // inner cutoff for dispersion table // custom data type for accessing Coulomb tables typedef union {int i; float f;} union_int_float_t; double THIRD; int vflag_fdotr; int maxeatom,maxvatom; virtual void ev_setup(int, int); void ev_unset(); void ev_tally_full(int, double, double, double, double, double, double); void ev_tally_xyz_full(int, double, double, double, double, double, double, double, double); void ev_tally4(int, int, int, int, double, double *, double *, double *, double *, double *, double *); void ev_tally_tip4p(int, int *, double *, double, double); void v_tally2(int, int, double, double *); void v_tally_tensor(int, int, int, int, double, double, double, double, double, double); void virial_fdotr_compute(); inline int sbmask(int j) { return j >> SBBITS & 3; } }; } #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Too many total bits for bitmapped lookup table Table size specified via pair_modify command is too large. Note that a value of N generates a 2^N size table. E: Cannot have both pair_modify shift and tail set to yes These 2 options are contradictory. E: Cannot use pair tail corrections with 2d simulations The correction factors are only currently defined for 3d systems. W: Using pair tail corrections with nonperiodic system This is probably a bogus thing to do, since tail corrections are computed by integrating the density of a periodic system out to infinity. E: All pair coeffs are not set All pair coefficients must be set in the data file or by the pair_coeff command before running a simulation. E: Pair style does not support pair_write The pair style does not have a single() function, so it can not be invoked by pair write. E: Invalid atom types in pair_write command Atom types must range from 1 to Ntypes inclusive. E: Invalid style in pair_write command Self-explanatory. Check the input script. E: Invalid cutoffs in pair_write command Inner cutoff must be larger than 0.0 and less than outer cutoff. E: Cannot open pair_write file The specified output file for pair energies and forces cannot be opened. Check that the path and name are correct. E: Bitmapped lookup tables require int/float be same size Cannot use pair tables on this machine, because of word sizes. Use the pair_modify command with table 0 instead. W: Table inner cutoff >= outer cutoff You specified an inner cutoff for a Coulombic table that is longer than the global cutoff. Probably not what you wanted. E: Too many exponent bits for lookup table Table size specified via pair_modify command does not work with your machine's floating point representation. E: Too many mantissa bits for lookup table Table size specified via pair_modify command does not work with your machine's floating point representation. E: Too few bits for lookup table Table size specified via pair_modify command does not work with your machine's floating point representation. */ diff --git a/src/pair_coul_wolf.cpp b/src/pair_coul_wolf.cpp index af9751fa1..8cf821c2d 100644 --- a/src/pair_coul_wolf.cpp +++ b/src/pair_coul_wolf.cpp @@ -1,324 +1,327 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Yongfeng Zhang (INL), yongfeng.zhang@inl.gov ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_coul_wolf.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" #include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairCoulWolf::PairCoulWolf(LAMMPS *lmp) : Pair(lmp) {} +PairCoulWolf::PairCoulWolf(LAMMPS *lmp) : Pair(lmp) +{ + writedata = 0; +} /* ---------------------------------------------------------------------- */ PairCoulWolf::~PairCoulWolf() { if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); } } /* ---------------------------------------------------------------------- */ void PairCoulWolf::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; double rsq,forcecoul,factor_coul; double prefactor; double r,rexp; int *ilist,*jlist,*numneigh,**firstneigh; double erfcc,erfcd,v_sh,dvdrr,e_self,e_shift,f_shift,qisq; ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; // self and shifted coulombic energy e_self = v_sh = 0.0; e_shift = erfc(alf*cut_coul)/cut_coul; f_shift = -(e_shift+ 2.0*alf/MY_PIS * exp(-alf*alf*cut_coul*cut_coul)) / cut_coul; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; qtmp = q[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; qisq = qtmp*qtmp; e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e; if (evflag) ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_coul = special_coul[sbmask(j)]; j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cut_coulsq) { r = sqrt(rsq); prefactor = qqrd2e*qtmp*q[j]/r; erfcc = erfc(alf*r); erfcd = exp(-alf*alf*r*r); v_sh = (erfcc - e_shift*r) * prefactor; dvdrr = (erfcc/rsq + 2.0*alf/MY_PIS * erfcd/r) + f_shift; forcecoul = dvdrr*rsq*prefactor; if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; fpair = forcecoul / rsq; f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } if (eflag) { if (rsq < cut_coulsq) { ecoul = v_sh; if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } else ecoul = 0.0; } if (evflag) ev_tally(i,j,nlocal,newton_pair, 0.0,ecoul,fpair,delx,dely,delz); } } } if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairCoulWolf::allocate() { allocated = 1; int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; memory->create(cutsq,n+1,n+1,"pair:cutsq"); } /* ---------------------------------------------------------------------- global settings unlike other pair styles, there are no individual pair settings that these override ------------------------------------------------------------------------- */ void PairCoulWolf::settings(int narg, char **arg) { if (narg != 2) error->all(FLERR,"Illegal pair_style command"); alf = force->numeric(arg[0]); cut_coul = force->numeric(arg[1]); } /* ---------------------------------------------------------------------- set cutoffs for one or more type pairs, optional ------------------------------------------------------------------------- */ void PairCoulWolf::coeff(int narg, char **arg) { if (narg != 2) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { setflag[i][j] = 1; count++; } } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairCoulWolf::init_style() { if (!atom->q_flag) error->all(FLERR,"Pair coul/wolf requires atom attribute q"); int irequest = neighbor->request(this); cut_coulsq = cut_coul*cut_coul; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairCoulWolf::init_one(int i, int j) { return cut_coul; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairCoulWolf::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) fwrite(&setflag[i][j],sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairCoulWolf::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairCoulWolf::write_restart_settings(FILE *fp) { fwrite(&alf,sizeof(double),1,fp); fwrite(&cut_coul,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairCoulWolf::read_restart_settings(FILE *fp) { if (comm->me == 0) { fread(&alf,sizeof(double),1,fp); fread(&cut_coul,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&alf,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- only the pair part is calculated here ------------------------------------------------------------------------- */ double PairCoulWolf::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double r6inv,r,prefactor,rexp; double forcecoul,forceborn,phicoul; double e_shift,f_shift,dvdrr,erfcc,erfcd; e_shift = erfc(alf*cut_coul) / cut_coul; f_shift = -(e_shift+ 2.0*alf/MY_PIS * exp(-alf*alf*cut_coul*cut_coul)) / cut_coul; if (rsq < cut_coulsq) { r = sqrt(rsq); prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; erfcc = erfc(alf*r); erfcd = exp(-alf*alf*r*r); dvdrr = (erfcc/rsq + 2.0*alf/MY_PIS * erfcd/r) + f_shift; forcecoul = dvdrr*rsq*prefactor; if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; } else forcecoul = 0.0; fforce = forcecoul / rsq; double eng = 0.0; if (rsq < cut_coulsq) { phicoul = prefactor * (erfcc-e_shift*r); if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; eng += phicoul; } return eng; } diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index 15e7e7b85..becf29452 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -1,719 +1,730 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Paul Crozier (SNL) ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_lj_cut.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "update.h" #include "integrate.h" #include "respa.h" #include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; /* ---------------------------------------------------------------------- */ PairLJCut::PairLJCut(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; writedata = 1; } /* ---------------------------------------------------------------------- */ PairLJCut::~PairLJCut() { if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(cut); memory->destroy(epsilon); memory->destroy(sigma); memory->destroy(lj1); memory->destroy(lj2); memory->destroy(lj3); memory->destroy(lj4); memory->destroy(offset); } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj; int *ilist,*jlist,*numneigh,**firstneigh; evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } if (eflag) { evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,0.0,fpair,delx,dely,delz); } } } if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_inner() { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; int *ilist,*jlist,*numneigh,**firstneigh; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = listinner->inum; ilist = listinner->ilist; numneigh = listinner->numneigh; firstneigh = listinner->firstneigh; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cut_out_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; jtype = type[j]; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; if (rsq > cut_out_on_sq) { rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw); } f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } } } } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_middle() { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; int *ilist,*jlist,*numneigh,**firstneigh; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = listmiddle->inum; ilist = listmiddle->ilist; numneigh = listmiddle->numneigh; firstneigh = listmiddle->firstneigh; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; double cut_out_on = cut_respa[2]; double cut_out_off = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; jtype = type[j]; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; if (rsq < cut_in_on_sq) { rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; fpair *= rsw*rsw*(3.0 - 2.0*rsw); } if (rsq > cut_out_on_sq) { rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); } f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } } } } } /* ---------------------------------------------------------------------- */ void PairLJCut::compute_outer(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; int *ilist,*jlist,*numneigh,**firstneigh; evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = 0; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = listouter->inum; ilist = listouter->ilist; numneigh = listouter->numneigh; firstneigh = listouter->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { if (rsq > cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; if (rsq < cut_in_on_sq) { rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; fpair *= rsw*rsw*(3.0 - 2.0*rsw); } f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; if (newton_pair || j < nlocal) { f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; } } if (eflag) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } if (vflag) { if (rsq <= cut_in_off_sq) { r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fpair = factor_lj*forcelj*r2inv; } else if (rsq < cut_in_on_sq) fpair = factor_lj*forcelj*r2inv; } if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,0.0,fpair,delx,dely,delz); } } } } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairLJCut::allocate() { allocated = 1; int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(cut,n+1,n+1,"pair:cut"); memory->create(epsilon,n+1,n+1,"pair:epsilon"); memory->create(sigma,n+1,n+1,"pair:sigma"); memory->create(lj1,n+1,n+1,"pair:lj1"); memory->create(lj2,n+1,n+1,"pair:lj2"); memory->create(lj3,n+1,n+1,"pair:lj3"); memory->create(lj4,n+1,n+1,"pair:lj4"); memory->create(offset,n+1,n+1,"pair:offset"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairLJCut::settings(int narg, char **arg) { if (narg != 1) error->all(FLERR,"Illegal pair_style command"); cut_global = force->numeric(arg[0]); // reset cutoffs that have been explicitly set if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; } } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairLJCut::coeff(int narg, char **arg) { if (narg < 4 || narg > 5) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); double epsilon_one = force->numeric(arg[2]); double sigma_one = force->numeric(arg[3]); double cut_one = cut_global; if (narg == 5) cut_one = force->numeric(arg[4]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut[i][j] = cut_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairLJCut::init_style() { // request regular or rRESPA neighbor lists int irequest; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; if (respa == 0) irequest = neighbor->request(this); else if (respa == 1) { irequest = neighbor->request(this); neighbor->requests[irequest]->id = 1; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respainner = 1; irequest = neighbor->request(this); neighbor->requests[irequest]->id = 3; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respaouter = 1; } else { irequest = neighbor->request(this); neighbor->requests[irequest]->id = 1; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respainner = 1; irequest = neighbor->request(this); neighbor->requests[irequest]->id = 2; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respamiddle = 1; irequest = neighbor->request(this); neighbor->requests[irequest]->id = 3; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respaouter = 1; } } else irequest = neighbor->request(this); // set rRESPA cutoffs if (strstr(update->integrate_style,"respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; else cut_respa = NULL; } /* ---------------------------------------------------------------------- neighbor callback to inform pair style of neighbor list to use regular or rRESPA ------------------------------------------------------------------------- */ void PairLJCut::init_list(int id, NeighList *ptr) { if (id == 0) list = ptr; else if (id == 1) listinner = ptr; else if (id == 2) listmiddle = ptr; else if (id == 3) listouter = ptr; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairLJCut::init_one(int i, int j) { if (setflag[i][j] == 0) { epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], sigma[i][i],sigma[j][j]); sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); cut[i][j] = mix_distance(cut[i][i],cut[j][j]); } lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); if (offset_flag) { double ratio = sigma[i][j] / cut[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; lj3[j][i] = lj3[i][j]; lj4[j][i] = lj4[i][j]; offset[j][i] = offset[i][j]; // check interior rRESPA cutoff if (cut_respa && cut[i][j] < cut_respa[3]) error->all(FLERR,"Pair cutoff < Respa interior cutoff"); // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce if (tail_flag) { int *type = atom->type; int nlocal = atom->nlocal; double count[2],all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); double sig2 = sigma[i][j]*sigma[i][j]; double sig6 = sig2*sig2*sig2; double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; double rc6 = rc3*rc3; double rc9 = rc3*rc6; etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9); } return cut[i][j]; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJCut::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&epsilon[i][j],sizeof(double),1,fp); fwrite(&sigma[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJCut::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) { fread(&epsilon[i][j],sizeof(double),1,fp); fread(&sigma[i][j],sizeof(double),1,fp); fread(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJCut::write_restart_settings(FILE *fp) { fwrite(&cut_global,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); fwrite(&tail_flag,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJCut::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { fread(&cut_global,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); fread(&tail_flag,sizeof(int),1,fp); } MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); MPI_Bcast(&tail_flag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- proc 0 writes to data file ------------------------------------------------------------------------- */ void PairLJCut::write_data(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); } +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairLJCut::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][i],sigma[i][i],cut[i][j]); +} + /* ---------------------------------------------------------------------- */ double PairLJCut::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double r2inv,r6inv,forcelj,philj; r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); fforce = factor_lj*forcelj*r2inv; philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; return factor_lj*philj; } /* ---------------------------------------------------------------------- */ void *PairLJCut::extract(const char *str, int &dim) { dim = 2; if (strcmp(str,"epsilon") == 0) return (void *) epsilon; if (strcmp(str,"sigma") == 0) return (void *) sigma; return NULL; } diff --git a/src/pair_lj_cut.h b/src/pair_lj_cut.h index 46a7bf560..cccd9d1c6 100644 --- a/src/pair_lj_cut.h +++ b/src/pair_lj_cut.h @@ -1,81 +1,82 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef PAIR_CLASS PairStyle(lj/cut,PairLJCut) #else #ifndef LMP_PAIR_LJ_CUT_H #define LMP_PAIR_LJ_CUT_H #include "pair.h" namespace LAMMPS_NS { class PairLJCut : public Pair { public: PairLJCut(class LAMMPS *); virtual ~PairLJCut(); virtual void compute(int, int); void settings(int, char **); void coeff(int, char **); void init_style(); void init_list(int, class NeighList *); double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); void write_data(FILE *); + void write_data_all(FILE *); double single(int, int, int, int, double, double, double, double &); void *extract(const char *, int &); void compute_inner(); void compute_middle(); void compute_outer(int, int); protected: double cut_global; double **cut; double **epsilon,**sigma; double **lj1,**lj2,**lj3,**lj4,**offset; double *cut_respa; void allocate(); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Incorrect args for pair coefficients Self-explanatory. Check the input script or data file. E: Pair cutoff < Respa interior cutoff One or more pairwise cutoffs are too short to use with the specified rRESPA cutoffs. */ diff --git a/src/read_data.cpp b/src/read_data.cpp index f03805cbc..56a7651e0 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -1,1675 +1,1717 @@ /* ---------------------------------------------------------------------- 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 "lmptype.h" #include "mpi.h" #include "math.h" #include "string.h" #include "stdlib.h" #include "ctype.h" #include "read_data.h" #include "atom.h" #include "atom_vec.h" #include "atom_vec_ellipsoid.h" #include "atom_vec_line.h" #include "atom_vec_tri.h" #include "comm.h" #include "update.h" #include "modify.h" #include "fix.h" #include "force.h" #include "pair.h" #include "domain.h" #include "bond.h" #include "angle.h" #include "dihedral.h" #include "improper.h" #include "special.h" #include "error.h" #include "memory.h" using namespace LAMMPS_NS; #define MAXLINE 256 #define LB_FACTOR 1.1 #define CHUNK 1024 #define DELTA 4 // must be 2 or larger #define MAXBODY 20 // max # of lines in one body, also in Atom class // customize for new sections -#define NSECTIONS 24 // change when add to header::section_keywords +#define NSECTIONS 25 // change when add to header::section_keywords /* ---------------------------------------------------------------------- */ ReadData::ReadData(LAMMPS *lmp) : Pointers(lmp) { MPI_Comm_rank(world,&me); line = new char[MAXLINE]; keyword = new char[MAXLINE]; buffer = new char[CHUNK*MAXLINE]; narg = maxarg = 0; arg = NULL; // customize for new sections // pointers to atom styles that store extra info nellipsoids = 0; avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); nlines = 0; avec_line = (AtomVecLine *) atom->style_match("line"); ntris = 0; avec_tri = (AtomVecTri *) atom->style_match("tri"); nbodies = 0; avec_body = (AtomVecBody *) atom->style_match("body"); } /* ---------------------------------------------------------------------- */ ReadData::~ReadData() { delete [] line; delete [] keyword; delete [] buffer; memory->sfree(arg); } /* ---------------------------------------------------------------------- */ void ReadData::command(int narg, char **arg) { if (narg < 1) error->all(FLERR,"Illegal read_data command"); if (domain->box_exist) error->all(FLERR,"Cannot read_data after simulation box is defined"); if (domain->dimension == 2 && domain->zperiodic == 0) error->all(FLERR,"Cannot run 2d simulation with nonperiodic Z dimension"); // fixes that process data file info nfix = 0; fix_index = NULL; fix_header = NULL; fix_section = NULL; int iarg = 1; while (iarg < narg) { if (strcmp(arg[iarg],"fix") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal read_data command"); memory->grow(fix_index,nfix+1,"read_data:fix_index"); fix_header = (char **) memory->srealloc(fix_header,(nfix+1)*sizeof(char *), "read_data:fix_header"); fix_section = (char **) memory->srealloc(fix_section,(nfix+1)*sizeof(char *), "read_data:fix_section"); fix_index[nfix] = modify->find_fix(arg[iarg+1]); if (fix_index[nfix] < 0) error->all(FLERR,"Fix ID for read_data does not exist"); int n = strlen(arg[iarg+2]) + 1; fix_header[nfix] = new char[n]; strcpy(fix_header[nfix],arg[iarg+2]); n = strlen(arg[iarg+3]) + 1; fix_section[nfix] = new char[n]; strcpy(fix_section[nfix],arg[iarg+3]); nfix++; iarg += 4; } else error->all(FLERR,"Illegal read_data command"); } // scan data file to determine max topology needed per atom // allocate initial topology arrays if (atom->molecular) { if (me == 0) { if (screen) fprintf(screen,"Scanning data file ...\n"); open(arg[0]); header(0); scan(atom->bond_per_atom,atom->angle_per_atom, atom->dihedral_per_atom,atom->improper_per_atom); if (compressed) pclose(fp); else fclose(fp); atom->bond_per_atom += atom->extra_bond_per_atom; } MPI_Bcast(&atom->bond_per_atom,1,MPI_INT,0,world); MPI_Bcast(&atom->angle_per_atom,1,MPI_INT,0,world); MPI_Bcast(&atom->dihedral_per_atom,1,MPI_INT,0,world); MPI_Bcast(&atom->improper_per_atom,1,MPI_INT,0,world); } else atom->bond_per_atom = atom->angle_per_atom = atom->dihedral_per_atom = atom->improper_per_atom = 0; // read header info if (me == 0) { if (screen) fprintf(screen,"Reading data file ...\n"); open(arg[0]); } header(1); domain->box_exist = 1; // problem setup using info from header update->ntimestep = 0; int n; if (comm->nprocs == 1) n = static_cast (atom->natoms); else n = static_cast (LB_FACTOR * atom->natoms / comm->nprocs); atom->allocate_type_arrays(); atom->avec->grow(n); n = atom->nmax; domain->print_box(" "); domain->set_initial_box(); domain->set_global_box(); comm->set_proc_grid(); domain->set_local_box(); // customize for new sections // read rest of file in free format int atomflag = 0; while (strlen(keyword)) { // allow special fixes first chance to match and process the section // if fix matches, continue to next section if (nfix) { for (n = 0; n < nfix; n++) if (strstr(line,fix_section[n])) { bigint nlines = modify->fix[fix_index[n]]->read_data_skip_lines(keyword); fix(n,keyword,nlines); parse_keyword(0,1); break; } if (n < nfix) continue; } if (strcmp(keyword,"Atoms") == 0) { atoms(); atomflag = 1; } else if (strcmp(keyword,"Velocities") == 0) { if (atomflag == 0) error->all(FLERR,"Must read Atoms before Velocities"); velocities(); } else if (strcmp(keyword,"Ellipsoids") == 0) { if (!avec_ellipsoid) error->all(FLERR,"Invalid data file section: Ellipsoids"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Ellipsoids"); bonus(nellipsoids,(AtomVec *) avec_ellipsoid,"ellipsoids"); } else if (strcmp(keyword,"Lines") == 0) { if (!avec_line) error->all(FLERR,"Invalid data file section: Lines"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Lines"); bonus(nlines,(AtomVec *) avec_line,"lines"); } else if (strcmp(keyword,"Triangles") == 0) { if (!avec_tri) error->all(FLERR,"Invalid data file section: Triangles"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Triangles"); bonus(ntris,(AtomVec *) avec_tri,"triangles"); } else if (strcmp(keyword,"Bodies") == 0) { if (!avec_body) error->all(FLERR,"Invalid data file section: Bodies"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Bodies"); bodies(); } else if (strcmp(keyword,"Bonds") == 0) { if (atom->avec->bonds_allow == 0) error->all(FLERR,"Invalid data file section: Bonds"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Bonds"); bonds(); } else if (strcmp(keyword,"Angles") == 0) { if (atom->avec->angles_allow == 0) error->all(FLERR,"Invalid data file section: Angles"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Angles"); angles(); } else if (strcmp(keyword,"Dihedrals") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: Dihedrals"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Dihedrals"); dihedrals(); } else if (strcmp(keyword,"Impropers") == 0) { if (atom->avec->impropers_allow == 0) error->all(FLERR,"Invalid data file section: Impropers"); if (atomflag == 0) error->all(FLERR,"Must read Atoms before Impropers"); impropers(); } else if (strcmp(keyword,"Masses") == 0) { mass(); } else if (strcmp(keyword,"Pair Coeffs") == 0) { if (force->pair == NULL) error->all(FLERR,"Must define pair_style before Pair Coeffs"); paircoeffs(); + } else if (strcmp(keyword,"PairIJ Coeffs") == 0) { + if (force->pair == NULL) + error->all(FLERR,"Must define pair_style before PairIJ Coeffs"); + pairIJcoeffs(); } else if (strcmp(keyword,"Bond Coeffs") == 0) { if (atom->avec->bonds_allow == 0) error->all(FLERR,"Invalid data file section: Bond Coeffs"); if (force->bond == NULL) error->all(FLERR,"Must define bond_style before Bond Coeffs"); bondcoeffs(); } else if (strcmp(keyword,"Angle Coeffs") == 0) { if (atom->avec->angles_allow == 0) error->all(FLERR,"Invalid data file section: Angle Coeffs"); if (force->angle == NULL) error->all(FLERR,"Must define angle_style before Angle Coeffs"); anglecoeffs(0); } else if (strcmp(keyword,"Dihedral Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: Dihedral Coeffs"); if (force->dihedral == NULL) error->all(FLERR,"Must define dihedral_style before Dihedral Coeffs"); dihedralcoeffs(0); } else if (strcmp(keyword,"Improper Coeffs") == 0) { if (atom->avec->impropers_allow == 0) error->all(FLERR,"Invalid data file section: Improper Coeffs"); if (force->improper == NULL) error->all(FLERR,"Must define improper_style before Improper Coeffs"); impropercoeffs(0); } else if (strcmp(keyword,"BondBond Coeffs") == 0) { if (atom->avec->angles_allow == 0) error->all(FLERR,"Invalid data file section: BondBond Coeffs"); if (force->angle == NULL) error->all(FLERR,"Must define angle_style before BondBond Coeffs"); anglecoeffs(1); } else if (strcmp(keyword,"BondAngle Coeffs") == 0) { if (atom->avec->angles_allow == 0) error->all(FLERR,"Invalid data file section: BondAngle Coeffs"); if (force->angle == NULL) error->all(FLERR,"Must define angle_style before BondAngle Coeffs"); anglecoeffs(2); } else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs"); if (force->dihedral == NULL) error->all(FLERR, "Must define dihedral_style before " "MiddleBondTorsion Coeffs"); dihedralcoeffs(1); } else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs"); if (force->dihedral == NULL) error->all(FLERR, "Must define dihedral_style before EndBondTorsion Coeffs"); dihedralcoeffs(2); } else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs"); if (force->dihedral == NULL) error->all(FLERR, "Must define dihedral_style before AngleTorsion Coeffs"); dihedralcoeffs(3); } else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs"); if (force->dihedral == NULL) error->all(FLERR, "Must define dihedral_style before " "AngleAngleTorsion Coeffs"); dihedralcoeffs(4); } else if (strcmp(keyword,"BondBond13 Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->all(FLERR,"Invalid data file section: BondBond13 Coeffs"); if (force->dihedral == NULL) error->all(FLERR,"Must define dihedral_style before BondBond13 Coeffs"); dihedralcoeffs(5); } else if (strcmp(keyword,"AngleAngle Coeffs") == 0) { if (atom->avec->impropers_allow == 0) error->all(FLERR,"Invalid data file section: AngleAngle Coeffs"); if (force->improper == NULL) error->all(FLERR,"Must define improper_style before AngleAngle Coeffs"); impropercoeffs(1); } else { char str[128]; sprintf(str,"Unknown identifier in data file: %s",keyword); error->all(FLERR,str); } parse_keyword(0,1); } // close file if (me == 0) { if (compressed) pclose(fp); else fclose(fp); } // error if natoms > 0 yet no atoms were read if (atom->natoms > 0 && atomflag == 0) error->all(FLERR,"No atoms in data file"); // create bond topology now that system is defined if (atom->molecular) { Special special(lmp); special.build(); } } /* ---------------------------------------------------------------------- read free-format header of data file if flag = 0, only called by proc 0 if flag = 1, called by all procs so bcast lines as read them 1st line and blank lines are skipped non-blank lines are checked for header keywords and leading value is read header ends with EOF or non-blank line containing no header keyword if EOF, line is set to blank line else line has first keyword line for rest of file ------------------------------------------------------------------------- */ void ReadData::header(int flag) { int n; char *ptr; // customize for new sections const char *section_keywords[NSECTIONS] = {"Atoms","Velocities","Ellipsoids","Lines","Triangles","Bodies", "Bonds","Angles","Dihedrals","Impropers", - "Masses","Pair Coeffs","Bond Coeffs","Angle Coeffs", + "Masses","Pair Coeffs","PairIJ Coeffs","Bond Coeffs","Angle Coeffs", "Dihedral Coeffs","Improper Coeffs", "BondBond Coeffs","BondAngle Coeffs","MiddleBondTorsion Coeffs", "EndBondTorsion Coeffs","AngleTorsion Coeffs", "AngleAngleTorsion Coeffs","BondBond13 Coeffs","AngleAngle Coeffs"}; // skip 1st line of file if (me == 0) { char *eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); } // customize for new header lines while (1) { // read a line and bcast length if flag is set if (me == 0) { if (fgets(line,MAXLINE,fp) == NULL) n = 0; else n = strlen(line) + 1; } if (flag) MPI_Bcast(&n,1,MPI_INT,0,world); // if n = 0 then end-of-file so return with blank line if (n == 0) { line[0] = '\0'; return; } // bcast line if flag is set if (flag) MPI_Bcast(line,n,MPI_CHAR,0,world); // trim anything from '#' onward // if line is blank, continue if (ptr = strchr(line,'#')) *ptr = '\0'; if (strspn(line," \t\n\r") == strlen(line)) continue; // allow special fixes first chance to match and process the line // if fix matches, continue to next header line if (nfix) { for (n = 0; n < nfix; n++) if (strstr(line,fix_header[n])) { modify->fix[fix_index[n]]->read_data_header(line); break; } if (n < nfix) continue; } // search line for header keyword and set corresponding variable if (strstr(line,"atoms")) sscanf(line,BIGINT_FORMAT,&atom->natoms); // check for these first // otherwise "triangles" will be matched as "angles" else if (strstr(line,"ellipsoids")) { if (!avec_ellipsoid && me == 0) error->one(FLERR,"No ellipsoids allowed with this atom style"); sscanf(line,BIGINT_FORMAT,&nellipsoids); } else if (strstr(line,"lines")) { if (!avec_line && me == 0) error->one(FLERR,"No lines allowed with this atom style"); sscanf(line,BIGINT_FORMAT,&nlines); } else if (strstr(line,"triangles")) { if (!avec_tri && me == 0) error->one(FLERR,"No triangles allowed with this atom style"); sscanf(line,BIGINT_FORMAT,&ntris); } else if (strstr(line,"bodies")) { if (!avec_body && me == 0) error->one(FLERR,"No bodies allowed with this atom style"); sscanf(line,BIGINT_FORMAT,&nbodies); } else if (strstr(line,"bonds")) sscanf(line,BIGINT_FORMAT,&atom->nbonds); else if (strstr(line,"angles")) sscanf(line,BIGINT_FORMAT,&atom->nangles); else if (strstr(line,"dihedrals")) sscanf(line,BIGINT_FORMAT, &atom->ndihedrals); else if (strstr(line,"impropers")) sscanf(line,BIGINT_FORMAT, &atom->nimpropers); else if (strstr(line,"atom types")) sscanf(line,"%d",&atom->ntypes); else if (strstr(line,"bond types")) sscanf(line,"%d",&atom->nbondtypes); else if (strstr(line,"angle types")) sscanf(line,"%d",&atom->nangletypes); else if (strstr(line,"dihedral types")) sscanf(line,"%d",&atom->ndihedraltypes); else if (strstr(line,"improper types")) sscanf(line,"%d",&atom->nimpropertypes); else if (strstr(line,"extra bond per atom")) sscanf(line,"%d",&atom->extra_bond_per_atom); else if (strstr(line,"xlo xhi")) sscanf(line,"%lg %lg",&domain->boxlo[0],&domain->boxhi[0]); else if (strstr(line,"ylo yhi")) sscanf(line,"%lg %lg",&domain->boxlo[1],&domain->boxhi[1]); else if (strstr(line,"zlo zhi")) sscanf(line,"%lg %lg",&domain->boxlo[2],&domain->boxhi[2]); else if (strstr(line,"xy xz yz")) { domain->triclinic = 1; sscanf(line,"%lg %lg %lg",&domain->xy,&domain->xz,&domain->yz); } else break; } // error check on total system size if (atom->natoms < 0 || atom->natoms > MAXBIGINT || atom->nbonds < 0 || atom->nbonds > MAXBIGINT || atom->nangles < 0 || atom->nangles > MAXBIGINT || atom->ndihedrals < 0 || atom->ndihedrals > MAXBIGINT || atom->nimpropers < 0 || atom->nimpropers > MAXBIGINT) { if (me == 0) error->one(FLERR,"System in data file is too big"); } // check that exiting string is a valid section keyword parse_keyword(1,flag); for (n = 0; n < NSECTIONS; n++) if (strcmp(keyword,section_keywords[n]) == 0) break; if (n == NSECTIONS && me == 0) { char str[128]; sprintf(str,"Unknown identifier in data file: %s",keyword); error->one(FLERR,str); } // error check on consistency of header values if ((atom->nbonds || atom->nbondtypes) && atom->avec->bonds_allow == 0 && me == 0) error->one(FLERR,"No bonds allowed with this atom style"); if ((atom->nangles || atom->nangletypes) && atom->avec->angles_allow == 0 && me == 0) error->one(FLERR,"No angles allowed with this atom style"); if ((atom->ndihedrals || atom->ndihedraltypes) && atom->avec->dihedrals_allow == 0 && me == 0) error->one(FLERR,"No dihedrals allowed with this atom style"); if ((atom->nimpropers || atom->nimpropertypes) && atom->avec->impropers_allow == 0 && me == 0) error->one(FLERR,"No impropers allowed with this atom style"); if (atom->nbonds > 0 && atom->nbondtypes <= 0 && me == 0) error->one(FLERR,"Bonds defined but no bond types"); if (atom->nangles > 0 && atom->nangletypes <= 0 && me == 0) error->one(FLERR,"Angles defined but no angle types"); if (atom->ndihedrals > 0 && atom->ndihedraltypes <= 0 && me == 0) error->one(FLERR,"Dihedrals defined but no dihedral types"); if (atom->nimpropers > 0 && atom->nimpropertypes <= 0 && me == 0) error->one(FLERR,"Impropers defined but no improper types"); } /* ---------------------------------------------------------------------- read all atoms ------------------------------------------------------------------------- */ void ReadData::atoms() { int i,m,nchunk; bigint nread = 0; bigint natoms = atom->natoms; while (nread < natoms) { if (natoms-nread > CHUNK) nchunk = CHUNK; else nchunk = natoms-nread; if (me == 0) { char *eof; m = 0; for (i = 0; i < nchunk; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); atom->data_atoms(nchunk,buffer); nread += nchunk; } // check that all atoms were assigned correctly bigint tmp = atom->nlocal; MPI_Allreduce(&tmp,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " atoms\n",natoms); if (logfile) fprintf(logfile," " BIGINT_FORMAT " atoms\n",natoms); } if (natoms != atom->natoms) error->all(FLERR,"Did not assign all atoms correctly"); // if any atom ID < 0, error // if all atom IDs = 0, tag_enable = 0 // if any atom ID > 0, error if any atom ID == 0 // not checking if atom IDs > natoms or are unique int nlocal = atom->nlocal; int *tag = atom->tag; int flag = 0; for (int i = 0; i < nlocal; i++) if (tag[i] < 0) flag = 1; int flag_all; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all(FLERR,"Invalid atom ID in Atoms section of data file"); flag = 0; for (int i = 0; i < nlocal; i++) if (tag[i] > 0) flag = 1; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world); if (flag_all == 0) atom->tag_enable = 0; if (atom->tag_enable) { flag = 0; for (int i = 0; i < nlocal; i++) if (tag[i] == 0) flag = 1; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); if (flag_all) error->all(FLERR,"Invalid atom ID in Atoms section of data file"); } // create global mapping if (atom->map_style) { atom->map_init(); atom->map_set(); } } /* ---------------------------------------------------------------------- read all velocities to find atoms, must build atom map if not a molecular system ------------------------------------------------------------------------- */ void ReadData::velocities() { int i,m,nchunk; int mapflag = 0; if (atom->map_style == 0) { mapflag = 1; atom->map_style = 1; atom->map_init(); atom->map_set(); } bigint nread = 0; bigint natoms = atom->natoms; while (nread < natoms) { if (natoms-nread > CHUNK) nchunk = CHUNK; else nchunk = natoms-nread; if (me == 0) { char *eof; m = 0; for (i = 0; i < nchunk; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); atom->data_vels(nchunk,buffer); nread += nchunk; } if (mapflag) { atom->map_delete(); atom->map_style = 0; } if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " velocities\n",natoms); if (logfile) fprintf(logfile," " BIGINT_FORMAT " velocities\n",natoms); } } /* ---------------------------------------------------------------------- read all bonus data to find atoms, must build atom map if not a molecular system ------------------------------------------------------------------------- */ void ReadData::bonus(bigint nbonus, AtomVec *ptr, const char *type) { int i,m,nchunk; int mapflag = 0; if (atom->map_style == 0) { mapflag = 1; atom->map_style = 1; atom->map_init(); atom->map_set(); } bigint nread = 0; bigint natoms = nbonus; while (nread < natoms) { if (natoms-nread > CHUNK) nchunk = CHUNK; else nchunk = natoms-nread; if (me == 0) { char *eof; m = 0; for (i = 0; i < nchunk; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); atom->data_bonus(nchunk,buffer,ptr); nread += nchunk; } if (mapflag) { atom->map_delete(); atom->map_style = 0; } if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " %s\n",natoms,type); if (logfile) fprintf(logfile," " BIGINT_FORMAT " %s\n",natoms,type); } } /* ---------------------------------------------------------------------- read all body data variable amount of info per body, described by ninteger and ndouble to find atoms, must build atom map if not a molecular system ------------------------------------------------------------------------- */ void ReadData::bodies() { int i,m,nchunk,nmax,ninteger,ndouble,tmp,onebody; char *eof; int mapflag = 0; if (atom->map_style == 0) { mapflag = 1; atom->map_style = 1; atom->map_init(); atom->map_set(); } // nmax = max # of bodies to read in this chunk // nchunk = actual # readr bigint nread = 0; bigint natoms = nbodies; while (nread < natoms) { if (natoms-nread > CHUNK) nmax = CHUNK; else nmax = natoms-nread; if (me == 0) { nchunk = 0; nlines = 0; m = 0; while (nchunk < nmax && nlines <= CHUNK-MAXBODY) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(&buffer[m],"%d %d %d",&tmp,&ninteger,&ndouble); m += strlen(&buffer[m]); onebody = 0; if (ninteger) onebody += (ninteger-1)/10 + 1; if (ndouble) onebody += (ndouble-1)/10 + 1; if (onebody+1 > MAXBODY) error->one(FLERR, "Too many lines in one body in data file - boost MAXBODY"); for (i = 0; i < onebody; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } nchunk++; nlines += onebody+1; } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&nchunk,1,MPI_INT,0,world); MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); atom->data_bodies(nchunk,buffer,avec_body); nread += nchunk; } if (mapflag) { atom->map_delete(); atom->map_style = 0; } if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " bodies\n",natoms); if (logfile) fprintf(logfile," " BIGINT_FORMAT " bodies\n",natoms); } } /* ---------------------------------------------------------------------- */ void ReadData::bonds() { int i,m,nchunk; bigint nread = 0; bigint nbonds = atom->nbonds; while (nread < nbonds) { nchunk = MIN(nbonds-nread,CHUNK); if (me == 0) { char *eof; m = 0; for (i = 0; i < nchunk; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); atom->data_bonds(nchunk,buffer); nread += nchunk; } // check that bonds were assigned correctly int nlocal = atom->nlocal; bigint sum; bigint n = 0; for (i = 0; i < nlocal; i++) n += atom->num_bond[i]; MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world); int factor = 1; if (!force->newton_bond) factor = 2; if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " bonds\n",sum/factor); if (logfile) fprintf(logfile," " BIGINT_FORMAT " bonds\n",sum/factor); } if (sum != factor*atom->nbonds) error->all(FLERR,"Bonds assigned incorrectly"); } /* ---------------------------------------------------------------------- */ void ReadData::angles() { int i,m,nchunk; bigint nread = 0; bigint nangles = atom->nangles; while (nread < nangles) { nchunk = MIN(nangles-nread,CHUNK); if (me == 0) { char *eof; m = 0; for (i = 0; i < nchunk; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); atom->data_angles(nchunk,buffer); nread += nchunk; } // check that ang int nlocal = atom->nlocal; bigint sum; bigint n = 0; for (i = 0; i < nlocal; i++) n += atom->num_angle[i]; MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world); int factor = 1; if (!force->newton_bond) factor = 3; if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " angles\n",sum/factor); if (logfile) fprintf(logfile," " BIGINT_FORMAT " angles\n",sum/factor); } if (sum != factor*atom->nangles) error->all(FLERR,"Angles assigned incorrectly"); } /* ---------------------------------------------------------------------- */ void ReadData::dihedrals() { int i,m,nchunk; bigint nread = 0; bigint ndihedrals = atom->ndihedrals; while (nread < ndihedrals) { nchunk = MIN(ndihedrals-nread,CHUNK); if (me == 0) { char *eof; m = 0; for (i = 0; i < nchunk; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); atom->data_dihedrals(nchunk,buffer); nread += nchunk; } // check that dihedrals were assigned correctly int nlocal = atom->nlocal; bigint sum; bigint n = 0; for (i = 0; i < nlocal; i++) n += atom->num_dihedral[i]; MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world); int factor = 1; if (!force->newton_bond) factor = 4; if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " dihedrals\n",sum/factor); if (logfile) fprintf(logfile," " BIGINT_FORMAT " dihedrals\n",sum/factor); } if (sum != factor*atom->ndihedrals) error->all(FLERR,"Dihedrals assigned incorrectly"); } /* ---------------------------------------------------------------------- */ void ReadData::impropers() { int i,m,nchunk; bigint nread = 0; bigint nimpropers = atom->nimpropers; while (nread < nimpropers) { nchunk = MIN(nimpropers-nread,CHUNK); if (me == 0) { char *eof; m = 0; for (i = 0; i < nchunk; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); atom->data_impropers(nchunk,buffer); nread += nchunk; } // check that impropers were assigned correctly int nlocal = atom->nlocal; bigint sum; bigint n = 0; for (i = 0; i < nlocal; i++) n += atom->num_improper[i]; MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world); int factor = 1; if (!force->newton_bond) factor = 4; if (me == 0) { if (screen) fprintf(screen," " BIGINT_FORMAT " impropers\n",sum/factor); if (logfile) fprintf(logfile," " BIGINT_FORMAT " impropers\n",sum/factor); } if (sum != factor*atom->nimpropers) error->all(FLERR,"Impropers assigned incorrectly"); } /* ---------------------------------------------------------------------- */ void ReadData::mass() { int i,m; char *buf = new char[atom->ntypes*MAXLINE]; char *original = buf; if (me == 0) { char *eof; m = 0; for (i = 0; i < atom->ntypes; i++) { eof = fgets(&buf[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buf[m]); if (buf[m-1] != '\n') strcpy(&buf[m++],"\n"); buf[m-1] = '\0'; } } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buf,m,MPI_CHAR,0,world); for (i = 0; i < atom->ntypes; i++) { atom->set_mass(buf); buf += strlen(buf) + 1; } delete [] original; } /* ---------------------------------------------------------------------- */ void ReadData::paircoeffs() { int i,m; char *buf = new char[atom->ntypes*MAXLINE]; char *original = buf; if (me == 0) { char *eof; m = 0; for (i = 0; i < atom->ntypes; i++) { eof = fgets(&buf[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buf[m]); if (buf[m-1] != '\n') strcpy(&buf[m++],"\n"); buf[m-1] = '\0'; } } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buf,m,MPI_CHAR,0,world); for (i = 0; i < atom->ntypes; i++) { m = strlen(buf) + 1; parse_coeffs(buf,NULL,1); force->pair->coeff(narg,arg); buf += m; } delete [] original; } /* ---------------------------------------------------------------------- */ +void ReadData::pairIJcoeffs() +{ + int i,j,m; + char *buf = new char[atom->ntypes*(atom->ntypes+1)/2 * MAXLINE]; + char *original = buf; + + if (me == 0) { + char *eof; + m = 0; + for (i = 0; i < atom->ntypes; i++) + for (j = i; j < atom->ntypes; j++) { + eof = fgets(&buf[m],MAXLINE,fp); + if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); + m += strlen(&buf[m]); + if (buf[m-1] != '\n') strcpy(&buf[m++],"\n"); + buf[m-1] = '\0'; + } + } + + MPI_Bcast(&m,1,MPI_INT,0,world); + MPI_Bcast(buf,m,MPI_CHAR,0,world); + + for (i = 0; i < atom->ntypes; i++) + for (j = i; j < atom->ntypes; j++) { + m = strlen(buf) + 1; + parse_coeffs(buf,NULL,0); + force->pair->coeff(narg,arg); + buf += m; + } + delete [] original; +} + +/* ---------------------------------------------------------------------- */ + void ReadData::bondcoeffs() { int i,m; char *buf = new char[atom->nbondtypes*MAXLINE]; char *original = buf; if (me == 0) { char *eof; m = 0; for (i = 0; i < atom->nbondtypes; i++) { eof = fgets(&buf[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buf[m]); if (buf[m-1] != '\n') strcpy(&buf[m++],"\n"); buf[m-1] = '\0'; } } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buf,m,MPI_CHAR,0,world); for (i = 0; i < atom->nbondtypes; i++) { m = strlen(buf) + 1; parse_coeffs(buf,NULL,0); force->bond->coeff(narg,arg); buf += m; } delete [] original; } /* ---------------------------------------------------------------------- */ void ReadData::anglecoeffs(int which) { int i,m; char *buf = new char[atom->nangletypes*MAXLINE]; char *original = buf; if (me == 0) { char *eof; m = 0; for (i = 0; i < atom->nangletypes; i++) { eof = fgets(&buf[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buf[m]); if (buf[m-1] != '\n') strcpy(&buf[m++],"\n"); buf[m-1] = '\0'; } } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buf,m,MPI_CHAR,0,world); for (i = 0; i < atom->nangletypes; i++) { m = strlen(buf) + 1; if (which == 0) parse_coeffs(buf,NULL,0); else if (which == 1) parse_coeffs(buf,"bb",0); else if (which == 2) parse_coeffs(buf,"ba",0); force->angle->coeff(narg,arg); buf += m; } delete [] original; } /* ---------------------------------------------------------------------- */ void ReadData::dihedralcoeffs(int which) { int i,m; char *buf = new char[atom->ndihedraltypes*MAXLINE]; char *original = buf; if (me == 0) { char *eof; m = 0; for (i = 0; i < atom->ndihedraltypes; i++) { eof = fgets(&buf[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buf[m]); if (buf[m-1] != '\n') strcpy(&buf[m++],"\n"); buf[m-1] = '\0'; } } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buf,m,MPI_CHAR,0,world); for (i = 0; i < atom->ndihedraltypes; i++) { m = strlen(buf) + 1; if (which == 0) parse_coeffs(buf,NULL,0); else if (which == 1) parse_coeffs(buf,"mbt",0); else if (which == 2) parse_coeffs(buf,"ebt",0); else if (which == 3) parse_coeffs(buf,"at",0); else if (which == 4) parse_coeffs(buf,"aat",0); else if (which == 5) parse_coeffs(buf,"bb13",0); force->dihedral->coeff(narg,arg); buf += m; } delete [] original; } /* ---------------------------------------------------------------------- */ void ReadData::impropercoeffs(int which) { int i,m; char *buf = new char[atom->nimpropertypes*MAXLINE]; char *original = buf; if (me == 0) { char *eof; m = 0; for (i = 0; i < atom->nimpropertypes; i++) { eof = fgets(&buf[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buf[m]); if (buf[m-1] != '\n') strcpy(&buf[m++],"\n"); buf[m-1] = '\0'; } } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buf,m,MPI_CHAR,0,world); for (i = 0; i < atom->nimpropertypes; i++) { m = strlen(buf) + 1; if (which == 0) parse_coeffs(buf,NULL,0); else if (which == 1) parse_coeffs(buf,"aa",0); force->improper->coeff(narg,arg); buf += m; } delete [] original; } /* ---------------------------------------------------------------------- read fix section, pass lines to fix to process n = index of fix ------------------------------------------------------------------------- */ void ReadData::fix(int ifix, char *line, bigint nlines) { int i,m,nchunk; bigint nread = 0; while (nread < nlines) { if (nlines-nread > CHUNK) nchunk = CHUNK; else nchunk = nlines-nread; if (me == 0) { char *eof; m = 0; for (i = 0; i < nchunk; i++) { eof = fgets(&buffer[m],MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); m += strlen(&buffer[m]); } if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n"); m++; } MPI_Bcast(&m,1,MPI_INT,0,world); MPI_Bcast(buffer,m,MPI_CHAR,0,world); modify->fix[ifix]->read_data_section(line,nchunk,buffer); nread += nchunk; } } /* ---------------------------------------------------------------------- proc 0 scans the data file for topology maximums ------------------------------------------------------------------------- */ void ReadData::scan(int &bond_per_atom, int &angle_per_atom, int &dihedral_per_atom, int &improper_per_atom) { int i,tmp1,tmp2,atom1,atom2,atom3,atom4; char *eof; if (atom->natoms > MAXSMALLINT) error->one(FLERR,"Molecular data file has too many atoms"); // customize for new sections int natoms = static_cast (atom->natoms); bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0; int ellipsoid_flag = 0; int line_flag = 0; int tri_flag = 0; int body_flag = 0; // customize for new sections // allocate topology counting vector // initially, array length = 1 to natoms // will grow via reallocate() if atom IDs > natoms int cmax = natoms + 1; int *count; memory->create(count,cmax,"read_data:count"); while (strlen(keyword)) { // allow special fixes first chance to match and process the section // if fix matches, continue to next section if (nfix) { for (i = 0; i < nfix; i++) { printf("LINE SECTION %s %s\n",line,fix_section[i]); if (strstr(line,fix_section[i])) { int n = modify->fix[fix_index[i]]->read_data_skip_lines(keyword); printf("NLINES SKIP %d\n",n); skip_lines(n); parse_keyword(0,0); break; } } if (i < nfix) continue; } if (strcmp(keyword,"Masses") == 0) skip_lines(atom->ntypes); else if (strcmp(keyword,"Atoms") == 0) skip_lines(natoms); else if (strcmp(keyword,"Velocities") == 0) skip_lines(natoms); else if (strcmp(keyword,"Ellipsoids") == 0) { if (!avec_ellipsoid) error->one(FLERR,"Invalid data file section: Ellipsoids"); ellipsoid_flag = 1; skip_lines(nellipsoids); } else if (strcmp(keyword,"Lines") == 0) { if (!avec_line) error->one(FLERR,"Invalid data file section: Lines"); line_flag = 1; skip_lines(nlines); } else if (strcmp(keyword,"Triangles") == 0) { if (!avec_tri) error->one(FLERR,"Invalid data file section: Triangles"); tri_flag = 1; skip_lines(ntris); } else if (strcmp(keyword,"Bodies") == 0) { if (!avec_body) error->one(FLERR,"Invalid data file section: Bodies"); body_flag = 1; skip_lines(nbodies); } else if (strcmp(keyword,"Pair Coeffs") == 0) { if (force->pair == NULL) error->one(FLERR,"Must define pair_style before Pair Coeffs"); skip_lines(atom->ntypes); + } else if (strcmp(keyword,"PairIJ Coeffs") == 0) { + if (force->pair == NULL) + error->one(FLERR,"Must define pair_style before Pair Coeffs"); + skip_lines(atom->ntypes*(atom->ntypes+1)/2); } else if (strcmp(keyword,"Bond Coeffs") == 0) { if (atom->avec->bonds_allow == 0) error->one(FLERR,"Invalid data file section: Bond Coeffs"); if (force->bond == NULL) error->one(FLERR,"Must define bond_style before Bond Coeffs"); skip_lines(atom->nbondtypes); } else if (strcmp(keyword,"Angle Coeffs") == 0) { if (atom->avec->angles_allow == 0) error->one(FLERR,"Invalid data file section: Angle Coeffs"); if (force->angle == NULL) error->one(FLERR,"Must define angle_style before Angle Coeffs"); skip_lines(atom->nangletypes); } else if (strcmp(keyword,"Dihedral Coeffs") == 0) { skip_lines(atom->ndihedraltypes); if (atom->avec->dihedrals_allow == 0) error->one(FLERR,"Invalid data file section: Dihedral Coeffs"); if (force->dihedral == NULL) error->one(FLERR,"Must define dihedral_style before Dihedral Coeffs"); } else if (strcmp(keyword,"Improper Coeffs") == 0) { if (atom->avec->impropers_allow == 0) error->one(FLERR,"Invalid data file section: Improper Coeffs"); if (force->improper == NULL) error->one(FLERR,"Must define improper_style before Improper Coeffs"); skip_lines(atom->nimpropertypes); } else if (strcmp(keyword,"BondBond Coeffs") == 0) { if (atom->avec->angles_allow == 0) error->one(FLERR,"Invalid data file section: BondBond Coeffs"); if (force->angle == NULL) error->one(FLERR,"Must define angle_style before BondBond Coeffs"); skip_lines(atom->nangletypes); } else if (strcmp(keyword,"BondAngle Coeffs") == 0) { if (atom->avec->angles_allow == 0) error->one(FLERR,"Invalid data file section: BondAngle Coeffs"); if (force->angle == NULL) error->one(FLERR,"Must define angle_style before BondAngle Coeffs"); skip_lines(atom->nangletypes); } else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->one(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs"); if (force->dihedral == NULL) error->one(FLERR, "Must define dihedral_style before " "MiddleBondTorsion Coeffs"); skip_lines(atom->ndihedraltypes); } else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->one(FLERR,"Invalid data file section: EndBondTorsion Coeffs"); if (force->dihedral == NULL) error->one(FLERR, "Must define dihedral_style before EndBondTorsion Coeffs"); skip_lines(atom->ndihedraltypes); } else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->one(FLERR,"Invalid data file section: AngleTorsion Coeffs"); if (force->dihedral == NULL) error->one(FLERR, "Must define dihedral_style before AngleTorsion Coeffs"); skip_lines(atom->ndihedraltypes); } else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->one(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs"); if (force->dihedral == NULL) error->one(FLERR, "Must define dihedral_style before " "AngleAngleTorsion Coeffs"); skip_lines(atom->ndihedraltypes); } else if (strcmp(keyword,"BondBond13 Coeffs") == 0) { if (atom->avec->dihedrals_allow == 0) error->one(FLERR,"Invalid data file section: BondBond13 Coeffs"); if (force->dihedral == NULL) error->one(FLERR,"Must define dihedral_style before BondBond13 Coeffs"); skip_lines(atom->ndihedraltypes); } else if (strcmp(keyword,"AngleAngle Coeffs") == 0) { if (atom->avec->impropers_allow == 0) error->one(FLERR,"Invalid data file section: AngleAngle Coeffs"); if (force->improper == NULL) error->one(FLERR,"Must define improper_style before AngleAngle Coeffs"); skip_lines(atom->nimpropertypes); } else if (strcmp(keyword,"Bonds") == 0) { for (i = 1; i < cmax; i++) count[i] = 0; if (force->newton_bond) for (i = 0; i < atom->nbonds; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(line,"%d %d %d %d",&tmp1,&tmp2,&atom1,&atom2); if (atom1 >= cmax) cmax = reallocate(&count,cmax,atom1); count[atom1]++; } else for (i = 0; i < atom->nbonds; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(line,"%d %d %d %d",&tmp1,&tmp2,&atom1,&atom2); int amax = MAX(atom1,atom2); if (amax >= cmax) cmax = reallocate(&count,cmax,amax); count[atom1]++; count[atom2]++; } for (i = 1; i < cmax; i++) bond_per_atom = MAX(bond_per_atom,count[i]); if (screen) fprintf(screen," %d = max bonds/atom\n",bond_per_atom); if (logfile) fprintf(logfile," %d = max bonds/atom\n",bond_per_atom); } else if (strcmp(keyword,"Angles") == 0) { for (i = 1; i < cmax; i++) count[i] = 0; if (force->newton_bond) for (i = 0; i < atom->nangles; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(line,"%d %d %d %d %d",&tmp1,&tmp2,&atom1,&atom2,&atom3); if (atom2 >= cmax) cmax = reallocate(&count,cmax,atom2); count[atom2]++; } else for (i = 0; i < atom->nangles; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(line,"%d %d %d %d %d",&tmp1,&tmp2,&atom1,&atom2,&atom3); int amax = MAX(atom1,atom2); amax = MAX(amax,atom3); if (amax >= cmax) cmax = reallocate(&count,cmax,amax); count[atom1]++; count[atom2]++; count[atom3]++; } for (i = 1; i < cmax; i++) angle_per_atom = MAX(angle_per_atom,count[i]); if (screen) fprintf(screen," %d = max angles/atom\n",angle_per_atom); if (logfile) fprintf(logfile," %d = max angles/atom\n",angle_per_atom); } else if (strcmp(keyword,"Dihedrals") == 0) { for (i = 1; i < cmax; i++) count[i] = 0; if (force->newton_bond) for (i = 0; i < atom->ndihedrals; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(line,"%d %d %d %d %d %d", &tmp1,&tmp2,&atom1,&atom2,&atom3,&atom4); if (atom2 >= cmax) cmax = reallocate(&count,cmax,atom2); count[atom2]++; } else for (i = 0; i < atom->ndihedrals; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(line,"%d %d %d %d %d %d", &tmp1,&tmp2,&atom1,&atom2,&atom3,&atom4); int amax = MAX(atom1,atom2); amax = MAX(amax,atom3); amax = MAX(amax,atom4); if (amax >= cmax) cmax = reallocate(&count,cmax,amax); count[atom1]++; count[atom2]++; count[atom3]++; count[atom4]++; } for (i = 1; i < cmax; i++) dihedral_per_atom = MAX(dihedral_per_atom,count[i]); if (screen) fprintf(screen," %d = max dihedrals/atom\n",dihedral_per_atom); if (logfile) fprintf(logfile," %d = max dihedrals/atom\n",dihedral_per_atom); } else if (strcmp(keyword,"Impropers") == 0) { for (i = 1; i < cmax; i++) count[i] = 0; if (force->newton_bond) for (i = 0; i < atom->nimpropers; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(line,"%d %d %d %d %d %d", &tmp1,&tmp2,&atom1,&atom2,&atom3,&atom4); if (atom2 >= cmax) cmax = reallocate(&count,cmax,atom2); count[atom2]++; } else for (i = 0; i < atom->nimpropers; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); sscanf(line,"%d %d %d %d %d %d", &tmp1,&tmp2,&atom1,&atom2,&atom3,&atom4); int amax = MAX(atom1,atom2); amax = MAX(amax,atom3); amax = MAX(amax,atom4); if (amax >= cmax) cmax = reallocate(&count,cmax,amax); count[atom1]++; count[atom2]++; count[atom3]++; count[atom4]++; } for (i = 1; i < cmax; i++) improper_per_atom = MAX(improper_per_atom,count[i]); if (screen) fprintf(screen," %d = max impropers/atom\n",improper_per_atom); if (logfile) fprintf(logfile," %d = max impropers/atom\n",improper_per_atom); } else { char str[128]; sprintf(str,"Unknown identifier in data file: %s",keyword); error->one(FLERR,str); } parse_keyword(0,0); } // free topology counting vector memory->destroy(count); // error check that topology was specified in file if ((atom->nbonds && !bond_per_atom) || (atom->nangles && !angle_per_atom) || (atom->ndihedrals && !dihedral_per_atom) || (atom->nimpropers && !improper_per_atom)) error->one(FLERR,"Needed topology not in data file"); // customize for new sections // error check that Bonus sections were speficied in file if (nellipsoids && !ellipsoid_flag) error->one(FLERR,"Needed bonus data not in data file"); if (nlines && !line_flag) error->one(FLERR,"Needed bonus data not in data file"); if (ntris && !tri_flag) error->one(FLERR,"Needed bonus data not in data file"); if (nbodies && !body_flag) error->one(FLERR,"Needed bonus data not in data file"); } /* ---------------------------------------------------------------------- reallocate the count vector from cmax to amax+1 and return new length zero new locations ------------------------------------------------------------------------- */ int ReadData::reallocate(int **pcount, int cmax, int amax) { int *count = *pcount; memory->grow(count,amax+1,"read_data:count"); for (int i = cmax; i <= amax; i++) count[i] = 0; *pcount = count; return amax+1; } /* ---------------------------------------------------------------------- proc 0 opens data file test if gzipped ------------------------------------------------------------------------- */ void ReadData::open(char *file) { compressed = 0; char *suffix = file + strlen(file) - 3; if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1; if (!compressed) fp = fopen(file,"r"); else { #ifdef LAMMPS_GZIP char gunzip[128]; sprintf(gunzip,"gunzip -c %s",file); fp = popen(gunzip,"r"); #else error->one(FLERR,"Cannot open gzipped file"); #endif } if (fp == NULL) { char str[128]; sprintf(str,"Cannot open file %s",file); error->one(FLERR,str); } } /* ---------------------------------------------------------------------- grab next keyword read lines until one is non-blank keyword is all text on line w/out leading & trailing white space read one additional line (assumed blank) if any read hits EOF, set keyword to empty if first = 1, line variable holds non-blank line that ended header if flag = 0, only proc 0 is calling so no bcast else flag = 1, bcast keyword line to all procs ------------------------------------------------------------------------- */ void ReadData::parse_keyword(int first, int flag) { int eof = 0; // proc 0 reads upto non-blank line plus 1 following line // eof is set to 1 if any read hits end-of-file if (me == 0) { if (!first) { 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(buffer,MAXLINE,fp) == NULL) eof = 1; } // if eof, set keyword empty and return if (flag) MPI_Bcast(&eof,1,MPI_INT,0,world); if (eof) { keyword[0] = '\0'; return; } // bcast keyword line to all procs if (flag) { 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]); } /* ---------------------------------------------------------------------- proc 0 reads N lines from file ------------------------------------------------------------------------- */ void ReadData::skip_lines(int n) { char *eof; for (int i = 0; i < n; i++) eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of data file"); } /* ---------------------------------------------------------------------- parse a line of coeffs into words, storing them in narg,arg trim anything from '#' onward word strings remain in line, are not copied if addstr != NULL, add addstr as extra arg for class2 angle/dihedral/improper if 2nd word starts with letter, then is hybrid style, add addstr after it else add addstr before 2nd word if dupflag, duplicate 1st word, so pair_coeff "2" becomes "2 2" ------------------------------------------------------------------------- */ void ReadData::parse_coeffs(char *line, const char *addstr, int dupflag) { char *ptr; if (ptr = strchr(line,'#')) *ptr = '\0'; narg = 0; char *word = strtok(line," \t\n\r\f"); while (word) { if (narg == maxarg) { maxarg += DELTA; arg = (char **) memory->srealloc(arg,maxarg*sizeof(char *),"read_data:arg"); } if (addstr && narg == 1 && !islower(word[0])) arg[narg++] = (char *) addstr; arg[narg++] = word; if (addstr && narg == 2 && islower(word[0])) arg[narg++] = (char *) addstr; if (dupflag && narg == 1) arg[narg++] = word; word = strtok(NULL," \t\n\r\f"); } } diff --git a/src/read_data.h b/src/read_data.h index 4788cdacf..b1d3bb48b 100644 --- a/src/read_data.h +++ b/src/read_data.h @@ -1,423 +1,424 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef COMMAND_CLASS CommandStyle(read_data,ReadData) #else #ifndef LMP_READ_DATA_H #define LMP_READ_DATA_H #include "stdio.h" #include "pointers.h" namespace LAMMPS_NS { class ReadData : protected Pointers { public: ReadData(class LAMMPS *); ~ReadData(); void command(int, char **); private: int me; char *line,*keyword,*buffer; FILE *fp; int narg,maxarg,compressed; char **arg; int nfix; // # of extra fixes that process/store info in data file int *fix_index; char **fix_header; char **fix_section; bigint nellipsoids; class AtomVecEllipsoid *avec_ellipsoid; bigint nlines; class AtomVecLine *avec_line; bigint ntris; class AtomVecTri *avec_tri; bigint nbodies; class AtomVecBody *avec_body; void open(char *); void scan(int &, int &, int &, int &); int reallocate(int **, int, int); void header(int); void parse_keyword(int, int); void skip_lines(int); void parse_coeffs(char *, const char *, int); void atoms(); void velocities(); void bonus(bigint, class AtomVec *, const char *); void bodies(); void bonds(); void angles(); void dihedrals(); void impropers(); void mass(); void paircoeffs(); + void pairIJcoeffs(); void bondcoeffs(); void anglecoeffs(int); void dihedralcoeffs(int); void impropercoeffs(int); void fix(int, char *, bigint); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Cannot read_data after simulation box is defined The read_data command cannot be used after a read_data, read_restart, or create_box command. E: Cannot run 2d simulation with nonperiodic Z dimension Use the boundary command to make the z dimension periodic in order to run a 2d simulation. E: Fix ID for read_data does not exist Self-explanatory. E: Must read Atoms before Velocities The Atoms section of a data file must come before a Velocities section. E: Invalid data file section: Ellipsoids Atom style does not allow ellipsoids. E: Must read Atoms before Ellipsoids The Atoms section of a data file must come before a Ellipsoids section. E: Invalid data file section: Lines Atom style does not allow lines. E: Must read Atoms before Lines The Atoms section of a data file must come before a Lines section. E: Invalid data file section: Triangles Atom style does not allow triangles. E: Must read Atoms before Triangles The Atoms section of a data file must come before a Triangles section. E: Invalid data file section: Bodies Atom style does not allow bodies. E: Must read Atoms before Bodies The Atoms section of a data file must come before a Bodies section. E: Invalid data file section: Bonds Atom style does not allow bonds. E: Must read Atoms before Bonds The Atoms section of a data file must come before a Bonds section. E: Invalid data file section: Angles Atom style does not allow angles. E: Must read Atoms before Angles The Atoms section of a data file must come before an Angles section. E: Invalid data file section: Dihedrals Atom style does not allow dihedrals. E: Must read Atoms before Dihedrals The Atoms section of a data file must come before a Dihedrals section. E: Invalid data file section: Impropers Atom style does not allow impropers. E: Must read Atoms before Impropers The Atoms section of a data file must come before an Impropers section. E: Must define pair_style before Pair Coeffs Must use a pair_style command before reading a data file that defines Pair Coeffs. E: Invalid data file section: Bond Coeffs Atom style does not allow bonds. E: Must define bond_style before Bond Coeffs Must use a bond_style command before reading a data file that defines Bond Coeffs. E: Invalid data file section: Angle Coeffs Atom style does not allow angles. E: Must define angle_style before Angle Coeffs Must use an angle_style command before reading a data file that defines Angle Coeffs. E: Invalid data file section: Dihedral Coeffs Atom style does not allow dihedrals. E: Must define dihedral_style before Dihedral Coeffs Must use a dihedral_style command before reading a data file that defines Dihedral Coeffs. E: Invalid data file section: Improper Coeffs Atom style does not allow impropers. E: Must define improper_style before Improper Coeffs Must use an improper_style command before reading a data file that defines Improper Coeffs. E: Invalid data file section: BondBond Coeffs Atom style does not allow angles. E: Must define angle_style before BondBond Coeffs Must use an angle_style command before reading a data file that defines Angle Coeffs. E: Invalid data file section: BondAngle Coeffs Atom style does not allow angles. E: Must define angle_style before BondAngle Coeffs Must use an angle_style command before reading a data file that defines Angle Coeffs. E: Invalid data file section: MiddleBondTorsion Coeffs Atom style does not allow dihedrals. E: Must define dihedral_style before MiddleBondTorsion Coeffs Must use a dihedral_style command before reading a data file that defines MiddleBondTorsion Coeffs. E: Invalid data file section: EndBondTorsion Coeffs Atom style does not allow dihedrals. E: Must define dihedral_style before EndBondTorsion Coeffs Must use a dihedral_style command before reading a data file that defines EndBondTorsion Coeffs. E: Invalid data file section: AngleTorsion Coeffs Atom style does not allow dihedrals. E: Must define dihedral_style before AngleTorsion Coeffs Must use a dihedral_style command before reading a data file that defines AngleTorsion Coeffs. E: Invalid data file section: AngleAngleTorsion Coeffs Atom style does not allow dihedrals. E: Must define dihedral_style before AngleAngleTorsion Coeffs Must use a dihedral_style command before reading a data file that defines AngleAngleTorsion Coeffs. E: Invalid data file section: BondBond13 Coeffs Atom style does not allow dihedrals. E: Must define dihedral_style before BondBond13 Coeffs Must use a dihedral_style command before reading a data file that defines BondBond13 Coeffs. E: Invalid data file section: AngleAngle Coeffs Atom style does not allow impropers. E: Must define improper_style before AngleAngle Coeffs Must use an improper_style command before reading a data file that defines AngleAngle Coeffs. E: Unknown identifier in data file: %s A section of the data file cannot be read by LAMMPS. E: No atoms in data file The header of the data file indicated that atoms would be included, but they were not present. E: Unexpected end of data file LAMMPS hit the end of the data file while attempting to read a section. Something is wrong with the format of the data file. E: No ellipsoids allowed with this atom style Self-explanatory. Check data file. E: No lines allowed with this atom style Self-explanatory. Check data file. E: No triangles allowed with this atom style Self-explanatory. Check data file. E: No bodies allowed with this atom style Self-explanatory. Check data file. E: System in data file is too big See the setting for bigint in the src/lmptype.h file. E: No bonds allowed with this atom style Self-explanatory. Check data file. E: No angles allowed with this atom style Self-explanatory. Check data file. E: No dihedrals allowed with this atom style Self-explanatory. Check data file. E: No impropers allowed with this atom style Self-explanatory. Check data file. E: Bonds defined but no bond types The data file header lists bonds but no bond types. E: Angles defined but no angle types The data file header lists angles but no angle types. E: Dihedrals defined but no dihedral types The data file header lists dihedrals but no dihedral types. E: Impropers defined but no improper types The data file header lists improper but no improper types. E: Did not assign all atoms correctly Atoms read in from a data file were not assigned correctly to processors. This is likely due to some atom coordinates being outside a non-periodic simulation box. E: Invalid atom ID in Atoms section of data file Atom IDs must be positive integers. E: Too many lines in one body in data file - boost MAXBODY MAXBODY is a setting at the top of the src/read_data.cpp file. Set it larger and re-compile the code. E: Bonds assigned incorrectly Bonds read in from the data file were not assigned correctly to atoms. This means there is something invalid about the topology definitions. E: Angles assigned incorrectly Angles read in from the data file were not assigned correctly to atoms. This means there is something invalid about the topology definitions. E: Dihedrals assigned incorrectly Dihedrals read in from the data file were not assigned correctly to atoms. This means there is something invalid about the topology definitions. E: Impropers assigned incorrectly Impropers read in from the data file were not assigned correctly to atoms. This means there is something invalid about the topology definitions. E: Molecular data file has too many atoms These kids of data files are currently limited to a number of atoms that fits in a 32-bit integer. E: Needed topology not in data file The header of the data file indicated that bonds or angles or dihedrals or impropers would be included, but they were not present. E: Needed bonus data not in data file Some atom styles require bonus data. See the read_data doc page for details. E: Cannot open gzipped file LAMMPS is attempting to open a gzipped version of the specified file but was unsuccessful. Check that the path and name are correct. E: Cannot open file %s The specified file cannot be opened. Check that the path and name are correct. */ diff --git a/src/write_data.cpp b/src/write_data.cpp index cfb1f7855..0f3de354e 100644 --- a/src/write_data.cpp +++ b/src/write_data.cpp @@ -1,602 +1,624 @@ /* ---------------------------------------------------------------------- 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 "lmptype.h" #include "mpi.h" #include "string.h" #include "write_data.h" #include "atom.h" #include "atom_vec.h" #include "group.h" #include "force.h" #include "pair.h" #include "bond.h" #include "angle.h" #include "dihedral.h" #include "improper.h" #include "update.h" #include "modify.h" #include "domain.h" #include "universe.h" #include "comm.h" #include "output.h" #include "thermo.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{IGNORE,WARN,ERROR}; // same as thermo.cpp +enum{II,IJ}; /* ---------------------------------------------------------------------- */ WriteData::WriteData(LAMMPS *lmp) : Pointers(lmp) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); } /* ---------------------------------------------------------------------- called as write_data command in input script ------------------------------------------------------------------------- */ void WriteData::command(int narg, char **arg) { if (domain->box_exist == 0) error->all(FLERR,"Write_data command before simulation box is defined"); - if (narg != 1) error->all(FLERR,"Illegal write_data command"); + + if (narg < 1) error->all(FLERR,"Illegal write_data command"); // if filename contains a "*", replace with current timestep char *ptr; int n = strlen(arg[0]) + 16; char *file = new char[n]; if (ptr = strchr(arg[0],'*')) { *ptr = '\0'; sprintf(file,"%s" BIGINT_FORMAT "%s",arg[0],update->ntimestep,ptr+1); } else strcpy(file,arg[0]); + // read optional args + + pairflag = II; + + int iarg = 1; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal write_data command"); + if (strcmp(arg[iarg+1],"ii") == 0) pairflag = II; + else if (strcmp(arg[iarg+1],"ij") == 0) pairflag = IJ; + else error->all(FLERR,"Illegal write_data command"); + iarg += 2; + } else error->all(FLERR,"Illegal write_data command"); + } + // init entire system since comm->exchange is done // comm::init needs neighbor::init needs pair::init needs kspace::init, etc if (comm->me == 0 && screen) fprintf(screen,"System init for write_data ...\n"); lmp->init(); // move atoms to new processors before writing file // do setup_pre_exchange to force update of per-atom info if needed // enforce PBC in case atoms are outside box // call borders() to rebuild atom map since exchange() destroys map modify->setup_pre_exchange(); if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); comm->exchange(); comm->borders(); if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); write(file); delete [] file; } /* ---------------------------------------------------------------------- called from command() later might let it be directly called within run/minimize loop ------------------------------------------------------------------------- */ void WriteData::write(char *file) { // special case where reneighboring is not done in integrator // on timestep data file is written (due to build_once being set) // if box is changing, must be reset, else data file will have // wrong box size and atoms will be lost when data file is read // other calls to pbc and domain and comm are not made, // b/c they only make sense if reneighboring is actually performed //if (neighbor->build_once) domain->reset_box(); // natoms = sum of nlocal = value to write into data file // if unequal and thermo lostflag is "error", don't write data file bigint nblocal = atom->nlocal; bigint natoms; MPI_Allreduce(&nblocal,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); if (natoms != atom->natoms && output->thermo->lostflag == ERROR) error->all(FLERR,"Atom count is inconsistent, cannot write data file"); // sum up bond,angle counts // may be different than atom->nbonds,nangles if broken/turned-off if (atom->nbonds || atom->nbondtypes) { nbonds_local = atom->avec->pack_bond(NULL); MPI_Allreduce(&nbonds_local,&nbonds,1,MPI_LMP_BIGINT,MPI_SUM,world); } if (atom->nangles || atom->nangletypes) { nangles_local = atom->avec->pack_angle(NULL); MPI_Allreduce(&nangles_local,&nangles,1,MPI_LMP_BIGINT,MPI_SUM,world); } // open data file if (me == 0) { fp = fopen(file,"w"); if (fp == NULL) { char str[128]; sprintf(str,"Cannot open data file %s",file); error->one(FLERR,str); } } // proc 0 writes header, ntype-length arrays, force fields if (me == 0) { header(); type_arrays(); force_fields(); } // per atom info if (natoms) atoms(); if (natoms) velocities(); if (atom->nbonds && nbonds) bonds(); if (atom->nangles && nangles) angles(); if (atom->ndihedrals) dihedrals(); if (atom->nimpropers) impropers(); // close data file if (me == 0) fclose(fp); } /* ---------------------------------------------------------------------- proc 0 writes out data file header ------------------------------------------------------------------------- */ void WriteData::header() { fprintf(fp,"LAMMPS data file via write_data, version %s, " "timestep = " BIGINT_FORMAT "\n", universe->version,update->ntimestep); fprintf(fp,"\n"); fprintf(fp,BIGINT_FORMAT " atoms\n",atom->natoms); fprintf(fp,"%d atom types\n",atom->ntypes); if (atom->nbonds || atom->nbondtypes) { fprintf(fp,BIGINT_FORMAT " bonds\n",nbonds); fprintf(fp,"%d bond types\n",atom->nbondtypes); } if (atom->nangles || atom->nangletypes) { fprintf(fp,BIGINT_FORMAT " angles\n",nangles); fprintf(fp,"%d angle types\n",atom->nangletypes); } if (atom->ndihedrals || atom->ndihedraltypes) { fprintf(fp,BIGINT_FORMAT " dihedrals\n",atom->ndihedrals); fprintf(fp,"%d dihedral types\n",atom->ndihedraltypes); } if (atom->nimpropers || atom->nimpropertypes) { fprintf(fp,BIGINT_FORMAT " impropers\n",atom->nimpropers); fprintf(fp,"%d improper types\n",atom->nimpropertypes); } fprintf(fp,"\n"); fprintf(fp,"%g %g xlo xhi\n",domain->boxlo[0],domain->boxhi[0]); fprintf(fp,"%g %g ylo yhi\n",domain->boxlo[1],domain->boxhi[1]); fprintf(fp,"%g %g zlo zhi\n",domain->boxlo[2],domain->boxhi[2]); if (domain->triclinic) fprintf(fp,"%g %g %g xy xz yz\n",domain->xy,domain->xz,domain->yz); } /* ---------------------------------------------------------------------- proc 0 writes out any type-based arrays that are defined ------------------------------------------------------------------------- */ void WriteData::type_arrays() { if (atom->mass) { double *mass = atom->mass; fprintf(fp,"\nMasses\n\n"); for (int i = 1; i <= atom->ntypes; i++) fprintf(fp,"%d %g\n",i,mass[i]); } } /* ---------------------------------------------------------------------- proc 0 writes out force field info ------------------------------------------------------------------------- */ void WriteData::force_fields() { if (force->pair && force->pair->writedata) { - fprintf(fp,"\nPair Coeffs\n\n"); - force->pair->write_data(fp); + if (pairflag == II) { + fprintf(fp,"\nPair Coeffs\n\n"); + force->pair->write_data(fp); + } else if (pairflag == IJ) { + fprintf(fp,"\nPairIJ Coeffs\n\n"); + force->pair->write_data_all(fp); + } } if (atom->avec->bonds_allow && force->bond && force->bond->writedata) { fprintf(fp,"\nBond Coeffs\n\n"); force->bond->write_data(fp); } if (atom->avec->angles_allow && force->angle && force->angle->writedata) { fprintf(fp,"\nAngle Coeffs\n\n"); force->angle->write_data(fp); } if (atom->avec->dihedrals_allow && force->dihedral && force->dihedral->writedata) { fprintf(fp,"\nDihedral Coeffs\n\n"); force->dihedral->write_data(fp); } if (atom->avec->impropers_allow && force->improper && force->improper->writedata) { fprintf(fp,"\nImproper Coeffs\n\n"); force->improper->write_data(fp); } } /* ---------------------------------------------------------------------- write out Atoms section of data file ------------------------------------------------------------------------- */ void WriteData::atoms() { // communication buffer for all my Atom info // max_size = largest buffer needed by any proc int ncol = atom->avec->size_data_atom + 3; int sendrow = atom->nlocal; int maxrow; MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world); double **buf; if (me == 0) memory->create(buf,maxrow,ncol,"write_data:buf"); else memory->create(buf,sendrow,ncol,"write_data:buf"); // pack my atom data into buf atom->avec->pack_data(buf); // write one chunk of atoms per proc to file // proc 0 pings each proc, receives its chunk, writes to file // all other procs wait for ping, send their chunk to proc 0 int tmp,recvrow; MPI_Status status; MPI_Request request; if (me == 0) { fprintf(fp,"\nAtoms\n\n"); for (int iproc = 0; iproc < nprocs; iproc++) { if (iproc) { MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_DOUBLE,iproc,0,world,&request); MPI_Send(&tmp,0,MPI_INT,iproc,0,world); MPI_Wait(&request,&status); MPI_Get_count(&status,MPI_DOUBLE,&recvrow); recvrow /= ncol; } else recvrow = sendrow; atom->avec->write_data(fp,recvrow,buf); } } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status); MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_DOUBLE,0,0,world); } memory->destroy(buf); } /* ---------------------------------------------------------------------- write out Velocities section of data file ------------------------------------------------------------------------- */ void WriteData::velocities() { // communication buffer for all my Atom info // max_size = largest buffer needed by any proc int ncol = atom->avec->size_velocity + 1; int sendrow = atom->nlocal; int maxrow; MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world); double **buf; if (me == 0) memory->create(buf,maxrow,ncol,"write_data:buf"); else memory->create(buf,sendrow,ncol,"write_data:buf"); // pack my velocity data into buf atom->avec->pack_vel(buf); // write one chunk of velocities per proc to file // proc 0 pings each proc, receives its chunk, writes to file // all other procs wait for ping, send their chunk to proc 0 int tmp,recvrow; MPI_Status status; MPI_Request request; if (me == 0) { fprintf(fp,"\nVelocities\n\n"); for (int iproc = 0; iproc < nprocs; iproc++) { if (iproc) { MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_DOUBLE,iproc,0,world,&request); MPI_Send(&tmp,0,MPI_INT,iproc,0,world); MPI_Wait(&request,&status); MPI_Get_count(&status,MPI_DOUBLE,&recvrow); recvrow /= ncol; } else recvrow = sendrow; atom->avec->write_vel(fp,recvrow,buf); } } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status); MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_DOUBLE,0,0,world); } memory->destroy(buf); } /* ---------------------------------------------------------------------- write out Bonds section of data file ------------------------------------------------------------------------- */ void WriteData::bonds() { // communication buffer for all my Bond info int ncol = 3; int sendrow = static_cast (nbonds_local); int maxrow; MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world); int **buf; if (me == 0) memory->create(buf,maxrow,ncol,"write_data:buf"); else memory->create(buf,sendrow,ncol,"write_data:buf"); // pack my bond data into buf int foo = atom->avec->pack_bond(buf); // write one chunk of atoms per proc to file // proc 0 pings each proc, receives its chunk, writes to file // all other procs wait for ping, send their chunk to proc 0 int tmp,recvrow; MPI_Status status; MPI_Request request; int index = 1; if (me == 0) { fprintf(fp,"\nBonds\n\n"); for (int iproc = 0; iproc < nprocs; iproc++) { if (iproc) { MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_INT,iproc,0,world,&request); MPI_Send(&tmp,0,MPI_INT,iproc,0,world); MPI_Wait(&request,&status); MPI_Get_count(&status,MPI_INT,&recvrow); recvrow /= ncol; } else recvrow = sendrow; atom->avec->write_bond(fp,recvrow,buf,index); index += recvrow; } } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status); MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_INT,0,0,world); } memory->destroy(buf); } /* ---------------------------------------------------------------------- write out Angles section of data file ------------------------------------------------------------------------- */ void WriteData::angles() { // communication buffer for all my Angle info int ncol = 4; int sendrow = static_cast (nangles_local); int maxrow; MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world); int **buf; if (me == 0) memory->create(buf,maxrow,ncol,"write_data:buf"); else memory->create(buf,sendrow,ncol,"write_data:buf"); // pack my angle data into buf atom->avec->pack_angle(buf); // write one chunk of atoms per proc to file // proc 0 pings each proc, receives its chunk, writes to file // all other procs wait for ping, send their chunk to proc 0 int tmp,recvrow; MPI_Status status; MPI_Request request; int index = 1; if (me == 0) { fprintf(fp,"\nAngles\n\n"); for (int iproc = 0; iproc < nprocs; iproc++) { if (iproc) { MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_INT,iproc,0,world,&request); MPI_Send(&tmp,0,MPI_INT,iproc,0,world); MPI_Wait(&request,&status); MPI_Get_count(&status,MPI_INT,&recvrow); recvrow /= ncol; } else recvrow = sendrow; atom->avec->write_angle(fp,recvrow,buf,index); index += recvrow; } } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status); MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_INT,0,0,world); } memory->destroy(buf); } /* ---------------------------------------------------------------------- write out Dihedrals section of data file ------------------------------------------------------------------------- */ void WriteData::dihedrals() { // communication buffer for all my Dihedral info // max_size = largest buffer needed by any proc int ncol = 5; int *tag = atom->tag; int *num_dihedral = atom->num_dihedral; int **dihedral_atom2 = atom->dihedral_atom2; int nlocal = atom->nlocal; int newton_bond = force->newton_bond; int i,j; int sendrow = 0; if (newton_bond) { for (i = 0; i < nlocal; i++) sendrow += num_dihedral[i]; } else { for (i = 0; i < nlocal; i++) for (j = 0; j < num_dihedral[i]; j++) if (tag[i] == dihedral_atom2[i][j]) sendrow++; } int maxrow; MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world); int **buf; if (me == 0) memory->create(buf,maxrow,ncol,"write_data:buf"); else memory->create(buf,sendrow,ncol,"write_data:buf"); // pack my dihedral data into buf atom->avec->pack_dihedral(buf); // write one chunk of atoms per proc to file // proc 0 pings each proc, receives its chunk, writes to file // all other procs wait for ping, send their chunk to proc 0 int tmp,recvrow; MPI_Status status; MPI_Request request; int index = 1; if (me == 0) { fprintf(fp,"\nDihedrals\n\n"); for (int iproc = 0; iproc < nprocs; iproc++) { if (iproc) { MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_INT,iproc,0,world,&request); MPI_Send(&tmp,0,MPI_INT,iproc,0,world); MPI_Wait(&request,&status); MPI_Get_count(&status,MPI_INT,&recvrow); recvrow /= ncol; } else recvrow = sendrow; atom->avec->write_dihedral(fp,recvrow,buf,index); index += recvrow; } } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status); MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_INT,0,0,world); } memory->destroy(buf); } /* ---------------------------------------------------------------------- write out Impropers section of data file ------------------------------------------------------------------------- */ void WriteData::impropers() { // communication buffer for all my Improper info // max_size = largest buffer needed by any proc int ncol = 5; int *tag = atom->tag; int *num_improper = atom->num_improper; int **improper_atom2 = atom->improper_atom2; int nlocal = atom->nlocal; int newton_bond = force->newton_bond; int i,j; int sendrow = 0; if (newton_bond) { for (i = 0; i < nlocal; i++) sendrow += num_improper[i]; } else { for (i = 0; i < nlocal; i++) for (j = 0; j < num_improper[i]; j++) if (tag[i] == improper_atom2[i][j]) sendrow++; } int maxrow; MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world); int **buf; if (me == 0) memory->create(buf,maxrow,ncol,"write_data:buf"); else memory->create(buf,sendrow,ncol,"write_data:buf"); // pack my improper data into buf atom->avec->pack_improper(buf); // write one chunk of atoms per proc to file // proc 0 pings each proc, receives its chunk, writes to file // all other procs wait for ping, send their chunk to proc 0 int tmp,recvrow; MPI_Status status; MPI_Request request; int index = 1; if (me == 0) { fprintf(fp,"\nImpropers\n\n"); for (int iproc = 0; iproc < nprocs; iproc++) { if (iproc) { MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_INT,iproc,0,world,&request); MPI_Send(&tmp,0,MPI_INT,iproc,0,world); MPI_Wait(&request,&status); MPI_Get_count(&status,MPI_INT,&recvrow); recvrow /= ncol; } else recvrow = sendrow; atom->avec->write_improper(fp,recvrow,buf,index); index += recvrow; } } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status); MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_INT,0,0,world); } memory->destroy(buf); } diff --git a/src/write_data.h b/src/write_data.h index 7b3044192..e879b00d7 100644 --- a/src/write_data.h +++ b/src/write_data.h @@ -1,54 +1,55 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef COMMAND_CLASS CommandStyle(write_data,WriteData) #else #ifndef LMP_WRITE_DATA_H #define LMP_WRITE_DATA_H #include "stdio.h" #include "pointers.h" namespace LAMMPS_NS { class WriteData : protected Pointers { public: WriteData(class LAMMPS *); void command(int, char **); void write(char *); private: int me,nprocs; + int pairflag; FILE *fp; bigint nbonds_local,nbonds; bigint nangles_local,nangles; void header(); void type_arrays(); void force_fields(); void atoms(); void velocities(); void bonds(); void angles(); void dihedrals(); void impropers(); }; } #endif #endif