water: from cube to sphere

Hello,

I’ve initialized and melted (while gradually increasing temp from 0-298) a cube of ice, using TIP3P charmm. The result is in close agreement w/ the literature when it comes to density and rdf. (final dump attached)

Now this is my problem: Next, I’ve cropped the cube into a sphere (data file attached), and tried to run that structure thru a similar minimization/melting procedure before using it in my final simulation. I’ve tried nve, nvt, and the npt ensembles but the sphere just shrinks unexplicably.

In the latter equilibration, I chose to expand the box so that the sphere is not too affected by its periodic images (the situation is similar to what my final objective simulation looks like, where there is fluid-solid interaction and the sphere will be placed on top of a periodic structure with the walls far away from it)

I’ve also tried using a procedure previously outlined in one of the threads of this forum, in which I fix the outer shell of the sphere and run the equilibration. I would then unfix a portion of that shell and run again, until I have unfixed the entire sphere. It doesn’t work well either. (input script attached). I’ve also tried different minimizations to proceed to sim to no avail.

I doubt that the physics of the simulation are wrong, so I suspect that I’m not using some of the lammps options right.

I was wondering if anyone has successfully equilibrated a water droplet in lammps? I’d appreciate any pointers.

(P.S. I realize I need to average the production data, but for the purposes of demonstrating my problem, the “last_cube_data” is actually the snapshot of the last water cube after melting.)

Many thanks
bs

last_cube_data (1.68 MB)

water_charmm_R35 (327 KB)

in.sphere.charmm.2 (1.4 KB)

Hello,

I've initialized and melted (while gradually increasing temp from 0-298) a
cube of ice, using TIP3P charmm. The result is in close agreement w/ the
literature when it comes to density and rdf. (final dump attached)\

it would be helpful, if you would provide the input to that as well,
or at least explain where you made changes relative to the droplet setup.

Now this is my problem: Next, I've cropped the cube into a sphere (data file
attached), and tried to run that structure thru a similar
minimization/melting procedure before using it in my final simulation. I've
tried nve, nvt, and the npt ensembles but the sphere just shrinks
unexplicably.

why that complicated effort? just expanding the cell of the
liquid cube and some equilibration would have resulted in
a droplet as well. i don't see, how you can do npt properly
for obtaining a droplet.

In the latter equilibration, I chose to expand the box so that the sphere is
not too affected by its periodic images (the situation is similar to what my
final objective simulation looks like, where there is fluid-solid
interaction and the sphere will be placed on top of a periodic structure
with the walls far away from it)

I've also tried using a procedure previously outlined in one of the threads
of this forum, in which I fix the outer shell of the sphere and run the
equilibration. I would then unfix a portion of that shell and run again,
until I have unfixed the entire sphere. It doesn't work well either. (input
script attached). I've also tried different minimizations to proceed to sim
to no avail.

I doubt that the physics of the simulation are wrong, so I suspect that I'm
not using some of the lammps options right.

I was wondering if anyone has successfully equilibrated a water droplet in
lammps? I'd appreciate any pointers.

yes, i haven't done it in LAMMPS and with other codes.
not a big deal, but you have to decide how to handle
evaporating water molecules and what kind of final
setup you are after.

axel.

Hello,

I've initialized and melted (while gradually increasing temp from 0-298) a
cube of ice, using TIP3P charmm. The result is in close agreement w/ the
literature when it comes to density and rdf. (final dump attached)

Now this is my problem: Next, I've cropped the cube into a sphere (data file
attached), and tried to run that structure thru a similar
minimization/melting procedure before using it in my final simulation. I've
tried nve, nvt, and the npt ensembles but the sphere just shrinks
unexplicably.

my guess is that you forgot to wrap back the molecules back
into the principal unit cell before you cut out the sphere. since
you seem to be using VMD, you'd have to execute the command:

pbc wrap -compound fragment -center origin

after loading the last equilibrated
did you wrap back all your molecules into the principal unit cell
before your cut out the sphere.

In the latter equilibration, I chose to expand the box so that the sphere is
not too affected by its periodic images (the situation is similar to what my
final objective simulation looks like, where there is fluid-solid
interaction and the sphere will be placed on top of a periodic structure
with the walls far away from it)

it is still worth pondering whether the (otherwise unscreened)
periodic interactions between the water droplets is desirable
compared to the error you are making, by not using PPPM
at all and a somewhat longer cutoff.

apropos cutoff, your choice of 5/6 anstrom for the lennard-jones
and realspace coulomb seems unreasonably short. typical is
about twice as much.

I've also tried using a procedure previously outlined in one of the threads
of this forum, in which I fix the outer shell of the sphere and run the
equilibration. I would then unfix a portion of that shell and run again,
until I have unfixed the entire sphere. It doesn't work well either. (input
script attached). I've also tried different minimizations to proceed to sim
to no avail.

this far to complex a scheme. what i would recommend in this case
is the following:

- take your equilibrated cube and wrap the positions into the
  principal unit cell as mentioned above

- set up a spherical bounding potential around that cube,
  so that water molecules cannot escape. e.g. with:

region container sphere 0.0 0.0 0.0 40.0 side in units box
fix bound all wall/region container harmonic 100.0 0.0 5.0

run as long as needed to have a reasonable droplet.
if you want to speed up dissipation of energy, you
could use fix nve plus fix langevin instead of fix nvt.
or use the former at the beginning and then switch
to the latter (or just plain nve).

I doubt that the physics of the simulation are wrong, so I suspect that I'm
not using some of the lammps options right.

nope. see above.

cheers,
     axel.

Thanks for the feedback, Dr. K.

I had tried all of what you mentioned, including the unwrapping, a 9/10 cutoff, etc. Some of these attempts were more of a random attempt to try and understand what’s happening than an educated guess.

I did not, however, try the fixed LJ wall, and it did somewhat accelerate the progress. After running for 100+ ps, the rdf plot is weird – I just don’t understand why I get the same profile as the previous cube but with strangely high probability values on the ordinate. (This misunderstanding may be due to my modest knowledge of rdfs, but I suppose it has to do with the periodicity and size of the box?)

The attached plot is obtained from the lammps rdf compute.

Am I good to go with what I have now?

sphere_rdf.pdf (9.96 KB)

upload (663 KB)

Thanks for the feedback, Dr. K.

I had tried all of what you mentioned, including the unwrapping, a 9/10
cutoff, etc. Some of these attempts were more of a random attempt to try and
understand what's happening than an educated guess.

you can't do proper MD simulations without knowing what you are doing.
if you are still just guessing, then you should go back and do it again,
step by step and come up with an explanation/justification for each step.
even though it may be a bit painful to move so slow right now, anything
you don't do well at this step, can hurt you much more later on.

I did not, however, try the fixed LJ wall, and it did somewhat accelerate
the progress. After running for 100+ ps, the rdf plot is weird -- I just

no, it isn't.

don't understand why I get the same profile as the previous cube but with
strangely high probability values on the ordinate. (This misunderstanding
may be due to my modest knowledge of rdfs, but I suppose it has to do with
the periodicity and size of the box?)

yes, you need to think more about how g(r) is computed *and* normalized.

The attached plot is obtained from the lammps rdf compute.

Am I good to go with what I have now?

i can't say and i won't. *you* have to be convinced that your results are
correct, since it will be *you* that has to defend the results against
future criticism. it won't help you a bit if you argue "but axel said
it was ok".

life is one of the hardest. sorry,
    axel.