Question regarding Napthalene C (atom type 92) in OPLSAA

[Nekkrad]

Dear moltemplate developers,

I am trying to simulate Anthracene with OPLS-AA force field in LAMMPS.

To generate the data files needed I assigned the following atom types:

• 90 (Aromatic C) for atoms 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
• 91 (Aromatic-H) for H atoms bonded to 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
• 92 (Naphthalene fusion C) for atoms 4a, 8a, 9a, 10a

I was checking the force field parameters against the one at boss4.9 documentation and I have a couple of questions.

I do not find in the boss documentation the naphthalene C fusion ring. The closest I was able to see is the C, R4C.
The non-bonded parameters are the following

Atom type	q	$\sigma$	$\epsilon$
90 Aromatic/Benzene C	-0.115	3.55	0.07
92 Napthalene Fusion C	0.0	3.55	0.07
R4C	0.0	3.5	0.066

Where do these Naphthalene fusion ring values come from? It looks like an Aromatic C without the charge. It makes sense to me chemically but I could not find it in the OPLSAA parameters.

The second question arises from the bond terms.

Bond	Coeff1	Coeff2
Bond_coeff 1	469.0	1.4
Bond_coeff 2	367.0	1.08
bond_coeff @bond:090_091	337.0	1.449

From the System.in.settings I found the first two bond terms which should reflect the C-C and C-H bonds, the last row shows what I think is the bond term for the C-H in an aromatic ring (taken from oplsaa.lt)

1. The fact that I have only 2 bond terms makes me think that indeed atom 92 is in fact an atom type 90 with 0 charge. Is this true?
2. Am I wrong in saying that the last row is the C-C bond term? I am confused because it says 090-091 which are atom types linked to C-H, but values which are consistent with C-C

My chemical knowledge tells me that 1.08 is indeed a C-H bond length and 1.4 should be the one for C-C in benzene.
Is there anything I am missing on how moltemplate works?

EDIT:
I generated the OPLSAA prm files for NAMD with libpargen . (OPLSAA_anthracene_Libpargen_NAMD_prm.txt)

I notice that there are extra terms compared to my settings for LAMMPS, for example:

Bond

• C810 C805 520.000 1.370. It appears only twice so I guess is the C-C bond term between fusion ring atoms.

Angle

• C802 C803 C804 85.000 134.900
• C804 C805 C810 85.000 117.300

I reckon there could be more in Dihedral and non bonded-term. I will investigate further

EDIT 2: I have modified the oplsaa.lt file in moltemplate to generate lammps files with the terms I think are missing except the improper dihedrals (I did not really understand how they are implemented).

My lammps simulations are crashing due to atoms out of range (bad dynamics probably coming from mistakes in the new ff parameters)

I will try to investigate further but I am kinda at my wits end to be fair.

Thank you again for your time and help.

[hothello]

Sorry for the super late reply!
I believe the R₄C type corresponds to @atom:84, which has the same Lennard-Jones parameters. Here is the list of carbon atom types that have zero charge:

    set type @atom:84 charge 0.0  # "Alkane >C<"
    set type @atom:86 charge 0.0  # "Alkene R2-C="
    set type @atom:92 charge 0.0  # "Naphthalene Fusion C"
    set type @atom:120 charge 0.0  # "Diene =CR-CR="
    set type @atom:159 charge 0.0  # "Methanethiol CH3-SH"
    set type @atom:200 charge 0.0  # "Imide >CH-CONHCO-"
    set type @atom:373 charge 0.0  # "Ethyl Anion CH3-CH2-"
    set type @atom:429 charge 0.0  # "Aryl Sulfonamide C-SO2N"
    set type @atom:460 charge 0.0  # "Biphenyl C1"
    set type @atom:654 charge 0.0  # "Cyclopropane -CR2-"
    set type @atom:733 charge 0.0  # "Amine CH3-NH2"
    set type @atom:769 charge 0.0  # "Alkyne RCCR"
    set type @atom:848 charge 0.0  # "Tertiary Amide -NRCH3"
    set type @atom:858 charge 0.0  # "B3-Pep CB GLY Main/C-Ter"
    set type @atom:897 charge 0.0  # "Triene R2-C= (mid C=C)"

Regarding the second question, I recently stumbled upon the same bond term and was equally confused. Actually, you have to first consider the following substitutions:

  replace{ @atom:90 @atom:90_b048_a048_d048_i048 }
  replace{ @atom:91 @atom:91_b049_a049_d049_i049 }
  replace{ @atom:92 @atom:92_b048_a048_d048_i048 }

The bond coefficients that apply to these three atom types are defined by the rules:

  write_once("Data Bonds By Type") {
    @bond:048_048 @atom:*_b048*_a*_d*_i* @atom:*_b048*_a*_d*_i*
    @bond:048_049 @atom:*_b048*_a*_d*_i* @atom:*_b049*_a*_d*_i*
  }

Therefore, the bond coefficients that apply to Naphthalene are:

    bond_coeff @bond:048_048  469.0 1.4
    bond_coeff @bond:048_049  367.0 1.08

To debug your system, please share the LT file with the molecular topology. Note that the simplest and most robust way to use the OPLS-AA force field is to let Moltemplate determine the molecular topology starting from the connectivity. You only need to assign the atom types and to specify the list of bonds with:

write('Data Bond List') {
  $bond:C1 $atom:C1 $atom:C1
  $bond:C2 $atom:C1 $atom:C9a
  ...
}

I made a tutorial on using Moltemplate plus OPLS-AA, which you can find on the newest version of the LAMMPS manual, section 8.6.4.
I hope you have figured it out by now. If not, here we go!

[Nekkrad]

Thank you for your reply!

In the end I have pivoted to a different system!

