read_data and tersoff potential

Dear Users/Developers

I am a newbie in LAMMPS, and I am stuck at one point. I searched the
documentation and the mailing list to no avail.

I would like to read a geometry from a file, and use the Tersoff
potential to do my calculations

When I use the prescribed "create_atom" route (i.e. in Tutorial 1,2),
an energy is calculated
however, when I try to input the system I actually want to calculate
using read_data, the calculated energy is 0, it is as if no atom is
interacting with each other.

The data file works with a Dreiding potential

I tried "full", "molecule", "atom" style formatted data files. Nothing
changes. Energy is still zero

Is there something special I need to do using a bond order potential
with read_data?

Please find the tests I've performed attached.

Best,
Baris

-Test 1- Fails, Energy is always 0
Input:
# ---------- Initialize Simulation ---------------------
clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array

# ---------- Create Atoms ---------------------
read_data rod001-7mrh8_5-relaxed-100nm-tersoff.lammps05
# ---------- Define Interatomic Potential ---------------------
pair_style tersoff
pair_coeff * * SiH-obm.tersoff Si H
mass 1 23.74
neighbor 2.0 bin
neigh_modify delay 10 check yes

# ---------- Define Settings ---------------------
compute eng all pe/atom
compute eatoms all reduce sum c_eng

# ---------- Run Minimization ---------------------
reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms
min_style cg
minimize 1e-25 1e-25 5000 10000

variable natoms equal "count(all)"
variable teng equal "c_eatoms"
variable length equal "lx"
variable ecoh equal "v_teng/v_natoms"

print "Total energy (eV) = \{teng\};" print "Number of atoms = {natoms};"
print "Lattice constant (Angstoms) = \{length\};" print "Cohesive energy \(eV\) = {ecoh};"

print "All done!"

Output:
...
Total energy (eV) = 0;
Number of atoms = 8280;
Lattice constant (Angstoms) = 50;
Cohesive energy (eV) = 0;
All done!

-----Test 2 -- Unphysical system to see if create atoms route works,
at least it calculates energy
# ---------- Initialize Simulation ---------------------
clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array

# ---------- Create Atoms ---------------------
lattice fcc 4
region box block 0 1 0 1 0 1 units lattice
create_box 1 box

lattice fcc 4 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 box
replicate 1 1 1

# ---------- Define Interatomic Potential ---------------------
pair_style tersoff
pair_coeff * * SiH-obm.tersoff Si
mass 1 23.74
neighbor 2.0 bin
neigh_modify delay 10 check yes

# ---------- Define Settings ---------------------
compute eng all pe/atom
compute eatoms all reduce sum c_eng

# ---------- Run Minimization ---------------------
reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms
min_style cg
minimize 1e-25 1e-25 5000 10000

variable natoms equal "count(all)"
variable teng equal "c_eatoms"
variable length equal "lx"
variable ecoh equal "v_teng/v_natoms"

print "Total energy (eV) = \{teng\};" print "Number of atoms = {natoms};"
print "Lattice constant (Angstoms) = \{length\};" print "Cohesive energy \(eV\) = {ecoh};"

print "All done!"

-Output-----
Total energy (eV) = -15.467621918130557646;
Number of atoms = 4;
Lattice constant (Angstoms) = 3.8975156250000040359;
Cohesive energy (eV) = -3.8669054795326394114;
All done!

Dear Users/Developers

I am a newbie in LAMMPS, and I am stuck at one point. I searched the
documentation and the mailing list to no avail.

I would like to read a geometry from a file, and use the Tersoff
potential to do my calculations

When I use the prescribed "create_atom" route (i.e. in Tutorial 1,2),
an energy is calculated
however, when I try to input the system I actually want to calculate
using read_data, the calculated energy is 0, it is as if no atom is
interacting with each other.

The data file works with a Dreiding potential

*yes* and that is the hint why you get what you get. your data file
has bonds (and angles etc.) defined, but the tersoff potential does
not use bonds. when you do have bonds define, the neighborlist of
atoms is subject to "non-bonded" exclusions according to the settings
of the "special_bonds" command. normally that would be 0.0 for bonds
and 0.0 for angles and 1.0 for dihedrals. the consequence is that such
atom pairs that have a 0.0 exclusion factor are pruned from the
neighbor list, and thus don't show in the tersoff calculation, even
though they should.

I tried "full", "molecule", "atom" style formatted data files. Nothing
changes. Energy is still zero

Is there something special I need to do using a bond order potential
with read_data?

yes. get rid of the bonds/angles/etc. or set special_bonds lj/coul 1.0 1.0 1.0.

axel.

Dear Axel,

My file for "atomic" style does not have anything but the atom positions. If I
remember correctly, the program forced me to remove everything else. The file
contains the minimum necessities as instructed by
http://lammps.sandia.gov/doc/read_data.html

