diff --git a/doc/src/fix_restrain.txt b/doc/src/fix_restrain.txt index 1a7c7b6ba5df11d6c9331de0963e2bc562131e42..9de63defb700fb0f713187d87f4077fd35b86e0f 100644 --- a/doc/src/fix_restrain.txt +++ b/doc/src/fix_restrain.txt @@ -187,10 +187,17 @@ 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 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 +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. diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp index 09cb8946e51d86692b201e4e5473f8839b0fc30a..6ad229fea7e70e16099b90cf9aee1e7b615e0ed3 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,6 +278,7 @@ void FixRestrain::restrain_bond(int m) if (r > 0.0) fbond = -2.0*rk/r; else fbond = 0.0; + ebond += rk*dr; energy += rk*dr; // apply force to each of 2 atoms @@ -378,6 +386,7 @@ void FixRestrain::restrain_angle(int m) dtheta = acos(c) - target[m]; tk = k * dtheta; + eangle += tk*dtheta; energy += tk*dtheta; a = -2.0 * tk * s; @@ -558,6 +567,7 @@ void FixRestrain::restrain_dihedral(int m) df1 *= -mult[m]; p += 1.0; + edihed += k * p; energy += k * p; fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm; @@ -635,3 +645,23 @@ 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) { + 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 14699ed5af38965d0e42768aec6e2fc114d3deea..4572905d46e1ae298d876fe91fb9c94f05e05d09 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,9 @@ class FixRestrain : public Fix { double *kstart,*kstop,*target; double *cos_target,*sin_target; double energy,energy_all; + double ebond,ebond_all; + double eangle,eangle_all; + double edihed,edihed_all; void restrain_bond(int); void restrain_angle(int);