[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