In any case I think that there are still some issues with Naphthalene and related systems.

1. Only two types of bond are assigned to Naphthalene (the ones you mentioned), but I think there should be atleast 1 additional bonded term (there 3 types of atom, Naphthalene fusion ring, Aromatic C and Aromatic H).
2. The same issue applies to Angles, Dihedrals and Improper Dihedrals

This implies that for Moltemplate the Naphthalene fusion ring C is considered to be an aromatic C with the non-bonded term being the only different part.

This didn't make much sense to me as a chemist and when I downloaded the FF parameters for Naphthalene from Jorgensen's Website (attached in the previous post), I did find extra terms that are missing from Moltemplate, or atleast I think so.

The possibility of having done a mistake setting up the .lt file for moltemplate is non-zero. I will attach my .lt to this post as soon as possible.

Thanks again!

[hothello]

The OPLS force field is designed to have universal atom types, therefore it is expected to trade off some chemical specificity. In the naphthalene example, the aromatic carbon atoms are all equivalent except for the charge. Therefore only the two C-C and C-H bonding terms are specified. The same logic applies to the higher bonding terms, as you can see from the renamed atom types.

Just try to set up the naphthalene molecule using the Data Bond List and let Moltemplate do the rest. I am quite sure this is the correct way to proceed.

[Nekkrad]

Yes and it works perfectly fine.

What I am confused about is that in both moltemplate and TINKER (https://github.com/TinkerTools/tinker/blob/release/params/oplsaa.prm) you only have two bond stretching terms , but in BOSS version 4.8 (actually 4.9, but according to the documentation no changes have been made to aromatics) and also on libpargen Naphthalene (Anthracene) has an extra bond term which refers to the Fusion C-Fusion C bond term (520.000 1.370).

I have tried to understand if this comes from moltemplate using an older OPLS-AA version compared to libpargen and BOSS 4.9, but it does not seem to be the case.

As I said I am not really working on Naphthelene anymore but I think it could be useful for others.
This are the lammps data, anthracene lt and settings file (do not have naphtalene but it should be fine)

system_data.txt
anthracene_lt.txt
system_settings.txt

The Anthracene param file I attached in my first post has an extra term. Where does this come from? That is my question. The file was taken using LibPargen and it is the same when using BOSS version 4.9.

[hothello]

Sorry, I have missed your point.

You are right, the LibParGen server does distinguish between aromatic and fusion C atoms. The OPLS-AA version which is distributed with Moltemplate is consistent with the one distributed with Tinker, which says that the

parameters are taken from those distributed with BOSS Version 4.8.

I believe the mismatch originates from using different versions of the OPLS force field. To force the distinction between the aromatic and fusion C atoms, you need to alter the renaming step to something like:

  replace{ @atom:92 @atom:92_b048F_a048_d048_i048 }

And then define a new bond by type rule:

  write_once("Data Bonds By Type") {
    @bond:60-60 @atom:*_b48_a*_d*_i* @atom:*_b48F_a*_d*_i*
  }

And do the same for the angle and dihedral terms, if there are specific terms. Thanks for raising this point: eventually, it will be necessary to assign a version to the OPLS-AA force field distributed with Moltemplate and perhaps update the parameters and rules to a more recent version.

[Nekkrad]

Yes,

This is what I tried and got stuck with adding the new improper dihedrals. I will attach the relevant files when I find them

On a side note Jorgensen just released a new OPLS-AA paper with all the parameters available in the SI. It might be a good time to implement them. I am very happy to help. (https://pubs.acs.org/doi/full/10.1021/acs.jpcb.3c06602)

[hothello]

Thank you very much for the suggestion: I would also involve @jewettaij in this work. I hope he is not too busy, but let's get this done soon.

[jewettaij]

Thank you Otello for reaching out to Nekkrad!

My apologies for my absence. I left science 2 years ago. I did not expect to be so busy that I could not respond to questions.

Truthfully, this is what I was afraid would happen. It's hard to keep up to date in the field, and that makes it hard to keep things updated. I've lost easy access to download papers. (I miss scihub.) I no longer even remember the procedure for calculating atomic partial charges. The only reason OPLS is available in moltemplate is that that Jason Lambert and I wrote a script to convert TINKER .prm files into moltemplate format. But the Ponder lab has not updated these files since Boss 4.8.

This is what needs to happen: I should write a new converter program which can read the force-field files directly from Jorgensen's web site:
https://traken.chem.yale.edu/
...and convert them to moltemplate format. That would make it much easier to keep moltemplate up-to-date with OPLS (and other force-fields too). But past experience teaches me that it would take a week. And I work "9-9-6" hours these days. I don't see it happening in the next 6 months.

I really hope I did not waste your time! The way things are now, I worry that moltemplate does more harm than good. I don't want people using out-of-date force fields. Since I can't keep these force fields up-to-date, I will clearly state exactly which version of OPLS is supported by moltemplate (the 2008 version). If you think I should remove support for OPLS in moltemplate, I can do that too. I don't want people to have to suffer they way you did.

I am grateful to both of you for your hard work, trying to use this software and keep it alive.
I haven't completely given up on it. Moltemplate is the most useful thing I have ever done. I will resume work on it eventually.

[hothello]

Hey Andrew! Moltemplate is a great tool, and you cannot be responsible for maintaining third-party force fields. It would be helpful to add a version to every force field shipped with the code, but ultimately it's a user's task to implement, check, and verify a force field.
As the nature of this project is open-source, it could also be worth considering opening the current force field library to the contribution of other users and developers, e.g. for the ATB force field.
Welcome back, we missed you!!

[jewettaij]

A new version of the OPLSAA force field (from late 2023) is being converted to moltemplate format. For details, see #104