diff --git a/user-doc/figs/trieste-5-colvars.png b/user-doc/figs/trieste-5-colvars.png new file mode 100644 index 0000000000000000000000000000000000000000..46efc381a2b2924844fb3734702fffd38fb3eb20 Binary files /dev/null and b/user-doc/figs/trieste-5-colvars.png differ diff --git a/user-doc/figs/trieste-5-demux-bad.png b/user-doc/figs/trieste-5-demux-bad.png new file mode 100644 index 0000000000000000000000000000000000000000..d4931f7e3be6e828bdaecc6947d3414d23baf2c2 Binary files /dev/null and b/user-doc/figs/trieste-5-demux-bad.png differ diff --git a/user-doc/figs/trieste-5-demux-good.png b/user-doc/figs/trieste-5-demux-good.png new file mode 100644 index 0000000000000000000000000000000000000000..b0a41fcb88f733c1509a045d47fff3f6e0465903 Binary files /dev/null and b/user-doc/figs/trieste-5-demux-good.png differ diff --git a/user-doc/tutorials/a-trieste-5.txt b/user-doc/tutorials/a-trieste-5.txt index 082c79e5e63374a8dfad0767ef19f49f90f93992..be0a2eff082a476bbee6f7c6722eeb120fd7dcfd 100644 --- a/user-doc/tutorials/a-trieste-5.txt +++ b/user-doc/tutorials/a-trieste-5.txt @@ -29,6 +29,9 @@ advantage of a special syntax for setting up multi-replica simulations that is o since version 2.4. Exercizes could be done also with version 2.3 but using a different syntax with respect to the one suggested. +Also notice that in the `.solutions` directory of the tarball you will find correct input files. +Please only look at these files after you have tried to solve the problems yourself. + \section trieste-5-introduction Introduction So far we always used PLUMED to run one simulation at a time. However, PLUMED can also @@ -414,9 +417,9 @@ two free-energy wells that correspond roughly to positive and negative \f$\phi\f This will tell you the unbiased population of the two minima. You should expect values close to these ones \verbatim -0.978318 0.0216822 +0.968885 0.0311148 \endverbatim -that tells you that the first minimum has a population of roughly 98%. +that tells you that the first minimum has a population of roughly 97%. \note This approach can be used to reweight any replica exchange simulation, irrespectively of how @@ -440,9 +443,33 @@ but using only three replicas: - two unbiased replicas. - one biased replica along psi. -Compute as before the relative population of the two minima. Which is the result? What can you learn from this example? +Create a file `plumed3.dat` aimed at this. If you are not able you can find a working one in the `solutions` directory. +Then use the following commands similarly to before: +\verbatim +> mpirun -np 3 gmx_mpi mdrun -plumed plumed3.dat -nsteps 1000000 -replex 120 -multi 3 +> gmx_mpi trjcat -f traj_comp{0,1,2}.xtc -cat -o fulltraj.xtc +> mpirun -np 3 plumed driver --mf_xtc fulltraj.xtc --plumed plumed3_reweight.dat --multi 3 +> paste COLVAR_WHAM.{0,1,2} | grep -v \# | awk '{for(i=1;i<=NF/2;i++) printf($(i*2)" ");printf("\n");}' > ALLBIAS +> python3 SCRIPTS/wham.py ALLBIAS 3 2.5 +\endverbatim + +Compute as before the relative population of the two minima. +\verbatim +> paste <(grep -v \# COLVAR_CONCAT.0) weights.dat | awk '{if($2>0)wb+=$NF; else wa+=$NF}END{print wa,wb}' +\endverbatim +Which is the result? What can you learn from this example? + +If you did things correclty, here you should find a population close to 2/3 for the first minimum and 1/3 for the second one. +Notice that this population just reflects the way we initialized the simulations (2 of them in one minimum and 1 in the other minimum) +and does not take into account which is the relative stability of the two minima according to the real potential energy function. +In other words, we are just obtaining the relative population that we used to initialize the simulation. + +Also plot the colvar files and look at them with attention. An excerpt is shown below. -Also plot the colvar files and look at them with attention. How many transitions between the two free-energy wells can you observe? +\image html trieste-5-colvars.png "Time series of phi from the simulations with a missing variable (exercise 3). Notice that there are appaernt transitions, but they are +always due to accepted exchanges between replicas." + +How many transitions between the two free-energy wells can you observe? Remember that replica exchange involves coordinate swaps that do not correspond to the real system dynamics. It is very useful to look at "demuxed" trajectories. @@ -464,9 +491,17 @@ continuous ones with the following commands You will now find new trajectories `0_trajout.xtc`, `1_trajout.xtc`, `2_trajout.xtc`, and `3_trajout.xtc`. These files can be in turn analyzed with `plumed driver`. Try to recompute collective variables on the trajectories -obtained in \ref trieste-5-ex-1 and with those obtained in \ref trieste-5-ex-3. Can you spot out the difference between -the two cases? +obtained in \ref trieste-5-ex-1 and with those obtained in \ref trieste-5-ex-3. + +\image html trieste-5-demux-good.png "Time series of phi obtained from the demuxed trajectories from exercise 2" +\image html trieste-5-demux-bad.png "Time series of phi obtained from the demuxed trajectories from exercise 3" + +As you can see, the trajectories from \ref trieste-5-ex-1 show a number of transitions. On the other hand, the trajectories from \ref trieste-5-ex-3 +are stuck. +As a general rule, notice that there is no way to estimate correctly the relative population +of two metastable states if there is not a continuous "demuxed" trajectory joining them, +possibly exhibiting multiple transitions. \section trieste-5-conclusion Conclusions diff --git a/user-doc/tutorials/trieste-5/.solutions/driver.dat b/user-doc/tutorials/trieste-5/.solutions/driver.dat new file mode 100644 index 0000000000000000000000000000000000000000..565bb19857680b0cbb25cfc5dda2278d30c7ec91 --- /dev/null +++ b/user-doc/tutorials/trieste-5/.solutions/driver.dat @@ -0,0 +1,9 @@ +# load the pdb to better select torsions later +MOLINFO STRUCTURE=reference.pdb +# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C +phi: TORSION ATOMS=@phi-2 +# Compute the backbone dihedral angle psi, defined by atoms N-CA-C-N +psi: TORSION ATOMS=@psi-2 + +# Just compute the collective variables +PRINT ARG=phi,psi FILE=COLVAR_DEMUX diff --git a/user-doc/tutorials/trieste-5/.solutions/plumed2_reweight.dat b/user-doc/tutorials/trieste-5/.solutions/plumed2_reweight.dat index d729a1a14475dc3872b2bdd6579f759cbc5c46b0..3b9ab7647a31c6ce6d0b600c5ff975cac30bd7d4 100644 --- a/user-doc/tutorials/trieste-5/.solutions/plumed2_reweight.dat +++ b/user-doc/tutorials/trieste-5/.solutions/plumed2_reweight.dat @@ -24,7 +24,7 @@ metad: METAD ... ... # Print both collective variables and the value of the bias potential on COLVAR_CONCAT file -PRINT ARG=phi,psi FILE=COLVAR_CONCAT STRIDE=10 +PRINT ARG=phi,psi FILE=COLVAR_CONCAT # Then write the bias potential (as we did for reweighting) # As a good practice, remember that there might be multiple bias diff --git a/user-doc/tutorials/trieste-5/.solutions/plumed3.dat b/user-doc/tutorials/trieste-5/.solutions/plumed3.dat new file mode 100644 index 0000000000000000000000000000000000000000..ef90a0b4ae53b9239fdd6a07acb1c50757fe67a0 --- /dev/null +++ b/user-doc/tutorials/trieste-5/.solutions/plumed3.dat @@ -0,0 +1,24 @@ +# load the pdb to better select torsions later +MOLINFO STRUCTURE=reference.pdb +# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C +phi: TORSION ATOMS=@phi-2 +# Compute the backbone dihedral angle psi, defined by atoms N-CA-C-N +psi: TORSION ATOMS=@psi-2 + +metad: METAD ... +# Activate well-tempered metadynamics +# Deposit a Gaussian every 500 time steps, with initial height equal +# to 1.2 kJoule/mol, biasfactor equal to 10.0 +# Remember: replica 0 and 1: no bias +# replica 2, bias on psi + ARG=@replicas:{phi,psi,psi} + HEIGHT=@replicas:{0.0,0.0,1.2} # make sure that replicas 0 and 1 feel no potential! + PACE=500 + BIASFACTOR=10.0 + SIGMA=0.35 +# Gaussians will be written to file and also stored on grid + FILE=HILLS GRID_MIN=-pi GRID_MAX=pi +... + +# Print both collective variables and the value of the bias potential on COLVAR file +PRINT ARG=phi,psi FILE=COLVAR STRIDE=10 diff --git a/user-doc/tutorials/trieste-5/.solutions/plumed3_reweight.dat b/user-doc/tutorials/trieste-5/.solutions/plumed3_reweight.dat new file mode 100644 index 0000000000000000000000000000000000000000..cdd52a5f90ca428ca4c0c30f96150e95a07ce322 --- /dev/null +++ b/user-doc/tutorials/trieste-5/.solutions/plumed3_reweight.dat @@ -0,0 +1,35 @@ +# load the pdb to better select torsions later +MOLINFO STRUCTURE=reference.pdb +# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C +phi: TORSION ATOMS=@phi-2 +# Compute the backbone dihedral angle psi, defined by atoms N-CA-C-N +psi: TORSION ATOMS=@psi-2 + +metad: METAD ... +# Activate well-tempered metadynamics +# Deposit a Gaussian every 500 time steps, with initial height equal +# to 1.2 kJoule/mol, biasfactor equal to 10.0 +# Remember: replica 0 and 1: no bias +# replica 2, bias on psi + ARG=@replicas:{phi,psi,psi} + HEIGHT=0 +# PACE=500 + BIASFACTOR=10.0 + SIGMA=0.35 + FILE=HILLS GRID_MIN=-pi GRID_MAX=pi + RESTART=YES + TEMP=300 + PACE=100000000 +... + +# Print both collective variables and the value of the bias potential on COLVAR_CONCAT file +PRINT ARG=phi,psi FILE=COLVAR_CONCAT + +# Then write the bias potential (as we did for reweighting) +# As a good practice, remember that there might be multiple bias +# potential applied (e.g. METAD and a RESTRAINT). +# to be sure that you include all of them, just combine them: +bias: COMBINE ARG=*.bias PERIODIC=NO +# here omit STRIDE, so that you will print the bias +# for all the frames in your concatenated trajectory +PRINT FILE=COLVAR_WHAM ARG=bias