Hello,
I am unclear how bonding works in periodic systems with lammps. I am
running calculations on a diamond lattice to determine the elastic
tensor following the ELASTIC example, but with pairwise morse functions
and also with bonding harmonic functions instead of the SW default
potential. My tensors for the same system with pairwise morse functions
and bonding harmonic functions are very different, with different
constants zero and non-zero, which is not what I would expect. I fear
that the bonds that are supposed to connect across a cell boundary to an
image are instead connecting across the whole length of the cell. Could
someone provide me with any insight into how this works? I have set my
cutoff distance to be ~10% larger than the equilibrium bond length.
Thanks a lot,
Michal
Hello,
I am unclear how bonding works in periodic systems with lammps. I am
running calculations on a diamond lattice to determine the elastic
tensor following the ELASTIC example, but with pairwise morse functions
and also with bonding harmonic functions instead of the SW default
potential. My tensors for the same system with pairwise morse functions
and bonding harmonic functions are very different, with different
constants zero and non-zero, which is not what I would expect. I fear
that the bonds that are supposed to connect across a cell boundary to an
image are instead connecting across the whole length of the cell. Could
someone provide me with any insight into how this works? I have set my
cutoff distance to be ~10% larger than the equilibrium bond length.
you should probably also set a suitable communicate cutoff.
if the system is set up properly, it should look as if you have bonds
across the box, but the neighborlist build will search for the closest
neighbor. having a sufficiently large ghost atom area is a
requirement.
i suggest you set up a system with only a few atoms and make some
tests that you can verify by hand and see, if you can make it work.
that is what i would do.
axel.
I have increased my cut-off distance to longer than the cell diagonal (50 in a 3.6^3 cell), and it made no difference to my results, while adding a lot to the runtime. And I am using only an 8 atom cell that I have solved previously, so I am confident of my input.
I spent the last six months writing my own code for this system to calculate these exact same properties before discovering lammps, so I am fairly familiar with the overall problem.
In the bonding system, I found only a non-zero value for C66, resistance to shear, while in the pairwise equivalent system I got non-zero values for C11,22,33,12,13,23,44,55 and 66. Since one is a using a morse potential while the other a harmonic potential, I would expect the values to be quantitatively different but qualitatively similar. (I am considering values of 10^-5 GPa and lower to be effectively zero, in the pairwise calculation I get values on the order or 10^0 - 10^2 GPa)
Thanks,
Michal
If you haven't yet, check to make sure that these bonds (that cross
periodic boundaries) are present in your data file. One way to do
this is to visualize the system in VMD/topotools. (I attached some
instructions for this to this post. If you already know how to do
this, then please excuse me.)
Anyway, you can visualize multiple periodic copies of the system by
selecting the "Graphics"->"Representations" menu option, and clicking
on the "Periodic" tab (middle-right), and checking the +x, +y, +z
boxes. If you are missing bonds across the boundary, you should
notice their absence when you view things this way.
Cheers!
Andrew
README_visualize.txt (2.85 KB)
Hello Andrew,
Thank you for the suggestions. I know I am listing the bonding
properly, but perhaps I am using too small a system and it thinks
they're in cell bonds? I am only using 8 atoms.
Also, when I visualize my input data in VMD, everything appears to be
zero in the z axis, when that is not what my input states. Could you
please have a look at my input data file and let me know if you can see
any issues with it?
Thanks,
Michal
1 Comments
2
3 8 atoms
4 16 bonds
5 0 angles
6 0 dihedrals
7 0 impropers
8
9 1 atom types
10 1 bond types
11 0 angle types
12 0 dihedral types
13 0 improper types
14
15 0 3.5699999999999998 xlo xhi
16 0 3.5699999999999998 ylo yhi
17 0 3.5699999999999998 ylo yhi
18 2.1859945364780254e-16 2.1859945364780254e-16 2.1859945364780254e-16
xy xz yz
19
20 Atoms
21
22 1 1 1 0 0 0
23 2 1 1 2.18599454e-16 1.785 1.785
24 3 1 1 1.785 1.785 0
25 4 1 1 1.785 1.09299727e-16 1.785
26 5 1 1 2.6775 0.8925 2.6775
27 6 1 1 0.8925 0.8925 0.8925
28 7 1 1 0.8925 2.6775 2.6775
29 8 1 1 2.6775 2.6775 0.8925
30
31 Bonds
32
33 1 1 1 5
34 2 1 1 6
35 3 1 1 7
36 4 1 1 8
37 5 1 2 5
38 6 1 2 6
39 7 1 2 7
40 8 1 2 8
41 9 1 3 5
42 10 1 3 6
43 11 1 3 7
44 12 1 3 8
45 13 1 4 5
46 14 1 4 6
47 15 1 4 7
48 16 1 4 8
To follow up on my last email, the issue with everything being flat was
that I used the wrong atom style in the tk console. However, all the
bonds that are supposed to cross the cell boundary are just being bonded
within the same cell. Is there a special format for specifying in-cell
vs out of cell bonds?
Thanks,
Michal
this data file has an error. it contains the line ylo yhi twice.
axel.
Ok, thanks, I can't believe I missed that...
With that fix and removing the 10^-16 tilt I had on the system, the
simulation appears to be working now.
I just have one last question, which I haven't been able to figure out
from the documentation. I want to use neigh_modify exclude molecule
group-ID to speed up my calculation, but I can't figure out what the
Group-ID for my atoms is. Could you let me know where I could find this
or what the default value would be? I'm using atom_style molecular.
More generally, how can I turn off 1-3 and 1-4 neighbour interactions
while maintaining 1-2 interactions?
Thanks a lot,
Michal
Ok, thanks, I can't believe I missed that...
With that fix and removing the 10^-16 tilt I had on the system, the
simulation appears to be working now.
I just have one last question, which I haven't been able to figure out
from the documentation. I want to use neigh_modify exclude molecule
group-ID to speed up my calculation, but I can't figure out what the
you need to speed up a calculation with 8 atoms???
Group-ID for my atoms is. Could you let me know where I could find this
http://lammps.sandia.gov/doc/group.html
or what the default value would be? I'm using atom_style molecular.
there is no default.
More generally, how can I turn off 1-3 and 1-4 neighbour interactions
while maintaining 1-2 interactions?
exclusion factors can be any value between 0 and 1; but your question
makes no sense,
based on the information you provided. please explain in more detail
what you are referring to and what your specific problem is.
axel.
When I increased my communicate cutoff value, the calculation slowed down to the point that every minimization took a second or two, and it was listing large numbers of 1-3 and even more 1-4 interactions per atom, even though I hadn't specified a potential between them. I only have 4 1-2 interactions per atom, and that's the bonded harmonic potential. I suspected that the additional calculation time may be related to these large amounts of long distance interactions, and want to turn them off explicitly. It's not a big deal running now, but as I scale these calculations up I want to make sure it's not doing more work than I need.
Does that make it more clear?
Thanks,
Michał
There is no need to have a large communicate cutoff. You lose time due to creation of useless ghost atoms.
If you have no nonbonded potential, then all those exclusion statistics have no meaning.
Remember that you made a lot of tests with a broken input.
Axel
Ok, thanks, I think I misinterpreted the exact purpose of the command
before. Everything works as expected now.
Thanks again for all the help,
Michal