tutorial-icon

Tutorial 1c: Analyzing the toy project

0) Download the finished toy project

If you didn’t follow all the setup and simulation run, you can download the completed toy project started in Tutorial 1a with the simulation results from here. Extract the folder to proceed and move to the project directory (pep_amber) to start with the analysis.

1) Check project status

Execute the following command inside the project folder to get a description of the replicas (independent simulations) contained and their simulation status. At this stage, we should find a True in Min, Eq and Prod tags:

> mdmix info project
 ==========================================================
 || pyMDMix User Interface ||
 ==========================================================
 || Author: Daniel Alvarez-Garcia ||
 || Version : 0.1 ||
 ==========================================================
------------------------------
PROJECT pep_amber INFO
------------------------------
SYSTEMS:
========
SYSTEM:pep
REPLICAS
========
REPLICA:ETA_1 system:pep_ETA nanos:1 restrMode:HA restrForce:0.010 Min:True Eq:True Prod:True Align:False
REPLICA:WAT_1 system:pep_WAT nanos:1 restrMode:HA restrForce:0.010 Min:True Eq:True Prod:True Align:False
------------------------------

2) Align trajectory and quality check

If all is correct, let’s align the trajectory for all the replicas in the project. The reference structure used for aligning the trajectory is a pdb file ending with _ref.pdb. In our example pep_ETA_ETA_1_ref.pdb. This file is exactly the same as pep_ref.pdb in main project folder.

The process involves the automatic detection of all residue numbers part of the protein. The backbone atoms of these residues will be used for the alignment process. If you wish to use different residues for the alignment process (e.g. align just a subdomain), read the alignment section in the documentation.

> mdmix analyze align all
==========================================================
 || pyMDMix User Interface ||
 ==========================================================
 || Author: Daniel Alvarez-Garcia ||
 || Version : 0.1 ||
 ==========================================================

Selected replica names: ['ETA_1', 'WAT_1']
INFO Running alignment for replicas [ETA_1 Replica, WAT_1 Replica]
INFO Running alignment of trajectory in replica WAT_1
INFO Running alignment of trajectory in replica ETA_1
DONE
Total execution time: 8.708s

The alignment of 1ns for each replica in this small system took just 8.7 seconds in my machine. If we check again project status, we should find a Align:True for each of the replicas.

Inside each replica folder (e.g. MD/ETA_1/), a new folder align will be written with the resulting trajectory files, input ptraj scripts, one pdb file for each job run containing the protein average, and two text files with rmsd data for each job run (one for the backbone atoms and one for all heavy atoms).If you remember, when we set up the project, we modified the defaults and each file was saving 500ps. We chose 1ns simulation and thus we will have two files for one nanosecond of simulation.

Eg. this is the content of MD/ETA_1/align folder:

md1_bb_rmsd.out
md1_ha_rmsd.out
md1.nc
md1.ptraj
md1_ptraj.log
md2_bb_rmsd.out
md2_ha_rmsd.out
md2.nc
md2.ptraj
md2_ptraj.log
prot_avg_1.pdb
prot_avg_2.pdb

3) Build density grids

Once trajectory has been aligned, we can calculate the density maps for each probe in the solvents we have simulated. We can query the solvent probe names. For instance:

> mdmix info solvents
 ==========================================================
 || pyMDMix User Interface ||
 ==========================================================
 || Author: Daniel Alvarez-Garcia ||
 || Version : 0.1 ||
 ==========================================================
--------------------------------------------------
SOLVENT DB in USE: /home/daniel/.local/lib/python2.7/site-packages/pyMDMix-0.1-py2.7.egg/pyMDMix/data/solventlib/SOLVENTS.db
Package SOLVENT DB: /home/daniel/.local/lib/python2.7/site-packages/pyMDMix-0.1-py2.7.egg/pyMDMix/data/solventlib/SOLVENTS.db
SOLVENTS Directory: /home/daniel/.local/lib/python2.7/site-packages/pyMDMix-0.1-py2.7.egg/pyMDMix/data/solventlib

Solvents:
 SOLVENT: MOH INFO: Methanol 20% mixture BOXUNIT: MOHWAT20 PROBES: MOH_WAT, MOH_OH, MOH_CT
 SOLVENT: ANT INFO: Acetonitrile 20% mixture BOXUNIT: ANTWAT20 PROBES: ANT_C, ANT_WAT, ANT_N
 SOLVENT: ION INFO: Ammonium Acetate mixture (Charged mixture) BOXUNIT: CN3COOWAT PROBES: ION_NEG, ION_WAT, ION_POS
 SOLVENT: ETA INFO: Ethanol 20% mixture BOXUNIT: ETAWAT20 PROBES: ETA_CT, ETA_WAT, ETA_OH
 SOLVENT: ISO INFO: Isopropanol 20% mixture BOXUNIT: ISOWAT20 PROBES: ISO_WAT, ISO_OH, ISO_CT
 SOLVENT: ISO5 INFO: Isopropanol 5% mixture BOXUNIT: ISOWAT5 PROBES: ISO5_WAT, ISO5_OH, ISO5_CT
 SOLVENT: WAT INFO: Water box BOXUNIT: TIP3PBOX PROBES: WAT_WAT
 SOLVENT: MAM INFO: Acetamide 20% mixture BOXUNIT: MAMWAT20 PROBES: MAM_N, MAM_WAT, MAM_O, MAM_CT
