Wrap HOOMD Force Computes as collective variables

hooks/hoomd/Driver.cpp

@@ -48,8 +48,11 @@ namespace SSAGES
         py::object exec_conf_obj = py::cast(exec_conf_);
         hoomd_module.attr("context").attr("exec_conf") = exec_conf_obj;
 
+        // Evaluate in scope of main module
+        py::object scope = py::module::import("__main__").attr("__dict__");
+
         // Evaluate the user script file
-        py::eval_file(rh_->GetInput().c_str());
+        py::eval_file(rh_->GetInput().c_str(),scope);
         py::object context = hoomd_module.attr("context").attr("current");
 
         // Get the current C++ SystemDefinition and update the HOOMDHook
@@ -73,6 +76,7 @@ namespace SSAGES
 
         // Run HOOMD
         hook_->update(0);  // this is a hack for initialization
+        hook_->setPyScope(scope);
         hook_->PreSimulationHook();
         hoomd_module.attr("run")(run_steps_);
         hook_->PostSimulationHook();
@@ -117,8 +121,9 @@ namespace SSAGES
             execution_mode = ExecutionConfiguration::executionMode::GPU;
             }
 
+        std::vector<int> gpu_id(1,-1);
         auto exec_conf = std::make_shared<ExecutionConfiguration>(
-          execution_mode, -1, false, false,
+          execution_mode, gpu_id, false, false,
           std::shared_ptr<Messenger>(), 0, rh->GetLocalComm());
 
         auto hook = std::make_shared<HOOMDHook>();

hooks/hoomd/HOOMDHook.cpp

@@ -6,7 +6,7 @@
 using namespace SSAGES;
 
 HOOMDHook::HOOMDHook() :
-Hook(), HalfStepHook()
+Hook(), HalfStepHook(), timestep_(0)
 {
   std::cerr << "Installing SSAGES HOOMD-blue hook." << std::endl;
 }

hooks/hoomd/HOOMDHook.h

@@ -2,11 +2,14 @@
 #define SSAGES_HOOMD_HOOK_H
 
 #include "Hook.h"
+#include "CVs/CVManager.h"
+
 #include "hoomd/ExecutionConfiguration.h"
 #include "hoomd/SystemDefinition.h"
 #include "hoomd/HalfStepHook.h"
 #include "hoomd/ComputeThermo.h"
 
+#include "hoomd/extern/pybind/include/pybind11/pybind11.h"
 
 // SSAGES Hook class for HOOMD
 class HOOMDHook : public SSAGES::Hook, public HalfStepHook
@@ -21,6 +24,9 @@ class HOOMDHook : public SSAGES::Hook, public HalfStepHook
         // HOOMD ComputeThermo
         std::shared_ptr<ComputeThermo> thermo_;
 
+        // Python interpreter scope
+        pybind11::object scope_;
+
     protected:
         // Implementation of the SyncToEngine interface.
         void SyncToEngine() override;
@@ -31,20 +37,61 @@ class HOOMDHook : public SSAGES::Hook, public HalfStepHook
     public:
         HOOMDHook();
 
+        //! Pre-simulation hook
+        /*!
+         * This should be called at the appropriate
+         * time by the Hook implementation.
+         */
+        virtual void PreSimulationHook() override
+        {
+            // call base class method
+            Hook::PreSimulationHook();
+
+            // Set simulation engine hook on CVs
+            for(auto& cv : cvmanager_->GetCVs())
+            {
+                cv->SetHook(*this);
+            }
+        }
+
         // Sets SystemDefinition to enable particle data access
-        void setSystemDefinition(std::shared_ptr<SystemDefinition> sysdef)
+        void setSystemDefinition(std::shared_ptr<SystemDefinition> sysdef) override
         {
             sysdef_ = sysdef;
         }
 
