Add energy decomposition (with PME)

configure

@@ -41,7 +41,7 @@ UsageFull() {
   echo "    -debugon   : Add -DDEBUG flag to activate additional internal debugging."
   echo "    -nolfs     : Do not enable large file support."
   echo "    -shared    : Configure for generating libcpptraj (implies -nosanderlib)."
-  echo "    -fftw3     : Use FFTW instead of pubfft for FFT."
+  echo "    -pubfft    : Use pubfft instead of FFTW for FFT (disables PME functionality)."
   echo "    -windows   : Set up for use with MinGW compilers for a native Windows build."
   echo "  LIBRARY OPTIONS"
   echo "    -<lib>             : Enable library."
@@ -1038,7 +1038,7 @@ EOF
       if [ "${LIB_STAT[$LFFTW3]}" = 'optional' ] ; then
         WrnMsg "FFTW test failed. CPPTRAJ will use the built-in PubFFT routines."
         if [ $BUILD_LIBS -eq 1 ] ; then
-          echo "To force CPPTRAJ to build its own FFTW, reconfigure and specify '-fftw3'."
+          echo "To force CPPTRAJ to build its own FFTW, reconfigure and specify '--buildlibs'."
         fi
         LIB_STAT[$LFFTW3]='off'
       elif [ "${LIB_STAT[$LFFTW3]}" = 'enabled' ] ; then
@@ -2605,6 +2605,7 @@ while [ ! -z "$1" ] ; do
     '-tune'          ) USE_OPT=2 ;;
     '-noc++11'       ) C11_SUPPORT='no' ;;
     '-windows'       ) PLATFORM='windows' ;;
+    '-pubfft'        ) LIB_STAT[$LFFTW3]='off' ;;
     # Cpptraj options
     '-nolfs'           ) LFS='' ;;
     '-single-ensemble' ) USE_SINGLEENSEMBLE=1 ;;

devtools/CreateMake.sh

@@ -38,6 +38,11 @@ DEL_FILE = /bin/rm -f
 # Objects
 OBJECTS=\$($VARNAME:.cpp=.o)
 
+# General rules
+.cpp.o:
+	\$(VB)echo CXX $<
+	\$(VB)\$(CXX) \$(DIRECTIVES) \$(CPPTRAJ_INC) \$(CXXFLAGS) -c -o \$@ $<
+
 # Default target: objects
 all: \$(OBJECTS)
 

doc/DocumentChecksums.txt

@@ -1,3 +1,3 @@
 b37726e7a841f6fc695ecd7fb040ffbf  CpptrajDevelopmentGuide.lyx
-736ee725ca9cc67d349f67d714f41c20  cpptraj.lyx
+151841b1d3283890c425aaae5a923e7f  cpptraj.lyx
 07c4039e732fc2eb1df0c1e0863cb949  CpptrajManual.lyx

doc/cpptraj.lyx

@@ -20478,7 +20478,7 @@ Mod
 \begin_layout Standard
 \align center
 \begin_inset Tabular
-<lyxtabular version="3" rows="92" columns="3">
+<lyxtabular version="3" rows="93" columns="3">
 <features islongtable="true" longtabularalignment="center">
 <column alignment="center" valignment="top" width="25text%">
 <column alignment="center" valignment="top" width="65text%">
@@ -21385,6 +21385,35 @@ Calculate secondary structure content using the DSSP algorithm
 <cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
 \begin_inset Text
 
+\begin_layout Plain Layout
+enedecomp
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Perform per-atom energy decomposition.
+\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
+
+\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
 energy
 \end_layout
@@ -27788,6 +27817,144 @@ nolink "false"
 .
 \end_layout
 
