PLUMED microtutorial ● Log into the machines ● Copy the plumed tutorial directory cp -r /group/dft-nao2012/tutorial_PLUMED ● You have four folders 2_how_to_enable_plumed 3_simple_bias 4_moving_bias 5_metadynamics sum_hills
The system ● Cl- CH3-Cl : undergo a SN2 reaction ● Studied by many (here ref is S. Yang, P. Fleurat-Lessard, I. Hristov, T. Ziegler, J Phys Chem A 108 (43) (2004) 9461–9468. )
ENABLING PLUMED Enter the dir 2_how_to_enable_plumed ● In the file control.in you find a plumed enabling line plumed .true. (requires a ● plumed.dat ) It is a molecular dynamics run (ensemble, length of the simulation, initialization ● specifications) MD_run 5.0 NVT_parrinello 300 0.1 MD_time_step 0.001 MD_clean_rotations .true. MD_restart .false. output_level MD_light MD_maxsteps -1 MD_MB_init 300 MD_RNG_seed 12345
plumed.dat DISTANCE LIST 1 <g1> g1-> 3 4 5 6 g1<- DISTANCE LIST 2 <g1> PRINT W_STRIDE 2 ENDMETA Distance in Bohr, energies in Hartree, angles in rad.
Running the calculation and outputs ● Source the environment (set the FHIAIMS and utilities in the path) source ../sourceme.sh ● Send the first calc on 4 processors mpiexec -np 4 $FHIAIMS ● Wait for some secs and you get PLUMED.OUT COLVAR
PLUMED.OUT ::::::::::::::::: READING PLUMED INPUT ::::::::::::::::: |- GROUP FOUND: g1 DISTANCE LIST 1 <g1> |- GROUP MEMBERS: 3 4 5 6 g1-> 3 4 5 6 1-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON g1<- |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 DISTANCE LIST 2 <g1> |- 1st SET MEMBERS: 1 PRINT W_STRIDE 2 |- 2nd SET MEMBERS: 3 4 5 6 ENDMETA 2-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 |- 1st SET MEMBERS: 2 |- 2nd SET MEMBERS: 3 4 5 6 |-PRINTING ON COLVAR FILE EVERY 2 STEPS |-INITIAL TIME OFFSET IS 0.000000 TIME UNITS |-ANALYSIS: YOU WILL ONLY MONITOR YOUR CVs DYNAMICS |- DIFFERENT COLLECTIVE VARIABLE WILL BE CALCULATED AT DIFFERENT TIMES |--CV 1 WILL BE EVALUATED ONLY WHEN NEEDED (OUTPUT OR EXCHANGE TRIAL) |--CV 2 WILL BE EVALUATED ONLY WHEN NEEDED (OUTPUT OR EXCHANGE TRIAL)
PLUMED.OUT ::::::::::::::::: READING PLUMED INPUT ::::::::::::::::: |- GROUP FOUND: g1 DISTANCE LIST 1 <g1> |- GROUP MEMBERS: 3 4 5 6 g1-> 3 4 5 6 1-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON g1<- |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 DISTANCE LIST 2 <g1> |- 1st SET MEMBERS: 1 PRINT W_STRIDE 2 |- 2nd SET MEMBERS: 3 4 5 6 ENDMETA 2-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 |- 1st SET MEMBERS: 2 |- 2nd SET MEMBERS: 3 4 5 6 |-PRINTING ON COLVAR FILE EVERY 2 STEPS |-INITIAL TIME OFFSET IS 0.000000 TIME UNITS |-ANALYSIS: YOU WILL ONLY MONITOR YOUR CVs DYNAMICS |- DIFFERENT COLLECTIVE VARIABLE WILL BE CALCULATED AT DIFFERENT TIMES |--CV 1 WILL BE EVALUATED ONLY WHEN NEEDED (OUTPUT OR EXCHANGE TRIAL) |--CV 2 WILL BE EVALUATED ONLY WHEN NEEDED (OUTPUT OR EXCHANGE TRIAL)
PLUMED.OUT ::::::::::::::::: READING PLUMED INPUT ::::::::::::::::: |- GROUP FOUND: g1 DISTANCE LIST 1 <g1> |- GROUP MEMBERS: 3 4 5 6 g1-> 3 4 5 6 1-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON g1<- |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 DISTANCE LIST 2 <g1> |- 1st SET MEMBERS: 1 PRINT W_STRIDE 2 |- 2nd SET MEMBERS: 3 4 5 6 ENDMETA 2-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 |- 1st SET MEMBERS: 2 |- 2nd SET MEMBERS: 3 4 5 6 |-PRINTING ON COLVAR FILE EVERY 2 STEPS |-INITIAL TIME OFFSET IS 0.000000 TIME UNITS |-ANALYSIS: YOU WILL ONLY MONITOR YOUR CVs DYNAMICS |- DIFFERENT COLLECTIVE VARIABLE WILL BE CALCULATED AT DIFFERENT TIMES |--CV 1 WILL BE EVALUATED ONLY WHEN NEEDED (OUTPUT OR EXCHANGE TRIAL) |--CV 2 WILL BE EVALUATED ONLY WHEN NEEDED (OUTPUT OR EXCHANGE TRIAL)
PLUMED.OUT ::::::::::::::::: READING PLUMED INPUT ::::::::::::::::: |- GROUP FOUND: g1 DISTANCE LIST 1 <g1> |- GROUP MEMBERS: 3 4 5 6 g1-> 3 4 5 6 1-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON g1<- |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 DISTANCE LIST 2 <g1> |- 1st SET MEMBERS: 1 PRINT W_STRIDE 2 |- 2nd SET MEMBERS: 3 4 5 6 ENDMETA 2-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 |- 1st SET MEMBERS: 2 |- 2nd SET MEMBERS: 3 4 5 6 |-PRINTING ON COLVAR FILE EVERY 2 STEPS |-INITIAL TIME OFFSET IS 0.000000 TIME UNITS |-ANALYSIS: YOU WILL ONLY MONITOR YOUR CVs DYNAMICS |- DIFFERENT COLLECTIVE VARIABLE WILL BE CALCULATED AT DIFFERENT TIMES |--CV 1 WILL BE EVALUATED ONLY WHEN NEEDED (OUTPUT OR EXCHANGE TRIAL) |--CV 2 WILL BE EVALUATED ONLY WHEN NEEDED (OUTPUT OR EXCHANGE TRIAL)
COLVAR #! FIELDS time cv1 cv2 Tells you which fields 0.0040 3.641223429 5.668572649 are there 0.0080 3.673303217 5.643192493 0.0120 3.709813875 5.611646432 0.0160 3.742140982 5.580829936 0.0200 3.769550265 5.551903013 0.0240 3.789024104 5.527580793 0.0280 3.802357539 5.504956398 0.0320 3.809773359 5.478647195 CV1 CV2 Time -COLVAR IS BACKED UP ONLY ONCE IN COLVAR.old. If you have another COLVAR.old it will be deleted -Distances are reported in Bohr, energies in Hartree, angles in rad.
PLOTTING COLVAR ● Use gnuplot to plot the timeline gnuplot gnuplot> plot “COLVAR” u 1:2 w lp gnuplot> plot “COLVAR” u 1:3 w lp (will be useful later)
A simple bias ● Put an umbrella potential to bias the simulation and “encourage” the reaction ● Enter the other directory: cd ../3_simple_bias ● plumed.dat: g1-> 3 4 5 6 g1<- DISTANCE LIST 1 <g1> DISTANCE LIST 2 <g1> # the distance between C-Cl' and C-Cl Use diff of distances Use diff of distances DISTANCE LIST 1 <g1> DIFFDIST 2 <g1> UMBRELLA CV 3 AT 0.0 KAPPA 0.01 Impose a 0-centered PRINT W_STRIDE 1 spring on CV3 ENDMETA
Run the calc The COLVAR file is: #! FIELDS time cv1 cv2 cv3 vwall XX XX RST3 WORK3 0.0020 3.604130111 5.694642933 -2.090512821 0.021851 RST 3 0.000 0.000 0.0040 3.664550915 5.645137333 -1.980586418 0.019613 RST 3 0.000 0.000 0.0060 3.702547964 5.609970760 -1.907422796 0.018191 RST 3 0.000 0.000 0.0080 3.755708366 5.558410133 -1.802701766 0.016248 RST 3 0.000 0.000 0.0100 3.815181937 5.498574383 -1.683392446 0.014169 RST 3 0.000 0.000 0.0120 3.873267229 5.437891473 -1.564624245 0.012240 RST 3 0.000 0.000 0.0140 3.934708487 5.371097388 -1.436388901 0.010316 RST 3 0.000 0.000 Work CV1 CV2 CV3 Bias from Time Center of (dist 1) (dist 2) (dist1- dist 2) Harmonic the bias spring
Plot the timeline and compare with the previous ● Call gnuplot and compare the distances now and before (without bias) gnuplot> p “COLVAR” u 1:2 w lp tit “d1+bias”, “” u 1:3 w lp tit “d2+bias”, “../how_to_enable_plumed/COLVAR” u 1:2 w lp tit “d1”,”” u 1:3 w lp tit “d2”
Now the two distances interchange and cross each others
Unbiasing the distribution ● Run the script to calculate the distribution and plot it: grep -v FIELDS COLVAR |awk '{print $1,$4}' | ./distribution.awk >distrib ● Plot it with gnuplot and retrieve the original unbiased free energy
Unbiasing the distribution ● Run the script to calculate the distribution and plot it: grep -v FIELDS COLVAR |awk '{print $1,$4}' | ./distribution.awk >distrib ● Plot it with gnuplot and retrieve the original unbiased free energy gnuplot> p "distrib" u 1:(-(0.597)*log($2)- 0.5*0.01*627.51*($1*$1)) w lp,"" u 1:(- (0.597)*log($2)) w lp ,"" u 1: (0.5*0.01*627.51*($1*$1)) w lp
Bias potential , biased fes and unbiased fes Anything wrong or suspicious???
Compare first and second half First and second half are different!
The timeline (plot COLVAR) 1) Initially the system is far from the eq. position of the harmonic spring: temperature is high-> it is like simulating a different ensemble 2) When the thermostat start having a role than a bistability appears: very hard to converge-> use harder springs and combine more segment of free energy
Produce a movie ● Make a vmd-loadable movie create_xyz_movie.pl out >movie.xyz ● Load in vmd: vmd movie.xyz ● Use CPK representation (or, better, dynamic bonds and CPK with zero bond radius)
A moving bias: steered-MD ● cd ../4_moving_bias ● The plumed input looks like g1-> 3 4 5 6 g1<- DISTANCE LIST 1 <g1> DISTANCE LIST 2 <g1> DISTANCE LIST 1 <g1> DIFFDIST 2 <g1> Steer the CV3 to 2.35 STEER CV 3 TO 2.35 VEL 2.50 KAPPA 0.1 Bohr, vel 2.50/kstep, PRINT W_STRIDE 1 kappa 0.1 Ha/bohr**2 ENDMETA ● A rather aggressive thermostatting should be used (ok the one we already set) ● Send the calculation mpirun -np 4 $FHIAIMS >out &
Recommend
More recommend
