Coulomb interaction calculations and time step for Ca and OH interactions

LAMMPS LAMMPS Mailing List Mirror
1

Dear All,

I’m simulating a jennite molecule (calcium silicate hydrate mineral) using COMPASS force field, but I have problems with the timestep which must be 1fs. Then, in order to understand the calculations behind in a simple system, I’m simulating HO and Ca, the code is showed below; but I have the same problem.

  1. I know that Coulombic interactions should be dominant in this system. Can somebody explain me how are these interactions calculated by LAMMPS (I’m using the default permeability 100%)? is correct to use partial charges to define the system?

  2. Using bigger time steps I have not numerical values, I believe that this is because of the higher values of the bond interactions. I would appreciate any help regarding to the possibility to increase this time step to 1fs.

Data file:

3 atoms
1 bonds
0 angles
0 dihedrals
0 impropers

3 atom types
1 bond types
0 angle types
0 dihedral types
0 improper types
-3.0 -1.0 xlo xhi
-1.0 1.0 ylo yhi
3.0 6.0 zlo zhi

Masses

1 40.078 # Ca
2 1.00794 # H
3 15.9994 # O

Atoms

1 1 2 0.410 -1.478802091 -0.767516298 5.706252374 #H
2 1 3 -0.410 -1.755135129 0.040436040 5.239274366 #O
3 1 1 2.000 -2.474013243 -0.360871658 3.149631413 #Ca

Bonds

1 1 1 2 #H: O

Input file:

#GENERAL PARAMETERS

log log.energymin_HO-Ca.csh
units real
boundary p p p
atom_style full

#STRUCTURE DEFINITION
read_data lmpdat.HO-Ca_compass
neighbor 2.0 bin

#INTERACTION DEFINITION
kspace_style ewald 0.0001
pair_style lj/class2/coul/long 12.5
pair_coeff 1 1 3.11000 1.56000 #Ca-Ca
pair_coeff 1 2 3.00514 0.18441 #Ca-H
pair_coeff 1 3 3.21200 0.34775 #Ca-O
pair_coeff 2 2 2.87800 0.02300 #H-H
pair_coeff 2 3 3.12418 0.03952 #H-O
pair_coeff 3 3 3.30000 0.08000 #O-O

bond_style class2
bond_coeff 1 0.9494 540.3633 -1311.8663 2132.4446 # H: O

dump dump1 all xyz 1 Ho-Ca-read.xyz
run 1
undump dump1

#ENERGY MINIMIZATION
velocity all set 0.0 0.0 0.0
min_style cg

thermo 1
thermo_style custom step ebond eangle edihed eimp evdwl ecoul emol epair etotal
thermo_modify norm no
minimize 1.0e-4 1.0e-6 120 1000

dump dump2 all xyz 1 Ho-Ca-dump2.xyz
run 1
undump dump2

NVT EQUILIBRATION

fix NVT_HOCa all nvt temp 298.0 298.0 0.01

#SETUP OUTPUT INFO
thermo 1

thermo_style custom step fmax temp etotal ke pe press lx ly lz vol
thermo_modify flush yes
thermo_style multi

dump dump3 all xyz 1000 OHCa-NVT.*.xyz

timestep 0.01
restart 5000 restart_NVT.8
run 10000
undump dump3
unfix NVT_HOCa

Thanks. Best wishes,

Ingrid Padilla

Joint School of Nanoscience and Nanoengineering
North Carolina A&T State University

2

Dear All,

I'm simulating a jennite molecule (calcium silicate hydrate mineral) using
COMPASS force field, but I have problems with the timestep which must be
1fs.

​why 1 fs?​

Then, in order to understand the calculations behind in a simple system,
I'm simulating HO and Ca, the code is showed below; but I have the same
problem.

1. I know that Coulombic interactions should be dominant in this system.
Can somebody explain me how are these interactions calculated by LAMMPS
(I'm using the default permeability 100%)? is correct to use partial
charges to define the system?

​LAMMPS implements regular force field models like many other MD codes.
whether the parameters you use are appropriate or not, should be
independent of​ the MD code.

​i don't see how you can get any kind of realistic information from your
input. your total system is charged and having a neutral OH component is
*very* doubtful, too.​ how did you obtain those parameters?

2. Using bigger time steps I have not numerical values, I believe that
this is because of the higher values of the bond interactions. I would
appreciate any help regarding to the possibility to increase this time step
to 1fs.

​i don't see how it should be possible in a system containing hydrogen
atoms to obtain a reliable time integration beyond with a timestep
0.1-0.2fs without constraining bonded interactions. please look up in your
favorite MD text book how the maximal possible time step depends on the
highest occurring vibration.

axel.