def contacts by filtering 1-4 bonding pairs

update the examples, this is simply a better choice of
the contacts. Otherwise it contains bonded pairs which
will always be there.

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,24 +14,87 @@
 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
 
 contact_cutoff = 0.5  # nanometer
 
 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)]

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