Skip to content

GitLab

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

update fortran NaCl example

parent 9187b702
Show whitespace changes
Inline Side-by-side
Showing with 422 additions and 248 deletions
examples/HelloSodiumChlorideInFortran.f90
View file @ fcd7b0c7
! OpenMM HelloSodiumChloride example in Fortran 95
program HelloSodiumChloride
! -----------------------------------------------------------------------------
! OpenMM(tm) HelloSodiumChloride example in Fortran 95 (June 2009)
! ------------------------------------------------------------------------------
! This is a complete, self-contained "hello world" example demonstrating
! GPU-accelerated constant temperature simulation of a very simple system with
! just nonbonded forces, consisting of several sodium (Na+) and chloride (Cl-)
! ions in implicit solvent. A multi-frame PDB file is written to stdout which
! can be read by VMD or other visualization tool to produce an animation of the
! resulting trajectory.
!
! Pay particular attention to the handling of units in this example. Incorrect
! handling of units is a very common error; this example shows how you can
! continue to work with Amber-style units like Angstroms, kCals, and van der
! Waals radii while correctly communicating with OpenMM in nm, kJ, and sigma.
!
! This example is written entirely in Fortran 95, using a Fortran interface
! module which is NOT official parts of the OpenMM distribution.
! ------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! ATOM, FORCE FIELD, AND SIMULATION PARAMETERS
!-------------------------------------------------------------------------------
! 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.
MODULE MyAtomInfo
! Simulation parameters:
real*8 Temperature, FrictionInPerPs, SolventDielectric, SoluteDielectric
parameter(Temperature = 300) !Kelvins
parameter(FrictionInPerPs = 91) !collisions per picosecond
parameter(SolventDielectric = 80) !typical for water
parameter(SoluteDielectric = 2) !typical for protein
real*8 StepSizeInFs, ReportIntervalInFs, SimulationTimeInPs
parameter(StepSizeInFs = 2) !integration step size (fs)
parameter(ReportIntervalInFs = 50) !how often for PDB frame (fs)
parameter(SimulationTimeInPs = 100) !total simulation time (ps)
! Currently energy calculation is not available in the GPU kernels so
! asking for it requires slow Reference Platform computation at
! reporting intervals.
logical, parameter :: WantEnergy = .true.
! Atom and force field information:
type Atom
character*4 pdb
real*8 mass, charge, vdwRadiusInAng, vdwEnergyInKcal
real*8 gbsaRadiusInAng, gbsaScaleFactor
real*8 initPosInAng(3)
real*8 posInAng(3) ! leave room for runtime state info
end type
integer, parameter :: NumAtoms = 6
type (Atom) :: atoms(NumAtoms) = (/ &
! pdb mass charge vdwRad vdwEnergy gbRad gbScale initPos pos
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/)),&
Atom(' CL ',35.45, -1, 2.4700, 0.1000, 1.735, 0.8, (/ 0,-9, 0/), (/0,0,0/)),&
Atom(' NA ',22.99, 1, 1.8680, 0.00277, 1.992, 0.8, (/ 0, 0,-10/), (/0,0,0/)),&
Atom(' CL ',35.45, -1, 2.4700, 0.1000, 1.735, 0.8, (/ 0, 0, 10/), (/0,0,0/)) &
/)
END MODULE
!-------------------------------------------------------------------------------
! MAIN PROGRAM
!-------------------------------------------------------------------------------
PROGRAM HelloSodiumChloride
use OpenMM
use MyAtomInfo
implicit none
!-------------------------------------------------------------------------------
! MODELING AND SIMULATION PARAMETERS
!-------------------------------------------------------------------------------
real*8 StepSizeInFs, ReportIntervalInFs, SimulationTimeInPs
parameter(StepSizeInFs = 2)
parameter(ReportIntervalInFs = 50)
parameter(SimulationTimeInPs = 100)
integer NumSilentSteps
integer NumReports, NumSilentSteps
parameter(NumReports = (SimulationTimeInPs*1000 / ReportIntervalInFs + 0.5))
parameter(NumSilentSteps = (ReportIntervalInFs / StepSizeInFs + 0.5))
type Atom
character*4 pdb
real*8 mass
real*8 charge
real*8 vdwRadiusInAng
real*8 vdwEnergyInKcal
real*8 initPosInAng(3)
end type
type (OpenMM_Objects) omm
character*10 platformName
real*8 timeInPs, energyInKcal
integer frame
integer NAtom
parameter(NAtom=6)
type(Atom) atoms(NAtom)
character*50 name
type(OpenMM_String) dir
call myInitializeOpenMM(omm, platformName)
atoms(1) = Atom(' NA ', 22.99, 1, 1.8680, 0.00277, (/ 8, 0, 0 /))
atoms(2) = Atom(' CL ', 35.45, -1, 2.4700, 0.1000, (/ -8, 0, 0 /))
atoms(3) = Atom(' NA ', 22.99, 1, 1.8680, 0.00277, (/ 0, 9, 0 /))
atoms(4) = Atom(' CL ', 35.45, -1, 2.4700, 0.1000, (/ 0,-9, 0 /))
atoms(5) = Atom(' NA ', 22.99, 1, 1.8680, 0.00277, (/ 0, 0,-10 /))
atoms(6) = Atom(' CL ', 35.45, -1, 2.4700, 0.1000, (/ 0, 0, 10 /))
print "('REMARK Using OpenMM platform ', A)", platformName
call myGetOpenMMState(omm, timeInPs, energyInKcal)
call myWritePDBFrame(0, timeInPs, energyInKcal)
call OpenMM_String_create(dir, '')
call OpenMM_Platform_getDefaultPluginsDirectory(dir)
call OpenMM_String_get(dir, name)
print *,'dir="',name,'"'
call OpenMM_Platform_loadPluginsFromDirectory(dir)
do frame = 1, NumReports
call myStepWithOpenMM(omm, NumSilentSteps)
call myGetOpenMMState(omm, timeInPs, energyInKcal)
call myWritePDBFrame(frame, timeInPs, energyInKcal)
end do
call simulateNaCl(atoms)
! Clean up top-level heap allocated objects that we're done with now.
call myTerminateOpenMM(omm)
! Allow subroutines to inherit type definitions from program level
CONTAINS
END PROGRAM
!-------------------------------------------------------------------------------
! NaCl SIMULATION
! PDB FILE WRITER
!-------------------------------------------------------------------------------
subroutine simulateNaCl(atoms)
implicit none
type(Atom) atoms(NAtom)
type(OpenMM_System) system
type(OpenMM_NonbondedForce) nonbond
type(OpenMM_Force) nonbondAsForce
type(OpenMM_Context) context
type(OpenMM_VerletIntegrator) verlet
type(OpenMM_Integrator) integrator
type(OpenMM_Vec3Array) initialPositionsInNm
character*10 platformName
real*8 posInNm(3),a(3),b(3),c(3), cutoff
! Periodic box size and cutoff in nm
parameter(a=(/5,0,0/), b=(/0,5,0/), c=(/0,0,5/), cutoff=2)
integer i
call OpenMM_System_create(system)
call OpenMM_NonbondedForce_create(nonbond)
call OpenMM_NonbondedForce_asForce(nonbond, nonbondAsForce)
call OpenMM_System_addForce(system, nonbondAsForce)
call OpenMM_NonbondedForce_setNonbondedMethod(nonbond, &
OpenMM_NonbondedForce_CutoffPeriodic)
! cutoff distance here is
call OpenMM_NonbondedForce_setCutoffDistance(nonbond, cutoff)
call OpenMM_NonbondedForce_setPeriodicBoxVectors(nonbond, a, b, c)
call OpenMM_Vec3Array_create(initialPositionsInNm, 0)
do i=1,NAtom
call OpenMM_System_addParticle(system, atoms(i)%mass)
call OpenMM_NonbondedForce_addParticle(nonbond, &
atoms(i)%charge, &
atoms(i)%vdwRadiusInAng * OpenMM_NmPerAngstrom &
* OpenMM_SigmaPerVdwRadius, &
atoms(i)%vdwEnergyInKcal * OpenMM_KJPerKcal)
! Convert this atom's initial position from Angstroms to nm
call OpenMM_Vec3_scale(atoms(i)%initPosInAng, OpenMM_NmPerAngstrom, posInNm)
call OpenMM_Vec3Array_append(initialPositionsInNm, posInNm)
end do
! 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
integer n
call OpenMM_VerletIntegrator_create(verlet, StepSizeInFs * OpenMM_PsPerFs)
call OpenMM_VerletIntegrator_asIntegrator(verlet, integrator)
call OpenMM_Context_create(context, system, integrator)
call OpenMM_Context_setPositions(context, initialPositionsInNm)
print "('MODEL',5X,I0)", frameNum
print "('REMARK 250 time=', F0.3, ' picoseconds; Energy=', F0.3, ' kcal/mole')", &
timeInPs, energyInKcal
do n = 1,NumAtoms
print "('ATOM ', I5, ' ', A4, ' SLT 1 ', 3F8.3, ' 1.00 0.00')", &
n, atoms(n)%pdb, atoms(n)%posInAng
end do
call OpenMM_Context_getPlatformName(context, platformName)
print "('ENDMDL')"
END SUBROUTINE
print "('REMARK Using OpenMM platform ', A)", platformName
call writePDB(atoms, context)
!-------------------------------------------------------------------------------
! OpenMM-USING CODE
!-------------------------------------------------------------------------------
do
call OpenMM_Integrator_step(integrator, NumSilentSteps)
call writePDB(atoms, context)
if (OpenMM_Context_getTime(context) >= SimulationTimeInPs) exit
end do
SUBROUTINE myInitializeOpenMM(omm, platformName)
use OpenMM; use MyAtomInfo
implicit none
type (OpenMM_Objects) omm
character*10 platformName
! Clean up top-level heap allocated objects that we're done with now.
! 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_LangevinIntegrator) langevin
real*8 posInNm(3)
integer n
character*50 name
type(OpenMM_String) dir
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_System_create(omm%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)
! Specify dielectrics for GBSA implicit solvation.
call OpenMM_GBSAOBCForce_setSolventDielectric(gbsa, SolventDielectric)
call OpenMM_GBSAOBCForce_setSoluteDielectric(gbsa, SoluteDielectric)
call OpenMM_Vec3Array_create(initialPosInNm, 0)
do n=1,NumAtoms
print *,'atom ',n,atoms(n)
call OpenMM_System_addParticle(omm%system, atoms(n)%mass)
call OpenMM_Vec3Array_destroy(initialPositionsInNm)
call OpenMM_Context_destroy(context)
call OpenMM_Integrator_destroy(integrator)
call OpenMM_NonbondedForce_addParticle(nonbond, &
atoms(n)%charge, &
atoms(n)%vdwRadiusInAng * OpenMM_NmPerAngstrom &
* OpenMM_SigmaPerVdwRadius, &
atoms(n)%vdwEnergyInKcal * OpenMM_KJPerKcal)
call OpenMM_GBSAOBCForce_addParticle(gbsa, &
atoms(n)%charge, &
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)
end do
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_Context_getPlatformName(omm%context, platformName)
END SUBROUTINE
end subroutine
!-------------------------------------------------------------------------------
! PDB FILE WRITER
! COPY STATE BACK TO CPU FROM OpenMM
!-------------------------------------------------------------------------------
subroutine writePDB(atoms, context)
SUBROUTINE myGetOpenMMState(omm, timeInPs, energyInKcal)
use OpenMM; use MyAtomInfo
implicit none
type(Atom) atoms(NAtom)
type(OpenMM_Context) context
integer, save :: modelFrameNumber = 0
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)
! Caution: at the moment asking for energy requires use of slow Reference
! platform calculation.
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, &
OpenMM_State_Positions + OpenMM_State_Velocities + OpenMM_State_Energy, &
state)
call OpenMM_Context_createState(omm%context, infoMask, state)
timeInPs = OpenMM_State_getTime(state) ! OpenMM time is in ps already.
energy = OpenMM_State_getPotentialEnergy(state) &
+ OpenMM_State_getKineticEnergy(state)
! 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)
npos = OpenMM_Vec3Array_size(posArray)
modelFrameNumber = modelFrameNumber + 1
print "('MODEL',5X,I0)", modelFrameNumber
print "('REMARK 250 time=', F0.3, ' picoseconds; Energy=', F0.3, ' kilojoules/mole')", &
OpenMM_State_getTime(state), energy
do i = 1,npos
call OpenMM_Vec3Array_get(posArray, i, posInNm)
call OpenMM_Vec3_scale(posInNm, OpenMM_AngstromsPerNm, posInAng)
print "('ATOM ', I5, ' ', A4, ' SLT 1 ', 3F8.3, ' 1.00 0.00 ')", &
i, atoms(i)%pdb, posInAng
do n = 1, NumAtoms
call OpenMM_Vec3Array_get(posArray, n, posInNm)
call OpenMM_Vec3_scale(posInNm, OpenMM_AngstromsPerNm, atoms(n)%posInAng)
end do
print "('ENDMDL')"
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
END SUBROUTINE
SUBROUTINE myStepWithOpenMM(omm, numSteps)
use OpenMM
implicit none
type (OpenMM_Objects) omm
integer numSteps
call OpenMM_Integrator_step(omm%integrator, numSteps)
END SUBROUTINE
END PROGRAM
SUBROUTINE myTerminateOpenMM(omm)
use OpenMM
implicit none
type (OpenMM_Objects) omm
call OpenMM_Objects_destroy(omm)
END SUBROUTINE
examples/wrappers/openmm.f90
View file @ fcd7b0c7
......@@ -29,6 +29,9 @@ module OpenMM_Types
type OpenMM_NonbondedForce
character, pointer :: handle => NULL()
end type
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
......@@ -42,16 +45,16 @@ module OpenMM_Types
end type
! OpenMM::State enumerations
integer OpenMM_State_Positions, OpenMM_State_Velocities
integer OpenMM_State_Forces, OpenMM_State_Energy
integer OpenMM_State_Parameters
integer*4 OpenMM_State_Positions, OpenMM_State_Velocities
integer*4 OpenMM_State_Forces, OpenMM_State_Energy
integer*4 OpenMM_State_Parameters
parameter(OpenMM_State_Positions=1, OpenMM_State_Velocities=2)
parameter(OpenMM_State_Forces=4, OpenMM_State_Energy=8)
parameter(OpenMM_State_Parameters=16)
!OpenMM::NonbondedForce enumerations
integer OpenMM_NonbondedForce_NoCutoff, OpenMM_NonbondedForce_CutoffNonPeriodic
integer OpenMM_NonbondedForce_CutoffPeriodic, OpenMM_NonbondedForce_Ewald
integer*4 OpenMM_NonbondedForce_NoCutoff, OpenMM_NonbondedForce_CutoffNonPeriodic
integer*4 OpenMM_NonbondedForce_CutoffPeriodic, OpenMM_NonbondedForce_Ewald
parameter(OpenMM_NonbondedForce_NoCutoff=0, OpenMM_NonbondedForce_CutoffNonPeriodic=1)
parameter(OpenMM_NonbondedForce_CutoffPeriodic=2, OpenMM_NonbondedForce_Ewald=3)
......@@ -66,6 +69,14 @@ 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
......@@ -77,17 +88,17 @@ module OpenMM
subroutine OpenMM_Vec3Array_create(array, n)
use OpenMM_Types
type (OpenMM_Vec3Array) array
integer n
integer*4 n
end
function OpenMM_Vec3Array_size(array)
use OpenMM_Types
integer OpenMM_Vec3Array_size
integer*4 OpenMM_Vec3Array_size
type (OpenMM_Vec3Array) array
end
subroutine OpenMM_Vec3Array_resize(array, n)
use OpenMM_Types
type (OpenMM_Vec3Array) array
integer n
integer*4 n
end
subroutine OpenMM_Vec3Array_destroy(array)
use OpenMM_Types
......@@ -119,7 +130,7 @@ module OpenMM
end
function OpenMM_String_length(string)
use OpenMM_Types
integer OpenMM_String_length
integer*4 OpenMM_String_length
type (OpenMM_String) string
end
subroutine OpenMM_String_get(string, fstring)
......@@ -181,7 +192,7 @@ module OpenMM
subroutine OpenMM_NonbondedForce_setNonbondedMethod(nonbond, method)
use OpenMM_Types
type (OpenMM_NonbondedForce) nonbond
integer method
integer*4 method
end
subroutine OpenMM_NonbondedForce_setCutoffDistance(nonbond, distanceInNm)
use OpenMM_Types
......@@ -193,12 +204,44 @@ module OpenMM
type (OpenMM_NonbondedForce) nonbond
real*8 a(3), b(3), c(3)
end
subroutine OpenMM_NonbondedForce_addParticle(nonbond, charge, sigmaInNm, vdwEnergyInKJ)
subroutine OpenMM_NonbondedForce_addParticle &
(nonbond, charge, sigmaInNm, vdwEnergyInKJ)
use OpenMM_Types
type (OpenMM_NonbondedForce) nonbond
real*8 charge, sigmaInNm, vdwEnergyInKJ
end
! OpenMM::GBSAOBCForce
subroutine OpenMM_GBSAOBCForce_create(gbsa)
use OpenMM_Types
type (OpenMM_GBSAOBCForce) gbsa
end
subroutine OpenMM_GBSAOBCForce_destroy(gbsa)
use OpenMM_Types
type (OpenMM_GBSAOBCForce) gbsa
end
subroutine OpenMM_GBSAOBCForce_asForce(gbsa, force)
use OpenMM_Types
type (OpenMM_GBSAOBCForce) gbsa
type (OpenMM_Force) force
end
subroutine OpenMM_GBSAOBCForce_setSolventDielectric(gbsa, d)
use OpenMM_Types
type (OpenMM_GBSAOBCForce) gbsa
real*8 d
end
subroutine OpenMM_GBSAOBCForce_setSoluteDielectric(gbsa, d)
use OpenMM_Types
type (OpenMM_GBSAOBCForce) gbsa
real*8 d
end
subroutine OpenMM_GBSAOBCForce_addParticle &
(gbsa, charge, radiusInNm, scalingFactor)
use OpenMM_Types
type (OpenMM_GBSAOBCForce) gbsa
real*8 charge, radiusInNm, scalingFactor
end
! OpenMM::Integrator
subroutine OpenMM_Integrator_step(integrator, numSteps)
use OpenMM_Types
......@@ -231,6 +274,28 @@ module OpenMM
integer*4 numSteps
end
! OpenMM::LangevinIntegrator
subroutine OpenMM_LangevinIntegrator_create &
(langevin, temperature, frictionInPerPs, stepSzInPs)
use OpenMM_Types
type (OpenMM_LangevinIntegrator) langevin
real*8 temperature, frictionInPerPs, stepSzInPs
end
subroutine OpenMM_LangevinIntegrator_destroy(langevin)
use OpenMM_Types
type (OpenMM_LangevinIntegrator) langevin
end
subroutine OpenMM_LangevinIntegrator_asIntegrator(langevin, integ)
use OpenMM_Types
type (OpenMM_LangevinIntegrator) langevin
type (OpenMM_Integrator) integ
end
subroutine OpenMM_LangevinIntegrator_step(langevin, numSteps)
use OpenMM_Types
type (OpenMM_LangevinIntegrator) langevin
integer*4 numSteps
end
! OpenMM::Context
subroutine OpenMM_Context_create(context, system, integrator)
use OpenMM_Types
......@@ -299,10 +364,19 @@ module OpenMM
type (OpenMM_State) state
type (OpenMM_Vec3Array) velocities
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)
......
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