OLD Minimization

From CHARMM tutorial
Jump to navigationJump to search

Description of Minimization

Description

Minimization reduces the energy of the system to the lowest possible point, simulating the system's energy as absolute zero.

Minimization algorithms

There are different methods used for calculating the lowest-energy conformation of a given structure. These descriptions here were taken from minimiz.doc, which is one of the standard documentation files included with CHARMM. The SD and ABNR methods are the most commonly utilized and will therefore be used in this tutorial when making minimization input scripts.

Steepest Descent (SD)
The simplest minimization algorithm is steepest descent (SD). In each step of this iterative procedure, the coordinates are adjusted in the negative direction of the gradient. It has one adjustable parameter, the step size, which determines how far to shift the coordinates at each step. The step size is adjusted depending on whether a step results in a lower energy. i.e., if the energy drops, the step size is increased by 20% to accelerate the convergence. If the energy rises, overshooting a minimum, the step size is halved. Steepest descent does not converge in general, but it will rapidly improve a very poor conformation.

Conjugate Gradient Technique (CONJ)
A second method is the conjugate gradient technique (CONJ) which has better convergence characteristics (R. Fletcher & C.M. Reeves, The Computer Journal 7:149 (1964)). The method is iterative and makes use of the previous history of minimization steps as well as the current gradient to determine the next step. It can be shown that the method converges to the minimum energy in N steps for a quadratic energy surface where N is the number of degrees of freedom in the energy. Since several terms in the potential are quadratic, it requires less evaluations of the energy and gradient to achieve the same reduction in energy in comparison to steepest descent. Its main drawback is that with very poor conformations, it is more likely to generate numerical overflows than steepest descent. The algorithm used in CHARMM has a slightly better interpolation scheme and automatic step size selection.

Conjugate Gradient Powell Method Minimizer (POWE)
A third method is the conjugate gradient powell method minimizer (POWE) (A. Brunger). Its efficiency is much improved over the Fletcher and Reeves method. The use of the POWELL minimizer is encouraged whenever ABNR is not possible. The POWELL minimizer has no INBFRQ or IHBFRQ feature. The use of CHARMM loops to mimic this important feature is encouraged. The CHARMM loop facilities allow harmonic constrained minimizations with periodic updates. In case of bad contacts or unlikely conformations the SHAKE method might give an error when the displacements become to large. Using a harmonic constraint setup with periodic updates solves this problem.

Newton-Raphson minimization (NRAP)
A fourth method involves solving the Newton-Raphson minimization equations iteratively (NRAP). This procedure requires the computation of the derivative of the gradient which is a matrix of size O(n2). The procedure here is to find a point where the gradient will be zero (hopefully a minimum in energy) assuming that the potential is quadratic. The Newton-Raphson equations can be solved by a number of means, but the method adopted for CHARMM involves diagonalizing the second derivative matrix and then finding the optimum step size along each eigenvector. When there are one or more negative eigenvalues, a blind application of the equations will find a saddle point in the potential. To overcome this problem, a single additional energy and gradient determination is performed along the eigenvector displacement for each small or negative eigenvalue. From this additional data, the energy function is approximated by a cubic potential and the step size that minimizes this function is adopted. The advantages of this method are that the geometry cannot remain at a saddle point, as sometimes occurs with the previous procedures, and that the procedure converges rapidly when the potential is nearly quadratic (or cubic). The major disadvantage is that this procedure can require excessive storage requirements, O(n2), and computation time, O(n3), for large molecules. Thus we are currently restricted to systems with about 200 atoms or less for this method. IMAGES and SHAKE are currently unavailable with this method.

Adopted Basis Newton-Raphson Method (ABNR)
The fifth method available is an adopted basis Newton-Raphson method (ABNR) (D. J. States). This routine performs energy minimization using a Newton-Raphson algorithm applied to a subspace of the coordinate vector spanned by the displacement coordinates of the last (MINDim) positions. The second derivative matrix is constructed numerically from the change in the gradient vectors, and is inverted by an eigenvector analysis allowing the routine to recognize and avoid saddle points in the energy surface. At each step the residual gradient vector is calculated and used to add a steepest descent step onto the Newton-Raphson step, incorporating new direction into the basis set. This method is the best for most circumstances.

Truncated-Newton Minimization Package (TN)
The sixth method is the truncated-Newton (TN) minimization package TNPACK, developed by T. Schlick and A. Fogelson. TNPACK is based on the preconditioned linear conjugate-gradient technique for solving the Newton equations. The structure of the problem --- sparsity of the Hessian --- is exploited for preconditioning. Thorough experience with the new version of TNPACK in CHARMM has been described in the paper: Journal of Computational Chemistry: 15, 532--552, 1994. Through comparisons among the minimization algorithms available in CHARMM, we find that TNPACK compares favorably to ABNR in terms of CPU time when curvature information is calculated by a finite-difference of gradients (the "NUMEric" option of TNPACK). With the analytic option, TNPACK can converge more rapidly than ABNR for small and medium systems (up to 400 atoms) as well as large molecules that have reasonably good starting conformations; for large systems that are poorly relaxed (i.e., the initial Brookhaven Protein Data Bank structures are poor approximations to the minimum), TNPACK performs similarly to ABNR. SHAKe is currently unavailable with this method.

