CHARMM Development Project Forums User Discussion & Questions Parameter Set Discussion 1-4 nonbonded interactions mixing rules in CGenFF

Hello,

I would like to ask a question regarding the mixing rules for Lennard-Jones interactions between 1-4 atoms. I am currently examining polyaromatic molecules grafted with aliphatic side chains and want to model them using CGenFF. As a result, in my simulations, two major atom types are the CG2R61 and CG321 types.

In the CGenFF prm file, one finds the following LJ entries regarding the aforementioned atom types:
CG2R61 0.0 -0.0700 1.9924
CG321 0.0 -0.0560 2.0100 0.0 -0.01 1.9

From what I understand, when two CG2R61 atoms are 1-4 neighbors, the values epsilon=-0.0700 and sigma=1.9924 are used without scaling and when two CG321 atoms are 1-4 neighbors, the values epsilon=-0.01 and sigma=1.9 are used; also unscaled.

What happens in the case of a CG2R61-CG321 1-4 pair? Are the parameters (-0.0700, 1.9924) and (-0.01, 1.9) combined via typical Lorentz Berthelot mixing?

I appreciate any input on this subject.

Thanks in advance.

When additional 1-4 terms are not given, it is assumed that they are the same as the primary epsilon and Rmin (not sigma!) values.

You may wish to devise a simple test, perhaps using SKIPE (energy.doc) to turn off most other energy terms.

Dear Dr. Venable,

Thank you very much for your reply.

Unfortunately, I cannot perform the test you propose, since I do not have access to the CHARMM program. I am currently implementing CGenFF for an in-house program that creates initial configurations to be used as input in atomistic simulations. So, my question is more force-field-oriented than CHARMM-oriented.

I will take into account your reply and apply it in our code.

I thank you again,
Orestis Ziogos

While I think your interpretation is right, the test Rick is suggesting is absolutely necessary to have faith in your code - too many things can go wrong. Coding stuff like this is like doing a very complicated calculation manually on paper; if you don't do a sanity test, you can't trust the results at all.

I'm a bit afraid to suggest it here, but if you don't have access to CHARMM, you might want to do the test with NAMD, which can read CHARMM parameter files directly and faithfully reproduces the energetics (at least at the short range, but that's all you need in this case). Obligatory disclaimer: we cannot support you on this forum with the practicalities of using NAMD, as the latter has its own support channels.

I managed to run some characteristic validation runs using NAMD and the latest CGenFF version on a series of isolated molecules that have additional 1-4 LJ terms (n-hexane and butylbenzene to name some).

It turns out that when secondary 1-4 parameters are not listed, the primary 1-4 parameters are used and are mixed via typical Lorentz-Berthelot rules, just as Dr. Venable suggested.

I would like to thank you again for all the input.

Orestis

This is as you originally thought, right? Also, am I right to assume your validation was a comparison between NAMD and your in-house implementation?

Yes. To be more precise, I cross checked the results from NAMD with calculations carried out by both our in-house program and LAMMPS.

The only problem with LAMMPS is that it does not support the shifting of Coulomb interactions, although it supports the typical CHARMM/Xplor switching function for LJ interactions. So, in order to be able to compare both interaction types, I deactivated the shifting in NAMD and made sure that the cutoff for every molecule was sufficient in order to account for all interactions.

The maximum differences between NAMD and our code and LAMMPS (our code yields exactly the same results with LAMMPS) for LJ and Coulomb interactions was 0.06% and 0.004% respectively. I believe the discrepancy is due to decimal rounding and minor programming differences between the packages.