[gmx-users] TIP4P molecules stuck together

John Whittaker johnwhittake at zedat.fu-berlin.de
Thu May 9 17:56:56 CEST 2019


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
-------------------------------------------------------------------------------
Potential                    503436         --          0          0 
(kJ/mol)

and the output from gmx dump:

traj.trr frame 0:
   natoms=        10  step=         0  time=7.1440000e+03  lambda=         0
   box (3x3):
      box[    0]={ 3.36795e+01,  0.00000e+00,  0.00000e+00}
      box[    1]={ 0.00000e+00,  7.54019e+00,  0.00000e+00}
      box[    2]={ 0.00000e+00,  0.00000e+00,  7.54019e+00}
   x (10x3):
      x[    0]={ 2.10910e+01,  3.64700e+00,  2.75200e+00}
      x[    1]={ 2.11150e+01,  3.60000e+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.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    3]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    4]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    5]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    6]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    7]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    8]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      v[    9]={ 0.00000e+00,  0.00000e+00,  0.00000e+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.00000e+00,  0.00000e+00,  0.00000e+00}
      f[    4]={ 0.00000e+00,  0.00000e+00,  0.00000e+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.00000e+00,  0.00000e+00,  0.00000e+00}
      f[    9]={ 0.00000e+00,  0.00000e+00,  0.00000e+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 <
> johnwhittake at 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-request at 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-request at gromacs.org.
>




More information about the gromacs.org_gmx-users mailing list