Convert itp files to ff files #327
Conversation
Co-authored-by: Peter C Kroon <pckroon@users.noreply.github.com>
|
@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. |
There was a problem hiding this comment.
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.
| parser_itp_ff.add_argument('-o', dest="outpath", type=Path) | ||
| parser_itp_ff.add_argument('-c', dest="charges", type=float, nargs='*') | ||
| parser_itp_ff.add_argument('-c', dest="res_charges", nargs='+', type=lambda s: s.split(':'),) |
There was a problem hiding this comment.
Needs a help with the format/syntax
| parser_itp_ff.add_argument('-f', dest='inpath', type=Path, required=False, default=[], | ||
| help='Input file (ITP|FF)', nargs='*') |
There was a problem hiding this comment.
https://docs.python.org/3/library/argparse.html#argparse.FileType
Maybe also in other places. Optional though
| self.resid += 1 | ||
|
|
||
| def label_fragments_from_graph(self, fragment_graphs): | ||
| def extract_unique_fragments(self, reference_graph): |
There was a problem hiding this comment.
I guess the docstring still needs updating
| # finally we simply collect one graph per restype | ||
| # which are the most centrail (i.e. avoid ends) |
| # which are the most centrail (i.e. avoid ends) | ||
| unique_fragments = {} | ||
| frag_centrality = {} | ||
| centrality = nx.betweenness_centrality(self.res_graph) |
There was a problem hiding this comment.
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.
| 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']: |
| if target_atoms == ref_inter.atoms and\ | ||
| target_inter.parameters != ref_inter.parameters: |
There was a problem hiding this comment.
| 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: |
| 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()} |
There was a problem hiding this comment.
| #link_to_mol_atoms = {value:key for key, value in mol_atoms_to_link_atoms.items()} |
| link_inter = Interaction(atoms=link_atoms, | ||
| parameters=target_inter.parameters, | ||
| meta={}) |
There was a problem hiding this comment.
| link_inter = Interaction(atoms=link_atoms, | |
| parameters=target_inter.parameters, | |
| meta={}) | |
| link_inter = Interaction(atoms=link_atoms, | |
| parameters=target_inter.parameters, | |
| meta={}) |
| molecule, | ||
| min(resids)) | ||
| link_atoms = mol_to_link.values() | ||
| link = vermouth.molecule.Link() |
There was a problem hiding this comment.
| link = vermouth.molecule.Link() |
| # a little dangerous but mostly ok; if there are no changes to | ||
| # the atoms we can continue | ||
| if len(replace_dict) == 0: | ||
| continue |
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:
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:
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
Known Issues