Buckingham potential

Hello,

I’m currently trying to module PVDF and obtain a working potential for further studies. I’m trying to utilize potentials that are stated in literature to compare results and in some of the references they use a buckingham potential for their pair style. When I utilize this function with their constants lammps produces

an error.

When I utilize an npt it gives me this error:

ERROR: Non-numeric pressure - simulation unstable (…/fix_nh.cpp:1029)

When I use a nvt, it gives me this error:

ERROR on proc 2: Bond atoms 2208 2211 missing on proc 2 at step 10 (…/ntopo_bond_all.cpp:63)

Here is the first part of my script for lammps:

#Settings

VARIABLES

variable T equal 463.15
variable V equal vol
variable p equal 400 # correlation length
variable s equal 5 # sample interval
variable d equal 2000 # dump interval
variable dt equal 1.0 # time step

#Converting Variables

variable kb equal 1.3806504e-23 #[J/K/** Boltzman Constant]
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal {atm2Pa}*{atm2Pa}{fs2s}*{A2m}{A2m}*{A2m}

variable fname index lammps_data_mod
variable simname index lammps_data_mod

Initialization

units real
boundary p p p
atom_style full
read_data lammps_data_mod

Potential Information

neighbor 0.4 bin
neigh_modify every 10 one 10000
bond_style harmonic
bond_coeff 1 327.6 1.085
bond_coeff 2 499.3 1.357
bond_coeff 3 308.9 1.534
angle_style harmonic
angle_coeff 1 42.9 108.45
angle_coeff 2 90.0 107.74
angle_coeff 3 38.5 109.27
angle_coeff 4 120.0 105.27
angle_coeff 5 62.0 107.1
angle_coeff 6 80.3 118.4
angle_coeff 7 80.3 118.4
dihedral_style fourier
dihedral_coeff 1 6 0.395 1 180 0.72 2 180 -0.38 3 180 -0.205 4 180 0.425 5 180 -0.025 6 180
dihedral_coeff 2 3 0 1 180 0 2 180 0 3 180
dihedral_coeff 3 3 0 1 180 0 2 180 0 3 180
dihedral_coeff 4 3 0 1 180 0 2 180 0 3 180
dihedral_coeff 5 6 0.355 1 180 0.345 2 180 -0.38 3 180 0.14 4 180 0.145 5 180 -0.025 6 180
dihedral_coeff 6 3 0 1 180 0 2 180 0 3 180
pair_style hybrid buck/coul/cut 10.5 lj/cut 10.5
dielectric 1.0
pair_coeff 1 1 buck/coul/cut 14976.0 0.323624595469 640.8 10.5
pair_coeff 1 2 buck/coul/cut 4320.0 0.292825768668 138.24 10.5
pair_coeff 1 3 buck/coul/cut 45094.0 0.26191037427 260.77 10.5
pair_coeff 1 4 buck/coul/cut 14976.0 0.323624595469 640.8 10.5
pair_coeff 2 2 buck/coul/cut 2649.6 0.267379679144 27.36 10.5
pair_coeff 2 3 buck/coul/cut 12300.0 0.241365161353 53.88 10.5
pair_coeff 2 4 buck/coul/cut 4320.0 0.292825768668 138.24 10.5
pair_coeff 3 3 buck/coul/cut 135782.0 0.221126417973 106.12 10.5
pair_coeff 3 4 buck/coul/cut 45094.0 0.26191037427 260.77 10.5
pair_coeff 4 4 buck/coul/cut 14976.0 0.323624595469 640.8 10.5
pair_coeff 4 5 lj/cut 0.173 3.600 10.5
pair_coeff 4 6 lj/cut 1.165 3.310 10.5
pair_coeff 1 5 lj/cut 0.173 3.600 10.5
pair_coeff 1 6 lj/cut 1.165 3.310 10.5
pair_coeff 2 5 lj/cut 0.063 2.600 10.5
pair_coeff 2 6 lj/cut 0.425 2.310 10.5
pair_coeff 3 5 lj/cut 0.126 3.087 10.5
pair_coeff 3 6 lj/cut 0.849 2.797 10.5
pair_coeff 5 5 lj/cut 0.200 3.200 10.5
pair_coeff 5 6 lj/cut 1.345 2.910 10.5
pair_coeff 6 6 lj/cut 9.052 2.620 10.5

There are more settings than just the force field parameters that determine the stability of an MD run: time step, neighbor list settings, thermostat parameters, initial geometry, choice of additional fixes.

Axel.

Maybe it is due to poor initial conditions. Have you visualized your initial system?

@Stephan Paquay: Yes I have visualized my assembly; what should I be looking for to see if its unstable? Also, when using all the same parameter but instead using a Lennard-Jones potential, I don’t get this error.

Thanks,

If the simulation is unstable you will observe some atoms acquiring an unrealistically high velocity or fly out of the simulation box. The fact that it works for a Lennard-Jones potential however seems to imply that the initial conditions should be OK. Does the simulation work for a smaller time step size?