[gmx-users] TIP4P molecules stuck together
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.
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.
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 simulations.