Free Energy Perturbation Calculations
The PERTurbe command allows the scaling of energy between PSFs for use in
energy analysis, comparisons, slow growth free energy simulations,
widowing free energy simulation, and for slow growth homology modelling.
This is a rather flexible implementation of free energy perturbation
that allows connectivity to change. Also, three energy restraint
terms (harmonic, dihedral and NOE) and the SKIP command flags are subject
to change which allows a flexible way in which to compute free energy
differences between different conformations. This code in implemented
in a non-intrusive manner and works with all minimizers and
integrators. SHAKE can now be applied to bonds which are mutated as
well; an appropriate constraint corrections is calculated
automatically in these cases.
* Menu:
* Syntax:: Syntax of PERT Commands
* Description:: Description of PERT Commands
* Restrictions:: Restrictions in PERT Command usage
* References:: References regarding Free Energy Perturbations
* Examples:: A Sample Input Files
* Constraints:: Special considerations regarding SHAKE
Syntax of Free Energy Perturbation Commands
[Syntax PERT]
PERTurb [OFF] [INBFrq int nonbond-specs] [RESEt] atom-selection
[INBFrq 0 ]
The PERT OFF command disables the free energy routines and the current
(lambda=1) PSF is used for subsequent commands.
When OFF is not specified, this command saves the current PSF as the lambda=0
state. The atom-selection indicated which atoms have changed. This is to
make the calculation run more efficiently. If only a small percent of the
atoms have changed, this doubles the performance. The nonbond specs are
included to make sure the nonbond exclusion lists are properly setup. This
then allows the connectivity to change during the simulation. INBFrq should
not be set to zero here unless the exclusion lists have already been setup
in a previous command.
----------------------------------------------------------------------------
ENERgy ... [ RESET ] [ free-energy-step-spec ]
DYNAmics ... [ PUNIt integer ] [WHAM integer]
MINImize ...
[ RESET ] ! Resets all all accumulation data and counters.
(automatic for the first PERT or after a PERT OFF command)
free-energy-step-spec::=
[PWINdow [LAMBda real] ] [PSTArt int] [PSTOp int] [LSTArt real] [LSTOp real] -
[PSLOwgrowth ]
[PINCrement int] [PEQUilibrate int] [LAVErage] [LINCrement real]
[PWINdow ] ! specifies the windowing algorithm (default)
[PSLOwgrowth] ! specifies the slow growth algorithm
[LAMBda real] ! specifies the lambda value for windowing methods or for
energy or minimization calculations.
[PSTArt int] ! starting dynamics step number for accumulation (default 1)
[PSTOp int] ! ending dynamics step number for accumulation (default 0)
[LSTArt real] ! specifies the starting lambda value (default 0.0)
[LSTOp real] ! specifies the final lambda value (default 1.0)
[PINCrement int] ! specifies number of steps to next window (auto mode).
[PEQUil int] ! specifies number of steps for equilibration (auto mode).
[LAVErage] ! specifies that lambda = (LSTART+LSTOP)/2 (auto mode).
[LINCrement] ! Specifies the lambda increment between windows (auto mode).
[ END ] ! Turns off the free energy code
The PSTArt and PSTOp values are relative to the number of dynamics steps
since PERT command was first enabled, or if a PERT RESET command is used
(regardless of whether the DYNAmics commands is run more than once or
whether the dynamic run involved the use of restart files).
By specifying the auto mode parameters (PINCrement, PEQUilibrate, LINCrement),
a new window will start at the conclusion of the current window with modified
parameters. Also, in auto mode, the run will terminate at the end of a window
where the LSTOP value is 1.0 or 0.0.
The PUNIt option allows the free-energy-step-spec to be specified
more than once and acts as a scheduler for a particular simulation.
The format of the PUNIt file is;
* title
*
repeat-lines-of(free-energy-step-specs)
END
The end is optional and terminates the free energy run.
Lack of an END (i.e. and end-of-file or blank lines) will put PERT into
auto mode which will continue until LSTOP becomes 1.0 or 0.0 (based on
the LINCrement value).
The WHAM option allows to write a formatted file to post-process with
the Weighted Histogram Analysis Method (see below).
Description of PERT Commands
The PERTurb command copies and saves the current PSF and restraint
data for harmonic, NOE and dihedral restraints to the initial (lambda=0)
saved state. The SKIP command flags are also saved to allow linear scaling
of entire energy terms.
The structure may then be modified or perturbed with patches,
SCALar commands, with the DELEte command, or by generating or reading
a new PSF.
The Basic mode of operation is;
....
PERTurb ! Define the lambda=0 state.
PATCH .... ! Define the lambda=1 state.
DYNAmics .... ! Run MD on intermediate surface...
....
The PSF in use when dynamics or energy minimization is invoked
becomes the final (lambda=1) state. The actual energy computed
is a linear combination of these two endpoints.
The PATCH command may be replaced with any other command that
modifies the PSF. Some examples which modify the PSF;
SCALAR CHARGE SET -0.55 SELE ATOM A 1 O* SHOW END ! change a charge
DELETE ATOM SELE ALL END
READ SEQUENCE ....
GENERATE ... ! generate a new different PSF.
DELETE CONNECTIVITY .... ! modify the PSF by changing the connectivity.
SCALAR TYPE SET 14 SELE ATOM A 1 O SHOW END ! change the vdw atom type
OPEN .... ! Read a new PSF
READ PSF ...
It is not required that the PSF be modified. If one wants to
carry out coordinate perturbation only, it is sufficient to modify
the harmonic restraints, the NOE restraints, or the dihedral restraints.
In this way, it is possible to calculate the free energy differences
between different conformers. (However, this option should not be
used with simulataneous change of SHAKE constraints)
Note that with this implementation, because two PSFs are used,
that the connectivity may change. The use of 1-4 interactions and
nonbond exclusions is fully supported. This allows this method to be
used for examining changes that involve bond changes, such as cystine
bridge formation.
The advanced mode of operation is;
....
PERTurb
PATCH ....
DYNAmics ....
....
PERTurb
PATCH ....
DYNAmics ....
....
PERTurb
PATCH ....
DYNAmics ....
....
In this way, several changes can be affected in a single CHARMM run.
For example, the first patch might be the removal of charges, and the
second patch could correspond to a change in atom size, and the third
patch could simply consist of modifying dihedral restraints so as to
affect a conformational change. The free energy differences and
fluctuations will be calculated for each window as well as the total
for all previous windows.
RESTRICTIONS:
The number of atoms in both sets must match! If the system of
interest has different numbers of atoms, then dummy atoms must be
used. The mapping of atoms between the first and last structure is
one to one.
The following CHARMM features are not currently supported for use
with free energy perturbation;
INTEraction_energy
These commands will continue to work, but will only use the final
(lambda=1) structure. Most other energy related CHARMM features
are supported. The IMAGE facility will be supported in the near
future.
The following CHARMM energy related features cannot be modified with
the PERT command (e.g. cannot be part of what is changing, and are
only determined by the final state).
HBON - hydrogen bond energy
ST2 - ST2 energy
CIC - internal coordinate constraint energy
CDRO - quartic droplet potential energy
USER - user supplied energy (USERLINK)
RXNF - Reaction field energy
IMNB - image van der Waal energy
IMEL - image electrostatic energy
IMHB - image hydrogen bond energy
IMST - image ST2 energy
SBOU - solvent boundary energy
UREY - Urey Bradley energy terms
XTLV - Crystal vdw terms
XTLE - Crystal electrostatics
EXTE - Extended electrostatics
Some References: M Mezei and D.L. Beveridge, in Annals of the NYAS, "Free Energy Simulations" 482 (1986) T. P. Straatsma Ph.D. Thesis "Free Energy Evaluation by Molecular Dynamics Simulations" Kollman, P. A.; et al. J. Am. Chem. Soc. 1987, 109, 1607. Kollman, P. A.; et al. J. Am. Chem. Soc. 1987, 109, 6283. Kollman, P. A.; et al. J. Chem. Phys. 1989, 91, 7831. Bell, C. D.; Harvey, S. C., J. Phys. Chem. 1986, 90, 6595. van Gunsteren, W.F. et al. in: Computer Simulation of Biomolecular Systems: Theoretical and Experimental Applications, vol. 2, eds. van Gunsteren W.F. and Weiner P.K. (Escom, Leiden, 1994), p. 349
Examples:
The input file:
* A SIMPLE TEST RUN FOR PERT
*
bomlev -1
OPEN READ FILE UNIT 1 NAME ~/c22pt/toph19.mod
READ RTF UNIT 1
OPEN READ FILE UNIT 2 NAME ~/c22pt/param19.mod
READ PARAMETER UNIT 2
READ SEQUENCE CARD
* FIRST SEQUENCE FOR SECOND DERIVATIVE TEST
*
2
AMN CBX
GENERATE A
GENERATE B DUPLICATE A
OPEN UNIT 3 READ CARD NAME perttest.crd
READ COOR CARD UNIT 3
! modify the charge for the lambda=0 state
SCALAR CHARGE SET -0.55 SELE ATOM A 1 O* SHOW END
! minimize initial state so initial forces will be small.
MINI ABNR NSTEP 100 CTOFNB 12.0 CUTNB 14.0
PERT ! save all PSF data for the lambda=0 state
! modify the charge again for the lambda=1 state
SCALAR CHARGE SET -0.15 SELE ATOM A 1 O* SHOW END
! carry out free energy run from first to final state
OPEN READ CARD UNIT 88 NAME perttest.punit
DYNA VERLET STRT NSTEP 12000 TIMESTEP 0.001 -
IPRFRQ 100 IHTFRQ 0 IEQFRQ 100 NTRFRQ 2000 -
IUNCRD 50 ISEED 314159 -
NPRINT 100 NSAVC 0 NSAVV 0 INBFRQ 25 IHBFRQ 0 -
CTOFNB 12.0 CUTNB 14.0 -
FIRSTT 300.0 FINALT 300.0 TEMINC 100.0 -
IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 1 TWINDH 20.0 TWINDL -20.0 -
PUNIT 88
PERT OFF
energy ! just a check at lamda=1
STOP
The punit file:
* PUNIT FILE FOR SIMPLE TEST CASE
* use window method with 2000 steps of equilibration
* and 8000 steps of analysis for each of 10 evenly spaces
* windows.
*
LSTART 0.0 LAMBDA 0.0 LSTOP 0.05 PSTART 12000 PSTOP 20000 PWIND
LSTART 0.05 LAMBDA 0.1 LSTOP 0.15 PSTART 22000 PSTOP 30000 PWIND
LSTART 0.15 LAMBDA 0.2 LSTOP 0.25 PSTART 32000 PSTOP 40000 PWIND
LSTART 0.25 LAMBDA 0.3 LSTOP 0.35 PSTART 42000 PSTOP 50000 PWIND
LSTART 0.35 LAMBDA 0.4 LSTOP 0.45 PSTART 52000 PSTOP 60000 PWIND
LSTART 0.45 LAMBDA 0.5 LSTOP 0.55 PSTART 62000 PSTOP 70000 PWIND
LSTART 0.55 LAMBDA 0.6 LSTOP 0.65 PSTART 72000 PSTOP 80000 PWIND
LSTART 0.65 LAMBDA 0.7 LSTOP 0.75 PSTART 82000 PSTOP 90000 PWIND
LSTART 0.75 LAMBDA 0.8 LSTOP 0.85 PSTART 92000 PSTOP 100000 PWIND
LSTART 0.85 LAMBDA 0.9 LSTOP 0.95 PSTART 102000 PSTOP 110000 PWIND
LSTART 0.95 LAMBDA 1.0 LSTOP 1.0 PSTART 112000 PSTOP 120000 PWIND
END
Or equivalently using auto mode:
* PUNIT FILE FOR SIMPLE TEST CASE
* use window method with 2000 steps of equilibration
* and 8000 steps of analysis for each of 10 evenly spaces
* windows.
*
LSTART 0.0 LAMBDA 0.0 LSTOP 0.05 PSTART 12000 PSTOP 20000 PWIND PEQUIL 2000 PINCR 10000 LINCR 0.1
Or also equivalently as:
* PUNIT FILE FOR SIMPLE TEST CASE
* use window method with 2000 steps of equilibration
* and 8000 steps of analysis for each of 10 evenly spaces
* windows.
*
LSTART 0.0 LAMBDA 0.0 LSTOP 0.05 PSTART 12000 PSTOP 20000 PWIND
LSTART 0.05 LAMBDA 0.1 LSTOP 0.15 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.15 LAMBDA 0.2 LSTOP 0.25 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.25 LAMBDA 0.3 LSTOP 0.35 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.35 LAMBDA 0.4 LSTOP 0.45 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.45 LAMBDA 0.5 LSTOP 0.55 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.55 LAMBDA 0.6 LSTOP 0.65 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.65 LAMBDA 0.7 LSTOP 0.75 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.75 LAMBDA 0.8 LSTOP 0.85 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.85 LAMBDA 0.9 LSTOP 0.95 PEQUIL 2000 PINCR 10000 PWIND
LSTART 0.95 LAMBDA 1.0 LSTOP 1.0 PEQUIL 2000 PINCR 10000 PWIND
END
------------------------------------------------------------------------------
An annotated output example:
This output is a short excerpt from perttest.out where the input
lines start with "|" and output lines start with ".....|".
The CHARMM command is:
| CHARMM> DYNA VERLET STRT NSTEP 12000 TIMESTEP 0.001 -
| CHARMM> IPRFRQ 100 IHTFRQ 0 IEQFRQ 100 NTRFRQ 2000 -
| CHARMM> IUNCRD 50 ISEED 314159 -
| CHARMM> NPRINT 100 NSAVC 0 NSAVV 0 INBFRQ 25 IHBFRQ 0 -
| CHARMM> CTOFNB 12.0 CUTNB 14.0 -
| CHARMM> FIRSTT 300.0 FINALT 300.0 TEMINC 100.0 -
| CHARMM> IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 1 TWINDH 20.0 TWINDL -20.0 -
| CHARMM> PUNIT 88
and the relevent punit data is in perttest.punit:
|* PUNIT FILE FOR SIMPLE TEST CASE
|* use window method with 2000 steps of equilibration
|* and 8000 steps of analysis for each of 10 evenly spaces
|* windows.
|*
| LSTART 0.0 LAMBDA 0.0 LSTOP 0.05 PSTART 1200 -
| PSTOP 2000 PWIND PEQUIL 200 PINCR 1000 LINCR 0.1
The output starting at line 1618 (in the middle of the dynamics command) is:
.....| PERTURBATION> Free energy perturbation results:
This indicates that a "window" was just completed.
**
.....| PERTURBATION> results, LSTART= 0.050000 LSTOP= 0.150000 LLAST= 0.100000 Number of steps used= 800
This says that the window started at lambda=0.05 and ended at 0.15.
The window "center" was at lambda=0.1
A total of 800 steps was used for collecting averages and fluctuations.
.....| PERTURBATION> result: EXPAVE=0.155456E+01 EXPFLC=0.225213E+00 DIFAVE= -0.256530 DIFFLC= 0.088960
The values:
EXPAVE is the time average of exp((ef(t)-ei(t)-ef(0)+ei(0))/kT)
EXPFLC is the fluctuation of this value about its average
DIFAVE is the time average of (ef(t)-ei(t)-ef(0)+ei(0))
DIFFLC is the fluctuation of this value about its average
Note: this value should not be much larger than kT for a
good window schedule. If this value is too large,
then smaller window lambda steps should be used.
The value here (0.09) indicates that a much larger
window would have been OK. ef(t) is the energy at
lambda=LSTOP, ei(t) is the energy at lambda=LSTART
ef(0) is the initial energy at LSTOP, ei(0) is the
initial energy at LSTART
.....| PERTURBATION> TP Windowing result, EPRTOT= 1.392400 EFORWARD= 0.914282 EPREF= 1.177303
This is the old format. In the new format the values EPRTOT,EFORWARD,
EBACKWARD,EPREFF,EPREFB are reported where the forward energy is from LLAST
to LSTOP and the backward from LLAST to LSTART. Separating the window into
two halves (double wide sampling) improves the accuracy of the the TP method.
.....| PERTURBATION> TI Windowing result, EPRTOT= 1.400218 EFORWARD= 0.920773 EPREF= 1.177303
EPRTOT is the total energy for this window and all previous (since a PERT
RESET).
EFORWARD is the energy for this current window.
EPREF is the initial energy difference (ef(0)-ei(0))
.....| PERTRES> LSTART= 0.05000 LSTOP= 0.15000 EPRTOT= 1.40022 EFORWARD= 0.92077 EPREF= 1.17730 DIFAVE= -0.25653 DIFFLC= 0.08896
This is the same data on a one line format. To use this, grep (search)
for "PERTRES".
.....| PERTURBATION> Averages for the last 800 steps:
.....|PAVE DYN: Step Time TOTEner TOTKe ENERgy TEMPerature
.....|PAVE PROP: GRMS VEREner VERKe EHFCor VIRKe
.....|PAVE EXTERN: VDWaals ELEC HBONds USER
.....|PAVE PRESS: VIRE VIRI PRESSE PRESSI VOLUme
.....| ---------- --------- --------- --------- --------- ---------
.....|PAVE> 800 0.00000 0.92077 -0.00039 0.92077 0.00000
<delF*v> <delE>
.....|PAVE PROP> 0.00000 0.00000 0.00000 0.00000 0.00000
.....|PAVE EXTERN> 0.00000 0.92077 0.00000 0.00000
.....|PAVE PRESS> 0.00000 0.00000 0.00000 0.00000 0.00000
This is the average values for this window (TI results).
Note: the <delF*v> is a "correction" term that shows how well the window is
equilibrated. This value should be close to zero and much smaller than
its fluctuation. If it is not, then the assumptions required for a
free energy calculation using windowing are not met.
This can occur if there is a "snapping" event with releases energy in
an irreversible manner, or if the system is not at equilibrium.
For a slow growth "window", this is a correction term which should be
multiplied by the estimated delay of the equilibration at the current
step and then added to the total.
For example, if the configuration distribution lags behind the energy
potential by 100 steps, this value should be scaled by 100.
When forcing a "continuous" change through slow growth, there tends to
be a delay since the structure does not respond instantly to the
potential.
.....| ---------- --------- --------- --------- --------- ---------
.....| PERTURBATION> Fluctuations for the last 800 steps:
.....|PFLC> 800 0.00000 0.08911 0.00985 0.08896 0.00000
.....|PFLC PROP> 0.00000 0.00000 0.00000 0.00000 0.00000
.....|PFLC EXTERN> 0.00000 0.08896 0.00000 0.00000
.....|PFLC PRESS> 0.00000 0.00000 0.00000 0.00000 0.00000
This is the fluctuation data for the current window (TI results).
.....| ---------- --------- --------- --------- --------- ---------
.....| PERTURBATION> TOTALS since last reset:
.....|PTOT> 800 0.00000 1.40022 -0.00038 1.40022 0.00000
.....|PTOT PROP> 0.00000 0.00000 0.00000 0.00000 0.00000
.....|PTOT EXTERN> 0.00000 1.40022 0.00000 0.00000
.....|PTOT PRESS> 0.00000 0.00000 0.00000 0.00000 0.00000
This is the total for this window and all previous (since the last PERT RESET).
.....| ---------- --------- --------- --------- --------- ---------
.....|
.....| PERTURBATION> EOF on punit file: PERT in auto mode.
This indicates that no more data was found on the PUNIT file.
.....| PERTURBATION> Free energy perturbation calculation continues.
Now we start a new window.
.....| PERTURBATION> PSTART= 3200 PSTOP= 4000
We will run equilibration until step 3200, then we will collect data until step 4000.
.....| PERTURBATION> LSTART= 0.150000 LSTOP= 0.250000 LAMBDA= 0.200000
This indicates the boundaries of the new window.
.....| PERTURBATION> Windowing will be used.
We are using windowing (fixed lambda values), instead of slow growth.
.....|
If SHAKE is applied to bond terms which are changed as the
result of an alchemical mutation a constraint correction is calculated
where required in slow-growth mode and TI in windowing mode. The
exponential formula in windowing mode is not supported. The user has
to beware of subtle problems regarding a possible "moment of inertia"
term that may be or may be not included in this correction (S. Boresch
& M. Karplus, to be published) In order for the constraint correction
to work properly, attention has to be given to the following points:
(1) SHAKE must not be applied to angle terms
(2) the PARA option has to be used (it's not clear, how to support
reference coordinates in the context of an alchemical mutation)
(3) the SHAKE command has to issued after the PERT command. (This
is similar to setting up IMAGEs in connection with PERT). A
typical input will look something like
PERT SELE ... END
!change psf; after ALL changes have been made
SHAKe BOND PARA
DYNA ... ! carry out MD simulations etc.
PERT OFF
(4) One should not mix situations where a constraint correction for
SHAKE is necessary with the use of harmonic, NOE and dihedral
restraints to calculate conformational free energy differences.
Items (2) and (3) can lead to error messages in situations where there
is actually no problem, e.g. you just want to apply SHAKe to your
solvent which is not affected by the mutation, so you specify SHAKe
before PERT and "bomb". Nevertheless, I thought better safe than
sorry and if wanted one can override the warnings with a BOMLEV -3.
Item (4) is simply untested.
The post-processing Wheigthed Histogram Analysis Method (WHAM) can
be used to help reaching better statistical convergence on free energy
perturbations calculations. The approach represents a generalization
of the histogram method developed by Ferrenberg and Swendsen (1989).
The central idea, which goes back to the maximum overlap method
developed by Bennet (1976) to estimate free energy differences,
consists in constructing an optimal estimate of the unbiased distribution
function as a weighted sum over the data extracted from all the
simulations and determining the functional form of the weight factors that
minimizes the statistical error. The WHAM approach can be used to calculate
the PMF along coordinates (Kollman, 1992; Boczko, 1993; Brooks, 1993;
Roux, 1995) and can also be used to post-process free energy perturbation
calculations in which no PMF is desired.
Assuming that the simulations are performed (as in PERT) with a potential
function given by a linear switching from E_0 to E_1, that is,
E_k = (Lambda_k-1)*E_1 + (Lambda_k)*E_0
The free energy constants F_k corresponding to any Lambda_k window are
given by,
exp(-F_k/kBT) = Sum_i { Sum_t ( Top_i(t) / Bot_i(t) ) },
where
Top_i(t) = Exp(-E_k[X_i(t)]/kBT)
and
Bot_i(t) = Sum_j { Ntime(j) * Exp(+F_j/kBT-E_j[x_i(t)]/kBT) }
where E_j[x_i(t)] = (lambda_j-1)*E_1[x_i(t)] + (lambda_j)*E_0[x_i(t)]
is the potential energy function of the j-th window evaluated with
the configuration taken from the i-th simulation. The WHAM equations
for the F_k must be solved interatively.
The syntax of the command is:
WHAM MAXWindow <integer> MAXTime <integer> unit <integer> -
tolerance <real> nstep <integer> [guess] -
ioffset <integer> nskip <integer> {lambda <real> lambda <real> ...}
where
MAXWindow is the total number of windows
MAXTime is the total number of time-step configuration for each window
unit is the unit of the formatted file containing all the information
tolerance is the tolerance on the F_k to reach convergence
nstep is the maximum number of iterations on the WHAM equations
guess to flag that an initial guess is provided for the F_k. Those
are read directly from the input stream with one line per
window in the format [ window <integer> F() <real> ]
nskip use only every nskip data point to reach convergence (faster)
lambda give any value of lambda for which you want the free energy (a list)
ioffset reference energy level at window number "ioffset".
The file containing the information can be written by PERT during dynamics
if the keyword WHAM is used in the dynamics command (see above). In
principle, the WHAM could also be used with non-linear perturbations,
but then the code in PERT would have to evaluate several energies since
those could not be obtained by lambda interpolations.
Some references on WHAM:
S. Kumar, D. Bouzida, R.H. Swendsen, P.A. Kollman, and J.M. Rosenberg.
J. Comp. Chem. 13, 1011--1021 (1992).
E.M. Boczko and C.L. Brooks III. J. Phys. Chem. 97, 4509--4513 (1993).
A.M. Ferrenberg and R.H. Swendsen. Phys. Rev. Lett. 63, 1195--1198 (1989).
C.M. Bennet. J. Comp. Phys. 22, 245--268 (1976).
C.L. Brooks III and L. Nilsson. J. Am. Chem. Soc. 115, 11034--11035 (1993).
B. Roux, Comp. Phys. Comm. 91, 275-282 (1995).
CHARMM Documentation / Rick_Venable@nih.gov