From 0edad83b2516669e8f038428c38f4cab0a1fa7d7 Mon Sep 17 00:00:00 2001 From: abbatux <31911482+abbatux@users.noreply.github.com> Date: Fri, 8 Dec 2017 11:29:12 +1100 Subject: [PATCH 01/19] Update atom_vec_smd.cpp --- src/USER-SMD/atom_vec_smd.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/USER-SMD/atom_vec_smd.cpp b/src/USER-SMD/atom_vec_smd.cpp index 2a0d65b642..07d526ba31 100644 --- a/src/USER-SMD/atom_vec_smd.cpp +++ b/src/USER-SMD/atom_vec_smd.cpp @@ -1143,10 +1143,14 @@ void AtomVecSMD::pack_data(double **buf) { buf[i][7] = x[i][0]; buf[i][8] = x[i][1]; buf[i][9] = x[i][2]; + + buf[i][10] = x0[i][0]; + buf[i][11] = x0[i][1]; + buf[i][12] = x0[i][2]; - buf[i][10] = ubuf((image[i] & IMGMASK) - IMGMAX).d; - buf[i][11] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d; - buf[i][12] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d; + buf[i][13] = ubuf((image[i] & IMGMASK) - IMGMAX).d; + buf[i][14] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d; + buf[i][15] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d; } } @@ -1167,9 +1171,9 @@ void AtomVecSMD::write_data(FILE *fp, int n, double **buf) { for (int i = 0; i < n; i++) fprintf(fp, TAGINT_FORMAT - " %d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n", (tagint) ubuf(buf[i][0]).i, + " %d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n", (tagint) ubuf(buf[i][0]).i, (int) ubuf(buf[i][1]).i, (int) ubuf(buf[i][2]).i, buf[i][3], buf[i][4], buf[i][5], buf[i][6], buf[i][7], buf[i][8], - buf[i][9], (int) ubuf(buf[i][7]).i, (int) ubuf(buf[i][8]).i, (int) ubuf(buf[i][9]).i); + buf[i][9], buf[i][10], buf[i][11], buf[i][12]); } /* ---------------------------------------------------------------------- From 2f857c6eda6fc1863896d549c4638e817d120bb6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 13 Dec 2017 15:12:14 -0500 Subject: [PATCH 02/19] correct fix neigh/history/omp to properly support -DLAMMPS_BIGBIG compilation --- src/USER-OMP/fix_neigh_history_omp.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/USER-OMP/fix_neigh_history_omp.cpp b/src/USER-OMP/fix_neigh_history_omp.cpp index 1b212732de..21aed6c7c0 100644 --- a/src/USER-OMP/fix_neigh_history_omp.cpp +++ b/src/USER-OMP/fix_neigh_history_omp.cpp @@ -98,7 +98,7 @@ void FixNeighHistoryOMP::pre_exchange_onesided() // nlocal can be larger if other fixes added atoms at this pre_exchange() // clear per-thread paged data structures - + MyPage &ipg = ipage_atom[tid]; MyPage &dpg = dpage_atom[tid]; ipg.reset(); @@ -540,8 +540,8 @@ void FixNeighHistoryOMP::post_neighbor() int *allflags; double *allvalues; - MyPage &ipg = ipage_atom[tid]; - MyPage &dpg = dpage_atom[tid]; + MyPage &ipg = ipage_neigh[tid]; + MyPage &dpg = dpage_neigh[tid]; ipg.reset(); dpg.reset(); From 9a71efc5d58fd83fecd6132a6426d2fd2b6592b8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 13 Dec 2017 15:19:46 -0500 Subject: [PATCH 03/19] fix neb bugfix from Emile Maras NEB was not working fine when using multiple proc per replica and the keywords last/efirst or last/efirst/middle I have corrected this in the enclosed fix_neb.cpp I also slightly modified the nudging for this free end so that it would be applied only when the target energy is larger than the energy. Anyway if the target energy is lower than the energy, the replica should relax toward the target energy without adding any nudging. I also modified the documentation according to this change. --- doc/src/fix_neb.txt | 17 ++++++------ src/REPLICA/fix_neb.cpp | 59 ++++++++++++++++++++++------------------- 2 files changed, 39 insertions(+), 37 deletions(-) diff --git a/doc/src/fix_neb.txt b/doc/src/fix_neb.txt index 73b3e31266..be4766eae9 100644 --- a/doc/src/fix_neb.txt +++ b/doc/src/fix_neb.txt @@ -138,17 +138,17 @@ By default, no additional forces act on the first and last replicas during the NEB relaxation, so these replicas simply relax toward their respective local minima. By using the key word {end}, additional forces can be applied to the first and/or last replicas, to enable -them to relax toward a MEP while constraining their energy. +them to relax toward a MEP while constraining their energy E to the +target energy ETarget. -The interatomic force Fi for the specified replica becomes: +If ETarget>E, the interatomic force Fi for the specified replica becomes: Fi = -Grad(V) + (Grad(V) dot T' + (E-ETarget)*Kspring3) T', {when} Grad(V) dot T' < 0 Fi = -Grad(V) + (Grad(V) dot T' + (ETarget- E)*Kspring3) T', {when} Grad(V) dot T' > 0 :pre -where E is the current energy of the replica and ETarget is the target -energy. The "spring" constant on the difference in energies is the -specified {Kspring3} value. +The "spring" constant on the difference in energies is the specified +{Kspring3} value. When {estyle} is specified as {first}, the force is applied to the first replica. When {estyle} is specified as {last}, the force is @@ -183,10 +183,9 @@ After converging a NEB calculation using an {estyle} of have a larger energy than the first replica. If this is not the case, the path is probably not a MEP. -Finally, note that if the last replica converges toward a local -minimum which has a larger energy than the energy of the first -replica, a NEB calculation using an {estyle} of {last/efirst} or -{last/efirst/middle} cannot reach final convergence. +Finally, note that the last replica may never reach the target energy +if it is stuck in a local minima which has a larger energy than the +target energy. [Restart, fix_modify, output, run start/stop, minimize info:] diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 13a6ad1d62..98e5539872 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -94,26 +94,26 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : if (iarg+3 > narg) error->all(FLERR,"Illegal fix neb command"); if (strcmp(arg[iarg+1],"first") == 0) { FreeEndIni = true; - kspringIni = force->numeric(FLERR,arg[iarg+2]); + kspringIni = force->numeric(FLERR,arg[iarg+2]); } else if (strcmp(arg[iarg+1],"last") == 0) { FreeEndFinal = true; FinalAndInterWithRespToEIni = false; FreeEndFinalWithRespToEIni = false; - kspringFinal = force->numeric(FLERR,arg[iarg+2]); + kspringFinal = force->numeric(FLERR,arg[iarg+2]); } else if (strcmp(arg[iarg+1],"last/efirst") == 0) { FreeEndFinal = false; FinalAndInterWithRespToEIni = false; FreeEndFinalWithRespToEIni = true; - kspringFinal = force->numeric(FLERR,arg[iarg+2]); + kspringFinal = force->numeric(FLERR,arg[iarg+2]); } else if (strcmp(arg[iarg+1],"last/efirst/middle") == 0) { FreeEndFinal = false; FinalAndInterWithRespToEIni = true; FreeEndFinalWithRespToEIni = true; - kspringFinal = force->numeric(FLERR,arg[iarg+2]); + kspringFinal = force->numeric(FLERR,arg[iarg+2]); } else error->all(FLERR,"Illegal fix neb command"); iarg += 3; - + } else error->all(FLERR,"Illegal fix neb command"); } @@ -298,12 +298,14 @@ void FixNEB::min_post_force(int vflag) if (ireplica == 0) vIni=veng; if (FreeEndFinalWithRespToEIni) { - if (me == 0) { + if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { int procFirst; procFirst=universe->root_proc[0]; MPI_Bcast(&vIni,1,MPI_DOUBLE,procFirst,uworld); - } - if (cmode == MULTI_PROC) { + }else { + if (me == 0) + MPI_Bcast(&vIni,1,MPI_DOUBLE,0,rootworld); + MPI_Bcast(&vIni,1,MPI_DOUBLE,0,world); } } @@ -313,8 +315,8 @@ void FixNEB::min_post_force(int vflag) // if (me == 0 ) if (update->ntimestep == 0) { EIniIni = veng; - // if (cmode == MULTI_PROC) - // MPI_Bcast(&EIniIni,1,MPI_DOUBLE,0,world); + // if (cmode == MULTI_PROC) + // MPI_Bcast(&EIniIni,1,MPI_DOUBLE,0,world); } }*/ @@ -360,7 +362,7 @@ void FixNEB::min_post_force(int vflag) tangent[i][2]=delzp; tlen += tangent[i][0]*tangent[i][0] + tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2]; - dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + + dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + f[i][2]*tangent[i][2]; } } @@ -383,9 +385,9 @@ void FixNEB::min_post_force(int vflag) tangent[i][0]=delxn; tangent[i][1]=delyn; tangent[i][2]=delzn; - tlen += tangent[i][0]*tangent[i][0] + + tlen += tangent[i][0]*tangent[i][0] + tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2]; - dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + + dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + f[i][2]*tangent[i][2]; } } @@ -431,15 +433,15 @@ void FixNEB::min_post_force(int vflag) } nlen += delxn*delxn + delyn*delyn + delzn*delzn; - tlen += tangent[i][0]*tangent[i][0] + + tlen += tangent[i][0]*tangent[i][0] + tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2]; gradlen += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2]; dotpath += delxp*delxn + delyp*delyn + delzp*delzn; - dottangrad += tangent[i][0]*f[i][0] + + dottangrad += tangent[i][0]*f[i][0] + tangent[i][1]*f[i][1] + tangent[i][2]*f[i][2]; - gradnextlen += fnext[i][0]*fnext[i][0] + + gradnextlen += fnext[i][0]*fnext[i][0] + fnext[i][1]*fnext[i][1] +fnext[i][2] * fnext[i][2]; - dotgrad += f[i][0]*fnext[i][0] + f[i][1]*fnext[i][1] + + dotgrad += f[i][0]*fnext[i][0] + f[i][1]*fnext[i][1] + f[i][2]*fnext[i][2]; springF[i][0] = kspringPerp*(delxn-delxp); @@ -514,9 +516,10 @@ void FixNEB::min_post_force(int vflag) MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world); dot=dotall/tlen; - if (dot<0) prefactor = -dot - kspringFinal*(veng-EFinalIni); - else prefactor = -dot + kspringFinal*(veng-EFinalIni); - + if (veng Date: Wed, 13 Dec 2017 16:20:49 -0500 Subject: [PATCH 04/19] change pair style airebo/intel to compile with -DLAMMPS_BIGBIG --- src/USER-INTEL/pair_airebo_intel.cpp | 30 ++++++++++++++-------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/USER-INTEL/pair_airebo_intel.cpp b/src/USER-INTEL/pair_airebo_intel.cpp index ad3c97c9df..cb8168de43 100644 --- a/src/USER-INTEL/pair_airebo_intel.cpp +++ b/src/USER-INTEL/pair_airebo_intel.cpp @@ -127,7 +127,7 @@ struct KernelArgsAIREBOT { struct NeighListAIREBO neigh_rebo; PairAIREBOIntelParam params; struct AtomAIREBOT * x; /* num_all */ - int * tag; /* num_all */ + tagint * tag; /* num_all */ flt_t * nC, * nH; /* num_all */ int * map; /* num_types+1 */ struct ResultForceT * result_f; /* num_all */ @@ -399,7 +399,7 @@ void PairAIREBOIntel::eval( const int * _noalias const numneigh = list->numneigh; const int * _noalias const cnumneigh = buffers->cnumneigh(list); const int * _noalias const firstneigh = buffers->firstneigh(list); - int * const tag = atom->tag; + tagint * const tag = atom->tag; const int ntypes = atom->ntypes + 1; const int eatom = this->eflag_atom; @@ -1912,10 +1912,10 @@ template void ref_torsion(KernelArgsAIREBOT * ka) { AtomAIREBOT * x = ka->x; int * map = ka->map; - int * tag = ka->tag; + tagint * tag = ka->tag; for (int ii = ka->frebo_from_atom; ii < ka->frebo_to_atom; ii++) { int i = ii; - int itag = tag[i]; + tagint itag = tag[i]; int itype = map[x[i].w]; if (itype != 0) continue; flt_t xtmp = x[i].x; @@ -1925,7 +1925,7 @@ void ref_torsion(KernelArgsAIREBOT * ka) { int jnum = ka->neigh_rebo.num[i]; for (int jj = 0; jj < jnum; jj++) { int j = REBO_neighs_i[jj]; - int jtag = tag[j]; + tagint jtag = tag[j]; if (itag > jtag) { if (((itag+jtag) & 1) == 0) continue; @@ -2039,9 +2039,9 @@ void ref_frebo_single_interaction(KernelArgsAIREBOT * ka, int i, template inline void ref_frebo_single_atom(KernelArgsAIREBOT * ka, int i) { AtomAIREBOT * x = ka->x; - int * tag = ka->tag; + tagint * tag = ka->tag; int jj; - int itag = tag[i]; + tagint itag = tag[i]; flt_t x_i = x[i].x; flt_t y_i = x[i].y; flt_t z_i = x[i].z; @@ -2049,7 +2049,7 @@ inline void ref_frebo_single_atom(KernelArgsAIREBOT * ka, int i) { int jnum = ka->neigh_rebo.num[i]; for (jj = 0; jj < jnum; jj++) { int j = neighs[jj]; - int jtag = tag[j]; + tagint jtag = tag[j]; if (itag > jtag) { if (((itag + jtag) & 1) == 0) continue; @@ -2165,9 +2165,9 @@ template void ref_lennard_jones_single_atom(KernelArgsAIREBOT * ka, int i, int morseflag) { AtomAIREBOT * x = ka->x; - int * tag = ka->tag; + tagint * tag = ka->tag; int jj; - int itag = tag[i]; + tagint itag = tag[i]; int * neighs = ka->neigh_lmp.entries + ka->neigh_lmp.offset[i]; int jnum = ka->neigh_lmp.num_half[i]; for (jj = 0; jj < jnum; jj++) { @@ -3640,7 +3640,7 @@ static void aut_frebo_batch_of_kind(KernelArgsAIREBOT * ka, int * i_buf, int * j_buf) { { // jump-scope for exceed_limits AtomAIREBOT * x = ka->x; - int * tag = ka->tag; + tagint * tag = ka->tag; int * map = ka->map; ResultForceT * result_f = ka->result_f; flt_t rcminij = ka->params.rcmin[itype][jtype]; @@ -3847,13 +3847,13 @@ exceed_limits: */ static void aut_frebo(KernelArgsAIREBOT * ka, int torflag) { AtomAIREBOT * _noalias x = ka->x; - int * _noalias tag = ka->tag; + tagint * _noalias tag = ka->tag; int * _noalias map = ka->map; int i_buf[2][2][fvec::VL]; int j_buf[2][2][fvec::VL]; int n_buf[2][2] = {0}; for (int i = ka->frebo_from_atom; i < ka->frebo_to_atom; i++) { - int itag = tag[i]; + tagint itag = tag[i]; int itype = map[x[i].w]; flt_t x_i = x[i].x; flt_t y_i = x[i].y; @@ -3862,7 +3862,7 @@ static void aut_frebo(KernelArgsAIREBOT * ka, int torflag) { int jnum = ka->neigh_rebo.num[i]; for (int jj = 0; jj < jnum; jj++) { int j = neighs[jj]; - int jtag = tag[j]; + tagint jtag = tag[j]; if (itag > jtag) { if (((itag + jtag) & 1) == 0) continue; @@ -4490,7 +4490,7 @@ exceed_limits: template static void aut_lennard_jones(KernelArgsAIREBOT * ka) { AtomAIREBOT * x = ka->x; - int * tag = ka->tag; + tagint * tag = ka->tag; int * map = ka->map; ResultForceT * result_f = ka->result_f; ivec c_i1 = ivec::set1(1); From 8a36cdc6bc6a1b537efdbaf1ef1ceb889a1fe068 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 13 Dec 2017 18:42:24 -0500 Subject: [PATCH 05/19] correct velocity output for write_data of atom style smd --- src/USER-SMD/atom_vec_smd.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/USER-SMD/atom_vec_smd.cpp b/src/USER-SMD/atom_vec_smd.cpp index 07d526ba31..8bffaf8a09 100644 --- a/src/USER-SMD/atom_vec_smd.cpp +++ b/src/USER-SMD/atom_vec_smd.cpp @@ -1215,8 +1215,7 @@ int AtomVecSMD::pack_vel_hybrid(int i, double *buf) { void AtomVecSMD::write_vel(FILE *fp, int n, double **buf) { for (int i = 0; i < n; i++) fprintf(fp, TAGINT_FORMAT - " %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n", (tagint) ubuf(buf[i][0]).i, buf[i][1], buf[i][2], buf[i][3], - buf[i][4], buf[i][5], buf[i][6]); + " %-1.16e %-1.16e %-1.16e\n", (tagint) ubuf(buf[i][0]).i, buf[i][1], buf[i][2], buf[i][3]); } /* ---------------------------------------------------------------------- From 10fa54b2fde854df152d2326ea3a44a5f6e4e8ae Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 13 Dec 2017 18:44:31 -0500 Subject: [PATCH 06/19] correct error messages. the atom style is called "smd" not "tlsph" --- src/USER-SMD/atom_vec_smd.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/USER-SMD/atom_vec_smd.cpp b/src/USER-SMD/atom_vec_smd.cpp index 8bffaf8a09..d2d72a855a 100644 --- a/src/USER-SMD/atom_vec_smd.cpp +++ b/src/USER-SMD/atom_vec_smd.cpp @@ -1099,7 +1099,7 @@ void AtomVecSMD::data_atom(double *coord, imageint imagetmp, char **values) { ------------------------------------------------------------------------- */ int AtomVecSMD::data_atom_hybrid(int nlocal, char **values) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style tlsph"); + error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); return -1; } @@ -1121,7 +1121,7 @@ void AtomVecSMD::data_vel(int m, char **values) { ------------------------------------------------------------------------- */ int AtomVecSMD::data_vel_hybrid(int m, char **values) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style tlsph"); + error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); return 0; } @@ -1159,7 +1159,7 @@ void AtomVecSMD::pack_data(double **buf) { ------------------------------------------------------------------------- */ int AtomVecSMD::pack_data_hybrid(int i, double *buf) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style tlsph"); + error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); return -1; } @@ -1181,7 +1181,7 @@ void AtomVecSMD::write_data(FILE *fp, int n, double **buf) { ------------------------------------------------------------------------- */ int AtomVecSMD::write_data_hybrid(FILE *fp, double *buf) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style tlsph"); + error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); return -1; } @@ -1204,7 +1204,7 @@ void AtomVecSMD::pack_vel(double **buf) { ------------------------------------------------------------------------- */ int AtomVecSMD::pack_vel_hybrid(int i, double *buf) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style tlsph"); + error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); return 0; } @@ -1223,7 +1223,7 @@ void AtomVecSMD::write_vel(FILE *fp, int n, double **buf) { ------------------------------------------------------------------------- */ int AtomVecSMD::write_vel_hybrid(FILE *fp, double *buf) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style tlsph"); + error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); return 3; } From 56e633a2cf417fc2a6da10b1fc39c1d98571b23d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 19 Dec 2017 14:54:35 -0500 Subject: [PATCH 07/19] error out on trying to compile USER-INTEL with -DLAMMPS_BIGBIG --- src/USER-INTEL/fix_intel.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/USER-INTEL/fix_intel.cpp b/src/USER-INTEL/fix_intel.cpp index 3e36c8f7a9..cac5e0b49c 100644 --- a/src/USER-INTEL/fix_intel.cpp +++ b/src/USER-INTEL/fix_intel.cpp @@ -34,6 +34,10 @@ #include #include +#ifdef LAMMPS_BIGBIG +#error "The USER-INTEL package is not compatible with -DLAMMPS_BIGBIG" +#endif + #ifdef _LMP_INTEL_OFFLOAD #ifndef INTEL_OFFLOAD_NOAFFINITY #include From 18acc6ae47b1378ab307eff8cffd7b388e8fd77d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 19 Dec 2017 15:01:11 -0500 Subject: [PATCH 08/19] remove some dead code --- src/USER-INTEL/pair_airebo_intel.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/USER-INTEL/pair_airebo_intel.cpp b/src/USER-INTEL/pair_airebo_intel.cpp index cb8168de43..54c5055812 100644 --- a/src/USER-INTEL/pair_airebo_intel.cpp +++ b/src/USER-INTEL/pair_airebo_intel.cpp @@ -2165,9 +2165,7 @@ template void ref_lennard_jones_single_atom(KernelArgsAIREBOT * ka, int i, int morseflag) { AtomAIREBOT * x = ka->x; - tagint * tag = ka->tag; int jj; - tagint itag = tag[i]; int * neighs = ka->neigh_lmp.entries + ka->neigh_lmp.offset[i]; int jnum = ka->neigh_lmp.num_half[i]; for (jj = 0; jj < jnum; jj++) { @@ -3640,7 +3638,6 @@ static void aut_frebo_batch_of_kind(KernelArgsAIREBOT * ka, int * i_buf, int * j_buf) { { // jump-scope for exceed_limits AtomAIREBOT * x = ka->x; - tagint * tag = ka->tag; int * map = ka->map; ResultForceT * result_f = ka->result_f; flt_t rcminij = ka->params.rcmin[itype][jtype]; @@ -4490,7 +4487,6 @@ exceed_limits: template static void aut_lennard_jones(KernelArgsAIREBOT * ka) { AtomAIREBOT * x = ka->x; - tagint * tag = ka->tag; int * map = ka->map; ResultForceT * result_f = ka->result_f; ivec c_i1 = ivec::set1(1); @@ -4521,7 +4517,6 @@ static void aut_lennard_jones(KernelArgsAIREBOT * ka) { int num_bo[2][2] = {0}; for (int i = ka->frebo_from_atom; i < ka->frebo_to_atom; i++) { - ivec itag_bc = ivec::set1(tag[i]); int itype = map[x[i].w]; fvec x_i = fvec::set1(x[i].x); fvec y_i = fvec::set1(x[i].y); From d4f45f4f85384dd57c7bc7275307830f60a20194 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 24 Dec 2017 17:45:58 -0500 Subject: [PATCH 09/19] correct set command example in fix property/atom --- doc/src/fix_property_atom.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_property_atom.txt b/doc/src/fix_property_atom.txt index aef534ad25..5a8601cf00 100644 --- a/doc/src/fix_property_atom.txt +++ b/doc/src/fix_property_atom.txt @@ -155,7 +155,7 @@ these commands could be used: fix prop all property/atom mol variable cluster atom ((id-1)/10)+1 -set id * mol v_cluster :pre +set atom * mol v_cluster :pre The "atom-style variable"_variable.html will create values for atoms with IDs 31,32,33,...40 that are 4.0,4.1,4.2,...,4.9. When the From a9e9a2046b49175d6db79a0f2439f551ed77122b Mon Sep 17 00:00:00 2001 From: Stefan Paquay Date: Mon, 25 Dec 2017 13:03:18 +0100 Subject: [PATCH 10/19] Fixes/clarifies the fix_property_atom docs. --- doc/src/fix_property_atom.txt | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_property_atom.txt b/doc/src/fix_property_atom.txt index aef534ad25..58a3c15cbf 100644 --- a/doc/src/fix_property_atom.txt +++ b/doc/src/fix_property_atom.txt @@ -111,6 +111,10 @@ need to communicate their new values to/from ghost atoms, an operation that can be invoked from within a "pair style"_pair_style.html or "fix"_fix.html or "compute"_compute.html that you write. +NOTE: If this fix is defined _after_ the simulation box is created, +a 'run 0' command should be issued to properly initialize the storage +created by this fix. + :line This fix is one of a small number that can be defined in an input @@ -155,7 +159,7 @@ these commands could be used: fix prop all property/atom mol variable cluster atom ((id-1)/10)+1 -set id * mol v_cluster :pre +set atom * mol v_cluster :pre The "atom-style variable"_variable.html will create values for atoms with IDs 31,32,33,...40 that are 4.0,4.1,4.2,...,4.9. When the From c333401e7250a9b787f27ef2291135e06d8c589a Mon Sep 17 00:00:00 2001 From: Stefan Paquay Date: Mon, 25 Dec 2017 13:27:14 +0100 Subject: [PATCH 11/19] Use bold font instead of underscores for emphasis. --- doc/src/fix_property_atom.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_property_atom.txt b/doc/src/fix_property_atom.txt index 58a3c15cbf..10bb89c94c 100644 --- a/doc/src/fix_property_atom.txt +++ b/doc/src/fix_property_atom.txt @@ -111,7 +111,7 @@ need to communicate their new values to/from ghost atoms, an operation that can be invoked from within a "pair style"_pair_style.html or "fix"_fix.html or "compute"_compute.html that you write. -NOTE: If this fix is defined _after_ the simulation box is created, +NOTE: If this fix is defined [after] the simulation box is created, a 'run 0' command should be issued to properly initialize the storage created by this fix. From 2896df2140b6993beed0f31956dd7b7356ab6104 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 4 Jan 2018 17:16:31 -0500 Subject: [PATCH 12/19] avoid that occasional neighbor lists requested from commands linger around for too long and thus cause segementation faults --- src/create_bonds.cpp | 4 ++++ src/delete_atoms.cpp | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/src/create_bonds.cpp b/src/create_bonds.cpp index c62c6df91c..2f6624f011 100644 --- a/src/create_bonds.cpp +++ b/src/create_bonds.cpp @@ -306,6 +306,10 @@ void CreateBonds::many() nadd_bonds,atom->nbonds); } } + // trigger clearing the list of available neighbor list requests + // and a full rebuild of them during the next run setup. + // otherwise the request from this command may linger around. + neighbor->init(); } /* ---------------------------------------------------------------------- */ diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp index 489c5bf5d5..a4219dbecb 100644 --- a/src/delete_atoms.cpp +++ b/src/delete_atoms.cpp @@ -407,6 +407,10 @@ void DeleteAtoms::delete_overlap(int narg, char **arg) break; } } + // trigger clearing the list of available neighbor list requests + // and a full rebuild of them during the next run setup. + // otherwise the request from this command may linger around. + neighbor->init(); } /* ---------------------------------------------------------------------- From 75f1a4f3f0ef30989c6ce3b5932cd7a4c251fa58 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 4 Jan 2018 20:43:53 -0500 Subject: [PATCH 13/19] use filelink.o consistently, since filelink does not work with some compilers, e.g. nvcc --- doc/src/Section_packages.txt | 2 +- lib/latte/README | 15 ++++++++------- src/LATTE/README | 2 +- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/doc/src/Section_packages.txt b/doc/src/Section_packages.txt index 912d371cd9..98028887ea 100644 --- a/doc/src/Section_packages.txt +++ b/doc/src/Section_packages.txt @@ -732,7 +732,7 @@ make lib-latte args="-b -m gfortran" # download and build in lib/latte and :pre Note that 3 symbolic (soft) links, "includelink" and "liblink" and -"filelink", are created in lib/latte to point into the LATTE home dir. +"filelink.o", are created in lib/latte to point into the LATTE home dir. When LAMMPS builds in src it will use these links. You should also check that the Makefile.lammps file you create is apporpriate for the compiler you use on your system to build LATTE. diff --git a/lib/latte/README b/lib/latte/README index bfa91e62cd..bdf67dd59c 100644 --- a/lib/latte/README +++ b/lib/latte/README @@ -28,14 +28,15 @@ Instructions: and finally, build the code via the make command % make -3. Create three soft links in this dir (lib/latte) +3. Create three symbolic links in this dir (lib/latte) E.g if you built LATTE in this dir: % ln -s ./LATTE-master/src includelink % ln -s ./LATTE-master liblink - % ln -s ./LATTE-master/src/latte_c_bind.o filelink + % ln -s ./LATTE-master/src/latte_c_bind.o filelink.o 4. Choose a Makefile.lammps.* file appropriate for your compiler - (GNU gfortran or Intel ifort) and copy it to Makefile.lammps. + (GNU gfortran with external BLAS, GNU gfortran with local liblinalg, + or Intel ifort with MKL) and copy it to Makefile.lammps. Note that you may need to edit Makefile.lammps for paths and compiler options appropriate to your system. @@ -49,7 +50,7 @@ with the LATTE package installed: % make g++ (or whatever target you wish) Note that if you download and unpack a new LAMMPS tarball, the -"includelink" and "liblink" and "filelink" files will be lost and you -will need to re-create them (step 3). If you built LATTE in this -directory (as opposed to somewhere else on your system), you will also -need to repeat steps 1,2,4. +"includelink" and "liblink" and "filelink.o" symbolic links will be +lost and you will need to re-create them (step 3). If you built LATTE +in this directory (as opposed to somewhere else on your system), you +will also need to repeat steps 1,2,4. diff --git a/src/LATTE/README b/src/LATTE/README index 2103074ad4..d7379f2461 100644 --- a/src/LATTE/README +++ b/src/LATTE/README @@ -17,7 +17,7 @@ Once you have successfully built LAMMPS with this package and the LATTE library you can test it using an input file from the examples latte dir, e.g. -lmp_serial < lammps/examples/latte/in.latte.water +lmp_serial -in lammps/examples/latte/in.latte.water This pair style was written in collaboration with the LATTE developers. From 5ecc3ce366c32128d8277e205b783ae91fb9242b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 4 Jan 2018 20:44:40 -0500 Subject: [PATCH 14/19] avoid division by zero when trying to run PPPM on a system without atoms --- src/KSPACE/pppm.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 6add8b58b7..28dda4abfc 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -1268,6 +1268,7 @@ double PPPM::compute_qopt() double PPPM::estimate_ik_error(double h, double prd, bigint natoms) { double sum = 0.0; + if (natoms == 0) return 0.0; for (int m = 0; m < order; m++) sum += acons[order][m] * pow(h*g_ewald,2.0*m); double value = q2 * pow(h*g_ewald,(double)order) * From 91993b236d52dbc9e75c9f7c45cf1438d14c51cc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 5 Jan 2018 19:52:51 -0500 Subject: [PATCH 15/19] avoid division by zero in PPPM for empty and uncharged systems. require kspace_modify gewald --- src/KOKKOS/pppm_kokkos.cpp | 3 +++ src/KSPACE/pppm.cpp | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/KOKKOS/pppm_kokkos.cpp b/src/KOKKOS/pppm_kokkos.cpp index cf6e2814c0..03de42c68c 100644 --- a/src/KOKKOS/pppm_kokkos.cpp +++ b/src/KOKKOS/pppm_kokkos.cpp @@ -969,6 +969,8 @@ void PPPMKokkos::set_grid_global() if (!gewaldflag) { if (accuracy <= 0.0) error->all(FLERR,"KSpace accuracy must be > 0"); + if (q2 == 0.0) + error->all(FLERR,"Must use 'kspace_modify gewald' for uncharged system"); g_ewald = accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*q2); if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/cutoff; else g_ewald = sqrt(-log(g_ewald)) / cutoff; @@ -1179,6 +1181,7 @@ double PPPMKokkos::final_accuracy() double yprd = domain->yprd; double zprd = domain->zprd; bigint natoms = atomKK->natoms; + if (natoms == 0) natoms = 1; // avoid division by zero double df_kspace = compute_df_kspace(); double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd); diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 28dda4abfc..53ab2e5a9d 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -1004,6 +1004,8 @@ void PPPM::set_grid_global() if (!gewaldflag) { if (accuracy <= 0.0) error->all(FLERR,"KSpace accuracy must be > 0"); + if (q2 == 0.0) + error->all(FLERR,"Must use kspace_modify gewald for uncharged system"); g_ewald = accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*q2); if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/cutoff; else g_ewald = sqrt(-log(g_ewald)) / cutoff; @@ -1346,6 +1348,7 @@ double PPPM::final_accuracy() double yprd = domain->yprd; double zprd = domain->zprd; bigint natoms = atom->natoms; + if (natoms == 0) natoms = 1; // avoid division by zero double df_kspace = compute_df_kspace(); double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd); From 6c058fb56c4240856a2545ada41b1cd253f01a83 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 5 Jan 2018 20:14:25 -0500 Subject: [PATCH 16/19] avoid division by zero in ewald for empty and uncharged systems. require kspace_modify gewald --- src/KSPACE/ewald.cpp | 3 +++ src/KSPACE/ewald_disp.cpp | 5 +++++ 2 files changed, 8 insertions(+) diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index d41443f02a..c7d242031f 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -154,6 +154,8 @@ void Ewald::init() if (!gewaldflag) { if (accuracy <= 0.0) error->all(FLERR,"KSpace accuracy must be > 0"); + if (q2 == 0.0) + error->all(FLERR,"Must use 'kspace_modify gewald' for uncharged system"); g_ewald = accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / (2.0*q2); if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/cutoff; else g_ewald = sqrt(-log(g_ewald)) / cutoff; @@ -340,6 +342,7 @@ void Ewald::setup() double Ewald::rms(int km, double prd, bigint natoms, double q2) { + if (natoms == 0) natoms = 1; // avoid division by zero double value = 2.0*q2*g_ewald/prd * sqrt(1.0/(MY_PI*km*natoms)) * exp(-MY_PI*MY_PI*km*km/(g_ewald*g_ewald*prd*prd)); diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp index 85e3da921b..e0424b0d93 100644 --- a/src/KSPACE/ewald_disp.cpp +++ b/src/KSPACE/ewald_disp.cpp @@ -193,6 +193,10 @@ void EwaldDisp::init() if (!gewaldflag) { if (function[0]) { + if (accuracy <= 0.0) + error->all(FLERR,"KSpace accuracy must be > 0"); + if (q2 == 0.0) + error->all(FLERR,"Must use 'kspace_modify gewald' for uncharged system"); g_ewald = accuracy*sqrt(natoms*(*cutoff)*shape_det(domain->h)) / (2.0*q2); if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/(*cutoff); else g_ewald = sqrt(-log(g_ewald)) / (*cutoff); @@ -298,6 +302,7 @@ double EwaldDisp::rms(int km, double prd, bigint natoms, double q2, double b2, double M2) { double value = 0.0; + if (natoms == 0) natoms = 1; // avoid division by zero // Coulombic From 6070182f06bb8dee5cf87485dc34707d482fe53a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 6 Jan 2018 11:03:48 -0500 Subject: [PATCH 17/19] correctly account for individually added bonds, angles, and dihedrals --- src/create_bonds.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/create_bonds.cpp b/src/create_bonds.cpp index 2f6624f011..ff62e86787 100644 --- a/src/create_bonds.cpp +++ b/src/create_bonds.cpp @@ -342,6 +342,7 @@ void CreateBonds::single_bond() bond_atom[m][num_bond[m]] = batom2; num_bond[m]++; } + ++atom->nbonds; if (force->newton_bond) return; @@ -389,6 +390,7 @@ void CreateBonds::single_angle() angle_atom3[m][num_angle[m]] = aatom3; num_angle[m]++; } + ++atom->nangles; if (force->newton_bond) return; @@ -452,6 +454,7 @@ void CreateBonds::single_dihedral() dihedral_atom4[m][num_dihedral[m]] = datom4; num_dihedral[m]++; } + ++atom->ndihedrals; if (force->newton_bond) return; From 46217db8a5e6d971e23f9261c7ded3b670b513d6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 6 Jan 2018 12:56:42 -0500 Subject: [PATCH 18/19] make python functions examples more complete by adding exception handling and initializing variables --- examples/python/funcs.py | 26 +++++++++++++++----------- examples/python/in.python | 1 + 2 files changed, 16 insertions(+), 11 deletions(-) diff --git a/examples/python/funcs.py b/examples/python/funcs.py index f38aca53f2..71865f62f4 100644 --- a/examples/python/funcs.py +++ b/examples/python/funcs.py @@ -8,17 +8,21 @@ def loop(N,cut0,thresh,lmpptr): from lammps import lammps lmp = lammps(ptr=lmpptr) natoms = lmp.get_natoms() - - for i in range(N): - cut = cut0 + i*0.1 - lmp.set_variable("cut",cut) # set a variable in LAMMPS + try: + for i in range(N): + cut = cut0 + i*0.1 - lmp.command("pair_style lj/cut ${cut}") # LAMMPS command - #lmp.command("pair_style lj/cut %d" % cut) # LAMMPS command option + lmp.set_variable("cut",cut) # set a variable in LAMMPS + + lmp.command("pair_style lj/cut ${cut}") # LAMMPS command + #lmp.command("pair_style lj/cut %d" % cut) # LAMMPS command option + + lmp.command("pair_coeff * * 1.0 1.0") # ditto + lmp.command("run 10") # ditto + pe = lmp.extract_compute("thermo_pe",0,0) # extract total PE from LAMMPS + print("PE",pe/natoms,thresh) + if pe/natoms < thresh: return + except Exception as e: + print("LOOP error:", e) - lmp.command("pair_coeff * * 1.0 1.0") # ditto - lmp.command("run 10") # ditto - pe = lmp.extract_compute("thermo_pe",0,0) # extract total PE from LAMMPS - print("PE",pe/natoms,thresh) - if pe/natoms < thresh: return diff --git a/examples/python/in.python b/examples/python/in.python index c5aa504d43..9372c684d6 100644 --- a/examples/python/in.python +++ b/examples/python/in.python @@ -28,6 +28,7 @@ python simple here """ from __future__ import print_function def simple(): + foo = 0 print("Inside simple function") try: foo += 1 From 3af389e6cff2ff41d8c1ad6367adc44db1e82dfd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 11 Dec 2017 08:08:49 -0500 Subject: [PATCH 19/19] Revert "no need to include library.o in the LAMMPS executable" This reverts commit 4a3a6b445597dc7d320fb59b910170bb5688377a. As it turns out, when using the LAMMPS python wrapper from inside code using the PYTHON package, the library symbols *are* needed. Thanks for Richard Berger (@rbberger) for pointing this out. --- src/Makefile | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Makefile b/src/Makefile index c3c84b3665..e0f0db77fe 100644 --- a/src/Makefile +++ b/src/Makefile @@ -17,12 +17,12 @@ SHLINK = liblammps.so OBJDIR = Obj_$@ OBJSHDIR = Obj_shared_$@ -SRC = $(filter-out library.cpp,$(wildcard *.cpp)) -INC = $(filter-out library.h,$(wildcard *.h)) +SRC = $(wildcard *.cpp) +INC = $(wildcard *.h) OBJ = $(SRC:.cpp=.o) -SRCLIB = $(filter-out main.cpp,$(SRC)) library.cpp -OBJLIB = $(filter-out main.o,$(OBJ)) library.o +SRCLIB = $(filter-out main.cpp,$(SRC)) +OBJLIB = $(filter-out main.o,$(OBJ)) # Command-line options for mode: exe (default), shexe, lib, shlib @@ -176,7 +176,7 @@ help: then cp Makefile.package.settings.empty Makefile.package.settings; fi @cp Makefile.package Makefile.package.settings $(objdir) @cd $(objdir); rm -f .depend; \ - $(MAKE) $(MFLAGS) "SRC = $(SRC) library.cpp" "INC = $(INC) library.h" depend || : + $(MAKE) $(MFLAGS) "SRC = $(SRC)" "INC = $(INC)" depend || : ifeq ($(mode),exe) @cd $(objdir); \ $(MAKE) $(MFLAGS) "OBJ = $(OBJ)" "INC = $(INC)" "SHFLAGS =" \