 |
 |
 |
| |
|
|
In this Masterclass, we will discuss how to run Hamiltonian replica exchange using PLUMED and GROMACS. We will also understand how to analyze the resulting trajectories. It is highly recommended to follow PLUMED Masterclass 21.5: Simulations with multiple replicas in advance.
Once you have completed this Masterclass you will be able to:
For this masterclass you will need versions of PLUMED and GROMACS that are compiled using the MPI library. In order to obtain the correct versions, please follow the instructions at this link.
The data needed to execute the exercises of this Masterclass can be found on GitHub. You can clone this repository locally on your machine using the following command:
git clone https://github.com/plumed/masterclass-22-10.git
Additional documentation about the replica exchange implementation discussed in this Masterclass can be found at this page: Using Hamiltonian replica exchange with GROMACS .
Throughout this Masterclass we will run simulations of alanine dipeptide in water using GROMACS and PLUMED. Whereas this system is too simple to be considered a proper benchmark for enhanced sampling methods, it is complex enough to be used in learning them. This Masterclass regards a specific implementation of Hamiltonian replica exchange that is only available when combining GROMACS and PLUMED.
Notice that simulations might take up to a couple of hours, depending on the hardware you have access to. The version of GROMACS that we provide on our conda channel is not optimized. Compiling GROMACS directly on your machine might lead to much better performance. In addition, using a GPU will also make your simulations significantly faster.
In PLUMED Masterclass 21.5: Simulations with multiple replicas we have learned how to run replica-exchange simulations where each replica was possibly feeling different biasing forces. A typical example is umbrella sampling with replica exchange, where each replica has a restraint located in a different position. Another example is bias-exchange metadynamics, where each replica is biased along a different collective variable. We have also seen parallel-tempering metadynamics, where replicas are kept at a different temperature and, at the same time, keep track of different history-dependent potentials. All these examples are a form of Hamiltonian replica exchange, since different replicas are subject to different potential energy functions, and thus different Hamiltonians.
In this masterclass we will consider a conceptually similar but technically different implementation of Hamiltonian replica exchange. Here, the different replicas will be simulated using different force field parameters. For most of the exercise, we will thus use an empty plumed.dat file (no bias added!). You will however need to include this file to ensure PLUMED can be enabled. We will also spend some time in learning how to edit the GROMACS topologies, so as to generate the modified force fields.
In order to use multiple-replica methods, you should run your simulation using MPI. This can be done prefixing your command with mpiexec -np N --oversubscribe, where N is the number of processes that you want to use and the --oversubscribe option is an OpenMPI option that is required to use more processes than the number of available processors. This is typically suboptimal, but we will need it in our lectures to run, e.g., simulations with 16 replicas even if we have a computer with 4 cores.
In brief, to run a GROMACS simulation where the individual replicas are in directories names dir0, dir1, etc and the (possibly empty) plumed.dat file is in the parent directory you will need a command such as
The option -replex 200 enables replica exchange and ensures exchanges are attempted every 200 steps. The option -hrex is specific for the implemenation discussed in this Masterclass, and informs GROMACS and PLUMED that there might be different force fields used in different replicas.
-hrex flag, no error will be issued, but the acceptance for exchanges will be incorrectly calculated.If you have random crashes on MacOS, try to set this environemnt variable:
export OMPI_MCA_btl="self,tcp"
In the root directory of this Masterclass you will find a conf.gro file, that can be used as a starting configuration for your simulations, a topol.top file, which contains the topology information, and a grompp.mdp file with reasonable simulation parameters. You have also a conf.pdb file, that is basically the conf.gro file converted to PDB format, and can be used for the MOLINFO keyword so as to facilitate atom selections in analysis.
Before swithing to solute tempering, we will play a bit with parallel tempering. Parallel tempering simulations require the highest replica to be hot enough for important energy barriers to be crossed in the simulation. We will first estimate how much we should raise the temperature so as to see forward and backward transitions between the two metastable states of our system.
To do so:
These are serial simulations, so you can run them with this command:
though it is recommended to use a script to run your simulations.
To analyze the resulting trajectories you can use PLUMED driver
with the following plumed.dat file:
MOLINFOSTRUCTURE=conf.pdb phi: TORSIONcompulsory keyword a file in pdb format containing a reference structure.ATOMS=@phi-2 PRINTthe four atoms involved in the torsional angleARG=phithe input for this action is the scalar output from one or more other actions.FILE=COLVARthe name of the file on which to output these quantities
You can then plot the second column of the resulting COLVAR file and see if there is a transition from the first basin (phi in range (-3,-1)) to the second basin (phi in range (0.5,1.5)), and back. Answer the following question:
The result will depend on stochastic factors, but also on the length of the simulation. I recommend running a 1 ns long simulation (500000 steps), but you can try with a longer or shorter trajectory.
Once you have identified this temperature, you can run a parallel tempering simulation. You will need a number of replicas to bridge from T=300 to the temperature you have identified. As a first guess, you can place them in a geometric series (that's the allocation that leads to uniform acceptance in a system with temperature-independent specific heat).
Let's say that you want to try with 16 replicas. You should create 16 directories named dir0, dir1, ... dir15 and place a file named topol.tpr in each. Each topol.tpr files will be created setting a different temperature in the grompp.mdp file. We will not need PLUMED for this exercise, so the simulation can be run with a command like this one:
mpiexec -np 32 --oversubscribe gmx_mpi mdrun -multidir dir? dir?? -replex 200 -nsteps 500000
The -replex option tells to GROMACS how frequently exchanges should be attempted. The average acceptance will be reported at the end of the md.log file.
The hottest replica will display transitions between the two metastable states. Thanks to the exchanges, both states will be observed also in the replica that is kept as T=300K.
Now answer the following questions:
In order to run Hamiltonian replica exchange simulations with multiple topologies, we will have first to learn how to generate the multiple topologies. Different tools could be used (including editing the topol.top file by hand!) but we will now learn how to use the partial_tempering tool available in PLUMED. This tool basically allows you to generate topologies where the energy of a subset of the atom have be scaled by a chosen factor. Remember that dividing the energy by a factor 2 is equivalent to multiplying the temperature by a factor 2.
Let's first have a look at the topol.top file. In this file there are a number of lines that look like #include ..., so that this file is not self-contained. The first thing that we have to do is to generate a self-contained topology file:
Have a look at the resulting processed.top file with a text editor. This file does not contain only information about alanine dipeptide, but also about the generic force-field parameters. It is thus self-contained.
Then you should edit the processed.top file to indicate which atoms you want to scale. To do so you have to add an underscore (_) to the atom name of the selected atoms. For instance, this line:
should be modified to this line:
To perform a solute tempering simulation, you should add the underscore to all the solute atoms (look in the [ atoms ] section). Once this is done you can use the following command:
This will scale the Hamiltonian of the selected atoms by a factor 1.0 (which means: no change!). Have a look at the resulting scaled.top file and find out what has changed. You can also try to apply two different scaling factors and check the difference:
Notice that these scaled topologies can be used to run GROMACS simulations. The parial_tempering script is far from perfect. We will now make some sanity check. To do these checks, we will either generate a new trajectory or simply take one of the trajectories that we generated in the previous exercise. We will then used a GROMACS tool named rerun, which allows recomputing the energy along a trajectory using a new topology.
First, create an energy file corresponding to the original topology:
The resulting ener.edr file can be converted to a text file with the gmx_mpi energy command.
Then, generate a new topology with scaling factor 1.0 and compute energies again:
The resulting energies should be identical!
As a second check, generate a new topology where you scaled all atoms (make sure to also select water!) by a factor 0.5:
This time, dihedral angles, LJ, and Coulomb energies should be multiplied by 0.5.
This is a very important check. The partial_tempering script is not compatible with some specific functional forms, e.g., CHARMM CMAP (search PLUMED mailing list for a fix). In case you are using an incompatible force field, you will find inconsistencies in these sanity checks. For the provided topol.top file, everything should work.
These tests are very important, please only proceed if you managed to pass them.
Another thing that we will have to check now is if acceptance is computed correctly. The code inserted in GROMACS to implement this calculation is non trivial and quite fragile. To do this check, you should run a short Hamiltonian replica exchange with two equivalent topology files. For instance, you can use the original topol.tpr file and the one that you obtained with scaling factor 1.0. Make sure that the two tpr files are using a different seed for randomizing the initial velocities or, even better, use two different conf.gro files to initialize the two simulations. For everything else, use the same settings you will use in production (ideally, same number of processes per replica, same GPU settings, etc).
Now run a short replica exchange simulation:
As written above, plumed.dat can be just an empty file. Check the resulting acceptance. Since the Hamiltonians are identical, the acceptance should be exactly 1.0.
Notice that this might be expected to fail if you use a GPU.**
You are now able to generate tpr files where the energy of a subset of the atoms is scaled. We will do something similar to Exercise 1a: Test with different temperatures, but playing with scaling factor instead of temperature. Make sure that you select all (and only) the atoms belonging to the solute (alanine dipeptide), so as to implement solute tempering. Try to run a set of 1 ns long simulations, with scaling factors decreasing (e.g., 1.0, 0.9, 0.8, etc). Now answer the following question:
Now run a Hamiltonian replica exchange simulation with multiple replicas, bridging from lambda=1.0 to the minimun value identified in the previous point. I would suggest using a linear distribution in lambda rather than a geometric one, but you can experiment. Similarly to parallel tempering, we can for now just analyze the reference replica (at lambda=1.0).
Now answer the following questions:
So far we only analized the reference replica (lambda=1.0). We can however do better and combine all replicas with WHAM. To this aim you should:
gmx_mpi trjcat -cat -f dir?/traj_comp.xtc dir??/traj_comp.xtc -o traj_multi.xtc).gmx_mpi mdrun -rerun traj_multi.xtc -s dir0/topol.tpr -e energy0.edr, gmx_mpi mdrun -rerun traj_multi.xtc -s dir1/topol.tpr -e ener1.edr, etc.).echo 11 | gmx_mpi energy -f ener0.edr -xvg no -e energy0.xvg -xvg no, etc.). 11 should select the potential energy.Then you can use the provided python wham script to compute the weight of all frames with the following python commands:
Now plot the weights array and answer the following questions:
The aim of this exercise is to optimize the list of values of lambda. As a target, we will try to make the product of the acceptances as large as possible. We will do this without running new simulations, i.e. just analyzing the simulation above.
This is a difficult exercise. It will be solved in the solution, but it is not necessary to complete it so as to proceed to the next step. I write some hints here.
sqrt(lambda). In other words, for each frame, you could write the energy as A*lambda+B*sqrt(lambda)+C, where A, B, and C are to be recomputed at each frame. Analyze the ene array that we generated in the last exercise so as to obtain these coefficients. The goal is to have a function that, for any frame, and for any value of lambda (including values not yet simulated!), returns the energy of the system.Once you have this function (let's call it predict_acceptance(lambda1,lambda2)) you can use the following script to find the ladder that maximizes the product of the acceptances. Let's say that the largest lambda is fixed to 1.0 and the lowest is fixed to lambda_min
You can then run a new simulation with the optimal lambdas and check if the predicted acceptances correspond to the observed ones!
As a very last point, repeat exercise Exercise 5: Parallel-tempering metadynamics using alanine dipeptide in explicit solvent and solute tempering with the lambda values optimized above. No further guidance is provided for this exercise, but it should not be difficult following the suggestions given for Exercise 5: Parallel-tempering metadynamics.
1.8.17