Implementation of the GEN_FF workflow

[fgrunewald]

This is the preliminary implementation of the gen_ff workflow. To generate ff-parameters from an itp file use:

polyply gen_ff -f <your_itp_file> -i <itp/top file> -s <cgsmiles_str> -o output.ff -c <resname>:<charge> 

You can then generate an arbitrary polymer with:

polyply gen_params -f output.ff -seq <cgsmiles_str> -o out.itp

LICENSE NOTICE

This code is currently licensed under PolyForm Noncommercial License 1.0.0, meaning it can be adopted shared and used as long as not for commercial purposes.

[tylerjereddy]

While working on the effort in gh-390 to build SEBS, using this command (via subprocess):

subprocess.CalledProcessError: Command '['polyply', 'gen_ff', '-i', 'smiles_molecule.acpype/smiles_molecule_GMX_OPLS.top', '-s', '{[#styrene]|2[#ethylene]|2[#butylene]|2[#styrene]|2}.{#styrene=[<]CC[>]c1ccccc1,#ethylene=[<]CCCC[>],#butylene=[<]CC[>](CC)}']' returned non-zero exit status 1.

I ended up with the traceback below the fold, affecting this line.

traceback

Traceback (most recent call last):
  File "/home/treddy/python_venvs/py_311_hemd/bin/polyply", line 292, in <module>
    main()
  File "/home/treddy/python_venvs/py_311_hemd/bin/polyply", line 288, in main
    subprogram(**args_dict)
  File "/home/treddy/python_venvs/py_311_hemd/lib/python3.11/site-packages/polyply/src/gen_ff.py", line 90, in gen_ff
    meta_mol = MetaMolecule.from_cgsmiles_str(force_field=force_field,
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/treddy/python_venvs/py_311_hemd/lib/python3.11/site-packages/polyply/src/meta_molecule.py", line 424, in from_cgsmiles_str
    atomname = meta_graph.nodes[node]['graph'].nodes[atom]['atomname']
               ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^
KeyError: 'atomname'

The error was triggered by this CGSmiles string originally: '{[#styrene]|2[#ethylene]|2[#butylene]|2[#styrene]|2}.{#styrene=[<]CC[>]c1ccccc1,#ethylene=[<]CCCC[>],#butylene=[<]CC[>](CC)}'

I can reproduce the error with a much simpler CGSmiles string, "{[#styrene]|2}.{#styrene=[<]CC[>]c1ccccc1}",
and even with: '{[#styrene]}.{#styrene=[<]CC[>]c1ccccc1}'.

Some debug printing shows the problematic atom in the context of the full CGSmiles input to be this one:

***** atom: 49
***** nodes: [45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56]
***** type(nodes): <class 'networkx.classes.reportviews.NodeView'>
***** atom index: {'charge': 0, 'aromatic': False, 'element': 'H', 'fragid': [3], 'fragname': 'ethylene', 'resname': 'ethylene'}

I'll see if I can dig a bit deeper..

[tylerjereddy]

It turns out that molecule.nodes[atom]['atomname'] is already populated in some cases, so the shim below seems to get past that particular issue:

--- a/polyply/src/meta_molecule.py
+++ b/polyply/src/meta_molecule.py
@@ -421,8 +421,12 @@ class MetaMolecule(nx.Graph):
                     molecule.nodes[atom]['resname'] = resname
                     molecule.nodes[atom]['resid'] = node + 1
                     #print(meta_graph.nodes[node]['graph'].nodes[atom])
-                    atomname = meta_graph.nodes[node]['graph'].nodes[atom]['atomname']
-                    molecule.nodes[atom]['atomname'] = atomname
+                    try:
+                        atomname = meta_graph.nodes[node]['graph'].nodes[atom]['atomname']
+                        molecule.nodes[atom]['atomname'] = atomname
+                    except KeyError:
+                        # the atomname tends to already be populated in some cases
+                        pass

But then there's yet another error:

Traceback (most recent call last):
  File "/home/treddy/python_venvs/py_311_hemd/lib/python3.11/site-packages/polyply/src/graph_utils.py", line 225, in find_one_ismags_match
    mapping = next(raw_matches)
              ^^^^^^^^^^^^^^^^^
StopIteration

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/treddy/python_venvs/py_311_hemd/bin/polyply", line 292, in <module>
    main()
  File "/home/treddy/python_venvs/py_311_hemd/bin/polyply", line 288, in main
    subprogram(**args_dict)
  File "/home/treddy/python_venvs/py_311_hemd/lib/python3.11/site-packages/polyply/src/gen_ff.py", line 97, in gen_ff
    unique_fragments, res_graph = FragmentFinder(target_mol).extract_unique_fragments(meta_mol.molecule)
                                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/treddy/python_venvs/py_311_hemd/lib/python3.11/site-packages/polyply/src/fragment_finder.py", line 141, in extract_unique_fragments
    mapping = find_one_ismags_match(self.molecule,
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/treddy/python_venvs/py_311_hemd/lib/python3.11/site-packages/polyply/src/graph_utils.py", line 228, in find_one_ismags_match
    raise IOError("no match_found")
OSError: no match_found

[fgrunewald]

@tylerjereddy I think this might be a dependency error, where the atomnames are different from what polyply expects. We've been changing quite a bit in the CGSmiles package over the past few weeks.I'll try taking a look to see if I can manage to figure it out.

[tylerjereddy]

@fgrunewald @pckroon I can see different errors at different pysmiles hashes alongside this branch, but it would be hard to bisect since it isn't clear which of the errors is really preferable.

Some pysmiles hash examples:

• lastest master branch at time of writing (1071fbf78929a5c5aed): AttributeError: module 'pysmiles.smiles_helper' has no attribute 'mark_aromatic_atoms'
• 24a15aee3557ad59ef2792437b75bf8de5dc2645 (Sept. 16th): AttributeError: module 'pysmiles.smiles_helper' has no attribute 'mark_aromatic_atoms'
• c43ff4455a6987df3bbe22cb02ffab7e2ae5cfa8 (Sept. 11th): KeyError: 'atomname'
• 8ad2702d48df4d0aab87c0f76827a04bd9141120 (Sept. 2nd): KeyError: 'atomname'
• f48710757c7d7c78d69a4389af95734a46525c7a (June 21): ImportError: cannot import name '_annotate_ez_isomers' from 'pysmiles.smiles_helper'

[tylerjereddy]

I also thought about checking the CI logs on this branch to see what got pulled in on the last push, but it has been a few months since the last commit and the logs are gone I think.

[pckroon]

Pysmiles has seen quite some changes in the past weeks. I think the issues you're running into appeared after pckroon/pysmiles#47 got merged. c43ff4455a6987df3bbe22cb02ffab7e2ae5cfa8 is the last commit before that PR. Note that after pckroon/pysmiles#49 gets merged I intend to release pysmiles 2.0, which should hopefully stabilize things quite a bit.

[tylerjereddy]

@pckroon I just tried that pysmiles commit alongside this branch in my workflow and ended up with same traceback

Traceback (most recent call last):
  File "/home/treddy/python_venvs/py_311_hemd/bin/polyply", line 292, in <module>
    main()
  File "/home/treddy/python_venvs/py_311_hemd/bin/polyply", line 288, in main
    subprogram(**args_dict)
  File "/home/treddy/python_venvs/py_311_hemd/lib/python3.11/site-packages/polyply/src/gen_ff.py", line 90, in gen_ff
    meta_mol = MetaMolecule.from_cgsmiles_str(force_field=force_field,
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/treddy/python_venvs/py_311_hemd/lib/python3.11/site-packages/polyply/src/meta_molecule.py", line 424, in from_cgsmiles_str
    atomname = meta_graph.nodes[node]['graph'].nodes[atom]['atomname']
               ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^
KeyError: 'atomname'

[pckroon]

That's probably an issue in cgsmiles. AFAIK that version of pysmiles is what cgsmiles currently understands, but it's all in motion.

[fgrunewald]

@tylerjereddy et al. This has been an issue in cgsmiles and is resolved upon merging CGSmiles PR33