MD with CP2K combined with PLUMED
161brun... at gmail.com
Thu Jul 6 13:48:08 UTC 2017
I have a question regarding CP2K MD simulations using PLUMED. This might
become lengthy to explain, but I hope that someone can perhaps give some
I will make corresponding files available through google drive.
The aim is to simulate the rearrangement of a substrate coordinated to a
metal centre. Specifically, this involves a "flip" of the coordination site
of a double bond.
>From previous static calculations I have found out a possible pathway for
this, which serves as a guidance for setting up the MD simulations. I have
defined 3 structures (start, transitions state, final) in the file
*path_flip.pdb**. *This should give you an idea what process I am
interested in, essentially the substrate flips around. As a note aside,
here I am interested in a molecular system, i.e. these simulations will use
a cubic supercell, and perhaps this will introduce some problems (see
below). Later on, I wish to return to the actual system, and embed this
molecular portion back into it's solid state environment (for various
reasons it is good to look at the molecular system first).
Now, since this process can be viewed as a conformational change, an
approach that appears reasonable here is the sampling method proposed by
Parinello (J. Chem. Phys., 2007, 126, 054103), which uses as the collective
variable an expression that measures the root-mean-square-deviation (RMSD)
of the sampled structures along a defined reaction path (as defined in
*path_flip.pdb*). This method appears to work quite nicely with PLUMED.
Indeed, a quick MD simulation (500 steps) using the movingconstraints
keyword in PLUMED and forcing the trajectory along the defined path shows
that the transition is smooth and follows the defined path
The reference point is always provided in the *path_flip.pdb* file. Note
also that only the atoms of the substrate are taken into account, in order
to remove "noise" from the movement of the other atoms (this is specified
in the .pdb file). Having established that this is a suitable collective
variable describing the rearrangement under study, I then moved on to
sample the trajectory using metadynamics (see
I have done 5 MD runs so far (4 restarts), and this is where things get
confusing. I have uploaded the files corresponding to the last restart. The
file *movie.xyz *file contains the complete trajectory, which I have
concatenated from the individual simulations. In each case I would follow
the following protocol:
(i) copy the CP2K .restart file from the previous md run to new folder
(ii) copy path_flip.pdb, plumed.dat, HILLS and COLVAR files to this folder
(iii) resubmit job.
*plumed.dat* contains the RESTART keyword, and I hope this is sufficient
information to continue with the simulation. The PLUMED.OUT file reports
that 192 Gaussians were read, which implies that the restart has worked.
Previous runs restarted from fewer numbers of hills, naturally.
What are the problems:
(i) From the COLVAR file it appears that the collective variable as defined
above (=path) is changing and the substrate is undergoing some
conformational changes. In the start the value is around 1.0, fluctuates
around this value for a while, then changes rapidly at around t=13-14 ps,
appears to reach the transition phase at 20 ps and is a t the final
structure thereafter (although there is some strong fluctuations going on
there). However, none of this is reflected in the actual trajectory, where
the substrate moves around but never undergoes the transition. There
appears to be a dichotomy between the values of the CV and the actual
(ii) It also appears that there are jumps in the sampling, as seen from the
COLVAR and HILL file. For example after 8 ps the next hills are only
spawned at 9.44 ps. How does this happen? Are my settings not right?
(iii) Another point of concern relates to the fact that this simulation is
molecular, and uses a cubic box (25 Ang). You can see that the molecule
rotates a lot, and I wonder if this is affects the estimation of the RMSD
value (i.e. the CV). I believe that this method can deal with translations,
but not so sure about rotations. HAs anyone experience with this, when the
molecule moves a lot during the simulation. I am also using the center
coordinates keyword (WAVELET!), how does this affect the results? Is it a
problem that the reference structures are in fact sitting on one corner of
the box, while the molecule for simulations is at the center (I think not,
since this should be accounted for through transnational symmetry).
Anyways, I would appreciate any reassurance here.
It appears that this is a neat way of doing these simulations, but it seems
to be also quite tricky to get it right. I am fairly new to this kind of
simulations, so any help would be very much appreciated. The initial
simulations hold promise, now I am trying to get the actual mtd right to
get some reliable results. Ultimately I would like to sample the system in
the solid state, i.e. using periodicity, in order to model this
rearrangement in dependence of the crystal structure. It would be nice to
work out the reason for variations in the barriers for this conformational
change in different systems, which we also study experimentally.
I am also using umbrella sampling as an alternative to the mtd part. I am
also aware that other collective variables could be defined to sample this
process (even without resorting to PLUMED), but I am simply interested in
making this work.
Thanks for your help, and apologies for this lengthy post. I hope that
there are some experts out there that can help.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the CP2K-user