Recognize Drude particles, improve handling of lone pairs

configure

@@ -168,7 +168,7 @@ PNETCDF_OPTS='--disable-fortran --disable-cxx'
 
 FFTW_SRCTAR='fftw-3.3.9.tar.gz'
 FFTW_SRCDIR='fftw-3.3.9'
-FFTW_URL="ftp://ftp.fftw.org/pub/fftw/$FFTW_SRCTAR"
+FFTW_URL="https://fftw.org/pub/fftw/$FFTW_SRCTAR"
 FFTW_OPTS='--enable-threads --with-combined-threads --with-pic --disable-fortran'
 
 # ----- Variables for downloading MPI ------------

doc/cpptraj.lyx

@@ -14247,7 +14247,7 @@ The following topology related commands are available:
 \begin_layout Standard
 \align center
 \begin_inset Tabular
-<lyxtabular version="3" rows="22" columns="2">
+<lyxtabular version="3" rows="23" columns="2">
 <features tabularvalignment="middle">
 <column alignment="center" valignment="top">
 <column alignment="center" valignment="top">
@@ -14345,6 +14345,26 @@ Print bond info for selected atoms.
 <cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
 \begin_inset Text
 
+\begin_layout Plain Layout
+bondparminfo
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Print the bond parameter table.
+\end_layout
+
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
 \begin_layout Plain Layout
 change
 \end_layout
@@ -15057,6 +15077,67 @@ If 2 masks are given instead of 1,
  print info for bonds with first atom in <mask1> and second atom in <mask2>.
 \end_layout
 
+\begin_layout Subsection
+bondparminfo
+\end_layout
+
+\begin_layout LyX-Code
+bondparminfo [parm <name> | crdset <set> | parmindex <#> | <#>] [out <file>]
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+[parm
+\begin_inset space ~
+\end_inset
+
+<name>
+\begin_inset space ~
+\end_inset
+
+|
+\begin_inset space ~
+\end_inset
+
+parmindex
+\begin_inset space ~
+\end_inset
+
+<#>
+\begin_inset space ~
+\end_inset
+
+|
+\begin_inset space ~
+\end_inset
+
+<#>] Name/tag or index of topology.
+ Default is first loaded topology.
+\end_layout
+
+\begin_layout Description
+[out
+\begin_inset space ~
+\end_inset
+
+<file>] File to print to (default STDOUT).
+\end_layout
+
+\end_deeper
+\begin_layout Standard
+Print the bond parameter table with format:
+\end_layout
+
+\begin_layout LyX-Code
+#Idx Rk Req
+\end_layout
+
+\begin_layout Standard
+Where Idx is the internal bond parameter index,
+ Rk is the bond force constant (in kcal/mol*Ang^2),
+ and Req is the bond equilibrium value (in Ang.).
+\end_layout
+
 \begin_layout Subsection
 change
 \end_layout
@@ -43023,7 +43104,7 @@ vector [<name>] <Type> [out <filename> [ptrajoutput]] [<mask1>] [<mask2>]
 \end_layout
 
 \begin_layout LyX-Code
-   	[magnitude] [ired] [gridset <grid>]
+   	[magnitude] [geom] [ired] [gridset <grid>] [debye]
 \end_layout
 
 \begin_layout LyX-Code
@@ -43103,6 +43184,11 @@ ptraj
 [magnitude] Store the magnitude of the vector with aspect [Mag].
 \end_layout
 
+\begin_layout Description
+[geom] If specified,
+ use geometric centers instead of centers of mass.
+\end_layout
+
 \begin_layout Description
 [ired] Mark this vector for subsequent IRED analysis with commands 'matrix ired' and 'ired'.
 \end_layout
@@ -43129,6 +43215,14 @@ ucell[x|y|z]
 .
 \end_layout
 
+\begin_layout Description
+[debye] (
+\series bold
+'dipole'
+\series default
+ vector only) Report dipole vector in units of Debye instead of e-*Ang.
+\end_layout
+
 \begin_layout Standard
 Data Sets Created:
 \end_layout
@@ -43186,8 +43280,13 @@ dipole Store the dipole and center of mass of the atoms specified in
 <mask1>
 \series default
 .
- The vector is not converted to appropriate units,
- nor is the value well-defined if the atoms in the mask are not overall charge neutral.
+ The dipole vector has units of e-*Ang unless 
+\series bold
+'debye'
+\series default
+ is specified for units of Debye.
+ The center is always stored as simply Ang (since it is just coordinates).
+ Note that the value may not be well-defined if the atoms in the mask are not overall charge neutral.
 \end_layout
 
 \begin_layout Description

src/Action_Angle.cpp

@@ -62,6 +62,12 @@ Action::RetType Action_Angle::Setup(ActionSetup& setup) {
     mprintf("Warning: angle: One or more masks contain 0 atoms.\n");
     return Action::SKIP;
   }
+  // Check if center of mass is possible
+  if (useMass_) {
+    if (setup.Top().MaskHasZeroMass(Mask1_)) useMass_ = false;
+    if (setup.Top().MaskHasZeroMass(Mask2_)) useMass_ = false;
+    if (setup.Top().MaskHasZeroMass(Mask3_)) useMass_ = false;
+  }
   return Action::OK;  
 }
 

src/Action_Dihedral.cpp

@@ -91,7 +91,13 @@ Action::RetType Action_Dihedral::Setup(ActionSetup& setup) {
     mprintf("Warning: One or more masks have no atoms.\n");
     return Action::SKIP;
   }
-
+  // Check if center of mass is possible
+  if (useMass_) {
+    if (setup.Top().MaskHasZeroMass(M1_)) useMass_ = false;
+    if (setup.Top().MaskHasZeroMass(M2_)) useMass_ = false;
+    if (setup.Top().MaskHasZeroMass(M3_)) useMass_ = false;
+    if (setup.Top().MaskHasZeroMass(M4_)) useMass_ = false;
+  }
   return Action::OK;  
 }
 

