# Did my Lennard-Jones potential LAMMPS simulation of a metal island on an FCC crystal surface find an energy minimum? (don't understand output)

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

orthogonal box = (-0.17677670 -0.10206207 -4.4907312) to (13.823223 12.022294 1.8177570)
1 by 1 by 1 MPI processor grid
1237 atoms

# ------------- 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

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

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

# ------------- 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

• Minimization is only for molecular statics, at each time step atoms position are changed to minimize over all energy of system. In molecular dynamics (greater than 0K), you don't do minimization. It is the ensemble average that system attains i.e equilibrium. It is very bad practice to add velocity and minimize. Think for while over this. your atoms are moving at some velocity and you are also trying to move again by minimization. Jul 22, 2022 at 17:16
• @Pranavkumar "Minimization is only for molecular statics" I disagree, finding lowest energy configurations is a fundamental approach in physics used widely. "It is very bad practice to add velocity and minimize" can you expand on that? Why can't a surface structure have a minimum energy configuration at a finite temperature? Isn't that in fact what happens in reality? Can't LAMMPS be used to model real situations?
– uhoh
Jul 22, 2022 at 22:17
• Minimize first followed by Use ensemble such as fix NPT, NVT to control on temperature. Let system go through ensemble average to be in equilibrium. Jul 23, 2022 at 2:51
• @Pranavkumar okay that's actionable help, thank you!
– uhoh
Jul 23, 2022 at 3:18
• Did you get an energy minimum? Do you need @Pranavkumar to turn those very useful comments into an answer? Jan 31 at 21:55