I want to calculate dynamic of CaF_2 bulk system.
The compound has FCC symmetry with Ca atom in (0,0,0) position
and F atoms in (1/4,1/4,1/4) and (-1/4,-1/4,-1/4) positions.
The input file follows bellow. The result looks reasonable, except
atoms close to the boundary jumps from one and of the box to
another. Is it normal or I did something wrong in the input file.
Thank you in advance,
German
# 3d metal
units metal
atom_style hybrid charge
boundary p p p
#lattice sc 5.46
lattice sc 5.4413
region box block 0 1 0 1 0 1
create_box 3 box
create_atoms 1 single 0 0 0
create_atoms 1 single 0.5 0.5 0
create_atoms 1 single 0.5 0 0.5
create_atoms 1 single 0 0.5 0.5
create_atoms 2 single 0.25 0.25 0.25
create_atoms 2 single 0.75 0.75 0.25
create_atoms 2 single 0.75 0.25 0.75
create_atoms 2 single 0.25 0.75 0.75
create_atoms 3 single 0.75 0.75 0.75
create_atoms 3 single 0.75 0.25 0.25
create_atoms 3 single 0.25 0.75 0.25
create_atoms 3 single 0.25 0.25 0.75
mass 1 40.0
mass 2 19.0
mass 3 19.0
set atom 1 charge 2.0
set atom 2 charge 2.0
set atom 3 charge 2.0
set atom 4 charge 2.0
set atom 5 charge -1.0
set atom 6 charge -1.0
set atom 7 charge -1.0
set atom 8 charge -1.0
set atom 9 charge -1.0
set atom 10 charge -1.0
set atom 11 charge -1.0
set atom 12 charge -1.0
1) It does not look like result of periodic boundaries.
In order to demonstrate what is going on I've attached
figure with atoms distribution at zero stem and at step 250.
Some of the atoms on the boundary removed and attached on the other
side out of convenient box.
The velocities are initialized at temperature equals 3 K.
2) How I can initialise charge for the type of atoms?
I've initialised it for each atom in the initial unit cell.
3) I need to calculate heat flux in this system.
the expression for the heat flux is
J = sum_i v_i epsilon_i + 1/2 sum_pairs F_ij r_ij . (v_i+v_j)
This problems already discussed by Philip Howell.
How I can extract pair forces F_{ij}
1) Why is that not PBC?
2) see the set command, or put the atoms in a data file with a charge
on every atom
3) you'd have to write your own fix to compute something like that, and
have it loop over the neighbor list and call pair->single() to get IJ force
I've used compute_group_group as an example for calculation of heat flux.
The modified files: style.h, compute_hf_atom.h and compute_hf_atom.cpp
(see attachments)
The potential energy calculated using force->pair->eatom() is in
agreement with result of
compute_pe_atom.cpp, this no wonder. There as then I try to calculate
potential energy
using return from eng =
pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair), the
result is
absolutely different. Is the value returned by pair->signle() single
pair contribution to
potential energy?
For pair buck/coul/long, it looks like single() should do
the same thing as the per-atom logging. So long as
you recognize that there will be no long-range component
to the energy, just like compute pe/atom. I.e. you can't
compare total energy to the output of "pe" by the thermo
command, which includes the long-range part (for the entire
system).
To debug, I would setup a 2-atom or 3-atom test, see what
compute pe/atom stores and your single() calls, and put
in print statements if they are different.
there last column corresponds to potential energy of each atom
there as pair->signle() single pair contribution to potential energy gives
-0.999227
0.00000
by the way 0.999227/2=0.4996135
I also do not understand why number of neighbors for the first atom in
the sum of pair->singe() equal 1 ,
but for the second atom 0?
This is right for more complicate systems too. I tried LJ fcc with
periodic boundary conditions.
the number of neighbors for different atoms is different. Sum of
energies over atoms in the unit cell
calculated from pair->single and from compute pe/atom is the same, but
contributions from separate atoms
are different.
no, sorry, I don't even understand what you're asking.
The total PE is the sum of the per-atom PE, and you
say that's correct. You talk about pair->single() but
neither of those quantities are computed via pair->single().
So what do you think LAMMPS is not doing correctly?
I would like to insert computation of heat flux to the code.
The expression for heat flux looks like
J = sum_i v_i epsilon_i + 1/2 sum_pairs F_ij r_ij . (v_i+v_j).
You recommended me to use pair->single() function.
In order to check if I use it right I compare potential energy per atom
calculated by
1) compute_PE_atom
and
2) run over list->firstneigh for each atom
and eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair)
The result is different.
As I can understand it because the list of neighbors I use is for i>j
so for the system with four atoms, single fcc cell,
first atom 0 has neighbors 1, 2, 3,
second atom 1 - 2,3,4 and the last one no neighbors.
If I summarize potential energy calculated for all atoms the sum is
the same for
both methods.
Question: there I can take "correct" list of atoms neighbors?
The standard neigh list in LAMMPS only stores the i,j pair once.
It might be with I, it might be with J. The compute pe/atom
uses energy tallied by pair and other styles that takes account
of this when eng is assigned to atoms.
If you use pair->single() with the neigh list, you will have to
worry about that also. E.g. assign half the result to each
atom, and worry about when atom J is on another proc.
