SMArt

Welcome to Simulation & Modelling Art

importing

import sys
import os
"""
SMArt_fd = os.path.abspath(os.getcwd())
for _ in range(2):
    SMArt_fd = os.path.split(smart_fd)[0]
sys.path.insert(0, smart_fd)
#"""
# if SMArt not in the path, uncomment the code above

import SMArt
data_fd = os.path.abspath(os.path.join(os.path.split(SMArt.__file__)[0], '..', 'doc', 'some_data'))

from SMArt.md import parse_top, parse_ff # functions to parse different files (e.g. topology or force field)
from SMArt.md import FF, Topology, Configuration # examples of some of the most important classes

#help(SMArt) # to get some general info on the package

parsing / writing

GROMOS

# just defining some paths to the data
out_fd = os.path.abspath(os.path.join(os.getcwd(), 'out_data'))
if not os.path.isdir(out_fd):os.mkdir(out_fd)
fd_gromos = os.path.join(data_fd, 'gromos')

Force field

# parsing
ifp_file = os.path.join(fd_gromos, '54a8.ifp')
ff_gr = parse_ff(ifp_file)
# same as
#ff_gr = FF(ifp_file)

# data of any object is stored in _containers
ff_gr._containers

['sys_title',
 'm_type',
 'gr_bonds',
 'gr_angles',
 'gr_impropers',
 'gr_dihedrals',
 'a_type',
 'vdw']

# mass types
ff_gr.m_type

OrderedDict([('1', 1 1.008),
             ('3', 3 13.019),
             ('4', 4 14.027),
             ('5', 5 15.035),
             ('6', 6 16.043),
             ('12', 12 12.011),
             ('14', 14 14.0067),
             ('16', 16 15.9994),
             ('19', 19 18.9984),
             ('23', 23 22.9898),
             ('24', 24 24.305),
             ('28', 28 28.08),
             ('31', 31 30.9738),
             ('32', 32 32.06),
             ('35', 35 35.453),
             ('39', 39 39.948),
             ('40', 40 40.08),
             ('56', 56 55.847),
             ('63', 63 63.546),
             ('65', 65 65.37),
             ('80', 80 79.904)])

ff_gr.m_type['1'] # mass type 1

1 1.008

ff_gr.a_type['2'] # atom type 2

2 OM

ff_gr.gr_bonds['1'] # bond type 2

BondType (15700000.0, 314000.0, 0.1)

Topology

top_file = os.path.join(fd_gromos, 'villin', 'wt.top')
#top_gr = parse_top(top_file)
# same as
top_gr = Topology(top_file)

top_gr.ff._containers # the topology object contains the force field data as well

['a_type', 'gr_bonds', 'gr_angles', 'gr_impropers', 'gr_dihedrals', 'vdw']

top_gr._containers # data stored in:

['sys_title',
 'ff',
 'residues',
 'cg',
 'excl_pair',
 'bonds',
 'angles',
 'impropers',
 'dihedrals']

# write out FF ifp file
out_file = os.path.join(out_fd, 'new_ff.ifp')
ff_gr.write_ff(out_file, format_type = 'gr')

# write out top file
out_file = os.path.join(out_fd, 'new_top.top')
top_gr.write_top(out_file)

some operations with the topology object

len(top_gr.atoms)

389

top_gr.residues['1'].atoms

[1 H1, 2 H2, 3 N, 4 H3, 5 CA, 6 CB, 7 CG, 8 SD, 9 CE, 10 C, 11 O]

# reduce topolgy to first 11 atoms (first residue)
red_top = top_gr.reduce_top([str(i) for i in range(1,12)])

red_top.residues

OrderedDict([('1', <SMArt.md.data_st.Residue at 0x7fdaa40b1fd0>)])

red_top.atoms

OrderedDict([('1', 1 H1),
             ('2', 2 H2),
             ('3', 3 N),
             ('4', 4 H3),
             ('5', 5 CA),
             ('6', 6 CB),
             ('7', 7 CG),
             ('8', 8 SD),
             ('9', 9 CE),
             ('10', 10 C),
             ('11', 11 O)])

out_file = os.path.join(out_fd, 'new_red_top.top')
red_top.write_top(out_file)

GROMACS

Force field

ff_gm = parse_ff('./gromos54a7.ff/forcefield.itp', format_type = 'gm') # default for format_type is 'gr' for GROMOS
# note that the force-field files can be automatically found from the gmx (if installed) - otherwise, provide the full path

list(ff_gm.defines.items())[:10] # list of definitions

