Dear Lammps users,
I have a simple question regarding bonded molecules and in particular a case involving a dumbbell molecule.
As a test I have set up a configuration of two LJ dumbell molecules which are read in via a read_data command. The parameters are therefore fully detailed within the input script and data file which I have attached therefore I shall not detail them here.
I calculate the pair interaction energy and harmonic energy from an external code as a test of principle that I understand the manor in which lammps reads and handles data. The problem however is that I get a different pair energy from that of lammps (The harmonic components match in both cases).
I have considered various components of the interactions which lammps might be omitting (obviously due to my ill understanding of it’s operation), however I have been unable to see my mistake.
Any help with this would be much appreciated.
Many thanks, Jeff.
in.dumb (470 Bytes)
data.dumb (556 Bytes)
Dear Lammps users,
I have a simple question regarding bonded molecules and in particular a case
involving a dumbbell molecule.
As a test I have set up a configuration of two LJ dumbell molecules which
are read in via a read_data command. The parameters are therefore fully
detailed within the input script and data file which I have attached
therefore I shall not detail them here.
I calculate the pair interaction energy and harmonic energy from an external
code as a test of principle that I understand the manor in which lammps
reads and handles data. The problem however is that I get a different pair
energy from that of lammps (The harmonic components match in both cases).
I have considered various components of the interactions which lammps might
be omitting (obviously due to my ill understanding of it's operation),
however I have been unable to see my mistake.
well, LAMMPS does what it says it does. but we cannot know
what your external code is doing. there are certain conventions
that are implemented differently in different codes. e.g. sometimes
the 4 is included in the epsilon, sometimes it is (like in LAMMPS).
axel.
Did you look at the special_bonds command? The pairwise
energy is affected for bonded atoms by those settings.
If I recall, the default is turn off pairwise LJ between bonded
atoms.
Steve
Thank you for the replies. Yes I am aware of the different conventions and was pretty careful when reading the lammps documentation to ensure I followed them when calculating the energy by hand.
I am also aware of the special_bonds command and read the defaults for my particular pair style and as you rightly point out this default is to not calculate pair terms for bonded atoms, as is the case in my manual calculation.
I also include the shifts in the same manner that lammps does i.e. subtracting the value of the pair potential for r=rcut.
One thing that isn't entirely clear to me from reading the documentation, and perhaps this is covered in a section which talks about the manipulation of extensive variables which I have not read yet, is whether the default thermo output is that of per-atom values or not.
Apologies for the vagueness of my initial "question". I guess I was hoping that I was not including something basic that would be very apparent to an experienced user.
Many thanks, Jeff.
Did you look at the special_bonds command? The pairwise
energy is affected for bonded atoms by those settings.
If I recall, the default is turn off pairwise LJ between bonded
atoms.
It's true. He is using the special_bonds command (with 0 0 0).
I don't know if this could be the problem, but take a look at your cutoffs.
You are using "lj/cut" with a 3.2 sigma cutoff:
pair_style lj/cut 3.2
When you calculate the Lennard-Jones energies yourself, did you
remember to throw-away contributions between particles separated
beyond 3.2 sigma? (Last I checked, your system seemed to be pretty
dilute.)
Also take a look at the "shift" pair correction.
http://lammps.sandia.gov/doc/pair_modify.html
If I understand correctly, then without the "shift" correction, there
is a discontinuous jump in the energy at the cutoff distance
threshold. The "shift" correction will elliminate this discontinuity
by raising the pairwise energies between particles (within the cutoff
distance) by a factor of approximately 1/(rcut^6) (approximately
0.001 for your system). There are other pair_styles which use
different cutoff schemes and don't need these particular corrections.
Cheers!
Andrew
(P.S. Regarding the "tail" correction: currently you're not using it
either. However I think this is fine. I don't think that the "tail"
correction was designed for use with molecular systems. And the
system you posted last time was so dilute, you probably don't need to
worry about it.)
Oops you posted your reply while I was posting mine.
I also include the shifts in the same manner that lammps does i.e. subtracting the value of the pair potential for r=rcut.
I'm under the impression that if you wan't lammps to do this, you need:
pair_modify shift yes
to your input script
One thing that isn't entirely clear to me from reading the documentation, and perhaps this is covered in a section which talks about the manipulation of extensive variables which I have not read yet, is whether the default thermo output is that of per-atom values or not.
I think the energies reported are the total energies (not per-atom).
One thing that isn't entirely clear to me from reading the documentation, and perhaps this is covered in a section which talks about the manipulation of extensive variables which I have not read yet, is whether the default thermo output is that of per-atom values or not.
If you use the "norm" option in thermo_modify then it will be a per
atom quantity for an "extensive" variable. The default values for the
keyword "norm" depends on the units style. norm = yes for unit style
of lj, norm = no for unit style of real and metal. A more detailed
discussion is found here
http://lammps.sandia.gov/doc/thermo_modify.html.
Thank you for the replies, I actually worked out the problem right after I sent my last message. You were correct Andrew. I was subtracting the cut-off value correctly however I forgot to include an "if" statement to neglect contributions above my cut-off.
A classic case of reading my code so many times that I had lost the ability to objectively look for errors 
In any case the replies regarding the general use of cut-offs and extensive variable data formats have been of much use.
Many thanks, Jeff.