src/Action_Distance.cpp

@@ -136,6 +136,13 @@ Action::RetType Action_Distance::Setup(ActionSetup& setup) {
   else
     mprintf(", imaging off");
   mprintf(".\n");
+  // Check if center of mass not possible
+  if (useMass_) {
+    if (setup.Top().MaskHasZeroMass( Mask1_ )) useMass_ = false;
+    if (mode_ == NORMAL) {
+      if (setup.Top().MaskHasZeroMass( Mask2_ )) useMass_ = false;
+    }
+  }
 
   return Action::OK;  
 }

src/Action_Vector.cpp

@@ -15,14 +15,16 @@ Action_Vector::Action_Vector() :
   vcorr_(0),
   ptrajoutput_(false),
   needBoxInfo_(false),
+  useMass_(true),
+  dipole_in_debye_(false),
   CurrentParm_(0),
   outfile_(0)
 {}
 
 // Action_Vector::Help()
 void Action_Vector::Help() const {
   mprintf("\t[<name>] <Type> [out <filename> [ptrajoutput]] [<mask1>] [<mask2>]\n"
-          "\t[magnitude] [ired] [gridset <grid>]\n"
+          "\t[magnitude] [geom] [ired] [gridset <grid>] [debye]\n"
           "\t<Type> = { mask     | minimage  | dipole | center   | corrplane | \n"
           "\t           box      | boxcenter | ucellx | ucelly   | ucellz    | \n"
           "\t           momentum | principal [x|y|z]  | velocity | force       }\n" 
@@ -38,7 +40,10 @@ void Action_Vector::Help() const {
           "    momentum         : Store total momentum vector of atoms in <mask1> (requires velocities).\n"
           "    principal [x|y|z]: X, Y, or Z principal axis vector for atoms in <mask1>.\n"
           "    velocity         : Store velocity of atoms in <mask1> (requires velocities).\n"
-          "    force            : Store force of atoms in <mask1> (requires forces).\n");
+          "    force            : Store force of atoms in <mask1> (requires forces).\n"
+          "  If 'magnitude' is specified, also calculate the vector magnitude.\n"
+          "  If 'geom' is specified, use geometric centers, otherwise use center of mass.\n"
+          "  If 'debye' is specified, report 'dipole' vector in Debye instead of e-*Ang.\n");
 }
 
 // DESTRUCTOR
@@ -97,6 +102,7 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d
     mprinterr("Error: 'ptrajoutput' and 'magnitude' are incompatible.\n");
     return Action::ERR;
   }
+  dipole_in_debye_ = actionArgs.hasKey("debye");
   needBoxInfo_ = false;
   // Deprecated: corrired, corr, ired
   if ( actionArgs.hasKey("principal") ) {
@@ -209,6 +215,12 @@ 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 (dipole_in_debye_)
+      mprintf("\tDipole vector units will be Debye.\n");
+    else
+      mprintf("\tDipole vector units will be e-*Ang.\n");
+  }
 
   return Action::OK;
 }
@@ -257,6 +269,11 @@ Action::RetType Action_Vector::Setup(ActionSetup& setup) {
       return Action::ERR;
     }
   }
+  // Check if center of mass is possible
+  if (useMass_) {
+    if (mask_.MaskStringSet() && setup.Top().MaskHasZeroMass( mask_ )) useMass_ = false;
+    if (mask2_.MaskStringSet() && setup.Top().MaskHasZeroMass( mask2_ )) useMass_ = false;
+  }
   CurrentParm_ = setup.TopAddress();
   return Action::OK;
 }
