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 < > 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. > -- 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.