BUG: branched dextran no longer possible to produce directly

[tylerjereddy]

Working on latest master branch of polyply.

I was trying to roughly follow the tutorial at https://github.com/marrink-lab/polyply_regression_tests/blob/main/examples/liquid_liquid_phase_separation/PolyplyTutorial.ipynb

Then, I noticed that dextran.martini3.ff doesn't even have alpha 1->3 linkages anymore, they were deleted for some reason in 0a76a13 -- was that an accident? Perhaps a branched build regression test is needed there? It looks like you already have an unbranched dextran regression test at polyply/tests/test_data/library_tests/martini3/DEX/polyply/command.

PEO polymer construction seems to work just fine for me locally, and so does unbranched dextran. For reproducer, this could be pared down quite a lot, but you'll get the idea--need to generate a branched JSON file with DEX 1-3 linkages/branches, then the failure to gen_params can be reproduced. I can help you produce a smaller regression test if you really want, but fixing the parameters in a MARTINI-3 compliant way is probably beyond the scope of my expertise, though you may just need to restore some deleted lines to fix it?

# modified from: https://github.com/marrink-lab/polyply_regression_tests/blob/main/examples/liquid_liquid_phase_separation/PolyplyTutorial.ipynb

import sys
import random
import json
import numpy as np
import networkx as nx
from networkx.readwrite import json_graph
from plot_polymers import plot_graph_nicley
import matplotlib.pyplot as plt

def write_graph_to_json(graph, filename):
    g_json = json_graph.node_link_data(graph)
    with open(filename, "w") as file_handle:
        json.dump(g_json, file_handle, indent=2)

def random_linear_comb(n_mon_tot, max_mon_side, p_side_chain):
    """
    Generate a graph that consists of a backbone and branches coming off
    the backbone such that the total number of monomers is `n_mon_tot`.
   
    Parameters
    ----------
    n_mon_tot: int
        total number of monomers in the polymer
    max_mon_side:
        the maximum number of monomers per branched off arm
    p_side_chain: list[float]
        probabilities to add at a level corrsponding to the
        position in the list i.e. probability of adding a monomer
        to a side chain beyond the back-bone is list[0], adding a monomer branched
        off 2 level from the back-bone is list[1].
        
    Returns
    --------
    nx.Graph
    """
    G = nx.Graph()
    # add the first edge to get started
    G.add_edge(0, 1, **{"linktype": "a16"})
    # backbone-nodes
    bb_nodes = [0, 1]
    ndx = 1
    prev_node = "BB"
    while ndx < n_mon_tot:
        # then we make the choice back-bone or side-chain
        # but only if we did not add a side-chain before
        # because we build a comb which has only 1 side chain
        # per back-bone otherwise it's a tree
        if prev_node == "BB":
            choice = random.choices(["BB", "SC"], weights=[1-p_side_chain[0], p_side_chain[0]])[0]
        else:
            choice = "BB"

        # if choice is BB we append to the list of BB nodes and continue
        if choice == "BB":
            G.add_edge(bb_nodes[-1], ndx+1, **{"linktype": "a16"})
            ndx += 1
            bb_nodes.append(ndx)
            prev_node = "BB"

        # if choice is SC we start growing a arm
        else:
            # add the first node of the arm
            G.add_edge(bb_nodes[-1], ndx+1, **{"linktype": "a13"})
            ndx += 1
            prev_node = "SC"
            # an arm can have at most max_side_chains monomers beyond first level
            # we add them until we draw a stop
            for sc_idxs in range(1, max_mon_side):
                choice = random.choices(["stop", "add"], weights=[
                                        1-p_side_chain[sc_idxs], p_side_chain[sc_idxs]])[0]
                if choice == "stop":
                    break
                else:
                    G.add_edge(ndx, ndx+1, **{"linktype": "a16"})
                    ndx += 1
                prev_node = "SC"
                # if we exceed the max number of monomers we terminate
                if ndx >= n_mon_tot:
                    return G

        # if we exceed the max number of monomers we terminate
        if ndx >= n_mon_tot:
            return G

if __name__ == "__main__":
    graph = random_linear_comb(61, 2, [0.5, 0.95, 0.95, 0.1])
    nx.set_node_attributes(graph, "DEX", "resname")
    resids = {node: node+1 for node in graph.nodes}
    nx.set_node_attributes(graph, resids, "resid")
    fig = plt.figure()
    ax_dextran = fig.add_subplot(111)
    plot_graph_nicley(graph, link_labels=True, ax=ax_dextran)
    ax_dextran.set_title("polyply branched dextran (10 kg/mol) prototype residue graph")
    fig.set_size_inches(20, 8)
    fig.savefig("dextran_residue_graph.png", dpi=300)
    write_graph_to_json(graph, 'dextran_polymer_1.json')

polyply gen_params -lib martini3 -seqf dextran_polymer_1.json -name DEX_polymer_1 -o dex_polymer_1.itp

<snip>
WARNING - general - Missing a link between residue 56 DEX and residue 57 DEX.
WARNING - general - Missing a link between residue 57 DEX and residue 58 DEX.
WARNING - general - Missing a link between residue 58 DEX and residue 59 DEX.
WARNING - general - Missing a link between residue 58 DEX and residue 61 DEX.
WARNING - general - Missing a link between residue 59 DEX and residue 60 DEX.
WARNING - general - Missing a link between residue 61 DEX and residue 62 DEX.
WARNING - general - Missing a link between residue 62 DEX and residue 63 DEX.
<snip>

[fgrunewald]

Dear @tylerjereddy,

I apologize for the delayed answer, however, I was at a conference last week. In a nutshell, the parameters for branched Dextran that are in the input files to be used with the jupyter-notebook are a preliminary Dextran model, which minimally diverges from the one published in JCTC. In order to be safe the branched links were therefore removed from the current Martini 3 repository (until someone complains). Not an ideal solution but we're really short on people.

Restoring the branched links should be as straightforward as adding the following directive to the existing links and parametrizing the dihedral angles, which in the previous versions were not present.

[ link ]
; alpha 1->3
resname "DEX"
[ atoms ]
A  { }
B  {"replace": {"atype": "SP1r"}}
C  { }
>A {"replace": {"atype": "SN6"}}
>B { }
>C { }

[ bonds ]
; i  j  funct  length  force constant
  B  >A   1      0.32  5500 {"group": "a13",  "edge": "false"}

[ angles ]
; i  j  k funct angle force con.
C   B   >A   1  180 180 {"edge": "flase", "group": "CBA - a13"}
B  >A   >B   1  115 180 {"edge": "flase", "group": "BAB - a13"}
A   B   >A   1  108 180 {"edge": "flase", "group": "ABA - a13"}
B  >A   >C   1  85  162 {"edge": "flase", "group": "BAC - a13"}

[ dihedrals ]
A   B   >A   >B   1 xxx

[ edges ]
B  >A {"linktype": "a13"}

[ link ]
resname "DEX"
[ dihedrals ]
B +A  +B ++A xxx {"version": "1", "group": "BABA - a13"}

If you want I can assist you in parametrizing them. I have the all-atom data and reference distributions. It just needs someone with a little time to optimize them. Let me know!