A bit of a strange situation… I have a simple OPLS-AA setup for bonded and nonbonded interactions for boron nitride (usually sheets or tubes). The topologies are typically generated with x2top and after a few minor tweaks everything is ready to go. I decided to migrate this setup to CHARMM (July 2022 version) and copied the corresponding parameters into ffbonded and ffnonbonded. There is really not much going on there: two atom types (nearly +/-1 charges, if that matters here), bonds, angles, dihedrals, all very simple.
Upon relaxing boron nitride nanotubes in vacuum, I find that some of the bonds are insanely compressed (1.3 A while the equilibrium distance is 1.446 A), which then leads to tube collapse within picoseconds of finite-temperature simulations. This is enormous strain and nothing like this is happening under OPLS-AA. Moreover, there is just no physical reason for any of this to happen.
Granted, casually moving interactions from one forcefield to another for a biomolecule is a terrible idea. However, this is a crystal structure and there is literally one bond type and a handful of angles. Any ideas about what might be going on here? Thanks!
Copying parameters like dihedrals won’t work because CHARMM uses a 1-4 scaling factor of 1.0 whereas OPLS uses 0.5, so the dihedral parameters are going to be completely wrong (different functional form notwithstanding).
Hi Justin! Different functional form isn’t an issue for the dihedrals, because I am manually making sure the function type is the same and the list of dihedrals is identical. I did however forget to remove the pairs section from the topology, thank you.
To be clear, the scaling factor you mention is unrelated to non-bonded business between molecules, correct? Just trying to make sure there’s nothing funny going on in addition to 1-4 and I don’t need to rescale any charges/epsilons in the nonbonded params.
The combination rules of the force fields are different, so I don’t know that you can assume that any ported values of ε and σ are valid. But what I was referring to was exclusively the 1-4 interactions. You certainly can’t assume the dihedral parameters work with what you’re doing since the scaling is totally different.
Yes, the non-bonded mixing rules are different, but I am not really worried about that at the moment in this particular case. Most of the issue now is structural and indeed I am able to prevent tubes from collapsing by cranking the dihedral energy term way up, but that’s just wrong. For graphene, everything is reasonably stable (all atoms are uncharged), but for hBN the entire model is broken, because these solid-state models are not really designed with 1-4 explicitly in mind.
By the way, could you suggest a reasonable aromatic carbon type to steal parameters from in ffnonbonded for something like graphene? I could run a graphene PDB through CHARMM-GUI, but I don’t trust it for non-bio…
The Nanomaterial Modeler should be fine, I don’t know why it wouldn’t be considered trustworthy. I have no experience with those types of systems, but I would expect a stock aromatic carbon type would be appropriate, e.g. CA.
I’ve posted about Nanomaterials Modeler in the past. All bonded interactions are completely wrong for graphene, unfortunately. Interestingly, CHARMM’s CA carbon LJ parameters are identical to those in OPLS-AA.
I wish nanomaterials modeler had something for hexagonal boron nitride, but there is nothing…