PLUMED microtutorial Log into the machines Copy the plumed tutorial - - PowerPoint PPT Presentation

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.

• utput_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 |- GROUP MEMBERS: 3 4 5 6 1-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 |- 1st SET MEMBERS: 1 |- 2nd SET MEMBERS: 3 4 5 6 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) DISTANCE LIST 1 <g1> g1-> 3 4 5 6 g1<- DISTANCE LIST 2 <g1> PRINT W_STRIDE 2 ENDMETA

PLUMED.OUT

::::::::::::::::: READING PLUMED INPUT ::::::::::::::::: |- GROUP FOUND: g1 |- GROUP MEMBERS: 3 4 5 6 1-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 |- 1st SET MEMBERS: 1 |- 2nd SET MEMBERS: 3 4 5 6 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) DISTANCE LIST 1 <g1> g1-> 3 4 5 6 g1<- DISTANCE LIST 2 <g1> PRINT W_STRIDE 2 ENDMETA

PLUMED.OUT

::::::::::::::::: READING PLUMED INPUT ::::::::::::::::: |- GROUP FOUND: g1 |- GROUP MEMBERS: 3 4 5 6 1-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 |- 1st SET MEMBERS: 1 |- 2nd SET MEMBERS: 3 4 5 6 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) DISTANCE LIST 1 <g1> g1-> 3 4 5 6 g1<- DISTANCE LIST 2 <g1> PRINT W_STRIDE 2 ENDMETA

PLUMED.OUT

::::::::::::::::: READING PLUMED INPUT ::::::::::::::::: |- GROUP FOUND: g1 |- GROUP MEMBERS: 3 4 5 6 1-DISTANCE: (1st SET: 1 ATOMS), (2nd SET: 4 ATOMS); PBC ON |- DISCARDING DISTANCE COMPONENTS (XYZ): 000 |- 1st SET MEMBERS: 1 |- 2nd SET MEMBERS: 3 4 5 6 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) DISTANCE LIST 1 <g1> g1-> 3 4 5 6 g1<- DISTANCE LIST 2 <g1> PRINT W_STRIDE 2 ENDMETA

COLVAR

#! FIELDS time cv1 cv2 0.0040 3.641223429 5.668572649 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 Tells you which fields are there Time CV1 CV2

• 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 DISTANCE LIST 1 <g1> DIFFDIST 2 <g1> UMBRELLA CV 3 AT 0.0 KAPPA 0.01 PRINT W_STRIDE 1 ENDMETA

Use diff of distances Use diff of distances Impose a 0-centered spring on CV3

Run the calc

#! 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 Time CV1 (dist 1) CV2 (dist 2) CV3 (dist1- dist 2) Bias from Harmonic spring Center of the bias Work The COLVAR file is:

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 CV 3 TO 2.35 VEL 2.50 KAPPA 0.1 PRINT W_STRIDE 1 ENDMETA Steer the CV3 to 2.35 Bohr, vel 2.50/kstep, kappa 0.1 Ha/bohr**2

• A rather aggressive thermostatting should be

used (ok the one we already set)

• Send the calculation

mpirun -np 4 $FHIAIMS >out &

• utputs
• COLVAR
• PLUMED.OUT

|-STEERING COLVAR 3 TO 2.350000: VELOCITY=2.500000 cvunit/kstep, SPRING=0.100000 #! FIELDS time cv1 cv2 cv3 XX XX RST3 WORK3 0.0020 3.604130111 5.694642 -2.090512 RST 3 -2.090512821 0.000000000 0.0040 3.640674604 5.669551 -2.028877 RST 3 -2.088012821 -0.000011548 0.0060 3.652244328 5.662087 -2.009843 RST 3 -2.085512821 -0.000028398 0.0080 3.658914636 5.658905 -1.999991 RST 3 -2.083012821 -0.000048235 0.0100 3.659760130 5.660854 -2.001094 RST 3 -2.080512821 -0.000068540 0.0120 3.655008550 5.667654 -2.012645 RST 3 -2.078012821 -0.000086638 Moving spring center Work is accumulated

Output: COLVAR

• In gnuplot:

