Add 'nobondstoh' keyword to energy. Minor change in angle energy calc.

doc/cpptraj.lyx

@@ -24682,7 +24682,7 @@ energy
 \end_layout
 
 \begin_layout LyX-Code
-energy [<name>] [<mask1>] [out <filename>]
+energy [<name>] [<mask1>] [out <filename>] [nobondstoh]
 \end_layout
 
 \begin_layout LyX-Code
@@ -24755,6 +24755,10 @@ energy [<name>] [<mask1>] [out <filename>]
 <filename>] 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

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[<name>] [<mask1>] [out <filename>]\n"
+  mprintf("\t[<name>] [<mask1>] [out <filename>] [nobondstoh]\n"
           "\t[bond] [angle] [dihedral] {[nb14] | [e14] | [v14]}\n"
           "\t{[nonbond] | [elec] [vdw]} [kinetic [ketype {vel|vv}] [dt <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;

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_;

src/Energy.cpp

@@ -7,18 +7,21 @@
 #include "Topology.h"
 #include "CharMask.h"
 #include "ExclusionArray.h"
+#include "Energy/Ene_Angle.h"
 
 const double Energy_Amber::QFAC = Constants::ELECTOAMBER * Constants::ELECTOAMBER;
 
 // CONSTRUCTOR
 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;
 }
 
@@ -43,8 +46,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 +86,9 @@ 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<double>(
+                     fIn.XYZ(a->A1()), fIn.XYZ(a->A2()), fIn.XYZ(a->A3()), ap.Teq(), ap.Tk());
+
       Eangle += ene;
 #     ifdef DEBUG_ENERGY
       mprintf("\tAngle %4u %4i -- %4i -- %4i: k= %12.5f  x0= %12.5f  t= %12.5f  E= %12.5e\n",

src/Energy.h

@@ -12,7 +12,7 @@ class Energy_Amber {
     typedef std::vector<double> 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&);

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 {
+/// \return Energy of the angle between xyz1, xyz2, and xyz3 (modeled after SANDER)
+template <typename T>
+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 theta = acos(ct0);
+
+  T da = theta - teq;
+  // for rms deviation from ideal angles:
+  // eadev = eadev + da*da
+  T df = tk * da;
+  T eaw = df * da;
+  return eaw;
+}
+}
+}
+#endif

src/Version.h

@@ -12,7 +12,7 @@
  * Whenever a number that precedes <revision> is incremented, all subsequent
  * numbers should be reset to 0.
  */
-#define CPPTRAJ_INTERNAL_VERSION "V6.11.0"
+#define CPPTRAJ_INTERNAL_VERSION "V6.12.0"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif

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