Energy decomposition with sander and igb8 solvent model¶

require¶

pytraj/cpptraj: http://amber-md.github.io/pytraj/latest/installation.html#from-source-code
```
  git clone https://github.com/Amber-MD/pytraj
  cd pytraj
  python setup.py install
```
libsander (http://ambermd.org/doc12/Amber15.pdf): If you have AmberTools15 and already make install, it is likely that you already have libsander

In [1]:

# load pytraj and load trajectories

import warnings
warnings.filterwarnings('ignore', category=DeprecationWarning)

import pytraj as pt

traj = pt.iterload('tz2.nc', 'tz2.parm7')
print(traj)

pytraj.TrajectoryIterator, 101 frames: 
Size: 0.000503 (GB)
<Topology: 223 atoms, 13 residues, 1 mols, non-PBC>

In [2]:

# compute the energy with igb=8 (GBneck2 solvation model)
data = pt.energy_decomposition(traj, igb=8)
data

Out[2]:

OrderedDict([('disp', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('angle',
              array([ 128.54514846,  105.06494493,  103.52028376, ...,   99.46042639,
                      103.49664937,   98.18016593])),
             ('polar', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('elec_14',
              array([ 720.65923256,  731.53353441,  717.35440566, ...,  736.47531336,
                      732.27034385,  726.61630671])),
             ('angle_ub', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('rism', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('vdw_14',
              array([ 52.55761103,  50.03499738,  54.29855009, ...,  50.42610919,
                      48.27053805,  51.9298654 ])),
             ('constraint', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('surf', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('vdw',
              array([-44.9811844 , -43.77030705, -52.33589852, ..., -70.01048987,
                     -66.01297848, -67.78679498])),
             ('les', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('tot',
              array([-225.72081333, -259.3103916 , -265.22164221, ..., -282.72357777,
                     -287.55610458, -293.18137453])),
             ('cmap', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('bond',
              array([ 0.01531366,  0.01358172,  0.01252056, ...,  0.0154962 ,
                      0.01510391,  0.0175098 ])),
             ('elec',
              array([-781.59559933, -807.4891344 , -741.17534022, ..., -835.16025964,
                     -805.08867184, -825.97031683])),
             ('pb', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('imp', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('dvdl', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('ct', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('dihedral',
              array([ 111.61132908,  105.39241304,   93.03084989, ...,  101.81038715,
                       94.63365271,   92.27451698])),
             ('amd_boost', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('scf', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('gb',
              array([-412.53266438, -400.09042163, -439.92701343, ..., -365.74056055,
                     -395.14074215, -368.44262754])),
             ('hbond', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('emap', array([ 0.,  0.,  0., ...,  0.,  0.,  0.])),
             ('noe', array([ 0.,  0.,  0., ...,  0.,  0.,  0.]))])

In [3]:

# get energy you want
print('potential energy', data['tot'])

potential energy [-225.72081333 -259.3103916  -265.22164221 ..., -282.72357777 -287.55610458
 -293.18137453]

In [4]:

# solvation energy (igb=8)
print('solvation energy', data['gb'])

solvation energy [-412.53266438 -400.09042163 -439.92701343 ..., -365.74056055 -395.14074215
 -368.44262754]

In [5]:

# you can also get other energies.
# for example: data['dihedral']
list(data.keys())

Out[5]:

['disp',
 'angle',
 'polar',
 'elec_14',
 'angle_ub',
 'rism',
 'vdw_14',
 'constraint',
 'surf',
 'vdw',
 'les',
 'tot',
 'cmap',
 'bond',
 'elec',
 'pb',
 'imp',
 'dvdl',
 'ct',
 'dihedral',
 'amd_boost',
 'scf',
 'gb',
 'hbond',
 'emap',
 'noe']

Parallel calculation¶

In [6]:

# serial: pt.energy_decomposition(traj, igb=8)

# parallel
data2 = pt.pmap(func=pt.energy_decomposition, traj=traj, igb=8, n_cores=4)
data2['gb']

Out[6]:

array([-412.53266438, -400.09042163, -439.92701343, ..., -365.74056055,
       -395.14074215, -368.44262754])

(energy_decomposition.ipynb; energy_decomposition_evaluated.ipynb; energy_decomposition.py)