How to contribute (coarse-grained) force field parameters with tabulated potentials?

[ricalessandri]

Just thought I would open first a discussion about this, since it might be useful for others. How to contribute coarse-grained force field parameters with tabulated potentials (I developed a couple of models in https://arxiv.org/abs/2209.02072)? Basically, what one needs are files that contain the potentials of all the pair interactions. Should they be uploaded/stored in the polyply data somewhere or uploaded on Zenodo/Figshare and then have a link to it in the .ff header?

[pckroon]

Good question. We should probably separate out bonded and non-bonded tab potentials, since vermouth doesn't deal with non-bonded parameters at all at the moment.
I think it would make sense to enable support for tabulated bonded potentials in vermouth, and I guess by extension non-bonded tab pot in polyply? This will require some (re)designing of file formats I'm afraid.

[fgrunewald]

We're tlaking gromacs herr? Are the interactions split in bonded and nonbonded contributions actually? Would it be possible for the bonded interactions to define them as types? If yes we can just offload them to the ff-files.

[ricalessandri]

Yes, I am talking about Gromacs. I am not familiar with how the bonded tabulated potentials are defined, actually, because I have parametrized those in the standard way - i.e., fitting standard harmonic, etc. potentials.

This is the stuff you need for the nonbonded:

1. The following lines in the mdp (with the bead types, in this case 5):

   energygrps               = VNL EST C1 C2 N4
   

   and (with all the pair interactions, in this case 15 pair interactions (5*5+1)/2):

   energygrp_table          = VNL VNL VNL EST VNL C1 VNL C2 VNL N4 EST EST EST C1 EST C2 EST N4 C1 C1 C1 C2 C1 N4 C2 C2 C2 N4 N4 N4
   

2. An index.ndx that contains which bead belongs to which bead type (somehow this is needed for how GMX works; this will of course vary depending on how many molecules one is simulating);

   [ VNL ]
   1    7   13   19 ...
   ...
   [ EST ]
   2    8   14   20  ...
   ...
   

3. Files called like table_VNL_C1.xvg for all the pair interactions and that contain the actual pair potential.

The .ff then actually looks like a "normal" .ff file; just the bead/atom types need to have the same names as the ones listed in the mdp and for which table_XX_YY.xvg files exist. Hence, at least for the nonbonded, probably integrating these in the .ff files is not that easy and probably does not make much sense... but maybe a description such as the above (with links to the files) can be included in a "Description/Info/Warning" section (that then gets also printed by polyply in the generated itp file)?

Let me know what your thoughts on this.

[pckroon]

IIRC (but it's been a while) it works similar for bonded potentials: you designate them the appropriate function type (forgot the number. 8?), and the parameter is an integer. That integer is then used to lookup the correct xvg table with energies, e.g. table_d0.xvg. I'm almost definitely wrong on the details, but general picture should be correct.

Do we want to limit this discussion to one of the two? There currently is a need for non-bonded tabulated potentials, so maybe just look at that?

I think it would make sense to keep the interaction tables (table_C1_VNL.xvg) as is in the data folder. The main issue will be detecting when non-bonded tables are required/wanted/.... I guess that's something that needs to be added to the .ff file format. After that writing the ndx file and mdp lines should be easy. For easy-of-use we can also copy the relevant xvg file(s) to the current working directory (or at least, next to the top file). All that would be left is deciding on how polyply calculates the non-bonded energy term between these kind of particles during the random walk. It shouldn't be too hard to interpolate the table using scipy?

[fgrunewald]

I think there are a couple of interesting but separate issues/additions/suggestions here:

1 - I like the idea of a warning directive in the ff format. That would be super useful for many polymers in polyply. Expect a PR this month. We can always add a warning independent from how we solve the problem in detail.

2 - Let's focus on the non bonded tabulated potentials first. There are two aspects in this topic as well:

2a: Are we going to provide the table parameters as part of the polyply library? I am against that because we also do not provide force-field parameters. Providing them means having to take care of those which I would rather not do. A warning or usage note would be better at least for the near future in my opinion.

2b: How are we generating coordinates? The random-walk uses a LJ potential independent of the underlying non bonded interactions in the actual force-field. So the tabulated potential only matters for estimating the volume of the residue. We have two options: (1) We can add a CLI argument on which the user provides the path to the tabulated potential files. We can generate the pair interactions automatically from the file names, interpolate the potential, find the zero energy and use that as a volume. Not sure that it will be worth it though. (2) We do the same as we plan to do for DNA. We provide a build file with a volume entry in the library directory. When provided with a library in gen_coords polyply reads the volume from either the build file provided via CLI or from the library. Thus we simply store an approximate volume with the ff-files in the library. That would be my favourite option, because we don't need to add any parsers or code. It already works to provide a volume via the build file of the CLI but the library option currently depends on PR #253

So my plan would be as follows:

• create a warnings directive like citations directive
• provide users with a link to tabulated potential through the warning entry; users find those files together with the force-field on their own
• provide custom volumes for the residues for coordinate generation in build file like for DNA
• patch the top-reader to not look for pair interactions when provided with a -custom_nb CLI flag

I think this should be the most easy option. @ricalessandri @pckroon do we want to divide the tasks among us?

[pckroon]

I'd be happy to implement the first TODO in vermouth (log triggers in the ff file format).

[pckroon]

There's an experimental and untested branch here: https://github.com/marrink-lab/vermouth-martinize/tree/ff_log_entry

[fgrunewald]

The branch looks good to me. So the entry would be:

[debug]
some lines

Perhaps we should rather call it something else than debug. But for now works for me. I'll create a temp polyply branch to see how this would pan out.

[fgrunewald]

@pckroon could we also have a link debug statement. This would be really cool because then you could issue a warning of the user does something not so good. For example, I could create a fake link that patches the debug when the user forgets to add the terminii for instance.

[pckroon]

Of course. My idea was to provide debug, info, warning and error sections for Blocks, Links and Modifications (And I just noticed I messed up, typoes everywhere :) ).
I'm still debating it's worth the effort to also add it to mappings. Probably not? Maybe later?

[fgrunewald]

hahaha yeah I saw them as well after playing with the code. Mappings I think are fine. There is no way to change a mapping anyways is there? For that we have the modifications. You could say if the mapping is strange a modification applies that induces a debug warning