(Apologies for submitting this again - I attached the wrong .colvars file the last time…)
Dear LAMMPS forum,
I have:
- two slabs of gamma-alumina (43.68 * 26.27 * 21.06Å; LWH) immersed in water
- initial separation in z of 25Å, with 30Å thickness of water between all other sides and the periodic box boundaries
- Total simulation box size: 45 * 86.27 * 127.12Å (x, y, z)
I am trying to pull the upper slab down towards the lower one, to a separation of 3Å, in increments of 1Å - the objective is eventually to extract the PMF using WHAM.
For the moment I’ve run quite a quick simulation, just to get an idea of what happens:
- 10k steps of 0.1fs, followed by 150k fs (1fs timestep) of just letting the system run.
- This is followed by (what should be) moving 22Å in increments of 1Å, spending 100k fs at each stage.
According to my .o file this seems to be happening (not able to attach in this message - 900kb attachment size limit). However, when I visualise the ouput.data file in VMD** it appears that this distance is closed in 9 frames (so ~90-100fs) which seems ridiculously fast.
**Attached screenshot composite from VMD: ‘transition.png’ (I’ve omitted the slab atoms for clarity)
I would be hugely grateful if anyone has any ideas on what’s possibly happening and why there seems to be this discrepancy. I’ve attached my run and .colvars files, if this helps.
I suspect it might be that I’ve made a mess with/misunderstood the idea of ‘centers’.
I’ve thought that centers, in the z dimension, would be midway through the thickness of the slab. Is it an absolute position, or relative?
Any explanation/suggestions/opinions on this very much appreciated.
Olivera
- (Please ignore the mention of SMD in my files; the name is inherited from my previous attempts to do this using SMD).
PMF-Olivera.colvars (756 Bytes)
run.in (1.83 KB)
Hello Olivera, thanks for uploading a file with actual parameter values (impossible to say anything from the previous file).
Looking at the run.in provided, it seems that you are outputting every 10000 steps, i.e. 1 ps. So shouldn’t 9 frames be 9 ps? Still fast, but I want to make sure the input file you attached is consistent with your results (seeing as the Colvars input in the previous email wasn’t).
I would recommend, before moving the restraint, to equilibrate first the system with a fixed restraint. This should be the procedure also if this was a production run. If the fixed restraint doesn’t work as you expect, that should tell you more clearly than when you are moving it.
Also, what is the value of the “dist” variable reported in the *.colvars.traj file?
Giacomo
Hi Giacomo,
thank you for replying. Yes, I meant 90k fs/9 ps for 9 frames - sorry!
I’ll try your suggestion with the fixed restraint and see how it goes.
I’m not quite sure about your last question - the colvars.traj file has a whole column of ‘dist’ values - which one are you referring to? The distance starts from -4.59974101180476e+01 and on the last stage is at -2.26273142414097e+01…I hope I understood your question correctly.
Kind regards,
Olivera
With an initial value of -46 Å for the variable “dist”, and a fairly stiff restraint centered at +46 Å, it is safe to say that “dist” will get quite a kick! The initial restraint energy would be 0.55.0(-46-46)^2 ~= 21,000 kcal/mol.
How about you switch the two groups to correct the sign of “dist”, or alternatively make the harmonic restraint consistent with it?
Giacomo
Hi Giacomo,
thank you for your help 
So I tried switching the two groups around: I think that was the problem! The upper slab now advances towards the lower one at the rate I’m expecting.
However: is there some way that I can constrain the movement of the upper slab in z only? At the moment it starts directly above the lower one but then starts hovering around moving in y (see attachment - last frame of simulation)…I’d like to keep it directly above the lower slab as it moves down, not tilting.
Maybe fix lineforce or fix setforce? I’ve tried both but they haven’t quite had the effect I was expecting…
Any help much appreciated.
Kind regards,
Olivera
You would need to also apply forces in the other two directions (X and Y I presume), but the distanceZ function is silent with respect to movements along those axes.
You could use, for example, the distanceXY function with a fixed restraint for that purpose. If you define one based on the two groups you have, you’ll restrain the relative movement. For the absolute movement, define two functions, each using one of the two groups and a dummyAtom group for the other group to give an absolute position.
Giacomo