[1]:
# importing and preparing some variables
import sys
import os
import glob
import re
"""
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'))
out_fd = os.path.abspath(os.path.join(os.getcwd(), 'PTP_out_data'))
if not os.path.isdir(out_fd):os.mkdir(out_fd)
Perturbation topologies¶
[2]:
from SMArt.md import parse_top
from SMArt import alchemy
[3]:
# GROMACS
# parse the topologies
top_wt_file = os.path.join(data_fd, 'gromacs', 'LYSH.top')
top_wt = parse_top(top_wt_file, format_type='gm')
top_PTM_file = os.path.join(data_fd, 'gromacs', 'K3C.top')
top_PTM = parse_top(top_PTM_file, format_type='gm')
mt_wt, mt_PTP = top_wt.molecule_types['mol_1'], top_PTM.molecule_types['mol_1']
# run the MCS algorithm
mcs = alchemy.point_mutation(mt_wt, mt_PTP)
mcs.solutions[0]
[3]:
top
0 1
atom 0 6 CA 6 CA
1 7 CB 7 CB
2 8 CG 8 CG
3 9 CD 9 CD
4 15 C 15 C
5 4 N 4 N
6 10 CE 10 CE
7 11 NZ 11 NZ
8 17 NTE 17 NTE
9 2 CN1 2 CN1
10 16 O 16 O
11 5 H 5 H
12 12 HZ1 14 CH3
13 14 HZ3 12 CH1
14 13 HZ2 13 CH2
15 19 H2 19 H2
16 18 H1 18 H1
17 3 ON2 3 ON2
18 1 CN2 1 CN2
[4]:
mcs.solutions[1] # an alternative solution (6 effectively identical solutions)
[4]:
top
0 1
atom 0 6 CA 6 CA
1 7 CB 7 CB
2 8 CG 8 CG
3 9 CD 9 CD
4 15 C 15 C
5 4 N 4 N
6 10 CE 10 CE
7 11 NZ 11 NZ
8 17 NTE 17 NTE
9 2 CN1 2 CN1
10 16 O 16 O
11 5 H 5 H
12 12 HZ1 14 CH3
13 14 HZ3 13 CH2
14 13 HZ2 12 CH1
15 19 H2 19 H2
16 18 H1 18 H1
17 3 ON2 3 ON2
18 1 CN2 1 CN2
[5]:
sol = alchemy.generate_2state_top(mcs, ff_dumm = top_wt.get_DUM_type) # solution = 0 by default
out_itp = os.path.join(out_fd, 'toptp.itp')
sol.toptp.write_itp(out_itp, flag_generate_excl_pairs = True)
[6]:
# GROMOS
# parse topologies
top_wt_file = os.path.join(data_fd, 'gromos', 'LYS_K3C_KAC', 'single_AA', 'LYS.top')
top_wt = parse_top(top_wt_file)
top_PTM_file = os.path.join(data_fd, 'gromos', 'LYS_K3C_KAC', 'single_AA', 'K3C.top')
top_PTM = parse_top(top_PTM_file)
# run the MCS algorithm
mcs = alchemy.point_mutation(top_wt, top_PTM)
sol = alchemy.generate_2state_top(mcs) # solution = 0 by default
sol.toptp.write_top(os.path.join(out_fd, 'toptp.top'))
sol.toptp.write_ptp(os.path.join(out_fd, 'toptp.ptp'))
#check_names(sol, at_map)
sol
[6]:
top
0 1
atom 0 6 CA 6 CA
1 7 CB 7 CB
2 8 CG 8 CG
3 9 CD 9 CD
4 4 N 4 N
5 15 C 15 C
6 10 CE 10 CE
7 11 NZ 11 NZ
8 2 CN1 2 CN1
9 17 NTE 17 NTE
10 5 H 5 H
11 16 O 16 O
12 13 HZ2 14 CH3
13 14 HZ3 13 CH2
14 12 HZ1 12 CH1
15 3 ON2 3 ON2
16 1 CN2 1 CN2
17 18 H1 19 H2
18 19 H2 18 H1
Alternative MCS search options¶
no bond perturbations allowed¶
[7]:
# for simplicity, let's reduce the topologies to the sidechains
top_wt = top_wt.reduce_top([str(i) for i in range(6,15)])
top_PTM = top_PTM.reduce_top([str(i) for i in range(6,15)])
mcs = alchemy.point_mutation(top_wt, top_PTM, flag_top_prune = 'bond') # makes sure no bonds are perturbed
mcs.solutions[0]
[7]:
top
0 1
atom 0 8 CG 8 CG
1 9 CD 9 CD
2 10 CE 10 CE
3 11 NZ 11 NZ
4 7 CB 7 CB
5 13 HZ2 DUM
6 12 HZ1 DUM
7 14 HZ3 DUM
8 6 CA 6 CA
9 DUM 12 CH1
10 DUM 14 CH3
11 DUM 13 CH2
add RMSD¶
[8]:
import numpy as np
coord = np.random.rand(len(top_wt.atoms) * 3).reshape(len(top_wt.atoms), 3) # random coordinates
coord2 = coord.copy()
mcs = alchemy.point_mutation(top_wt, top_PTM, add_RMSD='simple', coords = [coord, coord2])
mcs.solutions[0]
[8]:
top
0 1
atom 0 8 CG 8 CG
1 9 CD 9 CD
2 10 CE 10 CE
3 11 NZ 11 NZ
4 7 CB 7 CB
5 13 HZ2 13 CH2
6 12 HZ1 12 CH1
7 14 HZ3 14 CH3
8 6 CA 6 CA
[9]:
# let's swapp the the coordinates of last 2 atoms
print(coord2[-2:])
coord2[-2:] = coord2[-2:][::-1]
print('after swapping')
print(coord2[-2:])
[[0.40479395 0.89759386 0.03049702]
[0.15758507 0.15801852 0.00853602]]
after swapping
[[0.15758507 0.15801852 0.00853602]
[0.40479395 0.89759386 0.03049702]]
[10]:
mcs = alchemy.point_mutation(top_wt, top_PTM, add_RMSD='simple', coords = [coord, coord2])
mcs.solutions[0]
[10]:
top
0 1
atom 0 8 CG 8 CG
1 9 CD 9 CD
2 10 CE 10 CE
3 11 NZ 11 NZ
4 7 CB 7 CB
5 13 HZ2 14 CH3
6 12 HZ1 12 CH1
7 14 HZ3 13 CH2
8 6 CA 6 CA
other options include:
flag_partial_ring - defines if partial ring match is allowed or not (only for fused rings) - True by default
max_partial_ring_match - defines how many ring to non-ring atoms are allowed to be matched (2 by default)
dihedral_match_v - procedure to match proper dihedrals - 1 (default) based on the middle 2 atoms of a dihedral, or 2 match based on all atoms
use_ptp_flags - defines if the multiplicity of proper dihedrals is allowed to be perturbed (for GROMACS)
EDS topology¶
[11]:
TOPs_path = os.path.join(data_fd, 'gromos', 'GRA2', 'tops', '*top')
TOPs = glob.glob(TOPs_path)
tops = [parse_top(fpath) for fpath in TOPs]
mcs = alchemy.get_EDS(*tops[:5]) # first 5 tops; faster than all 16 topologies (line below)
#mcs = alchemy.get_EDS(*tops)
sol_EDS, state_names = alchemy.generate_EDS_top(mcs)
sol_EDS.toptp.write_top(os.path.join(out_fd, 'EDS.top'))
sol_EDS.toptp.write_EDS(os.path.join(out_fd, 'EDS.ptp'), state_names=state_names)
sol_EDS
[11]:
top
0 1 2 3 4
atom 0 3 C4 3 C4 3 C4 3 C4 3 C4
1 20 C25 20 C25 20 C25 20 C25 20 C25
2 1 C1 1 C1 1 C1 1 C1 1 C1
3 9 C10 9 C10 9 C10 9 C10 9 C10
4 5 C6 5 C6 5 C6 5 C6 5 C6
5 7 C8 7 C8 7 C8 7 C8 7 C8
6 17 S22 17 S22 17 S22 17 S22 17 S22
7 15 N20 15 N20 15 N20 15 N20 15 N20
8 10 N11 10 N11 10 N11 10 N11 10 N11
9 12 C13 12 C13 12 C13 12 C13 12 C13
10 11 C12 11 C12 11 C12 11 C12 11 C12
11 13 C14 14 C19 13 C14 13 C14 14 C19
12 4 H5 4 H5 DUM 4 H5 4 H5
13 2 H2 DUM 2 H2 2 H2 2 H2
14 6 H7 6 H7 6 H7 DUM 6 H7
15 8 H9 8 H9 8 H9 8 H9 DUM
16 19 O24 18 O23 18 O23 18 O23 19 O24
17 18 O23 19 O24 19 O24 19 O24 18 O23
18 16 H21 16 H21 16 H21 16 H21 16 H21
19 14 C19 13 C14 14 C19 14 C19 13 C14
20 DUM 2 F2 DUM DUM DUM
21 DUM DUM 4 F5 DUM DUM
22 DUM DUM DUM 6 F7 DUM
23 DUM DUM DUM DUM 8 F9
[ ]: