I'm starting to learn LAMMPs with the intention of replacing my simple Python model for equilibrium configurations of adatoms on a crystal surface. Here's an example of what I do now, using a Lennard-Jones or harmonic oscillator potential between bonded atoms, and a static hexagonal potential imitating the imposed periodic potential of a surface, and I currently use just a simple damped ODE solver to find equilibrium positions at zero temperature:
Since the lattice constant of the adatoms is about 20% larger than that of the substrate potential, they tend to seek various rotational angles where higher order coincidences can be found.
I've built a chunk of fcc crystal with the top face as the (111) surface and a nearest-neighbor distance of 1.0, using lj units and the lj/cut potential.
I set the sigma values at 21/6 for the substrate atoms and about 20% larger for the adatoms, and I understand that LAMMPS will us a geometrical average for the 1 2 pairings.
I've applied periodic boundary conditions in all three directions.
The .in file is below as well as the first bit of the .input file. The full input file with all the initial atom positions is long so I've put it here
Question: What happened when I ran the energy minimization at 10 degrees Kelvin? I understand that this is not an accurate model yet, I can see that the substrate atoms expanded vertically a little bit and the periodic boundary conditions prevented it from doing so laterally, so I need to tune my sigma or my atom spacing to a more relaxed situation.
But because I am new to LAMMPS it's hard to understand if this has really reached a minimum energy configuration. I don't see any big motion in the adatoms, which is the thing I'd like to study. Perhaps I need to fix the positions of the substrate atoms using something like fix freeze substrate 0 0 0 or similar? Would that help the minimizer focus more on the adatoms?
GIF generated by Ovito looking at the .dump file https://pastebin.com/rurbzFi9 I haven't been able to get Ovito to draw the atoms with smaller size yet.
log.lammps
LAMMPS (29 Sep 2021)
# ------------- Initialization
units lj
dimension 3
atom_style atomic
pair_style lj/cut 2.5
boundary p p p
read_data fccbox_v02_input.data
Reading data file ...
orthogonal box = (-0.17677670 -0.10206207 -4.4907312) to (13.823223 12.022294 1.8177570)
1 by 1 by 1 MPI processor grid
reading atoms ...
1237 atoms
read_data CPU = 0.003 seconds
# ------------- System definition
# region myregion block -0.1767767 13.8232233 -0.10206207 12.02229358 -4.4907312 2.25
# create_box 2 myregion
# ------------- Simulation settings
# mass 1 1
# mass 2 1
# pair_coeff 1 1 1.0 1.0 1.0 # these are in input.data
# pair_coeff 2 2 1.0 1.0 1.21
neighbor 0.3 bin # This is an extra skin for neighbors units of sigma
dump fundump all atom 100 fundump.dump
velocity all create 10.0 42 rot no dist gaussian
# ------------- Run
thermo 10
minimize 1.0e-4 1.0e-6 1000 10000
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (src/min.cpp:188)
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 1.6602422
ghost atom cutoff = 1.6602422
binsize = 0.83012109, bins = 17 15 8
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.338 | 5.338 | 5.338 Mbytes
Step Temp E_pair E_mol TotEng Press
0 10 -0.34005694 0 14.647817 59.679297
6 10 -2.5987468 0 12.389127 42.827225
Loop time of 0.0170679 on 1 procs for 6 steps with 1237 atoms
97.5% CPU use with 1 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
-0.340056939447598 -2.59874682375032 -2.59874682375032
Force two-norm initial, final = 1156.0012 320.30496
Force max component initial, final = 64.284628 33.671989
Final line search alpha, max atom move = 1.2865070e-10 4.3319251e-09
Iterations, force evaluations = 6 57
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.012193 | 0.012193 | 0.012193 | 0.0 | 71.44
Neigh | 0.0014319 | 0.0014319 | 0.0014319 | 0.0 | 8.39
Comm | 0.00075269 | 0.00075269 | 0.00075269 | 0.0 | 4.41
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.002691 | | | 15.76
Nlocal: 1237.00 ave 1237 max 1237 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1619.00 ave 1619 max 1619 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 7441.00 ave 7441 max 7441 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 7441
Ave neighs/atom = 6.0153597
Neighbor list builds = 2
Dangerous builds = 0
Total wall time: 0:00:00
.in
# ------------- Initialization
units lj
dimension 3
atom_style atomic
pair_style lj/cut 2.5
boundary p p p
read_data fccbox_v02_input.data
# ------------- System definition
# region myregion block -0.1767767 13.8232233 -0.10206207 12.02229358 -4.4907312 2.25
# create_box 2 myregion
# ------------- Simulation settings
# mass 1 1
# mass 2 1
# pair_coeff 1 1 1.0 1.0 1.0 # these are in input.data
# pair_coeff 2 2 1.0 1.0 1.21
neighbor 0.3 bin # This is an extra skin for neighbors units of sigma
dump fundump all atom 100 fundump.dump
velocity all create 10.0 42 rot no dist gaussian
# ------------- Run
thermo 10
minimize 1.0e-4 1.0e-6 1000 10000
beginning of.data
Data File created by fcc_box_v02.py
1237 atoms
2 atom types
-0.17677669529663687 13.823223304703363 xlo xhi
-0.10206207261596577 12.022293580366178 ylo yhi
-4.4907311951024935 1.8177570093457944 zlo zhi
Masses
1 1.0
2 1.9202965708989803
Pair Coeffs
1 1 1.0 1.122462048309373
2 2 1.0 1.360242170692667
Atoms
1 1 0.0 0.0 0.0
2 1 0.0 1.7320508075688774 0.0
3 1 0.0 3.464101615137755 0.0
4 1 0.0 5.196152422706632 0.0
5 1 0.0 6.92820323027551 0.0
6 1 0.0 8.660254037844387 0.0
7 1 0.0 10.392304845413264 0.0
8 1 1.0 0.0 0.0
9 1 1.0 1.7320508075688774 0.0
