[lammps-users] DREIDING force field+charge equilibration

Dear LAMMPS community:

I would appreciate it if anyone can help me to find an answer to the following question.

How charge equilibration can be performed along with DREIDING force field? As you may know, DREIDING force field itself does not calculate partial charges and not perform charge equilibration.

I have chosen a sample system containing 9 methane molecules in the unit cell with "atomic_style" of "full". I tried to use "fix qeq/reax" command to charge equilibration (QEq) method. Parts of data_file and input_file are attached. NVT-MD simulation at 300 K failed after several time steps, giving the following error:

ERROR on proc 16: Bond atoms 27 29 missing on proc 16 at step 154901
ERROR on proc 20: Bond atoms 17 19 missing on proc 20 at step 154901

Farshad

**************Input_File****************

# Initialization
units real
boundary p p p
atom_style full

# Dreiding potential information
neighbor 0.4 bin
neigh_modify every 10 delay 0 check no
bond_style harmonic
bond_coeff 1 350 1.09
angle_style harmonic
angle_coeff 1 50 109.5
pair_style buck/coul/cut 10.0
pair_coeff 1 1 88366.7 0.28 583.0
pair_coeff 1 2 17353.2 0.27 135.2
pair_coeff 2 2 3407.8 0.26 31.4

special_bonds dreiding

#Energy minimization-Steepest Descent method
minimize 1.0e-6 1.0e-6 1000 1000
min_style sd
write_restart min.sd.%

timestep 0.2

fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 param.qeq

#heating-up system from 0.1 K to 300 K (20 ps)
fix 3 all nvt temp 0.1 300.0 20.0
run 100000
unfix 3
write_restart nvt-0.1-300.%

#equilibrate system at 300 K (20 ps)
fix 5 all nvt temp 300.0 300.0 20.0
run 100000
unfix 5
write_restart nvt-equ-300-.%

************data_file************
LAMMPS data file. CGCMM style. generated by VMD/TopoTools v1.1 on Sat Jan 29 21:07:01 -0600 2011
45 atoms
36 bonds
54 angles
0 dihedrals
0 impropers
2 atom types
1 bond types
1 angle types
0 dihedral types
0 improper types
-7.750153 7.249847 xlo xhi
-6.925005 8.074995 ylo yhi
-7.278439 7.721561 zlo zhi

Masses

1 12.011000 # C
2 1.008000 # H

Atoms

1 1 2 0.000000 -2.086818 -1.097731 3.007030 # H 1
2 1 1 0.000000 -1.896039 -1.425871 1.976185 # C 1
3 1 2 0.000000 -1.204409 -0.723656 1.491671 # H 1
4 1 2 0.000000 -1.449813 -2.429589 1.987393 # H 1
5 1 2 0.000000 -2.841882 -1.450986 1.418153 # H 1

Bonds

1 1 1 2
2 1 2 3
3 1 2 4
4 1 2 5

Angles

1 1 1 2 3
2 1 1 2 4
3 1 1 2 5
4 1 3 2 4
5 1 3 2 5
6 1 4 2 5

Dear LAMMPS community:

I would appreciate it if anyone can help me to find an answer to the following question.

How charge equilibration can be performed along with DREIDING force field? As you may know, DREIDING force field itself does not calculate partial charges and not perform charge equilibration.

so what do you expect to get out of this then. the charges are part of the
parametrization and thus should not be modified.

I have chosen a sample system containing 9 methane molecules in the unit cell with "atomic_style" of "full". I tried to use "fix qeq/reax" command to charge equilibration (QEq) method. Parts of data_file and input_file are attached. NVT-MD simulation at 300 K failed after several time steps, giving the following error:

well, that fix is meant to be used with the reax/c pair style and
does operations required for the reaxx parameterization.

ERROR on proc 16: Bond atoms 27 29 missing on proc 16 at step 154901
ERROR on proc 20: Bond atoms 17 19 missing on proc 20 at step 154901

this is an error that often happens when something unphysical
happens in your simulation. have you visualized your trajectory?
does it happen, when you disable fix qeq/reax?

in general, i don't see how what you are doing can lead to any
meaningful simulation, so it doesn't surprise me much that
your simulation becomes unstable.

axel.