Simulation of a h2 box

Dear LAMMPS users

I am trying to simulate a box of H2 at room temperature. According to the literature that I found, the AIREBO potential should be used, however, it is used when there is graphene or similar and I do not understand how it should be used when only hydrogen is available. Is there any suggestion of a potential for lammps to realistically simulate a box of molecular hydrogen (H2)?

# ==================initial parameters ==========
variable dt equal 1
*variable Text equal 300.0 *
variable Pext equal 1.01325

units real
boundary p p p

atom_style full
bond_style harmonic

read_data data.hydrogen_H2
replicate 10 10 10
#=================Potential (Force Field) ===============

pair_style lj/cut/coul/long 9.0
pair_coeff 1 1 103.5 0.400 # H-H

kspace_style pppm 1e-6

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes one 11330 page 200000
timestep ${dt}

#============================= Thermo ===================
*thermo 10 *
thermo_style custom step temp press pe ke etotal ecoul epair lx ly lz vol
thermo_modify flush yes

#========== Minimization ===============================
*minimize 1.0e-4 1.0e-8 1000 10000 *

#=================== Equilibration ======================
*dump 1 all atom 10 H2.lammpstrj *
*fix 1 all npt temp {Text} {Text} (5.0*dt) iso {Pext} ${Pext} 10.0 *
run 4500

unfix 1

fix 2 all nvt temp 300 300 $(5.0dt)*
run 4500

unfix 2
fix 3 all nve
run 1000

write_restart restart.equil

I attach a script and the data.file.

data.hydrogen_H2 (447 Bytes)
in.H2 (1.3 KB)



It is quite possible to use AIREBO for just H2 molecules. The potential parameters for H-H interactions (inter- and intra-molecular) are present. The mapping of element to atom type is as for other manybody pair styles (please check in the documentation):

pair_style airebo 3.0 1 0
pair_coeff * *  CH.airebo H

However a few things need to be observed.

  • AIREBO uses metal units not real. That makes no difference for positions, but for time (picoseconds instead of femtoseconds)
  • there must be no explicit bonds. In AIREBO bond order and thus the bond topology is computed dynamically during a run. It is a reactive potential.
  • You cannot use a Kspace style with it (I don’t see the need in the first place and your data file with charged hydrogen atoms is bogus, too. LAMMPS warns about a total charge of 820. That is very bad).
  • You are also wasting some simulation time by not assigning an initial, randomized velocity, Since you have a perfectly ordered system, it takes significant simulation time until the symmetry is broken and the system mixes properly.

Whether AIREBO is a sufficiently realistic model for your purposes is for you to determine by analyzing the resulting trajectory after equilbration and comparing the results with expected results.

If you are looking for alternatives, you should continue your survey of the published literature. I am quite confident that there are more parameterizations available (and not only for plasmas).