We find the sander module to be a flexible way of incorporating a variety of restraints into a optimization procedure that includes energy minimization and dynamical simulated annealing. The "standard" sorts of NMR restraints, derived from NOE and J-coupling data, can be entered in a way very similar to that of programs like DISGEO, DIANA or X-PLOR; an aliasing syntax allows for definitions of pseudo-atoms, connections with peak numbers in spectra, and the use of "ambiguous" constraints from incompletely-assigned spectra. More "advanced" features include the use of time-averaged constraints, use of multiple copies (LES) in conjunction with NMR refinement, and direct refinement against NOESY intensities, paramagnetic and diamagnetic chemical shifts, or residual dipolar couplings. In addition, a key strength of the program is its ability to carry out the refinements (usually near the final stages) using an explicit-solvent representation that incorporates force fields and simulation protocols that are known to give pretty accurate results in many cases for unconstrained simulations; this ability should improve predictions in regions of low constraint density and should help reduce the number of places where the force field and the NMR constraints are in "competition" with one another.
Since there is no generally-accepted "recipe" for obtaining solution structures from NMR data, the comments below are intended to provide a guide to some commonly-used procedures. Generally speaking, the programs that need to be run to obtain NMR structures can be divided into three parts:
linkthick=1.5# second row:
PEAK: circle radius 0.4 "NOESY" "peak-list" right arrow 0.75 from PEAK.e "various" above "calibrations" below
DIST: box "7/8 column" "distance" "bounds" arrow 0.75 from DIST.e "makeDIST_RST" above
RDIS: ellipse "distance" "restraints"
arrow 0.75 from RDIS.e "makeVOL_RST" above
NOEEXP: ellipse "volume" "restraints"
# top row:
move up from PEAK.n RCS: ellipse "chemical" "shift" "restraints"
move up from DIST.n CS: circle radius 0.4 "chemical" "shifts" right arrow 0.75 from CS.e "analyze" above "ambiguities" below
MAP: box "MAP" "file"
move right 0.75 from MAP.e DMAP: ellipse "default" "MAP file"
arrow from CS.w to RCS.e "SHIFTS or" above "FANTASIAN" below
# third row:
move down from PEAK.s SPECTRUM: circle radius 0.4 "spectrum"
move right 0.75 from SPECTRUM.e COUPLING: box "5 column" "coupling" "constants" right arrow 0.75 from COUPLING.e "makeANG_RST" above
RANG: ellipse "J-coupling" "restraints"
move right 0.75 from RANG.e DISANG: ellipse "DISANG" "file"
# fourth row:
move down from COUPLING.s CHIR: ellipse "default" "chirality" "inform." right arrow 0.75 from CHIR.e "makeCHIR_RST" above
RCHIR: ellipse "chirality" "restraints"
arrow from RDIS.se to DISANG.nw arrow from RANG.e to DISANG.w arrow from RCHIR.ne to DISANG.sw arrow from MAP.s to RDIS.n arrow from DMAP.w to MAP.e arrow from SPECTRUM.n to PEAK.s arrow from SPECTRUM.e to COUPLING.w arrow from PEAK.ne to CS.sw spline from SPECTRUM.s down 1.6 then right 1.1 ->
# fifth row:
move down from CHIR.s DIP: box "direct" "dipolar" "couplings" right arrow 0.75 from DIP.e "makeDIP_RST" above
RDIP: ellipse "alignment" "restraints"
Notation: circles represent logical information, whose format might differ from one project to the next; solid rectangles are in a specific format (largely compatible with DIANA and other programs), and are intended to be read and edited by the user; ellipses are specific to Sander, and are generally not intended to be read or edited manually. The conversion of NOESY volumes to distance bounds can be carried out by a variety of programs such as mardigras or xpk2bound that are not included with Amber. Similarly, the analysis and partial assignment of ambiguous or overlapped peaks is a separate task; at TSRI, these are typically carried out using the programs xpkasgn and filter.pl.
Once the restraint files are made, here is a flow-chart of the general way in which sander refinements are carried out:
linkthick=1.5PDB: circle radius 0.4 "generic" "PDB-file" right arrow 0.75 from PDB.e "protonate" above
APDB: box "Amber" "PDB-file" right arrow 0.75 from APDB.e "LEaP" above
PRMTOP: box "prmtop &" "prmcrd" "files"
# second row:
move down from APDB.s MDIN: box "sander" "control" "file" right arrow 0.75 from MDIN.e
SANDER: ellipse "SANDER" move right 0.75 from SANDER.e RST: box "NMR" "restraints"
# third row:
move down from MDIN.s OUTPUTS: ellipse "sander" "output" left arrow 0.75 from OUTPUTS.w "sviol" above
VIOL: circle radius 0.4 "violation" "statistics"
move down from SANDER.s COORDS: ellipse "Amber" "coordinates" right arrow 0.75 from COORDS.e "ambpdb" above
PDBOUT: box "output" "pdb-files"
arrow from PRMTOP.s to SANDER.n arrow from SANDER.sw to OUTPUTS.n arrow from SANDER.s to COORDS.n arrow from RST.w to SANDER.e
Once families of PDB-format structures are available, there is a wide variety of tools available elsewhere to perform comparisons and analyses.
The basic ideas of this scheme owe a lot to the general experience of the nmr community over the past decade. Some good papers to look at are:
These papers outline procedures in the Scripps group, from which a lot of the NMR parts of SANDER are derived. They are by no means the only way to proceed. We hope that the flexibility incorporated into SANDER will encourage folks to experiment with refinement protocols.
The makeDIST_RST program converts a simplified description of distance bounds into a detailed input for sander. A variety of input and output filenames may be specified on the command line:
input:
-noe 8-col file of upper and lower bounds, OR
-ual 7-col file of NOESY volumes
-vol Brookhaven format file
-pdb MAP file (default:map.DG-AMBER)
-map LES atom mappings, made by addles
-les DGEOM95 restraint format
output:
-dgm SANDER restraint format
-rst Sander Volume Format
-svf (gives you this explanation, overrides other parameters)
other options:
-help (gives you short runtime diagnostic output)
-report (do not correct upper bound for r**-6 averaging)
-nocorr (use alternative form for the distance restraints)
-altdis
The 7/8 column distance bound file is essentially that used by the DIANA or DISGEO programs. It consists of one-line per restraint, which would typically look like the following:
If all peaks involved just single protons, and were fully assigned, this is all that one would need. In general, though, some peaks (especially methyl groups or fast-rotating aromatic rings) represent contributions from more than one proton, and many other peaks may not be fully assigned. Sander handles both of these situations in the same way, through the notion of an "ambiguous" peak, that may correspond to several assignments. These peaks are given two types of special names in the 7/8-column format file:
The default and molecule-specific MAP files are combined into a single file, which is used, along with the 7-column restraint file, the the program makeDIST_RST to construct the actual sander input files. You should consult the help file for makeDIST_RST for more information. For example, here are some lines added to the MAP file for a recent TSRI refinement:
Here the spectrum name and peak number were used to construct a label for each ambiguous peak. Then, an entry in the restraint file might look like this:
If the -les flag is set, the program will prepare distance restraints for multiple copies (LES) simulations. In this case, the input pdb file is one without LES copies, i.e. with just a single copy of the molecule. The "lesfile" specified by this flag is created by the addles program, and contains a mapping from original atom numbers into the copy numbers used in the multiple-copies simulation.
Distance constraints defined here cannot span two different copy sets, i.e., there cannot be one atom of a constraint torsion that is in one multiple copy set, where the other atom from the same restraint is in another copy set. It is OK to have some atoms with single copies, and others with multiple copies in the same restraint.
There are fewer "standards" for representing coupling constant information. We have followed the DIANA convention in the program makeANG_RST. This program takes as input a five-column torsion angle constraint file along with an AMBER pdb file of the molecule. It creates as output (to standard out) a list of constraints in RST format that is readable by AMBER.
The input torsion angle constraint file can be read from standard in or from a file specified by the -con option on the command line. The input constraint file should look something like this:
Note: It is assumed that the lower bound and the upper bound define a
region of allowed conformation on the unit circle that is swept out in a
clockwise direction from
. If the number in the lb column
is greater than the
the number in the ub column, 360o will successively be
subtracted from the lb
until
. This preserves the clockwise definition of the allowed
conformation space, while also making the number that specifies the lower
bound less than the number that specifies the upper bound, as is required by
AMBER. If this occurs, a warning message will be printed to stderr to
notify the user that the data has been modified.
The angles that one can constrain in this manner are defined in the library file that can be optionally specified on the command line with the -lib flag, or the default library "tordef.lib" (written by Garry P. Gippert) will be used. If you wish to specify your own nomenclature, or add angles that are not already defined in the default file, you should make a copy of this file and modify it to suit your needs. The general format for an entry in the library is:
If the first letter of the second field is "J", this torsion is assumed to be a J-coupling constraint. In that case, three additional floats are read at the end of the line, giving the A,B and C coefficients for the Karplus relation for this torsion. For example:
This program also supports pseudorotation phase angle constraints for prolines and nucleic acid sugars; each of these will generate restraints for the 5 component angles which correspond to the lb and ub values of the input psuedorotation constraint. In the torsion library, a pseudorotation definition looks like:
PPA stands for Pseudorotation Phase Angle and is the angle that should appear in the input constraint file when using pseudorotation constraints. The program then uses the definition of that PPA angle in the library file to look for the 5 other angles (NU0-NU4 in this case) which it then generates restraints for. PPA for proline residues is included in the standard library as well as for the DNA nucleotides.
If the -les flag is set, the program will prepare torsion angle restraints for multiple copies (LES) simulations. In this case, the input pdb file is one without LES copies, i.e. with just a single copy of the molecule. The "lesfile" specified by this flag is created by the addles program, and contains a mapping from original atom numbers into the copy numbers used in the multiple-copies simulation.
Torsion angle constraints defined here cannot span two different copy sets, i.e., there cannot be some atoms of a particular torsion that are in one multiple copy set, and other atoms from the same torsion that are in other copy sets. It is OK to have some atoms with single copies, and others with multiple copies in the same torsion. The program will create as many duplicate torsions as there are copies.
A good alternative to interpreting J-coupling constants in terms of torsion
angle restraints is to refine directly against the coupling constants
themselves, using a appropriate Karplus relation. See the discussion of the
variable RJCOEF, above.
We also find it useful to add chirality constraints and trans-peptide
constraints (where appropriate) to prevent chirality inversions or
peptide bond flips during the high-temperature portions of simulated
annealing runs. The program makeCHIR_RST will create these
constraints. Note that you may have to edit the output of this program to
change trans peptide constraints to cis, as appropriate.
Refinement directly against measured NOESY volume restraints is also possible in sander. The makeVOL_RST is used to prepare detailed input for this.
Defaults: cutoff is 5. Ang.,
On output, stdout contains a file that, with some minimal hand editing, should serve as input to SANDER. More explicitly, it creates groups and &noeexp namelists that inlcude many of the variables you will need. Any missing variables should be added by hand. [I do this since it seems easier to hand-edit the output file than to arrange for a lot of pass-throughs from this program into the output file.]
During refinement, measured NOESY volumes
are included in penalty functions that depend upon
where
is the experimentally measured value, and I is the value
corresponding the current conformation; the functional form of the
penalty depends upon the ipnlty variable.
Careful experimentation
will undoubtedly be required for each data set to define a reasonable
penalty function. Simply weighting each observed peak equally (with
the default values of awt and arange) is almost certainly
a bad idea, since this effectively gives too much influence to the
strong peaks at the expense of longer-range information.
For simulations with residual dipolar coupling restraints, makeDIP_RST is a simple code to prepare the input file. Type makeDIP_RST -help to obtain a more detailed description of the usage. For now, this code only handles backbone NH and CaH data. The header specifying values for various parameters needs to be manually added to the output of makeDIP_RST.
Use of residual dipolar coupling restraints is new both for AMBER and for the general NMR community. Refinement against these data should be carried out with care, and the optimal values for the force constant, penalty function, and initial guesses for the alignment tensor components are still under investigation. Here are some suggestions from the experiences so far:
If you specify LISTOUT=POUT when running sander, the output file
will contain a lot of detailed information about the remaining restraint
violations at the end of the run. When running a family of structures, it
can be useful to process these ouput files with sviol, which takes
a list of sander output files on the command line, and sends a summary
of energies and violations to STDOUT. If you have more than 20 or so
structures to analyze, the output from sviol becomes unwieldy. In
this case you may also wish to use sviol2, which prints out somewhat
less detailed information, but which can be used on larger families of
structures. The senergy script gives a more detailed view of
force-field energies from a series of structures. (We thank the TSRI NMR
community for helping to put these scripts together, and for providing many
useful suggestions.)
In carrying out NMR refinement simulations, it is common to modify the standard force field. The ifstrp variable allows one to treat nonbonded interactions as soft repulsions with no electrostatic contributions. Since bonds and angles are kept close to ideal values by the force constants inherent in the standard force field, and since the intrinsic dihedral barriers for single bonds are also quite small, this provides a "generic" or simplified representation of the allowed conformational space that may appeal to some users.
Other users will wish to use force fields that incorporate more of what we know about relative conformational energies, i.e. a more elaborate force field. Even here, though, some modifications may be advisable. For example, modifications we have found useful for peptides/proteins include increasing the torsional force constant for the peptide bond (to reduce the tendency of restained simulations to produce badly distorted bonds) and a reduction in the net charge of charged side chains (to compensate in part for neglect of explicit solvent). These changes are incorporated into a database and a frcmod file in the src/nmr_aux/forcefield subdirectory; see the README file in that directory.
The model of the previous sections involves the "simgle-average-structure" idea, and tries to fit all constraints to a single model, with minimal deviations. A generalization of this model treats distance constraints arising from from NOE crosspeaks (for example) as being the average distance determined from a trajectory, rather than as the single distance derived from an average structure. Time-averaged bonds and angles are calculated as

