Archive for February, 2010

Simple liquid simulations using OpenMD

In this example, we’ll build a simulation of a simple liquid (methanol) starting with a structure file (in XYZ format). Getting your molecule of choice into OpenMD is never a black box procedure, and will almost always require some hand adjustment of the input files. Assuming you’ve gotten a version of OpenMD built and installed (along with all of the prerequisite software), we’ll walk you through building the liquid simulation in the steps below.

  1. Start with a good structure for your molecule. If you are familiar with Avogadro, build the methanol structure (CH3OH), set up the force field, and optimize the structure. Save the structure as methanol.xyz. Alternatively, download the methanol.xyz structure to your working directory.
  2. Use the atom2omd program to convert the structure into a format that can be read by OpenMD:
    atom2omd -ixyz methanol.xyz

    This command will create an incomplete OpenMD file called methanol.omd that must be edited before it can be used.

  3. OpenMD can use a number of force fields, but in this example, we’ll use the Amber force field. If you are using this force field and are starting from an XYZ file or non-standard PDB file, you must edit the atom types. In the methanol.omd file:
    • change the atom typing for the methyl carbon from C3 to CT
    • change the O3 to OH
    • The hydrogens on the carbon should also be changed from HC to H1
    • The hydroxyl hydrogen can be left alone.
  4. At this point it is also a good idea to change the name of the molecule to something descriptive (perhaps “methanol”). This should be done in two places; once in the molecule description and another time in the component block.
  5. Before the simulation can run, add a forceField line after the component block:
    forceField = "Amber";

    At this stage, you should be able to run OpenMD on the file to check to make sure your hand-crafted atom typing can be matched up with types known by the force field:

    openmd methanol.omd

    If there are any problems, correct any unknown atom types, and repeat until you get an error about the “Integrator Factory”.

  6. Next, we’ll build a lattice of methanol molecules using this initial structure as a starting point. The density of liquid methanol is roughly 0.7918 g cm-3, so we’ll build a simple box of methanol molecules using the command:
    simpleBuilder -o liquid.omd --density=0.7918 --nx=3 --ny=3 --nz=3 methanol.omd

    This command creates a new system, liquid.omd which contains 108 copies of the methanol molecule arranged in a simple FCC lattice. FCC has 4 molecules in the unit cell, so the total number of molecules = 4 * 3 * 3 * 3 = 108. The molecules are packed at a distance commensurate with their liquid state density.

  7. To visualize what the system looks like at this stage, you can run:
    Dump2XYZ -i liquid.omd

    to create a file called liquid.xyz. This file can be viewed in VMD, Jmol, or any other chemical structure viewer.

  8. Add the following lines below the forceField line of the liquid.omd file:
    ensemble = NVT;
    cutoffMethod = "shifted_force";
    electrostaticScreeningMethod = "damped";
    cutoffRadius = 9;
    dampingAlpha = 0.18;
    targetTemp = 300;
    tauThermostat = 1000;
    dt = 1.0;
    runTime = 1e3;
    tempSet = "false";
    sampleTime = 100;
    statusTime = 10;
  9. Initial configurations that are created from bare structures typically have no velocity information. To give an initial kick to the atoms (i.e. to sample the velocities from a Maxwell-Boltzmann distribution), you can use the following command:
    thermalizer -o warm.omd -t 300 liquid.omd

    This creates a new OpenMD file called warm.omd which has initial velocity information.

  10. At this stage, a simple simulation can be run:
    openmd warm.omd
  11. This should complete relatively quickly, and should create a warm.stat file as well as a warm.dump file containing the actual trajectory.
  12. To view the contents of the trajectory file, you’ll need to convert the dump file into something another program can visualize:
    Dump2XYZ -i warm.dump

    will create a new file warm.xyz that can be viewed in VMD and many other chemical structure viewers.

  13. The “End-of-Run” file warm.eor can be re-purposed as the starting point for a new simulation:
    cp warm.eor  equilibration.omd

    Edit the equilibration.omd file, and change parameters you’d like to change before running openmd on the new file.

Converting a protein structure for use with OpenMD

