diff --git a/examples/openmm/cv_Q/abf.py b/examples/openmm/cv_Q/abf.py index 4929af75..73d6ba9b 100644 --- a/examples/openmm/cv_Q/abf.py +++ b/examples/openmm/cv_Q/abf.py @@ -1,5 +1,8 @@ #!/usr/bin/env python +from itertools import product + +import networkx as nx import numpy as np import openmm as mm import openmm.app as app @@ -11,6 +14,69 @@ from pysages.colvars.contacts import NativeContactFraction from pysages.methods import ABF, HistogramLogger +AUGC = ["A", "U", "G", "C"] + + +def init_graph(top): + graph = nx.Graph() + rna_indices = [] + full_id2rna_heavy_id = {} + index = 0 + for atom in top.atoms(): + if atom.residue.name in AUGC and atom.element.name != "hydrogen": + graph.add_node(index, name=atom.name, resid=atom.residue.index) + full_id2rna_heavy_id[atom.index] = index + rna_indices.append(atom.index) + index += 1 + + for bond in top.bonds(): + if ( + bond.atom1.residue.name in AUGC + and bond.atom1.element.name != "hydrogen" + and bond.atom2.residue.name in AUGC + and bond.atom2.element.name != "hydrogen" + ): + graph.add_edge( + full_id2rna_heavy_id[bond.atom1.index], full_id2rna_heavy_id[bond.atom2.index] + ) + + return graph, rna_indices + + +def gen_dihedral_indices(graph): + dihedral_indices = [] + for n1, n2 in nx.edge_dfs(graph): + neibor1 = list(graph.adj[n1]) + neibor1.remove(n2) + # we want to exclude neightbor n2 + # because we are trying to construct a dihedral pre-n1-n2-post, + # and we don't want pre is the same as n2 + neibor2 = list(graph.adj[n2]) + neibor2.remove(n1) + for pre, post in product(neibor1, neibor2): + if pre != post: # exlucde ring dihedrals + dihedral_indices.append([pre, n1, n2, post]) + return dihedral_indices + + +def calc_include_ij(graph, dihedral_indices): + """ + calculate all the atom pairs that are connected by 1 or 2 or 3 bonds. + Those pairs can't be included in the native contacts + return a matrix for all the atoms in the graph: + True means we can include the pair + False means we should exclude the pair + """ + n_atoms = len(graph) + exclude_ij = np.full((n_atoms, n_atoms), False) + for dihedral in dihedral_indices: + indices = np.array(np.meshgrid(dihedral, dihedral)).T.reshape(-1, 2) + i, j = indices[:, 0], indices[:, 1] + exclude_ij[i, j] = True + + return ~exclude_ij + + step_size = 2 * unit.femtosecond nsteps = 100 @@ -18,17 +84,17 @@ pdb = app.PDBFile("../../inputs/GAGA.box_0mM.pdb") positions = pdb.getPositions(asNumpy=True).value_in_unit(unit.nanometer) -rna_indices = [] -for i, residue in enumerate(pdb.topology.residues()): - if residue.name in ["A", "U", "G", "C"]: - for atom in residue.atoms(): - if atom.element.name != "hydrogen": - rna_indices.append(atom.index) + +rna_graph, rna_indices = init_graph(pdb.topology) +dihedral_indices = gen_dihedral_indices(rna_graph) +include_ij = calc_include_ij(rna_graph, dihedral_indices) rna_pos = positions.astype("float")[np.asarray(rna_indices)] contact_matrix = sd.squareform(sd.pdist(rna_pos)) < contact_cutoff contacts = np.transpose(np.nonzero(contact_matrix)) -rna_id_contacts = np.array([[rna_indices[i], rna_indices[j]] for i, j in contacts if i != j]) +rna_id_contacts = np.array( + [[rna_indices[i], rna_indices[j]] for i, j in contacts if i != j and include_ij[i, j]] +) # notice that we need to get rid of self-self contact! indices = np.unique(rna_id_contacts) references = positions.astype("float")[np.asarray(indices)] diff --git a/examples/openmm/cv_Q/contact_pairs.npy b/examples/openmm/cv_Q/contact_pairs.npy index f953d2fa..eaed2b7a 100644 Binary files a/examples/openmm/cv_Q/contact_pairs.npy and b/examples/openmm/cv_Q/contact_pairs.npy differ diff --git a/examples/openmm/cv_Q/contact_pairs_remapped.npy b/examples/openmm/cv_Q/contact_pairs_remapped.npy index 24cd2000..3be320a1 100644 Binary files a/examples/openmm/cv_Q/contact_pairs_remapped.npy and b/examples/openmm/cv_Q/contact_pairs_remapped.npy differ diff --git a/examples/openmm/cv_Q/log b/examples/openmm/cv_Q/log index cae3205f..420e8246 100644 --- a/examples/openmm/cv_Q/log +++ b/examples/openmm/cv_Q/log @@ -1,101 +1,101 @@ #"Step","Time (ps)","Speed (ns/day)","Elapsed Time (s)","Time Remaining" -1,0.002,0,0.00022292137145996094,-- -2,0.004,0.245,0.7055084705352783,1:09 -3,0.006,0.249,1.3880844116210938,1:07 -4,0.008,0.35,1.4817473888397217,0:47 -5,0.01,0.439,1.5761082172393799,0:37 -6,0.012,0.517,1.6709473133087158,0:31 -7,0.014,0.587,1.765103816986084,0:27 -8,0.016,0.632,1.9146194458007812,0:25 -9,0.018000000000000002,0.688,2.0096611976623535,0:22 -10,0.020000000000000004,0.739,2.1046254634857178,0:21 -11,0.022000000000000006,0.784,2.20461368560791,0:19 -12,0.024000000000000007,0.827,2.2998063564300537,0:18 -13,0.02600000000000001,0.866,2.3941025733947754,0:17 -14,0.02800000000000001,0.903,2.488455295562744,0:16 -15,0.030000000000000013,0.936,2.583986282348633,0:15 -16,0.032000000000000015,0.968,2.6787307262420654,0:15 -17,0.034000000000000016,0.977,2.828676462173462,0:14 -18,0.03600000000000002,1,2.925293207168579,0:14 -19,0.03800000000000002,1.03,3.021669626235962,0:13 -20,0.04000000000000002,1.05,3.116676092147827,0:13 -21,0.04200000000000002,1.08,3.211094856262207,0:12 -22,0.044000000000000025,1.1,3.3069169521331787,0:12 -23,0.04600000000000003,1.12,3.401900053024292,0:11 -24,0.04800000000000003,1.14,3.497349500656128,0:11 -25,0.05000000000000003,1.15,3.5924248695373535,0:11 -26,0.05200000000000003,1.17,3.6875064373016357,0:10 -27,0.054000000000000034,1.19,3.782982349395752,0:10 -28,0.056000000000000036,1.19,3.934417724609375,0:10 -29,0.05800000000000004,1.2,4.0313780307769775,0:10 -30,0.06000000000000004,1.21,4.127398252487183,0:09 -31,0.06200000000000004,1.23,4.22234582901001,0:09 -32,0.06400000000000004,1.24,4.3169190883636475,0:09 -33,0.06600000000000004,1.25,4.410840749740601,0:09 -34,0.06800000000000005,1.27,4.506067752838135,0:09 -35,0.07000000000000005,1.28,4.600708961486816,0:08 -36,0.07200000000000005,1.29,4.695460081100464,0:08 -37,0.07400000000000005,1.3,4.789996862411499,0:08 -38,0.07600000000000005,1.29,4.942267894744873,0:08 -39,0.07800000000000006,1.3,5.039082050323486,0:08 -40,0.08000000000000006,1.31,5.1348817348480225,0:07 -41,0.08200000000000006,1.32,5.230442762374878,0:07 -42,0.08400000000000006,1.33,5.325843095779419,0:07 -43,0.08600000000000006,1.34,5.420317649841309,0:07 -44,0.08800000000000006,1.35,5.515973329544067,0:07 -45,0.09000000000000007,1.35,5.612504482269287,0:07 -46,0.09200000000000007,1.36,5.708605527877808,0:06 -47,0.09400000000000007,1.37,5.803056955337524,0:06 -48,0.09600000000000007,1.36,5.95892071723938,0:06 -49,0.09800000000000007,1.37,6.054737567901611,0:06 -50,0.10000000000000007,1.38,6.149785280227661,0:06 -51,0.10200000000000008,1.38,6.244807243347168,0:06 -52,0.10400000000000008,1.39,6.340013265609741,0:05 -53,0.10600000000000008,1.4,6.434708833694458,0:05 -54,0.10800000000000008,1.4,6.529672622680664,0:05 -55,0.11000000000000008,1.41,6.624948740005493,0:05 -56,0.11200000000000009,1.41,6.720944166183472,0:05 -57,0.11400000000000009,1.42,6.816545009613037,0:05 -58,0.11600000000000009,1.41,6.974427938461304,0:05 -59,0.11800000000000009,1.42,7.071075916290283,0:04 -60,0.12000000000000009,1.42,7.167316675186157,0:04 -61,0.1220000000000001,1.43,7.262823104858398,0:04 -62,0.1240000000000001,1.43,7.357532739639282,0:04 -63,0.12600000000000008,1.44,7.452322006225586,0:04 -64,0.12800000000000009,1.44,7.54635763168335,0:04 -65,0.1300000000000001,1.45,7.640897512435913,0:04 -66,0.1320000000000001,1.45,7.735429525375366,0:04 -67,0.1340000000000001,1.46,7.830808162689209,0:03 -68,0.1360000000000001,1.45,7.9847142696380615,0:03 -69,0.1380000000000001,1.45,8.080800533294678,0:03 -70,0.1400000000000001,1.46,8.176419496536255,0:03 -71,0.1420000000000001,1.46,8.271305561065674,0:03 -72,0.1440000000000001,1.47,8.366270303726196,0:03 -73,0.1460000000000001,1.47,8.462299823760986,0:03 -74,0.1480000000000001,1.47,8.558221817016602,0:03 -75,0.1500000000000001,1.48,8.65354609489441,0:02 -76,0.1520000000000001,1.48,8.749397039413452,0:02 -77,0.1540000000000001,1.48,8.845061540603638,0:02 -78,0.1560000000000001,1.49,8.940609455108643,0:02 -79,0.1580000000000001,1.48,9.093991756439209,0:02 -80,0.16000000000000011,1.49,9.19190239906311,0:02 -81,0.16200000000000012,1.49,9.287736892700195,0:02 -82,0.16400000000000012,1.49,9.383140087127686,0:02 -83,0.16600000000000012,1.49,9.478439092636108,0:01 -84,0.16800000000000012,1.5,9.572929620742798,0:01 -85,0.17000000000000012,1.5,9.667733669281006,0:01 -86,0.17200000000000013,1.5,9.762798309326172,0:01 -87,0.17400000000000013,1.51,9.857316732406616,0:01 -88,0.17600000000000013,1.5,10.010178804397583,0:01 -89,0.17800000000000013,1.5,10.105587244033813,0:01 -90,0.18000000000000013,1.51,10.19968056678772,0:01 -91,0.18200000000000013,1.51,10.294670820236206,0:01 -92,0.18400000000000014,1.51,10.389363050460815,0:00 -93,0.18600000000000014,1.52,10.48458743095398,0:00 -94,0.18800000000000014,1.52,10.580000162124634,0:00 -95,0.19000000000000014,1.52,10.675233364105225,0:00 -96,0.19200000000000014,1.52,10.771239757537842,0:00 -97,0.19400000000000014,1.53,10.866724252700806,0:00 -98,0.19600000000000015,1.52,11.023080110549927,0:00 -99,0.19800000000000015,1.52,11.120671272277832,0:00 -100,0.20000000000000015,1.53,11.216019868850708,0:00 +1,0.002,0,0.00018405914306640625,-- +2,0.004,0.14,1.2383630275726318,2:01 +3,0.006,0.141,2.4549551010131836,1:59 +4,0.008,0.204,2.5401570796966553,1:21 +5,0.01,0.258,2.6760189533233643,1:03 +6,0.012,0.313,2.764267683029175,0:51 +7,0.014,0.364,2.850210666656494,0:44 +8,0.016,0.412,2.936511993408203,0:38 +9,0.018000000000000002,0.458,3.021488666534424,0:34 +10,0.020000000000000004,0.501,3.1063854694366455,0:31 +11,0.022000000000000006,0.537,3.2184159755706787,0:28 +12,0.024000000000000007,0.575,3.3063764572143555,0:26 +13,0.02600000000000001,0.611,3.391796112060547,0:24 +14,0.02800000000000001,0.637,3.528839349746704,0:23 +15,0.030000000000000013,0.669,3.618229866027832,0:21 +16,0.032000000000000015,0.7,3.70469069480896,0:20 +17,0.034000000000000016,0.73,3.7896921634674072,0:19 +18,0.03600000000000002,0.758,3.8752076625823975,0:18 +19,0.03800000000000002,0.785,3.959963798522949,0:17 +20,0.04000000000000002,0.812,4.044415712356567,0:17 +21,0.04200000000000002,0.837,4.129042148590088,0:16 +22,0.044000000000000025,0.861,4.213843584060669,0:15 +23,0.04600000000000003,0.884,4.299461603164673,0:15 +24,0.04800000000000003,0.896,4.43713116645813,0:14 +25,0.05000000000000003,0.916,4.525743722915649,0:14 +26,0.05200000000000003,0.937,4.611554384231567,0:13 +27,0.054000000000000034,0.957,4.696833610534668,0:13 +28,0.056000000000000036,0.976,4.782029390335083,0:12 +29,0.05800000000000004,0.994,4.866597890853882,0:12 +30,0.06000000000000004,1.01,4.952362775802612,0:11 +31,0.06200000000000004,1.03,5.038139581680298,0:11 +32,0.06400000000000004,1.05,5.1235737800598145,0:11 +33,0.06600000000000004,1.06,5.209028005599976,0:10 +34,0.06800000000000005,1.08,5.294060468673706,0:10 +35,0.07000000000000005,1.09,5.378709316253662,0:10 +36,0.07200000000000005,1.1,5.5171427726745605,0:10 +37,0.07400000000000005,1.11,5.607499361038208,0:09 +38,0.07600000000000005,1.12,5.6947526931762695,0:09 +39,0.07800000000000006,1.14,5.780440092086792,0:09 +40,0.08000000000000006,1.15,5.865738153457642,0:09 +41,0.08200000000000006,1.16,5.9505109786987305,0:08 +42,0.08400000000000006,1.17,6.036832332611084,0:08 +43,0.08600000000000006,1.19,6.122369050979614,0:08 +44,0.08800000000000006,1.2,6.208016395568848,0:08 +45,0.09000000000000007,1.2,6.345700979232788,0:07 +46,0.09200000000000007,1.21,6.436450719833374,0:07 +47,0.09400000000000007,1.22,6.523990631103516,0:07 +48,0.09600000000000007,1.23,6.610021352767944,0:07 +49,0.09800000000000007,1.24,6.696079730987549,0:07 +50,0.10000000000000007,1.25,6.78236722946167,0:06 +51,0.10200000000000008,1.26,6.86974573135376,0:06 +52,0.10400000000000008,1.27,6.956451654434204,0:06 +53,0.10600000000000008,1.28,7.042772054672241,0:06 +54,0.10800000000000008,1.28,7.128982305526733,0:06 +55,0.11000000000000008,1.29,7.215255260467529,0:06 +56,0.11200000000000009,1.29,7.354937314987183,0:05 +57,0.11400000000000009,1.3,7.445296287536621,0:05 +58,0.11600000000000009,1.31,7.532633543014526,0:05 +59,0.11800000000000009,1.32,7.619118928909302,0:05 +60,0.12000000000000009,1.32,7.706413745880127,0:05 +61,0.1220000000000001,1.33,7.792268514633179,0:05 +62,0.1240000000000001,1.34,7.878641843795776,0:04 +63,0.12600000000000008,1.35,7.964428663253784,0:04 +64,0.12800000000000009,1.35,8.050028324127197,0:04 +65,0.1300000000000001,1.36,8.135518312454224,0:04 +66,0.1320000000000001,1.37,8.221822261810303,0:04 +67,0.1340000000000001,1.37,8.30713152885437,0:04 +68,0.1360000000000001,1.37,8.44601321220398,0:04 +69,0.1380000000000001,1.38,8.535516738891602,0:03 +70,0.1400000000000001,1.38,8.623224973678589,0:03 +71,0.1420000000000001,1.39,8.710314750671387,0:03 +72,0.1440000000000001,1.39,8.796385526657104,0:03 +73,0.1460000000000001,1.4,8.882198572158813,0:03 +74,0.1480000000000001,1.41,8.968223094940186,0:03 +75,0.1500000000000001,1.41,9.054269313812256,0:03 +76,0.1520000000000001,1.42,9.139775276184082,0:02 +77,0.1540000000000001,1.42,9.279144048690796,0:02 +78,0.1560000000000001,1.42,9.368648052215576,0:02 +79,0.1580000000000001,1.43,9.45512843132019,0:02 +80,0.16000000000000011,1.43,9.54301142692566,0:02 +81,0.16200000000000012,1.44,9.629714012145996,0:02 +82,0.16400000000000012,1.44,9.717389106750488,0:02 +83,0.16600000000000012,1.45,9.803569078445435,0:02 +84,0.16800000000000012,1.45,9.88936734199524,0:01 +85,0.17000000000000012,1.46,9.975153923034668,0:01 +86,0.17200000000000013,1.46,10.061434745788574,0:01 +87,0.17400000000000013,1.46,10.202369928359985,0:01 +88,0.17600000000000013,1.46,10.292763471603394,0:01 +89,0.17800000000000013,1.47,10.379234552383423,0:01 +90,0.18000000000000013,1.47,10.464409589767456,0:01 +91,0.18200000000000013,1.47,10.549905061721802,0:01 +92,0.18400000000000014,1.48,10.635522603988647,0:00 +93,0.18600000000000014,1.48,10.721351146697998,0:00 +94,0.18800000000000014,1.49,10.807268857955933,0:00 +95,0.19000000000000014,1.49,10.89450216293335,0:00 +96,0.19200000000000014,1.5,10.980336904525757,0:00 +97,0.19400000000000014,1.5,11.065917491912842,0:00 +98,0.19600000000000015,1.5,11.152062177658081,0:00 +99,0.19800000000000015,1.5,11.291614294052124,0:00 +100,0.20000000000000015,1.5,11.381214141845703,0:00 diff --git a/examples/openmm/cv_Q/state.pkl b/examples/openmm/cv_Q/state.pkl index 93388932..cbd7e4ab 100644 Binary files a/examples/openmm/cv_Q/state.pkl and b/examples/openmm/cv_Q/state.pkl differ