+        // Get the SystemDefinition
+        std::shared_ptr<SystemDefinition> getSystemDefinition()
+            {
+            return sysdef_;
+            }
+
         // Sets ComputeThermo to enable temperature and energy computation
         void setComputeThermo(std::shared_ptr<ComputeThermo> thermo)
         {
             thermo_ = thermo;
         }
 
         // Synchronize snapshot with SSAGES after computing forces
-        void update(unsigned int timestep);
+        void update(unsigned int timestep) override;
+
+        // Return the current time step
+        unsigned int getTimeStep()
+            {
+            return timestep_;
+            }
+
+        // Set the python interpreter scope
+        void setPyScope(pybind11::object scope)
+            {
+            scope_ = scope;
+            }
+
+        // Set the python interpreter scope
+        pybind11::object getPyScope()
+            {
+            return scope_;
+            }
 
         ~HOOMDHook();
 };

hooks/hoomd/HOOMDWrapCV.h

@@ -0,0 +1,172 @@
+/**
+ * This file is part of
+ * SSAGES - Suite for Advanced Generalized Ensemble Simulations
+ *
+ * Copyright 2016 Yamil Colon <yamilcolon2015@u.northwestern.edu>
+ *                Hythem Sidky <hsidky@nd.edu>
+ *
+ * SSAGES is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * SSAGES is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with SSAGES.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#pragma once 
+
+#include "CVs/CollectiveVariable.h"
+#include "Validator/ObjectRequirement.h"
+#include "Drivers/DriverException.h"
+#include "Snapshot.h"
+#include "schema.h"
+
+#include "hoomd/ForceCompute.h"
+
+namespace SSAGES
+{
+    //! Collective variable to calculate angle.
+    class HOOMDWrapCV : public CollectiveVariable
+    {
+    private:
+        //! The ForceCompute we wrap
+        std::shared_ptr<ForceCompute> fc_;
+
+        //! A reference to the HOOMD hook
+        HOOMDHook& hook_;
+
+    public:
+
+        //! Constructor.
+        /*!
+         * \param fc The force compute to wrap
+         *
+         * Wrap a HOOMD ForceCompute into a CollectiveVariable
+         *
+         * \todo Bounds needs to be an input and periodic boundary conditions
+         */
+        HOOMDWrapCV(std::shared_ptr<ForceCompute> fc);
+        fc_(fc)
+        { }
+
+        //! Initialize necessary variables.
+        /*!
+         * \param snapshot Current simulation snapshot.
+         */
+        void Initialize(const Snapshot& snapshot) override
+        {
+        #if 0
+        using std::to_string;
+
+        std::vector<int> found;
+        snapshot.GetLocalIndices(atomids_, &found);
+        int nfound = found.size();
+        MPI_Allreduce(MPI_IN_PLACE, &nfound, 1, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
+
+        if(nfound != 3)
+            throw BuildException({
+                "TorsionalCV: Expected to find " + 
+                to_string(3) + 
+                " atoms, but only found " + 
+                to_string(nfound) + "."
+            });	
+        #endif
+        }
+
+        //! This sets a reference to the simulation engine hook.
+        /*!
+         * This method is called during initialization and can be used by the CV
+         * to access engine-specific data structures.
+         */
+        virtual void SetHook(const HOOMDHook& hook)
+            {
+            hook_ = hook;
+            }
+
+        //! Evaluate the CV.
+        /*!
+        * \param snapshot Current simulation snapshot.
+        */
+        void Evaluate(const Snapshot& snapshot) override
+        {
+            // evaluate the underlying ForceCompute
+            fc_->compute(hook_.getTimeStep());
+
+            // reduce CV
+            val_ = fc->calcEnergySum();
+        }
+
+        //! Compute the particle forces and virial by scaling the gradient with a constant
+        /*!
+         * \param bias The constant scaling factor to apply to all particle forces
+         * \param snapshot The snapshot that contains the forces
+         *
+         */
+        virtual void ApplyBias(double bias, const Snapshot& snapshot) const
+        {
+            // evaluate the underlying ForceCompute (only if necessary)
+            fc_->compute(hook_.getTimeStep());
+
+            ArrayHandle<Scalar4> h_force(fc_->getForce(), access_location::host, access_mode::readwrite);
+            ArrayHandle<Scalar> h_virial(fc_->getVirial(), access_location::host, access_mode::readwrite);
+            ArrayHandle<Scalar4> h_torque(fc_->getTorque(), access_location::host, access_mode::readwrite);
+
+            unsigned int virial_pitch = fc_->getVirial().getPitch();
+
+            auto pdata = hook_->getSystemDefinition()->getParticleData();
+            unsigned int n = pdata->getN();
+
+            // multiply particle forces with bias
+            for (unsigned int i = 0; i < n; ++i)
+                h_force.data[i] *= bias;
+
+            for (unsigned int j = 0; j < 6; ++j)
+                for (unsigned int i = 0; i < n; ++i)
+                    h_virial.data[j*virial_pitch + i] *= bias;
+
+            #if 0
+            // currently cannot set the external virial
+            for (unsigned int j = 0; j < 6; ++j)
+                fc_->setExternalVirial(j,fc_->getExternalVirial(j)*bias);
+            #endif
+        }
+
+        //! \copydoc CollectiveVariable::BuildCV()
+        static HOOMDWrapCV* Build(const Json::Value& json, const std::string& path)
+        {
+        Json::ObjectRequirement validator;
+        Json::Value schema;
+        Json::CharReaderBuilder rbuilder;
+        Json::CharReader* reader = rbuilder.newCharReader();
+
+        reader->parse(JsonSchema::HOOMDWrapCV.c_str(),
+                      JsonSchema::HOOMDWrapCV.c_str() + JsonSchema::HOOMDWrapCV.size(),
+                      &schema, NULL);
+        validator.Parse(schema, path);
+
+        // Validate inputs.
+        validator.Validate(json, path);
+        if(validator.HasErrors())
+            throw BuildException(validator.GetErrors());
+
+        std::vector<int> atomids;
+        auto force = json["force"].asString();
+
+        std::shared_ptr<ForceCompute> fc = context.attr(force).attr("cpp_force");
+        return new HOOMDWrapCV(fc);
+        }
+
+        //! Returns true if this CV can modify the particle forces in the snapshot
+        virtual bool ModifiesParticleForces() const
+        {
+        return false;
+        }
+
+
+    };
+}

