TraPPE-EH Lennard Jones Potential

Hi all,

Just wondering if anyone has been using TraPPE-EH force field in LAMMPS.

I used force field parameters from this paper: http://static.msi.umn.edu/rreports/2005/186.pdf
But my LAMMPS simualtion doesn’t give a correct CCCC torsional angle distribtuion.

One of the obvious problems in my modeling also comes from an incorrect system density at its equilibrium state, the volume is much higher than it’s supposed to be.

The TraPPE-EH used C-H bond length of 0.55, which is half of the the distance between carbon and hydrogen (1.1A).

I build a test system with 180 propane molecules and simulated at 273K, 1atm: keep the same parameter below, except that the bond_coeff( taken from CHARMM force field.) using 0.55 CH bond length and 1.1 bond length two for different cases, but none of them give a correct density.

pair_coeff 2 2 lj/cut 0.0306 3.31
pair_coeff 1 1 lj/cut 0.01 3.65

bond_coeff 1 harmonic 222.5 1.535
bond_coeff 2 harmonic 309.0 1.10

angle_coeff 1 harmonic 58.764 112.7
angle_coeff 2 harmonic 37.766 110.7
angle_coeff 3 harmonic 33.234 107.8

dihedral_coeff 1 harmonic 1.434 -1 3
dihedral_coeff 2 harmonic 1.708 -1 3
dihedral_coeff 3 opls 1.42012 -0.2728 3.16528 0.0

I suspect the problem is about the lj potential parameters. I tried the OPLS-AA LJ potential and it gives a closer result to the real propane gas density. I think I may implement the TraPPE-EH force field incorrectly, but I couldn’t figure out what the problem might be.

Any suggestions would be helpful!

Thanks,
Endian

How much is the density off? How long
did you equilibrate with NPT? What units are
you using?

Steve

There are different pathways to troubleshoot your problem. One is to create only one molecule and dump out individual energy contributions. Propane is quite small thus you could work it out by hand. The other way is to do make sure you can implement more simple cases of the TraPPE FF. TraPPE-UA for example.

Go to https://www.chem.umn.edu/groups/siepmann/trappe/intro.php pick butane (no dihedral for propane in TraPPE-UA) and see if you can succeed there instead. As always, make sure you are accounting for the Coulomb interactions, mixing rules, etc the right way.

Carlos

Hi Carlos and Steve,

Thanks for trying to help! I attached my simulation input files as well as moltemplate file *.lt that were used for building the input data.

I tested several simulation for Propane gas(180 Molecules) at 273K, 1atm with different C-H Bond length. As suggested in their publications, they allow to take bond potential parameters from CHARMM force filed for harmonic bond interactions.
I fixed TraPPE-EH LJ parameters in all cases of simulation, but change them with different bond length and bond potential, see here:

Bond E r0 Final Volume
CC 303 1.5
2.3*10^7

CH 80.9 0.55
CC 303 1.5
2.3*10^7

CH 309 0.55
CC 222 1.5
1.1*10^7

CH 309 1.1
CC 303 1.5
1.1*10^7

CH 309 1.1

For these four cases, change value of Potential E doesn’t produce much different results, but half the equilibrium CH bond distance from 1.1 to 0.55 will increase the final volume(nearly doubled the gas volume). But the expected volume for propane in this case should be 6.5*10^6.

I suspect that using 0.55 for CH bond will not give the system correct properties.

trappeEH.lt (1.52 KB)

system.lt (296 Bytes)

propane.lt (1.2 KB)

NPT.run (757 Bytes)

log.lammps (72.1 KB)

system.in.settings (396 Bytes)

system.in.init (328 Bytes)

system.data (233 KB)

Some comments,
1-) No time/interest in debugging your scripts. That is a struggle you have to experience yourself :wink: When faced with the right attitude good things come out of tough struggle.
2-) Certainly the non-bonded parameters are the ones responsible for the body of the liquid/gas. Its no so much where the interaction center gets placed as long as the non-bonded params are reset accordingly. I didn’t read the reference paper thus I ask, are you sure this particular parametrization can be applied to alkanes?
Carlos