--------------------------------------------------

Ethanol mixture (ETA) has three probes: ETA_CT, ETA_WAT and ETA_OH (being the methyl carbon, the water oxygen and the ethanol hydroxyl oxygen respectively). WAT mixture, which is actually not a mixture but just a box with water, has WAT_WAT probe.

For the current example we will calculate all probes for each solvent box plus a center of mass density map for each molecule type in the mixture:

> mdmix analyze density all --com
> mdmix analyze density all --com -C4 # same using 4 cpus

With this command we select all replicas in the project, calculate all density maps for all probes plus the center of mass density maps.

We could run just the calculation for all replicas with ETA solvent like this:

> mdmix analyze density bysolvent -s ETA --com -C3 # using 3 cpus
==========================================================
 || pyMDMix User Interface ||
 ==========================================================
 || Author: Daniel Alvarez-Garcia ||
 || Version : 0.1 ||
 ==========================================================

Selected replica names: ['ETA_1']
INFO Preparing analysis...
INFO Setting up density grids calculation for replica ETA_1
INFO Calculating density for probes: ['ETA_CT', 'ETA_WAT', 'ETA_OH', 'ETA_COM', 'WAT_COM']
INFO Ready to calculate
INFO --------------------------------------------------
INFO Running ETA_1 analysis...
INFO Working on file /home/daniel/Dropbox/WORK/mdmixDocAndMore/pep_amber/MD/ETA_1/align/md1.nc
INFO Working on file /home/daniel/Dropbox/WORK/mdmixDocAndMore/pep_amber/MD/ETA_1/align/md2.nc
INFO Processing density results for replica ETA_1
INFO Writing density grid dgrids/ETA_1_ETA_CT.dx
INFO Writing density grid dgrids/ETA_1_ETA_OH.dx
INFO Writing density grid dgrids/ETA_1_ETA_COM.dx
INFO Writing density grid dgrids/ETA_1_ETA_WAT.dx
INFO Writing density grid dgrids/ETA_1_WAT_COM.dx
INFO --------------------------------------------------
DONE
Total execution time: 66.731s

Alternatively, we could choose just some of the probes (e.g. ETA_CT and ETA_OH only):

> mdmix analyze density bysolvent -s ETA -P ETA_CT ETA_OH

Inside the replica folder (pep_amber/MD/ETA_1) a new directory dgrids/ will be created with one grid with format DX for each probe.

These grids can be visualized with VMD, PyMOL and the majority of molecular visualization programs. Use pep_ref.pdb in the main project folder as the reference structure for visualizing the results (the same used for aligning the trajectory).

Interpretation of the density grids obtained with this toy project are non-sense due to the short time of simulation and the simple peptide structure simulated.

4) Convert density maps to energy grids

Last step is converting the recently calculated density maps (or grids) to energy values using inverse Boltzmann equation as described in the documentation. There are two flavors of this function:

  • Convert single replicas: Each replica will have its own energy grids for the probes calculated.
  • Merge replicas and convert: Several replicas (same as independent simulations with same solvent) for which density maps have been calculated, can be merged to obtain more converged results. The process involves the addition of the density maps before conversion. This is done by selecting more than one replica of the same solvent. In this example, we have just one replica for each solvent and thus cannot do this analysis.

Converting single replicas

Again in project folder execute the following command to convert all density grids inside replica named ETA_1 to free energy maps (see Analysis Guide for theory and more examples).

> mdmix analyse energy byname -s ETA_1
 ==========================================================
 || pyMDMix User Interface ||
 ==========================================================
 || Author: Daniel Alvarez-Garcia ||
 || Version : 0.1 ||
 ==========================================================

Selected replica names: ['ETA_1']
INFO Converting grids for replica ETA_1
INFO Probe ETA_OH...
INFO Energy grid /home/daniel/Dropbox/WORK/mdmixDocAndMore/pep_amber/MD/ETA_1/egrids/ETA_OH_DG0.dx saved. OK.
INFO Probe ETA_CT...
INFO Energy grid /home/daniel/Dropbox/WORK/mdmixDocAndMore/pep_amber/MD/ETA_1/egrids/ETA_CT_DG0.dx saved. OK.
INFO Probe WAT_COM...
INFO Energy grid /home/daniel/Dropbox/WORK/mdmixDocAndMore/pep_amber/MD/ETA_1/egrids/WAT_COM_DG.dx saved. OK.
INFO Probe ETA_WAT...
INFO Energy grid /home/daniel/Dropbox/WORK/mdmixDocAndMore/pep_amber/MD/ETA_1/egrids/ETA_WAT_DG.dx saved. OK.
INFO Probe ETA_COM...
INFO Energy grid /home/daniel/Dropbox/WORK/mdmixDocAndMore/pep_amber/MD/ETA_1/egrids/ETA_COM_DG0.dx saved. OK.
INFO Done with all probes for replicas ['ETA_1']
Total execution time: 10.280s

Inside the replica folder (MD/ETA_1) a new folder egrids will be created containing the resulting grids after energy conversion.  These grids can be visualized with any molecular visualization program (pymol, VMD, etc.). When file name ends with DG0, it means that the volume correction has been applied to obtain the results. If the file name ends with DG, no volumetric correction has been done (as is the case for water grids which do not require such standard energy correction) – see Analysis Guide for details.

5) Otaining PDB files from grids

asdf

Posted in Tutorial and tagged .