Mod:Hunt Research Group/BmimCl equilibration
Equilibration of 216 pairs of [Bmim][Cl] (Notes by Ling)
1. NVT (20ps, 2000K)
2. NVT (29 ps, 20 ps equilibration, average over 9 ps, 58.5 hrs on 4 CPUs) => to get good T and P.
temperature 313
ensemble nvt hoover 0.1
timestep 0.00001 ps = 0.01 fs
cutoff 18.5000
delr width 1.5000
rvdw cutoff 18.0000
ewald precision 1.0E-6
shake tolerance 1.0E-8
Notes: SHAKE 10^-8, timestep = 0.5 fs (not working), 0.2 fs (not working), 0.1 fs (working). So am using 0.1fs for timestep in the next couple of restart runs. Hope I can increase the timestep to 0.5fs after that. People are using 2fs or 0.2 fs or 0.5 fs in the production run in the literature for [BMIM][Anion], and Leo is using 2fs for his water system in the production run. From the literature and Leo, they have used a much smaller timestep (can be as small as 1/1000 of the timestep in the production run) at the beginning of the simulation when the system is far from equilibrium.
The easiest way is to set 'ewald precision 1d-6' and cutoff ~ 1/2 cell size. And I have checked that 'eng_cou = -vir_cou' in the OUTPUT file so the Ewald sum is correctly converged in my new simulation. (Previously using 'ewald precision 1d-5' and cutoff < 1/2 cell size did NOT give 'eng_cou = -vir_cou'.)
T: In the first 20ps, it is simply rescaling the velocities, whenever the temperature goes outside a given tolerance from the target. After the equilibration time, it applies the Hoover thermostat. So it is doing proper NVT. Before 20 ps, it is still an NVT simulation, but the velocity rescaling may introduce unphysical effects, so those steps cannot be used reliably in average. And it only calculate averages in the last 9 ps. Final average temperature, after equilibration, was 313.66 - the equilibration was successful. Oscillations +/- 200 K are quite high. This suggests that the system may be too far from the equilibrium. Hope the oscillations will become smaller in the next restart run.
P: Can clearly see the end of the equilibration in the pressure plot too, and that the oscillations are sensible after that point too.
RDF => looks liquid. Simulation time is short, so see random spikes (possibly due to rdf 100 => change to rdf 5 next time, but also the simulation times must be much longer to get something smooth) within the broader oscillations coming from the shell structure. Aim for RDF: get an overall smoother curve, with oscillations in g(r) converging to 1 at large distance.
3. NPT (300 ps, 200 ps equilibration) => to equilibrate volume.
Strategy: steps 2800000, equilibration 2800000, after several runs until 200 ps, then start using equilibration 0. But take this into account when doing averages in the end (if it is necessary). Using 120 hrs per run. Once happy with the volume, can try NVE and take some averages (e.g. diffusion coefficients) and see if it compares to experiment/published data.
3.1 NPT (28.6 ps, 28.6 ps equilibration, 120 hrs on 4 CPUs). temperature 313 ensemble npt hoover 0.1 0.1 timestep 0.00001
Drift ~ 10 degrees over a period of 29 ps. => equilibrate well
Oscillations ~ ±20 katm.
At ~29ps, vol = 54770 A^3, => density = 6.549 mol/dm^3 (In Lop04, exp: density = 6.183 mol/dm^3; theo: 6.000±0.018 mol/dm^3).
Volume is steadily decreasing. This may be due to the initial volume was too large for this system. In any case, the volume is clearly not converged yet, the system is trying to shrink even more. But it is not exacly linear with time, so it is converging (slowly) to its correct value. The density will obviously become larger as the volume converges. Maybe the initial estimate of the density was not completely correct for this specific system? Density higher than expected. It looks as if the system is still undecided between solid and liquid phase. But it could also be that the volume itself is oscillating over very long periods of time, so would be the density.
If there are no thermostats or barostats, it should oscillate around an average. In principle, when computing physical avrages, one should run NVE (i.e. not thermostats or barostats and fixed volume).
3.2 NPT 4 CPUs, 27.6ps, 120 hrs.
T: Drift ~ 2 degrees over the last period of 10 ps. Temp ~ 310 K. => equilibrate well
P: Oscillations ~ ±10 katm over the last period of 10 ps. Press ~ -2.03 katm. => equilibrate well.
At ~57ps, vol = 52620 A^3, => density = 6.816 mol/dm^3.
Volume is steadily decreasing.
3.3 NPT 313K 4 CPUs, 100 ps, running.
To do: perhaps need to boost the temperature again to 1000K?
Trajectory looks frozen. Check with Ruth whether it is liquid.
Below are results for a test restart run for 1ps.
To do: check the position and height of the first peak and compare to Exp data. The spikes might be due to very short simulation (a test run over 1ps). The Cl are disordered (this gives many small peaks), but they had no time to move. When they move, the space between small spikes decreases, the height of the small spikes decreases, and the overall shell structure should become more apparent.
The volume is stabilising to a larger value than before. The temperature had that spike, but then is also stabilising. The structure looks OK. It might be that the system has become a liquind only now, and it has moved to a sensible density. Perhaps before it was just disordered, without being a liquid yet.