Add 'bonddipole' mode to 'vector' action (#1154)

* Start vector bond dipole calc

* Add missing printf var

* Average origin

* Enable bonddipole calc. Not yet working right

* Add Pauling electronegativity values for atoms

* Try to calculate bond dipoles using more chemical intuition

* Warn about unhandled numbers of bonds

* Use summation for bonds with multiple single substituents

* Generalize

* Try a hybrid method where the dipole vector is -SUM(R x Q) and origin is calculated from individual bonds, using electronegativity to orient them

* Shift the origin so that the center of the dipole vector will end up at the calculated bond average origin

* Hide some debug info

* Add a test

* Add warning

* 6.29.18. Revision bump for vector bonddipole

---------

Co-authored-by: Daniel R. Roe <daniel.roe@nih.gov>

Author: drroe

src/Action_Vector.cpp

@@ -2,6 +2,7 @@
 #include <algorithm> // sort
 #include "Constants.h" // PI
 #include "Action_Vector.h"
+#include "CharMask.h"
 #include "CpptrajStdio.h"
 #include "DistRoutines.h" // MinImagedVec, includes Matrix_3x3 for principal
 #include "DataSet_Vector.h"
@@ -18,9 +19,17 @@ Action_Vector::Action_Vector() :
   useMass_(true),
   dipole_in_debye_(false),
   CurrentParm_(0),
-  outfile_(0)
+  outfile_(0),
+  cmask_(0),
+  debug_(0)
 {}
 
+// DESTRUCTOR
+Action_Vector::~Action_Vector() {
+  if (vcorr_!=0) delete[] vcorr_;
+  if (cmask_ != 0) delete cmask_;
+}
+
 // Action_Vector::Help()
 void Action_Vector::Help() const {
   mprintf("\t[<name>] <Type> [out <filename> [ptrajoutput]] [<mask1>] [<mask2>]\n"
@@ -46,23 +55,20 @@ void Action_Vector::Help() const {
           "  If 'debye' is specified, report 'dipole' vector in Debye instead of e-*Ang.\n");
 }
 
-// DESTRUCTOR
-Action_Vector::~Action_Vector() {
-  if (vcorr_!=0) delete[] vcorr_;
-}
-
 const char* Action_Vector::ModeString_[] = {
   "NO_OP", "Principal X", "Principal Y", "Principal Z",
   "Dipole", "Box", "Mask",
   "CorrPlane", "Center", "Unit cell X", "Unit cell Y", "Unit cell Z",
-  "Box Center", "MinImage", "Momentum", "Velocity", "Force"
+  "Box Center", "MinImage", "Momentum", "Velocity", "Force",
+  "Bond Dipole"
 };
 
 const bool Action_Vector::NeedsOrigin_[] = {
   false, true, true, true,
   true, false, true,
   true, false, true, true, true,
-  false, true, false, false, false
+  false, true, false, false, false,
+  true
 };
 
 static Action::RetType WarnDeprecated() {
@@ -82,6 +88,7 @@ static inline Action::RetType DeprecatedErr(const char* key) {
 // Action_Vector::Init()
 Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
 {
+  debug_ = debugIn;
   DataFile* df = 0;
   std::string filename = actionArgs.GetStringKey("out");
   if (actionArgs.hasKey("trajout")) return DeprecatedErr( "trajout" );
@@ -120,6 +127,8 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d
     mode_ = FORCE;
   else if (actionArgs.hasKey("dipole"))
     mode_ = DIPOLE;
+  else if (actionArgs.hasKey("bonddipole"))
+    mode_ = BONDDIPOLE;
   else if (actionArgs.hasKey("box"))
     mode_ = BOX;
   else if (actionArgs.hasKey("corrplane"))
@@ -176,6 +185,9 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d
     }
     if (mask2_.SetMaskString( maskexpr )) return Action::ERR;
   }
+  // Allocate Char mask for bond dipole
+  if (mode_ == BONDDIPOLE)
+    cmask_ = new CharMask();
   // Set up vector dataset and IRED status
   MetaData md(actionArgs.GetStringNext(), MetaData::M_VECTOR);
   if (isIred) md.SetScalarType( MetaData::IREDVEC );
@@ -215,7 +227,9 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d
   mprintf("\n");
   if (gridSet_ != 0)
     mprintf("\tExtracting box vectors from grid set '%s'\n", gridSet_->legend());
-  if (mode_ == DIPOLE) {
+  if (mode_ == DIPOLE || mode_ == BONDDIPOLE) {
+    if (mode_ == BONDDIPOLE) //TODO remove when no longer testing
+      mprintf("Warning: 'bonddipole' mode is still experimental. Use with caution.\n");
     if (dipole_in_debye_)
       mprintf("\tDipole vector units will be Debye.\n");
     else
@@ -252,6 +266,10 @@ Action::RetType Action_Vector::Setup(ActionSetup& setup) {
       mprinterr("Error: First vector mask is empty.\n");
       return Action::ERR;
     }
+    // Set up char mask if needed
+    if (cmask_ != 0) {
+      *cmask_ = CharMask( mask_.ConvertToCharMask(), mask_.Nselected() );
+    }
   }
 
   // Allocate space for CORRPLANE. 
@@ -405,6 +423,187 @@ void Action_Vector::Mask(Frame const& currentFrame) {
   Vec_->AddVxyzo(VXYZ, CXYZ);
 }
 
+/** Calculate dipole */
+bool Action_Vector::calcBondDipole(Vec3& vBond, Vec3& vDipole, Vec3& vOrigin, double q0, Vec3 const& XYZ0, Vec3 const& XYZ1)
+{
+  bool hasCharge = true;
+  if (q0 < 0.0) {
+    vBond = XYZ0 - XYZ1;
+    vDipole = vBond * (-q0);
+    vOrigin = XYZ1;
+  } else if (q0 > 0.0) {
+    vBond = XYZ1 - XYZ0;
+    vDipole = vBond * q0;
+    vOrigin = XYZ0;
+  } else
+    hasCharge = false;
+  return hasCharge;
+}
+
+/// \return true if all atoms bonded to atomIn (ignoring ignoreIdx) have only the one bond back to atomIn
+static inline bool bondedAtomsHave1bond(Topology const& topIn, Atom const& atomIn, int ignoreIdx)
+{
+  for (int ii=0; ii < atomIn.Nbonds(); ii++) {
+    int bondedAtomIdx = atomIn.Bond(ii);
+    if (bondedAtomIdx != ignoreIdx) {
+      if (topIn[bondedAtomIdx].Nbonds() != 1) return false;
+    }
+  }
+  return true;
+}
+
+static inline double bondedAtomsCharge(Topology const& topIn, Atom const& atomIn, int ignoreIdx)
+{
+  double sumQ = 0.0;
+  for (int ii=0; ii < atomIn.Nbonds(); ii++) {
+    int bondedAtomIdx = atomIn.Bond(ii);
+    if (bondedAtomIdx != ignoreIdx) {
+      if (topIn[bondedAtomIdx].Nbonds() == 1)
+        sumQ += topIn[bondedAtomIdx].Charge();
+    }
+  }
+  return sumQ;
+}
+
+/** Calculate net dipole from individual bond dipoles. */
+void Action_Vector::BondDipole_individualBonds(Frame const& currentFrame) {
+  Vec3 VXYZ(0.0, 0.0, 0.0); // Negative, head
+  Vec3 OXYZ(0.0, 0.0, 0.0); // Positive, tail
+  Topology const& currentTop = *CurrentParm_;
+
+  unsigned int Nbonds = 0;
+  for (AtomMask::const_iterator iat = mask_.begin(); iat != mask_.end(); ++iat)
+  {
+    Atom const& atom0 = currentTop[*iat];
+    Vec3 XYZ0( currentFrame.XYZ(*iat) );
+    double q0 = atom0.Charge();
+    for (Atom::bond_iterator bit = atom0.bondbegin(); bit != atom0.bondend(); ++bit)
+    {
+      if (*bit > *iat && cmask_->AtomInCharMask( *bit )) {
+        // Both atoms in the bond are selected.
+        Atom const& atom1 = currentTop[*bit];
+        Vec3 XYZ1( currentFrame.XYZ(*bit) );
+        double q1 = atom1.Charge();
+        Vec3 vDipole(0.0, 0.0, 0.0);
+        Vec3 vBond(0.0, 0.0, 0.0);
+        Vec3 vOrigin(0.0, 0.0, 0.0);
+        bool hasCharge = false;
+        if (atom0.Element() != atom1.Element()) {
+          if (atom0.Element()==Atom::HYDROGEN || atom0.Nbonds() == 1) {
+            hasCharge = calcBondDipole( vBond, vDipole, vOrigin, q0, XYZ0, XYZ1 );
+          } else if (atom1.Element()==Atom::HYDROGEN || atom1.Nbonds() == 1) {
+            hasCharge = calcBondDipole( vBond, vDipole, vOrigin, q1, XYZ1, XYZ0 );
+          } else if (atom0.Nbonds() > 1 && bondedAtomsHave1bond(currentTop, atom0, *bit)) {
+            // atom0 is bonded to a single other atom, X-atom0-atom1. Sum up X and atom0 charge
+            double q0X = q0 + bondedAtomsCharge(currentTop, atom0, *bit);
+            if (debug_ > 1)
+              mprintf("DEBUG: Total charge on %s is %f\n", currentTop.AtomMaskName(*iat).c_str(), q0X);
+            hasCharge = calcBondDipole( vBond, vDipole, vOrigin, q0X, XYZ0, XYZ1 );
+          } else if (atom1.Nbonds() > 1 && bondedAtomsHave1bond(currentTop, atom1, *iat)) {
+            // atom1 is bonded to a single other atom, X-atom1-atom0. Sum up X and atom1 charge
+            double q1X = q1 + bondedAtomsCharge(currentTop, atom1, *iat);
+            if (debug_ > 1)
+            mprintf("DEBUG: Total charge on %s is %f\n", currentTop.AtomMaskName(*bit).c_str(), q1X);
+            hasCharge = calcBondDipole( vBond, vDipole, vOrigin, q1X, XYZ1, XYZ0 );
+          } else {
+            mprintf("Warning: Unhandled number of bonds: %s=%i, %s=%i\n",
+                     currentTop.AtomMaskName(*iat).c_str(), atom0.Nbonds(),
+                     currentTop.AtomMaskName(*bit).c_str(), atom1.Nbonds());
+          }
+        }
+
+        if (hasCharge) {
+          // Subtract half the dipole from the bond center to get the dipole origin
+          //Vec3 vBondCtr = vBond / 2.0;
+          //Vec3 vDipoleHalf = vDipole / 2.0;
+          //Vec3 vOrigin = vOrigin + vBondCtr - vDipoleHalf;
+          // DEBUG
+          if (debug_ > 1) {
+            mprintf("DEBUG: Bond %s(%f) -- %s(%f) bxyz={%f %f %f} oxyz={%f %f %f} vxyz={%f %f %f} mag=%f\n",
+                    currentTop.AtomMaskName(*iat).c_str(), q0,
+                    currentTop.AtomMaskName(*bit).c_str(), q1,
+                    vBond[0], vBond[1], vBond[2],
+                    vOrigin[0], vOrigin[1], vOrigin[2],
+                    vDipole[0], vDipole[1], vDipole[2],
+                    sqrt(vDipole.Magnitude2()));
+          }
+          // Sum
+          VXYZ += vDipole;
+          OXYZ += vOrigin;
+          Nbonds++;
+        }
+      }
+    } // END loop over bonded atoms
+  } // END loop over mask atoms
+  if (dipole_in_debye_)
+    VXYZ /= Constants::DEBYE_EA;
+  OXYZ /= (double)Nbonds;
+  Vec_->AddVxyzo( VXYZ, OXYZ );
+}
+
+/** Calculate net dipole from atoms, origin from individual bond dipoles. */
+void Action_Vector::BondDipole_net_bondOrigin(Frame const& currentFrame) {
+  Vec3 VXYZ(0.0, 0.0, 0.0); // Negative, head
+  Vec3 OXYZ(0.0, 0.0, 0.0); // Positive, tail
+  Topology const& currentTop = *CurrentParm_;
+
+  unsigned int Nbonds = 0;
+  for (AtomMask::const_iterator iat = mask_.begin(); iat != mask_.end(); ++iat)
+  {
+    Atom const& atom0 = currentTop[*iat];
+    Vec3 XYZ0( currentFrame.XYZ(*iat) );
+    double q0 = atom0.Charge();
+    double en0 = atom0.PaulingElectroNeg();
+    // Note that negative charge will decrease the vector (towards the origin)
+    // but by convention we want to point towards the negative charge, so at
+    // the end the vector will have to be negated.
+    VXYZ += ( XYZ0 * q0 );
+    // Store bond origins according to electronegativity.
+    for (Atom::bond_iterator bit = atom0.bondbegin(); bit != atom0.bondend(); ++bit)
+    {
+      if (*bit > *iat && cmask_->AtomInCharMask( *bit )) {
+        // Both atoms in the bond are selected.
+        Atom const& atom1 = currentTop[*bit];
+        Vec3 XYZ1( currentFrame.XYZ(*bit) );
+        double en1 = atom1.PaulingElectroNeg();
+        Vec3 vOrigin(0.0, 0.0, 0.0);
+        bool en_is_different = true;
+        if (en0 > en1) {
+          // atom0 is more negative
+          //vBond = XYZ0 - XYZ1;
+          vOrigin = XYZ1;
+        } else if (en1 > en0) {
+          // atom1 is more negative
+          //vBond = XYZ1 - XYZ0;
+          vOrigin = XYZ0;
+        } else
+          en_is_different = false;
+
+        if (en_is_different) {
+          // DEBUG
+          if (debug_ > 1) {
+            double q1 = atom1.Charge(); // DEBUG
+            mprintf("DEBUG: Bond %s(%f) -- %s(%f) oxyz={%f %f %f}\n",
+                    currentTop.AtomMaskName(*iat).c_str(), q0,
+                    currentTop.AtomMaskName(*bit).c_str(), q1,
+                    vOrigin[0], vOrigin[1], vOrigin[2]);
+          }
+          // Sum
+          OXYZ += vOrigin;
+          Nbonds++;
+        }
+      }
+    } // END loop over bonded atoms
+  } // END loop over mask atoms
+  if (dipole_in_debye_)
+    VXYZ /= Constants::DEBYE_EA;
+  OXYZ /= (double)Nbonds;
+  VXYZ.Neg();
+  // Shift it so the center of VXYZ will end up at the origin.
+  OXYZ -= (VXYZ / 2.0);
+  Vec_->AddVxyzo( VXYZ, OXYZ );
+}
+
 // Action_Vector::Dipole()
 void Action_Vector::Dipole(Frame const& currentFrame) {
   Vec3 VXYZ(0.0, 0.0, 0.0);
@@ -509,6 +708,8 @@ Action::RetType Action_Vector::DoAction(int frameNum, ActionFrame& frm) {
     case VELOCITY    : Vec_->AddVxyz( CalcCenter(frm.Frm().vAddress(), mask_) ); break;
     case FORCE       : Vec_->AddVxyz( CalcCenter(frm.Frm().fAddress(), mask_) ); break; 
     case DIPOLE      : Dipole(frm.Frm()); break;
+//    case BONDDIPOLE  : BondDipole_individualBonds(frm.Frm()); break;
+    case BONDDIPOLE  : BondDipole_net_bondOrigin(frm.Frm()); break;
     case PRINCIPAL_X :
     case PRINCIPAL_Y :
     case PRINCIPAL_Z : Principal(frm.Frm()); break;

src/Action_Vector.h

@@ -3,6 +3,7 @@
 #include "Action.h"
 class DataSet_Vector;
 class DataSet_3D;
+class CharMask;
 class Action_Vector : public Action {
   public:
     Action_Vector();
@@ -14,7 +15,8 @@ class Action_Vector : public Action {
       NO_OP=0,   PRINCIPAL_X, PRINCIPAL_Y, PRINCIPAL_Z,
       DIPOLE,    BOX,         MASK,
       CORRPLANE, CENTER,      BOX_X,       BOX_Y,       BOX_Z,
-      BOX_CTR,   MINIMAGE,    MOMENTUM,    VELOCITY,    FORCE
+      BOX_CTR,   MINIMAGE,    MOMENTUM,    VELOCITY,    FORCE,
+      BONDDIPOLE
     };
     static const char* ModeString_[];
     static const bool NeedsOrigin_[];
@@ -29,6 +31,9 @@ class Action_Vector : public Action {
     /// \return Center of mass or geometric center of atoms in given mask
     inline Vec3 GetVec(Frame const&, AtomMask const&) const;
     void Mask(Frame const&);
+    static inline bool calcBondDipole(Vec3&, Vec3&, Vec3&, double, Vec3 const&, Vec3 const&);
+    void BondDipole_individualBonds(Frame const&);
+    void BondDipole_net_bondOrigin(Frame const&);
     void Dipole(Frame const&);
     void Principal(Frame const&);
     void CorrPlane(Frame const&);
@@ -49,5 +54,7 @@ class Action_Vector : public Action {
     AtomMask mask_;
     AtomMask mask2_;
     CpptrajFile* outfile_;
+    CharMask* cmask_;
+    int debug_;
 };
 #endif

src/Atom.cpp

@@ -75,6 +75,32 @@ const double Atom::AtomicElementRadius_[NUMELEMENTS_] = { 1.0,
   0.000, 0.0 /* extra point/Drude has no radius */
 };
 
+/** Pauling electronegativity values from:
+  * A.L. Allred,
+  * Electronegativity values from thermochemical data,
+  * Journal of Inorganic and Nuclear Chemistry,
+  * Volume 17, Issues 3–4, 1961, Pages 215-221,
+  * https://doi.org/10.1016/0022-1902(61)80142-5.
+  * https://www.sciencedirect.com/science/article/pii/0022190261801425
+  * Atoms with no electronegativity value are assigned -1.0.
+  */
+const double Atom::PaulingElectroNeg_[NUMELEMENTS_] = { -1.0,
+  2.20, 2.04, 2.55, 3.04, 3.44, 3.98,
+  2.19, 2.58, 3.16, 2.96, 1.83, 1.00,
+  2.66, 1.31, 1.90, 0.98, 0.82, 0.82,
+  0.79, 1.65, 0.93, 1.61, -1.0, 2.18,
+  1.93, 2.54, 2.20, 1.57, 0.89, 2.02,
+  1.66, 1.88, 1.69, 0.70, 1.81, 2.01,
+  -1.0, 1.30, 2.00, 1.78, 2.20, 3.00,
+  1.55, 2.16, -1.0, 1.91, 1.60, 2.20,
+  2.20, 2.28, 2.33, 2.00, 2.20, 2.28,
+  1.90, -1.0, 0.90, 1.90, 1.36, 2.55,
+  0.95, 1.96, 2.05, 1.54, 1.90, 2.10,
+  1.50, 1.62, 1.63, 2.36, 2.60, 1.33,
+  1.22, 1.27, 1.17,
+  -1.0, -1.0 /* extra point/Drude has no EN value. TODO should that change? */
+};
+
 // CONSTRUCTOR
 Atom::Atom() : charge_(0.0), polar_(0.0), mass_(1.0), gb_radius_(0.0),
                gb_screen_(0.0), aname_(""), atype_(""), atype_index_(0),

src/Atom.h

@@ -81,6 +81,7 @@ class Atom {
     inline int AtomicNumber()          const { return AtomicElementNum_[element_];  }
     inline const char* ElementName()   const { return AtomicElementName_[element_]; }
     inline double ElementRadius()      const { return AtomicElementRadius_[element_]; }
+    double PaulingElectroNeg()         const { return PaulingElectroNeg_[element_]; }
     inline const NameType& Name()      const { return aname_; }
     inline const NameType& Type()      const { return atype_; }
     inline int TypeIndex()             const { return atype_index_; }
@@ -108,6 +109,7 @@ class Atom {
     static CPPTRAJ_EXPORT const char* AtomicElementName_[];
     static CPPTRAJ_EXPORT const double AtomicElementMass_[];
     static CPPTRAJ_EXPORT const double AtomicElementRadius_[];
+    static CPPTRAJ_EXPORT const double PaulingElectroNeg_[];
     double charge_;    ///< Charge in e-
     double polar_;     ///< Atomic polarizability in Ang^3
     double mass_;      ///< mass in amu

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.29.17"
+#define CPPTRAJ_INTERNAL_VERSION "V6.29.18"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif