From 58336dca9a279d7cfa35e521585d16b59b0b1ab1 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 14:13:43 +0000 Subject: [PATCH 01/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11773 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/fix_bond_create.html | 26 ++++++++++++-------------- doc/fix_bond_create.txt | 26 ++++++++++++-------------- 2 files changed, 24 insertions(+), 28 deletions(-) diff --git a/doc/fix_bond_create.html b/doc/fix_bond_create.html index 079bc359e8..3098617f90 100644 --- a/doc/fix_bond_create.html +++ b/doc/fix_bond_create.html @@ -117,22 +117,20 @@ structures that store this information must have space for it. When LAMMPS is initialized from a data file, the list of bonds is scanned and the maximum number of bonds per atom is tallied. If some atom will acquire more bonds than this limit as this fix operates, then the -"extra bonds per atom" parameter in the data file header must be set -to allow for it. See the read_data command for more -details. Note that if this parameter needs to be set, it means a data -file must be used to initialize the system, even if it initially has -no bonds. A data file with no atoms can be used if you wish to add -unbonded atoms via the create atoms command, -e.g. for a percolation simulation. +"extra bonds per atom" parameter must be set to allow for it. See the +read_data or create_box command for +more details. Note that a data file with no atoms can be used if you +wish to add unbonded atoms via the create atoms +command, e.g. for a percolation simulation.

IMPORTANT NOTE: LAMMPS also maintains a data structure that stores a -list of 1st, 2nd, and 3rd neighbors of each atom (in the bond topology -of the system) for use in weighting pairwise interactions for bonded -atoms. Adding a bond adds a single entry to this list. The "extra" -keyword of the special_bonds command should be -used to leave space for new bonds if the maximum number of entries for -any atom will be exceeded as this fix operates. See the -special_bonds command for details. +list of 1st, 2nd, and 3rd neighbors of each atom (within the bond +topology of the system) for use in weighting pairwise interactions for +bonded atoms. Adding a bond adds a single entry to this list. The +"extra special per atom" parameter must typically be set to allow for +it. There are 3 ways to do this. See the read_data +or create_box or "special_bonds extra" commands for +details.

Note that even if your simulation starts with no bonds, you must define a bond_style and use the diff --git a/doc/fix_bond_create.txt b/doc/fix_bond_create.txt index 89fb3a0c88..c8bbdab829 100755 --- a/doc/fix_bond_create.txt +++ b/doc/fix_bond_create.txt @@ -105,22 +105,20 @@ structures that store this information must have space for it. When LAMMPS is initialized from a data file, the list of bonds is scanned and the maximum number of bonds per atom is tallied. If some atom will acquire more bonds than this limit as this fix operates, then the -"extra bonds per atom" parameter in the data file header must be set -to allow for it. See the "read_data"_read_data.html command for more -details. Note that if this parameter needs to be set, it means a data -file must be used to initialize the system, even if it initially has -no bonds. A data file with no atoms can be used if you wish to add -unbonded atoms via the "create atoms"_create_atoms.html command, -e.g. for a percolation simulation. +"extra bonds per atom" parameter must be set to allow for it. See the +"read_data"_read_data.html or "create_box"_create_box.html command for +more details. Note that a data file with no atoms can be used if you +wish to add unbonded atoms via the "create atoms"_create_atoms.html +command, e.g. for a percolation simulation. IMPORTANT NOTE: LAMMPS also maintains a data structure that stores a -list of 1st, 2nd, and 3rd neighbors of each atom (in the bond topology -of the system) for use in weighting pairwise interactions for bonded -atoms. Adding a bond adds a single entry to this list. The "extra" -keyword of the "special_bonds"_special_bonds.html command should be -used to leave space for new bonds if the maximum number of entries for -any atom will be exceeded as this fix operates. See the -"special_bonds"_special_bonds.html command for details. +list of 1st, 2nd, and 3rd neighbors of each atom (within the bond +topology of the system) for use in weighting pairwise interactions for +bonded atoms. Adding a bond adds a single entry to this list. The +"extra special per atom" parameter must typically be set to allow for +it. There are 3 ways to do this. See the "read_data"_read_data.html +or "create_box"_create_box.html or "special_bonds extra" commands for +details. Note that even if your simulation starts with no bonds, you must define a "bond_style"_bond_style.html and use the From 9065c2c612913456ccc4b6cacd4fa79f37761565 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 14:29:55 +0000 Subject: [PATCH 02/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11774 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/fix_property_atom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp index 0c573fba10..f74e135125 100644 --- a/src/fix_property_atom.cpp +++ b/src/fix_property_atom.cpp @@ -429,7 +429,7 @@ int FixPropertyAtom::pack_border(int n, int *list, double *buf) int *ivector = atom->ivector[index[k]]; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = ubuf(ivector[j]).i; + buf[m++] = ubuf(ivector[j]).d; } } else if (style[k] == DOUBLE) { double *dvector = atom->dvector[index[k]]; From 1b6af8196fc1e058bdbf4e5256d8031cce029d8d Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 14:41:25 +0000 Subject: [PATCH 03/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11775 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MISC/compute_ti.cpp | 4 ++-- src/Make.py | 2 +- src/Makefile | 14 +++++++------- src/USER-OMP/fix_qeq_comb_omp.cpp | 6 +----- src/USER-OMP/neigh_full_omp.cpp | 6 +++--- src/USER-OMP/neigh_gran_omp.cpp | 8 ++++---- src/USER-OMP/neigh_half_bin_omp.cpp | 8 ++++---- src/USER-OMP/neigh_half_multi_omp.cpp | 6 +++--- src/USER-OMP/neigh_respa_omp.cpp | 18 +++++++++--------- src/compute_property_atom.cpp | 20 ++++++++++++++++++++ src/compute_property_atom.h | 2 ++ 11 files changed, 56 insertions(+), 38 deletions(-) diff --git a/src/MISC/compute_ti.cpp b/src/MISC/compute_ti.cpp index e47bbda556..90083b0666 100644 --- a/src/MISC/compute_ti.cpp +++ b/src/MISC/compute_ti.cpp @@ -127,7 +127,7 @@ void ComputeTI::init() for (int m = 0; m < nterms; m++) { ivar1[m] = input->variable->find(var1[m]); ivar2[m] = input->variable->find(var2[m]); - if (ivar1[m] < 0 || ivar2 < 0) + if (ivar1[m] < 0 || ivar2[m] < 0) error->all(FLERR,"Variable name for compute ti does not exist"); if (!input->variable->equalstyle(ivar1[m]) || !input->variable->equalstyle(ivar2[m])) @@ -185,7 +185,7 @@ double ComputeTI::compute_scalar() double *eatom = pptr[m]->eatom; - if (force->newton) npair += atom->nghost; + if (force->newton_pair) npair += atom->nghost; for (int i = 0; i < npair; i++) if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) eng += eatom[i]; MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/Make.py b/src/Make.py index d02620aa35..c704cb7037 100755 --- a/src/Make.py +++ b/src/Make.py @@ -73,7 +73,7 @@ support = ["Makefile","Make.sh","Makefile.package.empty", extlibs = {"USER-ATC": "atc", "USER-AWPMD": "awpmd", "USER-COLVARS": "colvars", "USER-CUDA": "cuda","GPU": "gpu","MEAM": "meam", "POEMS": "poems", - "REAX": "reax"} + "USER-QMMM": "qmmm", "REAX": "reax"} # help messages diff --git a/src/Makefile b/src/Makefile index 555b87e8bc..019bb4bed9 100755 --- a/src/Makefile +++ b/src/Makefile @@ -22,7 +22,7 @@ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \ user-omp user-phonon user-qmmm user-reaxc user-sph PACKLIB = gpu kim meam poems reax voronoi \ - user-atc user-awpmd user-colvars user-cuda user-molfile + user-atc user-awpmd user-colvars user-qmmm user-cuda user-molfile PACKALL = $(PACKAGE) $(PACKUSER) @@ -50,7 +50,7 @@ help: @echo 'make install-python install LAMMPS wrapper in Python' @echo '' @echo 'make package list available packages' - @echo 'make package-status status of all packages' + @echo 'make package-status (ps) status of all packages' @echo 'make yes-package install a single package in src dir' @echo 'make no-package remove a single package from src dir' @echo 'make yes-all install all packages in src dir' @@ -61,9 +61,9 @@ help: @echo 'make no-user remove all user packages' @echo 'make no-lib remove all packages with external libs' @echo '' - @echo 'make package-update replace src files with package files' + @echo 'make package-update (pu) replace src files with updated package files' @echo 'make package-overwrite replace package files with src files' - @echo 'make package-diff diff src files against package files' + @echo 'make package-diff (pd) diff src files against package files' @echo '' @echo 'make machine build LAMMPS where machine is one of:' @echo '' @@ -108,18 +108,18 @@ purge: Purge.list # Create a tarball of src dir and packages tar: - @cd STUBS; make clean + @cd STUBS; $(MAKE) clean @cd ..; tar cvzf src/$(ROOT)_src.tar.gz \ src/Make* src/Package.sh src/MAKE src/*.cpp src/*.h src/STUBS \ $(patsubst %,src/%,$(PACKAGEUC)) $(patsubst %,src/%,$(PACKUSERUC)) \ --exclude=*/.svn - @cd STUBS; make + @cd STUBS; $(MAKE) @echo "Created $(ROOT)_src.tar.gz" # Make MPI STUBS library stubs: - @cd STUBS; make clean; make + @cd STUBS; $(MAKE) clean; $(MAKE) # Create Makefile.lib, Makefile.shlib, and Makefile.list diff --git a/src/USER-OMP/fix_qeq_comb_omp.cpp b/src/USER-OMP/fix_qeq_comb_omp.cpp index f29ba5fb3d..e5de887cd3 100644 --- a/src/USER-OMP/fix_qeq_comb_omp.cpp +++ b/src/USER-OMP/fix_qeq_comb_omp.cpp @@ -54,7 +54,7 @@ void FixQEQCombOMP::init() if (!atom->q_flag) error->all(FLERR,"Fix qeq/comb/omp requires atom attribute q"); - if (NULL == force->pair_match("comb3",0)) + if (NULL != force->pair_match("comb3",0)) error->all(FLERR,"No support for comb3 currently available in USER-OMP"); comb = (PairComb *) force->pair_match("comb/omp",1); @@ -79,10 +79,6 @@ void FixQEQCombOMP::init() } int irequest = neighbor->request(this); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->fix = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; neighbor->requests[irequest]->omp = use_omp; } diff --git a/src/USER-OMP/neigh_full_omp.cpp b/src/USER-OMP/neigh_full_omp.cpp index f8036727dd..ee36e5610e 100644 --- a/src/USER-OMP/neigh_full_omp.cpp +++ b/src/USER-OMP/neigh_full_omp.cpp @@ -255,7 +255,7 @@ void Neighbor::full_bin_omp(NeighList *list) { // bin owned & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; @@ -367,7 +367,7 @@ void Neighbor::full_bin_ghost_omp(NeighList *list) { // bin owned & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = atom->nlocal; const int nall = nlocal + atom->nghost; @@ -510,7 +510,7 @@ void Neighbor::full_multi_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; diff --git a/src/USER-OMP/neigh_gran_omp.cpp b/src/USER-OMP/neigh_gran_omp.cpp index 36f07fa095..55ab604b02 100644 --- a/src/USER-OMP/neigh_gran_omp.cpp +++ b/src/USER-OMP/neigh_gran_omp.cpp @@ -96,7 +96,7 @@ void Neighbor::granular_nsq_no_newton_omp(NeighList *list) touchptr = ipage_touch->vget(); shearptr = dpage_shear->vget(); } - + xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -271,7 +271,7 @@ void Neighbor::granular_bin_no_newton_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; @@ -424,7 +424,7 @@ void Neighbor::granular_bin_newton_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; @@ -533,7 +533,7 @@ void Neighbor::granular_bin_newton_tri_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; diff --git a/src/USER-OMP/neigh_half_bin_omp.cpp b/src/USER-OMP/neigh_half_bin_omp.cpp index 3974d9f6ab..65fa07a48a 100644 --- a/src/USER-OMP/neigh_half_bin_omp.cpp +++ b/src/USER-OMP/neigh_half_bin_omp.cpp @@ -35,7 +35,7 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; @@ -152,7 +152,7 @@ void Neighbor::half_bin_no_newton_ghost_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = atom->nlocal; const int nall = nlocal + atom->nghost; @@ -299,7 +299,7 @@ void Neighbor::half_bin_newton_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; @@ -447,7 +447,7 @@ void Neighbor::half_bin_newton_tri_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; diff --git a/src/USER-OMP/neigh_half_multi_omp.cpp b/src/USER-OMP/neigh_half_multi_omp.cpp index ee5c3cf6f4..cc93bf6367 100644 --- a/src/USER-OMP/neigh_half_multi_omp.cpp +++ b/src/USER-OMP/neigh_half_multi_omp.cpp @@ -36,7 +36,7 @@ void Neighbor::half_multi_no_newton_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; @@ -158,7 +158,7 @@ void Neighbor::half_multi_newton_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; @@ -314,7 +314,7 @@ void Neighbor::half_multi_newton_tri_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; diff --git a/src/USER-OMP/neigh_respa_omp.cpp b/src/USER-OMP/neigh_respa_omp.cpp index 6ac70cf9d7..409090c9ea 100644 --- a/src/USER-OMP/neigh_respa_omp.cpp +++ b/src/USER-OMP/neigh_respa_omp.cpp @@ -143,7 +143,7 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list) tag[j]-tagprev); else which = 0; if (which == 0) neighptr[n++] = j; - else if (minchange = domain->minimum_image_check(delx,dely,delz)) + else if ((minchange = domain->minimum_image_check(delx,dely,delz))) neighptr[n++] = j; else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; @@ -328,7 +328,7 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list) tag[j]-tagprev); else which = 0; if (which == 0) neighptr[n++] = j; - else if (minchange = domain->minimum_image_check(delx,dely,delz)) + else if ((minchange = domain->minimum_image_check(delx,dely,delz))) neighptr[n++] = j; else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; @@ -390,7 +390,7 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; @@ -506,7 +506,7 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) tag[j]-tagprev); else which = 0; if (which == 0) neighptr[n++] = j; - else if (minchange = domain->minimum_image_check(delx,dely,delz)) + else if ((minchange = domain->minimum_image_check(delx,dely,delz))) neighptr[n++] = j; else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; @@ -569,7 +569,7 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; @@ -688,7 +688,7 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) tag[j]-tagprev); else which = 0; if (which == 0) neighptr[n++] = j; - else if (minchange = domain->minimum_image_check(delx,dely,delz)) + else if ((minchange = domain->minimum_image_check(delx,dely,delz))) neighptr[n++] = j; else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; @@ -732,7 +732,7 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) tag[j]-tagprev); else which = 0; if (which == 0) neighptr[n++] = j; - else if (minchange = domain->minimum_image_check(delx,dely,delz)) + else if ((minchange = domain->minimum_image_check(delx,dely,delz))) neighptr[n++] = j; else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; @@ -795,7 +795,7 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) { // bin local & ghost atoms - bin_atoms(); + if (binatomflag) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; @@ -919,7 +919,7 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) tag[j]-tagprev); else which = 0; if (which == 0) neighptr[n++] = j; - else if (minchange = domain->minimum_image_check(delx,dely,delz)) + else if ((minchange = domain->minimum_image_check(delx,dely,delz))) neighptr[n++] = j; else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); } else neighptr[n++] = j; diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp index c9e6b126e2..23d8fb8226 100644 --- a/src/compute_property_atom.cpp +++ b/src/compute_property_atom.cpp @@ -337,6 +337,11 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) : "atom property that isn't allocated"); pack_choice[i] = &ComputePropertyAtom::pack_corner3z; + } else if (strcmp(arg[iarg],"nbonds") == 0) { + if (!atom->molecule_flag) + error->all(FLERR,"Compute property/atom for " + "atom property that isn't allocated"); + pack_choice[i] = &ComputePropertyAtom::pack_nbonds; } else if (strstr(arg[iarg],"i_") == arg[iarg]) { int flag; index[i] = atom->find_custom(&arg[iarg][2],flag); @@ -1687,6 +1692,21 @@ void ComputePropertyAtom::pack_corner3z(int n) /* ---------------------------------------------------------------------- */ +void ComputePropertyAtom::pack_nbonds(int n) +{ + int *num_bond = atom->num_bond; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = num_bond[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void ComputePropertyAtom::pack_iname(int n) { int *ivector = atom->ivector[index[n]]; diff --git a/src/compute_property_atom.h b/src/compute_property_atom.h index a5d2c5832d..01e623011d 100644 --- a/src/compute_property_atom.h +++ b/src/compute_property_atom.h @@ -121,6 +121,8 @@ class ComputePropertyAtom : public Compute { void pack_corner3y(int); void pack_corner3z(int); + void pack_nbonds(int); + void pack_iname(int); void pack_dname(int); }; From d45cca6714a18e35e9f7b6796a8b0cca00d20581 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 14:42:32 +0000 Subject: [PATCH 04/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11776 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- examples/comb/in.comb.Cu2O.elastic | 8 ++++---- examples/comb/in.comb.Si.elastic | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/comb/in.comb.Cu2O.elastic b/examples/comb/in.comb.Cu2O.elastic index 491c0b9ed5..d9df9ccfc3 100644 --- a/examples/comb/in.comb.Cu2O.elastic +++ b/examples/comb/in.comb.Cu2O.elastic @@ -55,7 +55,7 @@ unfix 2 ## calculates C11, C12, C13 and C14 fix 2 all deform 1 x scale 1.0001 remap x -compute perfx all stress/atom pair +compute perfx all stress/atom NULL pair compute fx all reduce sum & c_perfx[1] c_perfx[2] c_perfx[3] & c_perfx[4] c_perfx[5] c_perfx[6] @@ -66,7 +66,7 @@ run 10 ## strain z, measure s_z: calculates C33 fix 2 all deform 1 z scale 1.0001 remap x -compute perfz all stress/atom pair +compute perfz all stress/atom NULL pair compute fz all reduce sum & c_perfz[1] c_perfz[2] c_perfz[3] & c_perfz[4] c_perfz[5] c_perfz[6] @@ -77,7 +77,7 @@ run 10 ## strain yz, measure s_yz: calculates C44 fix 2 all deform 1 yz erate 0.0001 remap x -compute perfyz all stress/atom pair +compute perfyz all stress/atom NULL pair compute fyz all reduce sum & c_perfyz[1] c_perfyz[2] c_perfyz[3] & c_perfyz[4] c_perfyz[5] c_perfyz[6] @@ -88,7 +88,7 @@ run 10 ## strain xy, measure s_xy: calculates C66 fix 2 all deform 1 xy erate 0.0001 remap x -compute perfxy all stress/atom pair +compute perfxy all stress/atom NULL pair compute fxy all reduce sum & c_perfxy[1] c_perfxy[2] c_perfxy[3] & c_perfxy[4] c_perfxy[5] c_perfxy[6] diff --git a/examples/comb/in.comb.Si.elastic b/examples/comb/in.comb.Si.elastic index a0764b19b4..ca80b34892 100644 --- a/examples/comb/in.comb.Si.elastic +++ b/examples/comb/in.comb.Si.elastic @@ -39,7 +39,7 @@ run 1 ## calculates C11, C12, C13 and C14 fix 2 all deform 1 x scale 1.0001 remap x -compute perfx all stress/atom pair +compute perfx all stress/atom NULL pair compute fx all reduce sum & c_perfx[1] c_perfx[2] c_perfx[3] & c_perfx[4] c_perfx[5] c_perfx[6] @@ -50,7 +50,7 @@ run 10 ## strain z, measure s_z: calculates C33 fix 2 all deform 1 z scale 1.0001 remap x -compute perfz all stress/atom pair +compute perfz all stress/atom NULL pair compute fz all reduce sum & c_perfz[1] c_perfz[2] c_perfz[3] & c_perfz[4] c_perfz[5] c_perfz[6] @@ -61,7 +61,7 @@ run 10 ## strain yz, measure s_yz: calculates C44 fix 2 all deform 1 yz erate 0.0001 remap x -compute perfyz all stress/atom pair +compute perfyz all stress/atom NULL pair compute fyz all reduce sum & c_perfyz[1] c_perfyz[2] c_perfyz[3] & c_perfyz[4] c_perfyz[5] c_perfyz[6] @@ -72,7 +72,7 @@ run 10 ## strain xy, measure s_xy: calculates C66 fix 2 all deform 1 xy erate 0.0001 remap x -compute perfxy all stress/atom pair +compute perfxy all stress/atom NULL pair compute fxy all reduce sum & c_perfxy[1] c_perfxy[2] c_perfxy[3] & c_perfxy[4] c_perfxy[5] c_perfxy[6] From 09168613150e2492ecaa14e0a99279e25ac8f4f5 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 14:46:49 +0000 Subject: [PATCH 05/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11777 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/ASPHERE/fix_nve_asphere_noforce.cpp | 1 - src/ASPHERE/pair_resquared.cpp | 4 +- src/COLLOID/fix_wall_colloid.h | 2 +- src/DIPOLE/pair_lj_cut_dipole_long.cpp | 1 - src/FLD/pair_lubricateU_poly.cpp | 15 +--- src/GPU/pair_eam_alloy_gpu.cpp | 2 +- src/GPU/pair_eam_fs_gpu.cpp | 2 +- src/GPU/pair_yukawa_colloid_gpu.cpp | 3 +- src/GRANULAR/fix_pour.cpp | 3 +- src/KSPACE/commgrid.cpp | 4 - src/KSPACE/pppm_disp.cpp | 4 +- src/MANYBODY/fix_qeq_comb.cpp | 1 - src/MANYBODY/pair_adp.cpp | 4 +- src/MANYBODY/pair_eam.cpp | 4 +- src/MISC/compute_ti.cpp | 15 +--- src/MISC/fix_deposit.cpp | 1 - src/POEMS/fix_poems.cpp | 2 +- src/REPLICA/neb.cpp | 2 +- src/RIGID/fix_rigid_small.cpp | 8 +- src/USER-CUDA/cuda.cpp | 17 ++-- src/USER-CUDA/fft3d_cuda.cpp | 8 +- src/USER-CUDA/fix_addforce_cuda.cpp | 7 +- src/USER-CUDA/fix_aveforce_cuda.cpp | 7 +- src/USER-CUDA/fix_gravity_cuda.cpp | 23 +++--- src/USER-CUDA/fix_nh_cuda.cpp | 74 ++++++++--------- src/USER-CUDA/fix_set_force_cuda.cpp | 7 +- src/USER-CUDA/fix_shake_cuda.cpp | 14 ++-- src/USER-CUDA/fix_temp_berendsen_cuda.cpp | 6 +- src/USER-CUDA/fix_temp_rescale_cuda.cpp | 10 +-- src/USER-CUDA/fix_temp_rescale_limit_cuda.cpp | 12 +-- src/USER-CUDA/pair_lj_sdk_coul_long_cuda.h | 5 +- src/USER-CUDA/pair_lj_sdk_cuda.h | 5 +- src/USER-CUDA/pppm_cuda.cpp | 4 +- src/USER-CUDA/verlet_cuda.cpp | 9 +-- src/USER-EFF/compute_temp_eff.h | 1 - src/USER-EFF/fix_langevin_eff.cpp | 10 --- src/USER-EFF/fix_nh_eff.cpp | 1 - src/USER-EFF/fix_nvt_sllod_eff.cpp | 2 +- src/USER-MISC/angle_dipole.cpp | 2 +- src/USER-MISC/improper_cossq.cpp | 2 +- src/USER-MISC/pair_cdeam.cpp | 4 +- src/USER-MOLFILE/molfile_interface.cpp | 41 ---------- src/USER-OMP/angle_class2_omp.cpp | 1 + src/USER-OMP/angle_cosine_delta_omp.cpp | 1 + src/USER-OMP/angle_cosine_omp.cpp | 1 + src/USER-OMP/angle_cosine_periodic_omp.cpp | 1 + src/USER-OMP/angle_cosine_shift_exp_omp.cpp | 1 + src/USER-OMP/angle_dipole_omp.cpp | 1 + src/USER-OMP/angle_fourier_simple_omp.cpp | 1 + src/USER-OMP/angle_quartic_omp.cpp | 1 + src/USER-OMP/angle_sdk_omp.cpp | 1 + src/USER-OMP/angle_table_omp.cpp | 1 + src/USER-OMP/bond_class2_omp.cpp | 1 + src/USER-OMP/bond_fene_expand_omp.cpp | 1 + src/USER-OMP/bond_fene_omp.cpp | 1 + src/USER-OMP/bond_harmonic_shift_cut_omp.cpp | 1 + src/USER-OMP/bond_morse_omp.cpp | 1 + src/USER-OMP/bond_table_omp.cpp | 1 + src/USER-OMP/dihedral_fourier_omp.cpp | 8 +- src/USER-OMP/fix_shear_history_omp.cpp | 1 - src/USER-OMP/pair_airebo_omp.cpp | 1 - src/USER-OMP/pair_lj_cut_dipole_cut_omp.cpp | 2 +- .../pair_lj_gromacs_coul_gromacs_omp.cpp | 1 + src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp | 2 +- src/USER-OMP/pair_resquared_omp.cpp | 1 + src/USER-QMMM/README | 2 +- src/USER-REAXC/reaxc_allocate.cpp | 18 ++--- src/USER-REAXC/reaxc_io_tools.cpp | 8 +- src/atom_vec_sphere.h | 2 +- src/balance.cpp | 80 ++++++++++--------- src/balance.h | 2 + src/compute_atom_molecule.h | 2 +- src/compute_bond_local.h | 2 +- src/compute_pair_local.h | 2 +- src/compute_reduce.cpp | 4 +- src/displace_atoms.cpp | 2 +- src/dump_cfg.h | 2 +- src/fix.cpp | 1 + src/fix_ave_correlate.h | 1 - src/fix_ave_histo.h | 2 +- src/fix_ave_spatial.cpp | 8 +- src/fix_ave_time.h | 2 +- src/fix_balance.cpp | 5 +- src/fix_indent.h | 2 +- src/fix_vector.h | 2 +- src/fix_wall_region.h | 1 - src/group.cpp | 2 +- src/input.cpp | 16 ++-- src/lammps.cpp | 13 ++- src/lattice.cpp | 4 +- src/pair_lj_cut_coul_dsf.cpp | 2 +- src/random_mars.h | 2 +- 92 files changed, 252 insertions(+), 318 deletions(-) diff --git a/src/ASPHERE/fix_nve_asphere_noforce.cpp b/src/ASPHERE/fix_nve_asphere_noforce.cpp index 2b93fb5a7b..1cbc8fb16c 100644 --- a/src/ASPHERE/fix_nve_asphere_noforce.cpp +++ b/src/ASPHERE/fix_nve_asphere_noforce.cpp @@ -73,7 +73,6 @@ void FixNVEAsphereNoforce::initial_integrate(int vflag) double **angmom = atom->angmom; double *rmass = atom->rmass; int *ellipsoid = atom->ellipsoid; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp index 2eb0a52820..9764fe61dd 100755 --- a/src/ASPHERE/pair_resquared.cpp +++ b/src/ASPHERE/pair_resquared.cpp @@ -36,8 +36,8 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ PairRESquared::PairRESquared(LAMMPS *lmp) : Pair(lmp), - b_alpha(45.0/56.0), - cr60(pow(60.0,1.0/3.0)) + cr60(pow(60.0,1.0/3.0)), + b_alpha(45.0/56.0) { single_enable = 0; diff --git a/src/COLLOID/fix_wall_colloid.h b/src/COLLOID/fix_wall_colloid.h index 61aa4078f9..f047818a16 100644 --- a/src/COLLOID/fix_wall_colloid.h +++ b/src/COLLOID/fix_wall_colloid.h @@ -32,7 +32,7 @@ class FixWallColloid : public FixWall { void wall_particle(int, int, double); private: - double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6]; + double coeff1[6],coeff2[6],coeff3[6],coeff4[6]; }; } diff --git a/src/DIPOLE/pair_lj_cut_dipole_long.cpp b/src/DIPOLE/pair_lj_cut_dipole_long.cpp index 2e85b0976a..2796b54ac1 100755 --- a/src/DIPOLE/pair_lj_cut_dipole_long.cpp +++ b/src/DIPOLE/pair_lj_cut_dipole_long.cpp @@ -101,7 +101,6 @@ void PairLJCutDipoleLong::compute(int eflag, int vflag) double **torque = atom->torque; int *type = atom->type; int nlocal = atom->nlocal; - int nall = atom->nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; diff --git a/src/FLD/pair_lubricateU_poly.cpp b/src/FLD/pair_lubricateU_poly.cpp index 894c778789..545f1e3207 100644 --- a/src/FLD/pair_lubricateU_poly.cpp +++ b/src/FLD/pair_lubricateU_poly.cpp @@ -150,7 +150,7 @@ void PairLubricateUPoly::compute(int eflag, int vflag) void PairLubricateUPoly::iterate(double **x, int stage) { - int i,j,ii,itype; + int i,j,ii; int inum = list->inum; int *ilist = list->ilist; @@ -160,14 +160,12 @@ void PairLubricateUPoly::iterate(double **x, int stage) double alpha,beta; double normi,error,normig; double send[2],recv[2],rcg_dot_rcg; - double mo_inertia,radi; double **v = atom->v; double **f = atom->f; double **omega = atom->omega; double **angmom = atom->angmom; double **torque = atom->torque; - double *radius = atom->radius; // First compute R_FE*E @@ -312,8 +310,6 @@ void PairLubricateUPoly::iterate(double **x, int stage) for (ii=0;iiv; double **f = atom->f; double **omega = atom->omega; - double **angmom = atom->angmom; double **torque = atom->torque; double *radius = atom->radius; @@ -616,10 +610,8 @@ void PairLubricateUPoly::compute_RU(double **x) double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz; double rsq,r,radi,radj,h_sep; - //double beta0,beta1; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; double vt1,vt2,vt3,wdotn,wt1,wt2,wt3; - double inv_inertia; double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3],pre[2]; double **v = atom->v; @@ -1123,8 +1115,6 @@ void PairLubricateUPoly::compute_RE(double **x) void PairLubricateUPoly::settings(int narg, char **arg) { - int itype; - if (narg < 5 || narg > 7) error->all(FLERR,"Illegal pair_style command"); mu = force->numeric(FLERR,arg[0]); @@ -1166,8 +1156,6 @@ void PairLubricateUPoly::settings(int narg, char **arg) Ef[2][0] = 0.0; Ef[2][1] = 0.0; Ef[2][2] = 0.0; - - } /* ---------------------------------------------------------------------- @@ -1188,7 +1176,6 @@ void PairLubricateUPoly::init_style() // for pair hybrid, should limit test to types using the pair style double *radius = atom->radius; - int *type = atom->type; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) diff --git a/src/GPU/pair_eam_alloy_gpu.cpp b/src/GPU/pair_eam_alloy_gpu.cpp index 48b3982107..74b43f1103 100644 --- a/src/GPU/pair_eam_alloy_gpu.cpp +++ b/src/GPU/pair_eam_alloy_gpu.cpp @@ -151,7 +151,7 @@ void PairEAMAlloyGPU::read_file(char *filename) char **words = new char*[file->nelements+1]; nwords = 0; strtok(line," \t\n\r\f"); - while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; + while ( (words[nwords++] = strtok(NULL," \t\n\r\f")) ) continue; file->elements = new char*[file->nelements]; for (int i = 0; i < file->nelements; i++) { diff --git a/src/GPU/pair_eam_fs_gpu.cpp b/src/GPU/pair_eam_fs_gpu.cpp index f1600a9a2b..a6e80547ee 100644 --- a/src/GPU/pair_eam_fs_gpu.cpp +++ b/src/GPU/pair_eam_fs_gpu.cpp @@ -151,7 +151,7 @@ void PairEAMFSGPU::read_file(char *filename) char **words = new char*[file->nelements+1]; nwords = 0; strtok(line," \t\n\r\f"); - while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; + while ( (words[nwords++] = strtok(NULL," \t\n\r\f")) ) continue; file->elements = new char*[file->nelements]; for (int i = 0; i < file->nelements; i++) { diff --git a/src/GPU/pair_yukawa_colloid_gpu.cpp b/src/GPU/pair_yukawa_colloid_gpu.cpp index dfff5465ed..df19b4bc54 100644 --- a/src/GPU/pair_yukawa_colloid_gpu.cpp +++ b/src/GPU/pair_yukawa_colloid_gpu.cpp @@ -183,7 +183,7 @@ void PairYukawaColloidGPU::cpu_compute(int start, int inum, int eflag, int **firstneigh) { int i,j,ii,jj,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,radi,radj; - double r,rsq,r2inv,rinv,screening,forceyukawa,factor; + double r,rsq,rinv,screening,forceyukawa,factor; int *jlist; double **x = atom->x; @@ -217,7 +217,6 @@ void PairYukawaColloidGPU::cpu_compute(int start, int inum, int eflag, radj = radius[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; r = sqrt(rsq); rinv = 1.0/r; screening = exp(-kappa*(r-(radi+radj))); diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 69f95eb92e..abcd5d7335 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -209,7 +209,8 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : // volume_one = volume of inserted particle (with max possible radius) // in 3d, insure dy >= 1, for quasi-2d simulations - double volume,volume_one; + double volume,volume_one=0.0; + dstyle = -1; if (domain->dimension == 3) { if (region_style == 1) { double dy = yhi - ylo; diff --git a/src/KSPACE/commgrid.cpp b/src/KSPACE/commgrid.cpp index 993a5eabdf..7e2db0b148 100644 --- a/src/KSPACE/commgrid.cpp +++ b/src/KSPACE/commgrid.cpp @@ -488,8 +488,6 @@ void CommGrid::setup() void CommGrid::forward_comm(KSpace *kspace, int which) { - int i,n; - for (int m = 0; m < nswap; m++) { if (swap[m].sendproc == me) kspace->pack_forward(which,buf2,swap[m].npack,swap[m].packlist); @@ -515,8 +513,6 @@ void CommGrid::forward_comm(KSpace *kspace, int which) void CommGrid::reverse_comm(KSpace *kspace, int which) { - int i,n; - for (int m = nswap-1; m >= 0; m--) { if (swap[m].recvproc == me) kspace->pack_reverse(which,buf2,swap[m].nunpack,swap[m].unpacklist); diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index f49f5c8bbb..a3bb9adf22 100755 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -258,7 +258,7 @@ void PPPMDisp::init() memset(function, 0, EWALD_FUNCS*sizeof(int)); for (int i=0; i<=EWALD_MAXORDER; ++i) // transcribe order if (ewald_order&(1< 1.0e-4) { char str[128]; - sprintf(str,"Error in splitting of dispersion coeffs is estimated %g%",err); + sprintf(str,"Estimated error in splitting of dispersion coeffs is %g",err); error->warning(FLERR, str); } // set B diff --git a/src/MANYBODY/fix_qeq_comb.cpp b/src/MANYBODY/fix_qeq_comb.cpp index 7a37dec974..c40147f467 100644 --- a/src/MANYBODY/fix_qeq_comb.cpp +++ b/src/MANYBODY/fix_qeq_comb.cpp @@ -206,7 +206,6 @@ void FixQEQComb::post_force(int vflag) double *q = atom->q; int *mask = atom->mask; - int nlocal = atom->nlocal; if (comb) { inum = comb->list->inum; diff --git a/src/MANYBODY/pair_adp.cpp b/src/MANYBODY/pair_adp.cpp index 0c06cc06eb..5a60a5f5dd 100644 --- a/src/MANYBODY/pair_adp.cpp +++ b/src/MANYBODY/pair_adp.cpp @@ -577,7 +577,7 @@ void PairADP::read_file(char *filename) char **words = new char*[file->nelements+1]; nwords = 0; strtok(line," \t\n\r\f"); - while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; + while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue; file->elements = new char*[file->nelements]; for (int i = 0; i < file->nelements; i++) { @@ -927,7 +927,7 @@ void PairADP::grab(FILE *fp, int n, double *list) fgets(line,MAXLINE,fp); ptr = strtok(line," \t\n\r\f"); list[i++] = atof(ptr); - while (ptr = strtok(NULL," \t\n\r\f")) list[i++] = atof(ptr); + while ((ptr = strtok(NULL," \t\n\r\f"))) list[i++] = atof(ptr); } } diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index 9a86a14c12..e38c2caa66 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -149,6 +149,7 @@ void PairEAM::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; int newton_pair = force->newton_pair; inum = list->inum; @@ -159,8 +160,7 @@ void PairEAM::compute(int eflag, int vflag) // zero out density if (newton_pair) { - m = nlocal + atom->nghost; - for (i = 0; i < m; i++) rho[i] = 0.0; + for (i = 0; i < nall; i++) rho[i] = 0.0; } else for (i = 0; i < nlocal; i++) rho[i] = 0.0; // rho = density at each atom diff --git a/src/MISC/compute_ti.cpp b/src/MISC/compute_ti.cpp index 90083b0666..f02277d6c6 100644 --- a/src/MISC/compute_ti.cpp +++ b/src/MISC/compute_ti.cpp @@ -160,6 +160,8 @@ double ComputeTI::compute_scalar() if (update->eflag_global != invoked_scalar) error->all(FLERR,"Energy was not tallied on needed timestep"); + const int nlocal = atom->nlocal; + const int * const type = atom->type; double dUdl = 0.0; for (int m = 0; m < nterms; m++) { @@ -171,18 +173,12 @@ double ComputeTI::compute_scalar() if (value1 == 0.0) continue; if (which[m] == PAIR) { - int ntypes = atom->ntypes; - int *mask = atom->mask; if (total_flag) { eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul; MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); } else { - int nlocal = atom->nlocal; int npair = nlocal; - int *mask = atom->mask; - int *type = atom->type; - double *eatom = pptr[m]->eatom; if (force->newton_pair) npair += atom->nghost; @@ -215,17 +211,10 @@ double ComputeTI::compute_scalar() dUdl += eng/value1 * value2; } else if (which[m] == KSPACE) { - int ntypes = atom->ntypes; - int *mask = atom->mask; if (total_flag) eng = force->kspace->energy; else { - int nlocal = atom->nlocal; - int npair = nlocal; - int *mask = atom->mask; - int *type = atom->type; double *eatom = force->kspace->eatom; - eng = 0; for(int i = 0; i < nlocal; i++) if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) eng += eatom[i]; diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp index 196b3b489b..ef146ea11b 100644 --- a/src/MISC/fix_deposit.cpp +++ b/src/MISC/fix_deposit.cpp @@ -248,7 +248,6 @@ void FixDeposit::pre_exchange() { int i,j,m,n,nlocalprev,flag,flagall; double coord[3],lamda[3],delx,dely,delz,rsq; - double alpha,beta,gamma; double r[3],vnew[3],rotmat[3][3],quat[4]; double *newcoord; diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp index e1cd2bae31..7f07accc57 100644 --- a/src/POEMS/fix_poems.cpp +++ b/src/POEMS/fix_poems.cpp @@ -926,7 +926,7 @@ void FixPOEMS::readfile(char *file) if (ptr == NULL || ptr[0] == '#') continue; ptr = strtok(NULL," ,\t\n\0"); - while (ptr = strtok(NULL," ,\t\n\0")) { + while ((ptr = strtok(NULL," ,\t\n\0"))) { id = atoi(ptr); i = atom->map(id); if (i < 0 || i >= nlocal) continue; diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index 5fd168a943..4a6649d1b6 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -412,7 +412,7 @@ void NEB::readfile(char *file, int flag) // loop over lines of atom coords // tokenize the line into values - for (int i = 0; i < nchunk; i++) { + for (i = 0; i < nchunk; i++) { next = strchr(buf,'\n'); values[0] = strtok(buf," \t\n\r\f"); diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index c4f3ccaab4..d82f17d0c0 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -384,7 +384,6 @@ void FixRigidSmall::setup_pre_neighbor() void FixRigidSmall::setup(int vflag) { int i,n,ibody; - double massone,radone; //check(1); @@ -394,7 +393,6 @@ void FixRigidSmall::setup(int vflag) double **x = atom->x; double **f = atom->f; - int *type = atom->type; imageint *image = atom->image; int nlocal = atom->nlocal; @@ -2033,7 +2031,7 @@ void FixRigidSmall::setup_bodies_static() void FixRigidSmall::setup_bodies_dynamic() { - int i,n,ibody; + int i,ibody; double massone,radone; // sum vcm, angmom across all rigid bodies @@ -2155,11 +2153,10 @@ void FixRigidSmall::readfile(int which, double **array, int *inbody) // key = mol ID of bodies my atoms own // value = index into local body array - tagint *molecule = atom->molecule; int nlocal = atom->nlocal; hash = new std::map(); - for (int i = 0; i < nlocal; i++) + for (i = 0; i < nlocal; i++) if (bodyown[i] >= 0) (*hash)[atom->molecule[i]] = bodyown[i]; // open file and read header @@ -2448,7 +2445,6 @@ void FixRigidSmall::set_molecule(int nlocalprev, tagint tagprev, int nlocal = atom->nlocal; if (nlocalprev == nlocal) return; - double **x = atom->x; tagint *tag = atom->tag; for (int i = nlocalprev; i < nlocal; i++) { diff --git a/src/USER-CUDA/cuda.cpp b/src/USER-CUDA/cuda.cpp index 1160c8cbf2..1d678144a8 100644 --- a/src/USER-CUDA/cuda.cpp +++ b/src/USER-CUDA/cuda.cpp @@ -131,7 +131,6 @@ Cuda::Cuda(LAMMPS* lmp) : Pointers(lmp) pinned = true; debugdata = 0; - new int[2 * CUDA_MAX_DEBUG_SIZE]; finished_setup = false; begin_setup = false; @@ -173,6 +172,7 @@ Cuda::~Cuda() delete cu_virial; delete cu_eng_vdwl; delete cu_eng_coul; + delete cu_extent; delete cu_eatom; delete cu_vatom; delete cu_radius; @@ -188,6 +188,9 @@ Cuda::~Cuda() delete cu_omega_rmass; delete [] omega_rmass; + delete cu_debugdata; + delete[] debugdata; + delete cu_map_array; std::map::iterator p = neigh_lists.begin(); @@ -213,14 +216,14 @@ void Cuda::accelerator(int narg, char** arg) if(++i == narg) error->all(FLERR, "Invalid Options for 'accelerator' command. Expecting a number after 'gpu/node' option."); - pppn = atoi(arg[i]); + pppn = force->inumeric(FLERR,arg[i]); } if(strcmp(arg[i], "gpu/node/special") == 0) { if(++i == narg) error->all(FLERR, "Invalid Options for 'accelerator' command. Expecting number of GPUs to be used per node after keyword 'gpu/node/special'."); - pppn = atoi(arg[i]); + pppn = force->inumeric(FLERR,arg[i]); if(pppn < 1) error->all(FLERR, "Invalid Options for 'accelerator' command. Expecting number of GPUs to be used per node after keyword 'gpu/node special'."); @@ -231,7 +234,7 @@ void Cuda::accelerator(int narg, char** arg) for(int k = 0; k < pppn; k++) { i++; - devicelist[k] = atoi(arg[i]); + devicelist[k] = force->inumeric(FLERR,arg[i]); } } @@ -239,7 +242,7 @@ void Cuda::accelerator(int narg, char** arg) if(++i == narg) error->all(FLERR, "Invalid Options for 'accelerator' command. Expecting a number after 'pinned' option."); - pinned = atoi(arg[i]) == 0 ? false : true; + pinned = force->inumeric(FLERR,arg[i]) == 0 ? false : true; if((pinned == false) && (universe->me == 0)) printf(" #CUDA: Pinned memory is not used for communication\n"); } @@ -263,7 +266,7 @@ void Cuda::accelerator(int narg, char** arg) if(++i == narg) error->all(FLERR, "Invalid Options for 'accelerator' command. Expecting a number after 'test' option."); - testatom = atof(arg[i]); + testatom = force->numeric(FLERR,arg[i]); dotestatom = true; } @@ -271,7 +274,7 @@ void Cuda::accelerator(int narg, char** arg) if(++i == narg) error->all(FLERR, "Invalid Options for 'accelerator' command. Expecting a number after 'override/bpa' option."); - shared_data.pair.override_block_per_atom = atoi(arg[i]); + shared_data.pair.override_block_per_atom = force->inumeric(FLERR,arg[i]); } } diff --git a/src/USER-CUDA/fft3d_cuda.cpp b/src/USER-CUDA/fft3d_cuda.cpp index bd1116e447..da52662bfb 100644 --- a/src/USER-CUDA/fft3d_cuda.cpp +++ b/src/USER-CUDA/fft3d_cuda.cpp @@ -243,7 +243,7 @@ struct fft_plan_3d *fft_3d_create_plan_cuda( remap_3d_create_plan(comm,in_ilo,in_ihi,in_jlo,in_jhi,in_klo,in_khi, first_ilo,first_ihi,first_jlo,first_jhi, first_klo,first_khi, - members,0,0,2); + members,0,0,2,0); if (plan->pre_plan == NULL) return NULL; } @@ -268,7 +268,7 @@ struct fft_plan_3d *fft_3d_create_plan_cuda( first_klo,first_khi, second_ilo,second_ihi,second_jlo,second_jhi, second_klo,second_khi, - 2,1,0,2); + 2,1,0,2,0); if (plan->mid1_plan == NULL) return NULL; // 1d FFTs along mid axis @@ -308,7 +308,7 @@ struct fft_plan_3d *fft_3d_create_plan_cuda( second_ilo,second_ihi, third_jlo,third_jhi,third_klo,third_khi, third_ilo,third_ihi, - 2,1,0,2); + 2,1,0,2,0); if (plan->mid2_plan == NULL) return NULL; // 1d FFTs along slow axis @@ -332,7 +332,7 @@ struct fft_plan_3d *fft_3d_create_plan_cuda( third_jlo,third_jhi, out_klo,out_khi,out_ilo,out_ihi, out_jlo,out_jhi, - 2,(permute+1)%3,0,2); + 2,(permute+1)%3,0,2,0); if (plan->post_plan == NULL) return NULL; } diff --git a/src/USER-CUDA/fix_addforce_cuda.cpp b/src/USER-CUDA/fix_addforce_cuda.cpp index 03519b80a9..9503e1d828 100644 --- a/src/USER-CUDA/fix_addforce_cuda.cpp +++ b/src/USER-CUDA/fix_addforce_cuda.cpp @@ -30,6 +30,7 @@ #include "update.h" #include "respa.h" #include "error.h" +#include "force.h" #include "domain.h" #include "cuda.h" #include "memory.h" @@ -58,9 +59,9 @@ FixAddForceCuda::FixAddForceCuda(LAMMPS *lmp, int narg, char **arg) : extscalar = 1; extvector = 1; - xvalue = atof(arg[3]); - yvalue = atof(arg[4]); - zvalue = atof(arg[5]); + xvalue = force->numeric(FLERR,arg[3]); + yvalue = force->numeric(FLERR,arg[4]); + zvalue = force->numeric(FLERR,arg[5]); // optional args diff --git a/src/USER-CUDA/fix_aveforce_cuda.cpp b/src/USER-CUDA/fix_aveforce_cuda.cpp index 8c332ee8a8..b6b9cb705a 100644 --- a/src/USER-CUDA/fix_aveforce_cuda.cpp +++ b/src/USER-CUDA/fix_aveforce_cuda.cpp @@ -38,6 +38,7 @@ #include "modify.h" #include "atom_masks.h" #include "error.h" +#include "force.h" using namespace LAMMPS_NS; @@ -72,7 +73,7 @@ FixAveForceCuda::FixAveForceCuda(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[3],"NULL") == 0) { xstyle = NONE; } else { - xvalue = atof(arg[3]); + xvalue = force->numeric(FLERR,arg[3]); xstyle = CONSTANT; } if (strstr(arg[4],"v_") == arg[4]) { @@ -82,7 +83,7 @@ FixAveForceCuda::FixAveForceCuda(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[4],"NULL") == 0) { ystyle = NONE; } else { - yvalue = atof(arg[4]); + yvalue = force->numeric(FLERR,arg[4]); ystyle = CONSTANT; } if (strstr(arg[5],"v_") == arg[5]) { @@ -92,7 +93,7 @@ FixAveForceCuda::FixAveForceCuda(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[5],"NULL") == 0) { zstyle = NONE; } else { - zvalue = atof(arg[5]); + zvalue = force->numeric(FLERR,arg[5]); zstyle = CONSTANT; } diff --git a/src/USER-CUDA/fix_gravity_cuda.cpp b/src/USER-CUDA/fix_gravity_cuda.cpp index 549f95edbf..a7234512d3 100644 --- a/src/USER-CUDA/fix_gravity_cuda.cpp +++ b/src/USER-CUDA/fix_gravity_cuda.cpp @@ -34,6 +34,7 @@ #include "cuda_modify_flags.h" #include "math_const.h" #include "error.h" +#include "force.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -53,31 +54,31 @@ FixGravityCuda::FixGravityCuda(LAMMPS *lmp, int narg, char **arg) : if (narg < 5) error->all(FLERR,"Illegal fix gravity command"); - magnitude = atof(arg[3]); + magnitude = force->numeric(FLERR,arg[3]); if (strcmp(arg[4],"chute") == 0) { if (narg != 6) error->all(FLERR,"Illegal fix gravity command"); style = CHUTE; phi = 0.0; - theta = 180.0 - atof(arg[5]); + theta = 180.0 - force->numeric(FLERR,arg[5]); } else if (strcmp(arg[4],"spherical") == 0) { if (narg != 7) error->all(FLERR,"Illegal fix gravity command"); style = SPHERICAL; - phi = atof(arg[5]); - theta = atof(arg[6]); + phi = force->numeric(FLERR,arg[5]); + theta = force->numeric(FLERR,arg[6]); } else if (strcmp(arg[4],"gradient") == 0) { if (narg != 9) error->all(FLERR,"Illegal fix gravity command"); style = GRADIENT; - phi = atof(arg[5]); - theta = atof(arg[6]); - phigrad = atof(arg[7]); - thetagrad = atof(arg[8]); + phi = force->numeric(FLERR,arg[5]); + theta = force->numeric(FLERR,arg[6]); + phigrad = force->numeric(FLERR,arg[7]); + thetagrad = force->numeric(FLERR,arg[8]); } else if (strcmp(arg[4],"vector") == 0) { if (narg != 8) error->all(FLERR,"Illegal fix gravity command"); style = VECTOR; - xdir = atof(arg[5]); - ydir = atof(arg[6]); - zdir = atof(arg[7]); + xdir = force->numeric(FLERR,arg[5]); + ydir = force->numeric(FLERR,arg[6]); + zdir = force->numeric(FLERR,arg[7]); } else error->all(FLERR,"Illegal fix gravity command"); degree2rad = MY_PI/180.0; diff --git a/src/USER-CUDA/fix_nh_cuda.cpp b/src/USER-CUDA/fix_nh_cuda.cpp index 43d06cc655..c0beac2182 100644 --- a/src/USER-CUDA/fix_nh_cuda.cpp +++ b/src/USER-CUDA/fix_nh_cuda.cpp @@ -99,9 +99,9 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (strcmp(arg[iarg],"temp") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); tstat_flag = 1; - t_start = atof(arg[iarg+1]); - t_stop = atof(arg[iarg+2]); - t_period = atof(arg[iarg+3]); + t_start = force->numeric(FLERR,arg[iarg+1]); + t_stop = force->numeric(FLERR,arg[iarg+2]); + t_period = force->numeric(FLERR,arg[iarg+3]); if (t_start < 0.0 || t_stop <= 0.0) error->all(FLERR,"Target T for fix nvt/npt/nph cannot be 0.0"); iarg += 4; @@ -109,9 +109,9 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg],"iso") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = XYZ; - p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); - p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); - p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = p_flag[1] = p_flag[2] = 1; if (dimension == 2) { p_start[2] = p_stop[2] = p_period[2] = 0.0; @@ -121,9 +121,9 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg],"aniso") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = NONE; - p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); - p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); - p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = p_flag[1] = p_flag[2] = 1; if (dimension == 2) { p_start[2] = p_stop[2] = p_period[2] = 0.0; @@ -133,13 +133,13 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg],"tri") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = NONE; - p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); - p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); - p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = p_flag[1] = p_flag[2] = 1; p_start[3] = p_start[4] = p_start[5] = 0.0; p_stop[3] = p_stop[4] = p_stop[5] = 0.0; - p_period[3] = p_period[4] = p_period[5] = atof(arg[iarg+3]); + p_period[3] = p_period[4] = p_period[5] = force->numeric(FLERR,arg[iarg+3]); p_flag[3] = p_flag[4] = p_flag[5] = 1; if (dimension == 2) { p_start[2] = p_stop[2] = p_period[2] = 0.0; @@ -153,25 +153,25 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg],"x") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - p_start[0] = atof(arg[iarg+1]); - p_stop[0] = atof(arg[iarg+2]); - p_period[0] = atof(arg[iarg+3]); + p_start[0] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = force->numeric(FLERR,arg[iarg+3]); p_flag[0] = 1; deviatoric_flag = 1; iarg += 4; } else if (strcmp(arg[iarg],"y") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - p_start[1] = atof(arg[iarg+1]); - p_stop[1] = atof(arg[iarg+2]); - p_period[1] = atof(arg[iarg+3]); + p_start[1] = force->numeric(FLERR,arg[iarg+1]); + p_stop[1] = force->numeric(FLERR,arg[iarg+2]); + p_period[1] = force->numeric(FLERR,arg[iarg+3]); p_flag[1] = 1; deviatoric_flag = 1; iarg += 4; } else if (strcmp(arg[iarg],"z") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - p_start[2] = atof(arg[iarg+1]); - p_stop[2] = atof(arg[iarg+2]); - p_period[2] = atof(arg[iarg+3]); + p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[2] = 1; deviatoric_flag = 1; iarg += 4; @@ -180,9 +180,9 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg],"yz") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - p_start[3] = atof(arg[iarg+1]); - p_stop[3] = atof(arg[iarg+2]); - p_period[3] = atof(arg[iarg+3]); + p_start[3] = force->numeric(FLERR,arg[iarg+1]); + p_stop[3] = force->numeric(FLERR,arg[iarg+2]); + p_period[3] = force->numeric(FLERR,arg[iarg+3]); p_flag[3] = 1; deviatoric_flag = 1; iarg += 4; @@ -190,9 +190,9 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); } else if (strcmp(arg[iarg],"xz") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - p_start[4] = atof(arg[iarg+1]); - p_stop[4] = atof(arg[iarg+2]); - p_period[4] = atof(arg[iarg+3]); + p_start[4] = force->numeric(FLERR,arg[iarg+1]); + p_stop[4] = force->numeric(FLERR,arg[iarg+2]); + p_period[4] = force->numeric(FLERR,arg[iarg+3]); p_flag[4] = 1; deviatoric_flag = 1; iarg += 4; @@ -200,9 +200,9 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); } else if (strcmp(arg[iarg],"xy") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - p_start[5] = atof(arg[iarg+1]); - p_stop[5] = atof(arg[iarg+2]); - p_period[5] = atof(arg[iarg+3]); + p_start[5] = force->numeric(FLERR,arg[iarg+1]); + p_stop[5] = force->numeric(FLERR,arg[iarg+2]); + p_period[5] = force->numeric(FLERR,arg[iarg+3]); p_flag[5] = 1; deviatoric_flag = 1; iarg += 4; @@ -219,7 +219,7 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg],"drag") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - drag = atof(arg[iarg+1]); + drag = force->numeric(FLERR,arg[iarg+1]); if (drag < 0.0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"dilate") == 0) { @@ -230,12 +230,12 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) iarg += 2; } else if (strcmp(arg[iarg],"tchain") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - mtchain = atoi(arg[iarg+1]); + mtchain = force->inumeric(FLERR,arg[iarg+1]); if (mtchain < 1) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"pchain") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - mpchain = atoi(arg[iarg+1]); + mpchain = force->inumeric(FLERR,arg[iarg+1]); if (mpchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"mtk") == 0) { @@ -246,17 +246,17 @@ FixNHCuda::FixNHCuda(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) iarg += 2; } else if (strcmp(arg[iarg],"tloop") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - nc_tchain = atoi(arg[iarg+1]); + nc_tchain = force->inumeric(FLERR,arg[iarg+1]); if (nc_tchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"ploop") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - nc_pchain = atoi(arg[iarg+1]); + nc_pchain = force->inumeric(FLERR,arg[iarg+1]); if (nc_pchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"nreset") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); - nreset_h0 = atoi(arg[iarg+1]); + nreset_h0 = force->inumeric(FLERR,arg[iarg+1]); if (nreset_h0 < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else error->all(FLERR,"Illegal fix nvt/npt/nph command"); diff --git a/src/USER-CUDA/fix_set_force_cuda.cpp b/src/USER-CUDA/fix_set_force_cuda.cpp index 1010c087ad..afaf5c7e88 100644 --- a/src/USER-CUDA/fix_set_force_cuda.cpp +++ b/src/USER-CUDA/fix_set_force_cuda.cpp @@ -28,6 +28,7 @@ #include "update.h" #include "respa.h" #include "error.h" +#include "force.h" #include "cuda.h" #include "memory.h" #include "cuda_modify_flags.h" @@ -55,11 +56,11 @@ FixSetForceCuda::FixSetForceCuda(LAMMPS *lmp, int narg, char **arg) : flagx = flagy = flagz = 1; if (strcmp(arg[3],"NULL") == 0) flagx = 0; - else xvalue = atof(arg[3]); + else xvalue = force->numeric(FLERR,arg[3]); if (strcmp(arg[4],"NULL") == 0) flagy = 0; - else yvalue = atof(arg[4]); + else yvalue = force->numeric(FLERR,arg[4]); if (strcmp(arg[5],"NULL") == 0) flagz = 0; - else zvalue = atof(arg[5]); + else zvalue = force->numeric(FLERR,arg[5]); force_flag = 0; foriginal[0] = foriginal[1] = foriginal[2] = 0.0; diff --git a/src/USER-CUDA/fix_shake_cuda.cpp b/src/USER-CUDA/fix_shake_cuda.cpp index 8391626615..6a2ea8d6e4 100644 --- a/src/USER-CUDA/fix_shake_cuda.cpp +++ b/src/USER-CUDA/fix_shake_cuda.cpp @@ -95,9 +95,9 @@ FixShakeCuda::FixShakeCuda(LAMMPS* lmp, int narg, char** arg) : if(narg < 8) error->all(FLERR, "Illegal fix shake command"); - tolerance = atof(arg[3]); - max_iter = atoi(arg[4]); - output_every = atoi(arg[5]); + tolerance = force->numeric(FLERR,arg[3]); + max_iter = force->inumeric(FLERR,arg[4]); + output_every = force->inumeric(FLERR,arg[5]); // parse SHAKE args for bond and angle types // will be used by find_clusters @@ -133,7 +133,7 @@ FixShakeCuda::FixShakeCuda(LAMMPS* lmp, int narg, char** arg) : atom->check_mass(); } else if(mode == 'b') { - int i = atoi(arg[next]); + int i = force->inumeric(FLERR,arg[next]); if(i < 1 || i > atom->nbondtypes) error->all(FLERR, "Invalid bond type index for fix shake"); @@ -141,7 +141,7 @@ FixShakeCuda::FixShakeCuda(LAMMPS* lmp, int narg, char** arg) : bond_flag[i] = 1; } else if(mode == 'a') { - int i = atoi(arg[next]); + int i = force->inumeric(FLERR,arg[next]); if(i < 1 || i > atom->nangletypes) error->all(FLERR, "Invalid angle type index for fix shake"); @@ -149,7 +149,7 @@ FixShakeCuda::FixShakeCuda(LAMMPS* lmp, int narg, char** arg) : angle_flag[i] = 1; } else if(mode == 't') { - int i = atoi(arg[next]); + int i = force->inumeric(FLERR,arg[next]); if(i < 1 || i > atom->ntypes) error->all(FLERR, "Invalid atom type index for fix shake"); @@ -157,7 +157,7 @@ FixShakeCuda::FixShakeCuda(LAMMPS* lmp, int narg, char** arg) : type_flag[i] = 1; } else if(mode == 'm') { - double massone = atof(arg[next]); + double massone = force->numeric(FLERR,arg[next]); if(massone == 0.0) error->all(FLERR, "Invalid atom mass for fix shake"); diff --git a/src/USER-CUDA/fix_temp_berendsen_cuda.cpp b/src/USER-CUDA/fix_temp_berendsen_cuda.cpp index 5f5e13dd3c..1e29de8db2 100644 --- a/src/USER-CUDA/fix_temp_berendsen_cuda.cpp +++ b/src/USER-CUDA/fix_temp_berendsen_cuda.cpp @@ -71,9 +71,9 @@ FixTempBerendsenCuda::FixTempBerendsenCuda(LAMMPS *lmp, int narg, char **arg) : nevery = 1; - t_start = atof(arg[3]); - t_stop = atof(arg[4]); - t_period = atof(arg[5]); + t_start = force->numeric(FLERR,arg[3]); + t_stop = force->numeric(FLERR,arg[4]); + t_period = force->numeric(FLERR,arg[5]); // error checks diff --git a/src/USER-CUDA/fix_temp_rescale_cuda.cpp b/src/USER-CUDA/fix_temp_rescale_cuda.cpp index 85684d9b46..36c75d446b 100644 --- a/src/USER-CUDA/fix_temp_rescale_cuda.cpp +++ b/src/USER-CUDA/fix_temp_rescale_cuda.cpp @@ -56,17 +56,17 @@ FixTempRescaleCuda::FixTempRescaleCuda(LAMMPS *lmp, int narg, char **arg) : if (narg < 8) error->all(FLERR,"Illegal fix temp/rescale/cuda command"); - nevery = atoi(arg[3]); + nevery = force->inumeric(FLERR,arg[3]); if (nevery <= 0) error->all(FLERR,"Illegal fix temp/rescale/cuda command"); scalar_flag = 1; global_freq = nevery; extscalar = 1; - t_start = atof(arg[4]); - t_stop = atof(arg[5]); - t_window = atof(arg[6]); - fraction = atof(arg[7]); + t_start = force->numeric(FLERR,arg[4]); + t_stop = force->numeric(FLERR,arg[5]); + t_window = force->numeric(FLERR,arg[6]); + fraction = force->numeric(FLERR,arg[7]); // create a new compute temp // id = fix-ID + temp, compute group = fix group diff --git a/src/USER-CUDA/fix_temp_rescale_limit_cuda.cpp b/src/USER-CUDA/fix_temp_rescale_limit_cuda.cpp index a7d4aaa4d7..3675ca0160 100644 --- a/src/USER-CUDA/fix_temp_rescale_limit_cuda.cpp +++ b/src/USER-CUDA/fix_temp_rescale_limit_cuda.cpp @@ -56,18 +56,18 @@ FixTempRescaleLimitCuda::FixTempRescaleLimitCuda(LAMMPS *lmp, int narg, char **a if (narg < 9) error->all(FLERR,"Illegal fix temp/rescale/limit/cuda command"); - nevery = atoi(arg[3]); + nevery = force->inumeric(FLERR,arg[3]); if (nevery <= 0) error->all(FLERR,"Illegal fix temp/rescale/limit/cuda command"); scalar_flag = 1; global_freq = nevery; extscalar = 1; - t_start = atof(arg[4]); - t_stop = atof(arg[5]); - t_window = atof(arg[6]); - fraction = atof(arg[7]); - limit = atof(arg[8]); + t_start = force->numeric(FLERR,arg[4]); + t_stop = force->numeric(FLERR,arg[5]); + t_window = force->numeric(FLERR,arg[6]); + fraction = force->numeric(FLERR,arg[7]); + limit = force->numeric(FLERR,arg[8]); if (limit <= 1.0) error->all(FLERR,"Illegal fix temp/rescale/limit/cuda command (limit must be > 1.0)"); diff --git a/src/USER-CUDA/pair_lj_sdk_coul_long_cuda.h b/src/USER-CUDA/pair_lj_sdk_coul_long_cuda.h index 8debb81666..78bc504d4f 100644 --- a/src/USER-CUDA/pair_lj_sdk_coul_long_cuda.h +++ b/src/USER-CUDA/pair_lj_sdk_coul_long_cuda.h @@ -23,12 +23,13 @@ #ifdef PAIR_CLASS +PairStyle(cg/cmm/coul/long/cuda,PairLJSDKCoulLongCuda) PairStyle(lj/sdk/coul/long/cuda,PairLJSDKCoulLongCuda) #else -#ifndef PAIR_CG_CMM_COUL_LONG_CUDA_H -#define PAIR_CG_CMM_COUL_LONG_CUDA_H +#ifndef PAIR_LJ_SDK_COUL_LONG_CUDA_H +#define PAIR_LJ_SDK_COUL_LONG_CUDA_H #include "pair_lj_sdk_coul_long.h" diff --git a/src/USER-CUDA/pair_lj_sdk_cuda.h b/src/USER-CUDA/pair_lj_sdk_cuda.h index f5737adadc..8afd6c0a73 100644 --- a/src/USER-CUDA/pair_lj_sdk_cuda.h +++ b/src/USER-CUDA/pair_lj_sdk_cuda.h @@ -24,11 +24,12 @@ #ifdef PAIR_CLASS PairStyle(lj/sdk/cuda,PairLJSDKCuda) +PairStyle(cg/cmm/cuda,PairLJSDKCuda) #else -#ifndef PAIR_CG_CMM_CUDA_H -#define PAIR_CG_CMM_CUDA_H +#ifndef PAIR_LJ_SDK_CUDA_H +#define PAIR_LJ_SDK_CUDA_H #include "pair_lj_sdk.h" #include "cuda_data.h" diff --git a/src/USER-CUDA/pppm_cuda.cpp b/src/USER-CUDA/pppm_cuda.cpp index f07310aa52..9754dd1f94 100644 --- a/src/USER-CUDA/pppm_cuda.cpp +++ b/src/USER-CUDA/pppm_cuda.cpp @@ -110,7 +110,7 @@ PPPMCuda::PPPMCuda(LAMMPS *lmp, int narg, char **arg) : #endif triclinic_support = 0; - accuracy_relative = atof(arg[0]); + accuracy_relative = fabs(force->numeric(FLERR,arg[0])); nfactors = 3; factors = new int[nfactors]; @@ -934,7 +934,7 @@ void PPPMCuda::allocate() remap = new Remap(lmp,world, nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in, nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft, - 1,0,0,2); + 1,0,0,2,0); pppm_device_init(cu_density_brick->dev_data(), cu_vdx_brick->dev_data(), cu_vdy_brick->dev_data(), cu_vdz_brick->dev_data(), cu_density_fft->dev_data(),cu_energy->dev_data(),cu_virial->dev_data() diff --git a/src/USER-CUDA/verlet_cuda.cpp b/src/USER-CUDA/verlet_cuda.cpp index 975b1349a4..d5ba68b7ec 100644 --- a/src/USER-CUDA/verlet_cuda.cpp +++ b/src/USER-CUDA/verlet_cuda.cpp @@ -116,10 +116,7 @@ void VerletCuda::setup() if(cuda->shared_data.me == 0) - printf("# CUDA: VerletCuda::setup: Allocate memory on device for maximum of %i atoms...\n", atom->nmax); - - if(cuda->shared_data.me == 0) - printf("# CUDA: Using precision: Global: %u X: %u V: %u F: %u PPPM: %u \n", CUDA_PRECISION == 1 ? 4 : 8, sizeof(X_FLOAT), sizeof(V_FLOAT), sizeof(F_FLOAT), sizeof(PPPM_FLOAT)); + printf("# CUDA: Using precision: Global: %u X: %u V: %u F: %u PPPM: %u \n", CUDA_PRECISION == 1 ? 4 : 8, (int) sizeof(X_FLOAT), (int) sizeof(V_FLOAT), (int) sizeof(F_FLOAT), (int) sizeof(PPPM_FLOAT)); cuda->allocate(); @@ -375,7 +372,7 @@ void VerletCuda::setup() test_atom(testatom, "post reverse comm"); if(cuda->shared_data.me == 0) - printf("# CUDA: Total Device Memory useage post setup: %lf MB\n", 1.0 * CudaWrapper_CheckMemUseage() / 1024 / 1024); + printf("# CUDA: Total Device Memory usage post setup: %lf MB\n", 1.0 * CudaWrapper_CheckMemUsage() / 1024 / 1024); MYDBG(printf("# CUDA: VerletCuda::setup: call modify setup\n");) modify->setup(vflag); @@ -565,7 +562,7 @@ void VerletCuda::setup_minimal(int flag) test_atom(testatom, "post reverse comm"); if(cuda->shared_data.me == 0) - printf("# CUDA: Total Device Memory useage post setup: %lf MB\n", 1.0 * CudaWrapper_CheckMemUseage() / 1024 / 1024); + printf("# CUDA: Total Device Memory usage post setup: %lf MB\n", 1.0 * CudaWrapper_CheckMemUsage() / 1024 / 1024); MYDBG(printf("# CUDA: VerletCuda::setup: call modify setup\n");) modify->setup(vflag); diff --git a/src/USER-EFF/compute_temp_eff.h b/src/USER-EFF/compute_temp_eff.h index bd41db6299..7996619f1a 100644 --- a/src/USER-EFF/compute_temp_eff.h +++ b/src/USER-EFF/compute_temp_eff.h @@ -36,7 +36,6 @@ class ComputeTempEff : public Compute { private: int fix_dof; double tfactor; - double *inertia; void dof_compute(); }; diff --git a/src/USER-EFF/fix_langevin_eff.cpp b/src/USER-EFF/fix_langevin_eff.cpp index 168b3b2ba1..826c606818 100644 --- a/src/USER-EFF/fix_langevin_eff.cpp +++ b/src/USER-EFF/fix_langevin_eff.cpp @@ -132,11 +132,6 @@ void FixLangevinEff::post_force_no_tally() double fran[4],fsum[4],fsumall[4]; fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0; - double boltz = force->boltz; - double dt = update->dt; - double mvv2e = force->mvv2e; - double ftm2v = force->ftm2v; - int particles = group->count(igroup); if (zeroflag) { if (particles == 0) @@ -306,11 +301,6 @@ void FixLangevinEff::post_force_tally() // test v = 0 since some computes mask non-participating atoms via v = 0 // and added force has extra term not multiplied by v = 0 - double boltz = force->boltz; - double dt = update->dt; - double mvv2e = force->mvv2e; - double ftm2v = force->ftm2v; - int particles = group->count(igroup); if (zeroflag) { if (particles == 0) diff --git a/src/USER-EFF/fix_nh_eff.cpp b/src/USER-EFF/fix_nh_eff.cpp index fb1a968ec5..1f996a565d 100644 --- a/src/USER-EFF/fix_nh_eff.cpp +++ b/src/USER-EFF/fix_nh_eff.cpp @@ -56,7 +56,6 @@ void FixNHEff::nve_v() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - int itype; double dtfm; for (int i = 0; i < nlocal; i++) { diff --git a/src/USER-EFF/fix_nvt_sllod_eff.cpp b/src/USER-EFF/fix_nvt_sllod_eff.cpp index 4aecdd4fde..95e52c1403 100644 --- a/src/USER-EFF/fix_nvt_sllod_eff.cpp +++ b/src/USER-EFF/fix_nvt_sllod_eff.cpp @@ -100,7 +100,7 @@ void FixNVTSllodEff::nh_v_temp() // calculate temperature since some computes require temp // computed on current nlocal atoms to remove bias - if (nondeformbias) double tmp = temperature->compute_scalar(); + if (nondeformbias) temperature->compute_scalar(); double **v = atom->v; double *ervel = atom->ervel; diff --git a/src/USER-MISC/angle_dipole.cpp b/src/USER-MISC/angle_dipole.cpp index 0762cf24fc..2d1b7eb3ef 100644 --- a/src/USER-MISC/angle_dipole.cpp +++ b/src/USER-MISC/angle_dipole.cpp @@ -52,7 +52,7 @@ void AngleDipole::compute(int eflag, int vflag) int iRef,iDip,iDummy,n,type; double delx,dely,delz; double eangle,tangle,f1[3],f3[3]; - double r,dr,cosGamma,deltaGamma,kdg,rmu; + double r,cosGamma,deltaGamma,kdg,rmu; eangle = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); diff --git a/src/USER-MISC/improper_cossq.cpp b/src/USER-MISC/improper_cossq.cpp index 53fd31e02c..e286277b60 100644 --- a/src/USER-MISC/improper_cossq.cpp +++ b/src/USER-MISC/improper_cossq.cpp @@ -56,7 +56,7 @@ ImproperCossq::~ImproperCossq() void ImproperCossq::compute(int eflag, int vflag) { - int i1,i2,i3,i4,m,n,type; + int i1,i2,i3,i4,n,type; double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z ; double eimproper,f1[3],f2[3],f3[3],f4[3]; double rjisq, rji, rlksq, rlk, cosphi, angfac; diff --git a/src/USER-MISC/pair_cdeam.cpp b/src/USER-MISC/pair_cdeam.cpp index b3ce1b18f4..cbb8e23f7c 100644 --- a/src/USER-MISC/pair_cdeam.cpp +++ b/src/USER-MISC/pair_cdeam.cpp @@ -342,7 +342,7 @@ void PairCDEAM::compute(int eflag, int vflag) x_j = rhoB[j]/rho[j]; ASSERT(x_j >= 0 && x_j<=1.0); - double D_j; + double D_j=0.0; if(cdeamVersion == 1) { // Calculate derivative of h(x_j) polynomial function. double h_prime_j = evalHprime(x_j); @@ -375,7 +375,7 @@ void PairCDEAM::compute(int eflag, int vflag) } else { // We have a concentration dependence for the i-j interaction. - double h; + double h=0.0; if(cdeamVersion == 1) { // Calculate h(x_i) polynomial function. double h_i = evalH(x_i); diff --git a/src/USER-MOLFILE/molfile_interface.cpp b/src/USER-MOLFILE/molfile_interface.cpp index f3b63e7eda..d0aff40f32 100644 --- a/src/USER-MOLFILE/molfile_interface.cpp +++ b/src/USER-MOLFILE/molfile_interface.cpp @@ -173,33 +173,6 @@ extern "C" { return pte_vdw_radius[idx]; } - static int get_pte_idx(const char *label) - { - int i; - char atom[3]; - - /* zap string */ - atom[0] = (char) 0; - atom[1] = (char) 0; - atom[2] = (char) 0; - /* if we don't have a null-pointer, there must be at least two - * chars, which is all we need. we convert to the capitalization - * convention of the table above during assignment. */ - if (label != NULL) { - atom[0] = (char) toupper((int) label[0]); - atom[1] = (char) tolower((int) label[1]); - } - /* discard numbers in atom label */ - if (isdigit(atom[1])) atom[1] = (char) 0; - - for (i=0; i < nr_pte_entries; ++i) { - if ( (pte_label[i][0] == atom[0]) - && (pte_label[i][1] == atom[1]) ) return i; - } - - return 0; - } - static int get_pte_idx_from_string(const char *label) { int i, ind; char atom[3]; @@ -296,15 +269,6 @@ extern "C" { return (void *)LoadLibrary(fname); } - // report error message from dlopen - static const char *my_dlerror(void) { - static CHAR szBuf[80]; - DWORD dw = GetLastError(); - - sprintf(szBuf, "my_dlopen failed: GetLastError returned %u\n", dw); - return szBuf; - } - // resolve a symbol in shared object static void *my_dlsym(void *h, const char *sym) { return (void *)GetProcAddress((HINSTANCE)h, sym); @@ -377,11 +341,6 @@ extern "C" { return dlopen(fname, RTLD_NOW); } - // report error message from dlopen - static const char *my_dlerror(void) { - return dlerror(); - } - // resolve a symbol in shared object static void *my_dlsym(void *h, const char *sym) { return dlsym(h, sym); diff --git a/src/USER-OMP/angle_class2_omp.cpp b/src/USER-OMP/angle_class2_omp.cpp index 6e61b0c9b4..bbe58ec48f 100644 --- a/src/USER-OMP/angle_class2_omp.cpp +++ b/src/USER-OMP/angle_class2_omp.cpp @@ -96,6 +96,7 @@ void AngleClass2OMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/angle_cosine_delta_omp.cpp b/src/USER-OMP/angle_cosine_delta_omp.cpp index 6fdd38acd4..a642694222 100644 --- a/src/USER-OMP/angle_cosine_delta_omp.cpp +++ b/src/USER-OMP/angle_cosine_delta_omp.cpp @@ -93,6 +93,7 @@ void AngleCosineDeltaOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/angle_cosine_omp.cpp b/src/USER-OMP/angle_cosine_omp.cpp index ed5f038b81..8aaddc528c 100644 --- a/src/USER-OMP/angle_cosine_omp.cpp +++ b/src/USER-OMP/angle_cosine_omp.cpp @@ -93,6 +93,7 @@ void AngleCosineOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/angle_cosine_periodic_omp.cpp b/src/USER-OMP/angle_cosine_periodic_omp.cpp index 0b13bfd34a..d63dfd8ed1 100644 --- a/src/USER-OMP/angle_cosine_periodic_omp.cpp +++ b/src/USER-OMP/angle_cosine_periodic_omp.cpp @@ -96,6 +96,7 @@ void AngleCosinePeriodicOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/angle_cosine_shift_exp_omp.cpp b/src/USER-OMP/angle_cosine_shift_exp_omp.cpp index 93fcfc7f8f..5abf621050 100644 --- a/src/USER-OMP/angle_cosine_shift_exp_omp.cpp +++ b/src/USER-OMP/angle_cosine_shift_exp_omp.cpp @@ -94,6 +94,7 @@ void AngleCosineShiftExpOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/angle_dipole_omp.cpp b/src/USER-OMP/angle_dipole_omp.cpp index 2f60f03bac..4e67801671 100644 --- a/src/USER-OMP/angle_dipole_omp.cpp +++ b/src/USER-OMP/angle_dipole_omp.cpp @@ -93,6 +93,7 @@ void AngleDipoleOMP::eval(int nfrom, int nto, ThrData * const thr) const int nlocal = atom->nlocal; const double f1[3] = {0.0, 0.0, 0.0}; const double f3[3] = {0.0, 0.0, 0.0}; + eangle = 0.0; for (n = nfrom; n < nto; n++) { iDip = anglelist[n][0]; // dipole whose orientation is to be restrained diff --git a/src/USER-OMP/angle_fourier_simple_omp.cpp b/src/USER-OMP/angle_fourier_simple_omp.cpp index 667dd8797f..9b27309e6c 100644 --- a/src/USER-OMP/angle_fourier_simple_omp.cpp +++ b/src/USER-OMP/angle_fourier_simple_omp.cpp @@ -94,6 +94,7 @@ void AngleFourierSimpleOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/angle_quartic_omp.cpp b/src/USER-OMP/angle_quartic_omp.cpp index a8119f51b7..26f783db79 100644 --- a/src/USER-OMP/angle_quartic_omp.cpp +++ b/src/USER-OMP/angle_quartic_omp.cpp @@ -94,6 +94,7 @@ void AngleQuarticOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/angle_sdk_omp.cpp b/src/USER-OMP/angle_sdk_omp.cpp index 60cb392f60..523f0836c7 100644 --- a/src/USER-OMP/angle_sdk_omp.cpp +++ b/src/USER-OMP/angle_sdk_omp.cpp @@ -96,6 +96,7 @@ void AngleSDKOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/angle_table_omp.cpp b/src/USER-OMP/angle_table_omp.cpp index d00e149634..b45956d54e 100644 --- a/src/USER-OMP/angle_table_omp.cpp +++ b/src/USER-OMP/angle_table_omp.cpp @@ -94,6 +94,7 @@ void AngleTableOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0]; const int nlocal = atom->nlocal; + eangle = 0.0; for (n = nfrom; n < nto; n++) { i1 = anglelist[n].a; diff --git a/src/USER-OMP/bond_class2_omp.cpp b/src/USER-OMP/bond_class2_omp.cpp index ab332299d6..47e684da1a 100644 --- a/src/USER-OMP/bond_class2_omp.cpp +++ b/src/USER-OMP/bond_class2_omp.cpp @@ -88,6 +88,7 @@ void BondClass2OMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int3_t * _noalias const bondlist = (int3_t *) neighbor->bondlist[0]; const int nlocal = atom->nlocal; + ebond = 0.0; for (n = nfrom; n < nto; n++) { i1 = bondlist[n].a; diff --git a/src/USER-OMP/bond_fene_expand_omp.cpp b/src/USER-OMP/bond_fene_expand_omp.cpp index e690aa215f..5a18461420 100644 --- a/src/USER-OMP/bond_fene_expand_omp.cpp +++ b/src/USER-OMP/bond_fene_expand_omp.cpp @@ -91,6 +91,7 @@ void BondFENEExpandOMP::eval(int nfrom, int nto, ThrData * const thr) const int3_t * _noalias const bondlist = (int3_t *) neighbor->bondlist[0]; const int nlocal = atom->nlocal; const int tid = thr->get_tid(); + ebond = 0.0; for (n = nfrom; n < nto; n++) { i1 = bondlist[n].a; diff --git a/src/USER-OMP/bond_fene_omp.cpp b/src/USER-OMP/bond_fene_omp.cpp index 164ba89195..bd7ed4a593 100644 --- a/src/USER-OMP/bond_fene_omp.cpp +++ b/src/USER-OMP/bond_fene_omp.cpp @@ -90,6 +90,7 @@ void BondFENEOMP::eval(int nfrom, int nto, ThrData * const thr) const int3_t * _noalias const bondlist = (int3_t *) neighbor->bondlist[0]; const int nlocal = atom->nlocal; const int tid = thr->get_tid(); + ebond = 0.0; for (n = nfrom; n < nto; n++) { i1 = bondlist[n].a; diff --git a/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp b/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp index 3d5d74acfa..1999912fae 100644 --- a/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp +++ b/src/USER-OMP/bond_harmonic_shift_cut_omp.cpp @@ -87,6 +87,7 @@ void BondHarmonicShiftCutOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int3_t * _noalias const bondlist = (int3_t *) neighbor->bondlist[0]; const int nlocal = atom->nlocal; + ebond = 0.0; for (n = nfrom; n < nto; n++) { i1 = bondlist[n].a; diff --git a/src/USER-OMP/bond_morse_omp.cpp b/src/USER-OMP/bond_morse_omp.cpp index 1d1af1b3db..2cae149e41 100644 --- a/src/USER-OMP/bond_morse_omp.cpp +++ b/src/USER-OMP/bond_morse_omp.cpp @@ -87,6 +87,7 @@ void BondMorseOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int3_t * _noalias const bondlist = (int3_t *) neighbor->bondlist[0]; const int nlocal = atom->nlocal; + ebond = 0.0; for (n = nfrom; n < nto; n++) { i1 = bondlist[n].a; diff --git a/src/USER-OMP/bond_table_omp.cpp b/src/USER-OMP/bond_table_omp.cpp index bd8fb5a3f3..3803eaa02b 100644 --- a/src/USER-OMP/bond_table_omp.cpp +++ b/src/USER-OMP/bond_table_omp.cpp @@ -88,6 +88,7 @@ void BondTableOMP::eval(int nfrom, int nto, ThrData * const thr) dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; const int3_t * _noalias const bondlist = (int3_t *) neighbor->bondlist[0]; const int nlocal = atom->nlocal; + ebond = 0.0; for (n = nfrom; n < nto; n++) { i1 = bondlist[n].a; diff --git a/src/USER-OMP/dihedral_fourier_omp.cpp b/src/USER-OMP/dihedral_fourier_omp.cpp index ef43313856..c8ab2bbf57 100644 --- a/src/USER-OMP/dihedral_fourier_omp.cpp +++ b/src/USER-OMP/dihedral_fourier_omp.cpp @@ -89,7 +89,7 @@ void DihedralFourierOMP::eval(int nfrom, int nto, ThrData * const thr) double edihedral,f1[3],f2[3],f3[3],f4[3]; double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; double df,df1_,ddf1_,fg,hg,fga,hgb,gaa,gbb; - double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; + double dtfx,dtfy,dtfz,dtgx,dtgy,dtgz,dthx,dthy,dthz; double c,s,p_,sx2,sy2,sz2; edihedral = 0.0; @@ -185,7 +185,7 @@ void DihedralFourierOMP::eval(int nfrom, int nto, ThrData * const thr) m = multiplicity[type][j]; p_ = 1.0; df1_ = 0.0; - + for (i = 0; i < m; i++) { ddf1_ = p_*c - df1_*s; df1_ = p_*s + df1_*c; @@ -196,13 +196,13 @@ void DihedralFourierOMP::eval(int nfrom, int nto, ThrData * const thr) df1_ = df1_*cos_shift[type][j] - ddf1_*sin_shift[type][j]; df1_ *= -m; p_ += 1.0; - + if (m == 0) { p_ = 1.0 + cos_shift[type][j]; df1_ = 0.0; } - if (EFLAG) edihedral += k[type][j] * p_; + if (EFLAG) edihedral += k[type][j] * p_; df += (-k[type][j] * df1_); } diff --git a/src/USER-OMP/fix_shear_history_omp.cpp b/src/USER-OMP/fix_shear_history_omp.cpp index 18c3c2af4e..0f5b82afa2 100644 --- a/src/USER-OMP/fix_shear_history_omp.cpp +++ b/src/USER-OMP/fix_shear_history_omp.cpp @@ -41,7 +41,6 @@ void FixShearHistoryOMP::pre_exchange() const int nlocal = atom->nlocal; const int nghost = atom->nghost; - const int nall = nlocal + nghost; const int nthreads = comm->nthreads; maxtouch = 0; diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp index bc785a2f81..5907460b7f 100644 --- a/src/USER-OMP/pair_airebo_omp.cpp +++ b/src/USER-OMP/pair_airebo_omp.cpp @@ -2659,7 +2659,6 @@ void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, void PairAIREBOOMP::REBO_neigh_thr() { - const int nlocal = atom->nlocal; const int nthreads = comm->nthreads; if (atom->nmax > maxlocal) { diff --git a/src/USER-OMP/pair_lj_cut_dipole_cut_omp.cpp b/src/USER-OMP/pair_lj_cut_dipole_cut_omp.cpp index 94900cca51..c7ec43af52 100755 --- a/src/USER-OMP/pair_lj_cut_dipole_cut_omp.cpp +++ b/src/USER-OMP/pair_lj_cut_dipole_cut_omp.cpp @@ -83,7 +83,7 @@ void PairLJCutDipoleCutOMP::eval(int iifrom, int iito, ThrData * const thr) double forcelj,factor_coul,factor_lj; int *ilist,*jlist,*numneigh,**firstneigh; - evdwl = 0.0; + evdwl = ecoul = 0.0; const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; diff --git a/src/USER-OMP/pair_lj_gromacs_coul_gromacs_omp.cpp b/src/USER-OMP/pair_lj_gromacs_coul_gromacs_omp.cpp index 1be1ba39eb..61fb667a69 100644 --- a/src/USER-OMP/pair_lj_gromacs_coul_gromacs_omp.cpp +++ b/src/USER-OMP/pair_lj_gromacs_coul_gromacs_omp.cpp @@ -122,6 +122,7 @@ void PairLJGromacsCoulGromacsOMP::eval(int iifrom, int iito, ThrData * const thr dely = ytmp - x[j].y; delz = ztmp - x[j].z; rsq = delx*delx + dely*dely + delz*delz; + tc = tlj = 0.0; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { diff --git a/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp b/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp index 2024a80dca..dbaf1a9639 100755 --- a/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp +++ b/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp @@ -87,7 +87,7 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) double rcutlj2inv, rcutcoul2inv,rcutlj6inv; int *ilist,*jlist,*numneigh,**firstneigh; - evdwl = 0.0; + evdwl = ecoul = 0.0; const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; diff --git a/src/USER-OMP/pair_resquared_omp.cpp b/src/USER-OMP/pair_resquared_omp.cpp index b35a2613c9..4e3bfded47 100644 --- a/src/USER-OMP/pair_resquared_omp.cpp +++ b/src/USER-OMP/pair_resquared_omp.cpp @@ -91,6 +91,7 @@ void PairRESquaredOMP::eval(int iifrom, int iito, ThrData * const thr) double fxtmp,fytmp,fztmp,t1tmp,t2tmp,t3tmp; + evdwl = 0.0; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; diff --git a/src/USER-QMMM/README b/src/USER-QMMM/README index 45f02e4b30..16bbdcc708 100644 --- a/src/USER-QMMM/README +++ b/src/USER-QMMM/README @@ -9,4 +9,4 @@ The person who created this package is Axel Kohlmeyer at Temple U --------------------------------- -Version: 2013-10-01 +Version: 2014-04-05 diff --git a/src/USER-REAXC/reaxc_allocate.cpp b/src/USER-REAXC/reaxc_allocate.cpp index 9352fd2c1d..7df5241a93 100644 --- a/src/USER-REAXC/reaxc_allocate.cpp +++ b/src/USER-REAXC/reaxc_allocate.cpp @@ -286,13 +286,11 @@ int Allocate_Workspace( reax_system *system, control_params *control, storage *workspace, int local_cap, int total_cap, MPI_Comm comm, char *msg ) { - int i, total_real, total_rvec, local_int, local_real, local_rvec; + int i, total_real, total_rvec, local_rvec; workspace->allocated = 1; total_real = total_cap * sizeof(real); total_rvec = total_cap * sizeof(rvec); - local_int = local_cap * sizeof(int); - local_real = local_cap * sizeof(real); local_rvec = local_cap * sizeof(rvec); /* communication storage */ @@ -551,10 +549,9 @@ int Estimate_GCell_Population( reax_system* system, MPI_Comm comm ) ivec c; grid *g; grid_cell *gc; - simulation_box *big_box, *my_ext_box; + simulation_box *my_ext_box; reax_atom *atoms; - big_box = &(system->big_box); my_ext_box = &(system->my_ext_box); g = &(system->my_grid); atoms = system->my_atoms; @@ -562,10 +559,6 @@ int Estimate_GCell_Population( reax_system* system, MPI_Comm comm ) for( l = 0; l < system->n; l++ ) { for( d = 0; d < 3; ++d ) { - //if( atoms[l].x[d] < big_box->min[d] ) - // atoms[l].x[d] += big_box->box_norms[d]; - //else if( atoms[l].x[d] >= big_box->max[d] ) - // atoms[l].x[d] -= big_box->box_norms[d]; c[d] = (int)((atoms[l].x[d]-my_ext_box->min[d])*g->inv_len[d]); @@ -757,7 +750,7 @@ void ReAllocate( reax_system *system, control_params *control, simulation_data *data, storage *workspace, reax_list **lists, mpi_datatypes *mpi_data ) { - int num_bonds, est_3body, nflag, Nflag, Hflag, ret; + int num_bonds, est_3body, Hflag, ret; int renbr, newsize; reallocate_data *realloc; reax_list *far_nbrs; @@ -792,14 +785,15 @@ void ReAllocate( reax_system *system, control_params *control, #endif // IMPORTANT: LOOSE ZONES CHECKS ARE DISABLED FOR NOW BY &&'ing with 0!!! - nflag = 0; + + int nflag = 0; if( system->n >= DANGER_ZONE * system->local_cap || (0 && system->n <= LOOSE_ZONE * system->local_cap) ) { nflag = 1; system->local_cap = MAX( (int)(system->n * safezone), mincap ); } - Nflag = 0; + int Nflag = 0; if( system->N >= DANGER_ZONE * system->total_cap || (0 && system->N <= LOOSE_ZONE * system->total_cap) ) { Nflag = 1; diff --git a/src/USER-REAXC/reaxc_io_tools.cpp b/src/USER-REAXC/reaxc_io_tools.cpp index b05a3fb8db..0b9d9b5ef9 100644 --- a/src/USER-REAXC/reaxc_io_tools.cpp +++ b/src/USER-REAXC/reaxc_io_tools.cpp @@ -613,7 +613,7 @@ void Print_Grid( grid* g, FILE *out ) } - +#if 0 void Print_GCell_Exchange_Bounds( int my_rank, neighbor_proc *my_nbrs ) { ivec r; @@ -653,6 +653,7 @@ void Print_GCell_Exchange_Bounds( int my_rank, neighbor_proc *my_nbrs ) fclose(f); } +#endif @@ -879,10 +880,9 @@ void Print_Symmetric_Sparse(reax_system *system, sparse_matrix *A, char *fname) void Print_Linear_System( reax_system *system, control_params *control, storage *workspace, int step ) { - int i, j; + int i; char fname[100]; - reax_atom *ai, *aj; - sparse_matrix *H; + reax_atom *ai; FILE *out; // print rhs and init guesses for QEq diff --git a/src/atom_vec_sphere.h b/src/atom_vec_sphere.h index fd82b7c5b1..c285267fff 100644 --- a/src/atom_vec_sphere.h +++ b/src/atom_vec_sphere.h @@ -73,7 +73,7 @@ class AtomVecSphere : public AtomVec { int *type,*mask; imageint *image; double **x,**v,**f; - double *radius,*density,*rmass; + double *radius,*rmass; double **omega,**torque; int radvary; }; diff --git a/src/balance.cpp b/src/balance.cpp index 010c2e8740..6073b25b64 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -11,6 +11,8 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +//#define BALANCE_DEBUG 1 + #include "lmptype.h" #include "mpi.h" #include "math.h" @@ -31,8 +33,6 @@ using namespace LAMMPS_NS; enum{NONE,UNIFORM,USER,DYNAMIC}; enum{X,Y,Z}; -//#define BALANCE_DEBUG 1 - /* ---------------------------------------------------------------------- */ Balance::Balance(LAMMPS *lmp) : Pointers(lmp) @@ -200,12 +200,13 @@ void Balance::command(int narg, char **arg) error->all(FLERR,"Illegal balance command"); if (dflag) { - for (int i = 0; i < strlen(bstr); i++) { + const int blen=strlen(bstr); + for (int i = 0; i < blen; i++) { if (bstr[i] != 'x' && bstr[i] != 'y' && bstr[i] != 'z') error->all(FLERR,"Balance dynamic string is invalid"); if (bstr[i] == 'z' && dimension == 2) error->all(FLERR,"Balance dynamic string is invalid"); - for (int j = i+1; j < strlen(bstr); j++) + for (int j = i+1; j < blen; j++) if (bstr[i] == bstr[j]) error->all(FLERR,"Balance dynamic string is invalid"); } @@ -433,7 +434,7 @@ void Balance::static_setup(char *str) ndim = strlen(str); bdim = new int[ndim]; - for (int i = 0; i < strlen(str); i++) { + for (int i = 0; i < ndim; i++) { if (str[i] == 'x') bdim[i] = X; if (str[i] == 'y') bdim[i] = Y; if (str[i] == 'z') bdim[i] = Z; @@ -866,11 +867,11 @@ void Balance::dumpout(bigint tstep, FILE *bfp) debug output for Idim and count only called by proc 0 ------------------------------------------------------------------------- */ - +#ifdef BALANCE_DEBUG void Balance::debug_output(int idim, int m, int np, double *split) { int i; - const char *dim; + const char *dim = NULL; double *boxlo = domain->boxlo; double *prd = domain->prd; @@ -878,41 +879,42 @@ void Balance::debug_output(int idim, int m, int np, double *split) if (bdim[idim] == X) dim = "X"; else if (bdim[idim] == Y) dim = "Y"; else if (bdim[idim] == Z) dim = "Z"; - printf("Dimension %s, Iteration %d\n",dim,m); + fprintf(stderr,"Dimension %s, Iteration %d\n",dim,m); - printf(" Count:"); - for (i = 0; i < np; i++) printf(" " BIGINT_FORMAT,count[i]); - printf("\n"); - printf(" Sum:"); - for (i = 0; i <= np; i++) printf(" " BIGINT_FORMAT,sum[i]); - printf("\n"); - printf(" Target:"); - for (i = 0; i <= np; i++) printf(" " BIGINT_FORMAT,target[i]); - printf("\n"); - printf(" Actual cut:"); + fprintf(stderr," Count:"); + for (i = 0; i < np; i++) fprintf(stderr," " BIGINT_FORMAT,count[i]); + fprintf(stderr,"\n"); + fprintf(stderr," Sum:"); + for (i = 0; i <= np; i++) fprintf(stderr," " BIGINT_FORMAT,sum[i]); + fprintf(stderr,"\n"); + fprintf(stderr," Target:"); + for (i = 0; i <= np; i++) fprintf(stderr," " BIGINT_FORMAT,target[i]); + fprintf(stderr,"\n"); + fprintf(stderr," Actual cut:"); for (i = 0; i <= np; i++) - printf(" %g",boxlo[bdim[idim]] + split[i]*prd[bdim[idim]]); - printf("\n"); - printf(" Split:"); - for (i = 0; i <= np; i++) printf(" %g",split[i]); - printf("\n"); - printf(" Low:"); - for (i = 0; i <= np; i++) printf(" %g",lo[i]); - printf("\n"); - printf(" Low-sum:"); - for (i = 0; i <= np; i++) printf(" " BIGINT_FORMAT,losum[i]); - printf("\n"); - printf(" Hi:"); - for (i = 0; i <= np; i++) printf(" %g",hi[i]); - printf("\n"); - printf(" Hi-sum:"); - for (i = 0; i <= np; i++) printf(" " BIGINT_FORMAT,hisum[i]); - printf("\n"); - printf(" Delta:"); - for (i = 0; i < np; i++) printf(" %g",split[i+1]-split[i]); - printf("\n"); + fprintf(stderr," %g",boxlo[bdim[idim]] + split[i]*prd[bdim[idim]]); + fprintf(stderr,"\n"); + fprintf(stderr," Split:"); + for (i = 0; i <= np; i++) fprintf(stderr," %g",split[i]); + fprintf(stderr,"\n"); + fprintf(stderr," Low:"); + for (i = 0; i <= np; i++) fprintf(stderr," %g",lo[i]); + fprintf(stderr,"\n"); + fprintf(stderr," Low-sum:"); + for (i = 0; i <= np; i++) fprintf(stderr," " BIGINT_FORMAT,losum[i]); + fprintf(stderr,"\n"); + fprintf(stderr," Hi:"); + for (i = 0; i <= np; i++) fprintf(stderr," %g",hi[i]); + fprintf(stderr,"\n"); + fprintf(stderr," Hi-sum:"); + for (i = 0; i <= np; i++) fprintf(stderr," " BIGINT_FORMAT,hisum[i]); + fprintf(stderr,"\n"); + fprintf(stderr," Delta:"); + for (i = 0; i < np; i++) fprintf(stderr," %g",split[i+1]-split[i]); + fprintf(stderr,"\n"); bigint max = 0; for (i = 0; i < np; i++) max = MAX(max,count[i]); - printf(" Imbalance factor: %g\n",1.0*max*np/target[np]); + fprintf(stderr," Imbalance factor: %g\n",1.0*max*np/target[np]); } +#endif diff --git a/src/balance.h b/src/balance.h index ee1d171546..013ed2189f 100644 --- a/src/balance.h +++ b/src/balance.h @@ -69,7 +69,9 @@ class Balance : protected Pointers { void tally(int, int, double *); int adjust(int, double *); int binary(double, int, double *); +#ifdef BALANCE_DEBUG void debug_output(int, int, int, double *); +#endif }; } diff --git a/src/compute_atom_molecule.h b/src/compute_atom_molecule.h index cd135e47d8..c93dcbee06 100644 --- a/src/compute_atom_molecule.h +++ b/src/compute_atom_molecule.h @@ -37,7 +37,7 @@ class ComputeAtomMolecule : public Compute { int nvalues,nmolecules; tagint idlo,idhi; - int *which,*argindex,*flavor,*value2index; + int *which,*argindex,*value2index; char **ids; int nstride,maxatom; diff --git a/src/compute_bond_local.h b/src/compute_bond_local.h index ede1bc6c17..9788ef9c8b 100644 --- a/src/compute_bond_local.h +++ b/src/compute_bond_local.h @@ -33,7 +33,7 @@ class ComputeBondLocal : public Compute { double memory_usage(); private: - int nvalues,dflag,eflag; + int nvalues; int ncount; int *bstyle; int singleflag; diff --git a/src/compute_pair_local.h b/src/compute_pair_local.h index 82851861c4..be4a23c614 100644 --- a/src/compute_pair_local.h +++ b/src/compute_pair_local.h @@ -34,7 +34,7 @@ class ComputePairLocal : public Compute { double memory_usage(); private: - int nvalues,dflag,eflag,fflag; + int nvalues; int ncount; int *pstyle; // style of each requested output diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index 9f87fcd09b..a0ee7ec01b 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -45,7 +45,7 @@ enum{PERATOM,LOCAL}; ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - int iarg; + int iarg = 0; if (strcmp(style,"reduce") == 0) { if (narg < 5) error->all(FLERR,"Illegal compute reduce command"); idregion = NULL; @@ -59,7 +59,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : idregion = new char[n]; strcpy(idregion,arg[3]); iarg = 4; - } + } else error->all(FLERR,"Unkown derived compute reduce style"); if (strcmp(arg[iarg],"sum") == 0) mode = SUM; else if (strcmp(arg[iarg],"min") == 0) mode = MINN; diff --git a/src/displace_atoms.cpp b/src/displace_atoms.cpp index 5701a04393..4a62647395 100644 --- a/src/displace_atoms.cpp +++ b/src/displace_atoms.cpp @@ -58,7 +58,7 @@ void DisplaceAtoms::command(int narg, char **arg) if (igroup == -1) error->all(FLERR,"Could not find displace_atoms group ID"); int groupbit = group->bitmask[igroup]; - int style; + int style=-1; if (strcmp(arg[1],"move") == 0) style = MOVE; else if (strcmp(arg[1],"ramp") == 0) style = RAMP; else if (strcmp(arg[1],"random") == 0) style = RANDOM; diff --git a/src/dump_cfg.h b/src/dump_cfg.h index 192074cd48..485b312b39 100755 --- a/src/dump_cfg.h +++ b/src/dump_cfg.h @@ -27,7 +27,7 @@ namespace LAMMPS_NS { class DumpCFG : public DumpCustom { public: DumpCFG(class LAMMPS *, int, char **); - ~DumpCFG(); + virtual ~DumpCFG(); private: char **auxname; // name strings of auxiliary properties diff --git a/src/fix.cpp b/src/fix.cpp index a2b772400d..e4b7e0d411 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -70,6 +70,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) // however, each fix that uses these values should explicitly set them nevery = 1; + global_freq = 1; maxvatom = 0; vatom = NULL; diff --git a/src/fix_ave_correlate.h b/src/fix_ave_correlate.h index d8b5510763..7928b3eda9 100644 --- a/src/fix_ave_correlate.h +++ b/src/fix_ave_correlate.h @@ -46,7 +46,6 @@ class FixAveCorrelate : public Fix { int type,ave,startstep,overwrite; double prefactor; - char *title1,*title2,*title3; long filepos; int firstindex; // index in values ring of earliest time sample diff --git a/src/fix_ave_histo.h b/src/fix_ave_histo.h index 5d3971014a..bb45e215a9 100644 --- a/src/fix_ave_histo.h +++ b/src/fix_ave_histo.h @@ -59,7 +59,7 @@ class FixAveHisto : public Fix { double *vector; int maxatom; - int ave,nwindow,nsum,startstep,mode; + int ave,nwindow,startstep,mode; char *title1,*title2,*title3; int iwindow,window_limit; diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index 000ae48c35..c9defd4078 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -730,10 +730,12 @@ void FixAveSpatial::end_of_step() MPI_DOUBLE,MPI_SUM,world); for (m = 0; m < nbins; m++) { if (count_sum[m] > 0.0) - for (j = 0; j < nvalues; j++) + for (j = 0; j < nvalues; j++) { if (which[j] == DENSITY_NUMBER) values_sum[m][j] /= repeat; - else if (which[j] == DENSITY_MASS) values_sum[m][j] *= mv2d/repeat; - else values_sum[m][j] /= count_sum[m]; + else if (which[j] == DENSITY_MASS) { + values_sum[m][j] *= mv2d/repeat; + } else values_sum[m][j] /= count_sum[m]; + } count_sum[m] /= repeat; } } else { diff --git a/src/fix_ave_time.h b/src/fix_ave_time.h index f35857ec57..4186d54580 100644 --- a/src/fix_ave_time.h +++ b/src/fix_ave_time.h @@ -47,7 +47,7 @@ class FixAveTime : public Fix { FILE *fp; int nrows; - int ave,nwindow,nsum,startstep,mode; + int ave,nwindow,startstep,mode; int noff,overwrite; int *offlist; char *title1,*title2,*title3; diff --git a/src/fix_balance.cpp b/src/fix_balance.cpp index 316140be3c..cf24cb9014 100644 --- a/src/fix_balance.cpp +++ b/src/fix_balance.cpp @@ -55,12 +55,13 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) : if (nevery < 0 || nitermax <= 0 || thresh < 1.0) error->all(FLERR,"Illegal fix balance command"); - for (int i = 0; i < strlen(bstr); i++) { + int blen = strlen(bstr); + for (int i = 0; i < blen; i++) { if (bstr[i] != 'x' && bstr[i] != 'y' && bstr[i] != 'z') error->all(FLERR,"Fix balance string is invalid"); if (bstr[i] == 'z' && dimension == 2) error->all(FLERR,"Fix balance string is invalid for 2d simulation"); - for (int j = i+1; j < strlen(bstr); j++) + for (int j = i+1; j < blen; j++) if (bstr[i] == bstr[j]) error->all(FLERR,"Fix balance string is invalid"); } diff --git a/src/fix_indent.h b/src/fix_indent.h index 7fbeaad1a7..f4639147a8 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -39,7 +39,7 @@ class FixIndent : public Fix { double compute_vector(int); private: - int istyle,scaleflag,thermo_flag,eflag_enable,side; + int istyle,scaleflag,side; double k,k3; char *xstr,*ystr,*zstr,*rstr,*pstr; int xvar,yvar,zvar,rvar,pvar; diff --git a/src/fix_vector.h b/src/fix_vector.h index 0d3511d695..819668afba 100644 --- a/src/fix_vector.h +++ b/src/fix_vector.h @@ -37,7 +37,7 @@ class FixVector : public Fix { private: int nvalues; - int *which,*argindex,*value2index,*offcol; + int *which,*argindex,*value2index; char **ids; int ncount; // # of values currently in growing vector or array diff --git a/src/fix_wall_region.h b/src/fix_wall_region.h index 92d7b393d3..23bec22d20 100644 --- a/src/fix_wall_region.h +++ b/src/fix_wall_region.h @@ -44,7 +44,6 @@ class FixWallRegion : public Fix { int eflag; double ewall[4],ewall_all[4]; int nlevels_respa; - double dt; char *idregion; double coeff1,coeff2,coeff3,coeff4,offset; diff --git a/src/group.cpp b/src/group.cpp index 2e93562b1a..4144108f24 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -187,7 +187,7 @@ void Group::assign(int narg, char **arg) strcmp(arg[2],"<=") == 0 || strcmp(arg[2],">=") == 0 || strcmp(arg[2],"<>") == 0)) { - int condition; + int condition = -1; if (strcmp(arg[2],"<") == 0) condition = LT; else if (strcmp(arg[2],"<=") == 0) condition = LE; else if (strcmp(arg[2],">") == 0) condition = GT; diff --git a/src/input.cpp b/src/input.cpp index e8f325b612..2c90363c30 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -893,7 +893,7 @@ void Input::partition() { if (narg < 3) error->all(FLERR,"Illegal partition command"); - int yesflag; + int yesflag=-1; if (strcmp(arg[0],"yes") == 0) yesflag = 1; else if (strcmp(arg[0],"no") == 0) yesflag = 0; else error->all(FLERR,"Illegal partition command"); @@ -1347,7 +1347,7 @@ void Input::neighbor_command() void Input::newton() { - int newton_pair,newton_bond; + int newton_pair=1,newton_bond=1; if (narg == 1) { if (strcmp(arg[0],"off") == 0) newton_pair = newton_bond = 0; @@ -1364,15 +1364,9 @@ void Input::newton() force->newton_pair = newton_pair; - if (newton_bond == 0) { - if (domain->box_exist && force->newton_bond == 1) - error->all(FLERR,"Newton bond change after simulation box is defined"); - force->newton_bond = 0; - } else { - if (domain->box_exist && force->newton_bond == 0) - error->all(FLERR,"Newton bond change after simulation box is defined"); - force->newton_bond = 1; - } + if (domain->box_exist && (newton_bond != force->newton_bond)) + error->all(FLERR,"Newton bond change after simulation box is defined"); + force->newton_bond = newton_bond; if (newton_pair || newton_bond) force->newton = 1; else force->newton = 0; diff --git a/src/lammps.cpp b/src/lammps.cpp index a4daffc34a..f37b4106dd 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -67,6 +67,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) screen = NULL; logfile = NULL; + infile = NULL; // parse input switches @@ -200,6 +201,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) if (iarg+1 > narg) error->universe_all(FLERR,"Invalid command-line argument"); helpflag = 1; + citeflag = 0; iarg += 1; } else error->universe_all(FLERR,"Invalid command-line argument"); } @@ -270,7 +272,6 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) screen = universe->uscreen; logfile = universe->ulogfile; world = universe->uworld; - infile = NULL; if (universe->me == 0) { if (inflag == 0) infile = stdin; @@ -421,7 +422,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) // error check on accelerator packages if (cudaflag == 1 && kokkosflag == 1) - error->all(FLERR,"Cannot use -cuda on and -kokkos on"); + error->all(FLERR,"Cannot use -cuda on and -kokkos on together"); // create Cuda class if USER-CUDA installed, unless explicitly switched off // instantiation creates dummy Cuda class if USER-CUDA is not installed @@ -518,13 +519,20 @@ LAMMPS::~LAMMPS() delete citeme; if (universe->nworlds == 1) { + if (screen && screen != stdout) fclose(screen); if (logfile) fclose(logfile); + logfile = NULL; + if (screen != stdout) screen = NULL; } else { if (screen && screen != stdout) fclose(screen); if (logfile) fclose(logfile); if (universe->ulogfile) fclose(universe->ulogfile); + logfile = NULL; + if (screen != stdout) screen = NULL; } + if (infile && infile != stdin) fclose(infile); + if (world != universe->uworld) MPI_Comm_free(&world); delete cuda; @@ -632,7 +640,6 @@ void LAMMPS::destroy() delete force; delete group; delete output; - delete modify; // modify must come after output, force, update // since they delete fixes delete domain; // domain must come after modify diff --git a/src/lattice.cpp b/src/lattice.cpp index ab184a04f7..e5081d75f7 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -153,12 +153,12 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) } else if (strcmp(arg[iarg],"orient") == 0) { if (iarg+5 > narg) error->all(FLERR,"Illegal lattice command"); - int dim; + int dim = -1; if (strcmp(arg[iarg+1],"x") == 0) dim = 0; else if (strcmp(arg[iarg+1],"y") == 0) dim = 1; else if (strcmp(arg[iarg+1],"z") == 0) dim = 2; else error->all(FLERR,"Illegal lattice command"); - int *orient; + int *orient = NULL; if (dim == 0) orient = orientx; else if (dim == 1) orient = orienty; else if (dim == 2) orient = orientz; diff --git a/src/pair_lj_cut_coul_dsf.cpp b/src/pair_lj_cut_coul_dsf.cpp index a18e9ccf6b..1fac1ef44f 100644 --- a/src/pair_lj_cut_coul_dsf.cpp +++ b/src/pair_lj_cut_coul_dsf.cpp @@ -456,7 +456,7 @@ double PairLJCutCoulDSF::single(int i, int j, int itype, int jtype, double rsq, eng += factor_lj*philj; } - if (r < cut_coulsq) { + if (rsq < cut_coulsq) { phicoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); eng += phicoul; } diff --git a/src/random_mars.h b/src/random_mars.h index b31999f68c..d3f6fa1887 100644 --- a/src/random_mars.h +++ b/src/random_mars.h @@ -26,7 +26,7 @@ class RanMars : protected Pointers { double gaussian(); private: - int seed,save; + int save; double second; double *u; int i97,j97; From d4518d2fd149d56b60c1955c551a909e1d349c25 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 20:23:37 +0000 Subject: [PATCH 06/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11779 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/fix_bond_break.html | 54 +++++++++++++++++------------------------ doc/fix_bond_break.txt | 54 +++++++++++++++++------------------------ 2 files changed, 44 insertions(+), 64 deletions(-) diff --git a/doc/fix_bond_break.html b/doc/fix_bond_break.html index 4cae77cb34..dd9d4d9983 100644 --- a/doc/fix_bond_break.html +++ b/doc/fix_bond_break.html @@ -47,12 +47,14 @@ specified criteria. This can be used to model the dissolution of a polymer network due to stretching of the simulation box or other deformations. In this context, a bond means an interaction between a pair of atoms computed by the bond_style command. -Once the bond is broken it will be permanently deleted. This is -different than a pairwise bond-order potential such -as Tersoff or AIREBO which infers bonds and many-body interactions -based on the current geometry of a small cluster of atoms and -effectively creates and destroys bonds from timestep to timestep as -atoms move. +Once the bond is broken it will be permanently deleted, as will all +angle, dihedral, and improper interactions that bond is part of. +

+

This is different than a pairwise bond-order +potential such as Tersoff or AIREBO which infers bonds and many-body +interactions based on the current geometry of a small cluster of atoms +and effectively creates and destroys bonds and higher-order many-body +interactions from timestep to timestep as atoms move.

A check for possible bond breakage is performed every Nevery timesteps. If two bonded atoms I,J are further than a distance Rmax @@ -80,17 +82,21 @@ A uniform random number between 0.0 and 1.0 is generated and the eligible bond is only broken if the random number < fraction.

When a bond is broken, data structures within LAMMPS that store bond -topology are updated to reflect the breakage. This can also affect -subsequent computation of pairwise interactions involving the atoms in -the bond. See the Restriction section below for additional -information. +topology are updated to reflect the breakage. Likewise, if the bond +is part of a 3-body (angle) or 4-body (dihedral, improper) +interaction, that interaction is removed as well.

-

Computationally, each timestep this fix operates, it loops over bond -lists and computes distances between pairs of bonded atoms in the -list. It also communicates between neighboring processors to -coordinate which bonds are broken. Thus it will increase the cost of -a timestep. Thus you should be cautious about invoking this fix too -frequently. +

Computationally, each timestep this fix operates, it loops over all +the bonds in the system and computes distances between pairs of bonded +atoms. It also communicates between neighboring processors to +coordinate which bonds are broken. Moreover, if any bonds are broken, +neighbor lists must be immediately updated on the same timestep. This +is to insure that any pairwise interactions that should be turned "on" +due to a bond breaking, because they are no longer excluded by the +presence of the bond and the settings of the +special_bonds command, will be immediately +recognized. All of these operations increase the cost of a timestep. +Thus you should be cautious about invoking this fix too frequently.

You can dump out snapshots of the current bond topology via the dump local command. @@ -132,22 +138,6 @@ minimization. built with that package. See the Making LAMMPS section for more info.

-

Currently, there are 2 restrictions for using this fix. We may relax -these in the future if there are new models that would be enabled by -it. -

-

When a bond is broken, you might wish to turn off angle and dihedral -interactions that include that bond. However, LAMMPS does not check -for these angles and dihedrals, even if your simulation defines an -angle_style or -dihedral_style. -

-

This fix requires that the pairwise weightings defined by the -special_bonds command be *,1,1 for 1-3 and 1-4 -neighbors within the bond topology (the 1-2 setting is not -restricted). This means that the pairwise interaction of I with J's -other bond partners is unaffected by the breaking of the bond. -

Related commands:

fix bond/create, fix diff --git a/doc/fix_bond_break.txt b/doc/fix_bond_break.txt index c3b1a502c2..556bfc1dcb 100755 --- a/doc/fix_bond_break.txt +++ b/doc/fix_bond_break.txt @@ -36,12 +36,14 @@ specified criteria. This can be used to model the dissolution of a polymer network due to stretching of the simulation box or other deformations. In this context, a bond means an interaction between a pair of atoms computed by the "bond_style"_bond_style.html command. -Once the bond is broken it will be permanently deleted. This is -different than a "pairwise"_pair_style.html bond-order potential such -as Tersoff or AIREBO which infers bonds and many-body interactions -based on the current geometry of a small cluster of atoms and -effectively creates and destroys bonds from timestep to timestep as -atoms move. +Once the bond is broken it will be permanently deleted, as will all +angle, dihedral, and improper interactions that bond is part of. + +This is different than a "pairwise"_pair_style.html bond-order +potential such as Tersoff or AIREBO which infers bonds and many-body +interactions based on the current geometry of a small cluster of atoms +and effectively creates and destroys bonds and higher-order many-body +interactions from timestep to timestep as atoms move. A check for possible bond breakage is performed every {Nevery} timesteps. If two bonded atoms I,J are further than a distance {Rmax} @@ -69,17 +71,21 @@ A uniform random number between 0.0 and 1.0 is generated and the eligible bond is only broken if the random number < fraction. When a bond is broken, data structures within LAMMPS that store bond -topology are updated to reflect the breakage. This can also affect -subsequent computation of pairwise interactions involving the atoms in -the bond. See the Restriction section below for additional -information. +topology are updated to reflect the breakage. Likewise, if the bond +is part of a 3-body (angle) or 4-body (dihedral, improper) +interaction, that interaction is removed as well. -Computationally, each timestep this fix operates, it loops over bond -lists and computes distances between pairs of bonded atoms in the -list. It also communicates between neighboring processors to -coordinate which bonds are broken. Thus it will increase the cost of -a timestep. Thus you should be cautious about invoking this fix too -frequently. +Computationally, each timestep this fix operates, it loops over all +the bonds in the system and computes distances between pairs of bonded +atoms. It also communicates between neighboring processors to +coordinate which bonds are broken. Moreover, if any bonds are broken, +neighbor lists must be immediately updated on the same timestep. This +is to insure that any pairwise interactions that should be turned "on" +due to a bond breaking, because they are no longer excluded by the +presence of the bond and the settings of the +"special_bonds"_special_bonds.html command, will be immediately +recognized. All of these operations increase the cost of a timestep. +Thus you should be cautious about invoking this fix too frequently. You can dump out snapshots of the current bond topology via the "dump local"_dump.html command. @@ -121,22 +127,6 @@ This fix is part of the MC package. It is only enabled if LAMMPS was built with that package. See the "Making LAMMPS"_Section_start.html#start_3 section for more info. -Currently, there are 2 restrictions for using this fix. We may relax -these in the future if there are new models that would be enabled by -it. - -When a bond is broken, you might wish to turn off angle and dihedral -interactions that include that bond. However, LAMMPS does not check -for these angles and dihedrals, even if your simulation defines an -"angle_style"_angle_style.html or -"dihedral_style"_dihedral_style.html. - -This fix requires that the pairwise weightings defined by the -"special_bonds"_special_bonds.html command be *,1,1 for 1-3 and 1-4 -neighbors within the bond topology (the 1-2 setting is not -restricted). This means that the pairwise interaction of I with J's -other bond partners is unaffected by the breaking of the bond. - [Related commands:] "fix bond/create"_fix_bond_create.html, "fix From 8871abb708f56e5217ee9125c4de1be1140015e7 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 20:23:51 +0000 Subject: [PATCH 07/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11780 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/comm.cpp | 55 ++++++++++++++++++++++++++++++++++++++++++++++++ src/comm.h | 2 ++ src/force.cpp | 15 +++++++++++++ src/force.h | 1 + src/neighbor.cpp | 19 +++++++++++++---- 5 files changed, 88 insertions(+), 4 deletions(-) diff --git a/src/comm.cpp b/src/comm.cpp index 1bb82bd8c4..bf5213af53 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -1518,6 +1518,61 @@ void Comm::forward_comm_array(int n, double **array) } } +/* ---------------------------------------------------------------------- + exchange info provided with all 6 stencil neighbors +------------------------------------------------------------------------- */ + +int Comm::exchange_variable(int n, double *inbuf, double *&outbuf) +{ + int nsend,nrecv,nrecv1,nrecv2; + MPI_Request request; + MPI_Status status; + + nrecv = n; + if (nrecv > maxrecv) grow_recv(nrecv); + memcpy(buf_recv,inbuf,nrecv*sizeof(double)); + + // loop over dimensions + + for (int dim = 0; dim < 3; dim++) { + + // no exchange if only one proc in a dimension + + if (procgrid[dim] == 1) continue; + + // send/recv info in both directions using same buf_recv + // if 2 procs in dimension, single send/recv + // if more than 2 procs in dimension, send/recv to both neighbors + + nsend = nrecv; + MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][0],0, + &nrecv1,1,MPI_INT,procneigh[dim][1],0,world,&status); + nrecv += nrecv1; + if (procgrid[dim] > 2) { + MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][1],0, + &nrecv2,1,MPI_INT,procneigh[dim][0],0,world,&status); + nrecv += nrecv2; + } else nrecv2 = 0; + + if (nrecv > maxrecv) grow_recv(nrecv); + + MPI_Irecv(&buf_recv[nsend],nrecv1,MPI_DOUBLE,procneigh[dim][1],0, + world,&request); + MPI_Send(buf_recv,nsend,MPI_DOUBLE,procneigh[dim][0],0,world); + MPI_Wait(&request,&status); + + if (procgrid[dim] > 2) { + MPI_Irecv(&buf_recv[nsend+nrecv1],nrecv2,MPI_DOUBLE,procneigh[dim][0],0, + world,&request); + MPI_Send(buf_recv,nsend,MPI_DOUBLE,procneigh[dim][1],0,world); + MPI_Wait(&request,&status); + } + } + + outbuf = buf_recv; + return nrecv; +} + /* ---------------------------------------------------------------------- communicate inbuf around full ring of processors with messtag nbytes = size of inbuf = n datums * nper bytes diff --git a/src/comm.h b/src/comm.h index e667ba2c32..2c72359552 100644 --- a/src/comm.h +++ b/src/comm.h @@ -63,6 +63,8 @@ class Comm : protected Pointers { virtual void reverse_comm_dump(class Dump *); // reverse comm from a Dump void forward_comm_array(int, double **); // forward comm of array + int exchange_variable(int, double *, // exchange of info on neigh stencil + double *&); void ring(int, int, void *, int, void (*)(int, char *), // ring comm void *, int self = 1); int read_lines_from_file(FILE *, int, int, char *); // read/bcast file lines diff --git a/src/force.cpp b/src/force.cpp index 04dea56633..7e8e189e84 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -485,6 +485,21 @@ Improper *Force::new_improper(const char *style, const char *suffix, int &sflag) return NULL; } +/* ---------------------------------------------------------------------- + return ptr to current improper class or hybrid sub-class if matches style +------------------------------------------------------------------------- */ + +Improper *Force::improper_match(const char *style) +{ + if (strcmp(improper_style,style) == 0) return improper; + else if (strcmp(improper_style,"hybrid") == 0) { + ImproperHybrid *hybrid = (ImproperHybrid *) bond; + for (int i = 0; i < hybrid->nstyles; i++) + if (strcmp(hybrid->keywords[i],style) == 0) return hybrid->styles[i]; + } + return NULL; +} + /* ---------------------------------------------------------------------- new kspace style ------------------------------------------------------------------------- */ diff --git a/src/force.h b/src/force.h index bd72567812..a9efb08fbf 100644 --- a/src/force.h +++ b/src/force.h @@ -93,6 +93,7 @@ class Force : protected Pointers { void create_improper(const char *, const char *suffix = NULL); class Improper *new_improper(const char *, const char *, int &); + class Improper *improper_match(const char *); void create_kspace(int, char **, const char *suffix = NULL); class KSpace *new_kspace(int, char **, const char *, int &); diff --git a/src/neighbor.cpp b/src/neighbor.cpp index c38d8db7aa..37305effca 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -565,11 +565,16 @@ void Neighbor::init() if (processed) continue; - // pair and half: if there is a full non-occasional non-skip list + // pair and half and newton != 2: + // if there is a full non-occasional non-skip list // change this list to half_from_full and point at the full list // parent could be copy list or pair or fix + // could remove newton != 2 check if added half_from_full_no_newton_ghost + // option in neigh_derive.cpp and below in choose_build() + // this would require full list had ghost info + // would be useful when reax/c used in hybrid mode, e.g. with airebo - if (requests[i]->pair && requests[i]->half) { + if (requests[i]->pair && requests[i]->half && requests[i]->newton != 2) { for (j = 0; j < nrequest; j++) { if (!lists[j]) continue; if (requests[j]->full && requests[j]->occasional == 0 && @@ -924,8 +929,14 @@ void Neighbor::choose_build(int index, NeighRequest *rq) else pb = &Neighbor::skip_from; } else if (rq->half_from_full) { - if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton; - else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton; + if (rq->newton == 0) { + if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton; + else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton; + } else if (rq->newton == 1) { + pb = &Neighbor::half_from_full_newton; + } else if (rq->newton == 2) { + pb = &Neighbor::half_from_full_no_newton; + } } else if (rq->half) { if (style == NSQ) { From 631aa15a321bb9760f3233d0112532f0584f2bae Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 23:35:35 +0000 Subject: [PATCH 08/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11781 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MC/fix_bond_break.cpp | 557 +++++++++++++++++++++++++++++++++++--- src/MC/fix_bond_break.h | 31 ++- src/thermo.cpp | 54 +++- src/thermo.h | 5 + 4 files changed, 598 insertions(+), 49 deletions(-) diff --git a/src/MC/fix_bond_break.cpp b/src/MC/fix_bond_break.cpp index 7f00bf37d5..99f78db824 100755 --- a/src/MC/fix_bond_break.cpp +++ b/src/MC/fix_bond_break.cpp @@ -19,6 +19,7 @@ #include "update.h" #include "respa.h" #include "atom.h" +#include "atom_vec.h" #include "force.h" #include "comm.h" #include "neighbor.h" @@ -30,6 +31,8 @@ using namespace LAMMPS_NS; using namespace FixConst; +#define DELTA 16 + /* ---------------------------------------------------------------------- */ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) : @@ -38,6 +41,7 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) : if (narg < 6) error->all(FLERR,"Illegal fix bond/break command"); MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); nevery = force->inumeric(FLERR,arg[3]); if (nevery <= 0) error->all(FLERR,"Illegal fix bond/break command"); @@ -50,7 +54,7 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) : extvector = 0; btype = force->inumeric(FLERR,arg[4]); - double cutoff = force->numeric(FLERR,arg[5]); + cutoff = force->numeric(FLERR,arg[5]); if (btype < 1 || btype > atom->nbondtypes) error->all(FLERR,"Invalid bond type in fix bond/break command"); @@ -86,8 +90,9 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) : random = new RanMars(lmp,seed + me); // set comm sizes needed by this fix + // forward is big due to comm of 1-2 neighbors - comm_forward = 2; + comm_forward = MAX(2,1+atom->maxspecial); comm_reverse = 2; // allocate arrays local to this fix @@ -96,6 +101,22 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) : partner = NULL; distsq = NULL; + maxbreak = maxbreakall = 0; + broken = brokenall = NULL; + inbuf = NULL; + + maxinfluenced = 0; + influenced = NULL; + + recvcounts = NULL; + displs = NULL; + + // copy = special list for one atom + // may contain 1-2 neighs of all 1-3 neighs before dedup() shrinks it + + int maxspecial = atom->maxspecial; + copy = new tagint[maxspecial*maxspecial + maxspecial]; + // zero out stats breakcount = 0; @@ -112,6 +133,13 @@ FixBondBreak::~FixBondBreak() memory->destroy(partner); memory->destroy(distsq); + memory->destroy(broken); + memory->destroy(brokenall); + memory->destroy(inbuf); + memory->destroy(influenced); + memory->destroy(recvcounts); + memory->destroy(displs); + delete [] copy; } /* ---------------------------------------------------------------------- */ @@ -128,33 +156,36 @@ int FixBondBreak::setmask() void FixBondBreak::init() { - // require special bonds = *,1,1 - // [0] can be anything b/c I,J are removed from each other's special list - // [1],[2] must be 1.0 b/c only special lists of I,J are updated when - // bond I-J is broken, not special lists of neighbors of I,J,etc - - int flag = 0; - if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) flag = 1; - if (force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) flag = 1; - if (flag) error->all(FLERR,"Fix bond/break requires special_bonds = *,1,1"); - - // warn if angles, dihedrals, impropers are being used - - if (force->angle || force->dihedral || force->improper) { - if (me == 0) - error->warning(FLERR,"Broken bonds will not alter angles, " - "dihedrals, or impropers"); - } - if (strstr(update->integrate_style,"respa")) nlevels_respa = ((Respa *) update->integrate)->nlevels; + + // commextent = 3*bondcutoff + // use 3 b/c 4 atom in 1-2-3-4 needs to know 1-2 bond has broken + // and that info could be known by atom 1 with + // atom 1 on one proc and atoms 2,3,4 on the other proc + // using "cutoff" as bond length is a guesstimate of what's OK + // see NOTE below for possible issues with this + + commextent = 3.0*cutoff; + + // improper class2 and ring styles not allowed for now + // due to different ordering of improper topology (not atom I centric) + + if (force->improper) { + if (force->improper_match("class2") || force->improper_match("ring")) + error->all(FLERR,"Cannot yet use fix bond/break with this " + "improper style"); + } + + // DEBUG + // print_bb(); } /* ---------------------------------------------------------------------- */ void FixBondBreak::post_integrate() { - int i,j,k,m,n,i1,i2,n1,n3,type; + int i,j,k,m,n,i1,i2,n1,n2,n3,type; double delx,dely,delz,rsq; tagint *slist; @@ -231,6 +262,7 @@ void FixBondBreak::post_integrate() if (partner[i]) probability[i] = random->uniform(); } + commflag = 1; comm->forward_comm_fix(this); // break bonds @@ -243,7 +275,7 @@ void FixBondBreak::post_integrate() int **nspecial = atom->nspecial; tagint **special = atom->special; - int nbreak = 0; + nbreak = 0; for (i = 0; i < nlocal; i++) { if (partner[i] == 0) continue; j = atom->map(partner[i]); @@ -278,17 +310,26 @@ void FixBondBreak::post_integrate() slist = special[i]; n1 = nspecial[i][0]; - n3 = nspecial[i][2]; for (m = 0; m < n1; m++) if (slist[m] == partner[i]) break; + n3 = nspecial[i][2]; for (; m < n3-1; m++) slist[m] = slist[m+1]; nspecial[i][0]--; nspecial[i][1]--; nspecial[i][2]--; - // count the broken bond once + // count the broken bond once and store in broken list - if (tag[i] < tag[j]) nbreak++; + if (tag[i] < tag[j]) { + if (nbreak == maxbreak) { + maxbreak += DELTA; + memory->grow(broken,maxbreak,2,"bond/break:broken"); + memory->grow(inbuf,2*maxbreak,"bond/break:inbuf"); + } + broken[nbreak][0] = tag[i]; + broken[nbreak][1] = tag[j]; + nbreak++; + } } // tally stats @@ -297,9 +338,370 @@ void FixBondBreak::post_integrate() breakcounttotal += breakcount; atom->nbonds -= breakcount; - // trigger reneighboring if any bonds were formed + // trigger reneighboring if any bonds were broken + // this insures neigh lists will immediately reflect the topology changes + // done if no bonds broken if (breakcount) next_reneighbor = update->ntimestep; + if (!breakcount) return; + + // communicate broken bonds to procs that need them + // local comm via exchange_variable() if commextent < prob sub-domain + // else global comm via MPI_Allgatherv() + // NOTE: not fully happy with this test, + // but want to avoid Allgather if at all possible + // issue is there is no simple way to guarantee that local comm is OK + // what is really needed is every 4 atom in 1-2-3-4 knows about 1-2 bond + // this is hard to predict dynamically b/c current lengths of 2-3,3-4 + // bonds are not known + // could use 3*(domain->maxbondall+BONDSTRETCH) for better estimate? + // am not doing the local comm via ghost atoms, so ghost cutoff + // is irrelevant + // this test is also for orthogonal boxes and equi-partitioning + + double *outbuf; + + int local = 1; + if (domain->xprd/comm->procgrid[0] < commextent) local = 0; + if (domain->yprd/comm->procgrid[1] < commextent) local = 0; + if (domain->dimension == 3 && domain->zprd/comm->procgrid[2] < commextent) + local = 0; + + if (local) { + m = 0; + for (i = 0; i < nbreak; i++) { + inbuf[m++] = broken[i][0]; + inbuf[m++] = broken[i][1]; + } + int ntotal = comm->exchange_variable(2*nbreak,inbuf,outbuf); + nbreakall = ntotal/2; + } else MPI_Allreduce(&nbreak,&nbreakall,1,MPI_INT,MPI_SUM,world); + + if (nbreakall > maxbreakall) { + while (maxbreakall < nbreakall) maxbreakall += DELTA; + memory->destroy(brokenall); + memory->grow(brokenall,maxbreakall,2,"bond/break:brokenall"); + } + + if (local) { + m = 0; + for (i = 0; i < nbreakall; i++) { + brokenall[i][0] = static_cast (outbuf[m++]); + brokenall[i][1] = static_cast (outbuf[m++]); + } + } else { + if (!recvcounts) { + memory->create(recvcounts,nprocs,"bond/break:recvcounts"); + memory->create(displs,nprocs,"bond/break:displs"); + } + int n = 2*nbreak; + MPI_Allgather(&n,1,MPI_INT,recvcounts,1,MPI_INT,world); + displs[0] = 0; + for (int iproc = 1; iproc < nprocs; iproc++) + displs[iproc] = displs[iproc-1] + recvcounts[iproc-1]; + + int *ptr = NULL; + if (nbreak) ptr = &broken[0][0]; + MPI_Allgatherv(ptr,2*nbreak,MPI_INT, + &brokenall[0][0],recvcounts,displs,MPI_INT,world); + } + + // communicate 1-2 special neighs of ghost atoms + // 1-2 neighs already reflect broken bonds + + commflag = 2; + comm->forward_comm_variable_fix(this); + + // update special neigh lists of all atoms affected by any broken bond + // also remove angles/dihedrals/impropers broken by broken bonds + + update_topology(); + + // DEBUG + // print_bb(); +} + +/* ---------------------------------------------------------------------- + double loop over my atoms and broken bonds + influenced = 1 if atom's topology is affected by any broken bond + yes if is one of 2 atoms in bond + yes if both atom IDs appear in atom's special list + else no + if influenced: + check for angles/dihedrals/impropers to break due to broken bond + rebuild the atom's special list of 1-2,1-3,1-4 neighs +------------------------------------------------------------------------- */ + +void FixBondBreak::update_topology() +{ + int i,j,k,n,influence,found; + tagint id1,id2; + tagint *slist; + + int angles_allow = atom->avec->angles_allow; + int dihedrals_allow = atom->avec->dihedrals_allow; + int impropers_allow = atom->avec->impropers_allow; + + tagint *tag = atom->tag; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + int nlocal = atom->nlocal; + + if (nlocal > maxinfluenced) { + maxinfluenced = atom->nmax; + memory->destroy(influenced); + memory->create(influenced,maxinfluenced,"bond/break:influenced"); + } + + nangles = 0; + ndihedrals = 0; + nimpropers = 0; + + for (i = 0; i < nlocal; i++) { + influenced[i] = 0; + slist = special[i]; + + for (j = 0; j < nbreakall; j++) { + id1 = brokenall[j][0]; + id2 = brokenall[j][1]; + + influence = 0; + if (tag[i] == id1 || tag[i] == id2) influence = 1; + else { + n = nspecial[i][2]; + found = 0; + for (k = 0; k < n; k++) + if (slist[k] == id1 || slist[k] == id2) found++; + if (found == 2) influence = 1; + } + if (!influence) continue; + influenced[i] = 1; + + if (angles_allow) break_angles(i,id1,id2); + if (dihedrals_allow) break_dihedrals(i,id1,id2); + if (impropers_allow) break_impropers(i,id1,id2); + } + } + + int newton_bond = force->newton_bond; + + int all; + if (angles_allow) { + MPI_Allreduce(&nangles,&all,1,MPI_INT,MPI_SUM,world); + if (!newton_bond) all /= 3; + atom->nangles -= all; + } + if (dihedrals_allow) { + MPI_Allreduce(&ndihedrals,&all,1,MPI_INT,MPI_SUM,world); + if (!newton_bond) all /= 4; + atom->ndihedrals -= all; + } + if (impropers_allow) { + MPI_Allreduce(&nimpropers,&all,1,MPI_INT,MPI_SUM,world); + if (!newton_bond) all /= 4; + atom->nimpropers -= all; + } + + for (i = 0; i < nlocal; i++) + if (influenced[i]) rebuild_special(i,id1,id2); +} + +/* ---------------------------------------------------------------------- +------------------------------------------------------------------------- */ + +void FixBondBreak::rebuild_special(int m, tagint id1, tagint id2) +{ + int i,j,n,n1,n2,n3,cn1,cn2,cn3; + tagint *slist; + + tagint *tag = atom->tag; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + + // new 1-2 neighs of atom M + + slist = special[m]; + n1 = nspecial[m][0]; + cn1 = 0; + for (i = 0; i < n1; i++) + copy[cn1++] = slist[i]; + + // new 1-3 neighs of atom M, based on 1-2 neighs of 1-2 neighs + // exclude self + // remove duplicates after adding all possible 1-3 neighs + + cn2 = cn1; + for (i = 0; i < cn1; i++) { + n = atom->map(copy[i]); + slist = special[n]; + n1 = nspecial[n][0]; + for (j = 0; j < n1; j++) + if (slist[j] != tag[m]) copy[cn2++] = slist[j]; + } + + cn2 = dedup(cn1,cn2,copy); + + // new 1-4 neighs of atom M, based on 1-2 neighs of 1-3 neighs + // exclude self + // remove duplicates after adding all possible 1-4 neighs + + cn3 = cn2; + for (i = cn1; i < cn2; i++) { + n = atom->map(copy[i]); + slist = special[n]; + n1 = nspecial[n][0]; + for (j = 0; j < n1; j++) + if (slist[j] != tag[m]) copy[cn3++] = slist[j]; + } + + cn3 = dedup(cn2,cn3,copy); + + // store new special list with atom M + + nspecial[m][0] = cn1; + nspecial[m][1] = cn2; + nspecial[m][2] = cn3; + memcpy(special[m],copy,cn3*sizeof(int)); +} + +/* ---------------------------------------------------------------------- +------------------------------------------------------------------------- */ + +void FixBondBreak::break_angles(int m, tagint id1, tagint id2) +{ + int j,found; + + int num_angle = atom->num_angle[m]; + int *angle_type = atom->angle_type[m]; + tagint *angle_atom1 = atom->angle_atom1[m]; + tagint *angle_atom2 = atom->angle_atom2[m]; + tagint *angle_atom3 = atom->angle_atom3[m]; + + int i = 0; + while (i < num_angle) { + found = 0; + if (angle_atom1[i] == id1 && angle_atom2[i] == id2) found = 1; + else if (angle_atom2[i] == id1 && angle_atom3[i] == id2) found = 1; + else if (angle_atom1[i] == id2 && angle_atom2[i] == id1) found = 1; + else if (angle_atom2[i] == id2 && angle_atom3[i] == id1) found = 1; + if (!found) i++; + else { + for (j = i; j < num_angle-1; j++) { + angle_type[j] = angle_type[j+1]; + angle_atom1[j] = angle_atom1[j+1]; + angle_atom2[j] = angle_atom2[j+1]; + angle_atom3[j] = angle_atom3[j+1]; + } + num_angle--; + nangles++; + } + } + + atom->num_angle[m] = num_angle; +} + +/* ---------------------------------------------------------------------- +------------------------------------------------------------------------- */ + +void FixBondBreak::break_dihedrals(int m, tagint id1, tagint id2) +{ + int j,found; + + int num_dihedral = atom->num_dihedral[m]; + int *dihedral_type = atom->dihedral_type[m]; + tagint *dihedral_atom1 = atom->dihedral_atom1[m]; + tagint *dihedral_atom2 = atom->dihedral_atom2[m]; + tagint *dihedral_atom3 = atom->dihedral_atom3[m]; + tagint *dihedral_atom4 = atom->dihedral_atom4[m]; + + int i = 0; + while (i < num_dihedral) { + found = 0; + if (dihedral_atom1[i] == id1 && dihedral_atom2[i] == id2) found = 1; + else if (dihedral_atom2[i] == id1 && dihedral_atom3[i] == id2) found = 1; + else if (dihedral_atom3[i] == id1 && dihedral_atom4[i] == id2) found = 1; + else if (dihedral_atom1[i] == id2 && dihedral_atom2[i] == id1) found = 1; + else if (dihedral_atom2[i] == id2 && dihedral_atom3[i] == id1) found = 1; + else if (dihedral_atom3[i] == id2 && dihedral_atom4[i] == id1) found = 1; + if (!found) i++; + else { + for (j = i; j < num_dihedral-1; j++) { + dihedral_type[j] = dihedral_type[j+1]; + dihedral_atom1[j] = dihedral_atom1[j+1]; + dihedral_atom2[j] = dihedral_atom2[j+1]; + dihedral_atom3[j] = dihedral_atom3[j+1]; + dihedral_atom4[j] = dihedral_atom4[j+1]; + } + num_dihedral--; + ndihedrals++; + } + } + + atom->num_dihedral[m] = num_dihedral; +} + +/* ---------------------------------------------------------------------- +------------------------------------------------------------------------- */ + +void FixBondBreak::break_impropers(int m, tagint id1, tagint id2) +{ + int j,found; + + int num_improper = atom->num_improper[m]; + int *improper_type = atom->improper_type[m]; + tagint *improper_atom1 = atom->improper_atom1[m]; + tagint *improper_atom2 = atom->improper_atom2[m]; + tagint *improper_atom3 = atom->improper_atom3[m]; + tagint *improper_atom4 = atom->improper_atom4[m]; + + int i = 0; + while (i < num_improper) { + found = 0; + if (improper_atom1[i] == id1 && improper_atom2[i] == id2) found = 1; + else if (improper_atom1[i] == id1 && improper_atom3[i] == id2) found = 1; + else if (improper_atom1[i] == id1 && improper_atom4[i] == id2) found = 1; + else if (improper_atom1[i] == id2 && improper_atom2[i] == id1) found = 1; + else if (improper_atom1[i] == id2 && improper_atom3[i] == id1) found = 1; + else if (improper_atom1[i] == id2 && improper_atom4[i] == id1) found = 1; + if (!found) i++; + else { + for (j = i; j < num_improper-1; j++) { + improper_type[j] = improper_type[j+1]; + improper_atom1[j] = improper_atom1[j+1]; + improper_atom2[j] = improper_atom2[j+1]; + improper_atom3[j] = improper_atom3[j+1]; + improper_atom4[j] = improper_atom4[j+1]; + } + num_improper--; + nimpropers++; + } + } + + atom->num_improper[m] = num_improper; +} + +/* ---------------------------------------------------------------------- + remove all ID duplicates in copy from Nstart:Nstop-1 + compare to all previous values in copy + return N decremented by any discarded duplicates +------------------------------------------------------------------------- */ + +int FixBondBreak::dedup(int nstart, int nstop, tagint *copy) +{ + int i; + + int m = nstart; + while (m < nstop) { + for (i = 0; i < m; i++) + if (copy[i] == copy[m]) { + copy[m] = copy[nstop-1]; + nstop--; + break; + } + if (i == m) m++; + } + + return nstop; } /* ---------------------------------------------------------------------- */ @@ -312,30 +714,61 @@ void FixBondBreak::post_integrate_respa(int ilevel, int iloop) /* ---------------------------------------------------------------------- */ int FixBondBreak::pack_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) + int pbc_flag, int *pbc) { - int i,j,m; + int i,j,k,m,ns; + + if (commflag == 1) { + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = ubuf(partner[j]).d; + buf[m++] = probability[j]; + } + return 2; + } + + int **nspecial = atom->nspecial; + int **special = atom->special; m = 0; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = partner[j]; - buf[m++] = probability[j]; + ns = nspecial[j][0]; + buf[m++] = ubuf(ns).d; + for (k = 0; k < ns; k++) + buf[m++] = ubuf(special[j][k]).d; } - return 2; + return m; } /* ---------------------------------------------------------------------- */ void FixBondBreak::unpack_comm(int n, int first, double *buf) { - int i,m,last; + int i,j,m,ns,last; - m = 0; - last = first + n; - for (i = first; i < last; i++) { - partner[i] = static_cast (buf[m++]); - probability[i] = buf[m++]; + if (commflag == 1) { + m = 0; + last = first + n; + for (i = first; i < last; i++) { + partner[i] = (tagint) ubuf(buf[m++]).i; + probability[i] = buf[m++]; + } + + } else { + + int **nspecial = atom->nspecial; + int **special = atom->special; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + ns = (int) ubuf(buf[m++]).i; + nspecial[i][0] = ns; + for (j = 0; j < ns; j++) + special[i][j] = (tagint) ubuf(buf[m++]).i; + } } } @@ -348,7 +781,7 @@ int FixBondBreak::pack_reverse_comm(int n, int first, double *buf) m = 0; last = first + n; for (i = first; i < last; i++) { - buf[m++] = partner[i]; + buf[m++] = ubuf(partner[i]).d; buf[m++] = distsq[i]; } return 2; @@ -364,12 +797,55 @@ void FixBondBreak::unpack_reverse_comm(int n, int *list, double *buf) for (i = 0; i < n; i++) { j = list[i]; if (buf[m+1] > distsq[j]) { - partner[j] = static_cast (buf[m++]); + partner[j] = (tagint) ubuf(buf[m++]).i; distsq[j] = buf[m++]; } else m += 2; } } + +/* ---------------------------------------------------------------------- */ + +void FixBondBreak::print_bb() +{ + for (int i = 0; i < atom->nlocal; i++) { + printf("TAG %i: %d nbonds: ",atom->tag[i],atom->num_bond[i]); + for (int j = 0; j < atom->num_bond[i]; j++) { + printf(" %d",atom->bond_atom[i][j]); + } + printf("\n"); + printf("TAG %i: %d nangles: ",atom->tag[i],atom->num_angle[i]); + for (int j = 0; j < atom->num_angle[i]; j++) { + printf(" %d %d %d,",atom->angle_atom1[i][j], + atom->angle_atom2[i][j],atom->angle_atom3[i][j]); + } + printf("\n"); + printf("TAG %i: %d ndihedrals: ",atom->tag[i],atom->num_dihedral[i]); + for (int j = 0; j < atom->num_dihedral[i]; j++) { + printf(" %d %d %d %d,",atom->dihedral_atom1[i][j], + atom->dihedral_atom2[i][j],atom->dihedral_atom3[i][j], + atom->dihedral_atom4[i][j]); + } + printf("\n"); + printf("TAG %i: %d %d %d nspecial: ",atom->tag[i], + atom->nspecial[i][0],atom->nspecial[i][1],atom->nspecial[i][2]); + for (int j = 0; j < atom->nspecial[i][2]; j++) { + printf(" %d",atom->special[i][j]); + } + printf("\n"); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixBondBreak::print_copy(const char *str, tagint m, + int n1, int n2, int n3, int *v) +{ + printf("%s %i: %d %d %d nspecial: ",str,m,n1,n2,n3); + for (int j = 0; j < n3; j++) printf(" %d",v[j]); + printf("\n"); +} + /* ---------------------------------------------------------------------- */ double FixBondBreak::compute_vector(int n) @@ -385,7 +861,8 @@ double FixBondBreak::compute_vector(int n) double FixBondBreak::memory_usage() { int nmax = atom->nmax; - double bytes = nmax * sizeof(int); + double bytes = nmax * sizeof(tagint); bytes += nmax * sizeof(double); + bytes += maxinfluenced * sizeof(int); return bytes; } diff --git a/src/MC/fix_bond_break.h b/src/MC/fix_bond_break.h index 009d516429..29c00c0e53 100755 --- a/src/MC/fix_bond_break.h +++ b/src/MC/fix_bond_break.h @@ -41,17 +41,44 @@ class FixBondBreak : public Fix { double memory_usage(); private: - int me; + int me,nprocs; int btype,seed; - double cutsq,fraction; + double cutoff,cutsq,fraction; + double commextent; int breakcount,breakcounttotal; int nmax; tagint *partner; double *distsq,*probability; + int nbreak,nbreakall; + int maxbreak,maxbreakall; + tagint **broken,**brokenall; + double *inbuf; + + int maxinfluenced; + int *influenced; + int *recvcounts,*displs; + tagint *copy; + class RanMars *random; int nlevels_respa; + + int commflag; + int nbroken; + int nangles,ndihedrals,nimpropers; + + void update_topology(); + void rebuild_special(int, tagint, tagint); + void break_angles(int, tagint, tagint); + void break_dihedrals(int, tagint, tagint); + void break_impropers(int, tagint, tagint); + int dedup(int, int, tagint *); + + // DEBUG + + void print_bb(); + void print_copy(const char *, tagint, int, int, int, int *); }; } diff --git a/src/thermo.cpp b/src/thermo.cpp index 3e5b5543ef..ce5a481d4b 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -52,6 +52,7 @@ using namespace MathConst; // evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail // vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, // xlat, ylat, zlat +// bonds, angles, dihedrals, impropers, // pxx, pyy, pzz, pxy, pxz, pyz // fmax, fnorm // cella, cellb, cellc, cellalpha, cellbeta, cellgamma @@ -764,6 +765,15 @@ void Thermo::parse_fields(char *str) } else if (strcmp(word,"zlat") == 0) { addfield("Zlat",&Thermo::compute_zlat,FLOAT); + } else if (strcmp(word,"bonds") == 0) { + addfield("Bonds",&Thermo::compute_bonds,BIGINT); + } else if (strcmp(word,"angles") == 0) { + addfield("Angles",&Thermo::compute_angles,BIGINT); + } else if (strcmp(word,"dihedrals") == 0) { + addfield("Diheds",&Thermo::compute_dihedrals,BIGINT); + } else if (strcmp(word,"impropers") == 0) { + addfield("Impros",&Thermo::compute_impropers,BIGINT); + } else if (strcmp(word,"pxx") == 0) { addfield("Pxx",&Thermo::compute_pxx,FLOAT); index_press_vector = add_compute(id_press,VECTOR); @@ -1266,14 +1276,16 @@ int Thermo::evaluate_keyword(char *word, double *answer) else if (strcmp(word,"xz") == 0) compute_xz(); else if (strcmp(word,"yz") == 0) compute_yz(); - else if (strcmp(word,"xlat") == 0) { - compute_xlat(); - } else if (strcmp(word,"ylat") == 0) { - compute_ylat(); - } else if (strcmp(word,"zlat") == 0) { - compute_zlat(); + else if (strcmp(word,"xlat") == 0) compute_xlat(); + else if (strcmp(word,"ylat") == 0) compute_ylat(); + else if (strcmp(word,"zlat") == 0) compute_zlat(); - } else if (strcmp(word,"pxx") == 0) { + else if (strcmp(word,"bonds") == 0) compute_bonds(); + else if (strcmp(word,"angles") == 0) compute_angles(); + else if (strcmp(word,"dihedrals") == 0) compute_dihedrals(); + else if (strcmp(word,"impropers") == 0) compute_impropers(); + + else if (strcmp(word,"pxx") == 0) { if (!pressure) error->all(FLERR,"Thermo keyword in variable requires " "thermo to use/init press"); @@ -1857,6 +1869,34 @@ void Thermo::compute_zlat() /* ---------------------------------------------------------------------- */ +void Thermo::compute_bonds() +{ + bivalue = atom->nbonds; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_angles() +{ + bivalue = atom->nangles; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_dihedrals() +{ + bivalue = atom->ndihedrals; +} + +/* ---------------------------------------------------------------------- */ + +void Thermo::compute_impropers() +{ + bivalue = atom->nimpropers; +} + +/* ---------------------------------------------------------------------- */ + void Thermo::compute_pxx() { dvalue = pressure->vector[0]; diff --git a/src/thermo.h b/src/thermo.h index 149f98dd16..6dbbdd7e05 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -169,6 +169,11 @@ class Thermo : protected Pointers { void compute_ylat(); void compute_zlat(); + void compute_bonds(); + void compute_angles(); + void compute_dihedrals(); + void compute_impropers(); + void compute_pxx(); void compute_pyy(); void compute_pzz(); From a939f8aab733429e00ede274bed42c877ac53821 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 23:37:01 +0000 Subject: [PATCH 09/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11782 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/thermo_style.html | 2 ++ doc/thermo_style.txt | 2 ++ 2 files changed, 4 insertions(+) diff --git a/doc/thermo_style.html b/doc/thermo_style.html index 2c43fde256..f8c8030406 100644 --- a/doc/thermo_style.html +++ b/doc/thermo_style.html @@ -29,6 +29,7 @@ emol, elong, etail, vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, xlat, ylat, zlat, + bonds, angles, dihedrals, impropers, pxx, pyy, pzz, pxy, pxz, pyz, fmax, fnorm, cella, cellb, cellc, cellalpha, cellbeta, cellgamma, @@ -68,6 +69,7 @@ xlo,xhi,ylo,yhi,zlo,zhi = box boundaries xy,xz,yz = box tilt for triclinic (non-orthogonal) simulation boxes xlat,ylat,zlat = lattice spacings as calculated by lattice command + bonds,angles,dihedrals,impropers = # of these interactions defined pxx,pyy,pzz,pxy,pxz,pyz = 6 components of pressure tensor fmax = max component of force on any atom in any dimension fnorm = length of force vector for all atoms diff --git a/doc/thermo_style.txt b/doc/thermo_style.txt index 0fbcbce8bb..f0adef0ab2 100644 --- a/doc/thermo_style.txt +++ b/doc/thermo_style.txt @@ -24,6 +24,7 @@ args = list of arguments for a particular style :l emol, elong, etail, vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, xlat, ylat, zlat, + bonds, angles, dihedrals, impropers, pxx, pyy, pzz, pxy, pxz, pyz, fmax, fnorm, cella, cellb, cellc, cellalpha, cellbeta, cellgamma, @@ -63,6 +64,7 @@ args = list of arguments for a particular style :l xlo,xhi,ylo,yhi,zlo,zhi = box boundaries xy,xz,yz = box tilt for triclinic (non-orthogonal) simulation boxes xlat,ylat,zlat = lattice spacings as calculated by "lattice"_lattice.html command + bonds,angles,dihedrals,impropers = # of these interactions defined pxx,pyy,pzz,pxy,pxz,pyz = 6 components of pressure tensor fmax = max component of force on any atom in any dimension fnorm = length of force vector for all atoms From 608c01bb8de12100e811487e21d8faef6dcf34af Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 23:37:40 +0000 Subject: [PATCH 10/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11783 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/version.h b/src/version.h index 7e739f872f..ee9cf1af8e 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "11 Apr 2014" +#define LAMMPS_VERSION "12 Apr 2014" From 5a26ea8238c2ef1f2f02a174eb3c3776793c13f4 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 11 Apr 2014 23:37:41 +0000 Subject: [PATCH 11/11] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11784 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Manual.html | 4 ++-- doc/Manual.txt | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/Manual.html b/doc/Manual.html index ee18255472..dce50bfbad 100644 --- a/doc/Manual.html +++ b/doc/Manual.html @@ -1,7 +1,7 @@ LAMMPS Users Manual - + @@ -22,7 +22,7 @@

LAMMPS Documentation

-

11 Apr 2014 version +

12 Apr 2014 version

Version info:

diff --git a/doc/Manual.txt b/doc/Manual.txt index a4642c9f19..6f58820859 100644 --- a/doc/Manual.txt +++ b/doc/Manual.txt @@ -1,6 +1,6 @@ LAMMPS Users Manual - + @@ -18,7 +18,7 @@

LAMMPS Documentation :c,h3 -11 Apr 2014 version :c,h4 +12 Apr 2014 version :c,h4 Version info: :h4