BN Potential in LAMMPS

Hi Lammps users, I am trying to parametrize Tersoff potential for BN using the values from

values are from Moon et.al Physica B 336 (2003) 329-334 as below, however, when I run my simulation using the potential values, I do not even get close to the minimum energy graph. I have included my input file as well. Anyone with some insight on what I maybe doing wrong, I would appreciate your input. Thanks

Jean Njoroge

these entries are in LAMMPS “metal” units:

A,B = eV; lambda1,lambda2,lambda3 = 1/Angstroms; R,D = Angstroms

format of a single entry (one or more lines):

element 1, element 2, element 3,

m, gamma, lambda3, c, d, costheta0,

n, beta, lambda2, B, R, D, lambda1, A

B B B 3.0 1.0 0.0 0.5269 0.001587 0.5
3.99290 0.000001601 1.38190 171.290 1.95 0.15 1.728 223.10

B B N 3.0 1.0 0.0 0.5269 0.001587 0.5
0.00000 0.000000 2.040950 0.000 1.95 0.15 3.582 0000.0

B N B 3.0 1.0 0.0 0.5269 0.001587 0.5
0.00000 0.00000 1.3819 0.0000 1.95 0.15 1.728 0000.0

B N N 3.0 1.0 0.0 0.5269 0.001587 0.5
1.330 0.0052938 2.040950 209.3237 1.95 0.15 3.582 1191.704

N B B 3.0 1.0 0.0 20312.0 25.510 -.562
1.330 0.0052938 2.040950 209.323789 1.95 0.150 3.582 1191.704

N B N 3.0 1.0 0.0 20312.0 25.510 -.562
0.00000 0.0000 2.70000 000.00 1.95 0.15 5.437 0000.0

N N B 3.0 1.0 0.0 20312.0 25.510 -.562
0.00000 0.0000000 2.040950 000.00 1.95 0.15 3.582 0000.0

N N N 3.0 1.0 0.0 20312.0 25.510 -.562
1.330 0.005293800 2.7000 511.760 1.95 0.15 5.437 6368.140

log log.cBN
#log file name

shell rm out.data
variable latt0 equal 3.49769 #Lattice parameter of Boron Nitride “Reference lattice constant 3.609 A” vary lattice betwwen
variable max_iter equal 1001
variable shift_i equal ceil({max_iter}/2.0) label loop variable i loop {max_iter}
if “i > {max_iter}” then “jump in.equillatt_bulk_mod”
variable j equal i-{shift_i}
variable latt equal ${latt0}+$j/1000.0 #Vary lattice constant in increments on 0.001

#-------------INITIALIZE SIMULATION---------------------------------------
clear
units metal
dimension 3
atom_style atomic
atom_modify map array
boundary p p p
newton on

#number of cpus in x y z dimensions
#processors 1 1 1

#-------------CREAT ATOMS-------------------------------------------------
lattice custom ${latt} a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 &
basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 &
basis 0.25 0.25 0.25 basis 0.75 0.75 0.25 &
basis 0.75 0.25 0.75 basis 0.25 0.75 0.75
region box block 0 5 0 5 0 5 #setup 1 region
create_box 2 box #create box with 1 region and 2 atoms
create_atoms 1 region box basis 5 2 basis 6 2 basis 7 2 basis 8 2

#-------------DEFINE INTERATOMIC POTENTIAL--------------------------------------------

mass 1 10.80
mass 2 14.00

include potential.mod

#------------DEFING AUXILLARY VARIABLES TO CONTAIN CHOESIVE ENERGY AND OUTPUT DATA--------------
variable EPOT equal pe
variable NP equal atoms
variable ELAT equal pe/{NP} variable VOL equal vol variable rho equal {NP}/{VOL} variable AtomVol equal {VOL}/${NP}

min_style cg
minimize 1e-25 1e-25 5000 10000

shell echo {latt} {rho} {AtomVol} {EPOT} {ELAT} {VOL} >> out.data
next i
jump in.equillatt_bulk_mod loop
label break
variable i delete

There is a factor X of the value 0.707. You need to scale A and B for
B-N pairs by this factor.

Ray

Ryan,
Thanks for the response. scale A & B using 0.707, even thought lammps documentation states only B? Perhaps I don’t understand, do you mind expounding on it.

Thanks
Jean

There is a factor X of the value 0.707. You need to scale A and B for
B-N pairs by this factor.

Ray

Perhaps, but you have to look at the paper closely to be sure. Also I
notice the following, " To keep N atoms stable inside the crystalline
structure, the value of wN2N is set to 0", which is not what you have
in your potential file.

Ray

Perhaps, but you have to look at the paper closely to be sure. Also I
notice the following, " To keep N atoms stable inside the crystalline
structure, the value of wN2N is set to 0", which is not what you have

That is X_{N-N}, from the paper you provided.

Ray

Ray,
I have looked at the paper, and I may be missing something, hence reason I am asking in this group. This is my first time trying to parametrize a potential in lammps. First you said to scale A & B using the scaling factor what did you mean? i.e Aij = (AiAj)^1/2 to multiply by the scaling factor as well?
Second, thanks for pointing out the part of X_{N-N}, the two body interactions in this calculations is wrong? Sorry, but I am getting confused by your answers?

N N N 3.0 1.0 0.0 20312.0 25.510 -.562

1.330 0.005293800 2.7000 511.760 1.95 0.15 5.437 6368.140

Perhaps, but you have to look at the paper closely to be sure. Also I
notice the following, " To keep N atoms stable inside the crystalline
structure, the value of wN2N is set to 0", which is not what you have

That is X_{N-N}, from the paper you provided.

Ray

Read the paper carefully to determine if A&B are both scaled. If
scaled, multiply A/B_{ij} by the factor.

X_{N-N} is zero meaning A/B values for the "N N N" entry should be zero.

Ray