Simulation of molecular dynamics in Python

Question

I am searching for a python package that I can use to simulate molecular dynamics in non-equilibrium situations. I need a setup that can handle a fairly large number of molecules in a primarily kinetic theory manner, and that can handle having solid surfaces present. With regards to the surfaces, I would need to be able to create arbitrary shapes and monitor pressure and other variables resulting from the molecular action. Alternatively, I could add the surface parts myself if I had molecules that could handle it.

Does anyone know of any packages that might be suitable?

Answer 1

Have you considered SimPy? SimPy is a rather generic Discrete Event Simulation package, but could feasibly meet your needs.

Better yet the Molecular Modelling ToolKit (MMTK) seems more specialized...

I have used neither, but this sounds like fun. Python, as a language, seems to be in privileged position for use in simulation software, whereby people can script the specific details of their model while relying on the framework for all the common logic, such as scheduling, visualization, monitoring etc. The unknown is how well such toolkits scale when fed with agent counts commensurate with biology models (BTW, how "big" is that?)

Answer 2

Lampps and gromacs are two well known molecular dynamics codes. These codes both have some python based wrapper stuff, but I am not sure how much functionality the wrappers expose. They may not give you enough control over the simulation.

Google for "GromacsWrapper" or google for "lammps" and "pizza.py"

Digital material and ASE are two molecular dynamics codes that expose a lot of functionality, but last time I looked, they were both fairly specialized. They may not allow you to use the force potentials that you want

Google for "digital material" and "cornell" or google for "ase" and dtu

Note to MJV: Normal MD-codes take one time step at a time, and they move all particles in each time step. Most of the time is spend calculating the total force on each atom. This involves iterating over a list of pairs of neighboring atoms. I think the best idea is to do the force calculation and a few more basics in c++ or fortran and then wrap that functionality in python. (But it could be fun to see how far one can get by using numpy matrices)