diff --git a/doc/Manual.html b/doc/Manual.html index cc9f910f05..6c68be1bc5 100644 --- a/doc/Manual.html +++ b/doc/Manual.html @@ -1,7 +1,7 @@ LAMMPS-ICMS Users Manual - + @@ -22,7 +22,7 @@

LAMMPS-ICMS Documentation

-

8 May 2014 version +

9 May 2014 version

Version info:

diff --git a/doc/Manual.txt b/doc/Manual.txt index 3261a9447c..c8c17b355d 100644 --- a/doc/Manual.txt +++ b/doc/Manual.txt @@ -1,6 +1,6 @@ LAMMPS-ICMS Users Manual - + @@ -18,7 +18,7 @@

LAMMPS-ICMS Documentation :c,h3 -8 May 2014 version :c,h4 +9 May 2014 version :c,h4 Version info: :h4 diff --git a/examples/VISCOSITY/in.gk.2d b/examples/VISCOSITY/in.gk.2d index 1438be5341..09fdef0a42 100644 --- a/examples/VISCOSITY/in.gk.2d +++ b/examples/VISCOSITY/in.gk.2d @@ -49,17 +49,28 @@ unfix 2 reset_timestep 0 +# Define distinct components of symmetric traceless stress tensor + variable pxy equal pxy +variable pxx equal pxx-press fix SS all ave/correlate $s $p $d & - v_pxy type auto file profile.gk.2d ave running + v_pxy v_pxx type auto file profile.gk.2d ave running + +# Diagonal components of SS are larger by factor 2-2/d, +# which is 4/3 for d=3, but 1 for d=2. +# See Daivis and Evans, J.Chem.Phys, 100, 541-547 (1994) variable scale equal 1.0/$t*vol*$s*dt -variable v11 equal trap(f_SS[3])*${scale} +variable diagfac equal 2-2/2 +variable vxy equal trap(f_SS[3])*${scale} +variable vxx equal trap(f_SS[4])*${scale}/${diagfac} -thermo_style custom step temp press pxy v_v11 +thermo_style custom step temp press pxy v_vxy v_vxx run 500000 -variable eta equal v_v11 +variable etaxy equal v_vxy +variable etaxx equal v_vxx +variable eta equal 0.5*(${etaxy}+${etaxx}) print "running average viscosity: ${eta}" diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp index 4c08097fa2..0f920a17ba 100644 --- a/src/ASPHERE/compute_temp_asphere.cpp +++ b/src/ASPHERE/compute_temp_asphere.cpp @@ -25,7 +25,6 @@ #include "force.h" #include "domain.h" #include "modify.h" -#include "fix.h" #include "group.h" #include "memory.h" #include "error.h" @@ -128,9 +127,7 @@ void ComputeTempAsphere::init() void ComputeTempAsphere::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -138,6 +135,8 @@ void ComputeTempAsphere::setup() void ComputeTempAsphere::dof_compute() { + if (fix_dof) adjust_dof_fix(); + // 6 dof for 3d, 3 dof for 2d // which dof are included also depends on mode // assume full rotation of extended particles diff --git a/src/ASPHERE/compute_temp_asphere.h b/src/ASPHERE/compute_temp_asphere.h index b746743823..d7335da4bd 100644 --- a/src/ASPHERE/compute_temp_asphere.h +++ b/src/ASPHERE/compute_temp_asphere.h @@ -39,7 +39,7 @@ class ComputeTempAsphere : public Compute { void restore_bias_thr(int, double *, double *); private: - int fix_dof,mode; + int mode; double tfactor; char *id_bias; class Compute *tbias; // ptr to additional bias compute diff --git a/src/DIPOLE/atom_vec_dipole.cpp b/src/DIPOLE/atom_vec_dipole.cpp index 08067d567d..07c8d203d8 100644 --- a/src/DIPOLE/atom_vec_dipole.cpp +++ b/src/DIPOLE/atom_vec_dipole.cpp @@ -24,8 +24,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecDipole::AtomVecDipole(LAMMPS *lmp) : AtomVec(lmp) @@ -48,13 +46,13 @@ AtomVecDipole::AtomVecDipole(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecDipole::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/MISC/fix_efield.cpp b/src/MISC/fix_efield.cpp index bb1aa73ab8..c472bf9f16 100644 --- a/src/MISC/fix_efield.cpp +++ b/src/MISC/fix_efield.cpp @@ -114,8 +114,8 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : force_flag = 0; fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0; - maxatom = -1; - efield = NULL; + maxatom = atom->nmax; + memory->create(efield,maxatom,4,"efield:efield"); } /* ---------------------------------------------------------------------- */ @@ -332,24 +332,16 @@ void FixEfield::post_force(int vflag) modify->clearstep_compute(); if (xstyle == EQUAL) ex = qe2f * input->variable->compute_equal(xvar); - else if (xstyle == ATOM) { - if (efield) input->variable->compute_atom(xvar,igroup,&efield[0][0],3,0); - else input->variable->compute_atom(xvar,igroup,NULL,3,0); - } + else if (xstyle == ATOM) + input->variable->compute_atom(xvar,igroup,&efield[0][0],3,0); if (ystyle == EQUAL) ey = qe2f * input->variable->compute_equal(yvar); - else if (ystyle == ATOM) { - if (efield) input->variable->compute_atom(yvar,igroup,&efield[0][1],3,0); - else input->variable->compute_atom(yvar,igroup,NULL,3,0); - } + else if (ystyle == ATOM) + input->variable->compute_atom(yvar,igroup,&efield[0][1],3,0); if (zstyle == EQUAL) ez = qe2f * input->variable->compute_equal(zvar); - else if (zstyle == ATOM) { - if (efield) input->variable->compute_atom(zvar,igroup,&efield[0][2],3,0); - else input->variable->compute_atom(zvar,igroup,NULL,3,0); - } - if (estyle == ATOM) { - if (efield) input->variable->compute_atom(evar,igroup,&efield[0][3],4,0); - else input->variable->compute_atom(evar,igroup,NULL,4,0); - } + else if (zstyle == ATOM) + input->variable->compute_atom(zvar,igroup,&efield[0][2],3,0); + if (estyle == ATOM) + input->variable->compute_atom(evar,igroup,&efield[0][3],4,0); modify->addstep_compute(update->ntimestep + 1); diff --git a/src/MOLECULE/atom_vec_angle.cpp b/src/MOLECULE/atom_vec_angle.cpp index bafac2391f..8241c693aa 100644 --- a/src/MOLECULE/atom_vec_angle.cpp +++ b/src/MOLECULE/atom_vec_angle.cpp @@ -23,8 +23,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecAngle::AtomVecAngle(LAMMPS *lmp) : AtomVec(lmp) @@ -47,13 +45,13 @@ AtomVecAngle::AtomVecAngle(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecAngle::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/MOLECULE/atom_vec_bond.cpp b/src/MOLECULE/atom_vec_bond.cpp index 030a86ac86..f53c71736e 100644 --- a/src/MOLECULE/atom_vec_bond.cpp +++ b/src/MOLECULE/atom_vec_bond.cpp @@ -23,8 +23,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecBond::AtomVecBond(LAMMPS *lmp) : AtomVec(lmp) @@ -47,13 +45,13 @@ AtomVecBond::AtomVecBond(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecBond::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/MOLECULE/atom_vec_full.cpp b/src/MOLECULE/atom_vec_full.cpp index f132fdc599..006917ce49 100644 --- a/src/MOLECULE/atom_vec_full.cpp +++ b/src/MOLECULE/atom_vec_full.cpp @@ -23,8 +23,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecFull::AtomVecFull(LAMMPS *lmp) : AtomVec(lmp) @@ -47,13 +45,13 @@ AtomVecFull::AtomVecFull(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecFull::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/MOLECULE/atom_vec_molecular.cpp b/src/MOLECULE/atom_vec_molecular.cpp index 766da9e354..3e347af5c9 100644 --- a/src/MOLECULE/atom_vec_molecular.cpp +++ b/src/MOLECULE/atom_vec_molecular.cpp @@ -23,8 +23,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecMolecular::AtomVecMolecular(LAMMPS *lmp) : AtomVec(lmp) @@ -47,13 +45,13 @@ AtomVecMolecular::AtomVecMolecular(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecMolecular::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/MOLECULE/atom_vec_template.cpp b/src/MOLECULE/atom_vec_template.cpp index a2bc632952..546b442825 100644 --- a/src/MOLECULE/atom_vec_template.cpp +++ b/src/MOLECULE/atom_vec_template.cpp @@ -25,8 +25,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecTemplate::AtomVecTemplate(LAMMPS *lmp) : AtomVec(lmp) @@ -90,13 +88,13 @@ void AtomVecTemplate::process_args(int narg, char **arg) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecTemplate::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/PERI/atom_vec_peri.cpp b/src/PERI/atom_vec_peri.cpp index c2791ce01f..165c9bc462 100644 --- a/src/PERI/atom_vec_peri.cpp +++ b/src/PERI/atom_vec_peri.cpp @@ -29,8 +29,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - static const char cite_peri_package[] = "PERI package for Peridynamics:\n\n" "@Article{Parks08,\n" @@ -66,13 +64,13 @@ AtomVecPeri::AtomVecPeri(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecPeri::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/USER-AWPMD/atom_vec_wavepacket.cpp b/src/USER-AWPMD/atom_vec_wavepacket.cpp index 9903e889fa..155ad88ca5 100644 --- a/src/USER-AWPMD/atom_vec_wavepacket.cpp +++ b/src/USER-AWPMD/atom_vec_wavepacket.cpp @@ -29,8 +29,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp) : AtomVec(lmp) @@ -58,13 +56,13 @@ AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom-electron arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecWavepacket::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; diff --git a/src/USER-CUDA/atom_vec_angle_cuda.cpp b/src/USER-CUDA/atom_vec_angle_cuda.cpp index 82a5b77dd2..2acc5d913e 100644 --- a/src/USER-CUDA/atom_vec_angle_cuda.cpp +++ b/src/USER-CUDA/atom_vec_angle_cuda.cpp @@ -52,7 +52,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 #define BUFFACTOR 1.5 #define BUFEXTRA 1000 #define NCUDAEXCHANGE 12 //nextra x y z vx vy vz tag type mask image molecule diff --git a/src/USER-CUDA/atom_vec_atomic_cuda.cpp b/src/USER-CUDA/atom_vec_atomic_cuda.cpp index bf8865b5cc..a2233fab6e 100644 --- a/src/USER-CUDA/atom_vec_atomic_cuda.cpp +++ b/src/USER-CUDA/atom_vec_atomic_cuda.cpp @@ -51,7 +51,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 #define BUFFACTOR 1.5 #define BUFEXTRA 1000 #define NCUDAEXCHANGE 11 //nextra x y z vx vy vz tag type mask image diff --git a/src/USER-CUDA/atom_vec_charge_cuda.cpp b/src/USER-CUDA/atom_vec_charge_cuda.cpp index f9cb3673e1..016e7fcca6 100644 --- a/src/USER-CUDA/atom_vec_charge_cuda.cpp +++ b/src/USER-CUDA/atom_vec_charge_cuda.cpp @@ -51,7 +51,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 #define BUFFACTOR 1.5 #define BUFEXTRA 1000 #define NCUDAEXCHANGE 12 //nextra x y z vx vy vz tag type mask image q diff --git a/src/USER-CUDA/atom_vec_full_cuda.cpp b/src/USER-CUDA/atom_vec_full_cuda.cpp index c430298cd6..a69c68c987 100644 --- a/src/USER-CUDA/atom_vec_full_cuda.cpp +++ b/src/USER-CUDA/atom_vec_full_cuda.cpp @@ -52,7 +52,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 #define BUFFACTOR 1.5 #define BUFEXTRA 1000 #define NCUDAEXCHANGE 13 //nextra x y z vx vy vz tag type mask image q molecule diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp index c68e348035..41837f1263 100644 --- a/src/USER-EFF/atom_vec_electron.cpp +++ b/src/USER-EFF/atom_vec_electron.cpp @@ -30,8 +30,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - static const char cite_user_eff_package[] = "USER-EFF package:\n\n" "@Article{Jaramillo-Botero11,\n" @@ -71,13 +69,13 @@ AtomVecElectron::AtomVecElectron(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom-electron arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecElectron::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp index 035f6aa90d..9b06f05a24 100644 --- a/src/USER-EFF/compute_temp_eff.cpp +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -23,8 +23,6 @@ #include "update.h" #include "force.h" #include "domain.h" -#include "modify.h" -#include "fix.h" #include "group.h" #include "error.h" @@ -58,9 +56,7 @@ ComputeTempEff::~ComputeTempEff() void ComputeTempEff::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -68,6 +64,7 @@ void ComputeTempEff::setup() void ComputeTempEff::dof_compute() { + if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; diff --git a/src/USER-MISC/compute_temp_rotate.cpp b/src/USER-MISC/compute_temp_rotate.cpp index 75950e2ca0..39f17d2d4f 100644 --- a/src/USER-MISC/compute_temp_rotate.cpp +++ b/src/USER-MISC/compute_temp_rotate.cpp @@ -23,8 +23,6 @@ #include "update.h" #include "force.h" #include "group.h" -#include "modify.h" -#include "fix.h" #include "domain.h" #include "lattice.h" #include "error.h" @@ -70,9 +68,7 @@ void ComputeTempRotate::init() void ComputeTempRotate::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -80,6 +76,7 @@ void ComputeTempRotate::setup() void ComputeTempRotate::dof_compute() { + if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = domain->dimension; dof = nper * natoms; diff --git a/src/USER-SPH/atom_vec_meso.cpp b/src/USER-SPH/atom_vec_meso.cpp index 808595bb35..3a57776d94 100644 --- a/src/USER-SPH/atom_vec_meso.cpp +++ b/src/USER-SPH/atom_vec_meso.cpp @@ -23,8 +23,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecMeso::AtomVecMeso(LAMMPS *lmp) : AtomVec(lmp) @@ -50,15 +48,14 @@ AtomVecMeso::AtomVecMeso(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ -void AtomVecMeso::grow(int n) { - if (n == 0) - nmax += DELTA; - else - nmax = n; +void AtomVecMeso::grow(int n) +{ + if (n == 0) grow_nmax(); + else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) error->one(FLERR,"Per-processor system is too big"); diff --git a/src/atom.cpp b/src/atom.cpp index cd471e581d..ba8283bc75 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -330,17 +330,12 @@ void Atom::create_avec(const char *style, int narg, char **arg, char *suffix) // create instance of AtomVec // use grow() to initialize atom-based arrays to length 1 // so that x[0][0] can always be referenced even if proc has no atoms - // but reset nmax = 0 - // so 2d arrays like bond_type will later be allocated correctly - // since currently, 2nd dimension bond_per_atom = 0 int sflag; avec = new_avec(style,suffix,sflag); avec->store_args(narg,arg); avec->process_args(narg,arg); avec->grow(1); - nmax = 0; - avec->reset(); if (sflag) { char estyle[256]; @@ -640,6 +635,48 @@ int Atom::count_words(const char *line) return n; } +/* ---------------------------------------------------------------------- + deallocate molecular topology arrays + done before realloc with (possibly) new 2nd dimension set to + correctly initialized per-atom values, e.g. bond_per_atom + needs to be called whenever 2nd dimensions are changed + and these arrays are already pre-allocated, + e.g. due to grow(1) in create_avec() +------------------------------------------------------------------------- */ + +void Atom::deallocate_topology() +{ + memory->destroy(atom->bond_type); + memory->destroy(atom->bond_atom); + atom->bond_type = NULL; + atom->bond_atom = NULL; + + memory->destroy(atom->angle_type); + memory->destroy(atom->angle_atom1); + memory->destroy(atom->angle_atom2); + memory->destroy(atom->angle_atom3); + atom->angle_type = NULL; + atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL; + + memory->destroy(atom->dihedral_type); + memory->destroy(atom->dihedral_atom1); + memory->destroy(atom->dihedral_atom2); + memory->destroy(atom->dihedral_atom3); + memory->destroy(atom->dihedral_atom4); + atom->dihedral_type = NULL; + atom->dihedral_atom1 = atom->dihedral_atom2 = + atom->dihedral_atom3 = atom->dihedral_atom4 = NULL; + + memory->destroy(atom->improper_type); + memory->destroy(atom->improper_atom1); + memory->destroy(atom->improper_atom2); + memory->destroy(atom->improper_atom3); + memory->destroy(atom->improper_atom4); + atom->improper_type = NULL; + atom->improper_atom1 = atom->improper_atom2 = + atom->improper_atom3 = atom->improper_atom4 = NULL; +} + /* ---------------------------------------------------------------------- unpack n lines from Atom section of data file call style-specific routine to parse line diff --git a/src/atom.h b/src/atom.h index 8138496052..5887b17c39 100644 --- a/src/atom.h +++ b/src/atom.h @@ -176,6 +176,8 @@ class Atom : protected Pointers { int parse_data(const char *); int count_words(const char *); + void deallocate_topology(); + void data_atoms(int, char *); void data_vels(int, char *); diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp index 24e9e2d713..0fd8043d85 100644 --- a/src/atom_vec.cpp +++ b/src/atom_vec.cpp @@ -21,6 +21,9 @@ using namespace LAMMPS_NS; +#define DELTA 16384 +#define DELTA_BONUS 8192 + /* ---------------------------------------------------------------------- */ AtomVec::AtomVec(LAMMPS *lmp) : Pointers(lmp) @@ -82,12 +85,24 @@ void AtomVec::init() } /* ---------------------------------------------------------------------- - reset nmax = 0, called by Atom::create_avec() when new atom style created + grow nmax so it is a multiple of DELTA ------------------------------------------------------------------------- */ -void AtomVec::reset() +void AtomVec::grow_nmax() { - nmax = 0; + nmax = nmax/DELTA * DELTA; + nmax += DELTA; +} + +/* ---------------------------------------------------------------------- + grow nmax_bonus so it is a multiple of DELTA_BONUS +------------------------------------------------------------------------- */ + +int AtomVec::grow_nmax_bonus(int nmax_bonus) +{ + nmax_bonus = nmax_bonus/DELTA_BONUS * DELTA_BONUS; + nmax_bonus += DELTA_BONUS; + return nmax_bonus; } /* ---------------------------------------------------------------------- diff --git a/src/atom_vec.h b/src/atom_vec.h index 241e84fadb..6597f55fa5 100644 --- a/src/atom_vec.h +++ b/src/atom_vec.h @@ -103,7 +103,6 @@ class AtomVec : protected Pointers { virtual void write_vel(FILE *, int, double **); virtual int write_vel_hybrid(FILE *, double *) {return 0;} - void reset(); int pack_bond(tagint **); void write_bond(FILE *, int, tagint **, int); int pack_angle(tagint **); @@ -139,6 +138,9 @@ class AtomVec : protected Pointers { ubuf(int64_t arg) : i(arg) {} ubuf(int arg) : i(arg) {} }; + + void grow_nmax(); + int grow_nmax_bonus(int); }; } diff --git a/src/atom_vec_atomic.cpp b/src/atom_vec_atomic.cpp index 5d96dc139e..c23d81e3cc 100644 --- a/src/atom_vec_atomic.cpp +++ b/src/atom_vec_atomic.cpp @@ -23,8 +23,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecAtomic::AtomVecAtomic(LAMMPS *lmp) : AtomVec(lmp) @@ -44,13 +42,13 @@ AtomVecAtomic::AtomVecAtomic(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecAtomic::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index 7058d1f34f..1461ecf501 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -28,9 +28,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 -#define DELTA_BONUS 10000 - /* ---------------------------------------------------------------------- */ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp) @@ -109,13 +106,13 @@ void AtomVecBody::process_args(int narg, char **arg) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecBody::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) @@ -158,7 +155,7 @@ void AtomVecBody::grow_reset() void AtomVecBody::grow_bonus() { - nmax_bonus += DELTA_BONUS; + nmax_bonus = grow_nmax_bonus(nmax_bonus); if (nmax_bonus < 0 || nmax_bonus > MAXSMALLINT) error->one(FLERR,"Per-processor system is too big"); diff --git a/src/atom_vec_charge.cpp b/src/atom_vec_charge.cpp index d8b4b7430c..da10ef42dd 100644 --- a/src/atom_vec_charge.cpp +++ b/src/atom_vec_charge.cpp @@ -23,8 +23,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecCharge::AtomVecCharge(LAMMPS *lmp) : AtomVec(lmp) @@ -46,13 +44,13 @@ AtomVecCharge::AtomVecCharge(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecCharge::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp index 512f3ba74b..81e4fbee20 100644 --- a/src/atom_vec_ellipsoid.cpp +++ b/src/atom_vec_ellipsoid.cpp @@ -31,9 +31,6 @@ using namespace LAMMPS_NS; using namespace MathConst; -#define DELTA 10000 -#define DELTA_BONUS 10000 - /* ---------------------------------------------------------------------- */ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp) : AtomVec(lmp) @@ -66,13 +63,13 @@ AtomVecEllipsoid::~AtomVecEllipsoid() /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecEllipsoid::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) @@ -115,7 +112,7 @@ void AtomVecEllipsoid::grow_reset() void AtomVecEllipsoid::grow_bonus() { - nmax_bonus += DELTA_BONUS; + nmax_bonus = grow_nmax_bonus(nmax_bonus); if (nmax_bonus < 0 || nmax_bonus > MAXSMALLINT) error->one(FLERR,"Per-processor system is too big"); diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp index b778e38782..88b4e6c43b 100644 --- a/src/atom_vec_hybrid.cpp +++ b/src/atom_vec_hybrid.cpp @@ -23,8 +23,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp) : AtomVec(lmp) {} @@ -135,13 +133,13 @@ void AtomVecHybrid::init() /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecHybrid::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/atom_vec_line.cpp b/src/atom_vec_line.cpp index bd5a5ea2b3..6760589ecb 100644 --- a/src/atom_vec_line.cpp +++ b/src/atom_vec_line.cpp @@ -26,8 +26,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 -#define DELTA_BONUS 10000 #define EPSILON 0.001 /* ---------------------------------------------------------------------- */ @@ -73,13 +71,13 @@ void AtomVecLine::init() /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecLine::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) @@ -124,7 +122,7 @@ void AtomVecLine::grow_reset() void AtomVecLine::grow_bonus() { - nmax_bonus += DELTA_BONUS; + nmax_bonus = grow_nmax_bonus(nmax_bonus); if (nmax_bonus < 0 || nmax_bonus > MAXSMALLINT) error->one(FLERR,"Per-processor system is too big"); diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp index acb1251981..de85a2eb15 100644 --- a/src/atom_vec_sphere.cpp +++ b/src/atom_vec_sphere.cpp @@ -29,8 +29,6 @@ using namespace LAMMPS_NS; using namespace MathConst; -#define DELTA 10000 - /* ---------------------------------------------------------------------- */ AtomVecSphere::AtomVecSphere(LAMMPS *lmp) : AtomVec(lmp) @@ -77,13 +75,13 @@ void AtomVecSphere::init() /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecSphere::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp index 20528ff0d6..f5f0f55f63 100644 --- a/src/atom_vec_tri.cpp +++ b/src/atom_vec_tri.cpp @@ -27,8 +27,6 @@ using namespace LAMMPS_NS; -#define DELTA 10000 -#define DELTA_BONUS 10000 #define EPSILON 0.001 /* ---------------------------------------------------------------------- */ @@ -74,13 +72,13 @@ void AtomVecTri::init() /* ---------------------------------------------------------------------- grow atom arrays - n = 0 grows arrays by DELTA + n = 0 grows arrays by a chunk n > 0 allocates arrays to size n ------------------------------------------------------------------------- */ void AtomVecTri::grow(int n) { - if (n == 0) nmax += DELTA; + if (n == 0) grow_nmax(); else nmax = n; atom->nmax = nmax; if (nmax < 0 || nmax > MAXSMALLINT) @@ -125,7 +123,7 @@ void AtomVecTri::grow_reset() void AtomVecTri::grow_bonus() { - nmax_bonus += DELTA_BONUS; + nmax_bonus = grow_nmax_bonus(nmax_bonus); if (nmax_bonus < 0 || nmax_bonus > MAXSMALLINT) error->one(FLERR,"Per-processor system is too big"); diff --git a/src/compute.cpp b/src/compute.cpp index 9367983656..fde11f6c7f 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -21,6 +21,8 @@ #include "domain.h" #include "comm.h" #include "group.h" +#include "modify.h" +#include "fix.h" #include "atom_masks.h" #include "memory.h" #include "error.h" @@ -131,6 +133,17 @@ void Compute::modify_params(int narg, char **arg) } } +/* ---------------------------------------------------------------------- + calculate adjustment in DOF due to fixes +------------------------------------------------------------------------- */ + +void Compute::adjust_dof_fix() +{ + fix_dof = 0; + for (int i = 0; i < modify->nfix; i++) + fix_dof += modify->fix[i]->dof(igroup); +} + /* ---------------------------------------------------------------------- reset extra_dof to its default value ------------------------------------------------------------------------- */ diff --git a/src/compute.h b/src/compute.h index 6d4bd3c0e3..d778564bfc 100644 --- a/src/compute.h +++ b/src/compute.h @@ -86,6 +86,7 @@ class Compute : protected Pointers { Compute(class LAMMPS *, int, char **); virtual ~Compute(); void modify_params(int, char **); + void adjust_dof_fix(); void reset_extra_dof(); virtual void init() = 0; @@ -124,6 +125,7 @@ class Compute : protected Pointers { protected: int extra_dof; // extra DOF for temperature computes + int fix_dof; // DOF due to fixes int dynamic; // recount atoms for temperature computes int thermoflag; // 1 if include fix PE for PE computes diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index c1e5e82eb0..3987c3d7ca 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -18,8 +18,6 @@ #include "update.h" #include "force.h" #include "domain.h" -#include "modify.h" -#include "fix.h" #include "group.h" #include "error.h" @@ -52,9 +50,7 @@ ComputeTemp::~ComputeTemp() void ComputeTemp::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -62,6 +58,7 @@ void ComputeTemp::setup() void ComputeTemp::dof_compute() { + if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= static_cast(extra_dof + fix_dof); diff --git a/src/compute_temp.h b/src/compute_temp.h index ef8a0465ec..cfe5012f81 100644 --- a/src/compute_temp.h +++ b/src/compute_temp.h @@ -34,7 +34,6 @@ class ComputeTemp : public Compute { void compute_vector(); protected: - int fix_dof; double tfactor; virtual void dof_compute(); diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp index 1d8b4c5304..536b11edb6 100644 --- a/src/compute_temp_com.cpp +++ b/src/compute_temp_com.cpp @@ -19,8 +19,6 @@ #include "update.h" #include "force.h" #include "group.h" -#include "modify.h" -#include "fix.h" #include "domain.h" #include "lattice.h" #include "error.h" @@ -62,9 +60,7 @@ void ComputeTempCOM::init() void ComputeTempCOM::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -72,6 +68,7 @@ void ComputeTempCOM::setup() void ComputeTempCOM::dof_compute() { + if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = domain->dimension; dof = nper * natoms; diff --git a/src/compute_temp_com.h b/src/compute_temp_com.h index 54c576d2d9..a8d1405cce 100644 --- a/src/compute_temp_com.h +++ b/src/compute_temp_com.h @@ -41,7 +41,6 @@ class ComputeTempCOM : public Compute { void restore_bias_thr(int, double *, double *); private: - int fix_dof; double tfactor,masstotal; void dof_compute(); diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index a54bd54b59..6171a818ce 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -86,9 +86,7 @@ void ComputeTempDeform::init() void ComputeTempDeform::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -96,6 +94,7 @@ void ComputeTempDeform::setup() void ComputeTempDeform::dof_compute() { + if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; diff --git a/src/compute_temp_deform.h b/src/compute_temp_deform.h index fd4d66a352..8f530b6332 100644 --- a/src/compute_temp_deform.h +++ b/src/compute_temp_deform.h @@ -42,7 +42,6 @@ class ComputeTempDeform : public Compute { double memory_usage(); protected: - int fix_dof; double tfactor; virtual void dof_compute(); diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp index 36c5dfacd1..90849ce7ed 100644 --- a/src/compute_temp_partial.cpp +++ b/src/compute_temp_partial.cpp @@ -18,8 +18,6 @@ #include "update.h" #include "force.h" #include "domain.h" -#include "modify.h" -#include "fix.h" #include "group.h" #include "memory.h" #include "error.h" @@ -63,9 +61,7 @@ ComputeTempPartial::~ComputeTempPartial() void ComputeTempPartial::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -76,6 +72,7 @@ void ComputeTempPartial::setup() void ComputeTempPartial::dof_compute() { + if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = xflag+yflag+zflag; dof = nper * natoms; diff --git a/src/compute_temp_partial.h b/src/compute_temp_partial.h index 5696336a9a..e43df3f6e5 100644 --- a/src/compute_temp_partial.h +++ b/src/compute_temp_partial.h @@ -44,7 +44,6 @@ class ComputeTempPartial : public Compute { protected: int xflag,yflag,zflag; - int fix_dof; double tfactor; void dof_compute(); diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index ecfad04046..75b26b95bf 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -161,9 +161,7 @@ ComputeTempProfile::~ComputeTempProfile() void ComputeTempProfile::init() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); // ptrs to domain data @@ -189,9 +187,7 @@ void ComputeTempProfile::init() void ComputeTempProfile::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -199,6 +195,7 @@ void ComputeTempProfile::setup() void ComputeTempProfile::dof_compute() { + if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = domain->dimension; dof = nper * natoms; diff --git a/src/compute_temp_profile.h b/src/compute_temp_profile.h index e43d538c63..c39974f746 100644 --- a/src/compute_temp_profile.h +++ b/src/compute_temp_profile.h @@ -46,7 +46,6 @@ class ComputeTempProfile : public Compute { int xflag,yflag,zflag,ncount,outflag; int nbinx,nbiny,nbinz,nbins; int ivx,ivy,ivz; - int fix_dof; double tfactor; int box_change,triclinic; diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index a4c9f4ec23..f5f2165501 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -117,9 +117,7 @@ ComputeTempRamp::~ComputeTempRamp() void ComputeTempRamp::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -127,6 +125,7 @@ void ComputeTempRamp::setup() void ComputeTempRamp::dof_compute() { + if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = domain->dimension; dof = nper * natoms; diff --git a/src/compute_temp_ramp.h b/src/compute_temp_ramp.h index a669bc48a5..9914645d23 100644 --- a/src/compute_temp_ramp.h +++ b/src/compute_temp_ramp.h @@ -46,7 +46,7 @@ class ComputeTempRamp : public Compute { double coord_lo,coord_hi; int v_dim; double v_lo,v_hi; - int scaleflag,fix_dof; + int scaleflag; double tfactor,xscale,yscale,zscale; void dof_compute(); diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index 2777463a11..d3fc4f0ce2 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -17,7 +17,6 @@ #include "atom.h" #include "update.h" #include "force.h" -#include "modify.h" #include "domain.h" #include "region.h" #include "memory.h" diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index 0bae77b3a7..821828a3de 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -20,7 +20,6 @@ #include "force.h" #include "domain.h" #include "modify.h" -#include "fix.h" #include "group.h" #include "error.h" @@ -109,9 +108,7 @@ void ComputeTempSphere::init() void ComputeTempSphere::setup() { - fix_dof = 0; - for (int i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); + fix_dof = -1; dof_compute(); } @@ -121,6 +118,8 @@ void ComputeTempSphere::dof_compute() { int count,count_all; + if (fix_dof) adjust_dof_fix(); + // 6 or 3 dof for extended/point particles for 3d // 3 or 2 dof for extended/point particles for 2d // which dof are included also depends on mode diff --git a/src/compute_temp_sphere.h b/src/compute_temp_sphere.h index 8645e3ea58..7ba74832e2 100644 --- a/src/compute_temp_sphere.h +++ b/src/compute_temp_sphere.h @@ -39,7 +39,7 @@ class ComputeTempSphere : public Compute { void restore_bias_thr(int, double *, double *); private: - int fix_dof,mode; + int mode; double tfactor; char *id_bias; Compute *tbias; // ptr to additional bias compute diff --git a/src/create_box.cpp b/src/create_box.cpp index 1e7ce40122..21e20aab56 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -164,11 +164,14 @@ void CreateBox::command(int narg, char **arg) } // problem setup using info from header - // no call to atom->grow since create_atoms or fixes will do it + // deallocate/grow insures any extra settings are used for topology arrays + // necessary in case no create_atoms is performed update->ntimestep = 0; atom->allocate_type_arrays(); + atom->deallocate_topology(); + atom->avec->grow(1); domain->print_box("Created "); domain->set_initial_box(); diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 3f04cb8254..d7bef2293b 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -370,7 +370,7 @@ void FixAdapt::change_settings() } /* ---------------------------------------------------------------------- - restore pair,kspace.atom parameters to original values + restore pair,kspace,atom parameters to original values ------------------------------------------------------------------------- */ void FixAdapt::restore_settings() diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp index f530731ecb..db2a78f2f6 100644 --- a/src/fix_addforce.cpp +++ b/src/fix_addforce.cpp @@ -104,8 +104,8 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) : force_flag = 0; foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0; - maxatom = -1; - sforce = NULL; + maxatom = atom->nmax; + memory->create(sforce,maxatom,4,"addforce:sforce"); } /* ---------------------------------------------------------------------- */ @@ -274,24 +274,16 @@ void FixAddForce::post_force(int vflag) modify->clearstep_compute(); if (xstyle == EQUAL) xvalue = input->variable->compute_equal(xvar); - else if (xstyle == ATOM) { - if (sforce) input->variable->compute_atom(xvar,igroup,&sforce[0][0],4,0); - else input->variable->compute_atom(xvar,igroup,NULL,4,0); - } + else if (xstyle == ATOM) + input->variable->compute_atom(xvar,igroup,&sforce[0][0],4,0); if (ystyle == EQUAL) yvalue = input->variable->compute_equal(yvar); - else if (ystyle == ATOM) { - if (sforce) input->variable->compute_atom(yvar,igroup,&sforce[0][1],4,0); - else input->variable->compute_atom(yvar,igroup,NULL,4,0); - } + else if (ystyle == ATOM) + input->variable->compute_atom(yvar,igroup,&sforce[0][1],4,0); if (zstyle == EQUAL) zvalue = input->variable->compute_equal(zvar); - else if (zstyle == ATOM) { - if (sforce) input->variable->compute_atom(zvar,igroup,&sforce[0][2],4,0); - else input->variable->compute_atom(zvar,igroup,NULL,4,0); - } - if (estyle == ATOM) { - if (sforce) input->variable->compute_atom(evar,igroup,&sforce[0][3],4,0); - else input->variable->compute_atom(evar,igroup,NULL,4,0); - } + else if (zstyle == ATOM) + input->variable->compute_atom(zvar,igroup,&sforce[0][2],4,0); + if (estyle == ATOM) + input->variable->compute_atom(evar,igroup,&sforce[0][3],4,0); modify->addstep_compute(update->ntimestep + 1); diff --git a/src/fix_move.cpp b/src/fix_move.cpp index 0169f4882d..0d80f5664f 100644 --- a/src/fix_move.cpp +++ b/src/fix_move.cpp @@ -242,8 +242,9 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : atom->add_callback(0); atom->add_callback(1); - maxatom = -1; - displace = velocity = NULL; + maxatom = atom->nmax; + if (displaceflag) memory->create(displace,maxatom,3,"move:displace"); + if (velocityflag) memory->create(velocity,maxatom,3,"move:velocity"); // xoriginal = initial unwrapped positions of atoms @@ -581,45 +582,27 @@ void FixMove::initial_integrate(int vflag) if (xvarstr) { if (xvarstyle == EQUAL) dx = input->variable->compute_equal(xvar); - else if (displace) - input->variable->compute_atom(xvar,igroup,&displace[0][0],3,0); - else - input->variable->compute_atom(xvar,igroup,NULL,3,0); + else input->variable->compute_atom(xvar,igroup,&displace[0][0],3,0); } if (yvarstr) { if (yvarstyle == EQUAL) dy = input->variable->compute_equal(yvar); - else if (displace) - input->variable->compute_atom(yvar,igroup,&displace[0][1],3,0); - else - input->variable->compute_atom(yvar,igroup,NULL,3,0); + else input->variable->compute_atom(yvar,igroup,&displace[0][1],3,0); } if (zvarstr) { if (zvarstyle == EQUAL) dz = input->variable->compute_equal(zvar); - else if (displace) - input->variable->compute_atom(zvar,igroup,&displace[0][2],3,0); - else - input->variable->compute_atom(zvar,igroup,NULL,3,0); + else input->variable->compute_atom(zvar,igroup,&displace[0][2],3,0); } if (vxvarstr) { if (vxvarstyle == EQUAL) vx = input->variable->compute_equal(vxvar); - else if (velocity) - input->variable->compute_atom(vxvar,igroup,&velocity[0][0],3,0); - else - input->variable->compute_atom(vxvar,igroup,NULL,3,0); + else input->variable->compute_atom(vxvar,igroup,&velocity[0][0],3,0); } if (vyvarstr) { if (vyvarstyle == EQUAL) vy = input->variable->compute_equal(vyvar); - else if (velocity) - input->variable->compute_atom(vyvar,igroup,&velocity[0][1],3,0); - else - input->variable->compute_atom(vyvar,igroup,NULL,3,0); + else input->variable->compute_atom(vyvar,igroup,&velocity[0][1],3,0); } if (vzvarstr) { if (vzvarstyle == EQUAL) vz = input->variable->compute_equal(vzvar); - else if (velocity) - input->variable->compute_atom(vzvar,igroup,&velocity[0][2],3,0); - else - input->variable->compute_atom(vzvar,igroup,NULL,3,0); + else input->variable->compute_atom(vzvar,igroup,&velocity[0][2],3,0); } modify->addstep_compute(update->ntimestep + 1); diff --git a/src/fix_setforce.cpp b/src/fix_setforce.cpp index 0d504f3b78..a6e1ecb972 100644 --- a/src/fix_setforce.cpp +++ b/src/fix_setforce.cpp @@ -99,8 +99,8 @@ FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) : force_flag = 0; foriginal[0] = foriginal[1] = foriginal[2] = 0.0; - maxatom = -1; - sforce = NULL; + maxatom = atom->nmax; + memory->create(sforce,maxatom,3,"setforce:sforce"); } /* ---------------------------------------------------------------------- */ @@ -257,20 +257,14 @@ void FixSetForce::post_force(int vflag) modify->clearstep_compute(); if (xstyle == EQUAL) xvalue = input->variable->compute_equal(xvar); - else if (xstyle == ATOM) { - if (sforce) input->variable->compute_atom(xvar,igroup,&sforce[0][0],3,0); - else input->variable->compute_atom(xvar,igroup,NULL,3,0); - } + else if (xstyle == ATOM) + input->variable->compute_atom(xvar,igroup,&sforce[0][0],3,0); if (ystyle == EQUAL) yvalue = input->variable->compute_equal(yvar); - else if (ystyle == ATOM) { - if (sforce) input->variable->compute_atom(yvar,igroup,&sforce[0][1],3,0); - else input->variable->compute_atom(yvar,igroup,NULL,3,0); - } + else if (ystyle == ATOM) + input->variable->compute_atom(yvar,igroup,&sforce[0][1],3,0); if (zstyle == EQUAL) zvalue = input->variable->compute_equal(zvar); - else if (zstyle == ATOM) { - if (sforce) input->variable->compute_atom(zvar,igroup,&sforce[0][2],3,0); - else input->variable->compute_atom(zvar,igroup,NULL,3,0); - } + else if (zstyle == ATOM) + input->variable->compute_atom(zvar,igroup,&sforce[0][2],3,0); modify->addstep_compute(update->ntimestep + 1); diff --git a/src/fix_store_state.cpp b/src/fix_store_state.cpp index 014c9a773d..e70eada845 100644 --- a/src/fix_store_state.cpp +++ b/src/fix_store_state.cpp @@ -467,11 +467,9 @@ void FixStoreState::end_of_step() // evaluate atom-style variable - } else if (which[m] == VARIABLE) - if (values) - input->variable->compute_atom(n,igroup,&values[0][m],nvalues,0); - else - input->variable->compute_atom(n,igroup,NULL,nvalues,0); + } else if (which[m] == VARIABLE) { + input->variable->compute_atom(n,igroup,&values[0][m],nvalues,0); + } } } diff --git a/src/modify.cpp b/src/modify.cpp index c454ff928c..8a193253f2 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -271,7 +271,7 @@ void Modify::init() void Modify::setup(int vflag) { // compute setup needs to come before fix setup - // b/c NH fixes need use DOF of temperature computes + // b/c NH fixes need DOF of temperature computes for (int i = 0; i < ncompute; i++) compute[i]->setup(); diff --git a/src/read_data.cpp b/src/read_data.cpp index 08551e99d1..0f54fa3394 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -206,9 +206,9 @@ void ReadData::command(int narg, char **arg) domain->box_exist = 1; update->ntimestep = 0; - // insure extra settings are applied before grow(), - // even if no topology in file - // if topology is in file, realloc and another grow() is done below + // apply extra settings before grow(), even if no topology in file + // deallocate() insures new settings are used for topology arrays + // if per-atom topology is in file, another grow() is done below atom->bond_per_atom = atom->extra_bond_per_atom; atom->angle_per_atom = atom->extra_angle_per_atom; @@ -220,6 +220,7 @@ void ReadData::command(int narg, char **arg) else n = static_cast (LB_FACTOR * atom->natoms / comm->nprocs); atom->allocate_type_arrays(); + atom->deallocate_topology(); atom->avec->grow(n); domain->print_box(" "); @@ -497,49 +498,11 @@ void ReadData::command(int narg, char **arg) firstpass = 0; // reallocate bond,angle,diehdral,improper arrays via grow() - // use new bond,angle,dihedral,improper per-atom values from 1st pass - // should leave other atom arrays unchanged, since already nmax in length - // if bonds/etc not in data file, initialize per-atom size - // with extra settings before grow() of these topology arrays - - if (bondflag) { - memory->destroy(atom->bond_type); - memory->destroy(atom->bond_atom); - atom->bond_type = NULL; - atom->bond_atom = NULL; - } - - if (angleflag) { - memory->destroy(atom->angle_type); - memory->destroy(atom->angle_atom1); - memory->destroy(atom->angle_atom2); - memory->destroy(atom->angle_atom3); - atom->angle_type = NULL; - atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL; - } - - if (dihedralflag) { - memory->destroy(atom->dihedral_type); - memory->destroy(atom->dihedral_atom1); - memory->destroy(atom->dihedral_atom2); - memory->destroy(atom->dihedral_atom3); - memory->destroy(atom->dihedral_atom4); - atom->dihedral_type = NULL; - atom->dihedral_atom1 = atom->dihedral_atom2 = - atom->dihedral_atom3 = atom->dihedral_atom4 = NULL; - } - - if (improperflag) { - memory->destroy(atom->improper_type); - memory->destroy(atom->improper_atom1); - memory->destroy(atom->improper_atom2); - memory->destroy(atom->improper_atom3); - memory->destroy(atom->improper_atom4); - atom->improper_type = NULL; - atom->improper_atom1 = atom->improper_atom2 = - atom->improper_atom3 = atom->improper_atom4 = NULL; - } + // will use new bond,angle,dihedral,improper per-atom values from 1st pass + // will also observe extra settings even if bond/etc topology not in file + // leaves other atom arrays unchanged, since already nmax in length + atom->deallocate_topology(); atom->avec->grow(atom->nmax); } diff --git a/src/read_restart.cpp b/src/read_restart.cpp index f72ced2c53..a224bfa2e3 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -151,6 +151,7 @@ void ReadRestart::command(int narg, char **arg) else n = static_cast (LB_FACTOR * atom->natoms / nprocs); atom->allocate_type_arrays(); + atom->deallocate_topology(); atom->avec->grow(n); n = atom->nmax; diff --git a/src/version.h b/src/version.h index 61569290c2..08a2177cc9 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "8 May 2014" +#define LAMMPS_VERSION "9 May 2014"