Merge pull request #412 from michellab/feature_pmefep_2023

Support for PME in SOMD FEP [ci skip]
michellab · Nov 16, 2022 · 13e5112 · 13e5112
2 parents 485c0c6 + 3d5296b
commit 13e5112
Show file tree

Hide file tree

Showing 20 changed files with 19,430 additions and 769 deletions.
diff --git a/corelib/src/libs/SireMove/CMakeLists.txt b/corelib/src/libs/SireMove/CMakeLists.txt
@@ -133,12 +133,14 @@ set ( SIREMOVE_OPENMM_HEADERS
       openmmmdintegrator.h
       openmmfrenergydt.h
       openmmfrenergyst.h
+      openmmpmefep.h
     )
 
 set ( SIREMOVE_OPENMM_SOURCES
       openmmmdintegrator.cpp
       openmmfrenergydt.cpp
       openmmfrenergyst.cpp
+      openmmpmefep.cpp
     )
 
 if ( ${SIRE_USE_OPENMM} )

diff --git a/corelib/src/libs/SireMove/openmmfrenergyst.cpp b/corelib/src/libs/SireMove/openmmfrenergyst.cpp
diff --git a/corelib/src/libs/SireMove/openmmfrenergyst.h b/corelib/src/libs/SireMove/openmmfrenergyst.h
@@ -52,8 +52,8 @@ SIREMOVE_EXPORT QDataStream& operator>>(QDataStream&, SireMove::OpenMMFrEnergyST
 
 namespace SireMove {
 
-    /** This class implements single topology a free energy method using OpenMM. 
- 
+    /** This class implements single topology a free energy method using OpenMM.
+
         @author Julien Michel,Gaetano Calabro and Antonia Mey
      */
     class SIREMOVE_EXPORT OpenMMFrEnergyST
@@ -108,7 +108,7 @@ namespace SireMove {
 
 	QString getCombiningRules(void);
 	void setCombiningRules(QString);
-	
+
         QString getCutoffType(void);
         void setCutoffType(QString);
 
@@ -200,17 +200,19 @@ namespace SireMove {
         int getRandomSeed(void);
         void setRandomSeed(int);
 
+	void setDebug(bool);
+
     private:
         void createContext(IntegratorWorkspace &workspace,
                 SireUnits::Dimension::Time timestep);
         void destroyContext();
-        void updateBoxDimensions(OpenMM::State &state_openmm, 
-        QVector<QVector<Vector>> &buffered_dimensions, bool Debug, 
-        AtomicVelocityWorkspace &ws);
-        
+        void updateBoxDimensions(OpenMM::State &state_openmm,
+				 QVector<QVector<Vector>> &buffered_dimensions,
+				 AtomicVelocityWorkspace &ws);
+
         double getPotentialEnergyAtLambda(double lambda);
         void updateOpenMMContextLambda(double lambda);
-        boost::tuples::tuple<double, double, double> calculateGradient(double increment_plus, 
+        boost::tuples::tuple<double, double, double> calculateGradient(double increment_plus,
         double increment_minus, double potential_energy_lambda, double beta);
         QVector<double> computeReducedPerturbedEnergies(double);
         void emptyContainers(void);
@@ -282,7 +284,7 @@ namespace SireMove {
         QVector<double> pot_energies;
 
         QVector<double> forward_Metropolis;
-        
+
         QVector<double> backward_Metropolis;
 
         QVector<QVector <double> > reduced_perturbed_energies;

diff --git a/corelib/src/libs/SireMove/openmmmdintegrator.cpp b/corelib/src/libs/SireMove/openmmmdintegrator.cpp
@@ -426,11 +426,11 @@ void OpenMMMDIntegrator::initialise()
         flag_cutoff = CUTOFFPERIODIC;
     else if (CutoffType == "ewald")
         flag_cutoff = EWALD;
-    else if (CutoffType == "pme")
+    else if (CutoffType == "PME")
         flag_cutoff = PME;
     else
         throw SireError::program_bug(QObject::tr(
-                                                 "The CutOff method has not been specified. Possible choices: nocutoff, cutoffnonperiodic, cutoffperiodic,ewal,pme"), CODELOC);
+                                                 "The CutOff method has not been specified. Possible choices: nocutoff, cutoffnonperiodic, cutoffperiodic,ewal,PME"), CODELOC);
 
     if (Debug)
         qDebug() << "\nCutoffType = " << CutoffType << "\n";
@@ -1485,8 +1485,7 @@ System OpenMMMDIntegrator::equilibrateSystem(System &system,
 
     const int nmols = ws.nMolecules();
     int k = 0;
-    double Ekin_sire = 0;
-    double Ekin_openmm = 0;
+
     for (int i = 0; i < nmols; i++)
     {
 
@@ -1713,7 +1712,7 @@ void OpenMMMDIntegrator::integrate(IntegratorWorkspace &workspace, const Symbol
 
     int k = 0;
     double Ekin_openmm=0;
-    double Ekin_sire=0;
+
     for (int i = 0; i < nmols; i++)
     {
         Vector *sire_coords = ws.coordsArray(i);