Some questions regarding bead selection in the Martini 3 force field

[zhangtrzhangtr]

In the image, the thiadiazole five-membered ring structure, where a sulfur atom is connected to two nitrogen atoms, forms a conjugated ring structure. How should I assign Martini 3 force field beads to this structure? It seems that I couldn't find a suitable match in the published Martini 3 force field bead mappings.

[ricalessandri]

@zhangtrzhangtr The first guess would be to use thiazole as a reference https://github.com/ricalessandri/Martini3-small-molecules/blob/main/models/itps/opt-mono/THAZ.itp; this would mean assigning TN6a to those 2 beads you highlighted.
Then, look into "An Analysis of the Lipophilicity of Furazan and Furoxan Derivatives using the CLOGP Algorithm" (by Rosella Calvino, Alberto Gasco, and Albert Leo, 1992) as it should contain the logP for 2,1,3-Benzothiadiazole. Map the rest (S --> TC6, the benzo moiety with two TC5 beads), compute the logP for the CG model and compare it to the exp. value to check if the choice of TN6a beads for the -N=C- fragments is good. Hope this helps.

[zhangtrzhangtr]

Okay, I will try to compare by calculating logP. Thank you for your guidance.

[zhangtrzhangtr]

@ricalessandri
I followed your initial suggestions, using TC6 beads to represent sulfur atoms, TN6a beads for the conjugated structure between nitrogen and carbon atoms, and TC5 beads for carbon atoms on the benzene ring. However, the calculated transfer free energy significantly deviates from the experimental results provided in the article, which lists a logP value of -1.98 for 2,1,3-Benzothiadiazole. Using the formula ΔG(O→W)=−RTln⁡P(O→W)=−2.303RTlog⁡P(O→W), the transfer free energy is calculated to be 11.113 kJ/mol. I then followed the thermodynamic integration method described in the tutorial @at http://cgmartini.nl/index.php/tutorials-general-introduction-gmx5/partitioning-techniques to calculate the desolvation energies in octanol and water, respectively. The results are ΔG(O→vacuum)=22.19 kJ/mol and ΔG(W→vacuum)=17.68 kJ/mol, leading to ΔG(S1→S2)=4.51kJ/mol.

I used ten λ points in my calculations. Below are my gro and itp files. Does this discrepancy suggest that the beads I selected are inappropriate? If so, could you suggest possible improvements? Thank you for your assistance.
;;;;;; benthiazole

[moleculetype]
; molname nrexcl
BTZ 1

[atoms]
; id type resnr residu atom cgnr charge
1 TC6 1 BTZ S1 1 0
2 TN6a 1 BTZ N2 2 0
3 TN6a 1 BTZ N3 3 0
4 TC5 1 BTZ C4 4 0
5 TC5 1 BTZ C5 5 0

[ bonds ]
; i j funct length force_k
1 2 1 0.192 25000 ;
1 3 1 0.192 25000 ;

#ifndef FLEXIBLE
[ constraints ]
#endif
; i j funct length
2 3 1 0.198 100000 ;
2 4 1 0.270 100000 ;
2 5 1 0.363 100000 ;
4 5 1 0.296 100000 ;
3 5 1 0.270 100000 ;

[dihedrals]
; improper
; i j k l funct ref.angle force_k
2 3 4 5 2 180 200 ; default for fused bicyclic compounds
1 2 3 5 2 180 200 ; default for fused bicyclic compounds

[exclusions]

1 2 3 4 5
2 3 4 5
3 4 5
4 5

Generated with cgbuilder
5
0MOL B0 1 0.225 0.000 0.000
0MOL B1 2 0.060 -0.099 0.000
0MOL B2 3 0.060 0.099 0.000
0MOL B3 4 -0.206 -0.148 0.000
0MOL B4 5 -0.206 0.148 0.000
10 10 10

[ricalessandri]

Thanks. To make sure that you're thermodynamic integration protocol is fine and hence remove that variable from the way, could you please compute the logP for a polycyclic compound like NAPH (https://github.com/ricalessandri/Martini3-small-molecules/blob/main/models/itps/opt-poly/NAPH.itp) for which we can compare then make a comparison with the DG I obtained in https://doi.org/10.1002/adts.202100391.

[zhangtrzhangtr]

@ricalessandri
Using the same method and the itp file you provided for naphthalene, I calculated its transfer free energy between octanol and water. The results were ΔG(O→vacuum)=13.97 kJ/molΔG(O→vacuum)=13.97kJ/mol and ΔG(W→vacuum)=−5.96 kJ/molΔG(W→vacuum)=−5.96kJ/mol, leading to a transfer free energy of 19.93 kJ/mol between octanol and water. In the article https://pubs.acs.org/doi/abs/10.1021/ac9902291, I found an experimental logP value of -3.40 for naphthalene, which corresponds to an experimental transfer free energy of approximately 19 kJ/mol. This result is close to the transfer free energy obtained from thermodynamic integration. Additionally, I calculated the transfer free energy for another small molecule, naphtho[1,2-b]thiophene, mentioned in the article, which resulted in 25.74 kJ/mol, compared to an experimental value of 25.6 kJ/mol. This suggests that my thermodynamic integration method is likely sound. I suspect the discrepancy in the transfer free energy calculation for 2,1,3-Benzothiadiazole might be related to the choice of bead types, although I'm not entirely sure—it's just a guess.

[ricalessandri]

@zhangtrzhangtr Sorry for the slow follow-up. Great, then there are no problems with your TI protocol. I have the following logP for NAPH: ~18.5 kJ/mol, so it's ~1.5 kJ/mol different from yours. I've used 20 lambda points during TI, so you may try that out to see if your value changes. However, I think the free energy you get is in the right ballpark for sure.

This means that the beads type choices must be revised for 2,1,3-Benzothiadiazole. Your data indicate that the Martini model is too hydrophilic w.r.t. the exp. data. I'd try to scale down the hydrophilic beads for the beads involving the nitrogens. Looking at Suppl. Table 19 of the M3 paper:

going from TN6a-TN6a to TN4a-TN4a should work. Let me know.

[zhangtrzhangtr]

@ricalessandri Yes, this approach was effective. After replacing the two nitrogen-containing beads with TN4a-TN4a, the calculated transfer free energy between octanol and water was 10 kJ/mol, which is very close to the experimental value of 11 kJ/mol. Thank you very much for your assistance and guidance.