Avogadro's packmol plugin

When running molecular dynamics (MD) simulations, one often needs to create an initial input structure. It is important that none of the atoms are too close to prevent the coordinates from "blowing up". Many MD programs provide their own utilities (e.g. gromacs' genbox) to create basic systems. For example, a protein with solvent (water) and counter ions. The counter ions are important to make the system charge neutral for algorithms like Ewald summation or the faster PME. While these utilities are often enough, creating complex input files is not always easy.

Luckily, there is a program dedicated to this task: Packmol. It uses simple text input files to define a system. The syntax is very easy to learn but writing the files by hand is a tedious task. Many of the required calculations can be automated and I started writing the Avogadro packmol plugin to do just that. The plugin should allow novice users to create complex systems and should also allow experienced users to make small adjustments (e.g. adding keywords, parameters).

The work flow of the plugin is as follows:
- Configure parameters in a wizard:

The plugin allows you to browse for all required files. This allows all files to be copied to a temporarily location to make sure all files are found. The screen shot shows the wizard for the most basic systems containing solvent and optionally solute.

- Clicking "Generate input file" shows the generated input file

If a new structure is added in the Text Mode tab, the plugin will ask for the file's location when the user clicks run. The syntax highlighting is there to prevent typing errors, it really helps :-)

- Clicking "Run Packmol" starts the program and shows the progress in the Output tab:

Once the program finishes, the resulting file is opened in Avogadro.

The result from the input file seen in the screen shot looks like:


One of the most important calculations required is the number of molecules required to fill a given volume. Given the density, molecular weight and volume, the number of molecules can easily be calculated using Avogadro's number. (This is the same calculation the volume guesser on packmol's site uses AFAIK)

The next step towards more complex systems would be to add solute and counter ions. I'm not completely sure how to subtract the volume taken by solute from the solvents volume. Computing the solute's volume might be an option. Ideas are welcome. :-)

Lipid bilayers are the next kind of system that can be created with packmol. This requires some more knowledge of the packmol syntax but it is still human readable. Before continuing, it might be good to review the bilayer input example from the packmol website. It has comments explaining how to the constraints are used to define the system. The wizard dialog currently looks like:


Currently only 1 "polar solvent" type and 1 "lipid" type are allowed but the idea is that this could allow for mixtures of lipids and so on.

To orient the lipids in the correct way, constraints are needed on specific atoms. This is currently a limitation since a user still has to specify the atom indexes for the polar head and lipophilic tail. In the screen shot below, the indexes have already been filled in.

The required indexes can easily be determined by opening the lipid structure file and enabling the label display type in Avogadro though.

The result from the input above (rendered with POV-Ray):


The next logical step would be to insert a membrane protein between the lipids. This requires the orientation and position of the protein in the membrane to be known. From the protein database files, there is no easy way to determine this. In many cases, the orientations of proteins in membrane (OPM) database provides exactly this. From their website:

The original PDB coordinates are transformed so that the calculated membrane normal coincides with the Z axis, and the origin of coordinates corresponds to the middle of membrane.

This is how the plugin generates the lipid bilayer which allows the protein to be inserted with only a few lines added:

Again, the volume of the protein isn't considered yet. Using slightly lower numbers should work for now. After waiting a while the result should look like:




The source code can be found here.

(The files from the OPM database contain dummy atoms marking where the membrane should be. Make sure to remove these if you try to replicate this)