Skip to content

Commit

Permalink
Add energy decomposition (with PME) (#1109)
Browse files Browse the repository at this point in the history
* Start EnergyDecomposer class

* Update dependencies

* Start init routine

* Add mask and options print

* Start energy decomp action

* Start adding the enedecomp command, start initialization

* Start doing setup

* Bond setup

* Setup output arrays

* Print atoms that energy will be saved for

* Just allocate an energy for each atom, will make bookkeeping easier

* Fix up the dataset dimension labels. Add calculation

* Add output file

* Add angle term. Need to sum before accumulation

* Decompose dihedrals

* LJ 6-12 energy template

* Create and use an LJ 6-12 template

* Start adding the 1-4 energy

* Add 1-4 elec

* Add debug for the 1-4 terms

* Start adding nonbond decomp

* Only do 1-4 calculation for NORMAL dihedrals

* Start trying to template the pairlist loop

* Add wrapper around Erfc

* Add Ewald adjust kernel

* Finish initial incarnation of the Ewald pairlist nonbond kernel. Make
ERFC non-const and put timing tracking directly into ErfcFxn

* Some cleanup

* Add code comments

* Change name

* Change class name

* Reorder

* Try as a separate class

* Remove for now

* Rename

* Finish renaming

* Remove unneeded depend.

* Rip out the new framework. Its really not at all ready for further
development until i come up with a plan

* Start a separate class for decomposing ewald energies. Bring back
ErfcFxn and LJswitch

* Start splitting out the Ewald parameters

* Update build files

* Just make a function template again

* Put in energy namespace and make erfc_func available

* Put Erfc and switching functions through the EwaldParams class

* Init class vars

* Add params for PME

* Second try for using pairlist template/engine model

* Do some Ewald setup

* Add routines to access Ewald parameters to the engine

* Start PME calc wrapper

* Add class for calculating long range VDW

* Use long range VDW calc

* PME calculation setup

* Add PME recip calc wrapper class

* Start using PME_Recip

* Allow PME_Recip to be used for coulomb or LJ

* Use COULOMB recip

* Add self energy calc

* Save recip coords in EwaldParams_PME so it can be used for multiple
recip objects. Start the nonbonded calc

* Finish up nonbonded pme calc with template

* Remove separate EwaldDecomp for now

* Add timing data

* Report pairlist timing data

* New test action for comparing the old and new PME methods

* Enable the new pmetest action (hidden)

* Ensure previous coords are cleared

* Add unit cell vector check to recip part

* Move LJ ewald coeff into a separate class

* Add C6 param and LJ self energy calcs

* Return the precalculated self energy

* First pass at the LJPME engine

* Use Ewald parameters

* Remove unneeded typedef

* Start the LJPME calc class

* Use the right variables

* Add Calc_LJPME to build

* Add Init routine

* Add Setup routine

* Add nonbonded calc routine

* Add timing routine

* Need to pass C6 params to LJPME recip calc, not charges

* Test new LJPME calc

* Report LJ PME recip timing in addition to regular recip timing

* Openmp parallelize the pairlist template loop

* Add timer to check that the recip calc itself is what takes the majority of the time, not setup (and it does). Comment it out.

* Rename

* const correctness

* Rename

* Finish renaming

* Start decomposable PME energy class

* Placeholder for decomp ewald engine

* Add engine placeholder

* Do engine init

* Add setup

* Add decomposed self energy

* Add to build

* Start incorporating ene decomp

* Self energy decomp seems to work

* Put lattice check in a separate routine

* Decomposable recip energy

* Decompose direct elec, vdw, and adjust energies

* Restore arrays needed for atom decomp of LR vdw correction

* Decompose the long range correction

* Actually save the energies back to EnergyDecomposer

* Add timing data. Hide some debug info.

* Put LJ PME vdw calc into a template

* Add LJPME exclude adjust kernel

* Remove old code

* Add Ewald LJPME decomp

* Create decomposed LJPME calc. Still need to do decomposed self6 energy.

* Add decomposed self6 calc

* Use decomposed self6 energies

* Update depends

* Doesnt appear like the recip term for the LJPME is working yet.

* Use cutoff from EwaldOptions instead of the engine

* Start Ecalc_Nonbond

* Start implementation

* Add pairlist setup

* Add pairlist setup to nonbond calc

* Start refactoring so that pairlist and exclusion array are not part of
the class

* Add .cpp.o rule so that subdir compile uses the correct flags

* Abstract base class for Ewald calcs

* Use abstract base class

* Use EwaldCalc ABC

* Add to build

* Add a nonbond kernel

* Add note about includes

* Add calc init and setup, and simple nonbond calc

* Use new nonbond calc

* Report pairlist timing

* Try adding a simple nonbond decomposable energy kernel

* Start ewald decomp ABC

* Fix includes

* Convert to EwaldCalc_Decomp class

* Convert decomposd LJPME to Ewald_Decomp class. Update depends

* Add energy decomp to Ecalc_Nonbond

* All ewald calcs need the pairlist

* Use Ecalc_Nonbond to do the decomp

* Start moving nfft and order to inside PME_Recip

* New class to hold recip coords

* Add routine to access coords array, update depends

* Remove files

* Move EwaldParams_PME Init, Setup, and coordsD routines to EwaldParams

* Use new EwaldParams

* Use EwaldParams

* Have recip print opts right after init

* Use new EwaldParams

* Update for new EwaldParams and PME_Recip

* Update depends

* Use EwaldParams intead of EwaldParams_PME. Update depends

* Start Ewald_Recip class

* Make 1 / sqrt(PI) constant available

* Add init and print routines

* Add setup

* Finish setup

* Finish regular recip calc

* Add timing

* Update docs

* No need to save ewald coeff, just calculate the factor in Init

* Start regular Ewald calc class

* Add setup

* Finish regular ewald

* Add EwaldCalc_Regular and Ewald_Recip to build

* Add regular ewald to Ecalc_Nonbond

* Add ability to test Ewald regular

* Add missing init

* Need to look for method == 4. Improve error message. Fix timing message.

* Fix reduction pragma

* Fix OpenMP compile

* Introduce new constant COULOMBFACTOR to replace all separate instances
of QFAC

* Start adding ifdefs

* Start splitting off the direct sum routines

* Start removing old Ewald

* These dont need to be separate yet

* Make vars available for GIST PME

* Use Ecalc_Nonbond

* Start converting to use EwaldParams

* Add PME_RecipParams class

* Use PME_RecipParams class

* Put PME_EwaldParams inside PME_Recip

* Dont call print inside the init

* Expose arrays needed by GIST_PME

* Finish converting GIST_PME, add VDW long range, recip, and adjust stuff.

* Add missing timing function

* Fix non-openmp compile

* Start adding LJPME OptType

* Use IsPmeType where appropriate. Ensure GIST PME is only used with PME options for now.

* Energy decomp with ljpme not yet allowed

* Allow LJPME to be directly set

* Remove debug message

* Make options printout more clear

* Make sure printed Ewald options line up

* Make print private

* Clean up printout of regular ewald recip opts

* Since fftw3 is always on now by default, change -fftw3 to -pubfft in case user wants to force using pubfft.

* PME recip needs LIBPME

* Add include of base Ewald decomp class when no LIBPME since no decomp for regular Ewald yet.

* Add files for enedecomp test

* Enable enedecomp test

* Remove all references to PairListLoop and associated files

* Update dependencies

* Add help text

* Dont repeat the same timing data. Should clean that up eventually

* Use Ecalc_Nonbond interal timer only when it is active

* Hide some debug info

* Remove old code

* Hide some debug behind ifdef

* Fix up info output

* Hide debug info behind ifdef

* Version 6.29.4. Revision bump for addition of the 'enedecomp' action and
rewrite of the Ewald classes.

* Add documentation on how to use the pairlist template

* Add longer test

* Add enedecomp to manual

* Put functionality to reduce an array of Stats<T> into a separate class
to be used by PairDist and EnergyDecomposer

* Use Stats_Reduce

* Ensure decomposed energy results are synced

* Add citation to code docs

* Update manual PDF and checksum

* Fix non-MPI compile (oops)

---------

Co-authored-by: Daniel R. Roe <[email protected]>
  • Loading branch information
drroe and Daniel R. Roe authored Oct 18, 2024
1 parent 8d1a163 commit addb723
Show file tree
Hide file tree
Showing 94 changed files with 7,465 additions and 1,446 deletions.
5 changes: 3 additions & 2 deletions configure
Original file line number Diff line number Diff line change
Expand Up @@ -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."
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 ;;
Expand Down
5 changes: 5 additions & 0 deletions devtools/CreateMake.sh
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Binary file modified doc/CpptrajManual.pdf
Binary file not shown.
2 changes: 1 addition & 1 deletion doc/DocumentChecksums.txt
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
b37726e7a841f6fc695ecd7fb040ffbf CpptrajDevelopmentGuide.lyx
736ee725ca9cc67d349f67d714f41c20 cpptraj.lyx
151841b1d3283890c425aaae5a923e7f cpptraj.lyx
07c4039e732fc2eb1df0c1e0863cb949 CpptrajManual.lyx
171 changes: 169 additions & 2 deletions doc/cpptraj.lyx
Original file line number Diff line number Diff line change
Expand Up @@ -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%">
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
51 changes: 51 additions & 0 deletions src/Action_EneDecomp.cpp
Original file line number Diff line number Diff line change
@@ -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();
}
32 changes: 32 additions & 0 deletions src/Action_EneDecomp.h
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit addb723

Please sign in to comment.