One of the major tasks facing a new user of OpenMD is getting their molecule of choice into a form that can be read by the program. This is never a black box procedure, and will almost always require some hand adjustment of the input files. Assuming you’ve gotten a version of OpenMD built and installed (along with all of the prerequisite software), we’ll try to walk you through the steps below.

  1. Start with a molecular structure in pdb or xyz formats. A good program for generating and optimizing starting structures is Avogadro, which is a wonderful open source molecular editor.
  2. Alternatively (and this is the route we’ll take here), you can download a protein structure directly from the Protein Data Bank. We’ve picked a short pentapeptide (the neurotransmitter met-enkephalin) to use as an example. The PDB-ID for this protein is 1PLW. Download the text version of the PDB file to your working directory.
  3. Use the atom2md program to convert the structure into a format that can be read by OpenMD:
    atom2md -ipdb 1PLW.pdb

    This command will create an incomplete OpenMD file called 1PLW.md that must be edited before it can be used.

  4. OpenMD can use the Amber force field for protein simulations. If you are using this force field, the N-terminal and C-terminal ends of the protein should be modified to use the correct atom types. This peptide has the sequence “Tyr-Gly-Gly-Phe-Met”, so Tyrosine is the N-terminal end and Methionine is the C-terminal end.

    In the MD file, the atom typing for the first residue must be changed from TYR-X to NTYR-X while for the last residue, it must be changed from MET-X to CMET-X. The special N- and C-terminal atom types usually only apply to the N, CA, C, O, OXT, HN, and HA base types. Other atom types can usually be left alone.

  5. Some PDB files specify more details on the atom name record than the Amber force field can use. In particular, the PHE-CD1 and PHE-CD2 types should be down-specified to PHE-CD. Similar down-specification can be done for PHE-CE1, PHE-CE2, PHE-HD1, PHE-HD2, PHE-HE1, PHE-HE2 atom types.
  6. At this point it is also a good idea to change the name of the molecule to something descriptive (perhaps “metenk”). This should be done in two places; once in the molecule description and another time in the component block.
  7. Before the simulation can run, add a forceField line after the component block:
    forceField = "Amber";

    At this stage, you should be able to run OpenMD on the file to check to make sure your hand-crafted atom typing can be matched up with types known by the force field:

    openmd 1PLW.md

    Correct any unknown atom types, and repeat until you get an error about the “Integrator Factory”.

  8. Add the following lines below the forceField line:
    ensemble = NVT;
    cutoffMethod = "shifted_force";
    electrostaticScreeningMethod = "damped";
    cutoffRadius = 10;
    dampingAlpha = 0.18;
    targetTemp = 300;
    tauThermostat = 1000;
    dt = 1.0;
    runTime = 1e3;
    tempSet = "false";
    sampleTime = 100;
    statusTime = 10;
  9. The default size for the periodic box or Hmat is typically a tight bounding box around the structure of the protein. Edit the Hmat line in the FrameData block to be at least twice the cutoff radius in each dimension. A good choice for a small bare protein might be:
            Hmat: {{ 30, 0, 0 }, { 0, 30, 0 }, { 0, 0, 30 }}

  10. Initial configurations that are created from crystal structures usually have no velocity information. To give an initial kick to the atoms (i.e. to sample the velocities from a Maxwell-Boltzmann distribution), you can use the following command:
    thermalizer -o warm.md -t 300 1PLW.md

    This creates a new OpenMD file called warm.md which has initial velocity information.

  11. At this stage, a simple simulation can be run:
    openmd warm.md

  12. This should complete relatively quickly, and should create a warm.stat file as well as a warm.dump file containing the actual trajectory.
  13. To view the contents of the trajectory file, you’ll need to convert the dump file into something another program can visualize:
    Dump2XYZ -i warm.dump

    will create a new file warm.xyz that can be viewed in VMD and many other chemical structure viewers.

  14. The “End-of-Run” file warm.eor can be re-purposed as the starting point for a new simulation:
    cp warm.eor  equilibration.md

    Edit the equilibration.md file, and change parameters you’d like to change before running openmd on the new file.