Convert itp files to ff files by fgrunewald · Pull Request #327 · marrink-lab/polyply_1.0

fgrunewald · 2023-06-16T13:21:39Z

The new utility allows users to extract ff-files automatically from itp files by defining fragments as SMILES. Currently the code only works for linear polymers. Let's look at some examples for how to use this Program.

Say we have a itp that describes any number of PEO oligomers then to extract the parameters for PEO we'd run:

polyply itp_to_ff -i <itp_file.ipt> -sm [CH2]O[CH2] -rn PEO -o new.ff -charge 0

Now of course with PEO there is the challenge of how we treat the termini. They can for example be terminated with an OH group. By default the algorithm will group all atoms that cannot be assigned to a fragment as given by the smile into the next connected fragment. Thus it creates a new bigger fragment. In this example, we would obtain the PEO fragment but also PEOter fragment that describes the termini. If they are asymmetric the code produces 2 new residues describing the termini.

Of course this might not be what we want, since we may want to build a modular ff format. In this case we can also specify the OH terminal explicitly as a fragment. For example, like so:

polyply itp_to_ff -i <itp_file.ipt> -sm [OH][CH2] [CH2]O[CH2] [OH][CH2] -rn OH PEO OH -o new.ff -charge 0

Note that we need to provide the terminal definitions on both ends, because the code works in blocks. Using this command we obtain definition for the termini as a separate block.

@ricalessandri

To Do

Remove print and draw statements
sort the atoms in a block definition according to their name
deal with chirality/tacticity
improve the charge equalisation; also round charges so they are at most 2 significant digits
test on charged molecules
make pysmiles optional dep. and raise Error if not installed
How many repeat units are needed???
add a warning that extract links only work with linear molecules atm
add help text to CLI

Known Issues

sometimes the first match picks a repeat unit that is not at the start of the chain. In that case, all other matches are likely messed up. The workaround is to provide the first smile with all the correct hydrogen atoms
if there are too few repeat units to extract all representative interactions between monomers the small fragment will work correctly but the larger ones will be incorrect
only works with all-atom at the moment

…ors are leftover

Co-authored-by: Peter C Kroon <pckroon@users.noreply.github.com>

fgrunewald · 2024-03-07T20:18:47Z

@ricalessandri that's it; except for the improvements on labeling edges the code is done and has all tests required. I ran it on my CHARMM database and will try some OPLS tomorrow. Please give it a try and see if there are any problems.

pckroon

Nice, this can be super useful!
I forsee some issues with the automatic handling of the termini, but it is a hard problem. You make a Link to deal with those, rather than a Modification? What's the reasoning here?

I'll finish this review after the BigSmiles PR is merged, and this branch has been updated.

pckroon · 2024-03-20T11:38:00Z

@@ -237,12 +237,12 @@ def main(): # pylint: disable=too-many-locals,too-many-statements
                                help='Enable debug logging output. Can be given '
                                'multiple times.', default=0)


Needs a help with the format/syntax

pckroon · 2024-03-20T11:39:22Z

@@ -237,12 +237,12 @@ def main(): # pylint: disable=too-many-locals,too-many-statements
                                help='Enable debug logging output. Can be given '


https://docs.python.org/3/library/argparse.html#argparse.FileType

Maybe also in other places. Optional though

pckroon · 2024-03-20T11:40:59Z

-                self.resid += 1
-
-    def label_fragments_from_graph(self, fragment_graphs):
+    def extract_unique_fragments(self, reference_graph):


I guess the docstring still needs updating

pckroon · 2024-03-20T11:42:42Z

+        # finally we simply collect one graph per restype
+        # which are the most centrail (i.e. avoid ends)


Add the "why"

pckroon · 2024-03-20T11:50:00Z

+        # which are the most centrail (i.e. avoid ends)
+        unique_fragments = {}
+        frag_centrality = {}
+        centrality = nx.betweenness_centrality(self.res_graph)


If you want the most central nodes, see also https://networkx.org/documentation/stable/reference/algorithms/generated/networkx.algorithms.distance_measures.center.html#center
You can experiment a little on what gives the best results.

pckroon · 2024-03-20T12:08:51Z

+        for node in target_block.nodes:
+            target_attrs = target_block.nodes[node]
+            ref_attrs = ref_block.nodes[target_attrs['atomname']]
+            for attr in ['atype', 'mass']:


Why limit to these attrs?

pckroon · 2024-03-20T12:09:21Z

+                    if target_atoms == ref_inter.atoms and\
+                    target_inter.parameters != ref_inter.parameters:


Suggested change

if target_atoms == ref_inter.atoms and\

target_inter.parameters != ref_inter.parameters:

if target_atoms == ref_inter.atoms and target_inter.parameters != ref_inter.parameters:

pckroon · 2024-03-20T12:10:32Z

+                         mol_atoms_to_link_atoms, edges, resnames = _extract_edges_from_shortest_path(target_inter.atoms,
+                                                                                                      molecule,
+                                                                                                      min(resids))
+                         #link_to_mol_atoms = {value:key for key, value in mol_atoms_to_link_atoms.items()}


