Ligand ff based on exsiting ligand ff (which way is faster)

GROMACS version: 2023.3
GROMACS modification: Yes/No

Dear All

introduction:

I have found a fantastic article “CHARMM-DYES: Parameterization of Fluorescent Dyes for Use with the CHARMM Force Field” with ff parameters for various dyes.
One can download charmm36_dyes.ff, place it in the gromacs directory and use it via
gmx pdb2gmx -f ‘pdb_of_your_dye’ - filename.gro, which will create topology files.
It is worth mentioning that the pdbs of all dyes from the article are given.

my goal:
Is to modify one of the existing dyes (let’s say to add an extra two benzene rings) and simulate it (charmm force field, but whether NAMD or gromacs does not matter ).

What I can do:

One way to go is to create force field parameters from scratch following the NAMD fftoolkit tutorial.
cons: because the molecule is too big, I have to compartmentalize it.
It is too long. I already have done it before, and I am not even sure that the parameters that I have gotten are really good.
Also as I am writing this I realize I can upload mol2 file into ccgenff get the parameters there and change it manually, but this is also tedious.

On the other hand, I want to use their parameters.
And the extra parameters from benzene rings I can get myself (using fftoolkit).

So my question is, after all, how can I tie all this together?
Should I do it manually, and if yes, which files in the charmm36_dyes.ff should I modify? Is there an algorithm that I can use?
Maybe you can recommend a tutorial or even steps on how to approach this kind of task.

Is it easier to create topology and parameters files in charmm or in gromacs?

Is there an easier way to do it that I am missing?

Thanks in advance
Vardan.

The residue templates are stored in .rtp files. If you want to introduce a new residue, follow the format of these files. Then, if you need hydrogens to be added automatically, you might also need to take care of .hdb.

If the modifications are one or two benzene rings, then perhaps you can get a good guess of atom types and e.g. extra dihedrals from CGenFF. You might still need to derive and symmetrize atomic charges from QM if you want to be precise about it. Then, the charges and types go into .rtp, while extra bonded parameters (such as dihedrals or bonds) should go into ffbonded.itp.

Overall, if the original .ff is well structured, it shouldn’t be much work in Gromacs.

Dear Milosz

Thanks for your reply; now I know how to approach this issue.
Also I started to read Gromacs documentation on how to add new residues.

best,
vardan.

Dear Milosz and all

I have one more follow-up question, although it might be slightly out of the gromacs topic but related to the same issue that I have already asked ( Ligand ff based on existing ligand ff (which way is faster)).

I have a new residue based on a residue cy5 (see attached image10.png or Figure 1, from " CHARMM-DYES: Parameterization of Fluorescent Dyes for Use with the CHARMM Force Field").

I have already tried what we have discussed before:
I took half of the residue, added hydrogen and carbon, made the corrections in the rtp files, and .ff worked fine during gmx pdb2gmx command.

The issue: how to obtain the optimized structure of the modified residue?

Let’s say I added a Benzin ring and took out the tail, following the same geometric optimization as in the article. However one of the dihedral angles in the ‘bridge’ between Nitrogens is 180 degrees rotated.

Do I understand correctly that if all the parameters are fine and only the dihedral is rotated, then I should not worry, or is that a mistake?

And if that is a mistake, I was thinking of compartmentalizing the molecule doing the geometric optimization for each part, and then connecting everything together.

The problem is I don’t know any software that connects the correct pdbs of the optimized sections. Any help on this topic?

Best,
Vardan.

The issue: how to obtain the optimized structure of the modified residue?

Through QM optimization, usually done with Gaussian or Orca on the DFT level (B3LYP/6-31+G* would be a reasonably cheap compromise). That’s just to get the geometry + charges though.

Do I understand correctly that if all the parameters are fine and only the dihedral is rotated, then I should not worry, or is that a mistake?

If your conjugated linker switches from all-trans to 12-cis or 13-cis (according to the numbering in the graphics), it might affect the assignment of charges. So I’d first make sure your isomers correspond well to what is found in the experiment.

I was thinking of compartmentalizing the molecule doing the geometric optimization for each part, and then connecting everything together.

In principle any localized changes to the molecule shouldn’t affect the chemistry far away, but your dye is extremely conjugated, so without specific knowledge of where the extra group goes, I wouldn’t discard the possibility that adding a phenyl ring might change things far away, even if slightly. Ask a local chemist, they are good at this :) Especially given that your scheme is just one of the possible resonance structures of the molecule.

The problem is I don’t know any software that connects the correct pdbs of the optimized sections. Any help on this topic?

For connecting two parts of the molecule, the easiest way is to keep a common substructure in both, and do a simple alignment on the resulting geometries, then remove the overlapping atoms from one half. But here again, I’m afraid partitioning the conjugated system might yield subsystems with somewhat different properties.