Adding Hydrogen Atoms to a Box of Aluminum

Hi,

I am looking to add 5-10 individual hydrogen atoms to the lattice of a slab of aluminum. My code for creating the aluminum is below:

generate the box and atom positions using a bcc lattice

variable a equal 4.05

boundary p p p

lattice fcc $a
region box block 0 5.0 0 5.0 0 5.0
create_box 1 box
create_atoms 1 box
change_box all triclinic

Need to set mass to something, just to satisfy LAMMPS

mass 1 26.98

What do I need to add to add Hydrogen atoms to the same box within the lattice?

Thanks!

Hi,

I am looking to add 5-10 individual hydrogen atoms to the lattice of a
slab of aluminum. My code for creating the aluminum is below:

# generate the box and atom positions using a bcc lattice
variable a equal 4.05

boundary p p p

lattice fcc $a
region box block 0 5.0 0 5.0 0 5.0
create_box 1 box
create_atoms 1 box
change_box all triclinic

# Need to set mass to something, just to satisfy LAMMPS
mass 1 26.98

What do I need to add to add Hydrogen atoms to the same box within the
lattice?

a) ​reserve space for two atom types when ​using create_box
b) use create_atoms with the hydrogen atom type (not aluminium)

axel.

Thank Axel, I appreciate your bluntness :).

I made these changes as per the manual. The box now has 2 atom types, and I believe the create_atoms 1 single should create one hydrogen atom between the lattice at (0 0 1.1).

This then runs and produces strains that increment properly and a stress strain curve that is plausible. Upon further inspection the only area that is not making sense is in the log file the length of the box in x, y and z is not going consistently up or down:

Step Lx Ly Lz Press Pxx Pyy Pzz PotEng Temp
0 8.55 8.55 8.55 77510.009 7717.2798 7421.5068 217391.24 -3489.7927 300
1000 8.6473147 8.6473147 8.6473147 7020.3404 -58049.98 -58126.86 137237.86 -5344.5882 3415.8198
2000 8.5768134 8.5768134 8.5768134 873.96864 -24860.661 -24310.971 51793.538 -5531.5891 1045.2548
3000 8.5708485 8.5708485 8.5708485 -269.41797 -20195.656 -20296.155 39683.557 -5487.9815 303.9941
4000 8.574311 8.574311 8.574311 -101.919 -20367.991 -20506.014 40568.248 -5445.5228 40.306878
5000 8.5720713 8.5720713 8.5720713 -192.47552 -19959.513 -20148.708 39530.794 -5490.3785 362.32949
6000 8.566135 8.566135 8.566135 -43.777323 -17486.306 -17778.212 35133.186 -5537.3381 517.68292
7000 8.5656566 8.5656566 8.5656566 -329.76228 -15574.279 -15251.177 29836.169 -5526.6592 318.58832
8000 8.5682278 8.5682278 8.5682278 -195.88622 -14456.345 -14205.151 28073.838 -5500.2921 158.12636
9000 8.5700727 8.5700727 8.5700727 -129.2463 -11750.85 -12540.185 23903.296 -5504.6193 236.67424
10000 8.570373 8.570373 8.570373 -46.138822 -7493.043 -8517.9112 15872.538 -5533.3915 422.03119

I have attached the new input script.

Cheers

Liam

in.liambox.txt (2.3 KB)

Thank Axel, I appreciate your bluntness :).

I made these changes as per the manual. The box now has 2 atom types, and
I believe the create_atoms 1 single should create one hydrogen atom between
the lattice at (0 0 1.1).

This then runs and produces strains that increment properly and a stress
strain curve that is plausible. Upon further inspection the only area that
is not making sense is in the log file the length of the box in x, y and z
is not going consistently up or down:

​i see several significant blunders in your input.

1) you have turned off the check for charge equilibration and i don't see a
fix qeq/reax command anywhere. charge equilibration is a crucial part of
the reaxff model.
2) when adding atoms at possibly high-energy locations, you should run a
minimization first to avoid having your added hydrogen atom being
catapulted through the system
3) you have obviously changed units from metal to real to accommodate the
reaxff parameter set units, but many time unit related input values look
like they were unchanged. so you've implicitly changed the time step from
1fs to 0.001fs, as metal units use ps for time units, while real units have
fs as time unit. considering the low mass of hydrogen, probably a choice of
0.1fs would be more reasonable than 1fs. that will require a 10x increase
of the time steps for the same duration of the respective runs. also check
out time constants for thermostats and barostats as well as the strain rate
parameter and output.

axel.

Hi,

Good catch on the units! Didn’t realize that. I have already written a minimization so I will add that in.

What I am unclear about is charge equilibration. I am see you make checkqeq yes, but then you need cut lo cut hi taper radius values. Are these just 0 and 10 as below or are they specific to the atom in question? Is there a resource that can provide the values needed to be put in the fix qeq/reax section

pair_style	reax/c NULL checkqeq yes
pair_coeff	* * ffield.reax.Fe_O_C_H H Fe
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

Hi,

Good catch on the units! Didn't realize that. I have already written a
minimization so I will add that in.

What I am unclear about is charge equilibration. I am see you make
checkqeq yes, but then you need cut lo cut hi taper radius values. Are
these just 0 and 10 as below or are they specific to the atom in question?
Is there a resource that can provide the values needed to be put in the fix
qeq/reax section

​yes, there is the LAMMPS manual, there are the examples bundled with
LAMMPS, there are the publications about the parameterization and
applications of the corresponding ​force field parameters. at the _very
least_ you should read and carefully study the publication that describes
the parameter set you are using. ReaxFF is not a black box force field. you
have to very carefully check, that you are using it correctly. this is best
done by reproducing published results and not trying to set up something
new.

pair_style reax/c NULL checkqeq yes
pair_coeff * * ffield.reax.Fe_O_C_H H Fe

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

fix 2 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

​this is *definitely* wrong. you don't seem to have spend sufficient time
to even understand the very basic logic of LAMMPS commands.
​why do you specific fix qeq/reax twice?

axel.​