Re: [gmx-users] TIP4P molecules stuck together
Dr. Shirts,
Thank you for the swift reply.
> If you have "unprotected" electrostatic sites (i.e. with nonzero repulsive
> terms directly on top of the charge), then there will always be some
> configurations with essentially infinitely negative energy.
That makes sense. Definitely something to think about, especially in these
simulations.
>Is your cap smoothly varying? If not, then your dynamics on hitting the
cap will be unphysical.
Indeed, our cap is implemented instantaneously and certainly introduces
non-physical dynamics when it is triggered. Our simulations consist of
non-interacting "tracer" particles that abruptly change resolution to
fully atomistic, interacting water molecules when they pass an interface
in the simulation box.
The force cap is a brute-force approach to ensure the simulation doesn't
explode when a particle crosses the boundary, gains atomistic features,
and finds itself in an unphysical configuration relative to an already
atomistically-resolved molecule (other groups have used a Monte Carlo move
to adjust the overlapping molecules... something we may consider in the
future).
We do not consider structural/dynamic properties within this interface
region where the force-cap is triggered by the immense energies due to
particle-particle overlap.
>How are the forces propagated into the energies (if grad U
>=/= F, then weird non-newtonian physics will also happen).
The forces are normalized to 2000 (iff > 2000) just before the velocities
are calculated in the first step of the stochastic dynamics integrator.
> What are the energies? Are they lower or higher than zero?
>From the single-point energy calculation of the dimer configuration, the
potential energy output is:
Energy Average Err.Est. RMSD Tot-Drift
---
Potential503436 -- 0 0
(kJ/mol)
and the output from gmx dump:
traj.trr frame 0:
natoms=10 step= 0 time=7.144e+03 lambda= 0
box (3x3):
box[0]={ 3.36795e+01, 0.0e+00, 0.0e+00}
box[1]={ 0.0e+00, 7.54019e+00, 0.0e+00}
box[2]={ 0.0e+00, 0.0e+00, 7.54019e+00}
x (10x3):
x[0]={ 2.10910e+01, 3.64700e+00, 2.75200e+00}
x[1]={ 2.11150e+01, 3.6e+00, 2.83200e+00}
x[2]={ 2.11700e+01, 3.69400e+00, 2.72700e+00}
x[3]={ 2.11040e+01, 3.64700e+00, 2.75900e+00}
x[4]={ 2.10960e+01, 3.64700e+00, 2.75500e+00}
x[5]={ 2.11700e+01, 3.72600e+00, 2.72800e+00}
x[6]={ 2.11700e+01, 3.65300e+00, 2.66500e+00}
x[7]={ 2.10770e+01, 3.73900e+00, 2.74900e+00}
x[8]={ 2.11580e+01, 3.71800e+00, 2.72200e+00}
x[9]={ 2.11650e+01, 3.72300e+00, 2.72500e+00}
v (10x3):
v[0]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[1]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[2]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[3]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[4]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[5]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[6]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[7]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[8]={ 0.0e+00, 0.0e+00, 0.0e+00}
v[9]={ 0.0e+00, 0.0e+00, 0.0e+00}
f (10x3):
f[0]={-3.67300e+07, -3.67266e+07, 1.11572e+07}
f[1]={-3.54829e+02, -4.32763e+01, 3.86598e+00}
f[2]={-4.24006e+04, 9.05089e+04, -1.34170e+04}
f[3]={ 0.0e+00, 0.0e+00, 0.0e+00}
f[4]={ 0.0e+00, 0.0e+00, 0.0e+00}
f[5]={ 3.67632e+07, 3.66652e+07, -1.11462e+07}
f[6]={ 3.81472e+03, -1.38062e+04, -2.40026e+02}
f[7]={ 5.73042e+03, -1.52270e+04, 2.66433e+03}
f[8]={ 0.0e+00, 0.0e+00, 0.0e+00}
f[9]={ 0.0e+00, 0.0e+00, 0.0e+00}
note: there are 5 atoms per TIP4P molecule in our case because we use a
virtual site constructed from the other atoms as the non-interacting
"tracer" particle (all 5 atoms exist at once as a hybrid molecule all the
time, interactions are just switched on for OW, HW1, HW2, and MW when they
transition to the atomistically-resolved region)
Sorry for the wall of text, I hope what I said makes sense. I appreciate
the help.
John
> On Thu, May 9, 2019 at 8:43 AM John Whittaker <
> johnwhitt...@zedat.fu-berlin.de> wrote:
>
>> Hi all,
>>
>> I have a rather strange question that I hope someone can shed some light
>> on.
>>
>> Before I begin, I want to note that I am pioneering some new
>> developments
>> of the Adaptive Resolution Simulation technique
>> (https://doi.org/10.1002/adts.201900014), so the simulations/techniques
>> I
>> am performing/implementing are fairly non-standard with respect to
>> normal
>> atomistic
Re: [gmx-users] TIP4P molecules stuck together
If you have "unprotected" electrostatic sites (i.e. with nonzero repulsive terms directly on top of the charge), then there will always be some configurations with essentially infinitely negative energy. For vanilla MD, these are essentially impossible to reach kinetically at any reasonable temperature or reasonable timestep, because the R^-12 of the neighboring atoms creates such an enormously large barrier. But with certain accelerated sampling approaches, you can skip over the barrier, and access these sites, which will either blow up the simulation, or get stuck forever. If you cap the forces, then weirder things will happen. Is your cap smoothly varying? If not, then your dynamics on hitting the cap will be unphysical. How are the forces propagated into the energies (if grad U =/= F, then weird non-newtonian physics will also happen). > the forces on each atom are massive (on the order of 10^7). What are the energies? Are they lower or higher than zero? On Thu, May 9, 2019 at 8:43 AM John Whittaker < johnwhitt...@zedat.fu-berlin.de> wrote: > Hi all, > > I have a rather strange question that I hope someone can shed some light > on. > > Before I begin, I want to note that I am pioneering some new developments > of the Adaptive Resolution Simulation technique > (https://doi.org/10.1002/adts.201900014), so the simulations/techniques I > am performing/implementing are fairly non-standard with respect to normal > atomistic simulations. > > With that in mind, I am simulating a box of TIP4P water and calculating > structural/static properties. My simulations utilize a force-cap of 2000 > kJ/(mol nm) at each time step - i.e., when the force on an atom is larger > than +/- 2000, the force is automatically normalized to +/- 2000 to > prevent explosive forces due to atomic overlaps. > > For the most part, this works for the purposes of my simulations but I > have observed some water molecules "sticking" together in the > configuration shown here: > > https://www.dropbox.com/s/p5rkximspp25flf/tip4pDimer.jpg?dl=0 > > with a corresponding O-H radial distribution function (unnormalized) shown > here: > > https://www.dropbox.com/s/ez56db4qggv1iii/rdf_OH_long.jpg?dl=0 > > where there is a clear (albeit, small) probability of finding a hydrogen > atom an extremely short distance from an oxygen. > > The molecules travel together like this for several ps and then, for > seemingly no reason, split apart and carry on perfectly fine for the rest > of the simulation. > > I have performed a single-point energy calculation on this configuration > in vacuum and have found, as one would expect, the forces on each atom are > massive (on the order of 10^7). Yet, the molecules do not repel and seem > to prefer this configuration for a short time. > > I have a feeling that this configuration is allowed when the forces are > normalized to 2000 and the molecules become trapped there. > > I am wondering if anyone may have some experience with TIP4P water > molecules taking on unphysical configurations for non-negligible times. I > have not tried this same simulation using TIP3P yet, so I'm unsure if it > has something to do with electrostatic interactions with the virtual site, > but I will test this tomorrow. > > Thank you for any information/speculation/guesses as to why this is > happening. > > - John > > -- > Gromacs Users mailing list > > * Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before > posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
