X2top creates non symmetric dihedrals for benzene (some C have 6 dihedrals, others 0)

Dear Gromacs Team,

first thanks for this great SW!

I am using:

GROMACS version: 2021.2
GROMACS modification: Yes (it says so, but this seems to be a bug in the WIndows build - I built from sources without modifications using CMake and VisualStudio).

I use x2top to convert Fullerene like structures to .top files using OPLS-AA/L with a modified atomname2type.n2t file but unmodified force field files.

An interesting effect is that for benzene x2top creates 6 dihedrals, all of which have the first C atom in common at position aj.

[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5
    1     2     4     3     3
    1     2     4     8     3
   12     2     4     3     3
   12     2     4     8     3
    1     2    12     6     3
    1     2    12    11     3

The even atoms are the C in order 2 4 8 10 6 12. So while C atom 2 appears in 6 dihedrals, C atom 10 doesn’t appear at all. One would expect that the top file for benzene is symmetric, so this is obviosuly wrong.

I wonder how many dihedrals benzene should have. Should there be 24 in total? 12x H-C-C-C (starting at each H going left and right), 6x H-C-C-H and 6xC-C-C-C? The corresponding dihedral types are available in OPLS-AA/L. In short should every bonded chain of 4 atoms have a dihedral? Or are there rules like one shouldn’t have more than a certain number of dihedrals per bond? E.g. how many dihedrals should ethan have? Only one around the C-C bound or 9 for each combination of H-C-C-H paths (3 H on each end make 3x3=9). I would expect that pencil and paper one uses just one dihedral, but one could certainly define a force field such that there should be 9 and the energies are divided accordingly. Is this a choice force fields can take? Please bear with me if I am silly - I am physicist / computer scientist.

Btw.: looking at x2top with computer scientist eyes I wonder if it wouldn’t be a better idea to create the top files for arbitrary structures using a SAT/SMT solver like Z3 rather than with the heuristics x2top uses. A SAT solver shouldn’t have difficulties to find an assignment of atom types such that for the complete structure the bond, angle and dihedral information exist in the force field, provides such an assignment exists. This should be pretty fast for up to a few 1000 atoms. As far as I can tell it is so that there are structures for which one cannot write an atomnames2types.n2t file so that x2top would select the right atom types - even if one would write the atomnames2types.n2t especially for this structure. E.g. in OPLS-AA/L one would choose a different type of nitrogen in case the second next atom is also nitrogen, but x2top only looks at direct neighbors (as far as I understand). One wouldn’t have this issue with a SAT solver, but then I don’t know if it is sufficient to have a complete set of bonds, angles and dihedrals to be chemically somehow reasonable (I understand that one shouldn’t use x2top if one wants to be chemically reasonable …). With a SMT solver one could add some weighting and minimize some cost. Things might get a bit slow then, though. A downside of a SAT solver approach would be that one can’t produce reasonable error messages - it would just say “constraints are unsatisfiable”. I am not sure, but I think in general the problem x2top tries to solve has a level of complexity which suggest one leverages the 1000s of man years people have put into SAT/SMT solvers to solve it. I just wonder if this has been tried. If you think it makes sense from a chemical point of view, I could do some quick prototype.

Best regards,

Michael

I meanwhile figured out that the answer likely is that bonds in rings shouldn’t have dihedrals (at least no proper ones) and x2top simpyl doesn’t detect rings. So I wrote a python script to do the work of x2top which does have a detection for closed rings.

I guess x2top assumes that each bond may at most have one dihedral and gets highly confused if it detects more than one possible dihedral for a bond, which I guess is a quite common situation. In my python script I count the dihedrals and scale the energy down accordingly. Usually all dihedrals are the same, so that the effect of this is an averaging of the fluctuations of atoms around the bond.

I have written a script, find it at GitHub - masrul/GenTopo: Gromacs topology template generator . Which detects rings and finds improper in addition to proper dihedrals.

1 Like

Thanks, nice! What I don’t understand is what force field it is intended for.

Did you think about putting this back into x2top? I regard my python code more as a prototype to experiment and to see what makes sense, but I think the goal should be to improve Gromacs.

It is not intended for any particular force-field. Main intention was to setup topology for non-standard/new force field, where user knows atom types and charges for each atom. All it does, create a template.

1 Like

Ah I see! This is indeed quite useful as a starting point when experimenting with new force fields.