[('_FF_GROMOS96', []),
 ('_FF_GROMOS54A7', []),
 ('gb_1', ['0.1000', '1.5700e+07']),
 ('gb_2', ['0.1000', '1.8700e+07']),
 ('gb_3', ['0.1090', '1.2300e+07']),
 ('gb_4', ['0.112', '3.7000e+07']),
 ('gb_5', ['0.1230', '1.6600e+07']),
 ('gb_6', ['0.1250', '1.3400e+07']),
 ('gb_7', ['0.1320', '1.2000e+07']),
 ('gb_8', ['0.1330', '8.8700e+06'])]

ff_gm._containers

['_segments',
 'sys_title',
 'a_type',
 'gm_vdw_normal',
 'gm_vdw_pairs',
 'gm_constraints',
 'gm_bonds',
 'gm_angles',
 'gm_dihedrals']

ff_gm.gm_bonds # bonde types

[BondType (0.204, 5030000.0), BondType (0.198, 640000.0)]

Topology

fd_gromacs = fd_gromos = os.path.join(data_fd, 'gromacs')
top_file = os.path.join(fd_gromacs, 'gromos', 'topol.top') # GROMACS topology using the GROMOS force field
top_gm = parse_top(top_file, format_type = 'gm')

top_gm._containers

['_segments', 'sys_title', 'ff', 'molecule_types', 'molecules']

top_gm.molecules

[Protein_chain_A Protein_chain_A 1]

top_gm.molecule_types

OrderedDict([('Protein_chain_A', Protein_chain_A Protein_chain_A),
             ('SOL', SOL SOL),
             ('CU1', CU1 CU1),
             ('CU', CU CU),
             ('ZN', ZN ZN),
             ('MG', MG MG),
             ('CA', CA CA),
             ('NA', NA NA),
             ('CL', CL CL)])

mol_type = top_gm.molecule_types['Protein_chain_A']
mol_type._containers

['residues',
 'cg',
 'bonds',
 'vdw_pairs',
 'angles',
 'dihedrals',
 'impropers',
 'excl_pair']

os.chdir(out_fd)

ff_gm.write_ff('gmx_ff.itp', from_str = 1, flag_segment_order=1, flag_split_non_bonded = 1, flag_if=1, format_type = 'gm')

top_gm.write_top('gm_top.top', format_type = 'gm')

g2g (converting between GROMOS and GROMACS formats)

# gromos to gromacs
ff_gr.gr2gm()
top_gr.gr2gm()
# gromacs to gromos
ff_gm.gm2gr()
top_gm.gm2gr()

# write out into files
ff_gr.write_ff('gr2gm_ff.itp', format_type = 'gm')
top_gr.write_top('gr2gm_top.top', from_str = 1, flag_segment_order=1, flag_if=1, flag_include = 1, sep_ff2itp='top_ff.itp', flag_defines_first_include = 1, sep_mol2itp=1, format_type='gm')

ff_gm.write_ff('gm2gr_ff.ifp', format_type = 'gr')
top_gm.write_top('gm2gr_top.top', format_type = 'gr')

some additional features

# find rings
mol_type.find_rings()

([177 CD1, 181 CE1, 185 CZ, 183 CE2, 179 CD2, 176 CG],
 [242 CE3,
  246 CZ3,
  248 CH2,
  244 CZ2,
  241 CE2,
  238 CD2,
  236 CD1,
  239 NE1,
  235 CG],
 [377 CD1, 381 CE1, 385 CZ, 383 CE2, 379 CD2, 376 CG],
 [216 CA, 217 CB, 218 CG, 219 CD, 215 N],
 [106 CD1, 110 CE1, 114 CZ, 112 CE2, 108 CD2, 105 CG],
 [62 CD1, 66 CE1, 70 CZ, 68 CE2, 64 CD2, 61 CG])

at = mol_type.atoms[1]
print(list(mol_type.BFS_l(at, 3))) # 3rd neighbours
el, pl = mol_type._generate_atom_excl_pair_list()
pl[at] # list of 1-4 pairs

[12 N, 11 O, 7 CG]

[7 CG, 11 O, 12 N]

print(list(mol_type.BFS_l(at, 1))) # first neighbours

[3 H2, 5 CA, 2 H1, 4 H3]

print(list(mol_type.BFS_d(at, 3))) # all up to 3rd neighbours (including itself)
print([at] + el[at] + pl[at])

[1 N, 3 H2, 5 CA, 2 H1, 4 H3, 10 C, 6 CB, 12 N, 11 O, 7 CG]
[1 N, 2 H1, 3 H2, 4 H3, 5 CA, 6 CB, 10 C, 7 CG, 11 O, 12 N]

#"""
# remove all backed-up files
import glob
for f_path in glob.glob('#bck*'):
    os.remove(f_path)
#"""