Suggested change

#link_to_mol_atoms = {value:key for key, value in mol_atoms_to_link_atoms.items()}

pckroon · 2024-03-20T12:10:58Z

+                         link_inter = Interaction(atoms=link_atoms,
+                                                  parameters=target_inter.parameters,
+                                                   meta={})


Suggested change

link_inter = Interaction(atoms=link_atoms,

parameters=target_inter.parameters,

meta={})

link_inter = Interaction(atoms=link_atoms,

parameters=target_inter.parameters,

meta={})

pckroon · 2024-03-20T13:05:55Z

+                                                                         molecule,
+                                                                         min(resids))
+        link_atoms = mol_to_link.values()
+        link = vermouth.molecule.Link()


Suggested change

link = vermouth.molecule.Link()

pckroon

To be continued after #358

pckroon · 2024-03-20T13:11:18Z

+        # a little dangerous but mostly ok; if there are no changes to
+        # the atoms we can continue
+        if len(replace_dict) == 0:
+            continue


Dangerous why?

f.grunewald added 17 commits June 13, 2023 16:02

init draft itp to ff

5968f83

imporve graph matching

c577052

fragment finder with prints

7eff22a

add tests for fragment finder

95c4b87

add test for 100% coverage

ae2794c

refactor graph matchin post isomorph check

101d2b7

add check on node naming

6261186

add pysmiles to tests

a8ce5a1

tests for ffoutput

b8dfa7b

use tmp-file for testing ffoutput

b3ea5ac

modify extract block and use in itp_to_ff

79c38fb

update test for generate templates accordingly

77dfe16

add isomorphism naming

214f5f2

properly check if interactions are equal

ef70012

read itp files

2410b0a

draft round robin tests

450ebc4

fix input types

888515b

ricalessandri added the enhancement New feature or request label Jun 24, 2023

f.grunewald added 2 commits June 26, 2023 11:28

add test print

44cc867

Merge branch 'master' into itp_to_ff

0877210

fgrunewald mentioned this pull request Aug 12, 2023

filter molecule for templates #339

Merged

2 tasks

Fabian Gruenewald added 9 commits November 22, 2023 15:47

clean up output

32cd8f8

methods to deal with charges

a8d1bb9

methods to deal with charges

37bad71

methods to deal with charges

7409098

adjust test

362372e

move extract block to molecule utils

4ed2979

small fix

fa32f76

allow for charged residues and make pysmiles optional import

5207801

make mass optional

737b45c

Fabian Gruenewald and others added 26 commits February 29, 2024 15:43

refactor fragment itp_to_ff

87510bb

change input for itp_to_ff to allow bigmsiles

05df2e5

take most central fragment

0ebfa6a

add special links for terminal modifications

a7cd590

type the charges to float in itp to ff

8e0c257

add H and ] as special characters in big smile parser

7cb3b4c

account for explicit hydrogen in the smiles string input

097ec84

test accounting for explicit hydrogen in the smiles string input

514ba1b

read provided ff file and use these blocks instead of making new ones

3e4a737

adjust doc string

d97632d

skip termini mods if none atoms are different

b6acc73

Merge branch 'big_smiles' into itp_to_ff

6095dc6

redo hydrogen based on valency not based on how many bonding descript…

9d9ee89

…ors are leftover

parse force-field in molprocessor, adjust hydrogen reconstruction

c4f1652

fix tests

7a5dd1f

Apply suggestions from code review

2b9e7a9

Co-authored-by: Peter C Kroon <pckroon@users.noreply.github.com>

allow nested branch expansion

b6d891f

test branch expansion

a867329

add comments all over residue expansion functions

b6f5cc0

address comments

f965e1d

allow for ionic bonds with . syntax

0335956

Merge branch 'big_smiles' into itp_to_ff

1f25c32

fix previous issue with link appending

47fef23

update itp_to_ff tests

7f7fe21

update tests for fragment finder

7268663

remove leftover files

15be6a6

pckroon requested changes Mar 20, 2024

View reviewed changes

ricalessandri mentioned this pull request Jul 8, 2024

Add [10.1021/jacsau.4c00276] polymers to the polyply library #378

Merged

3 tasks

		@@ -237,12 +237,12 @@ def main(): # pylint: disable=too-many-locals,too-many-statements
		help='Enable debug logging output. Can be given '
		'multiple times.', default=0)

		# finally we simply collect one graph per restype
		# which are the most centrail (i.e. avoid ends)

		if target_atoms == ref_inter.atoms and\
		target_inter.parameters != ref_inter.parameters:

	if target_atoms == ref_inter.atoms and\
	target_inter.parameters != ref_inter.parameters:
	if target_atoms == ref_inter.atoms and target_inter.parameters != ref_inter.parameters:

Conversation

fgrunewald commented Jun 16, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

fgrunewald commented Mar 7, 2024

Uh oh!

pckroon left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pckroon left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

fgrunewald commented Jun 16, 2023 •

edited

Loading