I am only familiar with the Morse potential in reduced LJ units, but I think an alpha of about 2/r0 is quite small, so the potential is very wide and overlap not that costly. This might contribute to your observations as well.
Aidan: R0 and Rmin are the same thing. This isn’t being squeezed in fact quite the opposite, which you can see from the fact that the RDF doesn’t level off to one, but in fact shows a depression in the density. I did run the squeezed simulation before sending you the current simulation. It looks completely different tall and spiky.
Stefan: I think in LJ units it would be. It corresponds to a force constant of 800kcal/Å-1. Why are you dividing it by R0?
Has anyone out there run this pair style successfully?
thanks
matthew
Because that was my natural length scale in reduced units, i.e., I usually put r0 to 1. Just plot the potential and see if it is reasonable. A harmonic approximation of the Morse potential around the minimum gives an effective spring constant of 2alpha^2D… For a Lennard-Jones potential it is 72*epsilon I believe. You can also use the pair_write command to see if the pair potential looks as expected or not.
For me in LJ units it works fine.
See attached. According to the function in the manual, the energy curve is plotted for 100,2,1.5. While reasonable in terms of the width, the minimum corresponds to the first minima not the first maxima. So again either I have done something really stupid or there is an error in the code.
matthew
morse_rdf.pdf (38 KB)
See attached. According to the function in the manual, the energy curve is
plotted for 100,2,1.5. While reasonable in terms of the width, the minimum
corresponds to the first minima not the first maxima. So again either I
have done something really stupid or there is an error in the code.
i think this is more a general problem. what is unusual about your
potential parameters the depth of the potential. 100 kcal/mol is at the
order of magnitude of a covalent bond. (C-H or C-C). dispersion energy is a
couple orders of magnitude smaller, i.e. less than 1 kcal/mol.
remember that the relation between log(g(r)) and a pairwise additive
potential only holds for dilute systems and correspondingly weak
interactions. that is not the case here. if you look at your plot, you can
rationalize what is going on: have a look at the _second_ peak in the g(r)
and recall that the number of atoms increases with O(r**3). thus you have
quite a few atoms comprising the second shell and those experience a strong
force toward the central atom, thus are also pushing the atoms in the
first shell toward it (until the repulsive part of the potential becomes
too repulsive).
axel.
Thanks Axel. That is basically what I said the first time. Matthew, your RDF clearly shows multiple coordination layers, a signature of a dense liquid. if you compute the RDF lower densities, you should see the first peak move outward asymptotically approaching R0 as density goes to zero. Conversely, if you keep the density the same, but jack up the stiffness parameter alpha, you will also see the peak move outwards in the same way.
Thank you Axel.
I find this phenomena most disturbing, but I think your description is correct. The critical values which appear to work are 10.0 3.0 1.5 while 10. 2.0 1.5 do not. I am now running an actual simulation of Liquid neon at 4K. Since the effect on behavior of D0 is strongly temperature dependent it suggests that there are temperatures under which the simulation will simply collapse in on itself for a given set of pair parameters. Most strange, but a very helpful answer.
matthew
Thank you Axel.
I find this phenomena most disturbing, but I think your description is
correct. The critical values which appear to work are 10.0 3.0 1.5 while
10. 2.0 1.5 do not. I am now running an actual simulation of Liquid neon at
4K. Since the effect on behavior of D0 is strongly temperature dependent
it suggests that there are temperatures under which the simulation will
simply collapse in on itself for a given set of pair parameters. Most
strange, but a very helpful answer.
i don't find it strange at all. it is just another manifestation of the
GIGO principle. the depth of your morse potential is that of a covalent
bond, but since it is also quite wide, atoms will "catch" other atoms that
fluctuate around whenever they come close enough, and then treat them as
bound atoms. so you'll be getting atoms with a large number of "bonded"
atoms, much more than you can accommodate in a single shell, and thus you
have a second shell resulting in the g(r) that you see. in a "normal"
setup, this would not happen, as covalent bonds are modeled from an
explicit bond topology, which will not allow additional bonds to happen.
also keep in mind, that non-bonded interactions between explicitly bonded
atoms are excluded (often beyond immediate pairs). thus using the same
potential and parameters for bonded interactions will work just fine,
similarly, as i have already pointed out, the interaction energy for
dispersion interactions is *much* smaller and thus the same effect cannot
happen with properly chosen parameters.
axel.