Hi:
I have a very simple simulation that is supposed to be of liquid density based upon the paper I am modelling. That would be rho of about 0.8. I get twice that consistently no matter how I permute things in my script.
I am hopeful that there might be some simple thing I don’t understand about my script.
Here are the high points:
– I have 10 chains of 10 beads each. timestep 0.005
– I have a harmonic bond, and a standard lj12/6 with a cutoff of 2.3. This cutoff causes the chains to collapse in upon themselves and crank up the density. If I use a cutoff of 1.12 I get no collapse. However the paper requires a 2.3 cutoff.
– I am using and npt ensemble:
fix f1 all npt temp 1.0 1.0 100.0 iso 0.0 0.0 100.0
In less than 1 million steps the box shrinks to the point of crystallization
If I go with a nvt ensemble things are reasonable since the volume is fixed, however the paper I model uses npt. Could the damping parameters be something here? I experimented a bit with them to either no avail or distinctly broken behavior.
– That is pretty much it. Very simple. Is there something obvious here to the more experienced eye?
Thanks for all the help!
Not a polymer guy and far from a liquid expert but here is a comment:
You seem to have a small system running at zero pressure thus one would expect large pressure fluctuations in the npt ensemble. Have you tried increasing the system size to reduce the potentially large stresses on the simulation box?
I assume no other problems in your script and things like particles interacting with themselves across boundaries due to large/small cutfoff/box size have been ruled out.
Carlos
I assume your bond length is sigma? Also you can
look at the special_bonds settings. By default there
is no 1/2 1/3 or 1/4 interactions. You may want
to turn off only 1/2 interactions.
Steve
my bond length is harmonic with a K=555.5, and r0=0.967. Sigma is 1 in the lj potential.
I will look into special bonds.
I have rerun with a ‘number of chains’ of 50 and 100. The box is considerably larger and away from my 2.3 cutoff. It can definitely be argued that my n=10 box is too small, but I think 50 and 100 are enough bigger from my cutoff. They both still clumped up and went to crystallization. Out of interest if I crank up the temperature to something between 2 and 3 I get more agreeable densities, but then I am still not modeling the paper.
Thanks for all the help!
I was going to mention the “low” temp before but figured you would have been careful enough to check for this. After covering all the obvious checks you can always try reaching the authors of the paper for some clarification.
On another front here is another Q: Does it matter what kind of initial config you start with for the chains? I mean, on your way to reaching equilibrium do each chain first collapses onto itself to give you a self-contain cluster entity with the collection of these “units” forming the liquid? or do you see different chains clustering to render the so called crystal state? Please take my Q above with a grain of salt as I did not read all the fine print in your original post.
Carlos
i don't see any comment on changing special_bonds from the default
setting of lj/coul 0.0 0.0 0.0 to special_bonds lj/coul 0.0 1.0 1.0
with the default you can run into a situation where a chain of
o-o-o-o-o-o-o-o-o-o
can turn into:
8-8-o-o-8-8 or similar
and then the atoms at the ends have two lj spheres in about the same
location (they are excluded and thus don't repel each other and there
is not angle and dihedral potential to keep the chain straight) and
thus this acts as if you have doubled the epsilon and that can explain
the shrinkage of the system.
axel.
Thanks to you all! The special_bonds command fixed me up. I am very close to their densities now. If I read their paper correctly they are using special_bonds lj 1 1 1, so I will use that for now.