Commit 5be64eae authored by Christopher Bruns's avatar Christopher Bruns
Browse files

Remove NML plugin subtree.

parent d20a82d1
...@@ -353,12 +353,6 @@ IF(OPENMM_BUILD_FREE_ENERGY_PLUGIN) ...@@ -353,12 +353,6 @@ IF(OPENMM_BUILD_FREE_ENERGY_PLUGIN)
ADD_SUBDIRECTORY(plugins/freeEnergy) ADD_SUBDIRECTORY(plugins/freeEnergy)
ENDIF(OPENMM_BUILD_FREE_ENERGY_PLUGIN) ENDIF(OPENMM_BUILD_FREE_ENERGY_PLUGIN)
# Normal mode langevin plugin
set(OPENMM_BUILD_NML_PLUGIN OFF CACHE BOOL "Build normal-mode langevin plugin")
if(OPENMM_BUILD_NML_PLUGIN)
add_subdirectory(plugins/normalModeLangevin)
endif(OPENMM_BUILD_NML_PLUGIN)
INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${SHARED_TARGET}) INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${SHARED_TARGET})
INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${STATIC_TARGET}) INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${STATIC_TARGET})
FILE(GLOB CORE_HEADERS include/*.h */include/*.h) FILE(GLOB CORE_HEADERS include/*.h */include/*.h)
......
#/* -------------------------------------------------------------------------- *
# * OpenMM *
# * -------------------------------------------------------------------------- *
# * This is part of the OpenMM molecular simulation toolkit originating from *
# * Simbios, the NIH National Center for Physics-Based Simulation of *
# * Biological Structures at Stanford, funded under the NIH Roadmap for *
# * Medical Research, grant U54 GM072970. See https://simtk.org. *
# * *
# * Portions copyright (c) 2009 Stanford University and the Authors. *
# * Authors: Chris Sweet, Christopher Bruns *
# * Contributors: *
# * *
# * Permission is hereby granted, free of charge, to any person obtaining a *
# * copy of this software and associated documentation files (the "Software"), *
# * to deal in the Software without restriction, including without limitation *
# * the rights to use, copy, modify, merge, publish, distribute, sublicense, *
# * and/or sell copies of the Software, and to permit persons to whom the *
# * Software is furnished to do so, subject to the following conditions: *
# * *
# * The above copyright notice and this permission notice shall be included in *
# * all copies or substantial portions of the Software. *
# * *
# * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
# * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
# * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
# * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
# * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
# * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
# * USE OR OTHER DEALINGS IN THE SOFTWARE. *
# * -------------------------------------------------------------------------- */
# Create Normal Mode Langevin (NML) plugin for OpenMM
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include)
file(GLOB nml_headers "include/*.h")
file(GLOB nml_sources "src/*.cpp")
# Reference platform implementation
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/platforms/reference/include)
include_directories(${CMAKE_SOURCE_DIR}/platforms/reference/src)
file(GLOB headers "${CMAKE_CURRENT_SOURCE_DIR}/platforms/reference/include/*.h")
file(GLOB sources "${CMAKE_CURRENT_SOURCE_DIR}/platforms/reference/src/*.cpp")
set(nml_headers ${nml_headers} ${headers})
set(nml_sources ${nml_sources} ${sources})
set(NML_PLUGIN NormalModeLangevin)
# On Unix or Cygwin we have to add the debug suffix manually
if(UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
set(NML_PLUGIN ${NML_PLUGIN}_d)
endif(UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
add_library(${NML_PLUGIN} SHARED ${nml_sources} ${nml_headers})
set_target_properties(${NML_PLUGIN} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY")
target_link_libraries( ${NML_PLUGIN} ${SHARED_TARGET} )
if(BUILD_TESTING)
# Copy plugin to test_plugin_dir
set(test_plugin_dir "${CMAKE_BINARY_DIR}/test_plugin_dir")
file(MAKE_DIRECTORY "${test_plugin_dir}")
# On Windows we need to copy the correct Release/Debug plugin for testing
if(MSVC)
add_custom_command(TARGET ${NML_PLUGIN} POST_BUILD
DEPENDS $(TargetPath)
COMMAND "${CMAKE_COMMAND}"
ARGS -E copy \"$\(TargetPath\)\" \"${test_plugin_dir}\"
COMMENT "Copying normal mode langevin plugin for testing (WIN32)")
else(MSVC)
get_target_property(old_loc ${NML_PLUGIN} LOCATION)
add_custom_command(TARGET ${NML_PLUGIN} POST_BUILD
COMMAND "${CMAKE_COMMAND}"
ARGS -E copy ${old_loc} ${test_plugin_dir}
COMMENT "Copying normal mode langevin plugin for testing")
endif(MSVC)
add_subdirectory(test)
endif(BUILD_TESTING)
if(WIN32)
# install DLL but not LIB
install(TARGETS ${NML_PLUGIN} RUNTIME DESTINATION lib/plugins)
else(WIN32)
install(TARGETS ${NML_PLUGIN} LIBRARY DESTINATION lib/plugins)
endif(WIN32)
/* -*- c++ -*- */
#ifndef EIGENVECTORINFO_H
#define EIGENVECTORINFO_H
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include <string>
#include <iostream>
#include <vector>
using std::cout;
using std::endl;
namespace ProtoMol
{
//_________________________________________________________________XYZ
/**
* Container holding coordinates/Vector3D and names
*/
struct EigenvectorInfo {
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Constructors, destructors, assignment
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
EigenvectorInfo() : myEigenvectorLength ( 0 ), myNumEigenvectors ( 0 ), myNumUsedEigenvectors( 0 ),
myEigenvectors ( 0 ), myOrigCEigval( 0.0 ), myNewCEigval( 0.0 ), myOrigTimestep( 0.0 ),
reDiagonalize( false ),
mySingleEigs( 0 ), myEigVecChanged( true ), myMinimumLimit( 0.5 ), currentMode( -1 ) {};
EigenvectorInfo( unsigned int n, unsigned int m ) : myEigenvectorLength( n ), myNumEigenvectors( m ), myNumUsedEigenvectors( m ),
myMaxEigenvalue( 0.0 ), myEigenvectors ( new double[n * m * 3] ),
myOrigCEigval( 0.0 ), myNewCEigval( 0.0 ), myOrigTimestep( 0.0 ), reDiagonalize( false ),
mySingleEigs( 0 ),
myEigVecChanged( true ), myMinimumLimit( 0.5 ), currentMode( -1 ) {}
~EigenvectorInfo() {
if ( myEigenvectors ) {
delete [] myEigenvectors;
myEigenvectors = 0;
}
if ( mySingleEigs ) {
delete [] mySingleEigs;
mySingleEigs = 0;
}
};
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// New methods of class EigenvectorInfo
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
bool initializeEigenvectors() {
try {
myEigenvectors = new double[myEigenvectorLength * myNumEigenvectors * 3];
} catch ( std::bad_alloc& ) {
return false;
}
return true;
}
float* getFloatEigPointer() {
//no eigenvectors assigned?
if(!myEigenvectors) return (float*)0;
const unsigned int arrayLen = myEigenvectorLength * myNumEigenvectors * 3;
//assign storage if required
if(!mySingleEigs) {
try {
mySingleEigs = new float[arrayLen];
} catch ( std::bad_alloc& ) {
return (float*)0;
}
}
//update values if double array updated
if(myEigVecChanged) {
for( unsigned int i=0; i<arrayLen; i++)
mySingleEigs[i] = (float)myEigenvectors[i];
myEigVecChanged = false;
}
return mySingleEigs;
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// My data members
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Eigenvector information
unsigned int myEigenvectorLength;
unsigned int myNumEigenvectors;
unsigned int myNumUsedEigenvectors;
double myMaxEigenvalue;
double *myEigenvectors;
//Current and original max subspace eigenvalue,
//for adaptive timestep
double myOrigCEigval, myNewCEigval;
double myOrigTimestep;
//re-diagonalization flag
bool reDiagonalize;
//OpenMM single precision interface
float *mySingleEigs;
bool myEigVecChanged;
double myMinimumLimit;
//Analytic integrator
std::vector< double > myEigenvalues;
int currentMode;
};
}
#endif /* EIGENVECTORINFO_H */
#ifndef INTEGRATE_NML_STEP_KERNEL_H
#define INTEGRATE_NML_STEP_KERNEL_H
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/KernelImpl.h"
#include "openmm/System.h"
#include "NMLIntegrator.h"
namespace OpenMM {
/**
* This kernel is invoked by NMLIntegrator to take one time step.
*/
class IntegrateNMLStepKernel : public KernelImpl {
public:
static std::string Name() {
return "IntegrateNMLStep";
}
IntegrateNMLStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param integrator the NMLIntegrator this kernel will be used for
*/
virtual void initialize(const System& system, const NMLIntegrator& integrator) = 0;
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the NMLIntegrator this kernel is being used for
*/
virtual void execute(ContextImpl& context, const NMLIntegrator& integrator, const double currentPE, const int stepType) = 0;
};
} /* namespace OpenMM */
#endif /* INTEGRATE_NML_STEP_KERNEL_H */
#ifndef OPENMM_NMLIntegrator_H_
#define OPENMM_NMLIntegrator_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/Integrator.h"
#include "openmm/Kernel.h"
#include "openmm/internal/windowsExport.h"
#include "EigenvectorInfo.h"
namespace OpenMM {
/**
* This is an Integrator which simulates a System using Langevin dynamics.
*/
class OPENMM_EXPORT NMLIntegrator : public Integrator {
public:
/**
* Create a NMLIntegrator.
*
* @param temperature the temperature of the heat bath (in Kelvin)
* @param frictionCoeff the friction coefficient which couples the system to the heat bath
* @param stepSize the step size with which to integrator the system (in picoseconds)
*/
NMLIntegrator(double temperature, double frictionCoeff, double stepSize, ProtoMol::EigenvectorInfo* projectionVectorInfo );
/**
* Get the temperature of the heat bath (in Kelvin).
*/
double getTemperature() const {
return temperature;
}
/**
* Set the temperature of the heat bath (in Kelvin).
*/
void setTemperature(double temp) {
temperature = temp;
}
/**
* Get the friction coefficient which determines how strongly the system is coupled to
* the heat bath.
*/
double getFriction() const {
return friction;
}
/**
* Set the friction coefficient which determines how strongly the system is coupled to
* the heat bath.
*/
void setFriction(double coeff) {
friction = coeff;
}
void setProjectionVectorInfo( ProtoMol::EigenvectorInfo* inVectorInfo ){
mProjectionVectorInfo = inVectorInfo;
}
ProtoMol::EigenvectorInfo* getProjectionVectorInfo() const{
return mProjectionVectorInfo;
}
unsigned int getNumProjectionVectors() const{
//return mProjectionVectorInfo->myNumEigenvectors;
return mProjectionVectorInfo->myNumUsedEigenvectors;
}
double getMinimumLimit() const{
return mProjectionVectorInfo->myMinimumLimit * 4.184; //Kcal->KJ
}
bool getProjVecChanged() const{
return mProjectionVectorInfo->myEigVecChanged;
}
double* getProjectionVectors() const{
return mProjectionVectorInfo->myEigenvectors;
}
void setProjVecChanged(bool inChanged){
mProjectionVectorInfo->myEigVecChanged = inChanged;
}
double getMaxEigenvalue() const{
//return max eigenvalue (f^2) in 1/ps from 1/charmm
return mProjectionVectorInfo->myMaxEigenvalue *
(1000.0 / 48.88821290839616) * (1000.0 / 48.88821290839616) //(1/charmm to 1/ps)^2
* 100.0; //(Angstron to nm [inverse])^2
}
/**
* Advance a simulation through time by taking a series of time steps.
*
* @param steps the number of time steps to take
*/
void step(int steps);
//Minimizer
void minimize(int maxsteps);
protected:
/**
* This will be called by the OpenMMContext when it is created. It informs the Integrator
* of what context it will be integrating, and gives it a chance to do any necessary initialization.
* It will also get called again if the application calls reinitialize() on the OpenMMContext.
*/
void initialize(ContextImpl& context);
/**
* Get the names of all Kernels used by this Integrator.
*/
std::vector<std::string> getKernelNames();
private:
ProtoMol::EigenvectorInfo* mProjectionVectorInfo;
double temperature, friction;
ContextImpl* context;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_NMLIntegrator_H_*/
#ifndef NML_KERNELS_H
#define NML_KERNELS_H
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Christopher Bruns *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "IntegrateNMLStepKernel.h"
#endif /* NML_KERNELS_H */
#ifndef OPENMM_NORMAL_MODE_LANGEVIN_H_
#define OPENMM_NORMAL_MODE_LANGEVIN_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Christopher Bruns *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program 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 Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
namespace OpenMM {
} // namespace OpenMM
#endif /*OPENMM_NORMAL_MODE_LANGEVIN_H_*/
#/* -------------------------------------------------------------------------- *
# * OpenMM *
# * -------------------------------------------------------------------------- *
# * This is part of the OpenMM molecular simulation toolkit originating from *
# * Simbios, the NIH National Center for Physics-Based Simulation of *
# * Biological Structures at Stanford, funded under the NIH Roadmap for *
# * Medical Research, grant U54 GM072970. See https://simtk.org. *
# * *
# * Portions copyright (c) 2009 Stanford University and the Authors. *
# * Authors: Chris Sweet, Christopher Bruns *
# * Contributors: *
# * *
# * Permission is hereby granted, free of charge, to any person obtaining a *
# * copy of this software and associated documentation files (the "Software"), *
# * to deal in the Software without restriction, including without limitation *
# * the rights to use, copy, modify, merge, publish, distribute, sublicense, *
# * and/or sell copies of the Software, and to permit persons to whom the *
# * Software is furnished to do so, subject to the following conditions: *
# * *
# * The above copyright notice and this permission notice shall be included in *
# * all copies or substantial portions of the Software. *
# * *
# * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
# * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
# * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
# * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
# * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
# * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
# * USE OR OTHER DEALINGS IN THE SOFTWARE. *
# * -------------------------------------------------------------------------- */
#ifndef REFERENCE_INTEGRATE_NML_STEP_KERNEL_H
#define REFERENCE_INTEGRATE_NML_STEP_KERNEL_H
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferencePlatform.h"
#include "ReferenceNMLDynamics.h"
#include "IntegrateNMLStepKernel.h"
namespace OpenMM {
/**
* This kernel is invoked by NMLIntegrator to take one time step.
*/
class ReferenceIntegrateNMLStepKernel : public IntegrateNMLStepKernel {
public:
ReferenceIntegrateNMLStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateNMLStepKernel(name, platform),
data(data), dynamics(0), constraints(0), masses(0), constraintDistances(0), constraintIndices(0), projectionVectors(0) {
}
~ReferenceIntegrateNMLStepKernel();
/**
* Initialize the kernel, setting up the particle masses.
*
* @param system the System this kernel will be applied to
* @param integrator the NMLIntegrator this kernel will be used for
*/
void initialize(const System& system, const NMLIntegrator& integrator);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the NMLIntegrator this kernel is being used for
*/
void execute(ContextImpl& context, const NMLIntegrator& integrator, const double currentPE, const int stepType);
private:
ReferencePlatform::PlatformData& data;
ReferenceNMLDynamics* dynamics;
ReferenceConstraintAlgorithm* constraints;
RealOpenMM* masses;
RealOpenMM* constraintDistances;
int** constraintIndices;
int numConstraints;
double prevTemp, prevFriction, prevStepSize;
//double prevTemp, prevFriction, prevErrorTol;
RealOpenMM* projectionVectors;
};
} /* namespace OpenMM */
#endif /* REFERENCE_INTEGRATE_NML_STEP_KERNEL_H */
#ifndef __ReferenceNMLDynamics_H__
#define __ReferenceNMLDynamics_H__
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns, Pande Group *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include "SimTKReference/ReferenceDynamics.h"
// ---------------------------------------------------------------------------------------
class ReferenceNMLDynamics : public ReferenceDynamics {
private:
//dt, myGamma, fdt, vdt, ndt, sqrtFCoverM from Langevin Leapfrog
//enum FixedParameters { GDT, EPH, EMH, EP, EM, B, C, D, V, X, Yv, Yx, MaxFixedParameters };
enum FixedParameters { DT, GAMMA, FDT, VDT, NDT, SQRTFCOVERM, MaxFixedParameters };
enum TwoDArrayIndicies { X2D, V2D, OldV, xPrime2D, vPrime2D, Max2DArrays };
enum OneDArrayIndicies { InverseMasses, Max1DArrays };
RealOpenMM _tau;
RealOpenMM _fixedParameters[MaxFixedParameters];
RealOpenMM* _projectionVectors;
unsigned int _numProjectionVectors;
RealOpenMM _minimumLimit, _maxEig;
RealOpenMM lastPE, lastSlope;
/**---------------------------------------------------------------------------------------
Set fixed values
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int _setFixedParameters( void );
public:
/**---------------------------------------------------------------------------------------
Constructor
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param tau viscosity
@param temperature temperature
--------------------------------------------------------------------------------------- */
ReferenceNMLDynamics( int numberOfAtoms, RealOpenMM deltaT, RealOpenMM tau, RealOpenMM temperature,
RealOpenMM* projectionVectors, unsigned int numProjectionVectors,
RealOpenMM minimumLimit, RealOpenMM maxEig
);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceNMLDynamics( );
/**---------------------------------------------------------------------------------------
Get tau
@return tau
--------------------------------------------------------------------------------------- */
RealOpenMM getTau( void ) const;
/**---------------------------------------------------------------------------------------
Get array of fixed parameters indexed by 'FixedParameters' enums
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* getFixedParameters( void ) const;
/**---------------------------------------------------------------------------------------
Print parameters
@param message message
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int printParameters( std::stringstream& message ) const;
/**---------------------------------------------------------------------------------------
Update
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int update( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses, const RealOpenMM currentPE, const int stepType );
/**---------------------------------------------------------------------------------------
Velocity update; Langevin Leapfrog
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param inverseMasses inverse atom masses
@param xVector xVector
--------------------------------------------------------------------------------------- */
void halfKick( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses, RealOpenMM* inverseMasses, RealOpenMM** xVector );
/**---------------------------------------------------------------------------------------
Position update; Langevin Leapfrog
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param inverseMasses inverse atom masses
@param xPrime xPrime
--------------------------------------------------------------------------------------- */
void drift( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses, RealOpenMM* inverseMasses,
RealOpenMM** xPrime );
/**---------------------------------------------------------------------------------------
Write state
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param state 0 if initial state; otherwise nonzero
@param baseFileName base file name
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int writeState( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses,
int state, const std::string& baseFileName ) const;
/**---------------------------------------------------------------------------------------
Find forces OR positions inside subspace (defined as the span of the 'eigenvectors' Q)
Take 'array' as input, 'outArray' as output (may be the same vector).
'projectionMode' determines:
a) if the projection is for force or positions, bit 1 set=force.
b) if the projection is for the sub-space or-complement space, bit 2 set=sub-space.
----------------------------------------------------------------------------------------- */
void subspaceProjection( RealOpenMM** arrayParam,
RealOpenMM** outArrayParam,
int numberOfAtoms,
RealOpenMM* masses,
RealOpenMM* inverseMasses,
int projectionMode);
};
// ---------------------------------------------------------------------------------------
#endif // __ReferenceNMLDynamics_H__
#ifndef OPENMM_REFERENCE_NMLKERNELFACTORY_H_
#define OPENMM_REFERENCE_NMLKERNELFACTORY_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/KernelFactory.h"
namespace OpenMM {
/**
* This KernelFactory creates all kernels for NormalModeLangevin.
*/
class ReferenceNMLKernelFactory : public KernelFactory {
public:
KernelImpl* createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const;
};
} // namespace OpenMM
#endif /*OPENMM_REFERENCE_NMLKERNELFACTORY_H_*/
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferenceIntegrateNMLStepKernel.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/internal/ContextImpl.h"
#include "SimTKReference/ReferenceCCMAAlgorithm.h"
#include <vector>
using namespace std;
namespace OpenMM {
// These methods lifted from ReferenceKernels.cpp, so NormalModeLangevin could be a plugin
// allocateIntArray
// disposeIntArray
// extractPositions
// extractVelocities
// extractForces
// findAnglesForCCMA
static int** allocateIntArray(int length, int width) {
int** array = new int*[length];
for (int i = 0; i < length; ++i)
array[i] = new int[width];
return array;
}
static void disposeIntArray(int** array, int size) {
if (array) {
for (int i = 0; i < size; ++i)
delete[] array[i];
delete[] array;
}
}
static RealOpenMM** extractPositions(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return (RealOpenMM**) data->positions;
}
static RealOpenMM** extractVelocities(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return (RealOpenMM**) data->velocities;
}
static RealOpenMM** extractForces(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return (RealOpenMM**) data->forces;
}
static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
for (int i = 0; i < system.getNumForces(); i++) {
const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
if (force != NULL) {
for (int j = 0; j < force->getNumAngles(); j++) {
int atom1, atom2, atom3;
double angle, k;
force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM)angle));
}
}
}
}
ReferenceIntegrateNMLStepKernel::~ReferenceIntegrateNMLStepKernel() {
if (dynamics)
delete dynamics;
if (constraints)
delete constraints;
if (masses)
delete[] masses;
if (constraintIndices)
disposeIntArray(constraintIndices, numConstraints);
if (constraintDistances)
delete[] constraintDistances;
if (projectionVectors)
delete projectionVectors;
}
void ReferenceIntegrateNMLStepKernel::initialize(const System& system, const NMLIntegrator& integrator) {
int numParticles = system.getNumParticles();
masses = new RealOpenMM[numParticles];
for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
numConstraints = system.getNumConstraints();
constraintIndices = allocateIntArray(numConstraints, 2);
constraintDistances = new RealOpenMM[numConstraints];
for (int i = 0; i < numConstraints; ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
constraintIndices[i][0] = particle1;
constraintIndices[i][1] = particle2;
constraintDistances[i] = static_cast<RealOpenMM>(distance);
}
//SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
}
void ReferenceIntegrateNMLStepKernel::execute(ContextImpl& context, const NMLIntegrator& integrator, const double currentPE, const int stepType) {
double temperature = integrator.getTemperature();
double friction = integrator.getFriction();
double stepSize = integrator.getStepSize();
double* dProjectionVectors = integrator.getProjectionVectors();
unsigned int numProjectionVectors = integrator.getNumProjectionVectors();
bool projVecChanged = integrator.getProjVecChanged();
double minimumLimit = integrator.getMinimumLimit();
double maxEig = integrator.getMaxEigenvalue();
RealOpenMM** posData = extractPositions(context);
RealOpenMM** velData = extractVelocities(context);
RealOpenMM** forceData = extractForces(context);
//projection vectors
if (projectionVectors == 0 || projVecChanged) {
std::cout << "Setting vectors " << projVecChanged << std::endl;
unsigned int arraySz = numProjectionVectors * context.getSystem().getNumParticles() * 3;
if(projectionVectors == 0) {
projectionVectors = new RealOpenMM[arraySz];
}
for(unsigned int i=0; i<arraySz; i++) {
projectionVectors[i] = static_cast<RealOpenMM>(dProjectionVectors[i]);
}
}
if (dynamics == 0 ){ //|| temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
// Recreate the computation objects with the new parameters.
if (dynamics) {
delete dynamics;
delete constraints;
}
RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
dynamics = new ReferenceNMLDynamics(context.getSystem().getNumParticles(),
static_cast<RealOpenMM>(stepSize),
static_cast<RealOpenMM>(tau),
static_cast<RealOpenMM>(temperature),
projectionVectors,
numProjectionVectors,
static_cast<RealOpenMM>(minimumLimit),
static_cast<RealOpenMM>(maxEig));
vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
findAnglesForCCMA(context.getSystem(), angles);
constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
dynamics->setReferenceConstraintAlgorithm(constraints);
prevTemp = temperature;
prevFriction = friction;
prevStepSize = stepSize;
}
//for minimization quadratic fit we need the potential energy
RealOpenMM currentPE2 = static_cast<RealOpenMM>(currentPE);
//
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, currentPE2, stepType );
//update at dynamic step 2
if(stepType == 2){
data.time += stepSize;
data.stepCount++;
}
}
} /* namespace OpenMM */
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns, Pande Group *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include <string.h>
#include <sstream>
#include "SimTKUtilities/SimTKOpenMMCommon.h"
#include "SimTKUtilities/SimTKOpenMMLog.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceNMLDynamics.h"
/**---------------------------------------------------------------------------------------
ReferenceNMLDynamics constructor
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param tau viscosity(?)
@param temperature temperature
--------------------------------------------------------------------------------------- */
ReferenceNMLDynamics::ReferenceNMLDynamics( int numberOfAtoms,
RealOpenMM deltaT, RealOpenMM tau,
RealOpenMM temperature,
RealOpenMM* projectionVectors,
unsigned int numProjectionVectors,
RealOpenMM minimumLimit, RealOpenMM maxEig ) :
ReferenceDynamics( numberOfAtoms, deltaT, temperature ), _tau( tau ),
_projectionVectors(projectionVectors), _numProjectionVectors(numProjectionVectors),
_minimumLimit(minimumLimit), _maxEig(maxEig) {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceNMLDynamics::ReferenceNMLDynamics";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
// insure tau is not zero -- if it is print warning message
if( _tau == zero ){
std::stringstream message;
message << methodName;
message << " input tau value=" << tau << " is invalid -- setting to 1.";
SimTKOpenMMLog::printError( message );
_tau = one;
}
_setFixedParameters( );
allocate2DArrays( numberOfAtoms, 3, Max2DArrays );
allocate1DArrays( numberOfAtoms, Max1DArrays );
}
/**---------------------------------------------------------------------------------------
ReferenceNMLDynamics destructor
--------------------------------------------------------------------------------------- */
ReferenceNMLDynamics::~ReferenceNMLDynamics( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceNMLDynamics::~ReferenceNMLDynamics";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Set fixed parameters
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceNMLDynamics::_setFixedParameters( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceNMLDynamics::_setFixedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM half = 0.5;
static const RealOpenMM eighty = 80.0;
// ---------------------------------------------------------------------------------------
//dt, myGamma, fdt, vdt, ndt, sqrtFCoverM from Langevin Leapfrog
_fixedParameters[DT] = (RealOpenMM) getTimeStep();
_fixedParameters[GAMMA] = (_tau == zero ? eighty : one/_tau);
_fixedParameters[FDT] = ( one - EXP( -half * _fixedParameters[GAMMA] * _fixedParameters[DT] ) ) / _fixedParameters[GAMMA];
_fixedParameters[VDT] = EXP(-half * _fixedParameters[GAMMA] * _fixedParameters[DT]);
_fixedParameters[NDT] = SQRT( ( one - EXP( -_fixedParameters[GAMMA] * _fixedParameters[DT] ) ) / (two * _fixedParameters[GAMMA]) );
_fixedParameters[SQRTFCOVERM] = SQRT( two * (RealOpenMM) BOLTZ * getTemperature() * _fixedParameters[GAMMA] );
return 0; // ReferenceDynamics::DefaultReturn;
};
/**---------------------------------------------------------------------------------------
Get tau
@return tau
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceNMLDynamics::getTau( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceNMLDynamics::getTau";
// ---------------------------------------------------------------------------------------
return _tau;
}
/**---------------------------------------------------------------------------------------
Get array of fixed parameters indexed by 'FixedParameters' enums
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* ReferenceNMLDynamics::getFixedParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceNMLDynamics::getFixedParameters";
// ---------------------------------------------------------------------------------------
return _fixedParameters;
}
/**---------------------------------------------------------------------------------------
Print parameters
@param message message
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceNMLDynamics::printParameters( std::stringstream& message ) const {
// ---------------------------------------------------------------------------------------
//dt, myGamma, fdt, vdt, ndt, sqrtFCoverM from Langevin Leapfrog
static const char* methodName = "\nReferenceNMLDynamics::printParameters";
static const char* parameterNames[MaxFixedParameters] = { "dt", "myGamma", "fdt", "ndt", "sqrtFCoverM" };
// ---------------------------------------------------------------------------------------
// print parameters
ReferenceDynamics::printParameters( message );
message << " tau=" << getTau();
message << " T=" << getTemperature();
int cut = 3;
for( int ii = 0; ii < MaxFixedParameters; ii++ ){
message << " " << parameterNames[ii] << "=" << _fixedParameters[ii];
if( cut++ > 5 ){
cut = 0;
message << std::endl;
}
}
return 0; // ReferenceDynamics::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Drift update
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param inverseMasses inverse atom masses
@param xPrime xPrime
@param oldVelocities previous velocities
@param xVector xVector
@param vVector vVector
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
void ReferenceNMLDynamics::halfKick( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses,
RealOpenMM* inverseMasses, RealOpenMM** xVector ){
//dt, myGamma, fdt, vdt, ndt, sqrtFCoverM from Langevin Leapfrog
static const RealOpenMM dt = (RealOpenMM) getTimeStep(); // in ps
static const RealOpenMM myGamma = (RealOpenMM) (_tau == 0.0 ? 80.0 : 1.0/_tau);
static const RealOpenMM fdt = (RealOpenMM) (( 1.0 - EXP( -0.5 * myGamma * dt ) ) / myGamma);
static const RealOpenMM vdt = (RealOpenMM) EXP(-0.5 * myGamma * dt);
static const RealOpenMM ndt = (RealOpenMM) SQRT( ( 1.0 - EXP( -myGamma * dt ) ) / (2.0 * myGamma) );
static const RealOpenMM sqrtFCoverM = (RealOpenMM) SQRT( 2.0 * (RealOpenMM) BOLTZ * getTemperature() * myGamma );
//project forces into sub space. TODO put this somewhere callable from NMLIntegrator, else projecting twice.
subspaceProjection(forces, forces, numberOfAtoms, masses, inverseMasses, 0);
for( int ii = 0; ii < numberOfAtoms; ii++ ){
//RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
for( int jj = 0; jj < 3; jj++ ){
xVector[ii][jj] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber() * SQRT( masses[ii] );
}
}
subspaceProjection(xVector, xVector, numberOfAtoms, masses, inverseMasses, 0);
for( int ii = 0; ii < numberOfAtoms; ii++ ){
//RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
for( int jj = 0; jj < 3; jj++ ){
xVector[ii][jj] *= inverseMasses[ii] ;
}
}
for( int i = 0; i < numberOfAtoms; i++ ) {
// semi-update velocities
for( int j = 0; j < 3; j++ ) {
velocities[i][j] = velocities[i][j] * vdt
+ forces[i][j] * fdt * inverseMasses[i]
+ xVector[i][j] * sqrtFCoverM * ndt;
//+ SimTKOpenMMUtilities::getNormallyDistributedRandomNumber() * sqrtFCoverM * ndt;
}
}
//project forces into complement space, put in xPrime
subspaceProjection(velocities, velocities, numberOfAtoms, masses, inverseMasses, 1); //in subspace but position/velocity projection
}
void ReferenceNMLDynamics::drift( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses, RealOpenMM* inverseMasses,
RealOpenMM** xPrime ){
static const RealOpenMM deltaT = (RealOpenMM) getTimeStep();
for( int i = 0; i < numberOfAtoms; i++ ){
for( int j = 0; j < 3; j++ ) {
xPrime[i][j] = atomCoordinates[i][j] + deltaT * velocities[i][j];
}
}
}
/**---------------------------------------------------------------------------------------
Update -- driver routine for performing stochastic dynamics update of coordinates
and velocities
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceNMLDynamics::update( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses,
const RealOpenMM currentPE, const int stepType ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceNMLDynamics::update";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static int debug = 0;
// ---------------------------------------------------------------------------------------
// get work arrays
RealOpenMM** xPrime = get2DArrayAtIndex( xPrime2D );
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
RealOpenMM** xVector = get2DArrayAtIndex( X2D );
RealOpenMM** vVector = get2DArrayAtIndex( V2D );
RealOpenMM* inverseMasses = get1DArrayAtIndex( InverseMasses );
// first-time-through initialization
if( getTimeStep() == 0 ){
std::stringstream message;
message << methodName;
int errors = 0;
// invert masses
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( masses[ii] <= zero ){
message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
errors++;
} else {
inverseMasses[ii] = one/masses[ii];
}
}
// set xVector
for( int ii = 0; ii < numberOfAtoms; ii++ ){
//RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
for( int jj = 0; jj < 3; jj++ ){
xVector[ii][jj] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber() * SQRT( masses[ii] );
}
}
subspaceProjection(xVector, xVector, numberOfAtoms, masses, inverseMasses, 0);
for( int ii = 0; ii < numberOfAtoms; ii++ ){
//RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
for( int jj = 0; jj < 3; jj++ ){
xVector[ii][jj] *= inverseMasses[ii] ;
}
}
// exit if errors
if( errors ){
SimTKOpenMMLog::printError( message );
}
}
std::stringstream messagea;
messagea << "Tau " << _tau << "\n";
//messagea << "Eigs " << _projectionVectors[0] << "," << _projectionVectors[1] << "," << _projectionVectors[2] << "," << _projectionVectors[3] << "," <<"\n";
//messagea << "Num vect " << _numProjectionVectors << "\n";
SimTKOpenMMLog::printMessage( messagea );
switch(stepType){
case 1:
halfKick( numberOfAtoms, atomCoordinates,
velocities,
forces, masses, inverseMasses, xVector );
drift( numberOfAtoms, atomCoordinates,
velocities,
forces, masses, inverseMasses,
xPrime );
break;
case 2:
halfKick( numberOfAtoms, atomCoordinates,
velocities,
forces, masses, inverseMasses, xVector );
break;
case 3:{
//save current PE in case quadratic required
lastPE = currentPE;
//project forces into complement space, put in xPrime
subspaceProjection(forces, xPrime, numberOfAtoms, masses, inverseMasses, 2);
//just equal to minus dot product of 'proposed position move' and forces (=-\nabla PE)
lastSlope = 0.0;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
for( int jj = 0; jj < 3; jj++ ){
lastSlope -= xPrime[ii][jj] * forces[ii][jj] * inverseMasses[ii]; //posTemp[k].dot((*myForcesP)[k]);
}
}
std::stringstream message1;
message1 << "Last slope " << lastSlope <<"\n";
SimTKOpenMMLog::printMessage( message1 );
// copy xPrime -> atomCoordinates
for( int ii = 0; ii < numberOfAtoms; ii++ ){
const RealOpenMM factor = inverseMasses[ii] / _maxEig;
xPrime[ii][0] = atomCoordinates[ii][0] + factor * xPrime[ii][0];
xPrime[ii][1] = atomCoordinates[ii][1] + factor * xPrime[ii][1];
xPrime[ii][2] = atomCoordinates[ii][2] + factor * xPrime[ii][2];
}
break;
}
case 4:
//project forces into complement space, put in xPrime
subspaceProjection(forces, xPrime, numberOfAtoms, masses, inverseMasses, 2);
//find slope of PE with /lambda here
RealOpenMM lambda = (RealOpenMM) (1.0 / _maxEig);
RealOpenMM oldLambda = lambda;
std::stringstream message;
message << "Lambda "<< lambda << " oldLambda " << oldLambda <<"\n";
SimTKOpenMMLog::printMessage( message );
//just equal to minus dot product of 'proposed position move' and forces (=-\nabla PE)
//RealOpenMM lambdaSlp = 0.0;
//for( int ii = 0; ii < numberOfAtoms; ii++ ){
// for( int jj = 0; jj < 3; jj++ ){
// lambdaSlp -= xPrime[ii][jj] * forces[ii][jj] * inverseMasses[ii]; //posTemp[k].dot((*myForcesP)[k]);
// }
//}
//solve for minimum for quadratic fit using two PE vales and the slope with /lambda=0
RealOpenMM a, b;
//a = -(((currentPE - lastPE) / lambda - lambdaSlp) / lambda);
//b = lambdaSlp - 2.0 * a * lambda;
a = (((currentPE - lastPE) / lambda - lastSlope) / lambda);
b = lastSlope;
if( a != 0.0){
lambda = -b / (2 * a);
}else{
lambda = (RealOpenMM) (oldLambda / 2.0);
}
//test if lambda negative, if so just use smaller lambda
if(lambda <= 0.0) lambda = (RealOpenMM) (oldLambda / 2.0);
//std::stringstream message;
//message << "Lambda "<< lambda << " oldLambda " << oldLambda << " slope " << lambdaSlp << " lastSlope " << lastSlope<<"\n";
message << "Lambda "<< lambda << " oldLambda " << oldLambda << " lastSlope " << lastSlope<<"\n";
SimTKOpenMMLog::printMessage( message );
// copy xPrime -> atomCoordinates
for( int ii = 0; ii < numberOfAtoms; ii++ ){
const RealOpenMM factor = inverseMasses[ii] * (lambda-oldLambda);
xPrime[ii][0] = atomCoordinates[ii][0] + factor * xPrime[ii][0];
xPrime[ii][1] = atomCoordinates[ii][1] + factor * xPrime[ii][1];
xPrime[ii][2] = atomCoordinates[ii][2] + factor * xPrime[ii][2];
}
break;
}
//New NML minimizer
std::stringstream message;
message << "Array "<<forces[0][0]<<" " <<forces[0][1]<<" " <<forces[0][2]<<" "<<xPrime[1][0]<<" " <<xPrime[1][1]<<" " <<xPrime[1][2]<<" " <<"\n";
RealOpenMM* mydata = &forces[0][0];
message << "Lin Array "<<mydata[0]<<" " <<mydata[1]<<" " <<mydata[2]<<" "<<mydata[3]<<" " <<mydata[4]<<" " <<mydata[5]<<" " <<"\n";
forces[0][0] = 1.0f; forces[0][1] = 1.1f; xPrime[1][0] = 2.0f;
message << "Lin Array "<<mydata[0]<<" " <<mydata[1]<<" " <<mydata[2]<<" "<<mydata[3]<<" " <<mydata[4]<<" " <<mydata[5]<<" " <<"\n";
SimTKOpenMMLog::printMessage( message );
//End NML Minimizer
if( debug ){
int maxPrint = 5;
std::stringstream message;
message << methodName << " Post SD2 atoms=" << numberOfAtoms << "\n";
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
for( int ii = 0; ii < maxPrint; ii++ ){
message << " x[";
SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
message << "] xp[";
SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
message << "] v[";
SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
message << "] ov[";
SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
message << "]\n";
}
SimTKOpenMMLog::printMessage( message );
}
#if 0
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm ){
/*
std::stringstream message;
message << methodName;
message << " calling constrain1\n";
SimTKOpenMMLog::printMessage( message );
*/
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
inverseMasses );
}
#endif
// copy xPrime -> atomCoordinates
for( int ii = 0; ii < numberOfAtoms; ii++ ){
atomCoordinates[ii][0] = xPrime[ii][0];
atomCoordinates[ii][1] = xPrime[ii][1];
atomCoordinates[ii][2] = xPrime[ii][2];
}
incrementTimeStep();
return 0; // ReferenceDynamics::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Write state
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param state 0 if initial state; otherwise nonzero
@param baseFileName base file name
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceNMLDynamics::writeState( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses,
int state, const std::string& baseFileName ) const {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceNMLDynamics::writeState";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const int threeI = 3;
// ---------------------------------------------------------------------------------------
std::stringstream stateFileName;
stateFileName << baseFileName;
stateFileName << "_Step" << getTimeStep();
// stateFileName << "_State" << state;
stateFileName << ".txt";
// ---------------------------------------------------------------------------------------
// open file -- return if unsuccessful
FILE* stateFile = NULL;
#ifdef WIN32
fopen_s( &stateFile, stateFileName.str().c_str(), "w" );
#else
stateFile = fopen( stateFileName.str().c_str(), "w" );
#endif
// ---------------------------------------------------------------------------------------
// diagnostics
if( stateFile != NULL ){
std::stringstream message;
message << methodName;
message << " Opened file=<" << stateFileName.str() << ">.\n";
SimTKOpenMMLog::printMessage( message );
} else {
std::stringstream message;
message << methodName;
message << " could not open file=<" << stateFileName.str() << "> -- abort output.\n";
SimTKOpenMMLog::printMessage( message );
return -1; // ReferenceDynamics::ErrorReturn;
}
// ---------------------------------------------------------------------------------------
StringVector scalarNameI;
IntVector scalarI;
StringVector scalarNameR;
RealOpenMMVector scalarR;
StringVector scalarNameR1;
RealOpenMMPtrVector scalarR1;
StringVector scalarNameR2;
RealOpenMMPtrPtrVector scalarR2;
scalarI.push_back( getNumberOfAtoms() );
scalarNameI.push_back( "Atoms" );
scalarI.push_back( getTimeStep() );
scalarNameI.push_back( "Timestep" );
if( state == 0 || state == -1 ){
scalarR.push_back( getDeltaT() );
scalarNameR.push_back( "delta_t" );
scalarR.push_back( getTemperature() );
scalarNameR.push_back( "T" );
scalarR.push_back( getTau() );
scalarNameR.push_back( "Tau" );
scalarR1.push_back( masses );
scalarNameR1.push_back( "mass" );
scalarR2.push_back( atomCoordinates );
scalarNameR2.push_back( "coord" );
scalarR2.push_back( velocities );
scalarNameR2.push_back( "velocities" );
scalarR2.push_back( forces );
scalarNameR2.push_back( "forces" );
if( state == -1 ){
RealOpenMM** xPrime = get2DArrayAtIndex( xPrime2D );
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
RealOpenMM** xVector = get2DArrayAtIndex( X2D );
RealOpenMM** vVector = get2DArrayAtIndex( V2D );
scalarR2.push_back( xPrime );
scalarNameR2.push_back( "xPrime" );
scalarR2.push_back( oldVelocities);
scalarNameR2.push_back( "vold" );
scalarR2.push_back( xVector );
scalarNameR2.push_back( "xVector" );
scalarR2.push_back( vVector );
scalarNameR2.push_back( "vVector" );
}
} else {
scalarR2.push_back( atomCoordinates );
scalarNameR2.push_back( "coord" );
scalarR2.push_back( velocities );
scalarNameR2.push_back( "velocities" );
}
writeStateToFile( stateFile, scalarNameI, scalarI, scalarNameR, scalarR, getNumberOfAtoms(), scalarNameR1, scalarR1, threeI, scalarNameR2, scalarR2 );
(void) fclose( stateFile );
return 0; // ReferenceDynamics::DefaultReturn;
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Find forces OR positions inside subspace (defined as the span of the 'eigenvectors' Q)
// Take 'array' as input, 'outArray' as output (may be the same vector).
// 'projectionMode' determines:
// a) if the projection is for force or positions, bit 1 set=force.
// b) if the projection is for the sub-space or-complement space, bit 2 set=sub-space.
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void ReferenceNMLDynamics::subspaceProjection( RealOpenMM** arrayParam,
RealOpenMM** outArrayParam,
int numberOfAtoms,
RealOpenMM* masses,
RealOpenMM* inverseMasses,
int projectionMode){
//move to linear arrays for Blas etc.
RealOpenMM* array = &arrayParam[0][0];
RealOpenMM* outArray = &outArrayParam[0][0];
//If 'array' and 'outArray are not the same array
//copy 'array' into outArray
const unsigned int _3N = numberOfAtoms * 3;
if (array != outArray) {
for (unsigned int i=0; i<_3N; i++)
outArray[i] = array[i];
}
//Temporary array for the 'mode' values, denoted 'c'
//the size is equal to the number of eigenvectors
RealOpenMM* tmpC = new RealOpenMM[_numProjectionVectors];
//~~~~We need to calculate M^{1/2}QQ^TM^{-1/2}force or M^{-1/2}QQ^TM^{1/2}positions~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//First we weight the array by multiplying by
//the square-root of the atomic masses for positions
//or the the inverse of the square-root of the atomic masses for forces.
//
//a'=M^{-1/2}*f for forces, OR a'= M^{1/2}*x for positions
//
if(!(projectionMode & 0x01)) { //Forces if zero so weight is sqrt(1/m)
//array is 3N long, for N atoms
for( int i=0; i < numberOfAtoms; i++) { //N loops
RealOpenMM massWt = SQRT( inverseMasses[i] );
for (unsigned int j=0; j<3; j++) //times 3 loops
outArray[i*3 + j] *= massWt;
}
} else { //positions if first bit set, so weight is sqrt(m)
//array is 3N long, for N atoms
for( int i=0; i < numberOfAtoms; i++) { //N loops
RealOpenMM massWt = SQRT( masses[i] );
for (unsigned int j=0; j<3; j++) //times 3 loops
outArray[i*3 + j] *= massWt;
}
}
//Project onto mode space by taking the matrix product of
//the transpose of the eigenvectors Q with the array.
//
//c=Q^T*a', a' from last algorithm step
//
#if defined(HAVE_LAPACK)
//~~~~Blas here~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
char *transA = "T"; // Transpose, LAPACK checks only first character N/V
int m = numberOfAtoms * 3; //eigenvector vector/array length
int n = _numProjectionVectors; //number of eigenvectors
int incxy = 1; //sizes
RealOpenMM alpha = 1.0; RealOpenMM beta = 0.0;
//if single precision use Blas 'sgemv'
#if RealOpenMMType == 1
//
sgemv_ (transA, &m, &n, &alpha, _projectionVectors, &m, outArray, &incxy, &beta, tmpC, &incxy);
//
//Now find projected force/positions a'' by matrix product with Eigenvectors Q
//or the complementary Eigenvectors
//a''=Qc
char *transB = "N"; /* LAPACK checks only first character N/V */
//if projectionMode has bit 2 set then complement space (I-Q^TQ) not (Q^TQ)
if(!(projectionMode & 0x02)) {
alpha = 1.0; beta = 0.0;
} else {
alpha = -1.0; beta = 1.0;
}
//
sgemv_ (transB, &m, &n, &alpha, _projectionVectors, &m, tmpC, &incxy, &beta, outArray, &incxy);
//if double precision use Blas 'dgemv'
#else
//
dgemv_ (transA, &m, &n, &alpha, _projectionVectors), &m, outArray, &incxy, &beta, tmpC, &incxy);
//
//Now find projected force/positions a'' by matrix product with Eigenvectors Q
//a''=Qc
char *transB = "N"; /* LAPACK checks only first character N/V */
//if projectionMode has bit 2 set then complement space (I-Q^TQ) not (Q^TQ)
if(!(projectionMode & 0x02)) {
alpha = 1.0; beta = 0.0;
} else {
alpha = -1.0; beta = 1.0;
}
//
dgemv_ (transB, &m, &n, &alpha, _projectionVectors, &m, tmpC, &incxy, &beta, outArray, &incxy);
#endif
#else
//~~~~NO Blas here~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//If no Blas is available we need to manually find the product c=A*b
//c_i=\sum_{j=1}^n A_{i,j} b_j
//c=Q^T*a', a' from last algorithm step
//Q is a linear array in column major format
//so tmpC_i = \sum_{j=1}^n Q_{j,i} outArray_j
//Q_{j,i}=_projectionVectors[j * numberOfAtoms * 3 + i]
//const int _3N = numberOfAtoms * 3; //eigenvector vector/array length
for(int i=0; i < (int) _numProjectionVectors; i++){ //over all eigenvectors
tmpC[i] = 0.0; //clear
for(int j=0; j< (int) _3N; j++){ //over each element in the vector
//tmpC[i] += _projectionVectors[j * _3N + i] * outArray[j];
tmpC[i] += _projectionVectors[j + i * _3N] * outArray[j];
}
}
//Now find projected force/positions a'' by matrix product with Eigenvectors Q
//a''=Qc
//so outArray_i = \sum_{j=1}^n Q_{i,j} tmpC_i
//if projectionMode has bit 2 set then complement space (I-Q^TQ) not (Q^TQ)
const bool pMode = !(projectionMode & 0x02);
//find product
for(int i=0; i< (int) _3N; i++){ //over each element in the vector
//if sub-space do Q*c
//else do a'-Q(Q^T a') = (I-QQ^T)a'
if(pMode){
outArray[i] = 0.0; //if not complement
for(int j=0; j < (int) _numProjectionVectors; j++){ //over all eigenvectors
//outArray[i] += _projectionVectors[i * _3N + j] * tmpC[j];
outArray[i] += _projectionVectors[i + j * _3N] * tmpC[j];
}
}else{
for(int j=0; j < (int) _numProjectionVectors; j++){ //over all eigenvectors
//outArray[i] -= _projectionVectors[i * _3N + j] * tmpC[j];
outArray[i] -= _projectionVectors[i + j * _3N] * tmpC[j];
}
}
}
//~~~~End of no-blas~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#endif
//Finally we weight the array by multiplying by
//the inverse of the square-root of the atomic masses for positions
//or the the square-root of the atomic masses for forces.
//
//a'''=M^{1/2}*a'' or a'''=M^{-1/2}*a''
//
if(!(projectionMode & 0x01)) { //Forces if zero so weight is sqrt(m)
for( int i=0; i < numberOfAtoms; i++) {
RealOpenMM massUnWt = SQRT( masses[i] );
for (unsigned int j=0; j<3; j++)
outArray[i*3 + j] *= massUnWt;
}
} else { //positions if first bit set, so weight is sqrt(1/m)
for( int i=0; i < numberOfAtoms; i++) {
RealOpenMM massUnWt = SQRT( inverseMasses[i] );
for (unsigned int j=0; j<3; j++)
outArray[i*3 + j] *= massUnWt;
}
}
//Clean up
//delete temporary array
delete [] tmpC;
}
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "ReferenceNMLKernelFactory.h"
#include "ReferenceIntegrateNMLStepKernel.h"
#include <iostream>
using namespace std;
namespace OpenMM {
KernelImpl* ReferenceNMLKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
ReferencePlatform::PlatformData& data = *static_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
if (name == IntegrateNMLStepKernel::Name())
return new ReferenceIntegrateNMLStepKernel(name, platform, data);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
}
} /* namespace OpenMM */
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "NMLIntegrator.h"
#include "IntegrateNMLStepKernel.h"
#include "openmm/Context.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/kernels.h"
#include <string>
using namespace OpenMM;
using std::string;
using std::vector;
NMLIntegrator::NMLIntegrator(double temperature, double frictionCoeff, double stepSize, ProtoMol::EigenvectorInfo* projectionVectorInfo ) {
setTemperature(temperature);
setFriction(frictionCoeff);
setStepSize(stepSize);
setProjectionVectorInfo(projectionVectorInfo);
setConstraintTolerance(1e-4);
}
void NMLIntegrator::initialize(ContextImpl& contextRef) {
context = &contextRef;
kernel = context->getPlatform().createKernel(IntegrateNMLStepKernel::Name(), contextRef);
dynamic_cast<IntegrateNMLStepKernel&>(kernel.getImpl()).initialize(contextRef.getSystem(), *this);
}
vector<string> NMLIntegrator::getKernelNames() {
std::vector<std::string> names;
names.push_back(IntegrateNMLStepKernel::Name());
return names;
}
void NMLIntegrator::step(int steps) {
std::cout << "One loop " << std::endl;
//updates per loop. Need it here? ####Moved within loop. Single v multiple steps error if not!
//context->updateContextState();
//context->calcForces();
//loop
for (int i = 0; i < steps; ++i) {
//updates per loop. Need it here?
context->updateContextState();
context->calcForces();
//do half kick and full drift
dynamic_cast<IntegrateNMLStepKernel&>(kernel.getImpl()).execute(*context, *this, 0.0, 1 ); //stepType 1 is halfKick/drift LangevinLeapfrog
//in case projection vectors changed, clear flag
setProjVecChanged(false);
//updates per loop, ####done in minimizer
//context->updateContextState();
//context->calcForces();
//minimize compliment space
minimize(50);
//do half kick
dynamic_cast<IntegrateNMLStepKernel&>(kernel.getImpl()).execute(*context, *this, 0.0, 2 ); //stepType 1 is halfKick LangevinLeapfrog
}
//update time
context->setTime(context->getTime()+getStepSize() * steps);
}
void NMLIntegrator::minimize(int maxsteps) {
//minimum limit
const double minlim = getMinimumLimit();
//loop
for (int i = 0; i < maxsteps; ++i) {
//updates per loop
context->updateContextState();
context->calcForces();
//get initial PE
const double initialPE = context->calcPotentialEnergy();
//minimize
dynamic_cast<IntegrateNMLStepKernel&>(kernel.getImpl()).execute(*context, *this, initialPE, 3 ); //stepType 3 is simple minimizer
setProjVecChanged(false);
//get PE difference
const double currentPE = context->calcPotentialEnergy();
const double peDifference = initialPE - currentPE;
//std::cout << "Loop " << i << ", Current " << currentPE << ", Old " << initialPE << ", lim " << minlim << std::endl;
//end condition met?
if(peDifference < minlim && peDifference >= 0.0){
//std::cout << "Breaking at " << i << ", diff " << peDifference << std::endl;
break;
}
//Require quadratic solution
if(peDifference < 0.0){
//std::cout << "Quadratic at " << i << std::endl;
//update for new positions
context->updateContextState();
context->calcForces();
//minimize, uses quadratic soultion as 'stepsize' not forced
dynamic_cast<IntegrateNMLStepKernel&>(kernel.getImpl()).execute(*context, *this, currentPE, 4 ); //stepType 4 is quadratic minimizer
setProjVecChanged(false);
std::cout << "Quadratic " << i << ", Current " << context->calcPotentialEnergy() << ", Old " << initialPE << ", lim " << minlim << std::endl;
const double quadraticPE = context->calcPotentialEnergy();
if(quadraticPE > initialPE){
break;
}
}
}
}
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Chris Sweet *
* Contributors: Christopher Bruns *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include <iostream>
#include "OpenMM.h"
#include "ReferenceNMLKernelFactory.h"
using namespace OpenMM;
using namespace std;
extern "C" void OPENMM_EXPORT registerPlatforms() {
cout << "calling NML registerPlatforms()..." << endl;
}
extern "C" void OPENMM_EXPORT registerKernelFactories() {
cout << "Initializing Normal Mode Langevin OpenMM plugin..." << endl;
for (int p = 0; p < Platform::getNumPlatforms(); ++p) {
cout << "Platform " << p << " name = " << Platform::getPlatform(p).getName() << endl;
}
// Only register cuda kernels if cuda platform is found
cout << "NML looking for Cuda plugin..." << endl;
try {
Platform& platform = Platform::getPlatformByName("Cuda");
cout << "NML found Cuda platform..." << endl;
// platform.registerKernelFactory("CudaNMLKernelFactory", new CudaNMLKernelFactory());
} catch (const std::exception& exc) { // non fatal
cout << "NML Cuda platform not found. " << exc.what() << endl;
}
cout << "NML looking for Reference plugin..." << endl;
try {
Platform& platform = Platform::getPlatformByName("Reference");
cout << "NML found Reference platform..." << endl;
ReferenceNMLKernelFactory* factory = new ReferenceNMLKernelFactory();
platform.registerKernelFactory("IntegrateNMLStep", factory);
} catch (const std::exception& exc) { // non fatal
cout << "NML Reference platform not found. " << exc.what() << endl;
}
}
add_executable(TestNormalModeLangevin TestNormalModeLangevin.cpp)
target_link_libraries(TestNormalModeLangevin ${SHARED_TARGET})
# pass OPENMM_PLUGIN_DIR as command line argument
add_test(TestNormalModeLangevin
"${EXECUTABLE_OUTPUT_PATH}/TestNormalModeLangevin"
"${test_plugin_dir}")
# does not work - nice try - maybe in cmake 2.8
set_tests_properties(TestNormalModeLangevin PROPERTIES
ENVIRONMENT "OPENMM_PLUGIN_DIR=${test_plugin_dir}")
#include "OpenMM.h"
#include <vector>
#include <string>
#include <iostream>
#include <cstdlib>
using namespace OpenMM;
using namespace std;
void mysetenv(const char* name, const char* value) {
#ifdef _MSC_VER
char buffer[1000];
sprintf(buffer, "%s=%s", name, value);
putenv(buffer);
#else
setenv(name, value, 1);
#endif
}
void testLoadNMLPlugin()
{
string pluginDir = Platform::getDefaultPluginsDirectory();
cout << "Default plugins directory = " << pluginDir << endl;
Platform::loadPluginsFromDirectory(pluginDir);
// Create a context, just to initialize all plugins
System system;
VerletIntegrator integrator(0.01);
Context context(system, integrator);
// Was NormalModeLangevin plugin loaded?
vector<string> kernelName;
kernelName.push_back("IntegrateNMLStep");
cout << "Searching for kernel IntegrateNMLStep" << endl;
Platform& platform = Platform::findPlatform(kernelName); // throws if no platform with kernel
}
int main(int argc, const char* argv[])
{
try
{
// Set OPENMM_PLUGIN_DIR from first command line argument
if (argc > 1) {
const char* plugin_dir = argv[1];
mysetenv("OPENMM_PLUGIN_DIR", plugin_dir);
cout << plugin_dir << endl;
cout << getenv("OPENMM_PLUGIN_DIR") << endl;
}
testLoadNMLPlugin();
cout << "tests passed" << endl;
return 0;
}
catch (const std::exception& exc)
{
cout << "FAILED: " << exc.what() << endl;
return 1;
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment