During equilibration of pure ionic liquid (TDTHP-DCA) system to estimate molar volume, self-diffusion constant etc, there are spikes arising in the energy plot. The spikes continued to exist even after repeating the simulation by re-optimizing the structure and reducing the time step from 1 fs to 0.5 fs, moreover the system energy shifted to positive. The respective energy plots are attached with the mail.
This time I save veldcd file to see if there are any abnormalities in velocities at the instants spikes are arising. I compared the velocities between the normal instants and instants where spikes are observed in terms of maximum value, mean and standard deviation. And these values are quite relatable, proving there are no abnormal velocities at the spiked time instants. What could be the reason for this behavior and how to solve?
Potential Energy_1 fs: https://www.dropbox.com/s/uk257az28xj8632/Potential%20Energy_1fs.jpg?dl=0
Potential Energy_0.5 fs: https://www.dropbox.com/s/ugmw83ws2gilhyg/Potential%20Energy_0.5fs.jpg?dl=0
Did you check the pressure? Is this caused by the fluctuation of pressure?
What is the source of your force field parameters? Have they been used in previous studies? What can you observe from the configurations in the frames where the energy spikes occur? Perhaps you just have some molecules that have moved around such that there is an unfavorable charge-charge interaction that causes the energy to rise. It also seems as if the problem is going away over time, with the spikes lowering in magnitude, but your simulations are very short; have you extended them at all? I would also note that whatever "reoptimization" you did before switching to a 0.5-fs time step dramatically increased the starting potential energy, so that system is further from likely reality than the original system. If you want to test the effect of time step, just do that - don't change anything else, because now you have two variables and you don't know which is doing what.
I doubt this is pressure-related; the instantaneous pressure in an MD simulation fluctuates wildly at every time step, not infrequently like these values.
I'm using OPLS-All Atom parameters reported in Lopes and Padua, J. Phys. Chem. B, vol. 110, no. 39, pp. 19586–19592, 2006. These parameters have been used in the literature, however not particularly on this system. To test the effect of time step, the simulations were performed at 1 fs, 0.5 fs and 0.2 fs (currently running) time steps starting with the same initial configuration. All systems are converging to similar energies and though the frequency of spikes have relatively reduced for 0.5 fs system, they continue to exist even with 0.2 fs. So large time steps are not the reason for such behavior. I checked the frames where spikes are arising once again, and could not find any anomalies under visual observation.
This is an equilibration stage, so I'm terminating the simulation as the system energy converges to a plateau, which is happening within 2 ns. Charge-charge interactions are a possibility as the ionic liquids have very active electrostatic interactions, but I would like to understand how normal is this behavior?
True that the instantaneous pressure varies wildly in any MD simulation, and probably is not the cause either.
Examine the individual components of the potential. One (or more) should stick out. If it's coming from unfavorable charged contacts, the electrostatic energy will show the same behavior.
As Justin said, extract one of the conformations that has a high potential energy and take a look at the components. That should point you in the right direction. My guess is the non-bonded LJ term. If you print out the gradients for that conformation, you can see which atoms have the highest forces acting on them. In a liquid with a thermostat resetting velocities (what thermostat are you using? what time constant? try rerunning one of your simulations with a longer time constant) , it's quite possible for two atoms to get too close: you probably wouldn't see it by just inspecting the system visually, but it might be enough to jack up the LJ term.
My guess is that you have problems with LJ. Have you checked the non-bonded interaction calculation setup?
Thanks everyone for the helpful suggestions. Breaking down the potential energy to individual components for the simulation run with 1 fs time step, surprisingly it is not the electrostatic energies but apparently the bonding and kinetic energies are responsible for this behavior. Particular to bond energies, I identified two bonding interactions in DCA anion showing spike behavior. The respective force field parameters were taken from the literature cited in the previous comment, and I have cross checked the numbers multiple times. After this development, I do not think these parameters can be considered trustworthy anymore.
What are the bond constants on those two bonds? Not in the literature but in the actual topology file you are feeding to your simulation program? But the main question is are those bond spikes he cause or the after-effect of your total energy spikes? In the total energy plots you showed, there were spikes of over 6000 kcal/mol but here you have spikes of at most 1600kcal/mol. Could it be that the bond spikes you see are reactions to something else (like bad vdW contacts) happening immedietly before?
N3-CZ: 502.62 kcal/mol; CZ-NZ: 925.6692 kcal/mol, the latter is a triple bond.
True that the spikes have a peak range of ~6000 kcal/mol. You can see from the figure that the kinetic energies are the reason for that. The nonbonded energies including vdW are absolutely fine.
- Badges
- Science topic