diff --git a/src/accelerator_kokkos.h b/src/accelerator_kokkos.h index dc8bf5d9b..2c542b74f 100644 --- a/src/accelerator_kokkos.h +++ b/src/accelerator_kokkos.h @@ -1,88 +1,94 @@ /* -*- 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_ACCELERATOR_KOKKOS_H #define LMP_ACCELERATOR_KOKKOS_H // true interface to KOKKOS // used when KOKKOS is installed #ifdef LMP_KOKKOS #include "kokkos.h" #include "atom_kokkos.h" #include "comm_kokkos.h" #include "domain_kokkos.h" #include "neighbor_kokkos.h" #include "modify_kokkos.h" #else // dummy interface to KOKKOS // needed for compiling when KOKKOS is not installed #include "atom.h" #include "comm_brick.h" #include "domain.h" #include "neighbor.h" #include "modify.h" namespace LAMMPS_NS { class KokkosLMP { public: int kokkos_exists; int num_threads; int numa; KokkosLMP(class LAMMPS *, int, char **) {kokkos_exists = 0;} ~KokkosLMP() {} void accelerator(int, char **) {} int neigh_list_kokkos(int) {return 0;} int neigh_count(int) {return 0;} }; class AtomKokkos : public Atom { public: AtomKokkos(class LAMMPS *lmp) : Atom(lmp) {} ~AtomKokkos() {} }; class CommKokkos : public CommBrick { public: CommKokkos(class LAMMPS *lmp) : CommBrick(lmp) {} ~CommKokkos() {} }; class DomainKokkos : public Domain { public: DomainKokkos(class LAMMPS *lmp) : Domain(lmp) {} ~DomainKokkos() {} }; class NeighborKokkos : public Neighbor { public: NeighborKokkos(class LAMMPS *lmp) : Neighbor(lmp) {} ~NeighborKokkos() {} }; class ModifyKokkos : public Modify { public: ModifyKokkos(class LAMMPS *lmp) : Modify(lmp) {} ~ModifyKokkos() {} }; +class DAT { + public: + typedef double tdual_xfloat_1d; + typedef int tdual_int_2d; +}; + } #endif #endif diff --git a/src/fix.cpp b/src/fix.cpp index 332e3c8ae..60e11de5a 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -1,203 +1,207 @@ /* ---------------------------------------------------------------------- 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 "string.h" #include "ctype.h" #include "fix.h" #include "atom.h" #include "group.h" #include "atom_masks.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; // allocate space for static class instance variable and initialize it int Fix::instance_total = 0; /* ---------------------------------------------------------------------- */ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) { instance_me = instance_total++; // fix ID, group, and style // ID must be all alphanumeric chars or underscores int n = strlen(arg[0]) + 1; id = new char[n]; strcpy(id,arg[0]); for (int i = 0; i < n-1; i++) if (!isalnum(id[i]) && id[i] != '_') error->all(FLERR,"Fix ID must be alphanumeric or underscore characters"); igroup = group->find(arg[1]); if (igroup == -1) error->all(FLERR,"Could not find fix group ID"); groupbit = group->bitmask[igroup]; n = strlen(arg[2]) + 1; style = new char[n]; strcpy(style,arg[2]); restart_global = restart_peratom = restart_file = 0; force_reneighbor = 0; box_change_size = box_change_shape = box_change_domain = 0; thermo_energy = 0; rigid_flag = 0; virial_flag = 0; no_change_box = 0; time_integrate = 0; time_depend = 0; create_attribute = 0; restart_pbc = 0; wd_header = wd_section = 0; dynamic_group_allow = 0; cudable_comm = 0; scalar_flag = vector_flag = array_flag = 0; peratom_flag = local_flag = 0; size_vector_variable = size_array_rows_variable = 0; comm_forward = comm_reverse = comm_border = 0; restart_reset = 0; // reasonable defaults // however, each fix that uses these values should explicitly set them nevery = 1; global_freq = 1; maxvatom = 0; vatom = NULL; // CUDA and KOKKOS per-fix data masks datamask = ALL_MASK; datamask_ext = ALL_MASK; execution_space = Host; datamask_read = ALL_MASK; datamask_modify = ALL_MASK; + + copymode = 0; } /* ---------------------------------------------------------------------- */ Fix::~Fix() { - delete [] id; - delete [] style; - memory->destroy(vatom); + if (!copymode) { + delete [] id; + delete [] style; + memory->destroy(vatom); + } } /* ---------------------------------------------------------------------- process params common to all fixes here if unknown param, call modify_param specific to the fix ------------------------------------------------------------------------- */ void Fix::modify_params(int narg, char **arg) { if (narg == 0) error->all(FLERR,"Illegal fix_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"energy") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command"); if (strcmp(arg[iarg+1],"no") == 0) thermo_energy = 0; else if (strcmp(arg[iarg+1],"yes") == 0) thermo_energy = 1; else error->all(FLERR,"Illegal fix_modify command"); iarg += 2; } else { int n = modify_param(narg-iarg,&arg[iarg]); if (n == 0) error->all(FLERR,"Illegal fix_modify command"); iarg += n; } } } /* ---------------------------------------------------------------------- setup for virial computation see integrate::ev_set() for values of vflag (0-6) ------------------------------------------------------------------------- */ void Fix::v_setup(int vflag) { int i,n; evflag = 1; vflag_global = vflag % 4; vflag_atom = vflag / 4; // reallocate per-atom array if necessary if (vflag_atom && atom->nlocal > maxvatom) { maxvatom = atom->nmax; memory->destroy(vatom); memory->create(vatom,maxvatom,6,"fix:vatom"); } // zero accumulators if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; if (vflag_atom) { n = atom->nlocal; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; vatom[i][1] = 0.0; vatom[i][2] = 0.0; vatom[i][3] = 0.0; vatom[i][4] = 0.0; vatom[i][5] = 0.0; } } } /* ---------------------------------------------------------------------- tally virial into global and per-atom accumulators v = total virial for the interaction involving total atoms n = # of local atoms involved, with local indices in list increment global virial by n/total fraction increment per-atom virial of each atom in list by 1/total fraction assumes other procs will tally left-over fractions ------------------------------------------------------------------------- */ void Fix::v_tally(int n, int *list, double total, double *v) { int m; if (vflag_global) { double fraction = n/total; virial[0] += fraction*v[0]; virial[1] += fraction*v[1]; virial[2] += fraction*v[2]; virial[3] += fraction*v[3]; virial[4] += fraction*v[4]; virial[5] += fraction*v[5]; } if (vflag_atom) { double fraction = 1.0/total; for (int i = 0; i < n; i++) { m = list[i]; vatom[m][0] += fraction*v[0]; vatom[m][1] += fraction*v[1]; vatom[m][2] += fraction*v[2]; vatom[m][3] += fraction*v[3]; vatom[m][4] += fraction*v[4]; vatom[m][5] += fraction*v[5]; } } } diff --git a/src/fix.h b/src/fix.h index b762f9902..56ff4d0fd 100644 --- a/src/fix.h +++ b/src/fix.h @@ -1,271 +1,274 @@ /* -*- 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_FIX_H #define LMP_FIX_H #include "pointers.h" namespace LAMMPS_NS { class Fix : protected Pointers { public: static int instance_total; // # of Fix classes ever instantiated char *id,*style; int igroup,groupbit; int restart_global; // 1 if Fix saves global state, 0 if not int restart_peratom; // 1 if Fix saves peratom state, 0 if not int restart_file; // 1 if Fix writes own restart file, 0 if not int force_reneighbor; // 1 if Fix forces reneighboring, 0 if not int box_change_size; // 1 if Fix changes box size, 0 if not int box_change_shape; // 1 if Fix changes box shape, 0 if not int box_change_domain; // 1 if Fix changes proc sub-domains, 0 if not bigint next_reneighbor; // next timestep to force a reneighboring int thermo_energy; // 1 if fix_modify enabled ThEng, 0 if not int nevery; // how often to call an end_of_step fix int rigid_flag; // 1 if Fix integrates rigid bodies, 0 if not int virial_flag; // 1 if Fix contributes to virial, 0 if not int no_change_box; // 1 if cannot swap ortho <-> triclinic int time_integrate; // 1 if fix performs time integration, 0 if no int time_depend; // 1 if requires continuous timestepping int create_attribute; // 1 if fix stores attributes that need // setting when a new atom is created int restart_pbc; // 1 if fix moves atoms (except integrate) // so write_restart must remap to PBC int wd_header; // # of header values fix writes to data file int wd_section; // # of sections fix writes to data file int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 int cudable_comm; // 1 if fix has CUDA-enabled communication int scalar_flag; // 0/1 if compute_scalar() function exists int vector_flag; // 0/1 if compute_vector() function exists int array_flag; // 0/1 if compute_array() function exists int size_vector; // length of global vector int size_array_rows; // rows in global array int size_array_cols; // columns in global array int size_vector_variable; // 1 if vec length is unknown in advance int size_array_rows_variable; // 1 if array rows is unknown in advance int global_freq; // frequency s/v data is available at int peratom_flag; // 0/1 if per-atom data is stored int size_peratom_cols; // 0 = vector, N = columns in peratom array int peratom_freq; // frequency per-atom data is available at int local_flag; // 0/1 if local data is stored int size_local_rows; // rows in local vector or array int size_local_cols; // 0 = vector, N = columns in local array int local_freq; // frequency local data is available at int extscalar; // 0/1 if global scalar is intensive/extensive int extvector; // 0/1/-1 if global vector is all int/ext/extlist int *extlist; // list of 0/1 int/ext for each vec component int extarray; // 0/1 if global array is intensive/extensive double *vector_atom; // computed per-atom vector double **array_atom; // computed per-atom array double *vector_local; // computed local vector double **array_local; // computed local array int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int comm_border; // size of border communication (0 if none) double virial[6]; // accumlated virial double **vatom; // accumulated per-atom virial int restart_reset; // 1 if restart just re-initialized fix // KOKKOS host/device flag and data masks ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; // USER-CUDA per-fix data masks unsigned int datamask; unsigned int datamask_ext; Fix(class LAMMPS *, int, char **); virtual ~Fix(); void modify_params(int, char **); virtual int setmask() = 0; virtual void post_constructor() {} virtual void init() {} virtual void init_list(int, class NeighList *) {} virtual void setup(int) {} virtual void setup_pre_exchange() {} virtual void setup_pre_neighbor() {} virtual void setup_pre_force(int) {} virtual void min_setup(int) {} virtual void initial_integrate(int) {} virtual void post_integrate() {} virtual void pre_exchange() {} virtual void pre_neighbor() {} virtual void pre_force(int) {} virtual void post_force(int) {} virtual void final_integrate() {} virtual void end_of_step() {} virtual void post_run() {} virtual void write_restart(FILE *) {} virtual void write_restart_file(char *) {} virtual void restart(char *) {} virtual void grow_arrays(int) {} virtual void copy_arrays(int, int, int) {} virtual void set_arrays(int) {} virtual void update_arrays(int, int) {} virtual void set_molecule(int, tagint, int, double *, double *, double *) {} virtual int pack_border(int, int *, double *) {return 0;} virtual int unpack_border(int, int, double *) {return 0;} virtual int pack_exchange(int, double *) {return 0;} virtual int unpack_exchange(int, double *) {return 0;} virtual int pack_restart(int, double *) {return 0;} virtual void unpack_restart(int, int) {} virtual int size_restart(int) {return 0;} virtual int maxsize_restart() {return 0;} virtual void setup_pre_force_respa(int, int) {} virtual void initial_integrate_respa(int, int, int) {} virtual void post_integrate_respa(int, int) {} virtual void pre_force_respa(int, int, int) {} virtual void post_force_respa(int, int, int) {} virtual void final_integrate_respa(int, int) {} virtual void min_setup_pre_exchange() {} virtual void min_setup_pre_neighbor() {} virtual void min_setup_pre_force(int) {} virtual void min_pre_exchange() {} virtual void min_pre_neighbor() {} virtual void min_pre_force(int) {} virtual void min_post_force(int) {} virtual double min_energy(double *) {return 0.0;} virtual void min_store() {} virtual void min_clearstore() {} virtual void min_pushstore() {} virtual void min_popstore() {} virtual int min_reset_ref() {return 0;} virtual void min_step(double, double *) {} virtual double max_alpha(double *) {return 0.0;} virtual int min_dof() {return 0;} virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} virtual double compute_scalar() {return 0.0;} virtual double compute_vector(int) {return 0.0;} virtual double compute_array(int,int) {return 0.0;} virtual int dof(int) {return 0;} virtual void deform(int) {} virtual void reset_target(double) {} virtual void reset_dt() {} virtual void reset_timestep(bigint) {} virtual void read_data_header(char *) {} virtual void read_data_section(char *, int, char *) {} virtual bigint read_data_skip_lines(char *) {return 0;} virtual void write_data_header(FILE *, int) {} virtual void write_data_section_size(int, int &, int &) {} virtual void write_data_section_pack(int, double **) {} virtual void write_data_section_keyword(int, FILE *) {} virtual void write_data_section(int, FILE *, int, double **, int) {} virtual void zero_momentum() {} virtual void zero_rotation() {} virtual int modify_param(int, char **) {return 0;} virtual void *extract(const char *, int &) {return NULL;} virtual double memory_usage() {return 0.0;} virtual unsigned int data_mask() {return datamask;} virtual unsigned int data_mask_ext() {return datamask_ext;} protected: int instance_me; // which Fix class instantiation I am int evflag; int vflag_global,vflag_atom; int maxvatom; + int copymode; // if set, do not deallocate during destruction + // required when classes are used as functors by Kokkos + void v_setup(int); void v_tally(int, int *, double, double *); // union data struct for packing 32-bit and 64-bit ints into double bufs // see atom_vec.h for documentation union ubuf { double d; int64_t i; ubuf(double arg) : d(arg) {} ubuf(int64_t arg) : i(arg) {} ubuf(int arg) : i(arg) {} }; }; namespace FixConst { static const int INITIAL_INTEGRATE = 1<<0; static const int POST_INTEGRATE = 1<<1; static const int PRE_EXCHANGE = 1<<2; static const int PRE_NEIGHBOR = 1<<3; static const int PRE_FORCE = 1<<4; static const int POST_FORCE = 1<<5; static const int FINAL_INTEGRATE = 1<<6; static const int END_OF_STEP = 1<<7; static const int THERMO_ENERGY = 1<<8; static const int INITIAL_INTEGRATE_RESPA = 1<<9; static const int POST_INTEGRATE_RESPA = 1<<10; static const int PRE_FORCE_RESPA = 1<<11; static const int POST_FORCE_RESPA = 1<<12; static const int FINAL_INTEGRATE_RESPA = 1<<13; static const int MIN_PRE_EXCHANGE = 1<<14; static const int MIN_PRE_NEIGHBOR = 1<<15; static const int MIN_PRE_FORCE = 1<<16; static const int MIN_POST_FORCE = 1<<17; static const int MIN_ENERGY = 1<<18; static const int POST_RUN = 1<<19; static const int FIX_CONST_LAST = 1<<20; } } #endif /* ERROR/WARNING messages: E: Fix ID must be alphanumeric or underscore characters Self-explanatory. E: Could not find fix group ID A group ID used in the fix command does not exist. 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. */ diff --git a/src/pair.cpp b/src/pair.cpp index 42ad10032..27811567c 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -1,1646 +1,1650 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "ctype.h" #include "float.h" #include "limits.h" #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair.h" #include "atom.h" #include "neighbor.h" #include "neigh_list.h" #include "domain.h" #include "comm.h" #include "force.h" #include "kspace.h" #include "update.h" #include "accelerator_cuda.h" #include "suffix.h" #include "atom_masks.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define EWALD_F 1.12837917 enum{NONE,RLINEAR,RSQ,BMP}; // allocate space for static class instance variable and initialize it int Pair::instance_total = 0; /* ---------------------------------------------------------------------- */ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) { instance_me = instance_total++; THIRD = 1.0/3.0; eng_vdwl = eng_coul = 0.0; comm_forward = comm_reverse = comm_reverse_off = 0; single_enable = 1; restartinfo = 1; respa_enable = 0; one_coeff = 0; no_virial_fdotr_compute = 0; writedata = 0; ghostneigh = 0; nextra = 0; pvector = NULL; single_extra = 0; svector = NULL; ewaldflag = pppmflag = msmflag = dispersionflag = tip4pflag = dipoleflag = 0; reinitflag = 1; // pair_modify settingsx compute_flag = 1; manybody_flag = 0; offset_flag = 0; mix_flag = GEOMETRIC; tail_flag = 0; etail = ptail = etail_ij = ptail_ij = 0.0; ncoultablebits = 12; ndisptablebits = 12; tabinner = sqrt(2.0); tabinner_disp = sqrt(2.0); allocated = 0; suffix_flag = Suffix::NONE; maxeatom = maxvatom = 0; eatom = NULL; vatom = NULL; // CUDA and KOKKOS per-fix data masks datamask = ALL_MASK; datamask_ext = ALL_MASK; execution_space = Host; datamask_read = ALL_MASK; datamask_modify = ALL_MASK; + + copymode = 0; } /* ---------------------------------------------------------------------- */ Pair::~Pair() { - memory->destroy(eatom); - memory->destroy(vatom); + if (!copymode) { + memory->destroy(eatom); + memory->destroy(vatom); + } } /* ---------------------------------------------------------------------- modify parameters of the pair style pair_hybrid has its own version of this routine to apply modifications to each of its sub-styles ------------------------------------------------------------------------- */ void Pair::modify_params(int narg, char **arg) { if (narg == 0) error->all(FLERR,"Illegal pair_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"mix") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); if (strcmp(arg[iarg+1],"geometric") == 0) mix_flag = GEOMETRIC; else if (strcmp(arg[iarg+1],"arithmetic") == 0) mix_flag = ARITHMETIC; else if (strcmp(arg[iarg+1],"sixthpower") == 0) mix_flag = SIXTHPOWER; else error->all(FLERR,"Illegal pair_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"shift") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); if (strcmp(arg[iarg+1],"yes") == 0) offset_flag = 1; else if (strcmp(arg[iarg+1],"no") == 0) offset_flag = 0; else error->all(FLERR,"Illegal pair_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"table") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); ncoultablebits = force->inumeric(FLERR,arg[iarg+1]); if (ncoultablebits > sizeof(float)*CHAR_BIT) error->all(FLERR,"Too many total bits for bitmapped lookup table"); iarg += 2; } else if (strcmp(arg[iarg],"table/disp") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); ndisptablebits = force->inumeric(FLERR,arg[iarg+1]); if (ndisptablebits > sizeof(float)*CHAR_BIT) error->all(FLERR,"Too many total bits for bitmapped lookup table"); iarg += 2; } else if (strcmp(arg[iarg],"tabinner") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); tabinner = force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"tabinner/disp") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); tabinner_disp = force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"tail") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); if (strcmp(arg[iarg+1],"yes") == 0) tail_flag = 1; else if (strcmp(arg[iarg+1],"no") == 0) tail_flag = 0; else error->all(FLERR,"Illegal pair_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"compute") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); if (strcmp(arg[iarg+1],"yes") == 0) compute_flag = 1; else if (strcmp(arg[iarg+1],"no") == 0) compute_flag = 0; else error->all(FLERR,"Illegal pair_modify command"); iarg += 2; } else error->all(FLERR,"Illegal pair_modify command"); } } /* ---------------------------------------------------------------------- */ void Pair::init() { int i,j; if (offset_flag && tail_flag) error->all(FLERR,"Cannot have both pair_modify shift and tail set to yes"); if (tail_flag && domain->dimension == 2) error->all(FLERR,"Cannot use pair tail corrections with 2d simulations"); if (tail_flag && domain->nonperiodic && comm->me == 0) error->warning(FLERR,"Using pair tail corrections with nonperiodic system"); // for manybody potentials // check if bonded exclusions could invalidate the neighbor list if (manybody_flag && atom->molecular) { int flag = 0; if (atom->nbonds > 0 && force->special_lj[1] == 0.0 && force->special_coul[1] == 0.0) flag = 1; if (atom->nangles > 0 && force->special_lj[2] == 0.0 && force->special_coul[2] == 0.0) flag = 1; if (atom->ndihedrals > 0 && force->special_lj[3] == 0.0 && force->special_coul[3] == 0.0) flag = 1; if (flag && comm->me == 0) error->warning(FLERR,"Using a manybody potential with " "bonds/angles/dihedrals and special_bond exclusions"); } // I,I coeffs must be set // init_one() will check if I,J is set explicitly or inferred by mixing if (!allocated) error->all(FLERR,"All pair coeffs are not set"); for (i = 1; i <= atom->ntypes; i++) if (setflag[i][i] == 0) error->all(FLERR,"All pair coeffs are not set"); // style-specific initialization init_style(); // call init_one() for each I,J // set cutsq for each I,J, used to neighbor // cutforce = max of all I,J cutoffs cutforce = 0.0; etail = ptail = 0.0; double cut; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { cut = init_one(i,j); cutsq[i][j] = cutsq[j][i] = cut*cut; cutforce = MAX(cutforce,cut); if (tail_flag) { etail += etail_ij; ptail += ptail_ij; if (i != j) { etail += etail_ij; ptail += ptail_ij; } } } } /* ---------------------------------------------------------------------- reset all type-based params by invoking init_one() for each I,J called by fix adapt after it changes one or more params ------------------------------------------------------------------------- */ void Pair::reinit() { // generalize this error message if reinit() is used by more than fix adapt if (!reinitflag) error->all(FLERR,"Fix adapt interface to this pair style not supported"); etail = ptail = 0.0; for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) { init_one(i,j); if (tail_flag) { etail += etail_ij; ptail += ptail_ij; if (i != j) { etail += etail_ij; ptail += ptail_ij; } } } } /* ---------------------------------------------------------------------- init specific to a pair style specific pair style can override this function if needs its own error checks if needs another kind of neighbor list request default neighbor list = half list ------------------------------------------------------------------------- */ void Pair::init_style() { neighbor->request(this,instance_me); } /* ---------------------------------------------------------------------- neighbor callback to inform pair style of neighbor list to use specific pair style can override this function ------------------------------------------------------------------------- */ void Pair::init_list(int which, NeighList *ptr) { list = ptr; } /* ---------------------------------------------------------------------- setup Coulomb force tables used in compute routines ------------------------------------------------------------------------- */ void Pair::init_tables(double cut_coul, double *cut_respa) { int masklo,maskhi; double r,grij,expm2,derfc,egamma,fgamma,rsw; double qqrd2e = force->qqrd2e; if (force->kspace == NULL) error->all(FLERR,"Pair style requres a KSpace style"); double g_ewald = force->kspace->g_ewald; double cut_coulsq = cut_coul * cut_coul; tabinnersq = tabinner*tabinner; init_bitmap(tabinner,cut_coul,ncoultablebits, masklo,maskhi,ncoulmask,ncoulshiftbits); int ntable = 1; for (int i = 0; i < ncoultablebits; i++) ntable *= 2; // linear lookup tables of length N = 2^ncoultablebits // stored value = value at lower edge of bin // d values = delta from lower edge to upper edge of bin if (ftable) free_tables(); memory->create(rtable,ntable,"pair:rtable"); memory->create(ftable,ntable,"pair:ftable"); memory->create(ctable,ntable,"pair:ctable"); memory->create(etable,ntable,"pair:etable"); memory->create(drtable,ntable,"pair:drtable"); memory->create(dftable,ntable,"pair:dftable"); memory->create(dctable,ntable,"pair:dctable"); memory->create(detable,ntable,"pair:detable"); if (cut_respa == NULL) { vtable = ptable = dvtable = dptable = NULL; } else { memory->create(vtable,ntable,"pair:vtable"); memory->create(ptable,ntable,"pair:ptable"); memory->create(dvtable,ntable,"pair:dvtable"); memory->create(dptable,ntable,"pair:dptable"); } union_int_float_t rsq_lookup; union_int_float_t minrsq_lookup; int itablemin; minrsq_lookup.i = 0 << ncoulshiftbits; minrsq_lookup.i |= maskhi; for (int i = 0; i < ntable; i++) { rsq_lookup.i = i << ncoulshiftbits; rsq_lookup.i |= masklo; if (rsq_lookup.f < tabinnersq) { rsq_lookup.i = i << ncoulshiftbits; rsq_lookup.i |= maskhi; } r = sqrtf(rsq_lookup.f); if (msmflag) { egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)* force->kspace->dgamma(r/cut_coul); } else { grij = g_ewald * r; expm2 = exp(-grij*grij); derfc = erfc(grij); } if (cut_respa == NULL) { rtable[i] = rsq_lookup.f; ctable[i] = qqrd2e/r; if (msmflag) { ftable[i] = qqrd2e/r * fgamma; etable[i] = qqrd2e/r * egamma; } else { ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); etable[i] = qqrd2e/r * derfc; } } else { rtable[i] = rsq_lookup.f; ctable[i] = 0.0; ptable[i] = qqrd2e/r; if (msmflag) { ftable[i] = qqrd2e/r * (fgamma - 1.0); etable[i] = qqrd2e/r * egamma; vtable[i] = qqrd2e/r * fgamma; } else { ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); etable[i] = qqrd2e/r * derfc; vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); } if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); } else { if (msmflag) ftable[i] = qqrd2e/r * fgamma; else ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); ctable[i] = qqrd2e/r; } } } minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); } tabinnersq = minrsq_lookup.f; int ntablem1 = ntable - 1; for (int i = 0; i < ntablem1; i++) { drtable[i] = 1.0/(rtable[i+1] - rtable[i]); dftable[i] = ftable[i+1] - ftable[i]; dctable[i] = ctable[i+1] - ctable[i]; detable[i] = etable[i+1] - etable[i]; } if (cut_respa) { for (int i = 0; i < ntablem1; i++) { dvtable[i] = vtable[i+1] - vtable[i]; dptable[i] = ptable[i+1] - ptable[i]; } } // get the delta values for the last table entries // tables are connected periodically between 0 and ntablem1 drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); dftable[ntablem1] = ftable[0] - ftable[ntablem1]; dctable[ntablem1] = ctable[0] - ctable[ntablem1]; detable[ntablem1] = etable[0] - etable[ntablem1]; if (cut_respa) { dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; dptable[ntablem1] = ptable[0] - ptable[ntablem1]; } // get the correct delta values at itablemax // smallest r is in bin itablemin // largest r is in bin itablemax, which is itablemin-1, // or ntablem1 if itablemin=0 // deltas at itablemax only needed if corresponding rsq < cut*cut // if so, compute deltas between rsq and cut*cut double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp; p_tmp = 0.0; v_tmp = 0.0; itablemin = minrsq_lookup.i & ncoulmask; itablemin >>= ncoulshiftbits; int itablemax = itablemin - 1; if (itablemin == 0) itablemax = ntablem1; rsq_lookup.i = itablemax << ncoulshiftbits; rsq_lookup.i |= maskhi; if (rsq_lookup.f < cut_coulsq) { rsq_lookup.f = cut_coulsq; r = sqrtf(rsq_lookup.f); if (msmflag) { egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)* force->kspace->dgamma(r/cut_coul); } else { grij = g_ewald * r; expm2 = exp(-grij*grij); derfc = erfc(grij); } if (cut_respa == NULL) { c_tmp = qqrd2e/r; if (msmflag) { f_tmp = qqrd2e/r * fgamma; e_tmp = qqrd2e/r * egamma; } else { f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); e_tmp = qqrd2e/r * derfc; } } else { c_tmp = 0.0; p_tmp = qqrd2e/r; if (msmflag) { f_tmp = qqrd2e/r * (fgamma - 1.0); e_tmp = qqrd2e/r * egamma; v_tmp = qqrd2e/r * fgamma; } else { f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); e_tmp = qqrd2e/r * derfc; v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); } if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); } else { if (msmflag) f_tmp = qqrd2e/r * fgamma; else f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); c_tmp = qqrd2e/r; } } } drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]); dftable[itablemax] = f_tmp - ftable[itablemax]; dctable[itablemax] = c_tmp - ctable[itablemax]; detable[itablemax] = e_tmp - etable[itablemax]; if (cut_respa) { dvtable[itablemax] = v_tmp - vtable[itablemax]; dptable[itablemax] = p_tmp - ptable[itablemax]; } } } /* ---------------------------------------------------------------------- setup force tables for dispersion used in compute routines ------------------------------------------------------------------------- */ void Pair::init_tables_disp(double cut_lj_global) { int masklo,maskhi; double rsq; double g_ewald_6 = force->kspace->g_ewald_6; double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2; tabinnerdispsq = tabinner_disp*tabinner_disp; init_bitmap(tabinner_disp,cut_lj_global,ndisptablebits, masklo,maskhi,ndispmask,ndispshiftbits); int ntable = 1; for (int i = 0; i < ndisptablebits; i++) ntable *= 2; // linear lookup tables of length N = 2^ndisptablebits // stored value = value at lower edge of bin // d values = delta from lower edge to upper edge of bin if (fdisptable) free_disp_tables(); memory->create(rdisptable,ntable,"pair:rdisptable"); memory->create(fdisptable,ntable,"pair:fdisptable"); memory->create(edisptable,ntable,"pair:edisptable"); memory->create(drdisptable,ntable,"pair:drdisptable"); memory->create(dfdisptable,ntable,"pair:dfdisptable"); memory->create(dedisptable,ntable,"pair:dedisptable"); union_int_float_t rsq_lookup; union_int_float_t minrsq_lookup; int itablemin; minrsq_lookup.i = 0 << ndispshiftbits; minrsq_lookup.i |= maskhi; for (int i = 0; i < ntable; i++) { rsq_lookup.i = i << ndispshiftbits; rsq_lookup.i |= masklo; if (rsq_lookup.f < tabinnerdispsq) { rsq_lookup.i = i << ndispshiftbits; rsq_lookup.i |= maskhi; } rsq = rsq_lookup.f; register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2); rdisptable[i] = rsq_lookup.f; fdisptable[i] = g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; edisptable[i] = g6*((a2+1.0)*a2+0.5)*x2; minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); } tabinnerdispsq = minrsq_lookup.f; int ntablem1 = ntable - 1; for (int i = 0; i < ntablem1; i++) { drdisptable[i] = 1.0/(rdisptable[i+1] - rdisptable[i]); dfdisptable[i] = fdisptable[i+1] - fdisptable[i]; dedisptable[i] = edisptable[i+1] - edisptable[i]; } // get the delta values for the last table entries // tables are connected periodically between 0 and ntablem1 drdisptable[ntablem1] = 1.0/(rdisptable[0] - rdisptable[ntablem1]); dfdisptable[ntablem1] = fdisptable[0] - fdisptable[ntablem1]; dedisptable[ntablem1] = edisptable[0] - edisptable[ntablem1]; // get the correct delta values at itablemax // smallest r is in bin itablemin // largest r is in bin itablemax, which is itablemin-1, // or ntablem1 if itablemin=0 // deltas at itablemax only needed if corresponding rsq < cut*cut // if so, compute deltas between rsq and cut*cut double f_tmp,e_tmp; double cut_lj_globalsq; itablemin = minrsq_lookup.i & ndispmask; itablemin >>= ndispshiftbits; int itablemax = itablemin - 1; if (itablemin == 0) itablemax = ntablem1; rsq_lookup.i = itablemax << ndispshiftbits; rsq_lookup.i |= maskhi; if (rsq_lookup.f < (cut_lj_globalsq = cut_lj_global * cut_lj_global)) { rsq_lookup.f = cut_lj_globalsq; register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2); f_tmp = g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; e_tmp = g6*((a2+1.0)*a2+0.5)*x2; drdisptable[itablemax] = 1.0/(rsq_lookup.f - rdisptable[itablemax]); dfdisptable[itablemax] = f_tmp - fdisptable[itablemax]; dedisptable[itablemax] = e_tmp - edisptable[itablemax]; } } /* ---------------------------------------------------------------------- free memory for tables used in Coulombic pair computations ------------------------------------------------------------------------- */ void Pair::free_tables() { memory->destroy(rtable); memory->destroy(drtable); memory->destroy(ftable); memory->destroy(dftable); memory->destroy(ctable); memory->destroy(dctable); memory->destroy(etable); memory->destroy(detable); memory->destroy(vtable); memory->destroy(dvtable); memory->destroy(ptable); memory->destroy(dptable); } /* ---------------------------------------------------------------------- free memory for tables used in pair computations for dispersion ------------------------------------------------------------------------- */ void Pair::free_disp_tables() { memory->destroy(rdisptable); memory->destroy(drdisptable); memory->destroy(fdisptable); memory->destroy(dfdisptable); memory->destroy(edisptable); memory->destroy(dedisptable); } /* ---------------------------------------------------------------------- mixing of pair potential prefactors (epsilon) ------------------------------------------------------------------------- */ double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2) { if (mix_flag == GEOMETRIC) return sqrt(eps1*eps2); else if (mix_flag == ARITHMETIC) return sqrt(eps1*eps2); else if (mix_flag == SIXTHPOWER) return (2.0 * sqrt(eps1*eps2) * pow(sig1,3.0) * pow(sig2,3.0) / (pow(sig1,6.0) + pow(sig2,6.0))); else return 0.0; } /* ---------------------------------------------------------------------- mixing of pair potential distances (sigma, cutoff) ------------------------------------------------------------------------- */ double Pair::mix_distance(double sig1, double sig2) { if (mix_flag == GEOMETRIC) return sqrt(sig1*sig2); else if (mix_flag == ARITHMETIC) return (0.5 * (sig1+sig2)); else if (mix_flag == SIXTHPOWER) return pow((0.5 * (pow(sig1,6.0) + pow(sig2,6.0))),1.0/6.0); else return 0.0; } /* ---------------------------------------------------------------------- */ void Pair::compute_dummy(int eflag, int vflag) { if (eflag || vflag) ev_setup(eflag,vflag); else evflag = 0; } /* ---------------------------------------------------------------------- setup for energy, virial computation see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) ------------------------------------------------------------------------- */ void Pair::ev_setup(int eflag, int vflag) { int i,n; evflag = 1; eflag_either = eflag; eflag_global = eflag % 2; eflag_atom = eflag / 2; vflag_either = vflag; vflag_global = vflag % 4; vflag_atom = vflag / 4; // reallocate per-atom arrays if necessary if (eflag_atom && atom->nmax > maxeatom) { maxeatom = atom->nmax; memory->destroy(eatom); memory->create(eatom,comm->nthreads*maxeatom,"pair:eatom"); } if (vflag_atom && atom->nmax > maxvatom) { maxvatom = atom->nmax; memory->destroy(vatom); memory->create(vatom,comm->nthreads*maxvatom,6,"pair:vatom"); } // zero accumulators // use force->newton instead of newton_pair // b/c some bonds/dihedrals call pair::ev_tally with pairwise info if (eflag_global) eng_vdwl = eng_coul = 0.0; if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; if (eflag_atom) { n = atom->nlocal; if (force->newton) n += atom->nghost; for (i = 0; i < n; i++) eatom[i] = 0.0; } if (vflag_atom) { n = atom->nlocal; if (force->newton) n += atom->nghost; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; vatom[i][1] = 0.0; vatom[i][2] = 0.0; vatom[i][3] = 0.0; vatom[i][4] = 0.0; vatom[i][5] = 0.0; } } // if vflag_global = 2 and pair::compute() calls virial_fdotr_compute() // compute global virial via (F dot r) instead of via pairwise summation // unset other flags as appropriate if (vflag_global == 2 && no_virial_fdotr_compute == 0) { vflag_fdotr = 1; vflag_global = 0; if (vflag_atom == 0) vflag_either = 0; if (vflag_either == 0 && eflag_either == 0) evflag = 0; } else vflag_fdotr = 0; if (lmp->cuda) lmp->cuda->evsetup_eatom_vatom(eflag_atom,vflag_atom); } /* ---------------------------------------------------------------------- set all flags to zero for energy, virial computation called by some complicated many-body potentials that use individual flags to insure no holdover of flags from previous timestep ------------------------------------------------------------------------- */ void Pair::ev_unset() { evflag = 0; eflag_either = 0; eflag_global = 0; eflag_atom = 0; vflag_either = 0; vflag_global = 0; vflag_atom = 0; vflag_fdotr = 0; } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators need i < nlocal test since called by bond_quartic and dihedral_charmm ------------------------------------------------------------------------- */ void Pair::ev_tally(int i, int j, int nlocal, int newton_pair, double evdwl, double ecoul, double fpair, double delx, double dely, double delz) { double evdwlhalf,ecoulhalf,epairhalf,v[6]; if (eflag_either) { if (eflag_global) { if (newton_pair) { eng_vdwl += evdwl; eng_coul += ecoul; } else { evdwlhalf = 0.5*evdwl; ecoulhalf = 0.5*ecoul; if (i < nlocal) { eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } if (j < nlocal) { eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } } } if (eflag_atom) { epairhalf = 0.5 * (evdwl + ecoul); if (newton_pair || i < nlocal) eatom[i] += epairhalf; if (newton_pair || j < nlocal) eatom[j] += epairhalf; } } if (vflag_either) { v[0] = delx*delx*fpair; v[1] = dely*dely*fpair; v[2] = delz*delz*fpair; v[3] = delx*dely*fpair; v[4] = delx*delz*fpair; v[5] = dely*delz*fpair; if (vflag_global) { if (newton_pair) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } else { if (i < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } if (j < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } } } if (vflag_atom) { if (newton_pair || i < nlocal) { vatom[i][0] += 0.5*v[0]; vatom[i][1] += 0.5*v[1]; vatom[i][2] += 0.5*v[2]; vatom[i][3] += 0.5*v[3]; vatom[i][4] += 0.5*v[4]; vatom[i][5] += 0.5*v[5]; } if (newton_pair || j < nlocal) { vatom[j][0] += 0.5*v[0]; vatom[j][1] += 0.5*v[1]; vatom[j][2] += 0.5*v[2]; vatom[j][3] += 0.5*v[3]; vatom[j][4] += 0.5*v[4]; vatom[j][5] += 0.5*v[5]; } } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators can use this version with full neighbor lists ------------------------------------------------------------------------- */ void Pair::ev_tally_full(int i, double evdwl, double ecoul, double fpair, double delx, double dely, double delz) { double v[6]; if (eflag_either) { if (eflag_global) { eng_vdwl += 0.5*evdwl; eng_coul += 0.5*ecoul; } if (eflag_atom) eatom[i] += 0.5 * (evdwl + ecoul); } if (vflag_either) { v[0] = 0.5*delx*delx*fpair; v[1] = 0.5*dely*dely*fpair; v[2] = 0.5*delz*delz*fpair; v[3] = 0.5*delx*dely*fpair; v[4] = 0.5*delx*delz*fpair; v[5] = 0.5*dely*delz*fpair; if (vflag_global) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } if (vflag_atom) { vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz ------------------------------------------------------------------------- */ void Pair::ev_tally_xyz(int i, int j, int nlocal, int newton_pair, double evdwl, double ecoul, double fx, double fy, double fz, double delx, double dely, double delz) { double evdwlhalf,ecoulhalf,epairhalf,v[6]; if (eflag_either) { if (eflag_global) { if (newton_pair) { eng_vdwl += evdwl; eng_coul += ecoul; } else { evdwlhalf = 0.5*evdwl; ecoulhalf = 0.5*ecoul; if (i < nlocal) { eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } if (j < nlocal) { eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } } } if (eflag_atom) { epairhalf = 0.5 * (evdwl + ecoul); if (newton_pair || i < nlocal) eatom[i] += epairhalf; if (newton_pair || j < nlocal) eatom[j] += epairhalf; } } if (vflag_either) { v[0] = delx*fx; v[1] = dely*fy; v[2] = delz*fz; v[3] = delx*fy; v[4] = delx*fz; v[5] = dely*fz; if (vflag_global) { if (newton_pair) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } else { if (i < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } if (j < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } } } if (vflag_atom) { if (newton_pair || i < nlocal) { vatom[i][0] += 0.5*v[0]; vatom[i][1] += 0.5*v[1]; vatom[i][2] += 0.5*v[2]; vatom[i][3] += 0.5*v[3]; vatom[i][4] += 0.5*v[4]; vatom[i][5] += 0.5*v[5]; } if (newton_pair || j < nlocal) { vatom[j][0] += 0.5*v[0]; vatom[j][1] += 0.5*v[1]; vatom[j][2] += 0.5*v[2]; vatom[j][3] += 0.5*v[3]; vatom[j][4] += 0.5*v[4]; vatom[j][5] += 0.5*v[5]; } } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz called when using full neighbor lists ------------------------------------------------------------------------- */ void Pair::ev_tally_xyz_full(int i, double evdwl, double ecoul, double fx, double fy, double fz, double delx, double dely, double delz) { double evdwlhalf,ecoulhalf,epairhalf,v[6]; if (eflag_either) { if (eflag_global) { evdwlhalf = 0.5*evdwl; ecoulhalf = 0.5*ecoul; eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } if (eflag_atom) { epairhalf = 0.5 * (evdwl + ecoul); eatom[i] += epairhalf; } } if (vflag_either) { v[0] = 0.5*delx*fx; v[1] = 0.5*dely*fy; v[2] = 0.5*delz*fz; v[3] = 0.5*delx*fy; v[4] = 0.5*delx*fz; v[5] = 0.5*dely*fz; if (vflag_global) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } if (vflag_atom) { vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators called by SW and hbond potentials, newton_pair is always on virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk ------------------------------------------------------------------------- */ void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul, double *fj, double *fk, double *drji, double *drki) { double epairthird,v[6]; if (eflag_either) { if (eflag_global) { eng_vdwl += evdwl; eng_coul += ecoul; } if (eflag_atom) { epairthird = THIRD * (evdwl + ecoul); eatom[i] += epairthird; eatom[j] += epairthird; eatom[k] += epairthird; } } if (vflag_either) { v[0] = drji[0]*fj[0] + drki[0]*fk[0]; v[1] = drji[1]*fj[1] + drki[1]*fk[1]; v[2] = drji[2]*fj[2] + drki[2]*fk[2]; v[3] = drji[0]*fj[1] + drki[0]*fk[1]; v[4] = drji[0]*fj[2] + drki[0]*fk[2]; v[5] = drji[1]*fj[2] + drki[1]*fk[2]; if (vflag_global) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } if (vflag_atom) { vatom[i][0] += THIRD*v[0]; vatom[i][1] += THIRD*v[1]; vatom[i][2] += THIRD*v[2]; vatom[i][3] += THIRD*v[3]; vatom[i][4] += THIRD*v[4]; vatom[i][5] += THIRD*v[5]; vatom[j][0] += THIRD*v[0]; vatom[j][1] += THIRD*v[1]; vatom[j][2] += THIRD*v[2]; vatom[j][3] += THIRD*v[3]; vatom[j][4] += THIRD*v[4]; vatom[j][5] += THIRD*v[5]; vatom[k][0] += THIRD*v[0]; vatom[k][1] += THIRD*v[1]; vatom[k][2] += THIRD*v[2]; vatom[k][3] += THIRD*v[3]; vatom[k][4] += THIRD*v[4]; vatom[k][5] += THIRD*v[5]; } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators called by AIREBO potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::ev_tally4(int i, int j, int k, int m, double evdwl, double *fi, double *fj, double *fk, double *drim, double *drjm, double *drkm) { double epairfourth,v[6]; if (eflag_either) { if (eflag_global) eng_vdwl += evdwl; if (eflag_atom) { epairfourth = 0.25 * evdwl; eatom[i] += epairfourth; eatom[j] += epairfourth; eatom[k] += epairfourth; eatom[m] += epairfourth; } } if (vflag_atom) { v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; } } /* ---------------------------------------------------------------------- tally ecoul and virial into each of atoms in list called by TIP4P potential, newton_pair is always on weight assignments by alpha, so contribution is all to O atom as alpha -> 0.0 key = 0 if neither atom = water O key = 1 if first atom = water O key = 2 if second atom = water O key = 3 if both atoms = water O ------------------------------------------------------------------------- */ void Pair::ev_tally_tip4p(int key, int *list, double *v, double ecoul, double alpha) { int i; if (eflag_either) { if (eflag_global) eng_coul += ecoul; if (eflag_atom) { if (key == 0) { eatom[list[0]] += 0.5*ecoul; eatom[list[1]] += 0.5*ecoul; } else if (key == 1) { eatom[list[0]] += 0.5*ecoul*(1-alpha); eatom[list[1]] += 0.25*ecoul*alpha; eatom[list[2]] += 0.25*ecoul*alpha; eatom[list[3]] += 0.5*ecoul; } else if (key == 2) { eatom[list[0]] += 0.5*ecoul; eatom[list[1]] += 0.5*ecoul*(1-alpha); eatom[list[2]] += 0.25*ecoul*alpha; eatom[list[3]] += 0.25*ecoul*alpha; } else { eatom[list[0]] += 0.5*ecoul*(1-alpha); eatom[list[1]] += 0.25*ecoul*alpha; eatom[list[2]] += 0.25*ecoul*alpha; eatom[list[3]] += 0.5*ecoul*(1-alpha); eatom[list[4]] += 0.25*ecoul*alpha; eatom[list[5]] += 0.25*ecoul*alpha; } } } if (vflag_either) { if (vflag_global) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } if (vflag_atom) { if (key == 0) { for (i = 0; i <= 5; i++) { vatom[list[0]][i] += 0.5*v[i]; vatom[list[1]][i] += 0.5*v[i]; } } else if (key == 1) { for (i = 0; i <= 5; i++) { vatom[list[0]][i] += 0.5*v[i]*(1-alpha); vatom[list[1]][i] += 0.25*v[i]*alpha; vatom[list[2]][i] += 0.25*v[i]*alpha; vatom[list[3]][i] += 0.5*v[i]; } } else if (key == 2) { for (i = 0; i <= 5; i++) { vatom[list[0]][i] += 0.5*v[i]; vatom[list[1]][i] += 0.5*v[i]*(1-alpha); vatom[list[2]][i] += 0.25*v[i]*alpha; vatom[list[3]][i] += 0.25*v[i]*alpha; } } else { for (i = 0; i <= 5; i++) { vatom[list[0]][i] += 0.5*v[i]*(1-alpha); vatom[list[1]][i] += 0.25*v[i]*alpha; vatom[list[2]][i] += 0.25*v[i]*alpha; vatom[list[3]][i] += 0.5*v[i]*(1-alpha); vatom[list[4]][i] += 0.25*v[i]*alpha; vatom[list[5]][i] += 0.25*v[i]*alpha; } } } } } /* ---------------------------------------------------------------------- tally virial into per-atom accumulators called by REAX/C potential, newton_pair is always on fi is magnitude of force on atom i ------------------------------------------------------------------------- */ void Pair::v_tally(int i, double *fi, double *deli) { double v[6]; v[0] = 0.5*deli[0]*fi[0]; v[1] = 0.5*deli[1]*fi[1]; v[2] = 0.5*deli[2]*fi[2]; v[3] = 0.5*deli[0]*fi[1]; v[4] = 0.5*deli[0]*fi[2]; v[5] = 0.5*deli[1]*fi[2]; vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; } /* ---------------------------------------------------------------------- tally virial into per-atom accumulators called by AIREBO potential, newton_pair is always on fpair is magnitude of force on atom I ------------------------------------------------------------------------- */ void Pair::v_tally2(int i, int j, double fpair, double *drij) { double v[6]; v[0] = 0.5 * drij[0]*drij[0]*fpair; v[1] = 0.5 * drij[1]*drij[1]*fpair; v[2] = 0.5 * drij[2]*drij[2]*fpair; v[3] = 0.5 * drij[0]*drij[1]*fpair; v[4] = 0.5 * drij[0]*drij[2]*fpair; v[5] = 0.5 * drij[1]*drij[2]*fpair; vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; } /* ---------------------------------------------------------------------- tally virial into per-atom accumulators called by AIREBO and Tersoff potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally3(int i, int j, int k, double *fi, double *fj, double *drik, double *drjk) { double v[6]; v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]); v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]); v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]); v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]); v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; } /* ---------------------------------------------------------------------- tally virial into per-atom accumulators called by AIREBO potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally4(int i, int j, int k, int m, double *fi, double *fj, double *fk, double *drim, double *drjm, double *drkm) { double v[6]; v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; } /* ---------------------------------------------------------------------- tally virial into global and per-atom accumulators called by pair lubricate potential with 6 tensor components ------------------------------------------------------------------------- */ void Pair::v_tally_tensor(int i, int j, int nlocal, int newton_pair, double vxx, double vyy, double vzz, double vxy, double vxz, double vyz) { double v[6]; v[0] = vxx; v[1] = vyy; v[2] = vzz; v[3] = vxy; v[4] = vxz; v[5] = vyz; if (vflag_global) { if (newton_pair) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } else { if (i < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } if (j < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } } } if (vflag_atom) { if (newton_pair || i < nlocal) { vatom[i][0] += 0.5*v[0]; vatom[i][1] += 0.5*v[1]; vatom[i][2] += 0.5*v[2]; vatom[i][3] += 0.5*v[3]; vatom[i][4] += 0.5*v[4]; vatom[i][5] += 0.5*v[5]; } if (newton_pair || j < nlocal) { vatom[j][0] += 0.5*v[0]; vatom[j][1] += 0.5*v[1]; vatom[j][2] += 0.5*v[2]; vatom[j][3] += 0.5*v[3]; vatom[j][4] += 0.5*v[4]; vatom[j][5] += 0.5*v[5]; } } } /* ---------------------------------------------------------------------- compute global pair virial via summing F dot r over own & ghost atoms at this point, only pairwise forces have been accumulated in atom->f ------------------------------------------------------------------------- */ void Pair::virial_fdotr_compute() { double **x = atom->x; double **f = atom->f; // sum over force on all particles including ghosts if (neighbor->includegroup == 0) { int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) { virial[0] += f[i][0]*x[i][0]; virial[1] += f[i][1]*x[i][1]; virial[2] += f[i][2]*x[i][2]; virial[3] += f[i][1]*x[i][0]; virial[4] += f[i][2]*x[i][0]; virial[5] += f[i][2]*x[i][1]; } // neighbor includegroup flag is set // sum over force on initial nfirst particles and ghosts } else { int nall = atom->nfirst; for (int i = 0; i < nall; i++) { virial[0] += f[i][0]*x[i][0]; virial[1] += f[i][1]*x[i][1]; virial[2] += f[i][2]*x[i][2]; virial[3] += f[i][1]*x[i][0]; virial[4] += f[i][2]*x[i][0]; virial[5] += f[i][2]*x[i][1]; } nall = atom->nlocal + atom->nghost; for (int i = atom->nlocal; i < nall; i++) { virial[0] += f[i][0]*x[i][0]; virial[1] += f[i][1]*x[i][1]; virial[2] += f[i][2]*x[i][2]; virial[3] += f[i][1]*x[i][0]; virial[4] += f[i][2]*x[i][0]; virial[5] += f[i][2]*x[i][1]; } } // prevent multiple calls to update the virial // when a hybrid pair style uses both a gpu and non-gpu pair style // or when respa is used with gpu pair styles vflag_fdotr = 0; } /* ---------------------------------------------------------------------- write a table of pair potential energy/force vs distance to a file ------------------------------------------------------------------------- */ void Pair::write_file(int narg, char **arg) { if (narg < 8) error->all(FLERR,"Illegal pair_write command"); if (single_enable == 0) error->all(FLERR,"Pair style does not support pair_write"); // parse arguments int itype = force->inumeric(FLERR,arg[0]); int jtype = force->inumeric(FLERR,arg[1]); if (itype < 1 || itype > atom->ntypes || jtype < 1 || jtype > atom->ntypes) error->all(FLERR,"Invalid atom types in pair_write command"); int n = force->inumeric(FLERR,arg[2]); int style = NONE; if (strcmp(arg[3],"r") == 0) style = RLINEAR; else if (strcmp(arg[3],"rsq") == 0) style = RSQ; else if (strcmp(arg[3],"bitmap") == 0) style = BMP; else error->all(FLERR,"Invalid style in pair_write command"); double inner = force->numeric(FLERR,arg[4]); double outer = force->numeric(FLERR,arg[5]); if (inner <= 0.0 || inner >= outer) error->all(FLERR,"Invalid cutoffs in pair_write command"); // open file in append mode // print header in format used by pair_style table int me; MPI_Comm_rank(world,&me); FILE *fp; if (me == 0) { fp = fopen(arg[6],"a"); if (fp == NULL) error->one(FLERR,"Cannot open pair_write file"); fprintf(fp,"# Pair potential %s for atom types %d %d: i,r,energy,force\n", force->pair_style,itype,jtype); if (style == RLINEAR) fprintf(fp,"\n%s\nN %d R %g %g\n\n",arg[7],n,inner,outer); if (style == RSQ) fprintf(fp,"\n%s\nN %d RSQ %g %g\n\n",arg[7],n,inner,outer); } // initialize potentials before evaluating pair potential // insures all pair coeffs are set and force constants // also initialize neighbor so that neighbor requests are processed // NOTE: might be safest to just do lmp->init() force->init(); neighbor->init(); // if pair style = any of EAM, swap in dummy fp vector double eamfp[2]; eamfp[0] = eamfp[1] = 0.0; double *eamfp_hold; Pair *epair = force->pair_match("eam",0); if (epair) epair->swap_eam(eamfp,&eamfp_hold); // if atom style defines charge, swap in dummy q vec double q[2]; q[0] = q[1] = 1.0; if (narg == 10) { q[0] = force->numeric(FLERR,arg[8]); q[1] = force->numeric(FLERR,arg[9]); } double *q_hold; if (atom->q) { q_hold = atom->q; atom->q = q; } // evaluate energy and force at each of N distances int masklo,maskhi,nmask,nshiftbits; if (style == BMP) { init_bitmap(inner,outer,n,masklo,maskhi,nmask,nshiftbits); int ntable = 1 << n; if (me == 0) fprintf(fp,"\n%s\nN %d BITMAP %g %g\n\n",arg[7],ntable,inner,outer); n = ntable; } double r,e,f,rsq; union_int_float_t rsq_lookup; for (int i = 0; i < n; i++) { if (style == RLINEAR) { r = inner + (outer-inner) * i/(n-1); rsq = r*r; } else if (style == RSQ) { rsq = inner*inner + (outer*outer - inner*inner) * i/(n-1); r = sqrt(rsq); } else if (style == BMP) { rsq_lookup.i = i << nshiftbits; rsq_lookup.i |= masklo; if (rsq_lookup.f < inner*inner) { rsq_lookup.i = i << nshiftbits; rsq_lookup.i |= maskhi; } rsq = rsq_lookup.f; r = sqrt(rsq); } if (rsq < cutsq[itype][jtype]) { e = single(0,1,itype,jtype,rsq,1.0,1.0,f); f *= r; } else e = f = 0.0; if (me == 0) fprintf(fp,"%d %g %g %g\n",i+1,r,e,f); } // restore original vecs that were swapped in for double *tmp; if (epair) epair->swap_eam(eamfp_hold,&tmp); if (atom->q) atom->q = q_hold; if (me == 0) fclose(fp); } /* ---------------------------------------------------------------------- define bitmap parameters based on inner and outer cutoffs ------------------------------------------------------------------------- */ void Pair::init_bitmap(double inner, double outer, int ntablebits, int &masklo, int &maskhi, int &nmask, int &nshiftbits) { if (sizeof(int) != sizeof(float)) error->all(FLERR,"Bitmapped lookup tables require int/float be same size"); if (ntablebits > sizeof(float)*CHAR_BIT) error->all(FLERR,"Too many total bits for bitmapped lookup table"); if (inner >= outer) error->warning(FLERR,"Table inner cutoff >= outer cutoff"); int nlowermin = 1; while (!((pow(double(2),(double)nlowermin) <= inner*inner) && (pow(double(2),(double)nlowermin+1.0) > inner*inner))) { if (pow(double(2),(double)nlowermin) <= inner*inner) nlowermin++; else nlowermin--; } int nexpbits = 0; double required_range = outer*outer / pow(double(2),(double)nlowermin); double available_range = 2.0; while (available_range < required_range) { nexpbits++; available_range = pow(double(2),pow(double(2),(double)nexpbits)); } int nmantbits = ntablebits - nexpbits; if (nexpbits > sizeof(float)*CHAR_BIT - FLT_MANT_DIG) error->all(FLERR,"Too many exponent bits for lookup table"); if (nmantbits+1 > FLT_MANT_DIG) error->all(FLERR,"Too many mantissa bits for lookup table"); if (nmantbits < 3) error->all(FLERR,"Too few bits for lookup table"); nshiftbits = FLT_MANT_DIG - (nmantbits+1); nmask = 1; for (int j = 0; j < ntablebits+nshiftbits; j++) nmask *= 2; nmask -= 1; union_int_float_t rsq_lookup; rsq_lookup.f = outer*outer; maskhi = rsq_lookup.i & ~(nmask); rsq_lookup.f = inner*inner; masklo = rsq_lookup.i & ~(nmask); } /* ---------------------------------------------------------------------- */ double Pair::memory_usage() { double bytes = comm->nthreads*maxeatom * sizeof(double); bytes += comm->nthreads*maxvatom*6 * sizeof(double); return bytes; } diff --git a/src/pair.h b/src/pair.h index d58b79de9..e7278a0d7 100644 --- a/src/pair.h +++ b/src/pair.h @@ -1,321 +1,327 @@ /* -*- 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" +#include "accelerator_kokkos.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: static int instance_total; // # of Pair classes ever instantiated 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 manybody_flag; // 1 if a manybody potential 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 dipoleflag; // 1 if compatible with dipole solver int reinitflag; // 1 if compatible with fix adapt and alike 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() // KOKKOS host/device flag and data masks ExecutionSpace execution_space; unsigned int datamask_read,datamask_modify; Pair(class LAMMPS *); virtual ~Pair(); // top-level Pair methods void init(); virtual 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_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} + virtual int pack_forward_comm_kokkos(int, DAT::tdual_int_2d, int, DAT::tdual_xfloat_1d&, int, int *) {return 0;}; + virtual void unpack_forward_comm_kokkos(int, int, DAT::tdual_xfloat_1d&) {} 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: int instance_me; // which Pair class instantiation I am enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // mixing options int special_lj[4]; // copied from force->special_lj for Kokkos 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; + int copymode; // if set, do not deallocate during destruction + // required when classes are used as functors by Kokkos + 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. W: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions This is likely not what you want to do. The exclusion settings will eliminate neighbors in the neighbor list, which the manybody potential needs to calculated its terms correctly. 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 requres a KSpace style Self-explanatory. 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_buck.cpp b/src/pair_buck.cpp index 096d4d0ff..5eaadd72c 100644 --- a/src/pair_buck.cpp +++ b/src/pair_buck.cpp @@ -1,406 +1,406 @@ /* ---------------------------------------------------------------------- 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 "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_buck.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neigh_list.h" #include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; /* ---------------------------------------------------------------------- */ PairBuck::PairBuck(LAMMPS *lmp) : Pair(lmp) { writedata = 1; } /* ---------------------------------------------------------------------- */ PairBuck::~PairBuck() { - if (allocated) { + if (allocated && !copymode) { memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(cut); memory->destroy(a); memory->destroy(rho); memory->destroy(c); memory->destroy(rhoinv); memory->destroy(buck1); memory->destroy(buck2); memory->destroy(offset); } } /* ---------------------------------------------------------------------- */ void PairBuck::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,forcebuck,factor_lj; double r,rexp; 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; r = sqrt(rsq); rexp = exp(-r*rhoinv[itype][jtype]); forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; fpair = factor_lj*forcebuck*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 = a[itype][jtype]*rexp - c[itype][jtype]*r6inv - 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(); } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairBuck::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_lj"); memory->create(a,n+1,n+1,"pair:a"); memory->create(rho,n+1,n+1,"pair:rho"); memory->create(c,n+1,n+1,"pair:c"); memory->create(rhoinv,n+1,n+1,"pair:rhoinv"); memory->create(buck1,n+1,n+1,"pair:buck1"); memory->create(buck2,n+1,n+1,"pair:buck2"); memory->create(offset,n+1,n+1,"pair:offset"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairBuck::settings(int narg, char **arg) { if (narg != 1) error->all(FLERR,"Illegal pair_style command"); cut_global = force->numeric(FLERR,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 PairBuck::coeff(int narg, char **arg) { if (narg < 5 || narg > 6) 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 a_one = force->numeric(FLERR,arg[2]); double rho_one = force->numeric(FLERR,arg[3]); if (rho_one <= 0) error->all(FLERR,"Incorrect args for pair coefficients"); double c_one = force->numeric(FLERR,arg[4]); double cut_one = cut_global; if (narg == 6) cut_one = force->numeric(FLERR,arg[5]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { a[i][j] = a_one; rho[i][j] = rho_one; c[i][j] = c_one; cut[i][j] = cut_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairBuck::init_one(int i, int j) { if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); rhoinv[i][j] = 1.0/rho[i][j]; buck1[i][j] = a[i][j]/rho[i][j]; buck2[i][j] = 6.0*c[i][j]; if (offset_flag) { double rexp = exp(-cut[i][j]/rho[i][j]); offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut[i][j],6.0); } else offset[i][j] = 0.0; a[j][i] = a[i][j]; c[j][i] = c[i][j]; rhoinv[j][i] = rhoinv[i][j]; buck1[j][i] = buck1[i][j]; buck2[j][i] = buck2[i][j]; offset[j][i] = offset[i][j]; // 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 rho1 = rho[i][j]; double rho2 = rho1*rho1; double rho3 = rho2*rho1; double rc = cut[i][j]; double rc2 = rc*rc; double rc3 = rc2*rc; etail_ij = 2.0*MY_PI*all[0]*all[1]* (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) - c[i][j]/(3.0*rc3)); ptail_ij = (-1/3.0)*2.0*MY_PI*all[0]*all[1]* (-a[i][j]*exp(-rc/rho1)* (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3); } return cut[i][j]; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairBuck::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(&a[i][j],sizeof(double),1,fp); fwrite(&rho[i][j],sizeof(double),1,fp); fwrite(&c[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairBuck::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(&a[i][j],sizeof(double),1,fp); fread(&rho[i][j],sizeof(double),1,fp); fread(&c[i][j],sizeof(double),1,fp); fread(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&a[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&rho[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&c[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairBuck::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 PairBuck::read_restart_settings(FILE *fp) { if (comm->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 PairBuck::write_data(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) fprintf(fp,"%d %g %g %g\n",i,a[i][i],rho[i][i],c[i][i]); } /* ---------------------------------------------------------------------- proc 0 writes all pairs to data file ------------------------------------------------------------------------- */ void PairBuck::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 %g\n",i,j, a[i][j],rho[i][j],c[i][j],cut[i][j]); } /* ---------------------------------------------------------------------- */ double PairBuck::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double r2inv,r6inv,r,rexp,forcebuck,phibuck; r2inv = 1.0/rsq; r6inv = r2inv*r2inv*r2inv; r = sqrt(rsq); rexp = exp(-r*rhoinv[itype][jtype]); forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; fforce = factor_lj*forcebuck*r2inv; phibuck = a[itype][jtype]*rexp - c[itype][jtype]*r6inv - offset[itype][jtype]; return factor_lj*phibuck; } /* ---------------------------------------------------------------------- */ void *PairBuck::extract(const char *str, int &dim) { dim = 2; if (strcmp(str,"a") == 0) return (void *) a; if (strcmp(str,"c") == 0) return (void *) c; return NULL; } diff --git a/src/pair_buck.h b/src/pair_buck.h index bdcb1560e..d74934de8 100644 --- a/src/pair_buck.h +++ b/src/pair_buck.h @@ -1,75 +1,75 @@ /* -*- 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. ------------------------------------------------------------------------- */ #ifdef PAIR_CLASS PairStyle(buck,PairBuck) #else #ifndef LMP_PAIR_BUCK_H #define LMP_PAIR_BUCK_H #include "pair.h" namespace LAMMPS_NS { class PairBuck : public Pair { public: PairBuck(class LAMMPS *); virtual ~PairBuck(); virtual void compute(int, int); void settings(int, char **); void coeff(int, char **); 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 &); protected: double cut_global; double **cut; double **a,**rho,**c; double **rhoinv,**buck1,**buck2,**offset; - void allocate(); + virtual 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: 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. */ diff --git a/src/pair_coul_dsf.cpp b/src/pair_coul_dsf.cpp index 0aba5548d..17b3dfaa7 100644 --- a/src/pair_coul_dsf.cpp +++ b/src/pair_coul_dsf.cpp @@ -1,331 +1,331 @@ /* ---------------------------------------------------------------------- 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: Trung Dac Nguyen (ORNL) References: Fennell and Gezelter, JCP 124, 234104 (2006) ------------------------------------------------------------------------- */ #include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_coul_dsf.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" #include "memory.h" #include "math_const.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; #define EWALD_F 1.12837917 #define EWALD_P 0.3275911 #define A1 0.254829592 #define A2 -0.284496736 #define A3 1.421413741 #define A4 -1.453152027 #define A5 1.061405429 /* ---------------------------------------------------------------------- */ PairCoulDSF::PairCoulDSF(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; } /* ---------------------------------------------------------------------- */ PairCoulDSF::~PairCoulDSF() { - if (allocated) { + if (allocated && !copymode) { memory->destroy(setflag); memory->destroy(cutsq); } } /* ---------------------------------------------------------------------- */ void PairCoulDSF::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; double r,rsq,r2inv,forcecoul,factor_coul; double prefactor,erfcc,erfcd,e_self,t; int *ilist,*jlist,*numneigh,**firstneigh; 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 nlocal = atom->nlocal; double *special_coul = force->special_coul; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; 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]; jlist = firstneigh[i]; jnum = numneigh[i]; if (eflag) { e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e; 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; if (rsq < cut_coulsq) { r2inv = 1.0/rsq; r = sqrt(rsq); prefactor = factor_coul * qqrd2e*qtmp*q[j]/r; erfcd = exp(-alpha*alpha*rsq); t = 1.0 / (1.0 + EWALD_P*alpha*r); erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + r*f_shift) * r; fpair = forcecoul * 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) { ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); } 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 PairCoulDSF::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 ------------------------------------------------------------------------- */ void PairCoulDSF::settings(int narg, char **arg) { if (narg != 2) error->all(FLERR,"Illegal pair_style command"); alpha = force->numeric(FLERR,arg[0]); cut_coul = force->numeric(FLERR,arg[1]); } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairCoulDSF::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 PairCoulDSF::init_style() { if (!atom->q_flag) error->all(FLERR,"Pair style coul/dsf requires atom attribute q"); neighbor->request(this,instance_me); cut_coulsq = cut_coul * cut_coul; double erfcc = erfc(alpha*cut_coul); double erfcd = exp(-alpha*alpha*cut_coul*cut_coul); f_shift = -(erfcc/cut_coulsq + 2.0/MY_PIS*alpha*erfcd/cut_coul); e_shift = erfcc/cut_coul - f_shift*cut_coul; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairCoulDSF::init_one(int i, int j) { return cut_coul; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairCoulDSF::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 PairCoulDSF::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 PairCoulDSF::write_restart_settings(FILE *fp) { fwrite(&alpha,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 PairCoulDSF::read_restart_settings(FILE *fp) { if (comm->me == 0) { fread(&alpha,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(&alpha,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); } /* ---------------------------------------------------------------------- */ double PairCoulDSF::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double r2inv,r,erfcc,erfcd,prefactor,t; double forcecoul,phicoul; r2inv = 1.0/rsq; double eng = 0.0; if (rsq < cut_coulsq) { r = sqrt(rsq); prefactor = factor_coul * force->qqrd2e * atom->q[i]*atom->q[j]/r; erfcd = exp(-alpha*alpha*rsq); t = 1.0 / (1.0 + EWALD_P*alpha*r); erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS*erfcd + r*f_shift) * r; phicoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); eng += phicoul; } else forcecoul = 0.0; fforce = forcecoul * r2inv; return eng; } /* ---------------------------------------------------------------------- */ void *PairCoulDSF::extract(const char *str, int &dim) { if (strcmp(str,"cut_coul") == 0) { dim = 0; return (void *) &cut_coul; } return NULL; } diff --git a/src/pair_coul_wolf.cpp b/src/pair_coul_wolf.cpp index 799b672a4..4687a75fc 100644 --- a/src/pair_coul_wolf.cpp +++ b/src/pair_coul_wolf.cpp @@ -1,323 +1,323 @@ /* ---------------------------------------------------------------------- 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) { single_enable = 0; } /* ---------------------------------------------------------------------- */ PairCoulWolf::~PairCoulWolf() { - if (allocated) { + if (allocated && !copymode) { memory->destroy(setflag); memory->destroy(cutsq); } } /* ---------------------------------------------------------------------- */ void PairCoulWolf::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; double rsq,forcecoul,factor_coul; double prefactor; double r; 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 nlocal = atom->nlocal; 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]; 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; 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(FLERR,arg[0]); cut_coul = force->numeric(FLERR,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"); neighbor->request(this,instance_me); 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 r,prefactor; double forcecoul,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; }