Dear LAMMPS users,
I’m having difficulty using the CVFF forcefield to model polymers, specifically crosslinked DGEBA and DETDA. My simulation is stable with only lj/cut applied, but when I take into account coulombic forces (lj/cut/coul/long or lj/cut/coul/cut) the hydroxyl and amine Hs overlap with N and O atoms causing the system to eventually crash. I made the data file using msi2lmp and I double-checked all the parameters and partial charges according to the cvff.frc file. I’m running version LAMMPS 64-bit 8Feb2019.
I can get around the problem by setting my hydroxyl O-H partial charges to zero, or by using ‘neigh_modify exclude’ to prevent ON-H interactions. But its obviously incorrect to exclude or modify the parameters given. Is it possible the values in cvff.frc aren’t the right units for LAMMPS, or maybe a special_bonds command needed for CVFF, or something specific to hydrogen bonding I’m not accounting for?
Below is my code. Any feedback is appreciated.
Thank you,
Doug Pratt
#-----Initialization------------------------------
units real
dimension 3
boundary p p p
atom_style full
timestep 1
neigh_modify delay 0 every 1 check yes
#-----Force Field (CVFF)--------------------------
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style cvff
#Coulomb and 12-6
#pair_style lj/cut/coul/cut 10.0
#Just 12-6
pair_style lj/cut 10.0
read_data data.16x8_cvff_13_5
#-----Exclude NO-H interactions-------------------
group ON type 10 20 21 4 8
group HnHo type 25 36 13
#neigh_modify exclude group ON HnHo
#-----Output--------------------------------------
thermo 50
thermo_style custom step pe ke etotal evdwl ecoul temp density press
#dump 1 all cfg 10 crosslinked*.dump mass type xs ys zs
fix 1 all nvt temp 300 300 $(dt*100)
run 500
Dear Doug,
something similar happened once to my simulation when using GAFF force field. Not sure about CVFF, but some force fields assume that such hydrogen atoms are in the van der Waals field of the bonded heteroatom and the pair coefficients for the hydrogens are zeros. I guess, if you set the LJ-coefficient for the hydrogen to a small number (say, 0.001 for epsilon and 0.55 for sigma) , this should not affect your results (NO GUARANTEE, you change the FF at your own risk) , but will prevent such collisions. Unfortunately, you didn’t append any data from the data file, so I cannot be completely sure here.
Alternatively, you could try to decrease your timestep.
I also noticed that you start directly at 300 K. Try to heat up first with
fix 1 all nvt temp 1 300 $(dt*100)
run 10000
Or set the initial velocities with
velocity all create 300.0 123456 rot yes dist gaussian
Good luck!
Arkady
Hi Arkady,
My O-H and N-H hydrogens were indeed 0 0 pair coefficients. It felt like these numbers were the problem since (to my understanding) with no LJ repulsion term there’s nothing to push against the coulombic attraction between the atoms. Changing them to your suggested values solved the problem. Do you happen to know if this overlapping is a LAMMPS-specific occurrence? It doesn’t make sense to me why all CVFF papers give these Hs zero value parameters when it causes this error, especially for common structures like hydroxyl and amine groups.
Thank you very much for your time,
Doug
Dear Doug,
as I mentioned, these hydrogens are supposed to be in the VdW field of the bonded heteroatom (N,O). So, in theory there must be an equilibrium between bond stretching of O-H/N-H and O-N Coulombic repulsion. It actually should work with zero coefficients.
So, double check your system before changing anything. I personally switched to a FF with a small repulsive LJ term for H (GAFF2) and everything is working just fine since then.
Arkady