Skip to content

GitLab

Commit acf79ffc authored by Michael Sherman's avatar Michael Sherman
Browse files

Fortran NaCl example is getting close now

parent 796b5603
Hide whitespace changes
Inline Side-by-side
Showing with 282 additions and 153 deletions
examples/HelloSodiumChlorideInFortran.f90
View file @ acf79ffc
......@@ -22,9 +22,11 @@
!-------------------------------------------------------------------------------
! We'll define this module as a simplified example of the kinds of data
! structures that may already be in an MD program that is to be converted
! to use OpenMM.
! to use OpenMM. Note that we're using data in Angstrom and kcal units; we'll
! show how to safely convert to and from OpenMM's internal units as we go.
MODULE MyAtomInfo
! Simulation parameters:
! Simulation parameters
! ---------------------
real*8 Temperature, FrictionInPerPs, SolventDielectric, SoluteDielectric
parameter(Temperature = 300) !Kelvins
parameter(FrictionInPerPs = 91) !collisions per picosecond
......@@ -38,10 +40,11 @@ MODULE MyAtomInfo
! Currently energy calculation is not available in the GPU kernels so
! asking for it requires slow Reference Platform computation at
! reporting intervals.
! reporting intervals. If you have a big system you'll want this off.
logical, parameter :: WantEnergy = .true.
! Atom and force field information:
! Atom and force field information
! --------------------------------
type Atom
character*4 pdb
real*8 mass, charge, vdwRadiusInAng, vdwEnergyInKcal
......@@ -51,7 +54,7 @@ MODULE MyAtomInfo
end type
integer, parameter :: NumAtoms = 6
type (Atom) :: atoms(NumAtoms) = (/ &
! pdb mass charge vdwRad vdwEnergy gbRad gbScale initPos pos
! pdb mass charge vdwRad vdwEnergy gbRad gbScale initPos runtime
Atom(' NA ',22.99, 1, 1.8680, 0.00277, 1.992, 0.8, (/ 8, 0, 0/), (/0,0,0/)),&
Atom(' CL ',35.45, -1, 2.4700, 0.1000, 1.735, 0.8, (/-8, 0, 0/), (/0,0,0/)),&
Atom(' NA ',22.99, 1, 1.8680, 0.00277, 1.992, 0.8, (/ 0, 9, 0/), (/0,0,0/)),&
......@@ -64,39 +67,52 @@ END MODULE
!-------------------------------------------------------------------------------
! MAIN PROGRAM
!-------------------------------------------------------------------------------
! This makes use of four subroutines that encapsulate all the OpenMM calls:
! myInitializeOpenMM
! myStepWithOpenMM
! myGetOpenMMState
! myTerminateOpenMM
! and one minimalist PDB file writer that has nothing to do with OpenMM:
! myWritePDBFrame
! All of these subroutines can be found later in this file. For use in a real
! MD code you would need to write your own interface routines along these lines.
! Note that the main program does NOT include the OpenMM module.
PROGRAM HelloSodiumChloride
use OpenMM
use MyAtomInfo
implicit none
!-------------------------------------------------------------------------------
! MODELING AND SIMULATION PARAMETERS
!-------------------------------------------------------------------------------
integer NumReports, NumSilentSteps
parameter(NumReports = (SimulationTimeInPs*1000 / ReportIntervalInFs + 0.5))
parameter(NumSilentSteps = (ReportIntervalInFs / StepSizeInFs + 0.5))
type (OpenMM_Objects) omm
character*10 platformName
real*8 timeInPs, energyInKcal
integer frame
call myInitializeOpenMM(omm, platformName)
print "('REMARK Using OpenMM platform ', A)", platformName
call myGetOpenMMState(omm, timeInPs, energyInKcal)
call myWritePDBFrame(0, timeInPs, energyInKcal)
use MyAtomInfo
do frame = 1, NumReports
call myStepWithOpenMM(omm, NumSilentSteps)
call myGetOpenMMState(omm, timeInPs, energyInKcal)
call myWritePDBFrame(frame, timeInPs, energyInKcal)
end do
! Calculate the number of PDB frames we want to write out and how
! many steps to take on the GPU in between.
integer NumReports, NumSilentSteps
parameter(NumReports = (SimulationTimeInPs*1000 / ReportIntervalInFs + 0.5))
parameter(NumSilentSteps = (ReportIntervalInFs / StepSizeInFs + 0.5))
character*10 platformName; real*8 timeInPs, energyInKcal; integer frame
! Clean up top-level heap allocated objects that we're done with now.
call myTerminateOpenMM(omm)
! This is an opaque handle to a container that holds the OpenMM runtime
! objects. You can use any type for this purpose as long as it is
! big enough to hold a pointer. (If you use a pointer type you'll have
! to declare the subroutine interfaces before calling them.)
integer*8 ommHandle
! Set up OpenMM data structures; returns platform name and handle.
call myInitializeOpenMM(ommHandle, platformName)
! Run the simulation:
! (1) Write the first line of the PDB file and the initial configuration.
! (2) Run silently entirely within OpenMM between reporting intervals.
! (3) Write a PDB frame when the time comes.
print "('REMARK Using OpenMM platform ', A)", platformName
call myGetOpenMMState(ommHandle, timeInPs, energyInKcal)
call myWritePDBFrame(0, timeInPs, energyInKcal)
do frame = 1, NumReports
call myStepWithOpenMM(ommHandle, NumSilentSteps)
call myGetOpenMMState(ommHandle, timeInPs, energyInKcal)
call myWritePDBFrame(frame, timeInPs, energyInKcal)
end do
! Clean up OpenMM data structures.
call myTerminateOpenMM(ommHandle)
END PROGRAM
!-------------------------------------------------------------------------------
......@@ -104,12 +120,10 @@ END PROGRAM
!-------------------------------------------------------------------------------
! Given state data, output a single frame (pdb "model") of the trajectory
SUBROUTINE myWritePDBFrame(frameNum, timeInPs, energyInKcal)
use MyAtomInfo
implicit none
integer frameNum
real*8 timeInPs, energyInKcal
use MyAtomInfo; implicit none
integer frameNum; real*8 timeInPs, energyInKcal
integer n
print "('MODEL',5X,I0)", frameNum
print "('REMARK 250 time=', F0.3, ' picoseconds; Energy=', F0.3, ' kcal/mole')", &
timeInPs, energyInKcal
......@@ -117,58 +131,94 @@ SUBROUTINE myWritePDBFrame(frameNum, timeInPs, energyInKcal)
print "('ATOM ', I5, ' ', A4, ' SLT 1 ', 3F8.3, ' 1.00 0.00')", &
n, atoms(n)%pdb, atoms(n)%posInAng
end do
print "('ENDMDL')"
END SUBROUTINE
!-------------------------------------------------------------------------------
! OpenMM-USING CODE
!-------------------------------------------------------------------------------
! The OpenMM Fortran interface module is included only at this point and below.
! Normally these subroutines would be in a separate compilation module; we're
! including them here for simplicity. We suggest that you write them in C++
! if possible, using extern "C" functions to make then callable from your
! Fortran main program. (See the C++ version of this example program for
! an implementation of very similar routines.) However, these routines are
! reimplemented entirely in Fortran 95 below in case you prefer.
!-------------------------------------------------------------------------------
! ------------------------------------------------------------------------------
! INITIALIZE OpenMM DATA STRUCTURES
! ------------------------------------------------------------------------------
! We take these actions here:
! (1) Load any available OpenMM plugins, e.g. Cuda and Brook.
! (2) Fill the OpenMM::System with the force field parameters we want to
! use and the particular set of atoms to be simulated.
! (3) Create an Integrator and a Context associating the Integrator with
! the System.
! (4) Select the OpenMM platform to be used.
! (5) Allocate a RuntimeObjects container to hang on to the System,
! Integrator, and Context.
! (6) Return an opaque handle to the RuntimeObjects and the name of the
! Platform in use.
!
! Note that this routine must understand the calling MD code's molecule and
! force field data structures so will need to be customized for each MD code.
SUBROUTINE myInitializeOpenMM(omm, platformName)
use OpenMM; use MyAtomInfo
implicit none
type (OpenMM_Objects) omm
character*10 platformName
SUBROUTINE myInitializeOpenMM(ommHandle, platformName)
use OpenMM; use MyAtomInfo; implicit none
integer*8, intent(out) :: ommHandle
character*10, intent(out) :: platformName
! This is the actual type of the opaque handle.
type (OpenMM_RuntimeObjects) omm
! These are the objects we'll create here and store in the RuntimeObjects
! container for later access.
type (OpenMM_System) system
type (OpenMM_Integrator) integrator
type (OpenMM_Context) context
! These are temporary OpenMM objects used and discarded here.
type(OpenMM_Vec3Array) initialPosInNm
type(OpenMM_NonbondedForce) nonbond
type(OpenMM_Force) nonbondAsForce
type(OpenMM_GBSAOBCForce) gbsa
type(OpenMM_Force) gbsaAsForce
type(OpenMM_Force) force ! generic Force
type(OpenMM_LangevinIntegrator) langevin
real*8 posInNm(3)
integer n
character*50 name
type(OpenMM_String) dir
! Create a string, get the name of the default plugins directory,
! and then load all the plugins found there.
call OpenMM_String_create(dir, '')
call OpenMM_Platform_getDefaultPluginsDirectory(dir)
call OpenMM_String_get(dir, name)
print *,'dir="',name,'"'
call OpenMM_Platform_loadPluginsFromDirectory(dir)
call OpenMM_String_destroy(dir)
call OpenMM_System_create(omm%system)
! Create a System and Force objects for it. The System will take
! over ownership of the Forces; don't destroy them yourself.
call OpenMM_System_create(system)
call OpenMM_NonbondedForce_create(nonbond)
call OpenMM_GBSAOBCForce_create(gbsa)
! Convert specific force types to generic OpenMM_Force so that we can
! add them to the OpenMM_System.
call OpenMM_NonbondedForce_asForce(nonbond, nonbondAsForce)
call OpenMM_GBSAOBCForce_asForce(gbsa, gbsaAsForce)
call OpenMM_System_addForce(omm%system, nonbondAsForce)
call OpenMM_System_addForce(omm%system, gbsaAsForce)
call OpenMM_NonbondedForce_asForce(nonbond, force)
call OpenMM_System_addForce(system, force)
call OpenMM_GBSAOBCForce_asForce(gbsa, force)
call OpenMM_System_addForce(system, force)
! Specify dielectrics for GBSA implicit solvation.
call OpenMM_GBSAOBCForce_setSolventDielectric(gbsa, SolventDielectric)
call OpenMM_GBSAOBCForce_setSoluteDielectric(gbsa, SoluteDielectric)
call OpenMM_Vec3Array_create(initialPosInNm, 0)
! Specify the atoms and their properties:
! (1) System needs to know the masses.
! (2) NonbondedForce needs charges,van der Waals properties (in MD units!).
! (3) GBSA needs charge, radius, and scale factor.
! (4) Collect default positions for initializing the simulation later.
call OpenMM_Vec3Array_create(initialPosInNm, NumAtoms)
do n=1,NumAtoms
print *,'atom ',n,atoms(n)
call OpenMM_System_addParticle(omm%system, atoms(n)%mass)
call OpenMM_System_addParticle(system, atoms(n)%mass)
call OpenMM_NonbondedForce_addParticle(nonbond, &
atoms(n)%charge, &
......@@ -181,84 +231,113 @@ SUBROUTINE myInitializeOpenMM(omm, platformName)
atoms(n)%gbsaRadiusInAng * OpenMM_NmPerAngstrom, &
atoms(n)%gbsaScaleFactor)
! Convert this atom's initial position from Angstroms to nm
call OpenMM_Vec3_scale(atoms(n)%initPosInAng, OpenMM_NmPerAngstrom, posInNm)
call OpenMM_Vec3Array_append(initialPosInNm, posInNm)
! Sets initPos(n) = atoms(n)%initPos * nm/Angstrom.
call OpenMM_Vec3Array_setScaled(initialPosInNm, n, &
atoms(n)%initPosInAng, &
OpenMM_NmPerAngstrom)
end do
! Choose an Integrator for advancing time, and a Context connecting the
! System with the Integrator for simulation. Let the Context choose the
! best available Platform. Initialize the configuration from the default
! positions we collected above. Initial velocities will be zero but could
! have been set here.
call OpenMM_LangevinIntegrator_create(langevin, &
Temperature, FrictionInPerPs, &
StepSizeInFs * OpenMM_PsPerFs)
call OpenMM_LangevinIntegrator_asIntegrator(langevin, omm%integrator)
call OpenMM_Context_create(omm%context, omm%system, omm%integrator)
call OpenMM_Context_setPositions(omm%context, initialPosInNm)
call OpenMM_LangevinIntegrator_asIntegrator(langevin, integrator)
call OpenMM_Context_create(context, system, integrator)
call OpenMM_Context_setPositions(context, initialPosInNm)
! Get the platform name to return.
call OpenMM_Context_getPlatformName(context, platformName)
call OpenMM_Context_getPlatformName(omm%context, platformName)
! Put the System, Integrator, and Context in the RuntimeObjects
! container and return an opaque reference to the container.
call OpenMM_RuntimeObjects_create(omm)
call OpenMM_RuntimeObjects_setSystem(omm, system)
call OpenMM_RuntimeObjects_setIntegrator(omm, integrator)
call OpenMM_RuntimeObjects_setContext(omm, context)
ommHandle = transfer(omm, ommHandle)
END SUBROUTINE
!-------------------------------------------------------------------------------
! COPY STATE BACK TO CPU FROM OpenMM
!-------------------------------------------------------------------------------
SUBROUTINE myGetOpenMMState(omm, timeInPs, energyInKcal)
use OpenMM; use MyAtomInfo
implicit none
type (OpenMM_Objects) omm
real*8 timeInPs, energyInKcal
type (OpenMM_State) state
type (OpenMM_Vec3Array) posArray
integer infoMask
integer n
integer npos, i, j
real*8 energy
real*8 posInAng(3), posInNm(3)
infoMask = OpenMM_State_Positions
if (WantEnergy) then
infoMask = infoMask + OpenMM_State_Velocities ! for KE (cheap)
infoMask = infoMask + OpenMM_State_Energy ! for PE (very expensive)
end if
! Forces are also available (and cheap).
! Don't forget to destroy this State when you're done with it.
call OpenMM_Context_createState(omm%context, infoMask, state)
timeInPs = OpenMM_State_getTime(state) ! OpenMM time is in ps already.
SUBROUTINE myGetOpenMMState(ommHandle, timeInPs, energyInKcal)
use OpenMM; use MyAtomInfo; implicit none
integer*8, intent(in) :: ommHandle
real*8, intent(out) :: timeInPs, energyInKcal
type (OpenMM_State) state
type (OpenMM_Vec3Array) posArray
integer infoMask, n
real*8 posInNm(3)
type (OpenMM_RuntimeObjects) omm
type (OpenMM_Context) context
! Positions are maintained as a Vec3Array inside the State. This will give
! us access, but don't destroy it yourself -- it will go away with the State.
call OpenMM_State_getPositions(state, posArray)
do n = 1, NumAtoms
call OpenMM_Vec3Array_get(posArray, n, posInNm)
call OpenMM_Vec3_scale(posInNm, OpenMM_AngstromsPerNm, atoms(n)%posInAng)
end do
omm = transfer(ommHandle, omm)
call OpenMM_RuntimeObjects_getContext(omm, context)
infoMask = OpenMM_State_Positions
if (WantEnergy) then
infoMask = infoMask + OpenMM_State_Velocities ! for KE (cheap)
infoMask = infoMask + OpenMM_State_Energy ! for PE (very expensive)
end if
! Forces are also available (and cheap).
! Don't forget to destroy this State when you're done with it.
call OpenMM_Context_createState(context, infoMask, state)
timeInPs = OpenMM_State_getTime(state) ! OpenMM time is in ps already.
! Positions are maintained as a Vec3Array inside the State. This will give
! us access, but don't destroy it yourself -- it will go away with the State.
call OpenMM_State_getPositions(state, posArray)
do n = 1, NumAtoms
call OpenMM_Vec3Array_get(posArray, n, posInNm)
call OpenMM_Vec3_scale(posInNm, OpenMM_AngstromsPerNm, atoms(n)%posInAng)
end do
energyInKcal = 0
if (WantEnergy) then
energyInKcal = ( OpenMM_State_getPotentialEnergy(state) &
+ OpenMM_State_getKineticEnergy(state)) &
* OpenMM_KcalPerKJ
end if
! Clean up the State memory
call OpenMM_State_destroy(state)
END SUBROUTINE
energyInKcal = 0
if (WantEnergy) then
energyInKcal = ( OpenMM_State_getPotentialEnergy(state) &
+ OpenMM_State_getKineticEnergy(state)) &
* OpenMM_KcalPerKJ
end if
! Clean up the State memory
call OpenMM_State_destroy(state)
END SUBROUTINE
!-------------------------------------------------------------------------------
! TAKE MULTIPLE STEPS USING OpenMM
!-------------------------------------------------------------------------------
SUBROUTINE myStepWithOpenMM(ommHandle, numSteps)
use OpenMM; implicit none
integer*8, intent(in) :: ommHandle
integer, intent(in) :: numSteps
SUBROUTINE myStepWithOpenMM(omm, numSteps)
use OpenMM
implicit none
type (OpenMM_Objects) omm
integer numSteps
type (OpenMM_RuntimeObjects) omm
type (OpenMM_Integrator) integrator
omm = transfer(ommHandle, omm)
call OpenMM_RuntimeObjects_getIntegrator(omm, integrator)
call OpenMM_Integrator_step(omm%integrator, numSteps)
call OpenMM_Integrator_step(integrator, numSteps)
END SUBROUTINE
!-------------------------------------------------------------------------------
! DEALLOCATE ALL OpenMM OBJECTS
!-------------------------------------------------------------------------------
SUBROUTINE myTerminateOpenMM(ommHandle)
use OpenMM; implicit none
integer*8, intent(in) :: ommHandle
type (OpenMM_RuntimeObjects) omm
omm = transfer(ommHandle, omm)
SUBROUTINE myTerminateOpenMM(omm)
use OpenMM
implicit none
type (OpenMM_Objects) omm
call OpenMM_Objects_destroy(omm)
call OpenMM_RuntimeObjects_destroy(omm)
END SUBROUTINE
examples/wrappers/openmm.f90
View file @ acf79ffc
......@@ -6,12 +6,28 @@
module OpenMM_Types
implicit none
! The System, Integrator, and Context must persist between calls.
! They can be conveniently grouped in a RuntimeObjects structure.
type OpenMM_System
character, pointer :: handle => NULL()
end type
! This is the generic Integrator type; it represents one of
! the concrete integrators like Verlet or Langevin.
type OpenMM_Integrator
character, pointer :: handle => NULL()
end type
type OpenMM_Context
character, pointer :: handle => NULL()
end type
! This data structure can be used to hold the set of OpenMM objects
! that must persist from call to call while running a simulation.
! It contains an OpenMM_System, _Integrator, and _Context.
type OpenMM_RuntimeObjects
character, pointer :: handle => NULL()
end type
type OpenMM_State
character, pointer :: handle => NULL()
end type
......@@ -21,8 +37,7 @@ module OpenMM_Types
type OpenMM_String
character, pointer :: handle => NULL()
end type
! This is the generic Force type. Each concrete Force type should
! be able to convert itself to this type.
! This is the generic Force type.
type OpenMM_Force
character, pointer :: handle => NULL()
end type
......@@ -32,11 +47,6 @@ module OpenMM_Types
type OpenMM_GBSAOBCForce
character, pointer :: handle => NULL()
end type
! This is the generic Integrator type. Each concrete Integrator type should
! be able to convert itself to this type.
type OpenMM_Integrator
character, pointer :: handle => NULL()
end type
type OpenMM_VerletIntegrator
character, pointer :: handle => NULL()
end type
......@@ -69,14 +79,6 @@ module OpenMM_Types
parameter(OpenMM_DegreesPerRadian=180.0/3.1415926535897932385)
parameter(OpenMM_SigmaPerVdwRadius=1.78179743628068)
! This data structure can be used to hold the set of OpenMM objects
! that must persist from call to call while running a simulation.
type OpenMM_Objects
type (OpenMM_System) system
type (OpenMM_Integrator) integrator
type (OpenMM_Context) context
end type
end module OpenMM_Types
module OpenMM
......@@ -113,7 +115,32 @@ module OpenMM
use OpenMM_Types
type (OpenMM_Vec3Array) array
integer*4 i
real*8 v3(3)
real*8, intent(out) :: v3(3)
end
subroutine OpenMM_Vec3Array_getScaled(array, i, s, v3)
use OpenMM_Types
type (OpenMM_Vec3Array) array
integer*4 i
real*8 s
real*8, intent(out) :: v3(3)
end
subroutine OpenMM_Vec3Array_set(array, i, v3)
use OpenMM_Types
type (OpenMM_Vec3Array) array
integer*4 i
real*8, intent(in) :: v3(3)
end
subroutine OpenMM_Vec3Array_setScaled(array, i, v3, s)
use OpenMM_Types
type (OpenMM_Vec3Array) array
integer*4 i
real*8, intent(in) :: v3(3)
real*8 s
end
subroutine OpenMM_Vec3_scale(v3in, s, v3out)
real*8, intent(in) :: v3in(3)
real*8 s
real*8, intent(out) :: v3out(3)
end
! OpenMM_String is an interface to std::string, with some
......@@ -364,25 +391,48 @@ module OpenMM
type (OpenMM_State) state
type (OpenMM_Vec3Array) velocities
end
subroutine OpenMM_RuntimeObjects_create(omm)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
end
subroutine OpenMM_RuntimeObjects_clear(omm)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
end
subroutine OpenMM_RuntimeObjects_destroy(omm)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
end
subroutine OpenMM_RuntimeObjects_setSystem(omm,sys)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
type (OpenMM_System) sys
end
subroutine OpenMM_RuntimeObjects_setIntegrator(omm,integ)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
type (OpenMM_Integrator) integ
end
subroutine OpenMM_RuntimeObjects_setContext(omm,context)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
type (OpenMM_Context) context
end
subroutine OpenMM_RuntimeObjects_getSystem(omm,sys)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
type (OpenMM_System) sys
end
subroutine OpenMM_RuntimeObjects_getIntegrator(omm,integ)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
type (OpenMM_Integrator) integ
end
subroutine OpenMM_RuntimeObjects_getContext(omm,context)
use OpenMM_Types
type (OpenMM_RuntimeObjects) omm
type (OpenMM_Context) context
end
end interface
CONTAINS
subroutine OpenMM_Objects_destroy(omm)
use OpenMM_Types
type (OpenMM_Objects) omm
call OpenMM_Context_destroy(omm%context)
call OpenMM_Integrator_destroy(omm%integrator)
call OpenMM_System_destroy(omm%system)
end subroutine
subroutine OpenMM_Vec3_scale(in, s, out)
real*8,intent(in) :: in(3)
real*8,intent(out) :: out(3)
real*8 s
out(1) = in(1) * s
out(2) = in(2) * s
out(3) = in(3) * s
end subroutine
end module OpenMM
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