Energy minimisation

Overview

Teaching: 10 min
Exercises: 0 min

Questions

• What should I do after generating my system’s topology and coordinates?

• How do I refine and remove bad contacts from the initial structure?

Objectives

• Run a sander command to minimise the initial structure.

• Explain briefly the options in the sander input file.

Before running any QM/MM simulation, we must ALWAYS equilibrate the system using only molecular mechanics (MM). This step will take approx. 20 min, so we are going to submit the calculation NOW on the short queue of ARCHER using the following script:

qsub send_mm_equil.pbs

We are going to minimise and equilibrate the system using the sander tool from AmberTools suite. We have provided commented input files and ARCHER submission scripts to make these steps easier to follow. You can find more information on the sander input options in Chapter 19 of the Amber20 manual.

To run sander, one needs to specify at least the initial coordinates (-c system.rst7) and the topology files of your system (-p system.parm7) and an input file (-i sander_min.in) providing the details of the simuilation to perform. Also, we can specify the name of output files such the log file (-o) and final coordinates (-r).

#!/bin/bash --login
#PBS -N MM_equil

# Select 1 nodes (maximum of 24 cores)
#PBS -q short
#PBS -l select=2
#PBS -l walltime=00:20:00

# Make sure you change this to your budget code
#PBS -A y14-amber2020

module swap PrgEnv-cray PrgEnv-gnu
module load amber-tools/20

# Move to directory that script was submitted from
export PBS_O_WORKDIR=$(readlink -f $PBS_O_WORKDIR)
cd $PBS_O_WORKDIR

date
echo "Minimisation"
aprun -n 8 sander.MPI -O -i sander_min.in -o min_classical.out -p system.parm7 -c system.rst7 -r system.min.rst7
date
echo "Thermalisation"
aprun -n 48 sander.MPI -O -i sander_heat.in -o heat_classical.out -p system.parm7 -c system.min.rst7 -r system.heat.rst7 -x system.heat.nc
date
echo "Pressure equilibration"
aprun -n 48 sander.MPI -O -i sander_equil.in -o equil_classical.out -p system.parm7 -c system.heat.rst7 -r system.equil.rst7 -x system.equil.nc
date

Energy minimisation with sander

We have devised a minimisation protocol that will perform 4000 steps using two different methods: 2000 of gradient descent and 2000 of conjugate gradient. The sander input file looks like this:

Minimisation of system 
 &cntrl
  imin=1,        ! Perform an energy minimization.
  maxcyc=4000,   ! The maximum number of cycles of minimization. 
  ncyc=2000,     ! The method will be switched from steepest descent to conjugate gradient after NCYC cycles.
 /

sander should have produced the following files:

• min_classical.out where it has written down the energy of the system every 50 steps.
• system.min.rst7 with the minimised coordinates.

If we open the min_classical.out file, the last step of the minimisation should look like this:

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
   4000      -1.7025E+05     2.0459E-01     9.4107E+00     OE1       169

 BOND    =    13362.4614  ANGLE   =      412.0010  DIHED      =     2157.5043
 VDWAALS =    34829.2930  EEL     =  -229750.9423  HBOND      =        0.0000
 1-4 VDW =      587.3195  1-4 EEL =     8150.8146  RESTRAINT  =        0.0000

We want the system to relax, so the ENERGY will decrease along the energy minimisation but we should also take the GMAX value into account. The GMAX is the slope on the Potential Energy Surface (PES) of the system of the current coordinates. GMAX should be close to 0. Other important information is the NAME and NUMBER, they point to the atom name and index, that is causing the major energy contribution.

Key Points

• After adding water molecules and combining coordinates of different system elements we have to minimise the system to fix any bad contacts (LEap had warned us already of some bad contacts in the structure).

• Removing bad contacts at this point will prevent our system from failing catastrophically later down the line.

• Combining different minimisation methods is a good practice and can help to avoid getting stuck into a local minima.