+\begin_layout Subsection
+enedecomp
+\end_layout
+
+\begin_layout LyX-Code
+enedecomp [<name>] [<mask>] [out <filename>]
+\end_layout
+
+\begin_layout LyX-Code
+          [ pme [cut <cutoff>] [dsumtol <dtol>] [ewcoeff <coeff>]
+\end_layout
+
+\begin_layout LyX-Code
+                [erfcdx <dx>] [skinnb <skinnb>] [ljswidth <width>]
+\end_layout
+
+\begin_layout LyX-Code
+                [order <order>] [nfft <nfft1>,<nfft2>,<nfft3>]
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+[<name>] Data set name.
+\end_layout
+
+\begin_layout Description
+[<mask>] Mask of atoms to calculate energy for.
+\end_layout
+
+\begin_layout Description
+[out
+\begin_inset space ~
+\end_inset
+
+<filename>] File to write results to.
+\end_layout
+
+\begin_layout Description
+[pme] Use particle mesh Ewald for electrostatics;
+ van der Waals energy will be calculated using a long-range correction for periodicity.
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+cut
+\begin_inset space ~
+\end_inset
+
+<cutoff> Direct space cutoff in Angstroms (default 8.0).
+\end_layout
+
+\begin_layout Description
+dsumtol
+\begin_inset space ~
+\end_inset
+
+<dtol> Direct sum tolerance (default 0.00001).
+ Used to determine Ewald coefficient.
+\end_layout
+
+\begin_layout Description
+ewcoeff
+\begin_inset space ~
+\end_inset
+
+<coeff> Ewald coefficient in 1/Ang.
+\end_layout
+
+\begin_layout Description
+erfcdx
+\begin_inset space ~
+\end_inset
+
+<dx> Spacing to use for the ERFC splines (default 0.0002 Ang.).
+\end_layout
+
+\begin_layout Description
+skinnb Used to determine pairlist atoms (added to cut,
+ so pairlist cutoff is 
+\series bold
+cut
+\series default
+ + 
+\series bold
+skinnb
+\series default
+);
+ included in order to maintain consistency with results from sander.
+\end_layout
+
+\begin_layout Description
+ljswidth
+\begin_inset space ~
+\end_inset
+
+<width> If specified,
+ use a force-switching form for the Lennard-Jones calculation from <cutoff>-<width> to <cutoff>.
+\end_layout
+
+\begin_layout Description
+order
+\begin_inset space ~
+\end_inset
+
+<order> Spline order for charges.
+\end_layout
+
+\begin_layout Description
+nfft
+\begin_inset space ~
+\end_inset
+
+<nfft1>,<nfft2>,<nfft3> Explicitly set the number of FFT grid points in each dimension.
+ Will be determined automatically from unit cell dimensions if not specified.
+\end_layout
+
+\end_deeper
+\begin_layout Standard
+DataSets created:
+\end_layout
+
+\begin_layout Description
+<name> Set containing atom index and the corresponding average energy over frames.
+\end_layout
+
+\end_deeper
+\begin_layout Standard
+Perform per-atom energy decomposition for selected atoms.
+ The energy is calculated for the entire system but only the energies for selected atoms will be reported.
+ The energy is composed of the regular bond,
+ angle,
+ torsion,
+ 1-4 nonbonded,
+ and nonbonded terms.
+ If 'pme' is specified the non-bonded terms will us PME for electrostatics and a long-range periodic correction for van der Waals,
+ otherwise a simple model with no cutoff will be used.
+\end_layout
+
 \begin_layout Subsection
 energy
 \end_layout
@@ -28154,7 +28321,7 @@ nfft
 \end_inset
 
 <nfft1>,<nfft2>,<nfft3> Explicitly set the number of FFT grid points in each dimension.
- Will be determined automatically if not specified.
+ Will be determined automatically from unit cell dimensions if not specified.
 \end_layout
 
 \begin_layout Description

src/Action_EneDecomp.cpp

@@ -0,0 +1,51 @@
+#include "Action_EneDecomp.h"
+#include "CpptrajStdio.h"
+
+// Action_EneDecomp::Help()
+void Action_EneDecomp::Help() const {
+  Cpptraj::Energy::EnergyDecomposer::HelpText(); 
+}
+
+// Action_EneDecomp::Init()
+Action::RetType Action_EneDecomp::Init(ArgList& actionArgs, ActionInit& init, int debugIn)
+{
+# ifdef MPI
+  trajComm_ = init.TrajComm();
+# endif
+  if (eneDecomp_.InitDecomposer( actionArgs, init.DSL(), init.DFL(), debugIn ))
+    return  Action::ERR;
+  mprintf("    ENEDECOMP: Decomposing energy for selected atoms.\n");
+  eneDecomp_.PrintOpts();
+
+  return Action::OK;
+}
+
+// Action_EneDecomp::Setup()
+Action::RetType Action_EneDecomp::Setup(ActionSetup& setup)
+{
+  int ret = eneDecomp_.SetupDecomposer( setup.Top(), setup.CoordInfo().TrajBox() );
+  if (ret == -1)
+    return Action::SKIP;
+  else if (ret == 1)
+    return Action::ERR;
+  return Action::OK;
+}
+
+// Action_EneDecomp::DoAction()
+Action::RetType Action_EneDecomp::DoAction(int frameNum, ActionFrame& frm)
+{
+  if (eneDecomp_.CalcEne( frm.Frm() ))
+    return Action::ERR;
+  return Action::OK;
+}
+
+#ifdef MPI
+int Action_EneDecomp::SyncAction() {
+  return eneDecomp_.ReduceToMaster(trajComm_);
+}
+#endif
+
+// Action_EneDecomp::Print()
+void Action_EneDecomp::Print() {
+  eneDecomp_.FinishCalc();
+}

src/Action_EneDecomp.h

@@ -0,0 +1,32 @@
+#ifndef INC_ACTION_ENEDECOMP_H
+#define INC_ACTION_ENEDECOMP_H
+#include "Action.h"
+#include "Energy/EnergyDecomposer.h"
+/// Used to decompose the pairwise additive energy of a system 
+class Action_EneDecomp : public Action {
+  public:
+    /// CONSTRUCTOR
+    Action_EneDecomp() {}
+    /// Allocator
+    DispatchObject* Alloc() const { return (DispatchObject*)new Action_EneDecomp(); }
+    /// Help text
+    void Help() const;
+  private:
+    /// Initialization
+    Action::RetType Init(ArgList&, ActionInit&, int);
+    /// Topology-based setup
+    Action::RetType Setup(ActionSetup&);
+    /// Do the action
+    Action::RetType DoAction(int, ActionFrame&);
+    /// Print results/finish calculations
+    void Print();
+#   ifdef MPI
+    int SyncAction();
+#   endif
+
+    Cpptraj::Energy::EnergyDecomposer eneDecomp_; ///< Do the actual decomposition
+#   ifdef MPI
+    Parallel::Comm trajComm_; ///< Across-trajectory communicator
+#   endif
+};
+#endif