! ----------------------------------------------------------------------------- ! 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 !------------------------------------------------------------------------------- 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) do frame = 1, NumReports call myStepWithOpenMM(omm, NumSilentSteps) call myGetOpenMMState(omm, timeInPs, energyInKcal) call myWritePDBFrame(frame, timeInPs, energyInKcal) end do ! Clean up top-level heap allocated objects that we're done with now. call myTerminateOpenMM(omm) END PROGRAM !------------------------------------------------------------------------------- ! PDB FILE WRITER !------------------------------------------------------------------------------- ! 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 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 print "('ENDMDL')" END SUBROUTINE !------------------------------------------------------------------------------- ! OpenMM-USING CODE !------------------------------------------------------------------------------- SUBROUTINE myInitializeOpenMM(omm, platformName) use OpenMM; use MyAtomInfo implicit none type (OpenMM_Objects) omm character*10 platformName ! 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_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 !------------------------------------------------------------------------------- ! 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. ! 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 SUBROUTINE myStepWithOpenMM(omm, numSteps) use OpenMM implicit none type (OpenMM_Objects) omm integer numSteps call OpenMM_Integrator_step(omm%integrator, numSteps) END SUBROUTINE SUBROUTINE myTerminateOpenMM(omm) use OpenMM implicit none type (OpenMM_Objects) omm call OpenMM_Objects_destroy(omm) END SUBROUTINE