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.