An anomalous result about using LAMMPS for the Tersoff potential.

LAMMPS LAMMPS Mailing List Mirror
1

Dear all LAMMPS users & developers.

First of all, I wish to send a gratitude to developers for sharing this nice molecular dynamics

program.

I like to discuss about an anomalous result that I found recently.

I typically do ‘bulk’ things, especially crystalline or liquid Si with the Tersoff potential

by using LAMMPS.

There is no problem at all about bulk things, but recently I found something weird about

Si nano-wire simulations (of course, with the Tersoff potential.).

What I found is that some atoms nearby in the surface of a nano-wire have more stable

potential energy values than the energy value of the ground-state of cubic diamond.

A value is little different with respect of a postion, but the difference is quite huge, for example,

the energy value of cubic diamond; E_cd = -4.63 eV, but some atoms in the nano-wire

has -4.71 eV.

At first, I thought this is due to a characteristic of the potential so I made a custom code to

analyze this phenomenon but only found the discrepancy between my code and LAMMPS

(a result of my code does not show such phenomenon.).

A next thing that I did is checking energy of a single isolated tetrahedron by LAMMPS

(molecular statics). Bond-angles and pair lengths are same as cubic diamond,

therefore atoms in the center of tetrahedron must have an energy value which is identical to

the value the cubic diamond. Nonetheless, it turned out that the center atom has a one-step

lower energy (-4.96 eV).

With the indentical LAMMPS scripts, the Stillinger-Weber and MEAM potential do not show

such anomaly.

I tested in two different LAMMPS version: 31Mar17 and 4Jan19.

The Tersoff potential that I used is Si© type in Physical Review B, vol 38, number 14 (1988).

I included my graphical summary pdf file and LAMMPS scripts with this e-mail.

Any comment about this problem will be helpful for me, personally I looked through

‘Tersoff.cpp & Tersoff.h’. The code itself seems fine but I suspect ev_tally() might be

the reason of this problem.

Si_tetrahedron_isolated.txt (341 Bytes)

Si89.tersoff (1.84 KB)

test.txt (1.57 KB)

anomaly.pdf (204 KB)

2

Hello,

I think the problem is related to what you have pointed out, that is, how the potential energy is partitioned among the atoms. The total energy of the system can be written as

U = \sum_i U_i,
where

U_i = 0.5 * \sum_{j != i} U_ij,
and

U_ij = [f_R(r_ij) - b_ij * f_A(r_ij)] * f_C(r_ij).
However, in LAMMPS, it seems that the total potential energy is not decomposed in this way. For example, it will attribute the part of U_ij that depends on a third atom k equally to atoms i, j, and k.

That said, the total potential energy, U, should not depend on the way of decomposition. Perhaps you can check if this is the case.

Usually, it is not very meaningful to consider the potential energy of a single atom. Some space-average should be able to eliminate the discrepancy between different decomposition schemes.

Bruce

3

yes, each 2-body interaction in Tersoff is assigned 1/2 to each of the I,J atoms

.Each 3-body interaction is assigned 1/3 to each of the I,J,K atoms.
Same for energy and virial (from the force on the 2 or 3 atoms).

Steve

4

Thanks for the kind comments, Bruce and Steve.

Definitely my custom code do ‘energy decomposition’ as the way only considering pair interaction participants; say ith & jth atoms. But considering only the enrgy part, ‘pair_tersoff.cpp’ is also only considering ith & jth atoms.

What I am not sure about is although the centered atom in the isolated tetrahedron & atoms in the nano-wire see the geometrically identical neighboring atoms like the case of cubic diamond, energy values are quite different with the value in cubic diamond.

I think that no matter which decomposition method is used, at least if an environment which an atom see is identical, the energy value of the atom must be same.

Therefore from a pair interaction,

let f_c - cutoff function, f_r - repulsive energy function, f_a - attractive energy function

V_ij = 0.5f_c(rij)[f_r(rij) - b_ij*f_a(rij)],

because in pair_tersoff.cpp, ev_tally() is only dedicating function for energy decomposition,

to ith atom (the centered atom in the isolated tetrahedron.), from 4 pairs, 40.5V_ij should be always stand for.

I also checked ‘pair_sw.cpp’, basically the general code structure is very similar to ‘pair_tersoff.cpp’, however simulations with the Stillinger-Weber potential do not show anomalies in the Tersoff case.

Jo.

2019년 5월 10일 (금) 오후 11:31, Steve Plimpton <[email protected]>님이 작성:

5

Ah, sorry for misunderstanding. Bruce and Steve.

I finally understand the reason of descrepancy, what a stupidity!

Basically, it is due to decomposition method, my code only do a loop ‘V_ij’ and not count ‘V_ji’.

And in the tetrahedral 5 atoms cluster, when i index stands for the centered atom, V_ij is totally identical with the bulk situation, but V_ji is quite different because on the behalf of jth atom, the only neighboring atom is the centered atom (ith atom) therefore b_ji = 1.0.

Consequently V_i = 4*(V_ij + V_ji) and the value V_ji is lower than V_ij because b_ji = 1.0.

2019년 5월 11일 (토) 오후 4:12, ‍조유환[학생](공과대학 기계공학과) <mangoman@…5003…>님이 작성:

6

I also want to correct my previous statement:
"
For example, it will attribute the part of U_ij that depends on a third atom k equally to atoms i, j, and k.
"

I think For Tersoff, the decomposition is only for pairs. The neighbor-atom k just affects the “environment” of atom i. For SW, the decomposition should be based on triplets.

Using your notation:
V_ij = 0.5f_c(rij)[f_r(rij) - b_ij*f_a(rij)],
The total energy V of the system should read (as written in Tersoff’s papers)
V = \sum_i V_i
V_i = \sum_{j != i} V_ij

Suppose the center atom in the tetrahedron is labeled as 1 and the others from 2 to 5, then the energy for atom 1 should read
V_1 = V_12 + V_13 + V14 + V15, (This gives -4.63 eV)
and the energies for the other atoms should be
V_i = V_i1 (i = 2, 3, 4, 5).

In LAMMPS, it might be calculated as
V_1 = (V_12 + V_21)/2 + (V_13+V_31)/2 + (V14 + V41)/2 + (V15 +V51)/2, (Smaller than -4.63 eV)

V_i = (V_i1 + V_1i)/2 (i = 2, 3, 4, 5).

That is, it seems in LAMMPS, the site potential is defined as
V_i = \sum_{j != i} (V_ij + V_ji)/2.

I prefer the original definition
V_i = \sum_{j != i} V_ij,
which can also be written as (as written in Brenner’s papers)
V_i = \sum_{j > i} (V_ij + V_ji) .

But the total energy is not affected by the decomposition method.

Bruce