From 69605295942797b0a9d1474f4a563b34ae845a5f Mon Sep 17 00:00:00 2001 From: robeme Date: Tue, 8 May 2018 10:02:10 +0200 Subject: [PATCH 1/5] set appropriate flags for computing a vector in constructor of fix restrain --- src/fix_restrain.cpp | 28 +++++++++++++++++++++++++--- src/fix_restrain.h | 2 ++ 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index 09cb8946e5..80775e0283 100644 --- a/src/fix_restrain.cpp +++ b/src/fix_restrain.cpp @@ -53,6 +53,9 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) : scalar_flag = 1; global_freq = 1; extscalar = 1; + vector_flag = 1; + size_vector = 3; + extvector = 1; respa_level_support = 1; ilevel_respa = 0; @@ -188,6 +191,10 @@ void FixRestrain::min_setup(int vflag) void FixRestrain::post_force(int vflag) { energy = 0.0; + + ebond = 0.0; + eangle = 0.0; + edihed = 0.0; for (int m = 0; m < nrestrain; m++) if (rstyle[m] == BOND) restrain_bond(m); @@ -271,7 +278,8 @@ void FixRestrain::restrain_bond(int m) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; - energy += rk*dr; + ebond += rk*dr; + energy += ebond; // apply force to each of 2 atoms @@ -378,7 +386,8 @@ void FixRestrain::restrain_angle(int m) dtheta = acos(c) - target[m]; tk = k * dtheta; - energy += tk*dtheta; + eangle += tk*dtheta; + energy += eangle; a = -2.0 * tk * s; a11 = a*c / rsq1; @@ -558,7 +567,8 @@ void FixRestrain::restrain_dihedral(int m) df1 *= -mult[m]; p += 1.0; - energy += k * p; + edihed += k * p; + energy += edihed; fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm; @@ -635,3 +645,15 @@ double FixRestrain::compute_scalar() MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); return energy_all; } + +/* ---------------------------------------------------------------------- + return individual energy contributions +------------------------------------------------------------------------- */ + +double FixRestrain::compute_vector(int n) +{ + if (n == 0) return ebond; + if (n == 1) return eangle; + if (n == 2) return edihed; + return 0.0; +} diff --git a/src/fix_restrain.h b/src/fix_restrain.h index 14699ed5af..c956e40982 100644 --- a/src/fix_restrain.h +++ b/src/fix_restrain.h @@ -36,6 +36,7 @@ class FixRestrain : public Fix { void post_force_respa(int, int, int); void min_post_force(int); double compute_scalar(); + double compute_vector(int); private: int ilevel_respa; @@ -46,6 +47,7 @@ class FixRestrain : public Fix { double *kstart,*kstop,*target; double *cos_target,*sin_target; double energy,energy_all; + double ebond,eangle,edihed; void restrain_bond(int); void restrain_angle(int); From 53fc9f1f0ffb6aeada893aff4d50e060b167891e Mon Sep 17 00:00:00 2001 From: robeme Date: Tue, 8 May 2018 10:11:00 +0200 Subject: [PATCH 2/5] updated doc for new fix restrain output --- doc/src/fix_restrain.txt | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_restrain.txt b/doc/src/fix_restrain.txt index 1a7c7b6ba5..6d3e20bdea 100644 --- a/doc/src/fix_restrain.txt +++ b/doc/src/fix_restrain.txt @@ -187,10 +187,16 @@ added forces to be included in the total potential energy of the system (the quantity being minimized), you MUST enable the "fix_modify"_fix_modify.html {energy} option for this fix. -This fix computes a global scalar, which can be accessed by various -"output commands"_Section_howto.html#howto_15. The scalar is the -potential energy for all the restraints as discussed above. The scalar -value calculated by this fix is "extensive". +This fix computes a global scalar and vector of length 3, which can be +accessed by various "output commands"_Section_howto.html#howto_15. The +scalar is the potential energy for all the restraints as discussed above. +The vector values are the following global quantities: + +1 = bond energy +2 = angle energy +3 = dihedral energy :ul + +The scalar and vector values calculated by this fix are "extensive". No parameter of this fix can be used with the {start/stop} keywords of the "run"_run.html command. From 16697fc4cb7d72de38442d552a9026586f2cf526 Mon Sep 17 00:00:00 2001 From: robeme Date: Tue, 8 May 2018 10:46:41 +0200 Subject: [PATCH 3/5] forgot to reduce energies --- src/fix_restrain.cpp | 16 ++++++++++++---- src/fix_restrain.h | 4 +++- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index 80775e0283..4c0d3e4a29 100644 --- a/src/fix_restrain.cpp +++ b/src/fix_restrain.cpp @@ -652,8 +652,16 @@ double FixRestrain::compute_scalar() double FixRestrain::compute_vector(int n) { - if (n == 0) return ebond; - if (n == 1) return eangle; - if (n == 2) return edihed; - return 0.0; + if (n == 0) { + MPI_Allreduce(&ebond,&ebond_all,1,MPI_DOUBLE,MPI_SUM,world); + return ebond_all; + } else if (n == 1) { + MPI_Allreduce(&eangle,&eangle_all,1,MPI_DOUBLE,MPI_SUM,world); + return eangle_all; + } else if (n == 2) { + MPI_Allreduce(&edihed,&edihed_all,1,MPI_DOUBLE,MPI_SUM,world); + return edihed_all; + } else { + return 0.0; + } } diff --git a/src/fix_restrain.h b/src/fix_restrain.h index c956e40982..4572905d46 100644 --- a/src/fix_restrain.h +++ b/src/fix_restrain.h @@ -47,7 +47,9 @@ class FixRestrain : public Fix { double *kstart,*kstop,*target; double *cos_target,*sin_target; double energy,energy_all; - double ebond,eangle,edihed; + double ebond,ebond_all; + double eangle,eangle_all; + double edihed,edihed_all; void restrain_bond(int); void restrain_angle(int); From 671876efefe8c6343f542657f2620a1a60709b7f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 8 May 2018 15:33:42 -0400 Subject: [PATCH 4/5] Make the description of global properties more explicit --- doc/src/fix_restrain.txt | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_restrain.txt b/doc/src/fix_restrain.txt index 6d3e20bdea..9de63defb7 100644 --- a/doc/src/fix_restrain.txt +++ b/doc/src/fix_restrain.txt @@ -187,10 +187,11 @@ added forces to be included in the total potential energy of the system (the quantity being minimized), you MUST enable the "fix_modify"_fix_modify.html {energy} option for this fix. -This fix computes a global scalar and vector of length 3, which can be -accessed by various "output commands"_Section_howto.html#howto_15. The -scalar is the potential energy for all the restraints as discussed above. -The vector values are the following global quantities: +This fix computes a global scalar and a global vector of length 3, which +can be accessed by various "output commands"_Section_howto.html#howto_15. +The scalar is the total potential energy for {all} the restraints as +discussed above. The vector values are the sum of contributions to the +following individual categories: 1 = bond energy 2 = angle energy From b4e5828a6085dcc4f54d7595b914ba6e14224d7d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 8 May 2018 15:36:32 -0400 Subject: [PATCH 5/5] Update fix_restrain.cpp --- src/fix_restrain.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index 4c0d3e4a29..6ad229fea7 100644 --- a/src/fix_restrain.cpp +++ b/src/fix_restrain.cpp @@ -279,7 +279,7 @@ void FixRestrain::restrain_bond(int m) else fbond = 0.0; ebond += rk*dr; - energy += ebond; + energy += rk*dr; // apply force to each of 2 atoms @@ -387,7 +387,7 @@ void FixRestrain::restrain_angle(int m) tk = k * dtheta; eangle += tk*dtheta; - energy += eangle; + energy += tk*dtheta; a = -2.0 * tk * s; a11 = a*c / rsq1; @@ -568,7 +568,7 @@ void FixRestrain::restrain_dihedral(int m) p += 1.0; edihed += k * p; - energy += edihed; + energy += k * p; fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm;