Using GMIN to generate endpoints (CHARMM)

We can use GMIN to explore the energy landscape of our peptide and find interesting structures (endpoints) to connect using OPTIM to form a pathway. After that, we can use PATHSAMPLE to calculate the rate of interconversion of these endpoints, taking into account multiple pathways.

Setting up CGMIN

Copy the metenk.crd file you have just created into a new directory, and rename it input.crd. CGMIN expects the input to always be called 'input.crd'.

mkdir CGMIN
cp metenk.crd CGMIN/input.crd
cd CGMIN

We now need to create a new file, data, which contains the GMIN keywords we would like to use. Copy the following into this file:

TARGET -157.082347
SLOPPYCONV 0.001
TIGHTCONV 0.00001
TRACKDATA
ACCEPTRATIO 0.3
EDIFF 0.01
DUMPINT 100
UPDATES 800
MAXIT 10000 10000
TEMPERATURE 0.8
STEPS 300 1.0
STEP 80.0 0.0
SAVE 10
CENTRE
CHPMAX  0.4
CHPMIN  0.2
CHNMAX  10
CHNMIN  0
CHARMMTYPE toph19_eef1_perm.inp param19_eef1_perm.inp
CHARMM
! Everything below the CHARMM line above is part of a CHARMM input file
set pardir "/home/csw34/svn/CHARMM31/toppar"

! BOMLev sets the level of warnings what do not cause the program to exit. -5 = very lax
BOMLev -5

! PRNLEV sets the ammount of output you get from CHARMM. 0 = small
PRNLEV 0

! Read standard topology and parameter files. The @top and ~par variables are set in the CHARMMTYPE line above
OPEN READ CARD UNIT 1 NAME @pardir/@top
READ RTF CARD UNIT 1
CLOSE UNIT 1

OPEN READ CARD UNIT 2 NAME @pardir/@par
READ PARAMETER CARD UNIT 2
CLOSE UNIT 2

! Generate the PSF for met-enk
READ SEQUence CARD
*
5
TYR GLY GLY PHE MET
GENErate FIRS NTER LAST CTER SETUp

! Read the initial coordinates from input.crd
OPEN UNIT 20 NAME input.crd READ CARD
READ COOR UNIT 20 CARD FREE
CLOSE UNIT 20

! Build the internal coordinate tables
IC FILL PRESERVE
IC PARAMETERS
IC PURGE
IC BUILD

! Set up the EEF1 solvent model
eef1 setup temp 298.15 unit 93 name "/home/csw34/svn/CHARMM31/toppar/solvpar.inp"
update ctonnb 7. ctofnb 9. cutnb 15. group rdie

As before, you may need to make some changes to this file, so that the paths point correctly to your CHARMM files. If you are working in CUC3 in Cambridge, this may just be as simple as changing my CRSID (csw34) for yours.

Everything above the CHARMM line is a GMIN keyword. It is highly recommended that you take a look at the GMIN documentation here, so that you know what each of them does. In brief, we are doing up to 300 basin-hopping steps (or until the TARGET - the global minimum - is found), and at each step generating new coordinates by changing between 0 and 10 dihedral angles (backbone or sidechain) by up to 80 degrees.

Everything below the CHARMM line is the CHARMM input that we use to set up the potential. If you'd like to know what this is doing, consult the CHARMM documentation! You may notice that is is very similar to the input we used to generate the metenk.crd file in the previous part of this tutorial.

Running CGMIN

Assuming that you have already compiled it - it's time to run CGMIN and find some minima!

CGMIN > charmm.out &

This should take about a minute to run, it's pretty quick! Once it is done, you should have quite a few output files in the directory. There are a few important ones you might want to take a look at:

GMIN_out

This is the primary GMIN output file, containing the energies of the minima found after each basin-hopping step, as well as the final 10 that have been saved in the dbase.x and dbase.x.pdb files. For example,

Qu          2 E=     -144.0873506553 steps=  518 RMS= 0.90078E-03 Markov E=    -143.4215906     t=        0.5

This is the line summarising the result of the second 'quench' (basin-hopping step). The energy of the minimum found after 518 L-BFGS minimisation steps is -144.0873506553, which is them compared to the current markov energy of -143.4215906. As this energy is lower, this minimum will replace the current markov minimum, and the next step will start from this conformation.

Final Quench      1 energy=     -157.0823478413 steps=   96 RMS force=      0.0000085 time=       28.85

This line shows the final quench (tight convergence) for the lowest energy minimum found, in this case - it is the global minimum. We know this from previous work, which is why we added the TARGET line to the data' file. The coordinates of this minimum are saved in dbase.1 in CHARMM .crd (CARD) format, and dbase.1.pdb in PDB format.

dbase.x and dbase.x.pdb

These files contain the coordinates of the lowest SAVE(GMIN keyword) minima found during the basin-hopping run. They are numbered in order of increasing energy, so dbase.1 and dbase.1.pdb contain the coordinates of the lowest energy minimum found in CHARMM .crd (CARD) and PDB format respectively. It is always worth checking that these structures are reasonable, so that you know nothing is really broken! For example, you could compare the lowest energy minimum in dbase.1.pdb with the minimised starting structure in initialmin.pdb. You can do this using VMD, or Pymol:

vmd -pdb initialmin.pdb -pdb dbase.1.pdb

This will load both sets of coordinates into a single molecule, that you can change between using the slider as if it was a trajectory. Frame 0 is the initial structure (initialmin.pdb) and frame 1 is the global minimum. You should see quite a difference!

NOTE: if you are working on a cluster, you might need to load the vmd module first using:

module load vmd

If you'd like to make sure that you have done this part right, you can find a complete set of output files here.

Now that we have a set of low energy minima, and a high energy extended structure, we can use OPTIM to connect them with pathways.