Regarding the issue of dihedrals parameters in the Martini 3 force field top files for constructing molecules.

[zhangtrzhangtr]

I have encountered some questions while establishing the Martini 3 force field top files for the molecule depicted in the image. I divided the molecule into 13 beads as illustrated. My goal is to fit the dihedral angle parameters between the central benzothiophene structure and the thiophene pentagon above it. Initially, I computed the torsion potential energy curve between the two structures through quantum chemical calculations. Following this, I obtained the molecular mechanics top file for this structure in the OPLS force field from the Lipargen server. I then calculated the molecular mechanics energy at different dihedral angles of the structure and subtracted the molecular mechanics energy from the quantum mechanics results to obtain the torsion potential energy curve. This curve was fitted using the Ryckaert-Bellemans function, which corresponds to dihedral function type 3 in GROMACS software.
For the Martini 3 force field, I attempted to represent the dihedral angle between benzothiophene and thiophene with beads 1-2-3-4, where 2 is the geometric center virtual site of the thiophene structure's bead, and 3 is the virtual site determined by the geometric center of the three beads that make up the benzene ring in the benzothiophene structure, as referenced in your previously published articles. May I directly apply the Ryckaert-Bellemans function fitting parameters, which were derived by subtracting molecular mechanics energy from quantum mechanics calculations, to the Martini 3 force field top files? That is, utilizing dihedral parameters of function type 3. However, in the models you have previously published, I observed that the top files predominantly utilized periodic types of dihedral parameters, specifically function types 1 or 9. I wish to know if it is appropriate to directly employ the Ryckaert-Bellemans function for fitting dihedral angle parameters in the Martini 3 force field top files and whether my methodology is sound.
Your guidance and advice would be greatly appreciated.

[ricalessandri]

1. The RB is a fine choice. Both the RB model and periodic dihedral potential (types 1 or 9 in Gromacs) are suitable for fitting dihedrals. In fact, the RB model offers more flexibility compared to the periodic dihedral function. My selection of Gromacs function types 1 or 9 was arbitrary.

2. Regarding your approach, directly using the RB parameters would be incorrect because the underlying AA nonbonded parameters are not present in M3. I believe the best approach is to directly implement the QM dihedral in M3, essentially using an RB-fit potential that directly matches the QM data instead of the QM-AA difference. Importantly, this also assumes that nonbonded interactions between the thiophene and central unit beads are excluded. To ensure accuracy, you should then conduct an AA simulation with the RB model (QM-AA as you described) and an M3 simulation using the procedure I outlined. Comparing the dihedral distributions from the 'mapped AA' and M3 simulations should yield consistent results.

[zhangtrzhangtr]

Thank you for your detailed guidance and explanation. I have basically understood your point. By directly using the parameters obtained from quantum chemical calculations and comparing them with the results of all-atom molecular dynamics simulations, excluding the non-bonded interactions between the thiophene structure and the central unit bead, I sincerely wish you smooth sailing in your research endeavors. Ha ha ha.

[yuqingxiong96]

I am very excited to see your discussion, I find it very interesting. I have also been fitting the proper dihedral parameters for a molecule, but I haven't found the tutorial. I wonder if you could give me with some help. I would appreciate it very much. May I have further communication with you? My profile has my email address, I hope you can contact me.

[ricalessandri]

Glad you find the discussion interesting. There's no tutorial on this topic unfortunately. I will add this to a tutorial wishlist, when the next workshop/occasion comes up, I'll try to make one. Feel free to post questions in this discussion thread in the meantime.

[yuqingxiong96]

For periodic dihedral potential, I don't know if I understand correctly: 1. The torsional potential energy curves between the structures of interest are obtained by QM calculation.
2. Directly fit the QM dihedral torsional potential energy curve as an object to get the k, n, ϕ value required by gromacs into the itp file
If this is done, the calculated value of QM is absolute energy, which will cause the calculated k to be very large, may I ask which step went wrong?

[zhangtrzhangtr]

I think you might not have subtracted the energy of each structure by the lowest energy point before fitting.

[yuqingxiong96]

Your suggestion is to use relative energy, right? So that the energy of the lowest point is set to 0?

[zhangtrzhangtr]

yes

[yuqingxiong96]

Thank you very much for your reply, which is very helpful to me.