@@ -371,10 +388,19 @@ Vec3 Action_Vector::leastSquaresPlane(int n, const double* vcorr) {
 }
 
 // -----------------------------------------------------------------------------
+Vec3 Action_Vector::GetVec(Frame const& currentFrame, AtomMask const& maskIn)
+const
+{
+  if (useMass_)
+    return currentFrame.VCenterOfMass(maskIn);
+  else
+    return currentFrame.VGeometricCenter(maskIn);
+}
+
 // Action_Vector::Mask()
 void Action_Vector::Mask(Frame const& currentFrame) {
-  Vec3 CXYZ = currentFrame.VCenterOfMass(mask_);
-  Vec3 VXYZ = currentFrame.VCenterOfMass(mask2_);
+  Vec3 CXYZ = GetVec(currentFrame, mask_);
+  Vec3 VXYZ = GetVec(currentFrame, mask2_);
   VXYZ -= CXYZ;
   Vec_->AddVxyzo(VXYZ, CXYZ);
 }
@@ -396,6 +422,8 @@ void Action_Vector::Dipole(Frame const& currentFrame) {
     VXYZ += ( XYZ );
   }
   CXYZ /= total_mass;
+  if (dipole_in_debye_)
+    VXYZ /= Constants::DEBYE_EA;
   Vec_->AddVxyzo( VXYZ, CXYZ );
 }
 
@@ -421,7 +449,7 @@ void Action_Vector::Principal(Frame const& currentFrame) {
 
 // Action_Vector::CorrPlane()
 void Action_Vector::CorrPlane(Frame const& currentFrame) {
-  Vec3 CXYZ = currentFrame.VCenterOfMass(mask_);
+  Vec3 CXYZ = GetVec(currentFrame, mask_);
   int idx = 0;
   for (AtomMask::const_iterator atom = mask_.begin();
                               atom != mask_.end(); ++atom)
@@ -454,8 +482,8 @@ void Action_Vector::BoxLengths(Box const& box) {
 
 // Action_Vector::MinImage()
 void Action_Vector::MinImage(Frame const& frm) {
-  Vec3 com1 = frm.VCenterOfMass(mask_);
-  Vec_->AddVxyzo( MinImagedVec(com1, frm.VCenterOfMass(mask2_), frm.BoxCrd().UnitCell(), frm.BoxCrd().FracCell()), com1 );
+  Vec3 com1 = GetVec(frm, mask_);
+  Vec_->AddVxyzo( MinImagedVec(com1, GetVec(frm, mask2_), frm.BoxCrd().UnitCell(), frm.BoxCrd().FracCell()), com1 );
 }
 
 /// \return The center of selected elements in given array.
@@ -476,7 +504,7 @@ static inline Vec3 CalcCenter(const double* xyz, AtomMask const& maskIn) {
 Action::RetType Action_Vector::DoAction(int frameNum, ActionFrame& frm) {
   switch ( mode_ ) {
     case MASK        : Mask(frm.Frm()); break;
-    case CENTER      : Vec_->AddVxyz( frm.Frm().VCenterOfMass(mask_) ); break;
+    case CENTER      : Vec_->AddVxyz( GetVec(frm.Frm(), mask_) ); break;
     case MOMENTUM    : Vec_->AddVxyz( frm.Frm().VMomentum(mask_) ); break;
     case VELOCITY    : Vec_->AddVxyz( CalcCenter(frm.Frm().vAddress(), mask_) ); break;
     case FORCE       : Vec_->AddVxyz( CalcCenter(frm.Frm().fAddress(), mask_) ); break; 

src/Action_Vector.h

@@ -26,6 +26,8 @@ class Action_Vector : public Action {
 
     static double solve_cubic_eq(double,double,double,double);
     static Vec3 leastSquaresPlane(int,const double*);
+    /// \return Center of mass or geometric center of atoms in given mask
+    inline Vec3 GetVec(Frame const&, AtomMask const&) const;
     void Mask(Frame const&);
     void Dipole(Frame const&);
     void Principal(Frame const&);
@@ -40,7 +42,9 @@ class Action_Vector : public Action {
     double* vcorr_;         ///< Temp. space for calculating CorrPlane
     vectorMode mode_;       ///< Vector calculation mode
     bool ptrajoutput_;      ///< If true output in ptraj format
-    bool needBoxInfo_;      ///< If true box info required. 
+    bool needBoxInfo_;      ///< If true box info required.
+    bool useMass_;          ///< If true, centers are mass-weighted 
+    bool dipole_in_debye_;  ///< If true, report dipole vector values in Debye
     Topology* CurrentParm_; ///< Current topology (for dipole)
     AtomMask mask_;
     AtomMask mask2_;