Dear LAMMPS users,
I am new to LAMMPS and Molecular Dynamics in general, so I am sorry if my questions are very basic.
For my masters thesis I want to study the diffusion of copper atoms on a copper (111) surface.
To start, I have created a Dimer on my copper surface in the hh (hcp-hcp) position. I want to obtain the potential energy of this position and of the Dimer in the ff (fcc-fcc) position, and the height of the energy barrier in between.
To obtain the hight of the energy barrier I tried to use the NEB command. I use a compute pe/atom command to get the potential energy of the Dimer and for every replica I create a dump file with the energy. While the potential energy difference of the start and end postition is more or less what you can find in literature, I get a potential well instead of a potential barrier in between. I guess there is something wrong with the way I think I can get the energy. Can you tell me the best way to obtain the height of the energy barrier when using the NEB command?
Thank you!
Regards,
Elisabeth
NEB is finding the energy as the entire system relaxes
around the dimer moving over the barrier. I'm not
clear on what you expect pe/atom to tell you unless
you are freezing everything but the dimer. NEB prints
out total energies for each replica as it proceeds thru
its 2 stages. Did those energies indicate it found a barrier
and not a well?
Steve
Dear Steve,
thank you for answering so fast.
In literature I find that there should be a potential energy barrier of about 40 meV between
those two positions. If I look at the total energies NEB prints out, I still find a well.
I guess there is some basic error in my code that I am just unable to find.
Here is my code:
boundary p p f
units metal
dimension 3
atom_style atomic
atom_modify map array sort 0 0.0
lattice fcc 3.615 orient x -1 -1 2 orient y 1 -1 0 orient z 1 1 1
region mybox block -3 3 -2 2 -4 4
create_box 1 mybox
create_atoms 1 box
pair_style eam
pair_coeff * * /mylammps/potentials/Cu_u3.eam
region void1 block INF INF INF INF 2.15 INF
region void2 block INF INF INF INF INF -1.85
region void union 2 void1 void2
delete_atoms region void
create_atoms 1 single -0.128837407 0.214056858 2.286883361
create_atoms 1 single 0.224956677 -0.022038463 2.286883361
region nebatoms block -1.0 1.0 -1.0 1.0 1.5 INF
group nebatoms region nebatoms
group nonneb subtract all nebatoms
fix 1 all nve
fix 2 nebatoms neb 10.0
fix 3 nonneb setforce 0.0 0.0 0.0
thermo 100
variable u uloop 20
dump 1 nebatoms atom 10 ../../neb/run1/eam/dump.neb.$u
dump 2 nonneb atom 10 ../../neb/run1/eam/dump.nonneb.$u
min_style quickmin
neb 1.0e-8 10000 10000 100 coord/2.coord
And the coordinates in the file 2.coord are:
1537 -1.31163 0.140507 14.3078
1538 0.776409 -1.06668 14.3077
Thanks for helping me!
Best,
Elisabeth
Don't know - I guess I would start with simpler problems
and see if you can get agreement. There are also lots
of EAM potentials.
Steve
Hey,
after I made an 'svn update' today, I get a potential barrier and the correct values for my code.
Thought you might want to know.
Best,
Elisabeth