Basic Minimization Script

Please read the Basic Input Script section if you have had no prior experience in creating a basic CHARMM input script.

Preparation
This tutorial for making a basic minimization input script continues from the Basic Input Script tutorial. The protein with the PDB.org ID of 1YJP will be used in this part of the tutorial. The PDB file and an associated water box are available from the main tutorial web site http://www.charmmtutorial.org/http://www.charmmtutorial.org. should have been combined into a single structure before minimization is attempted.

Create an input script, reading in the RTF and Parameter files. After this is done, the user must read in the appended PSF and CRD files that were created by the basic input script. This can be done with the following commands:

OPEN READ UNIT 2 CARD NAME 1yjp-full.psf
READ PSF CARD UNIT 2
CLOSe UNIT 2

OPEN UNIT 2 READ CARD NAME 1yjp-full.crd
READ COOR UNIT 2 CARD
CLOSe UNIT 2

CHARMM has the ability to work with and compare different coordinate sets. this tutorial, the root mean squared (RMS) difference between the structure before and after minimization will be computed. it is necessary to store the initial coordinates by copying them to the comparison set. "COOR" signifies that the user will retrieve the coordinates and "COPY" tells CHARMM what the users intent is with the coordinates. "COMP" is short for comparison set, so CHARMM knows the user wants to copy the coordinates to the comparison set. The line should look like:

COOR COPY COMP

SHAKe and non-bonded parameters
It is a good idea to fix hydrogen bond lengths in place to remove high frequency motion from the system, allowing the use of 1 fs time steps. can be accomplished with the "SHAKe BONH PARAmeters" command.

It is also necessary to set the non-bonded parameters to the same values that will be used for dynamics. full discussion of these parameters is given in the dynamics section.

SHAKe BONH PARAmeters
NBONd INBFrq -1 ELEC FSWItch VDW VSWItch CUTNb 16. CTOFnb 12. CTONnb 8.

Applying minimization
To tell CHARMM that a minimization is requested, start the command with "MINI". Next, choose what algorithm to use. This tutorial will first use steepest descent since for many structures it is likely to provide a substantial energy loss in a relatively short period of time. To inform CHARMM of this decision follow "MINI" with "SD". specify the number of steps, the "NSTEp" option is used followed by the number of steps. The line should look like the following:

MINI SD NSTEp 100

To use the ABNR algorithm, do the same thing as above except replace "SD" with "ABNR". The number of steps can be set to 100 the same way it was for steepest descent. CHARMM can print a report on minimization so with every set number of steps, CHARMM will print a notice in the output file. This can be done with the "NPRInt" command followed by the interval in which a notice should get printed.

With ABNR, minimization can end early if the root mean square of the gradient (GRMS) reaches a certain number. The smaller the number it is, the longer the minimization will run (unless it is terminated early by "NSTEp"), however a larger number will result in a less precise minimization. For this tutorial, the relatively large tolerance of 0.05 is used. the number of print intervals, use the command "TOLGrad" followed by "0.05". The next line should now look like:

MINI ABNR NSTEp 100 NPRInt 100 TOLG 0.05

To compare the RMS value to one previous to minimization, use the line:

COOR RMS

CHARMM will print the RMS between the main coordinate set (which now contains the minimized structure) and the comparison set (which still has the original coordinates).

To end the minimization process, write out a PDB file containing the new minimized structure, and then terminate the CHARMM script.

The final input file should resemble the following:

 * Minimize PDB
 *
 
 ! Read in Topology and  Parameter files
 OPEN UNIT 1 CARD READ NAME top_all27_prot_na.rtf
 READ RTF CARD UNIT 1
 CLOSE UNIT 1
 
 OPEN UNIT 1 CARD READ NAME par_all27_prot_na.prm
 READ PARA CARD UNIT 1
 CLOSE UNIT 1
 
 ! Read sequence from the PDB coordinate file
OPEN READ UNIT 27 CARD NAME 1yjp-full.psf
READ PSF CARD UNIT 27
CLOSE UNIT 27

OPEN READ UNIT 10 CARD NAME 1yjp-full.crd
READ COOR CARD UNIT 10
CLOSE UNIT 10

! constrain hydrogen bond lengths with SHAKe
SHAKe BONH PARAmeters

! non-bond setup (same values as will be used in dynamics)
 NBONd INBFrq -1 ELEC FSWItch VDW VSWItch CUTNb 16. CTOFnb 12. -
   CTONnb 8.
 
 ENERgy
 
 COOR COPY COMP
 MINI SD NSTEp 100
 MINI ABNR NSTEp 1000 NPRInt 100 TOLG 0.05
 COOR RMS
 
 OPEN WRITE UNIT 1 CARD NAME 1yjp-full-min.pdb
 WRITE COOR PDB UNIT 1
 * new_1yjp-11213-min.pdb
 *
 
 OPEN WRITE UNIT 1 CARD NAME 1yjp-full-min.crd
 WRITE COOR CARD UNIT 1
 * new_1yjp-11213-min.crd
 *
 
 STOP