CHARMM 24 TABLE OF CONTENTS commands.doc .......................................... 1 usage.doc .......................................... 1 charmm.doc .......................................... 1 analys.doc .......................................... 1 block.doc .......................................... 1 changelog.doc .......................................... 1 cons.doc .......................................... 1 corman.doc .......................................... 1 correl.doc .......................................... 1 crystl.doc .......................................... 1 developer.doc .......................................... 1 dynamc.doc .......................................... 1 energy.doc .......................................... 1 ewald.doc .......................................... 1 fourd.doc .......................................... 1 graphx.doc .......................................... 1 hbonds.doc .......................................... 1 hbuild.doc .......................................... 1 images.doc .......................................... 1 install.doc .......................................... 1 intcor.doc .......................................... 1 io.doc .......................................... 1 minimiz.doc .......................................... 1 miscom.doc .......................................... 1 mmfp.doc .......................................... 1 molvib.doc .......................................... 1 monitor.doc .......................................... 1 mts.doc .......................................... 1 nbonds.doc .......................................... 1 nmr.doc .......................................... 1 nose.doc .......................................... 1 parallel.doc .......................................... 1 parmfile.doc .......................................... 1 pdetail.doc .......................................... 1 pert.doc .......................................... 1 perturb.doc .......................................... 1 pimplem.doc .......................................... 1 pressure.doc .......................................... 1 qmmm.doc .......................................... 1 replica.doc .......................................... 1 rism.doc .......................................... 1 rtop.doc .......................................... 1 sbound.doc .......................................... 1 scalar.doc .......................................... 1 select.doc .......................................... 1 struct.doc .......................................... 1 support.doc .......................................... 1 test.doc .......................................... 1 testcase.doc .......................................... 1 travel.doc .......................................... 1 umbrel.doc .......................................... 1 vibran.doc .......................................... 1 CHARMM Element doc/commands.doc 1.1  File: Commands, Node: Top, Up: (chmdoc/charmm.doc), Previous: (chmdoc/parallel.doc), Next: (chmdoc/install.doc) CHARMM commands The commands available for use in CHARMM are classified in several groups. * Menu: * Analysis: (chmdoc/analys.doc ). Analysis facility * Block: (chmdoc/block.doc ). BLOCK free energy simulation * Cons: (chmdoc/cons.doc ). Harmonic and other constraints or SHAKE * Coordinates: (chmdoc/corman.doc ). Commands to manipulate coordinates * Correl: (chmdoc/correl.doc ). Time series and correlation functions * Crystl: (chmdoc/crystl.doc ). Crystal facility * Dynamics: (chmdoc/dynamc.doc ). Dynamics commands * Energy: (chmdoc/energy.doc ). Energy evaluation * Ewald: (chmdoc/ewald.doc ). Ewald summation * Graphx: (chmdoc/graphx.doc ). The graphics subsection for workstations * H-bond: (chmdoc/hbonds.doc ). Generation of hydrogen bonds * H-build: (chmdoc/hbuild.doc ). Construction of hydrogen positions * Images: (chmdoc/images.doc ). Use of periodic or crystal environment * Internal: (chmdoc/intcor.doc ). Manipulation of internal coordinates * I/O : (chmdoc/io.doc ). I/O of data structures and files * Minimiz: (chmdoc/minimiz.doc ). Description of the minimization methods * Miscellany: (chmdoc/miscom.doc ). Miscellaneous commands * MMFP: (chmdoc/mmfp.doc ). Miscelaneous Mean Field Potential * Molvib: (chmdoc/molvib.doc ). Molecular vibrational analysis facility * NMR: (chmdoc/nmr.doc ). NMR analysis facility * Non-bonded: (chmdoc/nbonds.doc ). Generation of the non-bonded interaction * Parameters: (chmdoc/parmfile.doc). CHARMM energy parameters * Perturb: (chmdoc/pert.doc ). Free energy perturbation simulations * Pressure: (chmdoc/pressure.doc). Pressure calculation and usage * Quantum: (chmdoc/qmmm.doc ). Quantum and Molecular Mechanical FF * Replica: (chmdoc/replica.doc ). REPLICA: molecular system replication * RISM: (chmdoc/rism.doc ). Reference Interaction Site Model * Sbound: (chmdoc/sbound.doc ). Stoichastic boundary * Scalar: (chmdoc/scalar.doc ). Scalar command for atom properties * Select: (chmdoc/select.doc ). Use of the atom selection facility * Structure: (chmdoc/struct.doc ). Structure manipulation (PSF generation) * Test: (chmdoc/test.doc ). Commands to test various things * Topology: (chmdoc/rtop.doc ). Residue Topology File * Travel: (chmdoc/travel.doc ). Reaction coordinate refinement command * TSM: (chmdoc/perturb.doc ). Thermodynamic Simulation Method * Umbrella: (chmdoc/umbrel.doc ). Umbrella Sampling * Vibration: (chmdoc/vibran.doc ). Vibrational analysis facility CHARMM Element doc/usage.doc 1.1  File: Usage, Node: Top, Up: (chmdoc/charmm.doc), previous: (chmdoc/install.doc), Next: (chmdoc/support.doc) How to use CHARMM The user of CHARMM controls its execution by executing commands sequentially from a command file or interactivly. In general the ordering of commands is limited only by the data required by the command. For example, the energy cannot be calculated unless the arrays holding the coordinates, the parameters, etc., have already been filled. This section deals with overall usage, as opposed to the detailed description of any given command. This is a good place to start when first learning CHARMM. * Menu: * Meta-Syntax:: Describing the Syntax of Commands * Command Syntax:: Rules for composing command input files. * Run Control:: Ways to modify control flow and stream switching. * I/O Units:: Correspondence between files and unit numbers used by CHARMM. * AKMA:: Units of Measurement used in CHARMM * Data Structures:: Data Structures used by CHARMM * Standard Files:: Descriptions of parameters, topologies, and coordinates available. * Examples:: Sample runs * Interface:: How to make your own private version of CHARMM * Syntactic Glossary:: Glossary of syntactic terms * Glossary:: Glossary of non-syntactic terms.  File: Usage, Node: Meta-Syntax, Up: Top, Next: Command Syntax, Previous: Top Rules for Describing the Syntax (The Meta-Syntax) The syntax of commands is described using the following rules: Capitalized words are keywords that must be specified as is. However, if the word is partially capitalized, it may be abbreviated to the capitalized part. Lower case words are to be replaced by a corresponding data entry. The symbol "::=" means "has the following syntactic form:". Anything enclosed in square brackets, "[]", is optional. If several things are stacked in square brackets, one may choose one optionally. Anything enclosed in curly brackets, "{}", specifies that a selection must be made of the choices stacked vertically inside. The syntactic entities which appear as an argument to "repeat" may be repeated any number (including zero) times. Defaults for optional parameters may be enclosed in apostrophes and placed under the entity they stand for. However, defaults are not specified in this manner if the rules for the default are complex. The syntactic glossary, see *note glossary: Syntactic Glossary, contains further syntactic entities which are used in the command descriptions. Finally, the options and operands in each command can usually be specified in any order except if otherwise noted.  File: Usage, Node: Command Syntax, Up: Top, Next: Run Control, Previous: Meta-Syntax Command language rules and lore A CHARMM run is controlled by a command file (or files). This section of the documentation describes the basic rules for the command file. Details of command level run control are described in the next node. A command file for CHARMM should begin with a specification of the title of the run. (See the syntactic glossary, *note syn: syntactic glossary, for the syntax of a title.) Then, any number of commands may be specified. Each command consists of a command line possibly followed by other data. The command line is scanned free field. This command line may be longer than one line in the file; to do this, one must place a hyphen at the end of line which is to be continued on the next line. Comments may be placed on a command line by preceding the comments by exclamation points. All lower case characters are converted to upper case. This format is identical to that used by the VAX command language interpreter. In addition, blank lines are permitted to separate blocks of commands for increased readability. The first word of every command line specifies the command. Generally, required operands of a command must follow in order. On the other hand, options may generally be specified in any order. Further, any number is always preceded by a key word so that any numeric operands, can be placed in arbitrary order. The command line is scanned in units of words and delimited strings. A word is defined by a sequence of non-blank characters, A delimited string consists of a keyword followed by a string of characters of variable length followed by a delimiter string. One example of where a delimeter string is used is in atom selection where the syntax is; SELE ...... END. Note, that the "END" is required and delimits the atom selection. Abbreviations are permitted in various contexts. The first word may be abbreviated to four characters and numerous options and operands may also be abbreviated to four characters. However, some key words which are used to mark numbers may not be abbreviated. See the processing for individual commands to see what can and cannot be abbreviated. Many of the various options and numeric values are maintained from one invocation of a command to the next. Once a value is specified, it is maintained until it is changed in any command. Therefore, if CUTNB is specified in a NBON command, that value will be used in the DYNA command unless it is changed therein. Usually, when a free field command line is read in, it is echoed onto a standard output. Each such echo will be prepended by a short marker, eg. "CHARMM>", which identifies the line of input as well as the command processor which is interpreting it. In general, as each of the command is interpreted, it is deleted from the command line. When command processing is finished, a check is made to see that nothing is left over. The presence of extraneous junk indicates that something was mistyped. For some commands, such as DYNAmics, where a mistake may be costly, extraneous characters result in a fatal error.  File: Usage, Node: Run Control, Up: Top, Next: I/O Units, Previous: Command Syntax Controlling a CHARMM Run IF command-parameter test-spec comparison-string command-spec GOTO label-string LABEL label-string STREAM [UNIT integer] [file-specification] RETURN SET command-parameter string INCRement command-parameter [BY real] DECRement command-parameter [BY real] This node describes commands that are used to modify the usual sequential interpretation of commands from the command file. Three methods are available to accomplish this: IF tests to conditionally execute a single command GOTO and LABEL transfers within a file STREAM and RETURN transfers to different command files. In addition commands can be modified by the use of command parameters. The command line reader scans input lines for parameters (specified by @n where n is an alphanumeric character) and will subsitute the appropriate parameter string. Command parameters are defined using the SET command to set one of the 36 command parameters, and their values (if numeric) can be modified by the INCRement command, which decodes the parameter string, does real arithmetic and encodes the result. The command parameters are identified by alphanumeric characters (0-9,A(a)-Z(z)(not case-sensitive)). IF compares the string in the specified parameter string to the comparison-string using the test-spec (GT GE EQ NE LE LT). If the comparison is true then the rest of the command line is executed (otherwise it is ignored). The EQ and NE comparisons are done as string comparisons, but the others require decoding of the two strings and comparison by real arithmetic. The command-spec can be any valid command line (including another IF test or a GOTO or STREAM specification). GOTO causes the current command file to be rewound and searched for a line containing the correct LABEL and label-string. The label-string is a single word. If multiple occurrences of a label are present, the first will be used. Command interpretation begins on the line following the LABEL (any information after the LABEL keyword and label-string is ignored). STREAM iunit begins reading commands from the specified fortran logical unit or from the stream file. The stream file is treated exactly as the main command file. It begins with a title and ends with a STOP or RETURN, the latter causing control to return to the previously active command file at the point where the stream switch occurred. The logical unit in OPEN, CLOSE, and REWIND commands are useful in working with streams see *note MISCOM:(chmdoc/miscom.doc). EXAMPLE: * This is a sample command file for CHARMM which calls a stream file * to build a structure and then maps out an adiabatic potential * surface defined by a pair of dihedrals * OPEN UNIT 10 READ FORM NAME makestruc.inp STREAM UNIT 10 SET 1 -180. SET 2 -180. LABEL LOOP CONS CLDH CONS DIHE first-dihedral-angle-spec FORCE 100.0 MIN @1 CONS DIHE second-dihedral-angle-spec FORCE 100.0 MIN @2 MINI minimization-spec INCR 1 BY 30.0 IF 1 LT 170. GOTO LOOP SET 1 -180. INCR 2 BY 30.0 IF 2 LT 170. GOTO LOOP STOP  File: Usage, Node: I/O units, Up: Top, Next: AKMA, Previous: Run Control Fortran I/O Units Usage by CHARMM In order to keep CHARMM as machine independent as possible, all specification of files is done through Fortran unit numbers. Two unit numbers have special signifigance, 5 and 6. Unit 5 is the command file interpreted by CHARMM. Unit 6 is the output file for all printed messages. As commands are read from unit 5, they are echoed on unit 6. All other unit numbers have no predefined meaning. The CHARMM OPEN command may be used to assign files to units. The tream file in "STREAM file-specfication" may be assigned to a logical unit between 100 and 119 (80 and 99 on Cray machines). Logical unit 0 through 9 may be used for CHARMM internal file handling. We recommend logical units 10 through 79 for user data files.  File: Usage, Node: AKMA, Up: Top, Next: Data Structures, Previous: I/O Units The CHARMM system of units: AKMA. CHARMM uses a distinct system of units, the AKMA system. I.e. Angstroms, Kilocalories/Mole, Atomic mass units. All distances are measured in Angstroms, energies in kcal/mole, mass in atomic mass units, and charge is in units of electron charge. Using this system, the AMKA unit of time is 4.888821E-14 seconds (based on the constants tabulated in Abramowitz and Stegun (1970)), however, for all input and output, the time is listed in picoseconds (20 AKMA time units is .978 picoseconds). In some places, the users may specify values in AKMA time units, and in some places both picosecond and AKMA time are output. Angles are given in degrees for the analysis and constraint sections. In parameter files, the minimum positions of angles are specified in degrees, but the force constants for angles, dihedrals, and dihedral constraints are specified in kcal/mole/radian/radian. Any numbers used in the documentation may be assumed to be in AKMA units unless otherwise noted.  File: Usage, Node: Data Structures, Up: Top, Next: Standard Files, Previous: AKMA Data Structures You Should Understand There are a number of data structures that CHARMM manipulates. Many of these data structures are important for most operations; others which are less important, are described with the commands that use them. Much more specific information is available in the various common blocks whose extension is .fcm in the source directory, ~/charmm/source/fcm ([...CHARMM.SOURCE.FCM] on VAX). The important data structures are given below: Each data structure name is followed by its abbreviation which is used as its name in commands. 1) Residue Topology File (RTF) The residue topology file stores the definitions of all residues. The atoms, atomic properties, bonds, bond angles, torsion angles, improper torsion angles, hydrogen bond donors and acceptors and antecedents, and non-bonded exclusions are all specified on a per residue basis. The term "residue" is somewhat historical, but can be any basic unit. 2) The Parameters (PARA or PARM) The parameters specify the force constants, equilibrium geometries, van der Waals radii, and other such data needed for calculating the energy. 3) Structure File (PSF) The structure file is the concatenation of information in the RTF. It specifies the information for the entire structure. It has a hierarchical organization wherein atoms are grouped into residues which are grouped into segments which comprise the structure. Each atom is uniquely identified within a residue by its IUPAC name, residue identifier, and its segment identifier. Identifiers may be up to 4 characters in length. 4) The Internal Coordinates (IC) The internal coordinates data structure contains information concerning the relative positions of atoms within a structure. This data structure is most commonly used to build or modify cartesian coordinates from known or desired internal coordinate values. It is also used in conjunction with the analysis of normal modes. Since there are complete editing facilities, it can be used as a simple but powerful method of examining or analyzing structures. 5) The Coordinates (COOR) The coordinates are the Cartesian coordinates for all the atoms in the PSF. There are two sets of coordinates provided. The main set is the default used for all operations involving the positions of the atoms. A comparison set (also called the reference set) is provided for a variety of purposes, such as a reference for rotation or operations which involve differences between coordinates for a particular molecule. Associated with each coordinate set is a general purpose weighting array (one element for each atom). 6) The Non-bonded List (NBON) The non-bonded list contains the list of non-bonded interactions to be used in calculating the energies as well as optional information about the charge, dipole moment, and quadrapole moments of the residues. This data structure depends on the coordinates for its construction and must be periodically updated if the coordinates are being modified. 7) The Hydrogen Bond List (HBON) The hydrogen bond list contains the list of hydrogen bonds. Like the non-bonded list, this data structure depends on the coordinates and must be periodically updated. 8) The Constraints (CONS) There is a variety of available constraints. All data pertaining to constraints reside in this data structure. 9) The Images data structure (IMAGES) The images data structure determines and defines the relative positions and orientations of any symmetric image of the primary molecule(s). The purpose of this data structure is to allow the simulation of crystal symmetry or the use of periodic boundary conditions. Also contined in this data structure is information concerning all nonbonded, H-bonds, and bonded interactions between primary and image atoms.  File: Usage, Node: Standard Files, Up: Top, Previous: Data Structures, Next: Examples, Files available for general use There are number of residue topology files, parameter files, coordinates files and files of other data structures available. The most important files generally available are residue topology and parameter files. Both such classes of files are stored for general use in the CnnPT: directories. The file names used for both these files consists of an alphabetic part followed by a number, e.g. PARAM7. There are two copies of each file; one with extension, .INP, which is a character files used as an command file to generate the binary file, with extension, .MOD. The .INP is meant for human eyes; the .MOD files is meant for CHARMM to read efficiently. The numeric part of each name is its version number. In general, one should use the highest version number of a file. Although parameter files and toplogy files are separate, they are usually associated, and they must be taken together when generating a structure (PSF). For example, a parameter set for proteins will not work with a DNA topology file. For information on the general use of directories, and the files they contain, see the following sections. * Menu: * Parameters: (chmdoc/parmfile.doc). Description of all the parameter files * Residue: (chmdoc/rtop.doc). Description of the topology files (RTF)  File: Usage, Node: Examples, Up: Top, Previous: Standard Files, Next: Interface Sample CHARMM Runs For an example of specification of a CHARMM run, examine a test case in ~/charmm/test. The file, TEST.INP, is an input to CHARMM which performs the test and contains examples of many commands. The file, TEST.OUT, contains the output from CHARMM produced on Fortran unit 6. Other test cases are found in the test directory.  File:Usage, Node: Interface, Up: Top, Next: Syntactic Glossary, Previous:Examples Interfacing to CHARMM A mechanism has been provided to allow users of the CHARMM to write their own special purpose subroutines which can be incorporated into the system without threatening its integrity. There are six "hooks" into the CHARMM which have been specially provided for casual modifiers. For detailed descriptions of each of these hooks, consult the routine in ~/charmm/source/main/usersb.src on UNIX machines or [...CHARMM.SOURCE.MAIN]USERSB.SRC under VAX/VMS. 1) USERSB The USER command invokes the subroutine, USERSB, and performs no other action. USERSB is a subroutine with no arguments. However, parameters may be passed to this subroutine via the COMMON blocks. These COMMON blocks store nearly all of the systems data. These common blocks may be obtained by including them from the directory containing the sources for the version of the program you are using. 2) USERE A user supplied energy routine may be provided that will be invoked on every energy evaluation. The force arrays should be modified accordingly. 3) USRSEL If one need to be able to select atoms in a manner not possible with the existing options, a user selection routine may be specified. One such example would be for for selecting atoms within a given rectangular solid, or other (nonsperical) solid. 4) USERNM Within VIBRAN, a user specified vector or mode may be generated with this routine. One command that appends this motion onto the existing set of vectors is "EDIT INCL USER integer". 5) USERF A user specified parameter fitting routine may be specified. 6) USRTIM A user specified time series routine may be provied for use in computing correlation functions. To simplify the use of these hooks and to allow users to replace subprograms in the CHARMM with their own versions of said subprograms, the command procedure BUILD has been provided. BUILD will produce a private version of the CHARMM in your default USER directory using your versions of USERSB and USERE. The procedure looks in your directory for USERSB.SRC and USERE.SRC. If either file (or both) is found, it is used in the make procedure of the CHARMM. BUILD command should always be used to generate a private version of the CHARMM as it will always use the correct files for linking. Before attempting to write your own USER functions, you should familiarize yourself with the information available onthe implementation of CHARMM. This interface procedure is designed for short, one time programs. If a user written subroutine is of general use, the routine should be rewritten to conform to parameter passing standards used in the system and then will be incorporated into the central CHARMM. There are several utility routines available to a user routine. Some of them are listed below. CALL GETE(X,Y,Z,...) will cause the energy and forces to be computed and values are saved in the appropriate common blocks. For this to work properly, NBONDS, HBONDS, and CODES must have been called. This can be done by executing both the NBONds and HBONds command, by the use of the UPDAte command, or by having previously found the energy (minimization, dynamics, etc..). CALL PRINTE(...) will write the current energy values (from common block values) to the specified unit (IUNIT). It will also write out the cycle or iteration number and optionally write out the standard header.  File:Usage, Node:Syntactic Glossary, Up:Top, Next:Glossary, Previous: Interface Glossary of Syntactic Terms char A character del The delimiter - a single character which is used to mark the end of a portion of a command. Initially, it is a dollar sign but can be changed using the DELIM command, see *note delim:(chmdoc/miscom.doc). It should be noted that the delimiter cannot be a character within any string it is supposed to delimit. deldel Two delimiters concatened together with no space in between. int or integer An integer iupac IUPAC name for an atom. Initially specified in the residue topology file. keyword A word, see below, serving to identify some option range equivalent to real real integer. The first real is the minimum value in the range, the second number is the maximum value in the range, and the third number gives the number of interval, i.e. lines or columns. real A real number. No decimal point is required for the number to be interpreted correctly resid Residue identifier (a string of upto 4 characters) resname Residue name (type of residue. e.g. GUA) segid Segment identifier (a string of upto 4 characters) string An ordered set of characters tag A string which is a tag, i.e. no embedded spaces. title A series of 1 to 32 lines of text (max 80 characters per line) each starting with a "*". The title is terminated by a line which an asterisk "*" as the first character. Used for commenting files. word A string with no blanks unit-number An integer which is a Fortran unit number.  File: Usage, Node: Glossary, Up: Top, Previous: Syntactic Glossary, Next: Top General Glossary data structure A collection of arrays, scalars, and possibly other data structures which are related by part of a larger entity. For example, a coordinate set is a data structure which hold the three dimensional positions of atoms. This data structure consists of 1 scalar and three arrays. The scalar is the number of coordinates; the three arrays are the X, Y, and Z components of the coordinates. Internal bonds, angles, torsions, improper torsions. coordinates Also, a data structure used for constructing coordinates. Iupac Name for The name of an atom with a residue. This name should be an atom unique within a residue and should conform to the IUPAC nomenclature, Biochemistry 9:3471 (1970) Hbonds hydrogen bonds Parameters constants in the energy expression ( force constants, minima of energy surfaces, charges, Lennard-Jones parameters, van der Waals radii, etc.) PSF structure file ( protein structure file ) : a list of the internal coordinates and related information Residue A string of four characters or less which uniquely specifies Identifier residue with in a segment. This value is currently set by CHARMM to be the character representation of the residue number in the segment starting from the first real monomer unit in it. RTF residue topology file : a list of standard internal coordinates, atom charges, atom types, excluded non-bonded interactions, etc. Segment A string of up to four characters uniquely designating Identifier a segment. Specified in the GENErate command, see *note gener: (chmdoc/struct.doc) Generate. Sequence list of residues CHARMM Element doc/charmm.doc 1.1  File: CHARMM, Node: Top Chemistry at HARvard Macromolecular Mechanics - --- - - Version 24b1 - August 15, 1995 Copyright(c) 1984,1987,1991,1994,1995 President and Fellows of Harvard College All rights reserved You are now using the INFO facility to view CHARMM 24 documentation. The paper; CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comp. Chem., Vol. 4, p187 (1983), is considered to be an integral part of this documentation. In places, this documentation and the paper will conflict. In all such cases, the documentation presented here should take precedence. * Menu: * Commands: (chmdoc/commands.doc). Discription and syntax of CHARMM commands * Install: (chmdoc/install.doc). Release notes How to install CHARMM on a user site * Usage: (chmdoc/usage.doc). How to use CHARMM * Support: (chmdoc/support.doc). Supporting data files and utilities * Testcase: (chmdoc/testcase.doc). CHARMM testcases * Develop: (chmdoc/developer.doc). Notes for CHARMM developers * News: (chmdoc/changelog.doc). New features and Modifications * Parallel: (chmdoc/parallel.doc). CHARMM on parallel platforms * Info: (Info). A description of the INFO facility. CHARMM Element doc/analys.doc 1.1  File: analys, Node: Top, Up: (chmdoc/commands.doc), Next: Description Analysis Commands In CHARMM23, new analysis commands are under development, but some features are currently functional. * Menu: * Description:: Description of analysis facility * Energy:: Energy partitioning  File: analys, Node: Description, Up: Top, Previous: Top, Next: Energy Description of the ANALysis Command ANALys {ON } Enable analysis and disable FAST routines. {OFF} Disable analysis and restore FAST option defaults. The ANALysis command is a new energy and structure table facility that is being developed to examine both static and dynamic properties. The current code only allows energy partition analysis and energy contribution analysis from free energy simulations. It also can produce a detailed printout of energy term contributions.  File: analys, Node: Energy, Up: Top, Previous: Description, Next: Top Energy option of the ANALysis Command Syntax: ANALys { ON } { TERM { [ALL] } { NONBond } [UNIT int] atom-selection } { { ANY } { [NONOnbond] } } { OFF } The ANALysis ON command enables energy partition analysis and disables the FAST routines. This will slow the calculation (especially on vector machines), but allow a detailed, atom by atom, energy analysis. Everytime the energy routine is invoked, the energy for each atom is stored in the ECONT array. During PERT dynamics, the EPCONT is filled with the time average energy difference on a atom by atom basis including every step of dynamics. This allows the free energy differences to be analyzed based on atom contributions. The ANALysis TERM command will cause all selected energy terms to be printed to the specified output unit (default: standard CHARMM output). The ALL keyword (default) will list elements where all atoms are selected. The ANY keyword will cause terms including any selected atom to be listed. The NONOnbond keyword (default) suppresses the listing of vdw and electrostatic energy terms. The NONBonded keywords also allows the analysis of vdw and electrostatic interactions. The ANALys OFF command enables the FAST routines and disables the resetting of the ECONT array (i.e. the ECONT array will not change, but may still be accessed. This command also disables the energy term analysis. The energy partition array can be accessed with the SCALar ECONt commands. *note Econt:(chmdoc/scalar.doc). The sum of all of the elements of the ECONT array is usually the total energy, but some energy terms, such as extended electrostatics, will not be included. The command: SCALar ECONT STATistics can be used to check the total energy and the command SCALar EPCONT .... can be used to examine atom contributions to energy differences for PERT. CHARMM Element doc/block.doc 1.1  File: BLOCK, Node: Top, Up: (chmdoc/commands.doc), Next: Syntax The commands described in this section are used to partition a molecular system into "blocks" and allow for the use of coefficients that scale the interaction energies (and forces) between these blocks. This has a number of applications, and specific commands to carry out free energy simulations with a component analysis scheme have been implemented. BLOCK was recently modified so that it works with the IMAGE module of CHARMM. As some changes to the documentation were necessary anyways, it was decided to also improve the existing documentation. The Syntax and Function section below are relatively unchanged; the added documentation is in the Hints section (READ IT if you are using BLOCK for the first time!). Comments/suggestions to boresch@tammy.harvard.edu. * Menu: * Syntax:: Syntax of the block commands * Function:: Purpose of each of the commands * Hints:: Some further explanations/hints * Limitations:: Some warnings...  File: BLOCK, Node: Syntax, Up: Top, Next: Function Syntax of BLOCK commands BLOCk [int] Subcommands: miscellaneous-command-spec ! see *note miscom:(chmdoc/miscom.doc). CALL int atom-selection LAMBda real COEFficient int int real NOFOrce FORCe FREE_energy_evaluation [OLDLambda real] [NEWLambda real] - FIRSt int [NUNIT int] [BEGIn int] [STOP int] [SKIP int] - [TEMPerature real] [CONTinuous int] [IHBF int] [INBF int] [IMGF int] INITialize CLEAr Energy_AVeraGe [OLDLambda real] [NEWLambda real] - FIRSt int [NUNIT int] [BEGIn int] [STOP int] [SKIP int] - [CONTinuous int] [IHBF int] [INBF int] [IMGF int] COMPonent_analysis DELL real NDEL int [TEMPerature real] - FIRSt int [NUNIT int] [BEGIn int] [STOP int] [SKIP int] [IHBF int] [INBF int] [IMGF] int AVERage {DISTance int int} {STRUcture} [PERT] [TEMPerature real] [OLDLambda real] [NEWLambda real] - FIRSt int [NUNIT int] [BEGIn int] [STOP int] [SKIP int] END  File: BLOCK, Node: Function, Up: Top, Previous: Syntax, Next: Hints 1) BLOCk [int] enters the block facility. The optional integer is only read when the block structure is initialized (usually the first call to block of a run) to specify the number of blocks for space allocation. If not specified, the default of three is assumed. 2) END exits the block facility. The assignment of blocks, the coefficient weighting of the energy function, the force/noforce option, etc. remain in place. For the terms of the energy function that are supported, each call to ENERGY (either directly or through MINIMIZE, DYNAMICS, etc. commands) results in an energy and force weighted as specified. The matrix of interaction coefficients is printed upon exiting. 3) CALL removes the atoms specified by "atom-selection" from their current block and assigns them to the block number specified by the integer. Initially all atoms are assigned to block 1. If atoms are removed from any block other than block 1, a warning message is issued. If blocks are assigned such that some energy terms (theta, phi, or imphi) are interactions between more than two blocks, a warning is issued when the END command is encountered. (Take such warnings seriously; this is a severe error and indicates that something is wrong. However, the problem might be not the CALL statement (or the atom selection) itself; quite possibly your hybrid molecule was generated impromperly) 4) LAMBda sets the value of lambda to "real". This command is only valid when there are three blocks active. Otherwise multiple COEF commands may be used to set the interaction coefficients manually. LAMBda x is equivalent to (let y=1.0-x) COEF 1 1 1.0 COEF 1 2 y COEF 1 3 x COEF 2 2 y COEF 2 3 0.0 COEF 3 3 x 5) COEF sets the interaction coefficient between two blocks (represented by the integers) to a value (the real number). When the block facility is invoked, all of the atoms are initially assigned to block 1 and all interaction coefficients are set to one. 6) NOFOrce specifies that in subsequent energy calculations, the forces are not required. This is economical when using the post-processing commands (FREE,EAVG,COMP). Forces may be turned back on with the FORCe command; this is necessary before running minimizations and dynamics if there was a prior NOFO command. 7) FREE calculates a free energy change using simple exponential averaging, i.e. the "exponential formula". If the old and new lambdas (OLDL,NEWL) are specified (can only be done when three blocks are active), the perturbation energy is calculated for these values (i.e. FREE gives you the free energy difference between NEWLambda and OLDLambda via perturbation from the lambda value at which your trajectory was calculated. If not, the current coefficient matrix is used (FREE should be used with three blocks, and the use of OLDL and NEWL is recommended). FIRSt_unit, NUNIt, BEGIn, STOP, and SKIP specify the trajectory/ies that is/are to be read (for a further description see the TRAJ command elsewhere in the CHARMM documentation). TEMPerature defaults to 300 K and gives the temperature value to be used in k_B*T. CONTinuous specifies the interval for writing cumulative free energies. A negative value causes binned (rather than cumulative average) values to be written. Be careful to make sure that you use correct non-bonded lists (see the hints section!) 8) INITialize is called automatically when the BLOCK facility is first entered and may also be called manually at some other point. All atoms are assigned to block one and all interaction coefficients are set to their initial value. 9) CLEAr removes all traces of the use of the BLOCK facility. The next command should generally be END, and then CHARMM will operate as if BLOCK had not ever been called. 10) EAVG The average value of the potential energy during a simulation can be calculated with the EAVG (Energy_AVeraGe) command. The parsing is very much like the FREE command above. The most frequent use of this command is to calculate the average value of dV/dlambda during the course of a simulation for use in thermodynamic integration. CONTinuous specifies the interval for writing cumulative free energies. A negative value causes binned (rather than cumulative average) values to be written. Be careful to make sure that you use correct non-bonded lists (see the hints section!) The command accepts the OLDL / NEWL option, similarly to FREE, but for EAVG it is recommended to set up the interaction matrix (using COEF commands) yourself -- see the hints section. 11) [COMP] The COMP module is essentially a modified version of the EAVG module which aside from calculating = at a given value of lambda l(i) will also give you expectation values of this quantity at l(i+-1), l(i+-2) etc. based on perturbation theory. COMP requires 4 blocks. Put the usual WT (reactant) in block 2 and MUT (product) in block 3. Put the portion of the environment whose contribution to the free energy change is desired into block 4 (this can be everything else, or just a subset) (Note that the same can be achieved easily with the EAVG command) You have to set up your own coefficient matrix. Much of the parsing is like the EAVG command. CONT is not supported. Two special subcommands (required) are DELL and NDEL. The normal output of COMP is evaluated at the lambda of the simulation. However, COMP also evaluates the same ensemble averages perturbed to lambda = lambda +/- {0,1,2,...NDEL}*DELL. This (sometimes) helps the quadrature in thermodynamic integration. Note that NDEL must be at least 1, and DELL should not be zero. (You have to specify these values; the default values will lead to an invalid input, i.e. you bomb...) Be careful to make sure that you use correct non-bonded lists (see the hints section!) A word of warning: If your initial ensemble average (at the lambda of the simulation) is not well converged, then your perturbed values are most likely random numbers. The approach taken by COMP is theoretically sound, but it should only be applied if convergence has been established! The output format of COMP is somewhat messy: COMP first prints = at lambda = lambda - NDEL*DELL lambda - (NDEL-1)*DELL ... lambda lambda + DELL ... lambda + NDEL*DELL; then it prints an average (integral) value over these results. The meaning of this last value is unclear to me. In earlier versions of this documentation, COMP has been recommended over EAVG. In my experience the opposite is true. There is little COMP can do which you can't do with EAVG (aside from obtaining expectation values for dU/dl). (Maybe the unclear output of the COMP module is the main reason why I don't like it). 12) [AVER] The AVERage command is used to extract ensemble average structural properties from a dynamics simulation. Features in this implementation allow averages taken over ensembles that are perturbed from that which the simulation corresponds to. This is particularly useful for calculating the average structure expected at lambda=0.0 from a simulation run at lambda=0.1, for example. One may calculate average structures [STRUcture] and average distances [DISTance int int; where the two integers are the atom numbers between which the average distance is requested], currently. The PERT keyword indicates that a perturbed ensemble from the dynamics trajectory is desired, with TEMPerature giving the temperature to use in the exponential for the perturbation (defaults to 300 K), OLDLambda and NEWLambda are the lambdas for which the simulation was run and for which the ensemble is requested, respectively (only valid if three blocks are active; if these are not specified, the perturbation energy is calculated with the current coefficient matrix), and the remaining keywords are used to specify the trajectory. NOTE: TO THE BEST OF MY KNOWLEDGE THIS COMMAND HAS NOT BE MAINTAINED (so you are on your own if you use it!)  File: BLOCK, Node: Hints, Up: Top, Previous: Function, Next: Limitations A warning is in order: the BLOCK module is quite user-unfriendly, AND the user (=you) has to know what he/she is doing, otherwise you won't get anywhere! (Of course, this could be a blessing in disguise) There are two applications for BLOCK: (i) Mere use as an energy partitioning facility, which may, e.g., very helpful as an alternative to the INTEraction energy command and (ii) use in free energy simulations. The focus here is on free energy applications. The following paragraphs assume that you are familiar with the theory of free energy difference simulations (e.g. Brooks et al. Advances in Chem. Physics, Vol. LXXI, 1988, chapter V); the emphasis here is to show how a rough tool as BLOCK can be used to implement the theory in a program and (of course) how to use it. Using BLOCK in order to calculate a free energy difference consists out of two rather dissimilar parts (as far as practical problems are concerned): (i) Run your system at various values of lambda and save trajectories. (ii) Postprocess the trajectories with the FREE or the EAVG command (possibly COMP), use the quantities which these modules give you to calculate the free energy difference. (i) The actual simulations ========================== It's probably easiest to use a concrete example, and the free energy difference between ethane and methanol in aqueous solution is used for that purpose. BLOCK is a so-called dual topology method (D. Pearlman, JPC 1994, 98, 1487) i.e. one has to duplicate any atom that is different with respect to any of its parameters. In the ethane/methanol case this means that you have to run with a solute which looks something like H1 \ /H4 \ C1E ---- C2-H5 H2 = { } \H6 / C1M --- OG / \HG1 H3 (and there is water.) Conceptually, this system is divided into three regions: environment: water, H1, H2, H3 (the region where nothing changes) reactant: C1E, C2, H4, H5, H6 (ethane half) product: C1M, OG, HG1 (methanol half), where of course the role of reactant and product is interchangeable. The steps involved to start running dynamics are as follows: (1) set up the hybrid (generate psf). In principle straightforward, but there is a practical pitfall: The autogenerate angles and dihedrals option(s) may produce artificial dihedrals/angles between the two/three parts of the system, e.g. you don't want angles H1-C1E-OG etc. or dihedrals H3-C1M-C2-H4 etc. Also, make sure to specify nonbonded exclusions between the reactant and product part, otherwise you'll get endless distance warnings and may even bomb if two atom positions coincide. (2) Place the hybrid into water (stochastic or periodic boundary conditions -- yes, IMAGE is now supported) as usual (3) Partition the system, i.e. enter BLOCK The following script fragment will do the trick: block 3 call 2 sele end call 3 sele end end (reactant and product have to be defined according to your system). BLOCK 3 initializes the block module with 3 blocks, all atoms are in block 1. The two CALL commands bring the reactant and product part of the system into block 2 and 3 respectively. (4) Run the necessary MD simulations. Let's assume that you decide to use the following values of lambda, lambda = 0.1, 0.3, 0.5, 0.7, 0.9. You want to start your simulation at lambda = 0.1 and you have already partitioned your system as shown in (3). (This information is kept within the same script between calls to block, but it is not saved in restart files or the psf, i.e. you have to repeat this step (as well as step (3)) in every input file). Enter block again, e.g. block lamb 0.1 end From now on interactions between the 3 blocks will be scaled according to the following matrix (lambda = l = 0.1 ==> 1-l = 0.9): block | 1 2 3 ------|-------------------- 1 | 1.0 1-l l 2 | 1-l 1-l 0. 3 | l 0. l Please note that BLOCK will first calculate an interaction, then check to which block the two atoms belong and scale the energy (and forces) appropriately. Therefore, if the distance between 2 atoms is zero (e.g. in the ethane/methanol example I would define C1M and C1E on top of each other!) then you need non-bonded exclusions, otherwise you encounter a division by 0 error! The LAMB command is a shortcut for the following sequence of COEF commands, the following code fragment should be self-explanatory: block coef 1 1 1.0 coef 1 2 0.9 coef 1 3 0.1 coef 2 2 0.9 coef 2 3 0.0 coef 3 3 0.1 end BLOCK only accepts and uses symmetric matrices, i.e. it doesn't matter whether you specify COEF 1 2 or COEF 2 1. Whenever you now call the energy routines, the energies/forces returned from them will be scaled according to the matrix you have set up. Minimizers and Dynamcis can be used as always. So you are ready to run dynamics, and for arguments sake say that you run at every value of lambda 10,000 steps equilibration and 20,000 steps production (i.e. you save coordinates to trajectories) You don't need to save every step, every 5th to 20th step is probably more than enough. (If you saved every step you'd obtain highly correlated data, i.e. you have larger trajectories, but you won't gain anything in terms of convergence.) (ii) Post-processing -- how to obtain a free energy difference ============================================================== At this point in our example, you would have five trajectories corresponding to lambda = 0.1, 0.3, ..., 0.9 The BLOCK module now has to be used to obtain the average quantities you need for either the exponential formula (FREE) or thermodynamic integration (EAVG,COMP) from the trajectories you generated in step (i) (1) At this point, some issues regarding the non-bonded list have to be considered. No special considerations were necessary while running dynamics (aside from having some non-bonded exclusions where necessary); you just set up list updates as usual. During post-processing there are two considerations: (a) efficiency -- you just want to calculate the necessary subset of interactions (otherwise your post-processing run will take about as much time as the simulation itself), and (b) proper list-updating. (a) Efficiency: In none of the post-processing routines do you need the interactions between particles that belong to the environment; therefore you should avoid calculating them. This can be done easily by specifying cons fix sele end Note that this is not necessary, but it will reduce the CPU time necessary from hours to minutes (and results are identical!) However, if you had atoms belonging to reactant or product or both FIXed during the simulations in step (i), you MUST NOT FIX them now; otherwise you'll omit contributions. (b) List updating: While the efficiency considerations in principle are optional, you have to follow one of the two strategies below otherwise you'll get erroneous results. If you used IMAGE, you have to use the second protocol! Originally, the BLOCK post-processing commands would not do any list updating. This meant that you had to have a nonbonded list which would include all possible interactions before starting post-processing -- don't forget that you post-process over, e.g., 20 ps and particles will move quite far. You can easily create such a nonbonded list by specifying a CUTNB value of, e.g. 99. or 999. Ang (surely, all possible interactions will be included). A CHARMM script looks approximately as follows: !set up system (psf, initial coordinates) block !partition system end cons fix sele end ==> energy cutnb 99. !open trajectories block !postprocessing end In this case, do not use the inbf, ihbf and imgf options of the post-processing commands, they will default to 0, i.e. no update. This approach, however, CANNOT work with IMAGES! Proper use of IMAGEs requires that the minimum image convention is checked periodically, i.e. particles have to be repartitioned between primary and image region. As the BLOCK post-processing commands now understand INBF, IHBF and IMGF, this doesn't pose a problem. However, the automated update is not supported (if you specify a negative value, you'll get a mild warning and the system will default to +1), and I recommend that you use 1 for all frequencies (don't forget, the frames in your trajectory are several steps apart, i.e. in general an update may be necessary) The above scheme now looks like: !set up system (psf, initial coordinates) block !partition system end cons fix sele end ! set up images if needed ==> energy !open trajectories block eavg inbf 1 ihbf ? (imgf 1) end Unless you have explicit hbond terms, ihbf can of course be 0! (Please note that there may or may not be problems with CRYSTAL, see Limitations section) (2) The actual post-processing commands. In the following I'll explain how to set things up for FREE, EAVG and COMP (as well as why). To speed up things further, you'll also want to specify the NOFOrce option at some point. FREE: This module allows you to calculate the necessary ensemble average for the exponential formula. Using our example, you can for example estimate the free energy difference between l=0.1 (a value at which you ran a trajectory) and l=0.0, or, based on your l=0.1 trajectory the free energy difference between l=0.0 and 0.2 (double wide sampling), i.e. A(0.0)-A(0.1) = -k_B*T*ln _(l=0.1) or A(0.2)-A(0.0) = -k_B*T*ln _(l=0.1) You should set up your system with 3 blocks and the usual environment, reactant and product partitions. Before entering block to issue the free command, you have to open the trajectory/ies. ! all the stuff shown above for non-bond lists open file unit 10 read name dat01.trj block free oldl 0.1 newl 0.0 first 10 nunit 1 [temp 300. - inbf 1 imgf 1] end or, for double wide sampling, the free line would be replaced by free oldl 0.0 newl 0.2 first 10 nunit 1 [temp 300. - inbf 1 imgf 1] Here dat01.trj is the trajectory which contains your 20 ps of dynamics at lambda = 0.1. Based on the oldl/newl values (which correspond to A(newl) - A(oldl)), FREE generates the appropriate interaction matrix, which it prints; I recommend that you try to understand why it generates this matrix! FIRST is the unit number of the first trajectory file (10 in our example), NUNIT is the number of trajectories (1 in our example). These (and the other options regarding the trajectories work as in any other post-processing command in CHARMM, see e.g. the TRAJ command) The update frequencies are optional depending on how you decided to handle your non-bonded updates. temp defaults to T=300 K, cf. equations above. If you specify CONT +n, you'll get a cumulative average every n steps; in this case the last value equals the final result; if you specify CONT -n, you'll get the average over every n frames, plus of course the final result at the end. Note that trajectories are not rewound after use; i.e. before any subsequent FREE (or EAVG,COMP) command you have to rewind (or reopen) them! Once you have all the free energy pieces you need, you simply add them up to obtain the free energy difference (beware of sign errors depending on how you defined oldl/newl) EAVG: The main use of this module lies in obtaining the required ensemble averages for thermodynamic integration. The most significant difference to EAVG is that you have to specify your own interactions matrix. BLOCK uses linear coupling in lambda in the potential energy function, i.e. V(l) = V0 + (1-l)*V_reac + l*V_prod, where V0 contains all the intra-environment terms, V_reac are the intra-reactant and reactant-env. interactions, and V_prod are the intra-product and product-env. interactions, respectively. The quantity of interest in TI is dV/dl; for the above potential energy function we have dV/dl = V_prod - V_reac It's very easy to obtain this quantity from EAVG. Use 3 blocks, partition the system as before. ! all the stuff shown above for non-bond lists open file unit 10 read name dat01.trj block coef 1 1 0. coef 1 2 -1. coef 2 2 -1. coef 1 3 1. coef 2 3 0. coef 3 3 1. eavg first 10 nunit 1 [inbf 1 imgf 1 cont +-n] end You will calculate the average interaction energy over all the frames in the trajectory according to the following (symmetric) matrix 0.0 -1.0 -1.0 1.0 0.0 1.0; i.e. it's easy to see that the above script will give you _(l=0.1). If you post-process the other trajectories (l=0.3, 0.5, ..,0.9) in an analogous fashion, you just have to approximate the TI integral by the trapezoidal formula (for basic Newton Cotes formulae (open and closed) see, e.g., Numerical Recipes), i.e. in this case you would have dA = 0.2 * (dV(0.1)+dV(0.3)+...+dV(0.9)), where dV(0.1) = _(l=0.1), etc. The above is an example of the basic use of EAVG. You automatically get the formal components according to interaction type. Cont +-n works similarly to the FREE case. If you wanted to exclude the intramolecular contributions from ethane and methanol you could set up a slightly different coefficient matrix, i.e. coef 1 1 0. coef 1 2 -1. coef 2 2 0. coef 1 3 1. coef 2 3 0. coef 3 3 0. and you'll get only the solute-solvent contributions. You can use more blocks (m > 3) to extract only a subset of interactions, e.g. block 1: environment - x block 2: reactant block 3: product block 4: x, where x is the region of interest, e.g. a specific sidechain in a protein (but not the one that is mutated!) Using EAVG with an appropriate coefficient matrix, e.g. coef 1 1 0. coef 1 2 0. coef 1 3 0. coef 1 4 0. coef 2 2 0. coef 2 3 0. coef 2 4 -1. coef 3 3 0. coef 3 4 1. coef 4 4 0. will give you (after integration over lambda) the free energy contribution of the interaction of sidechain x with the mutation site. Note that such formal free energy components may be (strongly) path-dependent. These last two examples have hopefully provided a flavor of what can be done with the EAVG module. COMP: This module is also used for thermodynamic integration. It always operates with four (and only four) blocks, just as the advanced example last given for EAVG, so it facilitates COMPonent analysis. Here I want to focus on the second unique aspect of COMP, it's capability to extrapolate additional datapoints, and so I consider in the framework of our ethane/methanol example the "special" case where I want the total free energy difference (as before in EAVG). In order to do this, the system needs to be partitioned as follows block 1: -- block 2: reactant block 3: product block 4: environment Whereas EAVG gave us _l only for those lambda values at which we had actually done the simulations, COMP gives us additional values via perturbation (see Bruce Tidor's thesis). Using ! all the stuff shown above for non-bond lists open file unit 10 read name dat01.trj block coef 1 1 0. coef 1 2 0. coef 1 3 0. coef 1 4 0. coef 2 2 -1. coef 2 3 0. coef 2 4 -1. coef 3 3 1. coef 3 4 1. coef 4 4 0. comp first 10 nunit 1 [inbf 1 imgf 1] dell 0.06667 ndel 1 end will now give us _l at l=0.03334, l=0.1 and l=0.16667. If we use the same script on the other trajectories, we have 15 instead of 5 datapoints for the integration, i.e. we can obtain dA as dA = 0.06667 * (dV'(0.03334)+dV(0.1)+...+dV'(0.96667)), where dV(0.1) = _(l=0.1), etc. and the ' indicates that this is a perturbed quantity. In principle, this should give a better numerical integration; however, in practice everything depends on how well your actual data (l=0.1, 0.3, ...,0.9) are converged. There is no check whether your ndel/dell combination is meaningful; and you cannot run COMP without using the perturbation feature, i.e. NDEL should be set to at least 1 (valid values are 1 through 5). The defaults (if you don't specify ndel/dell) lead to an invalid input (This should be fixed...)  File: BLOCK, Node: Limitations, Up: Top, Previous: Hints (1) Please be advised (again) that the AVERage command is unsupported, and I would not be surprised if it does not work (anymore). Unless someone who understands this module better than I do maintains it, I recommend that we remove it. (2) BLOCK now coexists with IMAGE "peacefully" and essentially transperantly to the user. It works correctly for the case of a periodic water-box (cf. the block3.inp testcase). I would, however, check carefully whether things really work before I would use it on something fancier like infinite alpha helices. Similarly, it is not clear to me whether things work with the CRYSTAL facility. If one modifies block3 as to use CRYSTAL instead of IMAGE things (seem to) work. On the other hand, I know that I didn't support XTLFRQ in the post-processing routines as I don't understand its meaning. I'll fix things if someone is willing to help me with the bits and pieces I don't understand. (3) Bond and bond angle terms (including Urey-Bradleys). Be advised that if you run a simulation at lambda = 0 or lambda = 1 you may effectively remove bond (and bond angle terms) as they get scaled by zero. In other words, you would have ghost particles that can move freely through your systems, and this leads to all sorts of nasty side-effects. Furthermore, this approach is not sound theoretically (S. Boresch & M. Karplus, unpublished). So in general, avoid running at lambda = 0 and 1. If you have your bonds constrained you're safe as the constraint will keep things together (that won't take care of angles however!) In order to avoid artifacts from noisy, diverging bond and bond angle contributions throw them out during post-processing, e.g. by using the SKIP BOND ANGL UREY command before starting block post-processing. If you want to see what can go wrong, look at the block2 test-case... CHARMM Element doc/changelog.doc 1.1  File: ChangeLog, Node: Top, Up: (chmdoc/charmm.doc), Previous: (chmdoc/developer.doc), Next: (chmdoc/parallel.doc) CHARMM Developer's Change Log Entries in C21-C22 are recorded by CHARMM developers to indicate new and modified features before the CHARMM 22.0b release (January 1, 1992). C20-C22 and C22-C23 summarize new features and modifications carried out during C22 development (till c22g2 release, August 15, 1992) and C23 development (till c23f3 release, February 1, 1994), respectively. C23-C24 is the record of developments during the period of c24a1 (February 15, 1994) through c24b1 (August 15, 1995). ------------------------------------------------------ CHARMM22.0.b Release April 22, 1991 CHARMM22.0.b1 Release September 30, 1991 CHARMM22 Release January 1, 1992 c22g1 Release February 15, 1992 c22g2 Release July 7, 1992 c22g3 Release November 3, 1992 c22g4 Release March 1, 1993 c22g5 Release August 1, 1993 CHARMM23.0 c23a1 Developmental August 15, 1992 c23a2 Developmental October 25, 1992 c23f Developmental March 1, 1993 c23f1 Developmental March 15, 1993 c23f2 Developmental August 15, 1993 c23f3 Release February 1, 1994 c23f4 Release August 15, 1994 c23f5 Release March 15, 1995 CHARMM24.0 c24a1 Developmental February 15, 1994 c24x1 Evaluation February 15, 1994 c24a2 Developmental August 15, 1994 c24a3 Developmental March 15, 1995 c24b1 Release August 15, 1995 ------------------------------------------------------ * Menu: * C21-C22:: Modifications of Developmental CHARMM21 to CHARMM22 * C20-C22:: Major enhancements and developments in CHARMM22 * C22-C23:: Major enhancements and developments in CHARMM23 * C23-C24:: Major enhancements and developments in CHARMM24  File: ChangeLog, Node: C21-C22, Up: Top, Previous: Top, Next: C20-C22 Summary of Modifications of Developmental CHARMM21 to CHARMM22 ------------------------------------------------------------------------------ Linear pressure ramping added to CPT code (see pressure.doc) ------------------------------------------------------------------------------ Frequency based crystal update is now supported Relevent new keyword is IXTFrq (see image.doc) ------------------------------------------------------------------------------ Constant Pressure and Temperature (CPT) dynamics (See PRESSURE.DOC) TRICLINIC unit cell is now supported. ------------------------------------------------------------------------------ Miscellaneous commands: UPPEr and LOWEr keywords added (see miscom.doc) ------------------------------------------------------------------------------ Minimization: new keyword (FMEM) for ABNER minimizer (see minimiz.doc) ------------------------------------------------------------------------------ Internal coordinates (see INTCOR.DOC) New commands: IC SAVE IC RESTore IC RANDom [iseed] Internal coordinates converted to double precision. ------------------------------------------------------------------------------ Coordinate Manipulation (See MISCOM.DOC and CORMAN.DOC) New inline command varaibles added: ?THETa , ?XMOVe , ?YMOVe , ?ZMOVe , ?RMS New CORMAN commands added: COOR HELIx COOR PUCKer COOR COVAraince COOR SEARch ... RBUFF ... ------------------------------------------------------------------------------ Energy, Angles Urey-Bradley 1-3 terms have been added as an option. Format of parameter file affected. (See IO.DOC) Energy analysis code added (ANALysis ON command). (See ANALYS.DOC) ------------------------------------------------------------------------------ NOE distance restraints (See CONS.DOC) Overhaulled to become a general distance restraint term. Commands syntax overhaulled as well. ------------------------------------------------------------------------------ PSF common structure modified Unused PSF arrays removed. All size limits increased. Binary file format changed to INTEGER*4 and REAL*8 PSF numbers added to ?variable list (See MISCOM.DOC). ------------------------------------------------------------------------------ Output redirecting implemented. (See MISCOM.DOC) OUTU replaces all writes to unit 6. ------------------------------------------------------------------------------ ATLIM modified to allow a limit of several days. PASMID has been changed to an integer which points the current day. See MISCOM.DOC ------------------------------------------------------------------------------ Free energy perturbation commands added. (See PERT.DOC) Several new commands and features have been modified to allow free energy perturbation simulations to be performed. ------------------------------------------------------------------------------ Partition function and classical free energy codee added to the vibrational analysis code. (See VIBRAN.DOC) Atom selection added for EDIT commands. Atom selection added for WRITE SECOnd-derivatives CARD command. ------------------------------------------------------------------------------ New time series commands and options (See CORREL.DOC) ENTER PUCKer ENTER HELIx ENTER RMS ENTER ENERgy ENTER RMS [MASS] atom-selection ENTER ATOM CROSsproduct ENTER FLUC CROSsproduct ENTER VECT CROSsproduct ENTER HBOND ENTER MODE ENTER RMS [MASS] [ORIEnt] ... TRAJ ... atom-selection MANTIME SQUARE (vectors now allowed) MANTIME ABS (vectors now allowed) MANTIME ACOS Off-by-one error removed in time series data (time series now do not start at time zero, but at time DELTA*SKIP). ------------------------------------------------------------------------------ Langevin dynamics modified. An improved algorithm has been incorporated which gives a more accurate integration at low gamma values as well as the proper brownian dynamics limiting values in the large gamma limit (and is more efficient). The gaussian random generator has been replaced to give a much more accurate distribution and uses only one random number call per atom by using an error function lookup table. ------------------------------------------------------------------------------ Miscellaneous commands added. (See MISCOM.DOC) DIVIde, EXONent, RANDom, and SHOW New miscellaneous variables added. ?RAND, ------------------------------------------------------------------------------ Precision and index limits improved. The entire program (except for the graphics section) has been converted to REAL*8 and INTEGER*4 from REAL*4 and INTEGER*2. ------------------------------------------------------------------------------ Constant Pressure and Temperature (CPT) dynamics added. (See PRESSURE.DOC) Pressure analysis code added. NTRFRQ usage modified so that it works for IMAGES and CRYSTAL. ------------------------------------------------------------------------------ Heuristic nonbond update feature added. (See NBONDS.DOC) ------------------------------------------------------------------------------ New (consistent) energy print format with search line indicators. ------------------------------------------------------------------------------ Graphics subsection added for workstations. ------------------------------------------------------------------------------ New GRADient option added for most minimization methods for searching for saddle points. ------------------------------------------------------------------------------ FAST option is now the default. It is no longer necessary to have the command "FAST 1" in order to use the efficient energy routines. ------------------------------------------------------------------------------ Constrained reference now only set for selected atoms for the CONS HARMonic command (the old method limited versatility). (See CONS.DOC) ------------------------------------------------------------------------------ Parallelization for shared memory multi-processor machines has been implemented. Functionality for the fast energy routines has been increased. The vector/parallel routines will now to no electrostatics and novdw as well as simple cut-offs. ------------------------------------------------------------------------------ SPECIfy command. Controls various options such as I/O buffer flushing maximum number of processors to be used and whether to use the fast nonbond list generator. ------------------------------------------------------------------------------ SYSTem "unix bourne shell commands" This command permits the user to issue Unix shell commands from the program. The command string must be enclosed in double quotes to prevent the CHARMm parser from converting the string to uppercase. ------------------------------------------------------------------------------ SHAKE FAST This command specifies the use of the new vector/parallel SHAKE ------------------------------------------------------------------------------ Deleted Features: The old VAX analysis facility has been removed. Sigma van der Waal switching and shifting options has been removed. BARRI command removed. ------------------------------------------------------------------------------  File: ChangeLog, Node: C20-C22, Up: Top, Previous: C21-C22, Next: C22-C23 Major Enhancements and Developments in CHARMM22 As CHARMM20 is not clearly defined, it is not straightforward to sort out major differences between the current version of CHARMM (CHARMM22.0) and a previous version (CHARMM20 or CHARMm21). The VAX version CHARMM on HUCHE1 turns out to be a "developmental" version towards CHARMM21 and contains the crystal facility, BLOCK, etc. The following is prepared by comparing the developmental VAX version CHARMM21 source code and that of CHARMM22.0. Obsolete Modules Deleted from CHARMM20 -------------------------------------- [1] GRAMPS It is supported only in the VAX version CHARMM20. TH:[MK.PROT.SOURCE.VAX]GRAMPS.FLX contains an interactive routine that writes several files for the command language interpreter for producing computer graphics on the Evans & Sutherland Multi-Picture-System called GRAMPS. This obsolete feature is no longer supported in CHARMM22. [2] PARAmeter Optimization PARMOP is not incorporated in the VAX version CHARMM20 either except at the point of command parsing. It seems that the feature has never been included in the central version. New Features in CHARMM22 ------------------------ [1] BLOCK The developmental CHARMM21 VAX version supports some BLOCK commands. The BLOCK commands are used to partition the molecular system into blocks and allows for the use of coefficients that scale the interaction energies between the blocks. Specific commands to carry out free energy simulations with a component analysis scheme have been implemented. [2] CRYStal The CRYStal commands are used to build a crystal with any space group symmetry, to optimize its lattice parameters and molecular coordinates and to carry out a vibrational analysis. The CRYSTAL program is incorporated into the IMAGE module. The VAX developmental version has a separate CRYSTL module. [3] COOR COVAri The new COORdinate subcommand COVAriance is added. It computes covariances of the spatial atom displacements of a dynamics trajectory for selected pairs of atoms. [4] CORR HELIx / CORR PUCKer The New CORRelation commands HELIx and PUCKer are introduced. The HELIx command computes time series of the helical axis orientation and PUCKer computes that of the sugar pucker phase and amplitude. [5] DRAW, GRAP The new module GRAPHICS provides CHARMM the capability of displaying molecular structures when run on a graphics workstation. (Currently works only on Apollo machines.) [6] HBTRim The HBTRim command deletes hydrogen bonds that have an energy of interaction that is higher than the specified cutoff. This command is used to reduce a list of hydrogen bonds to that of important hydrogen bonds. [7] MOLVIB MOLVIB is a general purpose vibrational analysis program, suitable for small to medium sized molecules (less than 50 atoms). It performs canonic force field calculations (KANO), crystal normal mode analysis for k=0 (CRYS) and other vibrational analyses in internal coordinates or in Cartesian coordinates. Details are documented in molvib.doc. [8] PERT The PERTurbe command allows the scaling between PSFs for use in energy analysis, comparisons, slow growth free energy simulations, and widowing free energy simulations. This is a rather flexible implementation of free energy perturbation that allows connectivity to change. Also, three energy restraint terms (harmonic, dihedral and NOE) are subject to change which allows a flexible way in which to compute free energy differences between different conformations. [9] QUANTUM Quantum mechanical and molecular mechanical combined force field method is implemented by employing the semi-empirical SCF method of the MOPAC program. This module has not been tested nor documented. The code does not confirm CHARMM coding standards. The future of the code is not certain at the time of the current release. [10] RMSD The new RMSDyn routine is a modified CORMAN routine by William D. Laidig, which computes the RMS difference between two trajectory files and make a matrix of results. [11] RXNCOR The RXNCor command is used for defining a reaction coordinate for any molecule based on its structure and impose an umbrella potential along that reaction coordinate (i.e., to run activated dynamics along this coordinate) in order to trace out the free energy profile during the structural change along the coordinate. [12] SOLANA The solvent analysis facility computes solvent averaged properties, e.g., the solvent velocity autocorrelation function, mean-square displacement function, solvent-solvent radial distribution functions, solvent-reference site radial distribution function, and the solvent - reference site deformable boundary force. [13] TRAJ The new TRAJectory command is used to merges or to break up a dynamics coordinate or velocity trajectory into different numbers of units. [14] TSM The Thermodynamics Simulation Method module performs the free energy simulation. [15] Urey-Bradley Energy Term Urey-Bradley 1-3 terms have been added. The developmental CHARMM21 also includes U-B terms. [16] Update Two new non-bonded neighbour list updating schemes are introduced; one has something to do with an automated updating procedure and the other with the list generation algorithm. When INBFRQ is set to -1 (which is the default), heuristic testing is performed every time ENERGY is called and a list update is done if necessary. A new routine NBNDGC (nbndgc.src), a modification of NBONDG, is introduced. NBNDGC is based on a cubical grid searching algorithm and generates the nonbonded list in linear time, as opposed to quadratic. On the Convex C220, which is a vector machine, it is faster than NBONDG for any system larger than a few hundred atoms. [17] Integrator The leap-frog integrator has been implemented. While the "old" Verlet integrator is still available via the DYNA VERLet command (and is the default), the new integrator can be accessed by DYNA LEAP. The velocity Verlet integrator is also added in CHARMM. This new velocity Verlet integrator can be called by DYNA VVER. [18] Constant Pressure & Temperature Dynamics (DYNCPT) The constant pressure/temperature dynamics algorithm is implemented following the paper by Berendsen et al. (J. Chem. Phys. (1984) 81(8) p.3684). Modification of CHARMM20 to CHARMM22 ------------------------------------ [1] ANALysis The VAX version analysis facility is replaced by an energy contribution array (ECONT). All evaluated energy terms are partitioned into each atomic contribution and collected in the array, which is accessible through the SCALAR command. [2] XRAY The XRAY command of CHARMM20 is replaced by the READ XRAY command in CHARMM22. In CHARMM22, all I/O functions are parsed in mainio.src. The subroutine XRAY is changed to RDXRAY, which generates a card file compatible with Richard Feldmann's XRAY display program. [3] NOE NOE constraint has been overhauled. It now handles general distance restraint terms. [4] MISCOM The miscellaneous command parser (miscom.src in CHARMM22) is modified. (1) The SKIPE command is parsed in MISCOM. (2) New command parameter (@x) handling commands are added: DIVIde, EXPOnentiate, GET, MULTiply and SHOW. (3) The RANDOM command is added to set random number specifications. (4) The STOP command is parsed in MISCOM. (5) The QUICk (or Q) command is added to carry out a quick coordinate analysis. [5] HANDLE The subroutine HANDLE is improved to accept command line arguments given with the CHARMM command issued to an operating system. It works on most UNIX, UNICOS and VAX/VMS versions. [6] Command Parameters In CHARMM20, we have ten command parameters @n, where n is a single digit, 0 through 9. It is expanded to support any single alpha-numeric character so that one can use upto 36 command parameters (0-9, a-z). [7] Dynamic Memory Allocation Most of UNIX versions now support VEHEAP. VEHEAP was originally implemented by employing VAX/VMS system calls. It expands the HEAP common block when more HEAP space is needed. In UNIX versions, we use the UNIX system library routine malloc(), if available (the availability depends on the machine), to perform the same function. [8] File Format / Compatibility All binary files except dynamics trajectory are written in double precision format and not compatible with old versions. For PSF, topology, parameter, etc. one should use CARD format to transfer previous version files to CHARMM22. Trajectory files are written in single precision and compatible with all CHARMM versions and QUANTA. Old version dynamics restart files are not compatible with CHARMM22. [9] Random Number Generator All random number routines are implemented in double precision (64-bit words). Box-Muller algorithm is used for generating a Gaussian random deviat. A machine specific random number routine (RANV of CONVEX VECLIB) is used in a CMU version.  File: ChangeLog, Node: C22-C23, Up: Top, Previous: C20-C22, Next: C23-C24 Major Enhancements and Developments in CHARMM23 As an on-going project, CHARMM development has been carried out with CHARMM version 23 series. CHARMM development entails two objectives. First, we maintain an integrated macromolecular science package running on a wide range of computing devices. Second, we incorporate and exploit molecular simulation methodologies at the frontier of current research. In order to establish the first objective, we maintain all source and support files under CVS (Concurrent Versions System) control. The ROOT repository is tammy.harvard.edu:/prog/chmgr/CVS. CHARMM23 is stored in /prog/chmgr/CVS/c23a. A particular version is retrieved with the version name as the rivision tag (e.g., c23f3). Since we branched out from the CHARMM22 release version c22g2, we have made two alpha versions and four FORTRAN versions. c23a1 Developmental August 15, 1992 c23a2 Developmental October 25, 1992 c23f Developmental March 1, 1993 c23f1 Developmental March 15, 1993 c23f2 Developmental August 15, 1993 c23f3 Release February 1, 1994 c23f3 is the current release version. As the "f" in c23f stands for FORTRAN version, we converted FLECS source into FORTRAN. The conversion task had been completed as of c23f2. Now CHARMM is written in full FORTRAN except several machine dependant codes written in C. The universal languages (C and FORTRAN) make it easier to port to new machines in a broad range of architectural designs and to incorporate new methodologies into a research version of CHARMM. During the c23 development cycle, we have added and tested several new features as described below. We have also ported c23 to new machines and supported c23f versions on the following platforms. Platforms Supported ------------------- RREFX key Platforms ----------- ------------------------------------- ALLIANT Alliant ALPHA DEC alpha workstation APOLLO HP-Apollo, both AEGIS and UNIX ARDENT Stardent CONVEX Convex Computer CRAY Cray Research Inc. DEC DEC ULTRIX HPUX Hewlett-Packard series 700 IBM IBM-3090 running AIX IBMMVS IBM's MVS platform IBMRS IBM RS/6000 IBMVM IBM's VM platform IRIS Silicon Graphics MACINTOSH Apple Macintosh computers (system 7) SUN Sun Microsystems VAX Digital Equipment Corp. VAX VMS New Features in CHARMM23 ------------------------ [1] Cray Fast Code - Douglas J. Tobias Vector/parallel code for energy calculation, shake, and nonbonded list generation on the Cray was implemented. Dynamic heap and stack allocation on the Cray was added. [2] PARALLEL - Bernard R. Brooks General code for support of CHARMM on MIMD machines is completed. This includes control of the I/O levels for all file I/O. For parallel machines or workstation clusters, only node zero performs I/O and it broadcasts are to other nodes. All compuationally intensive code exercised in MD is now fully parallel which includes: DYNAMC, ENERGY (and most subsections), SHAKE, PRSSRE, DYNLNG, IMAGES,... Almost all comutationally intensive code in the first order minimizers is fully parallel. Other usage of the energy routines are parallel (such as the energy time series in CORREL). [3] Dynamics Integrator 3.1 Leap-Frog Integrator - Bernard R. Brooks Berendsen's method was modified so that it would work for very small systems and for very weak coupling constants. Now it is possible to use SHAKE with CPT and get correct pressures and temperatures. Another change is to calculate the change in potential energy due to the constant pressure algorithm. The energy lost due to the changes in box size is now added to the kinetic energy during the constant temperature procedure. This allows the constant presure code to nearly conserve energy and allows the constant temperature code to be used with weak coupling times. This correction was made when we found that water box simulations with the Berendsen's method were running about 10 degrees too cold when both temperature and pressure coupling times of 1ps were used. Now the correct target temperature is achieved, even in the limit of very weak couplings. 3.2 EULER Dynamics Integrator - Bernard R. Brooks The incorporation of of the Langevin/Implicit Euler dynamics integrator has been achieved. The effect is to remove the energy in the high frequency degrees of freedom which eliminates the noise in free energy studies where bonds are being modified. To support the Implicit Euler integration, a Truncated Newton Minimizer has been added. This minimizer may be used directly using the MINI TN command. The minimizer is not yet fully implemented (it works, but is not as efficient as it will be), but it is already very competitive relative to existing minimization methods. MINI TN does not work with SHAKE. This code has been developed by Tamar Schlick at NYU. It has been integrated within CHARMM with some modifications. 3.3 EHFC: High Freequency Correction - Bernard R. Brooks The leap-frog dynamics integrator has been modified to have an improved high frequency correction (HFC) term. With the old term, energy was conserved within a harmonic degree of freedom, but total energy would drift as energy exchanged between high and low frequency degrees of freedom. The new code avoids this problem. The total energy and kinetic energy that is printed in the first line of dynamics energy printout has reverted to the standard Verlet energies, and these match the output of the old integrator. The HFC terms (total energy, and kinetic energy) are now printed on the second line. The fluctuation of the HFC total energy is usually an order of magnitude smaller than that of the total energy. The HCF total energy is a good indicator of problems with NVE dynamics because small changes in total energy are not lost in the noise of high frequency oscillations. 3.4 Velocity Verlet Integrator - Masa Watanabe Velocity Verlet method has been implemented. Two integrator (Verlet and Leap-frog) methods presented in CHARMM have their own flavors, but Verlet method handles velocities rather awkward and may introduce some numerical imprecision. On the other hand, the Leap-frog integrator minimizes loss of precision on a computer, but it does not handle the velocities in a satisfactory manner. Velocity Verlet integrator can store positions, velocities, and accelerations all at the same time and minimizes round-off error. 3.5 Nose-Hoover Constant Temperature Method - Masa Watanabe The constant temperature method has been implemented based on S. Nose, JCP 81, 511 (1984) and W.G. Hoover, Phy. Rev. A 31, 1695 (1985). This is an another type of constant temperature method, but an equilibration time in the vicinity of the desired temperature is faster than other routines which are available in CHARMM. Also multi-temperature controls are also developed in order to equilibrate the system faster and keep the system in the desired temperature well. This method works with Verlet and Velocity Verlet integrators. 3.6 Multiple Time-Scaled Method - Masa Watanabe Tuckerman et al proposed a reversible RESPA algorithm recently (Tuckerman, Berne, Martyna, JCP 97, 1990 (1992)). Previous MTS methods have the disadvantages of loosing accuracy due to the approximation of holding the slow variables fixed while integrating the equations for the fast variables. But in this reversible RESPA equations of motions are derived from Liouville operators and Trotter theorem. The method gives more accurate dynamics than previous methods. In this implementation, one can specify up to three different time steps in dynamic simulation run. [4] RISM (Reference Interaction Site Model) - Georgios Archontis The RISM module allows the user to calculate the site-site radial distribution functions g(r) and pair correlation functions c(r) for a multi-component molecular liquid. These functions can then be used to determine quantities such as the potential of mean force or the cavity interaction term between two solute molecules into a solvent, and the excess chemical potential of solvation of a solute into a solvent. The change in the solvent g(r) upon solvation can be determined and this allows for the decomposition of the excess chemical potential into the energy and entropy of solvation. [5] MMFP (Miscellaneous Mean Field Potential) - Benoit Roux The MMFP Commands are primarily used for setting up special restraining potentials on some or all of the atoms. The key word MMFP is used to enter the MMFP environement. In the MMFP environment, all miscelaneous commands (label, goto, if, etc...), and string substitutions (with @1, @2, etc...) are supported. The key word END returns to the main parser. The restraining potentials are used in all energy calculations, unless SKIP is used. The subcommand RESET clears the potential. This module is still under developement and only the subcommand GEO is released. The subcommand GEO (standing for geometrical) is used to setup various restraining potential (spherical, planar or cyclindrical restraints) on some or all atoms. The selection specification should be at the end of the command. The default atom selection includes all atoms. Future subcommands will include continuum electrostatic reaction field and solvent mean field potentials. Expected date of release is Spring 1994. [6] NMR Analysis - Benoit Roux The NMR commands may be used to obtain a set of time series for a number of NMR properties from a trajectory. Among the possible properties are relaxation rates due to dipole-dipole fluctuations (T1, T2, NOE, ROE), chemical shift anisotropy and Deuterium order parameters for oriented samples. [7] REPLICA - Leo Caves Tool to support LES and MCSS calculations. Performs replication of arbitrary regions of PSF. Data structure interfaces to non-bond list generation routines, to perform appropriate exclusions. In association with BLOCK can provide appropriate energy/force normalizations for various classes of methods employing replicas. Introduced REPLICA and REPDEB preprocessor directives. Code for cray multi-tasking list generation routine used inference and has not been tested. Convex parallel code works fine. Added miscellaneous parameters to report number of atom/group pairs from non-bonded routines: ?NNBA, ?NNBG, ?NNBI for atom/group/images respectively. For replica-based exclusions from the list there are ?NRXA and ?NRXG for atom and group exclusions. [8] Clustr code integrated into CORREL - Charles L. Brooks III The CLUSTER command clusters time series data obtained within the CORREL facility. The data are grouped into sets with similar time series values, using euclidean distance as the dissimilarity measure between different time frames of a set of time series. It is useful, for example, for grouping together similar conformations or energy levels. [9] GRAPHICS - Richard M. Venable Graphics code converted to FORTRAN and overhauled. Versions that work with Xwindows and GL are in progress. A new preflx keyword, NODISPLAY, builds a version which produces HPGL, PLUTO FDAT, and LIGHT.atm files without requiring any screen display capabilities. The SG (IRIS) code incorporation is relatively untested. Postscript file output similar to HPGL (but much nicer looking, hopefully) is also implemented. Major Modifications ------------------- [1] Command Line Handling 1.1 Extension of Command Line Parameter Handling - Leo Caves A command line parameter token can now be a string rather than just one of the single characters 0-9 and A(a)-Z(z). For substitution, a token is indicated by the use of the @ character as before. The token is end-delimited by any non-alphanumeric character. In the case that the token is not found in the parameter table, a check is made to see if the first character of the token is itself a token in the parameter table. If this single character token is in the table, the corresponding value is substituted -- this is the necessary scheme to allow backwards compatibilty with the old parameter substitution, which allowed parameters embedded in strings. For unambiguous token detection, "protect" the token with brackets {} --- this allows for the use of non alphanumerics in tokens such as -, _. 1.2 New Parsing Options - Bernard R. Brooks The IF command will be expanded to allow commands such as: IF ?ENER .GT. ?VDW THEN GOTO label or IF ?NSEL .LT. 8 THEN GOTO label 1.3 MSCNUM - Bernard R. Brooks New code for flexible miscellaneous command substitutions has been fully incoporated. Additional types were needed to make this code more flexible. Three types are supported, REAL(*8), INTEGER, CHARACTER. There are three subroutines which can be called; integer (SETMSI), character (SETMSC), and real (SETMSR) to specify a command substitution variable. Now it is possible for ?NATOM to return an integer, ?RSM to return a real number, and ?SEGID to return the segment identifier of the first selected atom. [2] QUANTUM Quantum mechanical and molecular mechanical combined force field method was implemented by employing the semi-empirical SCF method of the MOPAC program in the CHARMM version 22. The QUANTUM code has been modified extensively to meet CHARMM standards. There were several problems with the quantum code that have been fixed. The van der Waal group nonbond list was missing due to an improper interpretation of the group-group exclusion list in CHARMM (It's a two state list, not a 3 state as in the atom-atom exclusion list). All vdw interactions between QM and MM group where any QM atom had an exclusion or a 1-4 interaction with any MM atom were not computed. This caused major problems in certain situations where there was a strong electrostatic attraction with no compensating vdw interaction. New code to add link and place link atoms has been written. [3] Frequency Based Crystal Update - Ryszard Czerminski The modification allowes for automated, frequency based, crystal update. New variable (IXTFRQ) is introduced which controls frequency of the crystal update. [4] Ability to Linearly Increase/Decrease Pressure - Ryszard Czerminski The goal was to allow for linear increase (decrease) of the pressure during single dynamic run. New variables/keywords were introduced (PIXX - initial value of XX component of pressure tensor, PFXX - final value etc... for other components). [5] Atom Selection 5.1 Atom Parse - Bernard R. Brooks A new atom name parsing subroutine has been developed. This makes the code simpler and facilitates further advancements in atom parsing. One new feature allows an atom selection to be used to select a series of atoms. This is very useful in CORREL for specifying clusters of atoms for analysis. When the atom selection feature is used to specify 4 atoms of a dihedral, the first 4 selected atoms will be chosen. 5.2 New Tokens - Bernard R. Brooks new operator; .BYGROUP. new token; IGROup : have been added to allow the selection of atoms based on electrostatic groupings. Several keynames have been added to allow the query of the characterstics of selected atoms; ?SELATOM - number of first atom selected ?SELIRES - number of first residue selected ?SELISEG - number of first segment selected ?SELTYPE - name of first atom selected ?SELRESI - resid of first residue selected ?SELSEGI - segid of first residue selected ?SELRESN - residue type of first atom selected ?SELCHEM - chemical type of first atom selected These new keywords are in addition to the existing keyword; ?NSEL - Number of atoms selected [6] Correlation 6.1 New MANTim Options in CORREL - Bernard R. Brooks A histogram option to time series manipulation has been developed. This is executed by the command; MANTime time-series-name HISTogram min-value max-value num-steps The selected time series is replaced with a histogram which contains the probability of finding the time series within a given value range. Also, new options (RATIo and KMULt) added to the CORREL MANTIME command. 6.2 Dihedral Time Series in CORREL. - Bernard R. Brooks Fixed problems with the diheral code in correl to account for torsional timeseries. The correct fluctuation is now determined. The extra processing has been removed from the SHOW command because the data may no longer be valid for this processing when MANTIME commands are present in a script. A new command option "MANTime CONTinuous-dihedral" has been added to allow a dihedral timeseries to be unfolded to a continuous function. 6.3 Extension of Solanal ANALysis command - Arnaud Blondel A command -CROSs- was added to allow a cross analysis on two selected subsets of atoms. For the moment the exclusion of the couple of atoms belonging to the same SEGId is not implemented. The keyword CROSs cannot be selected with the following options: WATer, SITE, IKIRkg, ISDIst, IFDBf. IVAC, IMSD and IFMIn have not been tested with CROSs. [7] SCALAR Command Enhancement - Bernard R. Brooks The ASP arrays (IGNOre, ASPV and VDWS) are now accessible. There is a sort option for the SHOW command. There is a new MASS keyword for the STATistics and AVERage commands A new SCALAR READ option has been added. It allows values to be entered from a file. The use is: OPEN READ CARD UNIT 12 NAME file.dat SCALar WMAIn READ 12 SELE ... END which will read selected entries to the weighting array. [8] SURFACE - Bernard R. Brooks New analytic surface area code and energy terms for ASP (Atomic Solvation Parameters) energy and forces have been fully integrated (and parallelized for multi-machines). This has been achieved by the incorporation and adaptation of the code from Wesson and Eisenberg. The default for the COOR SURFace command is now the analytic surface area. The anaylitic answer is less expensive and more accurate. The older Lee and Richard's algorithm may still be invoked by specifying a nonzero RPRObe value. The maximum number of contacts that a sphere may have has been increased from 15 to 35. [9] QAUGMENT - Bernard R. Brooks It is desirable for a patch to be able to augment the charge of an atom. The current code could only set a charge. The new code can add or subtract a value from the charge. This is done by using a patch charge value near 100.0. For example, a charge of 100.15 will add 0.15 to the current charge. A charge value of -101.0 will subtract 1.0 from the current charge. Charge values less than -90.0 or larger than 90.0 are no longer allowed for generate or patch without charge augment. It allows more flexible patches to be developed where the prior charge on modified atoms need not be known. [10] COORdinate Commands 10.1 VACUUM_OP: COOR SEARCH Subcommand - Bernard R. Brooks The ability to manipulate pixel bitmaps generated from the COOR SEARCH command has been developed. The new syntax for the COOR SEARCH command is; COOR SEARch {PRINt [UNIT int]} { } {[VACUum]} {[RESEt]} [SAVE] {[NOPRint] } {[RCUT real]} { FILLed } { AND } {[RBUFf real]} { HOLES } { OR } { XOR } The new keywords are; SAVE - save the resultant bitmap for subsequent operations AND - logical AND the new bitmap with the previously saved map OR - logical OR the new bitmap with the previously saved map XOR - logical XOR the new bitmap with the previously saved map HOLES - search for holes (vacuum points surrounded by filled points) 10.2 New COOR DIST command - Bernard R. Brooks The COOR DISTance command has been overhauled and has additional features. One such feature is the ability to get g(r) plots from trajectory files using atom selections. It has several other features. The new syntax is: COOR DISTance { WEIGhting vector-spec atom-selection } { } { [UNIT int] [CUT real] [ENERGy [CLOSe]] 2X(atom-selection) - } { [Nonbonds] } { [NO14exclusions] } { [NOEXclusions] } - { NONOnbonds } { 14EXclusions } { EXCLusions } [TRIAngle] [ HISTogram HMIN real HMAX real HNUM integer - [HSAVe] [HPRInt] [HNORm real] [HDENsity real] ] [11] JOIN/RENUMBER Command - Bernard R. Brooks A "JOIN segid RENUMBER" feature is added in the JOIN command. This allows resid's to be made sequential within a single segment. [12] PREFX.SRC overhauled. - Bernard R. Brooks The PREFX program has been overhauled. The new code has the following features: - It allows "!" comments at the end of valid FORTRAN statements. - Conversion to single precision is performed ONLY if the SINGLE keyword is present. - It allows the use of identifier comments in ## statements. For example: ##IF PERT (pertprint) ... ##ELSE (pertprint) ... ##ENDIF (pertprint) This makes the code easier to read and allows ##ENDIF statements to be uniquely identified. A fatal error is flagged if the identifiers do not match.  File: ChangeLog, Node: C23-C24, Up: Top, Previous: C22-C23, Next: Top Major Enhancements and Developments in CHARMM24 During the C24 development cycle, February 15, 1994 to February 15, 1996, we made two bugfix-updates in the c23 releases and three alpha versions and one beta version in the c24 development line. c24x1 is the MMFF implementation in CHARMM developed at the Molecular Simulations Inc. CHARMM23.0 c23f4 Release August 15, 1994 c23f5 Release March 15, 1995 CHARMM24.0 c24a1 Developmental February 15, 1994 c24x1 Evaluation February 15, 1994 c24a2 Developmental August 15, 1994 c24a3 Developmental March 15, 1995 c24b1 Release August 15, 1995 Only bugfixes are incorporated into CHARMM23 and all new developments and enhancements have been carried out with the CHARMM24 developmental versions. All modifications are thoroughly recorded in the ChangeLog.c24 file and the following is the summary of new features and major enhancements in CHARMM 24. New Features in CHARMM24 ------------------------ [1] New Ports and Parallel Versions 1.1 Enhancement to Parallel Code - Bernard R. Brooks and Milan Hodoscek There has been continued development of the parallel code for CHARMM. This includes new features run in parallel, new machine types supported, new parallelization methods, and code made to run more efficiently. Due to conflict in routine names with library routines, the subroutines: WRITEC and READC had to be renamed. Initial code to allow the use of the Terra parallel computer has been added. Added preflx keyword SGIMP for multiprocessor SG machines using PVM massage passing library. The difference between PVM and (SGIMP, PVM) is that all the processes are spawned on one host and some communication parameters are not supported on MP machines. It can be used on a single processor SG for testing purpose. Use PVM only on a cluster of any type of workstation. 1.2 Convex Exemplar SPP-100 and generic PVM Ports - Charles L. Brooks, III and Stephen H. Fleischman A port of CHARMM version 24a2 to general PVM based parallelism using existing parallel code as well as a port to the Convex parallel machine are included. 1.3 Cray T3D Port - Charles L. Brooks, III and Barry C. Bolding A port of CHARMM version 24a2 to the Cray T3D parallel computer using existing parallel code is included. 1.4 Port of parallel CHARMM to Convex Exemplar SPP-1000 and generic MPI - Charles L. Brooks, III and Stephen H. Fleischman A port of CHARMM version 24a3 to general MPI based parallelism using existing parallel code as well as a port to the Convex parallel machine are included. 1.5 Thinking Machine's CM5 Port - Robert Nagle Previous communication scheme was based on a simple send and receive model. By using TMC's active message layer, communication bandwith can be increased by anywhere from 50% to 5X. 1.6 OS/2 Port - Stefan Boresch CHARMM (c23f4 and c24a3) has been ported to the OS/2 operating system, version 2.x and higher. The Watcom Fortran compiler (v. 9.5, patch-level (c)) has been used. A new pre-processor keyword, OS2, has been introduced, and all OS/2 related changes hide behind the OS2 keyword. There is currently no install script. Please contact me if you want to build an OS/2 version of CHARMM (boresch@tammy.harvard.edu). [2] Fast Multipole Code for Electrostatic interactions - Robert Nagle This is an initial implementation of a fast multipole method, based on John Board's work. A new non-bond option (FMA) has been added. This replaces cut-off parameters with a no cut-off hierarchical technique. The advantages of this method are that you can control the error and that it is amenable to parallelization. FMA is an O(N) technique but the constant is large and so FMA will, in general, be slower for systems of less that 5000 atoms, for the same accuracy. Two options, LEVEL and TERMS, govern how many hierarchical levels are used and how many terms are retained in the expansion, respectively. In the method, each box at every level is subdivided into 8 sub-boxes - you should select LEVEL so that the boxes at the lowest (i.e. finest) level contain 10-20 atoms on average: 3 or 4 will be typical choices. You then select TERMS to control the accuracy that you require: 4 will often suffice but I would generally recommend 6 or even 8. See the references in fma.doc for a detailed description of the error bounds. NOFMA is the nonbond option which turns off the multipole method. Compilation of FMA is controlled by the flag, FMA, in pref.dat. FAST ON is required for this initial implementation. This implementation is not yet parallelized. [3] Energy Embedding by the Addition of a Higher Spatial Dimension - Elan Z. Eisenmesser / Carol Post The energy embedding technique entails placing a molecule into a higher spatial dimension [Crippen, G. M. & Havel, T. F. (1990) J. Chem. Inf. Comput. Sci. Vol 30, 222-227]. The possibility of surmounting energy barriers with these added degrees of freedom may lead to lower energy minima. With the recent success of using four dimensions in the GROMOS force field [Van Schaik, R. C., Berendsen, H. J. C., Torda, A. E., & van Gunsteren, W. F. (1993) J. Mol. Biol. Vol 234, 751-762], creating a similar option in CHARMM should also prove advantageous. Specifically, another cartesian coordinate was added to the usual X, Y, and Z coordinates and was appropriately named FDIM for Fourth DIMension. This implementation has led to alterations in some existing code along with the addition of several algorithms. [4] DIMB (Diagonalization In a Mixed Basis) Method - David Perahia, Liliane Mouawad, Herman van Vlijmen The DIMB (Diagonalization In a Mixed Basis) method (see L. Mouawad and D. Perahia (1993), Biopolymers, 33, 599) is an iterative method to calculate the N lowest normal modes of molecules. It is especially targeted to do large molecules, since it does not require the full Hessian to be stored in memory or on disk. In short, the method does repetitive reduced-basis diagonalizations in bases that consist partially of the approximate eigenvectors, and partially of Cartesian coordinates. Eigenvectors are saved to file during the process. Before that is done, a new basis is again created, which consists of the approximate eigenvectors at that point + the residual vectors (Lanczos vectors). This accelerates the convergence. A very good property of this method is that the final eigenvectors are as accurate as the user wants them to be, so the results are no different from a full-blown diagonalization. Because the method is iterative, it takes longer to converge than a regular diagonalization. Sizewise it can handle almost anything on a moderately sized computer. David Perahia calculated a few dozen modes of Hemoglobin (~600 residues = ~6000 atoms = ~18000 d.o.f.) on a SGI workstation with 90 Mb memory. I have done several calculations on 900 residue systems. The actual time to reach convergence depends on the available memory, the desired accuracy, and the number of requested normal modes. One other area where the method saves memory is in the storage of the original Hessian. Since this matrix is usually sparse for large systems, a compressed Hessian is set up, which contains all non-zero elements. In addition, I added the option to used this compressed Hessian in the reduced-basis diagonalization option of VIBRAN. Before, the same size limits applied to full diagonalizations and reduced-basis diagonalizations. This should not be: people usually want to do reduced-basis calculations because the molecule is too big for the Hessian to be stored in memory. The option VIBRAn REDUce CMPAct will fill the compact Hessian and form the reduced-basis Hessian from this compact Hessian. Overall, this is a big saving on memory space. [5] Arithmetic Expression Interpreter - Benoit Roux An interpretor of arithmetic expression has been added to the CHARMM command parser. It is called at the level of the miscellaneous command handling using simply by the word CALC (for calculator). It can be used to evaluate algebraic numerical expression. The command supports all mathematical numerical expression with arbitrary number of nesting of recursive parentheses, e.g., exp[1.0-cos(2*(log(2*pi))**2)/0.5] The parsing is actually very crude since the expression is translated back and forth between character string and a real variable to handle the logic (there is no real subroutine recursion). [6] TNPACK Update - Tamar Schlick, Phillipe Derreumaux and Eric Barth The truncated-Newton minimization package TNPACK, developed by T. Schlick and A. Fogelson, has been incorporated into CHARMM and adopted for biomolecular energy minimization. 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 a paper now in press in the Journal of Computational Chemistry: Applications are reported for a series of molecular systems including Alanine Dipeptide (N-Methyl-Alanyl-Acetamide), a dimer of N-Methyl-Acetamide, Deca-Alanine, Mellitin (26 residues), Avian Pancreatic Polypeptide (36 residues), Rubredoxin (52 residues), Bovine Pancreatic Trypsin Inhibitor (58 residues), a dimer of Insulin (99 residues), and Lysozyme (130 residues). Through comparisons among the minimization algorithms available in CHARMM, we find that TNPACK performs significantly better than ABNR in terms of CPU time when curvature information is calculated by a finite-difference of gradients (the "numeric" option of TNPACK). The CPU gain is 50% or more (speedup factors of 1.5 to 2.5) for the largest molecular systems tested and even greater for smaller systems (CPU factors of 1 to 4 for small systems and 1 to 5 for medium systems). With the analytic option, TNPACK converges 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. TNPACK uses curvature information to escape from undesired configurational regions and to ensure the identification of true local minima. It converges rapidly once a convex region is reached and achieves very low final gradient norms, such as of order 10E-8, with little additional work. Even greater overall CPU gains are expected for large-scale minimization problems by making the architectures of CHARMM and TNPACK more compatible with respect to the second-derivative calculations. This work should be the focus of future developments. Such work involves sparse storage of the Hessian, efficient sparse Hessian/vector multiplications, and separation of the gradient and Hessian calculations. [7] X-window graphics extensively modified. - Richard M. Venable Several new features have been added to the X-window version of CHARMM graphics. This code has also been tested on a wider variety of hardware platforms (for example: SGI). Changes include: double-buffering, clipping, StaticColor, symbol fonts, window title, modified colormap calls, and a misc. Bug fixes in the labeling of the X axis. A NODISPLAY compile option has been added to the X windows version of CHARMM graphics in which only derivative files are produced. The GRAPhics NOWIndow option can be used to generate the same effect at run time. [8] Minimum Image Periodic Boundary Code - Charles L. Brooks, III, William A. Shirley and Stephen H. Fleischman Simple minimum periodic boundary conditions are added for cubic, truncated octahedra and rhomboidal (dodecahedra) periodicities which augments the image facility and enhances parallel scaling on scalar parallel machines as well as significantly reducing the memory requirements. This code is developed and fully tested for the simulation cells described above when the cell edgelength is the same in all dimensions. The (trivial) extension to non-identical cell sides will be added. However, it is critical to see reasonable performance on all scalar parallel platforms where simulations using images are currently employed that this enhancement be added now. [9] GAMESS Code - Bernard R. Brooks and Milan Hodoscek The CHARMM-GAMMES interface is under development. The interface part is completed and testing is in progress. Major Enhancements in CHARMM24 ------------------------------ [1] New Dihedral / Improper Dihedral Energy Routines - Arnaud Blondel The previous energy routines used the derivatives d(cos(phi))/dr to calculate the forces and the second derivatives. This choice introduced an artificial singularity at sin(phi)=0. The new routines use the derivative d(phi)/dr and thus have no singularities. This removes the tests to avoid numerical overflow or the switch functions in the vector improper routines. The new dihedral routines now support cases where planar conformation is not an extremum. Thus a value other than 0 or 180 can be specified in the dihedral parameters. The dihedral constraints can also use the dihedral functional form using the key word PERIod and giving a non-zero number. [2] Extended Pressure System, Langevin Piston Code - Bernard R. Brooks, Scott E. Feller and Yuhong Zhang The constant pressure code has been overhauled. The old method based on Berendsen's method has been replaced with a Langevin Piston Method. When no friction is applied, this method becomes the standard method based on Nose and Klein (adapted from Andersen). At the limit of infinite friction with no random force, this reverts to the Berendsen method. The unit cell information has been added to the trajectory file format. This implementation required an update to the image and crystal code which cleaned up some ancient problems. Options for including the surface tension (gamma-Area) term is also completed and tested. This has been developed for the accurate simulation of interfacial systems. [3] Anisotropic Harmonic Restraints - Bernard R. Brooks The global scale factors: "XSCAle", "YSCAle", and "ZSCAle" have been added to the "CONS HARM" command. This allows using the CONS HARM to enforce a planar or linear restraint. This feature is also useful for use in conjunction with our COORPLAS program (for generating 3-D coordinates from plastic models). [4] New RESDistance Facility - Bernard R. Brooks A new facility, RESD, has been created to allow general distance restraints based on a linear combination of distances. This is useful for searching reaction pathways. [5] New READ PARAm APPEnd Option - Bernard R. Brooks An append option has been added to the READ PARAM CARD command. This allows just a few parameters to be modified without editing an entire parameter file. A modification to the binary parameter file format was necessary. Old binary files may not be appended, but they are still supported. [6] New READ PSF APPEnd Option - Bernard R. Brooks An append option has been added to the READ PSF command. This allows PSFs to be easily merged to make a larger PSF. No modification to the binary parameter file format was necessary. This option works with both FILE and CARD options. [7] Best Fit Option to CORREL TRAJectory Command - Bernard R. Brooks The TRAJectory command in correl now accepts an ORIENt keyword with an optional [MASS] qualifier in conjunction with a second atom selection that will best fit selected atoms with respect to the rms deviation from the reference structure (in the comparison coordinate set). This operation is done prior to the determination of any time series value. This operation will not affect any time series value that is based only on relative distances and angles. [8] QM/MM Exclude Group Option - Bernard R. Brooks An option EXGRoup has been added which causes all atoms in the group of the link atom host to be excluded from the QM/MM electrostatic interaction terms. Code for specifying the charge of link atoms and their placement has also been added. [9] Enhancements to the Ewald Code - Bernard R. Brooks, Scott E. Feller and Steve Bogusz The EWALD electrostatic option now runs efficiently for parallel architectures. Also, the maximum K-space values can be specified independently for each direction. Several bugs were fixed. Additional ways to compute ERFC() were added, including a lookup table. [10] MMFP/SSBP Upgrade - Benoit Roux and Dmitrii Beglov The Miscellaneous Mean-Field Potentials (MMFP) has been upgraded. The spherical solvent boundary potential (SSBP) has also been incorporated into EPERT. A new "membrane-like" planar potential has been introduced using Gaussians to provide a smooth free energy function based on hydropathy profile of individual amino acids and solvent exposure. This is useful to orient membrane proteins. A new primary shell of hydration has been added to the MMFP facility to provide one layer of solvent around a flexible polypeptide. For more information, see Beglov & Roux, Biopolymers 35: 171-178 (1995). A solvent boundary potential for the simulation of water at constant pressure is also added to the Miscellaneous Mean Field Potential module. The boundary potential is an approximation but follows from a rigorous statistical mechanical treatment of the boundary. In light of the difficulties raised by the previous treatments, a different route was chosen to formulate and develop the solvent boundary potential for computer simulations of a finite representation of an infinite bulk system. The present theoretical formulation is based on a separation of the multidimensional solute-solvent configurational integral in terms of n "inner" solvent molecules nearest to an arbitrary solute, and the remaining "outer" bulk solvent molecules. This formulation, which differs significantly from previous treatments, provides further insight into the statistical mechanical basis of the solvent boundary potential and is helpful in constructing useful approximations for computer simulations in dense liquids. An approximation to the solvent boundary potential is constructed for simulations of bulk water at constant pressure, including the influence of van der Waals (done with RISM) and electrostatic interactions (done with a Kirkwood-like multipole expansion). The approach has been tested with success on several typical systems (water, ions, n-butane and alanine dipeptide). [11] Upgrade of the NMR module - Benoit Roux The NMR module is upgraded to have better output style. The old version used the value of PRNLEV to choose the printed quantities. Since this was a non-standard style in CHARMM, a series of logical flags have been included in the command calls to print some chosen quantities. In addition, the chemical shift anisotropy (CSA, used in solid state NMR of membrane proteins in oriented samples) has been redefined in term of a zmatrix to prevent confusion. The deuterium quadrupolar splittings (DQS) command is also upgraded. A bug in a call to NORMAL was fixed. [12] New Options to CORREL - Lennart Nilsson Two new MANTime options have been added to CORREL: CROS and DOTP. CROSsprod name Q(T) = Q(T) x Q2(T) produces the 3D crossproduct of the two 3D vectors formed by the selected and named timeseries and DOTProd name Q(T) = x-comp of Q(T)= Q(T) . Q2(T) gives x-comp of Q2(T) angle in degrees between the two vectors. [13] The COOR HBONd Command - Lennart Nilsson An option for the analysis of H-bond patterns from trajectories has been added to corman. COORdinates HBONd 2X(atom-selection) [CUT ] [CUTA ] [IUNIt ] [BRIDge ] [FIRSt int] [NUNIts int] [NSKIp int] [BEGIn int] [STOP int] The HBONd command analyses a trajectory for hydrogen bonding patterns. For each acceptor/donor in the first selection the average number and average lifetime of hydrogen bonds to any atom in the second selection is calculated. A hydrogen bond is assumed to exist when two candidate atoms are closer than the value specified by CUT (default 2.4A, (reasonable criterion, DeLoof et al. (1992) JACS 114, 4028), and if a value for CUTAngle is given the angle formed by D-H..A is greater than this CUTAngle (in degrees, 180 is a linear H-bond); the default is to allow all angles. The current implementation assumes that hbonding hydrogens are present in the PSF and also uses ACCEptor and DONOr information from the PSF to determine what pairs are possible. If output is wanted to a separate file the IUNIt option can be used. If the BRIDge option is used the routine calculates average number and lifetime of bridges formed between all pairs of atoms in the two selections; a bridge is counted a residue of the type specified with the BRIDge hydrogen bonds (using same criteria as for direct hbonding) to at least one atom in each selection. The typical use of this would be to find water bridges. Here again, results are presented for each atom in the first selection. In order not to find hbonds between bonded atoms UPDATE is called, which requires coordinates to be present when invoking this module. Since this is done just to get the non-bond exclusion lists, the cut-offs are set to very small values, and could influence subsequent energy evaluations if the non-bond cutoffs are not then respecified. [14] NORESET Option for SHAKE - Lennart Nilsson The NORESET option is added to allow multiple shake commands. It is useful to be able to define shake on bonds, bonh or so on several different sets of atoms, with different shake options. The NORESET keyword to shake command allows this by not zeroing counter. [15] Trajectory Reading - Lennart Nilsson READCV is modified to read coordinates at multiples of skip FROM the actual first coordinate set in a trajectory file. [16] Make BLOCK work with IMAGE/CRYSTAL and vice versa - Stefan Boresch In order to make BLOCK work / coexist with the IMAGE module two things had to be changed: (1) A memory allocation problem in the BLOCK datastructure and (2) the post-processing modules needed to be overhauled to allow for nonbonded list updates while reading frames from the trajectory. Ad (1), memory allocation: BLOCK uses two data-structures, one containing the interaction matrix between blocks, and one containing the block number for each atom (IBLCKP). This array was allocated so far as INTEG4(NATOM) on the heap. However, when IMAGE atoms are present, the energy routines attempt to find out to which block an IMAGE atom belongs. This at one point or the other causes a memory access violation. The solution consists out of two parts. (i) The IBLCKP data-structure is now allocated as INTEG4(MAXAIM) on the heap; therefore there is always enough space provided. (ii) The entries for the IMAGE atoms have to be initialized, and this has to be done at EVERY image update. However, similar things are already done for a number of other quantities like masses, vdW params, charges etc. All this is done among a number of other things in subroutine MKIMAT in upimag.src, where I have added an appropriate statement. Ad (2), changes to post-processing routines: Real/Image atoms leave/enter the simulation box/system dynamically. Therefore, the nonbonded/image interaction lists have to be updated during post-processing. The hooks were already in the program, subroutine BLUPLST. The real changes hide in this routine, most changes in BLFREE, BLEAVG and BLCOMP are either cosmetic or ensure proper printout. Post-processing routines FREE, EAVG and COMP will actually print IMAGE terms if present. The routine BLUPLST is a sibling of routine updeci in heurist.src. The heuristic update scheme itself is removed, as I feel that one should update the lists at every frame. Also, the CRYSTAL specific section of UPDECI is not present in BLUPLST as I don't understand it. Therefore, care should be exercised when using BLOCK with CRYSTAL! Negative values of INBFRQ/IMGFRQ are trapped, in this case they are set to 1; Printout from the update / list generation routines is suppressed by temporarily raising the PRNLEV to 1. The BLOCK documentation (block.doc) has been revised and reflects these modifications. A new testcase block3.inp has been added to test/c24test. [17] Constraint correction for PERT - Stefan Boresch The current version of PERT cannot handle situations where SHAKE is applied to bonds which change in length due to an alchemical mutation as SHAKE and PERT do not "communicate". Furthermore, in such cases a constraint correction has to be computed and added to the free energy difference. Two steps are required to fix this problem: (1) The constraint list needs