How to calculate the simple Keating valence force field model with LAMMPS

Dear users,

I am trying to set up a simple atomistic Keating valence force field model with LAMMPS for zincblende material heterostructures (like e.g. InAs quantum dots in GaAs). The aim of the Keating model is to compute the change of distance between an atom and its nearest neighbours (4 for zincblende lattice) and the angles of the bonds. Do you, please, know of any procedures in LAMMPS which can do that? Specifically, the equations needed to be computed are in equation (2.47) in the following thesis: https://elib.suub.uni-bremen.de/edocs/00103899-1.pdf (the results of which I actually try to reproduce in the end).

I attach below the LAMMPS infile which I made and where I try to do that.

Can you please help me?

Thank you very much,

Petr

infile:

units metal
atom_style atomic
boundary p p p

variable no_steps equal 2000
variable ba equal 5.65325
variable da equal 6.0583
variable bx equal 200
variable by equal 200
variable bz equal 100
variable lz equal 10
variable dlz equal 20
variable duz equal 40
variable dr equal 100

region box block -\{bx\} {bx} -\{by\} {by} 0 \{bz\} \#region box block 0 {bx} -\{by\} {by} 0 ${bz}
#create_box 3 box
create_box 4 box

region box1 block -\{bx\} {bx} -\{by\} {by} \{lz\} {dlz}
region dot cone z 0 0 \{dr\} 0 {dlz} ${duz}

#region box2 block -\{bx\} {bx} -\{by\} {by} \{lz\} {dlz}
#region dot2 cone z \{da4\} {da4} \{dr\} 0 {dlz} ${duz}

mass 1 69.723
mass 2 74.921595
mass 3 114.818
mass 4 74.921595

region del1 intersect 2 box box1
region del2 intersect 2 box dot

lattice fcc ${ba} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0
create_atoms 1 box

lattice fcc ${ba} origin 0.25 0.25 0.25 orient x 1 0 0 orient y 0 1 0
create_atoms 2 box

delete_atoms region del1
delete_atoms region del2

lattice fcc ${da} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0
create_atoms 3 region box1
create_atoms 3 region dot

lattice fcc ${da} origin 0.25 0.25 0.25 orient x 1 0 0 orient y 0 1 0
create_atoms 4 region box1
create_atoms 4 region dot

region show union 2 box1 dot
group show1 region show

#min_style hftn
min_style cg
min_modify dmax 0.0005
minimize 1e-12 1e-1 5000 5000
#fix 1 all nvt temp 0.01 0.001 0.01
fix 1 all nvt temp 20.01 10.001 0.1 #0.01

compute dsp all displace/atom
dump d all custom ${no_steps} dists.eim id c_dsp[1] c_dsp[2] c_dsp[3] c_dsp[4]

dump id all custom ${no_steps} dump.eim id ix iy iz fx fy fz

dump 2 all image ${no_steps} testim.*.jpg type type &
                         axes yes 1 0.05 view 89 -30 zoom 3
dump_modify 2 region box acolor 4/3/2/1 yellow/blue/green/red
dump_modify 2 region box adiam 3 2
dump_modify 2 region box adiam 4 1 #1.5
dump_modify 2 region box adiam 1 0.2 #5 #0.1
dump_modify 2 region box adiam 2 0.2 #5 #0.1

run ${no_steps} #2000

Dear users,

I am trying to set up a simple atomistic Keating valence force field
model with LAMMPS for zincblende material heterostructures (like e.g.
InAs quantum dots in GaAs). The aim of the Keating model is to compute
the change of distance between an atom and its nearest neighbours (4 for
zincblende lattice) and the angles of the bonds. Do you, please, know of
any procedures in LAMMPS which can do that? Specifically, the equations

that would have to be implemented as a pair style. there are several pair styles in the MANYBODY package that look at neighbors of atoms to compute interactions, but I am not aware of any implementation that is compatible with your needs. You would have to program it (in C++) yourself or perhaps find it implemented by somebody, that has not contributed it to the LAMMPS community.

needed to be computed are in equation (2.47) in the following thesis:
https://elib.suub.uni-bremen.de/edocs/00103899-1.pdf (the results of
which I actually try to reproduce in the end).

I attach below the LAMMPS infile which I made and where I try to do that.

that input does not contain any command enabling computation of forces, so it won’t do anything.
as mentioned before, you need to find or program such a model and activate it via a pair_style command.

Can you please help me?

not much beyond what I wrote above.

axel.

Dear Alex and LAMMPS users,

first of all, thank you very much for your reply. Please see also below my response to your comments and suggestions.

    Dear users,

    I am trying to set up a simple atomistic Keating valence force field
    model with LAMMPS for zincblende material heterostructures (like e.g.
    InAs quantum dots in GaAs). The aim of the Keating model is to compute
    the change of distance between an atom and its nearest neighbours (4
    for
    zincblende lattice) and the angles of the bonds. Do you, please,
    know of
    any procedures in LAMMPS which can do that? Specifically, the equations