Please find the beginning of the file attached.

With your suggestion, I tried adding special_bonds lj/coul 1.0 1.0
1.0, to the previously mentioned test file. This doesn't have any
effect as well. The calculated energy is still zero.

Is there a problem with my compilation?

C.
Baris

   8280 atoms
   2 atom types

     0.000000000 50.000000000 xlo xhi
     0.000000000 50.000000000 ylo yhi
     0.000000000 1002.228700000 zlo zhi

Masses

   1 28.086000 #Si
   2 1.007970 #H_

Atoms

      1 1 32.226430169 26.467452091 4.073018235
      2 1 26.467451546 32.226430833 4.073018235
      3 1 20.776541513 32.226369363 1.373746281
....

Impossible to say anything with such incomplete information, which in addition contradicts the description in your previous mail.

provide a complete(!) and small test case or sort it out yourself.

Axel

Dear Axel,

You are right, maybe there is more than one thing that is resulting
this error. I used the "test1" input file attached in the first
communication with an atomic style data file for a silicon cluster and
the tersoff potential provided with the lammps distribution. It didn't
work. In case this might be due to another problem I am attaching the
exact input file (which is the test 1 of previous communication,
modified with your suggested special_bonds command), the tersoff
potential, and the data file I use.

Best,
Baris

Test 1:

# ---------- Initialize Simulation ---------------------
clear
units metal
dimension 3
boundary p p p
atom_style atomic

# ---------- Create Atoms ---------------------
read_data rod001-7mrh8_5-relaxed-100nm-tersoff.lammps05
# ---------- Define Interatomic Potential ---------------------
special_bonds lj/coul 1.0 1.0 1.0
pair_style tersoff
pair_coeff * * SiH-obm.tersoff Si H
mass 1 23.74
neighbor 2.0 bin
neigh_modify delay 10 check yes

# ---------- Define Settings ---------------------
compute eng all pe/atom
compute eatoms all reduce sum c_eng

# ---------- Run Minimization ---------------------
reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms
min_style cg
minimize 1e-25 1e-25 5000 10000

variable natoms equal "count(all)"
variable teng equal "c_eatoms"
variable length equal "lx"
variable ecoh equal "v_teng/v_natoms"

print "Total energy (eV) = \{teng\};" print "Number of atoms = {natoms};"
print "Lattice constant (Angstoms) = \{length\};" print "Cohesive energy \(eV\) = {ecoh};"

print "All done!"

rod001-7mrh8_5-relaxed-100nm-tersoff.lammps05.gz (68.1 KB)

SiH-obm.tersoff.gz (732 Bytes)

Dear Axel,

You are right, maybe there is more than one thing that is resulting
this error. I used the "test1" input file attached in the first
communication with an atomic style data file for a silicon cluster and
the tersoff potential provided with the lammps distribution. It didn't
work. In case this might be due to another problem I am attaching the
exact input file (which is the test 1 of previous communication,
modified with your suggested special_bonds command), the tersoff
potential, and the data file I use.

the explanation is simple: your atom positions are too far apart.
tersoff models are fairly shortranged. could it be that you have
confused atomic units with angstroms somewhere?

i just scaled down the whole thing by a factor of 0.6 and it ran fine.
even with the zero energy, LAMMPS was just doing what you asked it to
do.

change_box all x scale 0.6 y scale 0.6 z scale 0.6 remap

apart from that, there are a number of issues in your input

- you needlessly compute the total energy via a compute reduce. just
use the "pe" thermo keyword.
- what is the point of using fix box/relax while forcing it to update
the box isotropically.?
- why do a minimization on such a huge, yet perfectly symmetrically
replicated system? why no minimize a small subunit and then use the
replicate command?

axel.

Dear Axel,

Thank you very much. I now see clearly that there is something wrong
with the units. Apparently I was not careful enough about the
documentation of the read_data and using the msi2lmp converter.

Thank you very much for the pointers. I will surely be mindful of them
in the future. The script you are seeing is not something I use for my
calculations, but something directly copied and pasted from the lammps
tutorial I have found browsing the lammps website, in order to
minimize a possible error I might be doing with my own scripts
(https://icme.hpc.msstate.edu/mediawiki/index.php/LAMMPS_Help2). Since
the Dreiding uses bonds etc, it was compensating for the error, and I
missed the strangeness of the units.

My current goal is: I am trying to have an idea about the scale of
impact of self interaction arising from fast moving acoustic phonons
due to periodic boundary conditions on the heat current auto
correlation function. After that I'll try and see if I can do a
rolling sampling or a model fit to remove any pbc related unfeasible
features. Length is an important factor for my quest, and there are
some supercell effects which I would like to be sure they are
included.

Best,
Baris