schema/CVs/HOOMD.cv.json

@@ -0,0 +1,23 @@
+{
+    "type" : "object",
+    "varname" : "HOOMDCV",
+    "properties" : {
+        "type" : {
+            "type" : "string",
+            "enum" : ["HOOMDCV"]
+        },
+        "force" : {
+            "type" : "string"
+        },
+        "bounds" : {
+        "type" : "array",
+        "minItems" : 2,
+        "maxItems" : 2,
+        "items" : {
+            "type" : "number"
+            }
+        }
+    },
+    "required" : ["type", "force"],
+    "additionalProperties" : false
+}

src/CVs/CollectiveVariable.cpp

@@ -36,6 +36,10 @@
 #include "json/json.h"
 #include <stdexcept>
 
+#ifdef ENABLE_HOOMD
+#include "../../hooks/hoomd/HOOMDWrapCV.h"
+#endif
+
 namespace SSAGES
 {
 	CollectiveVariable* CollectiveVariable::BuildCV(const Json::Value &json, const std::string& path)
@@ -67,6 +71,12 @@ namespace SSAGES
 			return ParallelBetaRMSDCV::Build(json, path);
 		else if (type == "AntiBetaRMSD")
 			return AntiBetaRMSDCV::Build(json, path);
+        else if (type == "HOOMDCV")
+        #ifdef ENABLE_HOOMD
+            return HOOMDWrapCV::Build(json, path);
+        #else
+            throw std::invalid_argument(path + ": HOOMDCV not available. Compile with HOOMD Hook.");
+        #endif
 		else
 			throw std::invalid_argument(path + ": Unknown CV type specified.");
 	}