where



Time-averaged torsions are calculated as

where
is the torsion, and
and
are calculated using the equation above with
or
substituted for r(t').
Forces for time-averaged restraints can be calculated either of two ways. This option is chosen with the DISAVI / ANGAVI / TORAVI commands (Section 1). In the first (the default),

(and analogously for y and z).
The forces then correspond to the standard flat-bottomed well functional
form, with the instantaneous value of the internal replaced by the
time-averaged value. For example, when
,

and similarly for other ranges of
.
When the second option for calculating forces is chosen (IINC = 1 on a DISAVI, ANGAVI or TORAVI card), forces are calculated as

For example, when
,

Integration of this equation does not give Equation (4), but rather
a non-intuitive expression for the energy (although one that still
forces the bond to the target range). The reason that it may
sometimes be preferable to use this second option is that the term
,
which occurs in the exact expression [Eq. (3)],
varies as
. When i=3, this means the
forces can be varying with the fourth power the distance, which
can possibly lead to very large transient forces and instabilities
in the molecular dynamics trajectory. [Note that this will not be the case
when linear scaling is performed, i.e. when i=-1, as is generally
the case for valence and torsion angles. Thus, for linear scaling,
the default (exact) force calculation should be used].
It should be noted that forces calculated using Equation (5) are not conservative forces, and would cause the system to gradually heat up, if no velocity rescaling were performed. The temperature coupling algorithm should act to maintain the average temperature near the target value. At any rate, this heating tendency should not be a problem in simulations, such as fitting NMR data, where MD is being used to sample conformational space rather than to extract thermodynamic data.
This section has described the methods of time-averaged restraints. For more discussion, the interested user is strongly urged to consult recent studies:
NMR restraints can be made compatible with the multiple copies (LES) facility; see the following chapter for more information about LES. To use NMR constraints with LES, you need to do two things:
file wnmr name=(lesnmr) wovr" to you input to
addles. The filename (lesnmr in this example) may be whatever
you wish. This will cause addles to output an additional file that is
needed at the next step.
-les lesnmr" to the command line arguments to makeDIST_RST.
This will read in the file created by addles containing information
about the copies. All NMR restraints will then be interpreted as
"ambiguous" restraints, so that if any of the copies satisifies the
restraint, the penalty goes to zero.
Note that although this scheme has worked well on small peptide test cases, we have yet not used it extensively for larger problems. This should be treated as an experimental option, and users should use caution in applying or interpreting the results.
The next few pages contain excerpts from some sample NMR refinement files used at TSRI. The first example just sets up a simple (but often effective) simulated annealing run. You may have to adjust the length, temperature maximum, etc. somewhat to fit your problem, but these values work well for many "ordinary" NMR problems.
+---------------------------------------------------------------------------+ | 1. Simulated annealing NMR refinement | +---------------------------------------------------------------------------+ || | 15ps simulated annealing protocol | | &cntrl | | nstlim=15000, ntt=1, (time limit, temp. control) | | scee=1.2, (scee must be set - 1-4 scale factor) | | ntpr=500, pencut=0.1, (control of printout) | | ipnlty=1, nmropt=1, (NMR penalty function options) | | vlimit=10, (prevent bad temp. jumps) | | ntb=0, (non-periodic simulation) | | igb=3, (don't use pme, do use r dielectric) | | &end | |# | |# Simple simulated annealing algorithm: | |# | |# from steps 0 to 1000: raise target temperature 10->1200K | |# from steps 1000 to 3000: leave at 1200K | |# from steps 3000 to 15000: re-cool to low temperatures | |# | | &wt type='TEMP0', istep1=0,istep2=1000,value1=10., | | value2=1200., &end | | &wt type='TEMP0', istep1=1001, istep2=3000, value1=1200., | | value2=1200.0, &end | | &wt type='TEMP0', istep1=3001, istep2=15000, value1=0., | | value2=0.0, &end | |# | |# Strength of temperature coupling: | |# steps 0 to 3000: tight coupling for heating and equilibration | |# steps 3000 to 11000: slow cooling phase | |# steps 11000 to 13000: somewhat faster cooling | |# steps 13000 to 15000: fast cooling, like a minimization | |# | | &wt type='TAUTP', istep1=0,istep2=3000,value1=0.2, | | value2=0.2, &end | | &wt type='TAUTP', istep1=3001,istep2=11000,value1=4.0, | | value2=2.0, &end | | &wt type='TAUTP', istep1=11001,istep2=13000,value1=1.0, | | value2=1.0, &end | | &wt type='TAUTP', istep1=13001,istep2=14000,value1=0.5, | | value2=0.5, &end | | &wt type='TAUTP', istep1=14001,istep2=15000,value1=0.05, | | value2=0.05, &end | | | | (continued on next page) | +---------------------------------------------------------------------------+
+---------------------------------------------------------------+ | 1. Simulated annealing NMR refinement (continued) | +---------------------------------------------------------------+ |# | |# "Ramp up" the restraints over the first 3000 steps: | |# | | &wt type='REST', istep1=0,istep2=3000,value1=0.1, | | value2=1.0, &end | | &wt type='REST', istep1=3001,istep2=15000,value1=1.0, | | value2=1.0, &end | | | | &wt type='END' &end | |LISTOUT=POUT (get restraint violation list) | |DISANG=RST.f (file containing NMR restraints) | +---------------------------------------------------------------+
The next example just shows some parts of the actual RST file that
sander would read. This file would ordinarily not be made or edited
by hand; rather, run the programs makeDIST_RST, makeANG_RST and
makeCHIR_RST, combining the three outputs together to construct the
RST file.
+--------------------------------------------------------------------+ | 2. Part of the RST.f file referred to above | +--------------------------------------------------------------------+ || |# first, some distance constraints prepared by makeDIST_RST: | |# (comment line is input to makeRST, &rst namelist is output) | |# | |#( proton 1 proton 2 upper bound) | |#--------------------------------------------- | |# | |# 2 ILE HA 3 ALA HN 4.00 | |# | | &rst iat= 23, 40, r3= 4.00, r4= 4.50, | | r1 = 1.3, r2 = 1.8, rk2=0.0, rk3=32.0, ir6=1, &end | |# | |# 3 ALA HA 4 GLU HN 4.00 | |# | | &rst iat= 42, 50, r3= 4.00, r4= 4.50, &end | |# | |# 3 ALA HN 3 ALA MB 5.50 | |# | | &rst iat= 40, -1, r3= 6.22, r4= 6.72, | | igr1= 0, 0, 0, 0, igr2= 44, 45, 46, 0, &end | |# | |# .......etc...... | |# | |# next, some dihedral angle constraints, from makeANG_RST: | |# | | &rst iat= 213, 215, 217, 233, r1=-190.0, | | r2=-160.0, r3= -80.0, r4= -50.0, &end | | | | &rst iat= 233, 235, 237, 249, r1=-190.0, | | r2=-160.0, r3= -80.0, r4= -50.0, &end | | | |# .......etc....... | | | +--------------------------------------------------------------------+
+--------------------------------------------------------------------+ | 2. Part of the RST.f file referred to above (continued) | +--------------------------------------------------------------------+ |# | |# | |# next, chirality and omega constraints prepared by makeCHIR_RST: | |# | |# | |# chirality for residue 1 atoms: CA CG HB2 HB3 | | &rst iat= 3 , 8 , 6 , 7 , | | r1=10., r2=60., r3=80., r4=130., rk2 = 10., rk3=10., &end | |# | |# chirality for residue 1 atoms: CB SD HG2 HG3 | | &rst iat= 5 , 11 , 9 , 10 , &end | |# | |# chirality for residue 1 atoms: N C HA CB | | &rst iat= 1 , 18 , 4 , 5 , &end | |# | |# chirality for residue 2 atoms: CA CG2 CG1 HB | | &rst iat= 22 , 26 , 30 , 25 , &end# | |# | | ......etc........ | | | |# trans-omega constraint for residue 2 | | &rst iat= 22 , 20 , 18 , 3 , | | r1=155., r2=175., r3=185., r4=205., rk2 = 80., rk3=80., &end | |# | |# trans-omega constraint for residue 3 | | &rst iat= 41 , 39 , 37 , 22 , &end | |# | |# trans-omega constraint for residue 4 | | &rst iat= 51 , 49 , 47 , 41 , &end | |# | |# ......etc........ | |# | +--------------------------------------------------------------------+
The next example is an input file for volume-based NOE refinement. As with
the distance/angle RST file shown above, the user would generally not
construct this file, but create it from a "7-column" file using the
makeVOL_RST program. Hand-editing might be used at the top of the file, to
change the correlation times, etc.
+---------------------------------------------------------------------------------+ | 3. Sample NOESY intensity input file | +---------------------------------------------------------------------------------+ |# A part of a NOESY intensity file, made by makeVOL_RST : | | | | &noeexp | | id2o=1, (exchangeable protons removed) | | oscale=6.21e-4, (scale between exp. and calc. intensity units) | | taumet=0.04, (correlation time for methyl rotation, in ns.) | | taurot=4.2, (protein tumbling time, in ns.) | | NPEAK = 13*3, (three peaks, each with 13 mixing times) | | EMIX = 2.0E-02, 3.0E-02, 4.0E-02, 5.0E-02, 6.0E-02, | | 8.0E-02, 0.1, 0.126, 0.175, 0.2, 0.25, 0.3, 0.35, | | (mixing times, in sec.) | | IHP(1,1) = 13*423, IHP(1,2) = 13*1029, IHP(1,3) = 13*421, | | (number of the first proton) | | JHP(1,1) = 78*568, JHP(1,2) = 65*1057, JHP(1,3) = 13*421, | | (number of the second proton) | | AEXP(1,1) = 5.7244, 7.6276, 7.7677, 9.3519, | | 10.733, 15.348, 18.601, | | 21.314, 26.999, 30.579, | | 33.57, 37.23, 40.011, | | (intensities for the first cross-peak) | | AEXP(1,2) = 8.067, 11.095, 13.127, 18.316, | | 22.19, 26.514, 30.748, | | 39.438, 44.065, 47.336, | | 54.467, 56.06, 60.113, | | AEXP(1,3) = 7.708, 13.019, 15.943, 19.374, | | 25.322, 28.118, 35.118, | | 40.581, 49.054, 53.083, | | 56.297, 59.326, 62.174, | | &end | |SUBMOL1 | |RES 27 27 29 29 39 41 57 57 70 70 72 72 82 82 (residues in this submol) | |END | |END | | | +---------------------------------------------------------------------------------+
Next, we illustrate the form of the file that holds residual dipolar coupling restraints. Again, this would generally be created from a human-readable input using the program makeDIP_RST.
+-----------------------------------------------------------------+ | 5. Residual dipolar restraints, prepared by makeDIP_RST: | +-----------------------------------------------------------------+ | &align | | ndip=91, dcut=-1.0, fj=0.1, dfac_al=37*-0.46623, 54*-0.21341, | | s11=3.883, s22=53.922, s12=33.855, s13=-4.508, s23=-0.559, | | id(1)=188, jd(1)=189, dobs(1)= 6.24, | | id(2)=208, jd(2)=209, dobs(2)= -10.39, | | id(3)=243, jd(3)=244, dobs(3)= -8.12, | | | | .... | | id(91)=1393, jd(91)=1394, dobs(91)= -19.64, | | &end | +-----------------------------------------------------------------+
Finally, we show how the detailed input to sander could be used to
generate a more complicated restraint. Here is where the user would have to
understand the details of the RST file, since there are no "canned"
programs to create this sort of restraint. This illustrates, though, the
potential power of the program.
+---------------------------------------------------------------------------------------+
| 5. A more complicated constraint |
+---------------------------------------------------------------------------------------+
| |
|# 1) Define two centers of mass. COM1 is defined by |
|# {C1 in residue 1; C1 in residue 2; N2 in residue 3; C1 in residue 4}. |
|# COM2 is defined by {C4 in residue 1; O4 in residue 1; N* in residue 1}. |
|# (These definitions are effected by the igr1/igr2 and grnam1/grnam2 |
|# variables; You can use up to 200 atoms to define a center-of-mass |
|# group) |
|# |
|# 2) Set up a distance restraint between COM1 and COM2 which goes from a |
|# target value of 5.0A to 2.5A, with a force constant of 1.0, over steps 1-5000. |
|# |
|# 3) Set up a distance restraint between COM1 and COM2 which remains fixed |
|# at the value of 2.5A as the force slowly constant decreases from |
|# 1.0 to 0.01 over steps 5001-10000. |
|# |
|# 4) Sets up no distance restraint past step 10000, so that free (unrestrained) |
|# dynamics takes place past this step. |
|# |
| |
| &rst iat=-1,-1, nstep1=1,nstep2=5000, |
| iresid=1,irstyp=0,ifvari=1,ninc=0,imult=0,ir6=0,ifntyp=0, |
| r1=0.00000E+00,r2=5.0000,r3=5.0000, |
| r4=99.000,rk2=1.0000,rk3=1.0000, |
| r1a=0.00000E+00,r2a=2.5000,r3a=2.5000, |
| r4a=99.000,rk2a=1.0000,rk3a=1.0000, |
| igr1 = 2,3,4,5,0, |
| grnam1(1)='C1',grnam1(2)='C1',grnam1(3)='N2',grnam1(4)='C1', |
| igr2 = 1,1,1,0, |
| grnam2(1)='C4',grnam2(2)='O4',grnam2(3)='N*', |
| &end |
| &rst iat=-1,-1, nstep1=5001,nstep2=10000, |
| iresid=1,irstyp=0,ifvari=1,ninc=0,imult=0,ir6=0,ifntyp=0, |
| r1=0.00000E+00,r2=2.5000,r3=2.5000, |
| r4=99.000,rk2=1.0000,rk3=1.0000, |
| r1a=0.00000E+00,r2a=2.5000,r3a=2.5000, |
| r4a=99.000,rk2a=1.0000,rk3a=0.0100, |
| igr1 = 2,3,4,5,0, |
| grnam1(1)='C1',grnam1(2)='C1',grnam1(3)='N2',grnam1(4)='C1', |
| igr2 = 1,1,1,0, |
| grnam2(1)='C4',grnam2(2)='O4',grnam2(3)='N*', |
| &end |
| |
+---------------------------------------------------------------------------------------+
\