plot “COLVAR” u 1:4,”” u 1:7

Output: work

• In gnuplot:

plot “COLVAR” u 1:4,”” u 1:7

The work is not the free energy!

• Work and free energy are related through

Jarzynski equation

• They are very similar whenever the orthogonal

degrees of freedom are few and rapidly averaged (precisely this case here)

Metadynamics (homework)

• Here we will test metadynamics in 2 dimensions

assuming the two distances as independent

• Metadynamics requires the hills width for the

two CVs-> we take them from the first run (dir 2_how_to_enable_plumed)

• Make the distributions

grep -v FIELDS COLVAR.lda | awk '{print $1,$2}' | awk -f ./distribution.awk >distrib_cv1 grep -v FIELDS COLVAR.lda | awk '{print $1,$3}' | awk -f ./distribution.awk >distrib_cv2

Fit the widths in gnuplot

f(x)=a*exp(-b*(x-c)*(x-c)) g(x)=d*exp(-e*(x-f)*(x-f)) a = 0.06 b = 40.0 c = 3.6 d = 0.05 e = 2.5 f = 6.0 fit g(x) "distrib_cv2" via d,e,f fit f(x) "distrib_cv1" via a,b,c print 1/(0.5*b)**(0.5) print 1/(0.5*e)**(0.5) Plot “distrib_cv1” w lp ,f(x) w lp, “distrib_cv2” w lp , g(x) w lp Sigma for first cv Sigma for second cv

The distributions

Metadynamics input

• Calculations in 2d require lots of time. Here at

least 50 ps (i.e. 25K steps with 2fs ts).

DISTANCE LIST 1 <g1> SIGMA 0.21 g1-> 3 4 5 6 g1<- UWALL CV 1 LIMIT 7. KAPPA 0.5 DISTANCE LIST 2 <g1> SIGMA 0.21 UWALL CV 2 LIMIT 7. KAPPA 0.5 HILLS HEIGHT 0.00047 W_STRIDE 50 PRINT W_STRIDE 2 ENDMETA Each CV has a sigma Metadynamics params Confining restraining corral

Metadynamics input

• Calculations in 2d require lots of time. Here at

least 50 ps (i.e. 25K steps with 2fs ts).

DISTANCE LIST 1 <g1> SIGMA 0.21 g1-> 3 4 5 6 g1<- UWALL CV 1 LIMIT 7. KAPPA 0.5 DISTANCE LIST 2 <g1> SIGMA 0.21 UWALL CV 2 LIMIT 7. KAPPA 0.5 HILLS HEIGHT 0.00047 W_STRIDE 50 PRINT W_STRIDE 2 ENDMETA Each CV has a sigma Metadynamics params Confining restraining corral

Metadynamics: output

• COLVAR
• HILLS

0.100 3.595283559 5.277905326 0.210 0.210 0.000470000 0.000 0.200 3.755987525 5.732841610 0.210 0.210 0.000470000 0.000 0.300 3.515045105 5.671404412 0.210 0.210 0.000470000 0.000 time cv1 cv2 Sigma cv1 Sigma cv2 height #! FIELDS time cv1 cv2 vbias vwall 0.1000 3.595283559 5.277905326 0.000470000 0.000000000 0.1040 3.548153486 5.307082166 0.000453911 0.000000000 0.1080 3.521752306 5.318650588 0.000433813 0.000000000 0.1120 3.517271798 5.316846599 0.000431188 0.000000000 0.1160 3.529824480 5.311633765 0.000441977 0.000000000

Postprocess metadynamics output

• Use sum_hills tool

../sum_hills/sum_hills.x -NDIM 2 -NDW 1 2 -file HILLS

and get a fes.dat

• In gnuplot:

set pm3d splot “fes.dat” w pm3d (evtl: set contour surf, unset clabel, set cntrparams level incremental -10,1.,0 )

...you should have something like this

U_WALL LIMIT 7.0 U _ W A L L L I M I T 7 .

A bit more on the convergence of metadynamics

• One should have many RECROSSINGS

And average many free energy realizations OR use well-tempered metadynamics

• Many other information, example,

documentation on www.plumed-code.org