Skip to content

Generating a periodic system for FHI-aims (full system)

Objective: Turn the full periodic MD system into a geometry file suitable for FHI-aims to be simulated with periodic boundary conditions.

Tip

The current tutorial exports the complete classical MD system and introduces the basic steps that we will also use when we are exporting a subsystem.

Protocol

The following commands can be used interactively in the python or ipython interpreter or a Jupyter notebook. Alternatively, combine them in a script extract_all_pbc.py and run the script with python extract_all_pbc.py.

Set-up

Import MDAnalysis and the transformations module: 

import MDAnalysis as mda
from MDAnalysis import transformations

Load the topology (TPR) and trajectory (TRR) of the classical MD simulation into the universe u: 

u = mda.Universe("./files/MD/NPT/md.tpr", "./files/MD/NPT/md.trr")

Change atom names to elements because these names are written to the output file and FHI-aims wants elements, not force field atom names:

u.atoms.names = u.atoms.elements

Note

With TPR files, information is available for elements. With other input files, MDAnalysis will have to guess the elements.

Create AtomGroup instances for the water molecules (water) and the capped lysine (solute) by using selection by residue name

water = u.atoms.select_atoms("resname SOL")
solute = u.atoms.select_atoms("resname ACE LYSN NME")

We can check that we partitioned the whole system into solute and water with a simple assertion that the combined AtomGroup of solute + water is the same as all the atoms in the system: 

assert u.atoms == solute + water

If the assertion fails then you should check your selections.

Image of the classical MD system on the last frame of the trajectory. Water molecules are shown in CPK and the solute in VDW representation

Classical MD system on the last frame of the trajectory. Water molecules (AtomGroup `water`) are shown in CPK and the solute (AtomGroup `solute`) in VDW representation; the unit cell is shown as a blue cubic box. (Note that PBC artifacts have been fixed in this image to show the selections more clearly. See the next section for what this system *actually* looks like.)

Centering

The objective is to generate a trajectory where

  • the peptide is at the center of the simulation cell
  • water is packed around the peptide
  • water molecules are within periodic boundaries.

We use on-the-fly transformations to rewrite the coordinates of the trajectory while reading the TRR file from the disk. We create a workflow of functions that is run between reading the coordinates from file and storing them in the Timestep instance for the rest of MDAnalysis to see.

workflow = [transformations.unwrap(u.atoms),
            transformations.center_in_box(solute, center='geometry'),
            transformations.wrap(water, compound='atoms')]
u.trajectory.add_transformations(*workflow)
  1. unwrap(u.atoms) unwraps all molecules over the box boundaries if necessary.
  2. center_in_box(solute, center='geometry') shifts all coordinates so that the center of the box coincides with the center of geometry of the AtomGroup solute (which we selected previously).
  3. wrap(water, compound='atoms') ensures that all water atoms are inside the box that was centered on the solute in transformation 2.

We then add the workflow to the u.trajectory instance and activate it for all future use of the universe instance.

Image of the classical MD system on the last frame of the trajectory with molecules broken across periodic boundaries Image of the classical MD system on the last frame with the system centered on the lysine solute and all atoms packed inside periodic boundaries.

Top: System with molecules broken across the periodic boundary, indicated by long bonds. Bottom: System after centering on the solute (lysine in VDW representation) and wrapping all atoms of water molecules back into the box of the periodic unitcell. Water molecules are broken across the periodic boundary, which is the correct behavior.

Select trajectory frame

MDAnalysis keeps track of the current frame of the trajectory. At the beginning, the frame "cursor" points to frame 0.

Info

MDAnalysis uses zero-based indexing for frames.

You can move to any frame between 0 and u.trajectory.n_frames - 1 by indexing the u.trajectory instance, i.e., u.trajectory[frame].

To move to the last frame of trajectory use 

u.trajectory[-1]

All AtomGroups will now read data from the last trajectory frame (suitably transformed).

Output to geometry.in

We want to write the our system to a FHI-aims geometry.in file (MDAnalysis supports reading and writing the FHI-aims file format). We simply take all atoms of the system from the u.atoms attribute of the Universe:

u.atoms.write("geometry.in")

We can also write a file in PDB format for visualization in tools that may not understand the FHI-aims format (ignore warnings about attributes that are set to default values to satisfy the PDB standard):

u.atoms.write("geometry.pdb")
Suppressing warnings in Python

You can temporarily silence the warnings with the warnings filter:

import warnings
with warnings.catch_warnings():
    warnings.simplefilter("ignore", UserWarning)
    smaller_system.write("geometry.pdb")

Complete script

Save this script as extract_all_pbc.py:

import MDAnalysis as mda
from MDAnalysis import transformations

u = mda.Universe("./files/MD/NPT/md.tpr", "./files/MD/NPT/md.trr")
u.atoms.names = u.atoms.elements
water = u.atoms.select_atoms("resname SOL")
solute = u.atoms.select_atoms("resname ACE LYSN NME")
assert u.atoms == solute + water

workflow = [transformations.unwrap(u.atoms),
            transformations.center_in_box(solute, center='geometry'),
            transformations.wrap(water, compound='atoms')]
u.trajectory.add_transformations(*workflow)

u.trajectory[-1]

u.atoms.write("geometry.in")
u.atoms.write("geometry.pdb")

Adjust the input files and the selection criteria and then run it with

python extract_all_pbc.py