that would have to be implemented as a pair style. there are several pair styles in the MANYBODY package that look at neighbors of atoms to compute interactions, but I am not aware of any implementation that is compatible with your needs. You would have to program it (in C++) yourself or perhaps find it implemented by somebody, that has not contributed it to the LAMMPS community.

Yes, in the meantime I realized that I would need some kind of pair_style defined in LAMMPS to have the atoms "move" during minimization. Actually, the Keating potential itself is not really critical for my application for several reasons:

1. for example the Tersoff potential would suit very well too and its even more precise than Keating and I only need to minimize reasonably the energy of the lattice with quantum dot from different zincblende crystal. The problem is that I have not found parameters for other zincblende lattices apart from GaAs that is provided in the examples of LAMMPS. Do you know about a repository of such params for zinclblende crystals, please?

2. since my main aim is the tight-binding code itself, after the minimization of the energy I need to compute the change of bond length relative to that before the minimization and, moreover, to find the change in the directional cosines, i.e., angles of the bonds. This then used to factorize the offdiagonal params in tight binding, see e.g. Eq. (2.51), (2.52), and (2.53) in the thesis https://elib.suub.uni-bremen.de/edocs/00103899-1.pdf I think that can be in some way computed in LAMMPS as I have seen on the program web documentation.

    needed to be computed are in equation (2.47) in the following thesis:
    https://elib.suub.uni-bremen.de/edocs/00103899-1.pdf (the results of
    which I actually try to reproduce in the end).

    I attach below the LAMMPS infile which I made and where I try to do
    that.

that input does not contain *any* command enabling computation of forces, so it won't do anything.
as mentioned before, you need to find or program such a model and activate it via a pair_style command.

I was recently experimenting with the following combination of commands added to my code below. Those were:

pair_style zero 2.5 #5.0
pair_coeff * *
neighbor 2.0 bin
group anion type 2 4
group cation type 1 3
compute gf anion group/group cation pair yes
dump frc all custom ${no_steps} force.eim c_gf

The atom types 2 and 4 are for Arsenic, 1 is Gallium and 3 is Indium. By the neighbor command I was aiming at choosing only the nearest neighbours which are always As for Ga or In atoms and vice versa in zincblende lattice. Am I aiming in the right direction, please? Moreover, I get error with output of the group/group compute command, where is the error please?

    Can you please help me?

not much beyond what I wrote above.

I think you already have, at least by confirming some of my conclusions as seen above, so I thank you very much indeed.

axel.

Petr

  1. for example the Tersoff potential would suit very well too and its
    even more precise than Keating and I only need to minimize reasonably
    the energy of the lattice with quantum dot from different zincblende
    crystal. The problem is that I have not found parameters for other
    zincblende lattices apart from GaAs that is provided in the examples of
    LAMMPS. Do you know about a repository of such params for zinclblende
    crystals, please?

no. you have to search the published literature.

  1. since my main aim is the tight-binding code itself, after the
    minimization of the energy I need to compute the change of bond length
    relative to that before the minimization and, moreover, to find the
    change in the directional cosines, i.e., angles of the bonds. This then
    used to factorize the offdiagonal params in tight binding, see e.g. Eq.
    (2.51), (2.52), and (2.53) in the thesis
    https://elib.suub.uni-bremen.de/edocs/00103899-1.pdf I think that can be
    in some way computed in LAMMPS as I have seen on the program web
    documentation.

if you have seen it, you should try it. i have no time to lookup stuff in big pdfs.

I was recently experimenting with the following combination of commands
added to my code below. Those were:

pair_style zero 2.5 #5.0
pair_coeff * *

pair style zero is a pointless exercise. it does not compute any forces either. please see the documentation.

neighbor 2.0 bin
group anion type 2 4
group cation type 1 3
compute gf anion group/group cation pair yes

you cannot compute pairwise energies and forces from a pair style that does not compute any forces and energies.

dump frc all custom ${no_steps} force.eim c_gf

The atom types 2 and 4 are for Arsenic, 1 is Gallium and 3 is Indium. By
the neighbor command I was aiming at choosing only the nearest
neighbours which are always As for Ga or In atoms and vice versa in
zincblende lattice. Am I aiming in the right direction, please?

hard to say. you seem to be quite far from grasping the basics of force field calculations.
also, it is not clear what this has to do with tight-binding calculations. if your purpose is to parameterize one of those, you seem to be coming from the wrong direction. i would have expected that you would do some DFT calculation (and thus avoid the whole mess of finding pair-wise potentials and parameters.

Moreover, I get error with output of the group/group compute command,

what error? you are not making it easy to help you, by not providing crucial information.

where is the error please?

my guess is that you get ‘Pair style does not support compute group/group’, which is appropriate for pair style zero, although one could argue, that it should have such a function and return force = 0.0 and energy = 0.0

axel.