MD with CP2K combined with PLUMED

[Tobias Kraemer]

Dear all, 


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 
advice.
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 
(see https://drive.google.com/open?id=0B3k7Eu2qqsBBNUl2dFlneWROajA). 

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 
https://drive.google.com/open?id=0B3k7Eu2qqsBBQ21hczVpbUt1azQ).

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 
as well
(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 
structure. 

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


Tobias
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20170706/91dcf116/attachment.htm>