From 2650e1ea69479af5652660108cda9a6fcf6cc575 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 20 Jul 2022 08:19:29 -0400 Subject: [PATCH 1/7] Make angle energy calc exactly what is in sander --- src/Energy.cpp | 14 ++++++++++---- src/Energy/Ene_Angle.h | 43 ++++++++++++++++++++++++++++++++++++++++++ src/cpptrajdepend | 2 +- 3 files changed, 54 insertions(+), 5 deletions(-) create mode 100644 src/Energy/Ene_Angle.h diff --git a/src/Energy.cpp b/src/Energy.cpp index 4e71df345f..2adb31b2b4 100644 --- a/src/Energy.cpp +++ b/src/Energy.cpp @@ -7,6 +7,7 @@ #include "Topology.h" #include "CharMask.h" #include "ExclusionArray.h" +#include "Energy/Ene_Angle.h" const double Energy_Amber::QFAC = Constants::ELECTOAMBER * Constants::ELECTOAMBER; @@ -43,8 +44,10 @@ double Energy_Amber::CalcBondEnergy(Frame const& fIn, BondArray const& Bonds, double rdiff = r - bp.Req(); double ene = bp.Rk() * (rdiff * rdiff); Ebond += ene; + mprintf("EBOND %4li %4i -- %4i: k= %12.5f x0= %12.5f r= %12.5f E= %12.5f\n", + b - Bonds.begin(), b->A1()+1, b->A2()+1, bp.Rk(), bp.Req(), r, ene); # ifdef DEBUG_ENERGY - mprintf("\tBond %4u %4i -- %4i: k= %12.5f x0= %12.5f r= %12.5f E= %12.5e\n", + mprintf("\tBond %4li %4i -- %4i: k= %12.5f x0= %12.5f r= %12.5f E= %12.5e\n", b - Bonds.begin(), b->A1()+1, b->A2()+1, bp.Rk(), bp.Req(), r, ene); # endif } @@ -81,9 +84,12 @@ double Energy_Amber::CalcAngleEnergy(Frame const& fIn, AngleArray const& Angles, continue; } AngleParmType const& ap = APA[apidx]; - double theta = CalcAngle(fIn.XYZ(a->A1()), fIn.XYZ(a->A2()), fIn.XYZ(a->A3())); - double tdiff = theta - ap.Teq(); - double ene = ap.Tk() * (tdiff * tdiff); + double ene = Cpptraj::Energy::Ene_Angle( + fIn.XYZ(a->A1()), fIn.XYZ(a->A2()), fIn.XYZ(a->A3()), ap.Teq(), ap.Tk()); + +// double theta = CalcAngle(fIn.XYZ(a->A1()), fIn.XYZ(a->A2()), fIn.XYZ(a->A3())); +// double tdiff = theta - ap.Teq(); +// double ene = ap.Tk() * (tdiff * tdiff); Eangle += ene; # ifdef DEBUG_ENERGY mprintf("\tAngle %4u %4i -- %4i -- %4i: k= %12.5f x0= %12.5f t= %12.5f E= %12.5e\n", diff --git a/src/Energy/Ene_Angle.h b/src/Energy/Ene_Angle.h new file mode 100644 index 0000000000..c36453ade3 --- /dev/null +++ b/src/Energy/Ene_Angle.h @@ -0,0 +1,43 @@ +#ifndef INC_ENERGY_ENE_ANGLE_H +#define INC_ENERGY_ENE_ANGLE_H +namespace Cpptraj { +namespace Energy { + +template +T Ene_Angle(T const* xyz1, T const* xyz2, T const* xyz3, T const& teq, T const& tk) +{ + static const T pt999 = 0.9990; + T ij[3], kj[3]; + + ij[0] = xyz1[0] - xyz2[0]; + ij[1] = xyz1[1] - xyz2[1]; + ij[2] = xyz1[2] - xyz2[2]; + kj[0] = xyz3[0] - xyz2[0]; + kj[1] = xyz3[1] - xyz2[1]; + kj[2] = xyz3[2] - xyz2[2]; + + T rij = ij[0]*ij[0] + ij[1]*ij[1] + ij[2]*ij[2]; + T rkj = kj[0]*kj[0] + kj[1]*kj[1] + kj[2]*kj[2]; + + T rik = sqrt(rij*rkj); + T ct0 = (ij[0]*kj[0] + ij[1]*kj[1] + ij[2]*kj[2]) / rik; + if (ct0 < -pt999) + ct0 = -pt999; + else if (ct0 > pt999) + ct0 = pt999; + //T ct1 = std::max(-pt999, ct0); + //T ct2 = std::min( pt999, ct1); + + //T cst = ct2; + T ant = acos(ct0); + + T da = ant - teq; + // for rms deviation from ideal angles: + // eadev = eadev + da*da + T df = tk * da; + T eaw = df * da; + return eaw; +} +} +} +#endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index c68e2f9e36..b5309e59f2 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -253,7 +253,7 @@ DataSet_unsignedInt.o : DataSet_unsignedInt.cpp AssociatedData.h CpptrajFile.h D Deprecated.o : Deprecated.cpp CpptrajStdio.h Deprecated.h DispatchObject.h DihedralSearch.o : DihedralSearch.cpp ArgList.h Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h DihedralSearch.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h DistRoutines.o : DistRoutines.cpp Box.h DistRoutines.h ImageOption.h Matrix_3x3.h Parallel.h Vec3.h -Energy.o : Energy.cpp Atom.h AtomMask.h AtomType.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h Energy.h ExclusionArray.h FileName.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h +Energy.o : Energy.cpp Atom.h AtomMask.h AtomType.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h Energy.h Energy/Ene_Angle.h ExclusionArray.h FileName.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h EnergyArray.o : EnergyArray.cpp CpptrajFile.h CpptrajStdio.h EnergyArray.h FileIO.h FileName.h Parallel.h Energy_Sander.o : Energy_Sander.cpp ArgList.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h Energy_Sander.h FileName.h FileTypes.h File_TempName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ParmFile.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h EnsembleIn.o : EnsembleIn.cpp Atom.h AtomMask.h Box.h CoordinateInfo.h CpptrajStdio.h EnsembleIn.h FileName.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h Timer.h TrajFrameCounter.h Unit.h Vec3.h From a5de48d0ca61538829dbc577e5a19832ba3fd857 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 20 Jul 2022 08:32:50 -0400 Subject: [PATCH 2/7] Add nobondstoh keyword --- src/Action_Energy.cpp | 11 ++++++++--- src/Action_Energy.h | 1 + src/Energy.cpp | 6 ++++-- src/Energy.h | 2 +- 4 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/Action_Energy.cpp b/src/Action_Energy.cpp index 9c1ec1f8b6..6e5d1da93d 100644 --- a/src/Action_Energy.cpp +++ b/src/Action_Energy.cpp @@ -13,7 +13,8 @@ Action_Energy::Action_Energy() : EW_(0), dt_(0), need_lj_params_(false), - needs_exclList_(false) + needs_exclList_(false), + bondsToH_(true) {} /// DESTRUCTOR @@ -22,7 +23,7 @@ Action_Energy::~Action_Energy() { } void Action_Energy::Help() const { - mprintf("\t[] [] [out ]\n" + mprintf("\t[] [] [out ] [nobondstoh]\n" "\t[bond] [angle] [dihedral] {[nb14] | [e14] | [v14]}\n" "\t{[nonbond] | [elec] [vdw]} [kinetic [ketype {vel|vv}] [dt
]]\n" "\t[ etype { simple |\n" @@ -116,6 +117,8 @@ Action::RetType Action_Energy::Init(ArgList& actionArgs, ActionInit& init, int d if (KEtype_ != KE_VEL) dt_ = actionArgs.getKeyDouble("dt", 0.001); } + // Do we want bonds to H calculated? + bondsToH_ = !(actionArgs.hasKey("nobondstoh")); // Electrostatics type. std::string etypearg = actionArgs.GetStringKey("etype"); elecType_ = NO_ELE; @@ -235,6 +238,8 @@ Action::RetType Action_Energy::Init(ArgList& actionArgs, ActionInit& init, int d for (int i = 0; i != (int)TOTAL+1; i++) if (termEnabled[i]) mprintf(" '%s'", EtypeStr[i]); mprintf("\n"); + if (termEnabled[BOND] && !bondsToH_) + mprintf("\tNot calculating energy of bonds to hydrogen.\n"); if (elecType_ != NO_ELE) mprintf("\tElectrostatics method: %s\n", ElecStr[elecType_]); if (elecType_ == DIRECTSUM) { @@ -344,7 +349,7 @@ Action::RetType Action_Energy::DoAction(int frameNum, ActionFrame& frm) { switch (*calc) { case C_BND: time_bond_.Start(); - ene = ENE_.E_bond(frm.Frm(), *currentParm_, Mask1_); + ene = ENE_.E_bond(frm.Frm(), *currentParm_, Mask1_, bondsToH_); time_bond_.Stop(); Energy_[BOND]->Add(frameNum, &ene); Etot += ene; diff --git a/src/Action_Energy.h b/src/Action_Energy.h index fcdd2b0a0f..99eb861c0e 100644 --- a/src/Action_Energy.h +++ b/src/Action_Energy.h @@ -54,6 +54,7 @@ class Action_Energy: public Action { double dt_; ///< Time step for estimating kinetic energy (leapfrog) bool need_lj_params_; ///< True if LJ parameters needed. bool needs_exclList_; ///< True if Excluded_ needs to be set up. + bool bondsToH_; ///< True if we want to calculate energy of bonds with H Timer time_total_; Timer time_bond_; Timer time_angle_; diff --git a/src/Energy.cpp b/src/Energy.cpp index 2adb31b2b4..0c26a1a65a 100644 --- a/src/Energy.cpp +++ b/src/Energy.cpp @@ -15,11 +15,13 @@ const double Energy_Amber::QFAC = Constants::ELECTOAMBER * Constants::ELECTOAMBE Energy_Amber::Energy_Amber() : debug_(0) {} /** Bond energy */ -double Energy_Amber::E_bond(Frame const& fIn, Topology const& tIn, CharMask const& mask) +double Energy_Amber::E_bond(Frame const& fIn, Topology const& tIn, CharMask const& mask, bool bondsToH) { // Heavy atom bonds double Ebond = CalcBondEnergy(fIn, tIn.Bonds(), tIn.BondParm(), mask); - Ebond += CalcBondEnergy(fIn, tIn.BondsH(), tIn.BondParm(), mask); + // Bonds to hydrogen + if (bondsToH) + Ebond += CalcBondEnergy(fIn, tIn.BondsH(), tIn.BondParm(), mask); return Ebond; } diff --git a/src/Energy.h b/src/Energy.h index 3bb88a3250..08d5e4c674 100644 --- a/src/Energy.h +++ b/src/Energy.h @@ -12,7 +12,7 @@ class Energy_Amber { typedef std::vector Darray; Energy_Amber(); - double E_bond(Frame const&, Topology const&, CharMask const&); + double E_bond(Frame const&, Topology const&, CharMask const&, bool); double E_angle(Frame const&, Topology const&, CharMask const&); double E_torsion(Frame const&, Topology const&, CharMask const&); double E_14_Nonbond(Frame const&, Topology const&, CharMask const&, double&); From ccec9794915c9c47a8e3324e8ade587fc62dc30a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 20 Jul 2022 08:59:03 -0400 Subject: [PATCH 3/7] Change var name to theta --- src/Energy/Ene_Angle.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Energy/Ene_Angle.h b/src/Energy/Ene_Angle.h index c36453ade3..61987f6e8e 100644 --- a/src/Energy/Ene_Angle.h +++ b/src/Energy/Ene_Angle.h @@ -29,9 +29,9 @@ T Ene_Angle(T const* xyz1, T const* xyz2, T const* xyz3, T const& teq, T const& //T ct2 = std::min( pt999, ct1); //T cst = ct2; - T ant = acos(ct0); + T theta = acos(ct0); - T da = ant - teq; + T da = theta - teq; // for rms deviation from ideal angles: // eadev = eadev + da*da T df = tk * da; From 1764a0d5a64e214aa6c9d1f7e78a12b7b8e8f77d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 20 Jul 2022 09:06:46 -0400 Subject: [PATCH 4/7] Add nobondstoh option to energy --- doc/cpptraj.lyx | 6 +++++- src/Energy.cpp | 3 --- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index b99e66b481..394dc72e2e 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -24682,7 +24682,7 @@ energy \end_layout \begin_layout LyX-Code -energy [] [] [out ] +energy [] [] [out ] [nobondstoh] \end_layout \begin_layout LyX-Code @@ -24755,6 +24755,10 @@ energy [] [] [out ] ] File to write results to. \end_layout +\begin_layout Description +[nobondstoh] Skip calculating the energy of bonds to hydrogen. +\end_layout + \begin_layout Description [bond] Calculate bond energy. \end_layout diff --git a/src/Energy.cpp b/src/Energy.cpp index 0c26a1a65a..41d6a1bc6d 100644 --- a/src/Energy.cpp +++ b/src/Energy.cpp @@ -89,9 +89,6 @@ double Energy_Amber::CalcAngleEnergy(Frame const& fIn, AngleArray const& Angles, double ene = Cpptraj::Energy::Ene_Angle( fIn.XYZ(a->A1()), fIn.XYZ(a->A2()), fIn.XYZ(a->A3()), ap.Teq(), ap.Tk()); -// double theta = CalcAngle(fIn.XYZ(a->A1()), fIn.XYZ(a->A2()), fIn.XYZ(a->A3())); -// double tdiff = theta - ap.Teq(); -// double ene = ap.Tk() * (tdiff * tdiff); Eangle += ene; # ifdef DEBUG_ENERGY mprintf("\tAngle %4u %4i -- %4i -- %4i: k= %12.5f x0= %12.5f t= %12.5f E= %12.5e\n", From 6b7a9e6e68a68e04fdfd330042e6bbbeaa58f4d0 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 20 Jul 2022 09:08:53 -0400 Subject: [PATCH 5/7] Add code doc --- src/Energy/Ene_Angle.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Energy/Ene_Angle.h b/src/Energy/Ene_Angle.h index 61987f6e8e..2347cd7352 100644 --- a/src/Energy/Ene_Angle.h +++ b/src/Energy/Ene_Angle.h @@ -2,7 +2,7 @@ #define INC_ENERGY_ENE_ANGLE_H namespace Cpptraj { namespace Energy { - +/// \return Energy of the angle between xyz1, xyz2, and xyz3 (modeled after SANDER) template T Ene_Angle(T const* xyz1, T const* xyz2, T const* xyz3, T const& teq, T const& tk) { From b795e3ad1c70ee78c40ae81417b9ba6f2751c075 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 20 Jul 2022 09:11:27 -0400 Subject: [PATCH 6/7] 6.11.1. Revision bump for adding nobondstoh keyword to energy. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index 965a5cf0d2..216db18a32 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V6.11.0" +#define CPPTRAJ_INTERNAL_VERSION "V6.11.1" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif From 9b691303a2108604444eaa107e900349e5db82d8 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 20 Jul 2022 09:14:57 -0400 Subject: [PATCH 7/7] Should actually be 6.12.0 (minor bump) since change in energy angle calc can slightly change angle energy output. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index 216db18a32..fbdf1b5160 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V6.11.1" +#define CPPTRAJ_INTERNAL_VERSION "V6.12.0" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif