Overprediction of Density by LAMMPS for R134 (1,1,1,2 Polytetraflouroethylene)

Please check that you intended to have the C-C bond flexible instead of constrained-length and that its bond constant is correct.

Also:

  1. Take out the line with xy yz xz if you don’t intend to deform the box. That line is forcing LAMMPS into triclinic instead of orthogonal box mode, which slows down the calculation for no benefit (and might even be the cause of the discrepancy, which would be valuable to know if true)
  2. It would’ve been good if we could have seen all this earlier!

I don’t intend to keep the C-C bond flexible. But, in lammps using shake, I can at most keep 3 bonds rigid. That is why bond #02 is missing. Besides, I thought that bond constants values that I have chosen (600 kcal/mol) would be high enough to constrain the bond anyways. I did not know how else to use fix shake to constrain all 4 bonds.

More things to consider:

  1. GOMC doesn’t sample bond length changes (…why?) but your LAMMPS script has flexible C-C bonds. I imagine your GROMACS script has fixed all bond lengths (which is easier there than in LAMMPS). But given the very high bond strength I can’t imagine this being significant.

  2. More importantly both GOMC and GROMACS use full 1-4 non-bonded interactions by default, while LAMMPS ignores them by default – that is, the default special_bonds setting is 0 0 0 in LAMMPS and (equivalent to) 0 0 1 in GOMC and GROMACS.

2 Likes
  1. Yes it was easy in GROMACS via LINCS to constrain all bonds easily. But still, yet I am not sure how to constrain more than 3 bonds in lammps via shake.

  2. Ohh!! I hope this is the prime reason for the deviation !!! Let me recheck and rerun the simulations if desired.

Best,
Asad

1 Like

I ran lammps with 1-4 interaction turned on and error as compared to GOMC reduced to 0.4%.
I think this overprediction was due to 1-4 interaction turned off.

On the other hand, just to confirm, I am also running GOMC with 1-4 interactions turned off to see if it gives a greater density as lammps was giving before.

Thank you for your help all along.

All the best! Since you know GROMACS reproduces the GOMC result with all bond lengths fixed, you might get the clearest comparison by running LAMMPS with 1-4 interactions turned off and GROMACS with the C-C bond length changing instead of held constant.

Also bear in mind that the model with fixed C-C bond is technically a different model from having a flexible C-C bond – but with the GROMACS runs you shouls be able to bound the difference.

1 Like

Hello there Shern,

So, I have been struggling with the density difference still now and mainly focusing on dihedrals. I have checked all angles and they are sampled correctly by both lammps and gomc.
There are 2 dihedrals defined in the forcefield FCCF & FCCH. As a quick reminder, my system is R134 (1,1,1,2 PTFE). According to the forcefield, both dihedrals should sample -180, -60, +60, +180 degrees.
When I checked GOMC (.dcd analysis via VMD), all dihedrals were sampling all angles (as they should).
But when I checked Lammps dcd file, I found that each dihedral is sampling separate angles only. Details are in this pdf (if statement is not clear). There are 02 types of dihedrals and total 09 dihedrals in one molecule.
Any idea why would this be? I have checked and double checked the dihedral constants in the input lammps files which are correct. Any insights are appreciated.
Probably this anomaly could be the reason of density difference between GOMC and lammps.

Dihedral Comprehensive Results for Lammps.pdf (78.6 KB)

Best,
Asad

See the note dihedral_style command — LAMMPS documentation :

When comparing the formulas and coefficients for various LAMMPS dihedral styles with dihedral equations defined by other force fields, note that some force field implementations divide/multiply the energy prefactor K by the multiple number of torsions that contain the J–K bond in an I–J–K–L torsion. LAMMPS does not do this (i.e., the listed dihedral equation applies to each individual dihedral). Thus, you need to define K appropriately via the dihedral_coeff command to account for this difference if necessary.

It is quite likely that your LAMMPS dihedrals are stiffened by a factor of three relative to the GOMC dihedrals, which would explain the lack of torsional exploration about the central C-C bond. I do not know GOMC or its documentation well enough to say from the GOMC end whether this is the case.

Thank you @srtee
I have looked into this possibility of different pre-factors and made sure that I am taking care of the correct pre-factor by comparing original forcefield dihedral equation with lammps equation. But maybe, diherals are stuck at any one minima and that is why only 1 sampling with each dihedral.
Thanks again.

Asad