CHARMM 27 B1 Documentation files CHARMM Element doc/ace.doc $Revision: 1.1 $  File: ACE, Node: Top, Up: (chmdoc/commands.doc), Next: Syntax Analytical Continuum Solvent (ACS) Potential Purpose: calculate solvation free energy and forces based on a continuum description of the solvent, in particular the analytical continuum electrostatics (ACE) potential. Please report problems to Michael Schaefer at schaefer@brel.u-strasbg.fr WARNING: The module is still being developed and may change in future versions. REFERENCES: M. Schaefer & M. Karplus (1996) J. Phys. Chem. 100, 1578-1599. M. Schaefer, C. Bartels & M. Karplus (1998) J. Mol. Biol., in press. * Menu: * Syntax:: Syntax of the ACE specifications * Defaults:: Defaults and Recommended values * Function:: Purpose of each of the specifications * Examples:: Usage examples of the ACE module  File: ACE, Node: Syntax, Up: Top, Previous: Top, Next: Defaults Syntax [SYNTAX ACE functions] Syntax: The ACE specifications can be specified any time the nbond specification parser is invoked, e.g., ENERgy [other-spec] [ace-spec] ace-spec::= [ ACE ] [ IEPS real ] [ SEPS real ] [ ALPHa real ] [ SIGMa real ]  File: ACE, Node: Defaults, Up: Top, Next: Function, Previous: Syntax The defaults for the ACE potential are IEPS 1.0 SEPS 80.0 ALPHa 1.2 SIGMa 3.0 In the current implementation, ACE should be used with united atom parameters, ALPHa set equal to 1.2 and the PARAM19 parameter file param19-1.2.inp; the param19-1.2.inp file inludes a table of atom volumes at the end which is compatible with ALPHa 1.2.  File: ACE, Node: Function, Up: Top, Previous: Syntax, Next: Examples 0. Introduction The analytical continuum solvent (ACS) potential is introduced to perform molecular dynamics/minimization calculations with a continuum approximation of the solvent. Two solvent contributions to the effective (free) energy of a solute are included: the electrostatic solvation free energy, and the non-polar (i.e., non-electrostatic) solvation free energy. The first (electrostatic) contribution (G_el) is calculated using an analytical approximation to the solution of the Poisson-equation called ACE (from: analytical continuum electrostatics). The non-polar solvation free energy (G_np) is approximated by a pairwise potential which yields results that are very similar to the well-known surface area approximations of the hydrophobic (solvation) energy (e.g., Wesson and Eisenberg, Prot. Sci. 1 (1992), 227--235; see the ASP potential in CHARMM). Restriction: The ACE solvation potential has to be used together with no cutoff or with atom based switching. Compatibility: ACE can be used with BLOCK (but: the diagonal elements of the BLOCK matrix MUST NOT be zero). Meaning of the ACE parameters: 1. IEPS Dielectric constant for the space occupied by the atoms that are treated explicitly, e.g., the space occupied by the protein. 2. SEPS Dielectric constant for the space occupied by the solvent that is treated as a continuum (i.e., the complement of the space occupied by the protein). 3. ALPHa The volumes occupied by individual (protein) atoms are described by Gaussian density distributions. The factor ALPHa controls the width of these Gaussians. The net volume of the individual atom Gaussian distributions is defined in the volume table at the end of the parameter file param19-1.2.inp. The width of the atom volume distributions and the volume table have to be consistent -- currently, the volumes in the param19-1.2.inp file are optimal for an ALPHa of 1.2. Changing ALPHa without adapting the volume table is expected to reduce the precision of the results. 4. SIGMa The ACE solvation potential includes a hydrophobic contribution which is roughly proportional to the solvent accessible surface area. The factor SIGMa scales the hydrophobic contribution. For peptides with about 10-15 residues, a SIGMa factor of 3 results in hydrophobic contributions that are approximately equal to the solvent accessible surface area multiplied by 8 cal/(mol*A*A).  File: ACE, Node: Examples, Up: Top, Previous: Function, Next: Top Examples To set up simulations/minimizations with the ACE solvation potential, the following energy call is expected to be adequate in most situations: ENERgy ATOM ACE IEPS 1.0 SEPS 80.0 ALPHa 1.2 SIGMa 3 SWITch - VDIS VSWI CUTNB 13.0 CTONNB 8.0 CTOFNB 12.0 When you run molecular dynamics or minimization with ACE, you get two more lines in the log file printout with energy terms, e.g., DYNA DYN: Step Time TOTEner TOTKe ENERgy TEMPerature DYNA PROP: GRMS HFCTote HFCKe EHFCor VIRKe DYNA INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers DYNA EXTERN: VDWaals ELEC HBONds ASP USER DYNA PRESS: VIRE VIRI PRESSE PRESSI VOLUme DYNA ACE1: HYDRophobic SELF SCREENing COULomb DYNA ACE2: SOLVation INTERaction ---------- --------- --------- --------- --------- --------- DYNA> 0 0.00000 -3423.29671 0.00000 -3423.29671 0.00000 DYNA PROP> 4.45310 -3423.12228 0.52327 0.17442 -532.70519 DYNA INTERN> 6.58717 60.43092 0.00000 56.00750 7.32144 DYNA EXTERN> -380.26218 -3173.38156 0.00000 0.00000 0.00000 DYNA PRESS> 0.00000 355.13679 0.00000 0.00000 0.00000 DYNA ACE1> 109.04469 -3829.20991 2750.59427 -2203.81062 DYNA ACE2> -1078.61564 546.78365 ---------- --------- --------- --------- --------- --------- and the same during minimization (MINI...). The terms in lines with ACE1 and ACE2 are: HYDRophobic: Hydrophobic potential, equivalent to a surface based solvation term proportional to the sigma input parameter; SELF: Self contribution to electrostatic solvation free energy, Delta-E_self, first term of eq(8) (i.e., sum over all atomic solvation energies, Delta-E_self_i, eq(28)); SCREENing: Interaction contribution to electrostatic solvation free energy, i.e., screening of Coulomb interactions, eq(38) (sum over all atom pairs, including bonded and 1-3 atom pairs!); COULomb: Coulomb energy with constant dielectric of EPSI (sum over all atom pairs for the first term in eq(36) -- excluding bonded and 1-3 atom pairs, and 1-4 atom pair contributions scaled with E14FAC); SOLVation: Electrostatic (!) solvation free energy, sum of SELF and SCREENing; INTERaction: Electrostatic interaction, sum of SCREENing and COULomb (eq(36), but taking account of the bonded, 1-3, and 1-4 exclusion in the Coulomb term, see above). The term "ELEC" in line "DYNA EXTERN>..." is the total of the electrostatic energy PLUS the hydrophobic solvation energy (may change this in the future to avoid confusion): ELEC: Sum of HYDRophobic, SELF, SCREENing, COULomb. Equation numbers refer to Schaefer & Karplus, J. Phys. Chem. 100 (1996), 1578. see also: test cases c26test/ace1.inp and c26test/ace2.inp. =============================================================================== CHARMM Element doc/adumb.doc 1.1  File: ADUMB, Node: Top, Up: (chmdoc/commands.doc), Next: Syntax Adaptive Umbrella Sampling Module Setting up of adaptive umbrella potentials. Currently supported types of umbrella potentials are functions of dihedral angles and functions of the potential energy of the system (energy sampling). WARNING: The module is still being developed and some details are likely to change in future versions. Please report problems to Christian Bartels at cb@brel.u-strasbg.fr REFERENCES: C. Bartels & M. Karplus, J. Comp. Chem. 18 (1997) 1450- C. Bartels & M. Karplus, J. Phys. Chem. 102 (1998) 865- M. Schaefer, C. Bartels, & M. Karplus, J. Mol. Biol. (1998) * Menu: * Syntax:: Syntax of the ADUMB commands * Function:: Purpose of each of the commands * Examples:: Usage examples of the ADUMB module  File: ADUMB, Node: Syntax, Up: Top, Previous: Top, Next: Function Syntax [SYNTAX ADUMB functions] Syntax: ADUMb DIHE NRES int TRIG int POLY int 4X(atom-spec) ADUMb ENER NRES int TRIG int POLY int MAXE real MINE real [MAXT real] [MINT real] ADUMb INIT NSIM int [UPDA int] [EQUI int] [TEMP real] [AGIN real] [NEXT int] [THRE real] [UCUN int] [WUNI int] [RUNI int] ADUMb PROB UCUN int [TEMP real] [PUNI int] [TUNI int] ADUMb STON ADUMb STOFf where: atom-spec ::= { segid resid iupac } { resnumber iupac }  File: ADUMB, Node: Function, Up: Top, Previous: Syntax, Next: Examples 0. Introduction The module provides commands to define degrees of freedom along which adaptive umbrella potentials are applied in molecular dynamics simulations. Statistics on the sampling of the degrees of freedom are recorded during the md simulations and periodically used to update the umbrella potential such that uniform sampling of the degrees of freedom can be expected. Currently, dihedral angles and the potential energy are supported as degrees of freedom. If several degrees of freedom are defined, multidimensional adaptive umbrella sampling is performed. Two sorts of input/output files are used by the module. The "umbrella" files contain the umbrella potentials that were used in the simulations together with the statistics of the sampling of the bins during the simulations. Based on this information the potential of mean force can be calculated and the umbrella potential expected to lead to uniform sampling can be determined. The second sort of files contains the values of the umbrella coordinates (=degree of freedom for adaptive umbrella sampling) for each time step in which coordinates were saved to the trajectory files. The umbrella coordinates are normalized to the range 0 to 1, independent of the degrees of freedom used. From the umbrella coordinates saved, weighting factors can be calculated which are needed to calculate average properties of the unbiased system. 1. ADUMb DIHE Define a dihedral angle as degree of freedom for adaptive umbrella sampling. To record the statistics the degree of freedom is partitioned into NRES bins. The umbrella potentials are represented as a linear combination of two times TRIG trigonometric functions and polynomial functions of degree 0 to POLY - 1. Repeating the command results in a multidimensional adaptive umbrella potential. The coordinates written to the umbrella coordinates file are normalized to the range 0 to 1 with 0 corresponding to -180 degrees and 1 corresponding to +180 degrees. 2. ADUMb ENER Define the potential energy as degree of freedom for adaptive umbrella sampling. NRES, TRIG and POLY have the same meaning as in ADUMb DIHE. MINE and MAXE specify the potential energy range: Statistics on the sampling are recorded in the range MINE-0.5*(MAXE-MINE) to MAXE+0.5*(MAXE-MINE). In the range outside of MINE to MAXE the umbrella potential is kept constant to prevent the system from leaving the range in which statistics are recorded. MINT and MAXT (default values: 273 K and 1000 K, respectively) are minimal and maximal temperatures to restrict sampling in the relevant temperature range. To set up a system, get a rough estimate of the potential energy of the system at the desired TMIN and TMAX (from short unbiased simulations at TMIN and TMAX). Set EMIN and EMAX to the values determined minus/plus a small tolerance, respectively. The coordinates written to the umbrella coordinates file are normalized to the range 0 to 1 with 0 corresponding to MINE-0.5*(MAXE-MINE) and 1 corresponding to MAXE+0.5*(MAXE-MINE). 3. ADUMb INIT Defines or redefines the parameters for adaptive umbrella sampling and initializes the umbrella potential. The umbrella potential is updated every UPDAte steps. After each update, no statistics are recorded for EQUI steps. For the remaining UPDA - EQUI steps, statistics on the sampling of the umbrella coordinates are recorded and stored separately from previous statistics and together with the umbrella potential active when recording the statistics. NSIM separate statistics can be kept in memory. If the number of updates performed in a run exceeds NSIM, the oldest statistics are discarded to make space for the most recent statistics. After each update the umbrella potential and the statistics are written to standard output (the log file). The written table contains, from left to right, the number of the bin, the number of integration time steps in which the system was in the bin since the last update, the potential of mean force calculated with the WHAM equations, the negative of the updated umbrella potential (potential of mean force modified to restrict sampling if necessary and fitted to the set of trigonometric and polynomial functions), the total number of times the bin was visited in the entire simulation, and the umbrella coordinates of the center of the bin. The temperature TEMP should be set to the temperature used in the simulations. It is used to calculate the umbrella potentials from the sampling statistics and to restrict sampling if potential energy sampling is performed. Umbrella coordinates are written to unit UCUN. At each update, the statistics are written to unit WUNI together with the umbrella potential active when recording the statistics. Statistics from previous runs can be read from unit RUNI. The statistics read must be from adaptive umbrella sampling simulations with the same parameters as the present one, in particular, the same degrees of freedom have to be used as umbrella coordinates. If adaptive umbrella sampling of the potential energy is used, umbrella potentials from runs at different temperatures can be read by repeating the ADUMb INIT command with RUNI set to the unit containing the statistics of each of the runs and TEMP set to the temperature of the run. To define the umbrella potential of bins for which no statistics have been acquired so far, the umbrella potential has to be extrapolated. In the current implementation (might change in future implementations), the umbrella potential of the bins that were not sampled is set to the same value (ext-cons). To determine ext-cons, the potential of the bins that were sampled is linearly extrapolated for NEXT bins, and the maximal value (max-extrapolated) of the linearly extrapolated potentials is determined. Then, the minimal value (min-sampled) of the potentials of the bins that were sampled is determined and ext-cons is set to min-sampled or max-extrapolated whatever value is smaller. A few statistics that differ significantly from the rest of the statistics can be due to problems with the convergence caused by the extrapolation or due to the occurrence of rare events. In the former case, outliers should occur only in the first few simulations and it is advantageous to eliminate them. By default, the module eliminates statistics that differ from the averaged statistics by THRE times the average deviation. If one wants to prevent statistics from being eliminated THRE has to be set to a value larger than NSIM. At each update, the deviations of the statistics from the averaged statistics is printed to standard output (log file), e.g., 0 Deviation of simulation 1 : 0.955 0 Deviation of simulation 2 : 0.513E-01 0 Deviation of simulation 3 : 0.787E-01 0 Deviation of simulation 4 : 0.292 0 Deviation of simulation 5 : 0.170 0 Deviation of simulation 6 : 0.201 0 Deviation of simulation 7 : 0.933 0 Deviation of simulation 8 : 0.208 0 Deviation of simulation 9 : 0.270 0 Deviation of simulation 10 : 0.131 0 Deviation of simulation 11 : 0.394 0 Deviation of simulation 12 : 1.52 0 Deviation of simulation 13 : 0.969 0 Deviation of simulation 14 : 0.502 0 Deviation of simulation 15 : 1.47 0 Deviation of simulation 16 : 2.97 -1 Deviation of simulation 17 : 210. 0 Deviation of simulation 18 : 0.695E-01 0 Deviation of simulation 19 : 0.160 0 Deviation of simulation 20 : 0.450 The 0 or -1 on each line indicates whether the statistics of a particular simulation are used (0) or were discarded (-1) based on the THRE criterion. For complex systems, there might exist no umbrella potential that enables the system to diffuse rapidly along the umbrella coordinate. In such cases it has been found to be advantageous to give a higher weight to the most recent statistics. This is implemented using the AGINg factor. For an umbrella potential calculated from n statistics, the i'th statistics (i=1,2,..,n) are weighted by AGINg**(n-i). 4. ADUMb PROB Average properties of the unbiased system can be obtained by weighting the conformations of an adaptive umbrella sampling run by appropriate factors. The ADUMb PROB command calculates these weighting factors from the umbrella coordinates read from unit UCUN and writes them to unit PUNI. For the command to work the umbrella potentials and statistics from the run must have been read with the ADUMb INIT command. If the potential energy was used as umbrella coordinate, the TEMP specifies the temperature at which properties of the unbiased system should be calculated. 5. ADUMb STON ADUMb STOFf By default statistics on the sampling of the umbrella coordinates are recorded in each call to the energy routines. The ADUMb STOFf command prevents that statistics are recorded. This might be useful when doing a minimization or running a md simulation with an umbrella potential that should not change during the simulation.  File: ADUMB, Node: Examples, Up: Top, Previous: Function, Next: Top Examples This examples are meant to be a partial guide in setting up an input file for ADUMB. There are three test files, adumb-phichi.inp, adumb-enum.inp and ace2.inp. Example (1) ----------- Set up and run an adaptive umbrella sampling simulation using two dihedral angles as umbrella coordinates. ! define the phi and chi1 dihedral angle as the two umbrella coordinates umbrella dihe nresol 36 trig 6 poly 1 pept 1 N pept 1 CA pept 1 CB pept 1 OG1 umbrella dihe nresol 36 trig 6 poly 1 pept 1 CY pept 1 N pept 1 CA pept 1 C umbrella init nsim 100 update 10000 equi 1000 thresh 10 temp 300 - ucun 10 wuni 11 ! perform adaptive umbrella sampling md simulation dynamics nose tref 300 qref 20 start - nstep 20000 timestep 0.001 - ihbfrq 0 inbfrq 10 ilbfrq 5 - iseed 12 - nprint 1000 iprfreq 1000 - isvfrq 1000 iunwrite -1 iunread -1 - wmin 1.2 Example(2) ---------- Set up and run an adaptive umbrella sampling simulation using the potential energy as umbrella coordinate (=energy sampling, multicanonical simulation, entropic sampling). ! set up umbrella; the range of relevant potential energies is assumed to ! extend form -50 kcal/mol to 100 kcal/mol. umbrella ener nresol 200 trig 20 poly 5 mine -50 maxe 100.0 mint 280 maxt 2000 open write formatted unit 9 name @9enum.umb open write formatted unit 10 name @9enum.uco open write unformatted unit 11 name @9enum.cor umbrella init nsim 100 update 10000 equi 1000 temp 1000 thres 100 - wuni 9 ucun 10 ! energy sampling simulation dynamics langevin start - nstep 50000 timestep 0.001 - inbfrq 10 ilbfrq 10 rbuffer 0.0 tbath 1000 - iseed 12 - nprint 1000 iprfreq 1000 - isvfrq 1000 iunwrite -1 iunread -1 - nsavc 100 iuncrd 11 - wmin 1.2 Example(3) ---------- Determine the weighting factors to calculate properties of the unbiased system. ! define the umbrella coordinates umbrella ener nresol 200 trig 20 poly 5 mine -50 maxe 100.0 mint 280 maxt 2000 open read formatted unit 10 name ../scr/@n.umb umbrella init nsim 100 update 10000 equi 1000 runi 10 temp 1000 thres 200 ! translate umbrella coordinates into probability factors at 300K open read formatted unit 11 name ../scr/@n.uco open write formatted unit 12 name ../scr/@nT300K.pfa umbrella prob ucun 11 puni 12 temp 300 ! translate umbrella coordinates into probability factors at 1000K open read formatted unit 11 name ../scr/@n.uco open write formatted unit 12 name ../scr/@nT1000K.pfa umbrella prob ucun 11 puni 12 temp 1000 CHARMM Element doc/analys.doc 1.1  File: analys, Node: Top, Up: (chmdoc/commands.doc), Next: Description Analysis Commands The ANALysis command is an energy and structure analysis facility that has been developed to examine both static and dynamic properties. The current code allows energy partition analysis and energy contribution analysis from free energy simulations. It also can produce a detailed printout of structural and energy term contributions for selected atoms * Menu: * Description:: Description of analysis facility * Energy:: Energy partitioning  File: analys, Node: Description, Up: Top, Previous: Top, Next: Energy Description of the ANALysis Command Syntax: ANALys { ON } { TERM { [ALL] } { NONBond } [UNIT int] atom-selection } { { ANY } { [NONOnbond] } } { OFF } ON Enable energy partition analysis and disable FAST routines. OFF Disable analysis and restore FAST option defaults. TERM Setup energy term print data and disable FAST routines. ALL (default) Print energy terms involving only selected atoms ANY Print energy terms when any of the atoms is selected NONBond In addition to internal terms, also print nonbond terms NONOnbond (default) Do not print electrostatic and vdw energy data UNIT integer Write the energy term printout data to a formatted file Otherwise, write data to the output file.  File: analys, Node: Energy, Up: Top, Previous: Description, Next: Top Energy term option of the ANALysis Command The ANALysis ON command enables energy partition analysis and disables the FAST routines. This will slow the calculation (especially on vector machines), but allows 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 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. 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. CHARMM Element doc/aspenr.doc $Revision: 1.1 $  File: ASPENR, Node: Top, Up: (chmdoc/commands.doc), Next: Syntax Atomic Solvation Parameter Based Energy Purpose: calculate solvation free energy and forces based on the exposed surface area of each atom using Atomic Solvation Parameters. Please report problems to brbrooks@helix.nih.gov REFERENCES: M. Wesson and D. Eisenberg, 19??. * Menu: * Syntax:: Syntax of ASP input * Structure:: Structure of the .surf file containing ASP data * Examples:: Usage examples of the ASP module  File: ASPENR, Node: Syntax, Up: Top, Previous: Top, Next: Structure Syntax [SYNTAX ASP functions] Syntax: The ASP specifications can be specified any time prior to an energy calculation and can be input either through reading a file or parsed directly off the command line - although the file route is more usual. Once turned on, the ASP energy term is in place during the course of the CHARMM run, i.e., it cannot be turned off except using the skipe command, see *note Skipe (chmdoc/energy.doc). Reading surf file: open unit 1 read vap_to_wat_kd.surf read surf unit 1 close unit 1  File: ASP, Node: Structure, Up: Top, Next: Examples, Previous: Syntax This module computes solvation energies and forces based on the surface area model proposed by Wesson and Eisenberg, i.e., E_solv = Sum (Gamma_i * ASA_i + Eref_i), where Gamma_i is a parameter describing the free energy cost of burying atom i (in units of cal/mol/A^2), ASA_i is the surface area of atom i with radius RvdW_i and probe radius Rprobe and Eref_i is a reference solvation energy. The analytic expressions for atomic surface areas and corresponding cartesian derivitives are used in these calculations. The values of the required parameters are read from a "surf" file which has the following syntax: * file: vap_to_wat_kd.surf * ! Note: These are asp's from Wolfendon water to vapor numbers, ! adjusted for standard state by Kyte and Doolittle. ! They are in units of cal/(mol*A**2). ! Table of ASP's ! 1.400000 ! the probe radius ! ! residue-type atom-name asp-value radius reference-area swap-pairs ANY H* 00.0 -1.0 0.0 ! ignore hydrogens ANY C 04.0 1.90 0.00 ANY OT1 -112.0 1.40 0.00 OT2 ANY OT2 -112.0 1.40 0.00 OT1 ANY N -112.0 1.70 0.00 . . . . TRP CZ2 04.0 1.90 0.00 TRP CZ3 04.0 1.90 0.00 TRP CH2 04.0 1.90 0.00 ASN OD1 -112.0 1.40 0.00 ASP OD1 -112.0 1.40 0.00 OD2 ASP OD2 -166.0 1.40 0.00 OD1 END Notes: -ANY refers to any residue type -A negative radius causes the atom to be ignored (such as hydrogens,...) -Atom name can use CHARMM wildcard rules (not residue names). -These commands ar eprocessed sequentially. If an atom is matched by more then one line the LAST line is used. -This file is free field format.  File: ACE, Node: Examples, Up: Top, Previous: Structure, Next: Top Examples To set up energy calculations/simulations/minimizations with the ASP potential, the following call is expected to be adequate in most situations: open unit 1 read form name vap_to_wat_kd.surf read surf unit 1 close unit 1 When you do an energy calculation, dynamics or minimization with ASP, you get columns in the log file printout with energy terms for ASP, e.g., ENER ENR: Eval# ENERgy Delta-E GRMS ENER INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers ENER EXTERN: VDWaals ELEC HBONds ASP USER ---------- --------- --------- --------- --------- --------- ENER> 0 -44.02560 7.74091 6.01738 ENER INTERN> 0.00000 0.04160 0.00000 0.00000 0.04556 ENER EXTERN> 5.95140 -42.32325 0.00000 -7.74091 0.00000 ---------- --------- --------- --------- --------- --------- and the same during minimization and dynamics. see also: test cases c27test/aspenr.inp 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. The lambda-dynamics, an alternative way of performing free energy calculations and screening binding molecules, has also been implemented. Subcommands related to BLOCK will be described here. To see how to output the results of a dynamics run, please see DYNAMICS documentation (keywords are IUNLDM, NSAVL, and LDTITLE). Please refer to PDETAIL.DOC for detailed description of the lambda dynamics and its implementation. 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 - [BOND real] [ANGL real] [DIHEdral real] [ELEC real] [VDW 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] LDINitialize int real real real real RMBOnd RMANgle LDMAtrix LDBI int LDBV int int int int real real int LDRStart LDWRite IUNL int NSAVL 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. The required real value (first specified) scales all energy terms expect those specific terms which are named with alternative corresponding scale factors. 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!) 13) LDINitialize specifies input parameters for running lambda dynamics. It sets up the value of lambda**2, the velocity of the lambda, its mass and reference free energy (or biasing potential). E.g, the following input lines set up parameters for the third lambda with [lambda(3)]**2 = 0.4, lambdaV(3) = 0.0, lambdaM(3) = 20.0, and lambdaF(3)=5.0 (note that lambdaF(1) should always be set to zero). LDIN 3 0.4 0.0 20.0 5.0 For more details, see Node Hints, section "lambda-dynamics simulations". 14) LDMAtrix will automatically map the input lambda**2 values onto the coefficient matrix of the interaction energies (and forces) between blocks. 15) LDBI provides an option on applying biasing potentials on lambda variables. The integer value specifies the total number of biasing potentials to be used. E.g, LDBI 3 will include total of 3 biasing potentials in the simulation. 16) LDBV sets up the specific form of the biasing potentials. At the moment, the functional form is of power law and allows three different classes (for details see "the actual simulations"). The input format is LDBV INDEX I J CLASS REF CFORCE NPOWER e.g. LDBV 2 2 3 3 0.0 50.0 4 will assign the second biasing potential acting between lambda(2) and lambda(3). The potential form belongs to the third class with reference value of zero, the force constant of 50 kcal/mol and the power of four. 17) LDRStart is used to restart the lambda dynamics runs. 18) LDWRite specifies the FORTRAN output unit No. and the frequency for writing lambda histogram by assigning an integer to IUNL and an integer to NSAVL. (IUNL and NSAVL can be reset in DYNAmic command, see *note dynamc:(chmdoc/dynamc.doc) ) 19) RMBOnd and RMANgle are used when no scaling of bond and angle energy terms is desired. END  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...) (iii) Lambda-dynamics simulations ================================= In an efforts to make the transition from using previous subcommands to running the lambda dynamics as smoothly as possible, we purposely parallel new syntax after the COEF subcommand. There are total of eight new keywords for setting up new dynamics. They are classified according to their functionality. (a) LDINitialize and LDMAtrix These two keywords are basic commands for starting the lambda dynamics run. The correct use of them is tied together with the BLOCK and CALL commands. Using the same example as the one given in "the actual simulations", the input script fragment will be as following: block 3 call 2 sele end call 3 sele end LDIN 1 1.0 0.0 20.0 0.0 LDIN 2 0.9 0.0 20.0 0.0 LDIN 3 0.1 0.0 20.0 0.0 LDMA end Here, the LDINitialize command models after the COEF command with the format LDIN INDEX LAMBDA**2 LAMBDAV LAMBDAM LAMBDAF Several comments are in order. First, notice that [lambda(1)**2] = 1.0 and [lambda(2)]**2 + [lambda(3)]**2 = 1.0. They are quite similar to the inputs of COEF subcommand. However, since one index instead of a pair is required here, only diagonal elements of the interaction coefficient matrix are specified. To fill up the matrix, LDMA is provided to finish the job automatically. In general, if there is total of N blocks, the first one is by default assumed to be the region where nothing changes. Therefore, [lambda(1)**2] = 1.0 is always true. The condition N ____ \ / [lambda(i)**2] = 1.0 (1) ---- i = 2 has to be satisfied for the partion of the system Hamiltonian. Due to some technical reasons in our implementation (details see PDETAIL.DOC), we have used [lambda(i)**2] instead of lambda(i) in our partion of the system Hamiltonian. Next, to make sure the above condition is met at any given simulation step, we have also enforced a condition containing velocities of the lambda variables N ____ \ / lambda(i)*lambdaV(i) = 0.0 (2) ---- i = 2 We used lambdaV(i) = 0.0 in the above script just to simplify the input. As far as the mass parameter lambdaM is concerned, the minimum requirement is that the value of mass has to be chosen such that the time step (or frequency) of lambda variables is consistent with that used for spatial coordinates x, y, z. Since the lambda variable is introduced into the system by using extended Lagrangian, considerations gone into the similar quantities, such as the adjustable parameter Q in a Nose thermostat are applicable to the choice of lambdaM. Some crude estimation can be made by examining the derivative of the system Hamiltonian with respect to the lambda, the curvature (simple harmonic approximation) or energy difference between two end-point states (0 and 1). Our experience has indicated that a conservative choice of the mass, i.e. a little bit heavier mass than that of the crude estimate, serves us well so far. The biasing potential LAMBDAF has two functions: (1) In the screening calculations LAMBDAF corresponds to the free energy difference of the ligands in the unbound state. Such calculations can identify ligands with favorable binding free energy and a ranking of the ligands can be obtained from the probability of each ligand in the lambda=1 state; (2) In precise free energy calculations, LAMBDAF corresponds to the best estimate of free energy from previous calculations. Therefore the estimate of free energy can be improved iteratively. (b) LDBI and LDBV In order to provide better control over simulation efficiency and sampling space, an option of applying biasing (or umbrella) potentials is furnished. LDBI specifies how many biasing potentials will be applied and LDBV supplies all the details. The general input format is LDBV INDEX I J CLASS REF CFORCE NPOWER Let us look at the following script block LDBI 3 LDBV 1 2 2 1 0.2 40.0 2 LDBV 2 3 3 2 0.6 50.0 2 LDBV 3 2 3 3 0.0 20.0 2 end It states that there is total of 3 biasing potentials. The first one (INDEX = 1) is acting on lambda(2) itself (I = J = 2), the second one on lambda(3) and the third one is coupling lambda(2) and lambda(3) together. At the moment, three different classes of functional forms are supported: CLASS 1: __ | CFORCE*(lambda - REF)**NPOWER if lambda < REF V =| | 0 otherwise |__ CLASS 2: __ | CFORCE*(lambda - REF)**NPOWER if lambda > REF V =| | 0 otherwise |__ CLASS 3: V = CFORCE*[lambda(I) - lambda(J)]**NPOWER (c) LDRStart LDRStart is used only if for some reason, e.g. execution of EXIT command, the logical variable QLDM for the lambda dynamics has been set to false. In this case, to restart the dynamics, LDRStart can be used to reset QLDM = .TRUE.. However, if LDIN is also being used in restarting the dynamics, it will automatically reset QLDM. Therefore, LDRS does not need to be called in this case. (d) LDERite LDWRite provides specifications for writing out lambda dynamics, i.e. the histogram of the lambda variables, the biasing potential etc. The integer variable IUNLdm is the FORTRAN unit on which the output data (unformatted) are to be saved. The value of the integer NSAVL sets step frequency for writing lambda histograms. IUNLdm is defaulted to -1 and NSAVL is defaulted to 0. Both IUNLdm and NSAVl can be reset in DYNAmics command (Please refer to *note dynamc:(chmdoc/dynamc.doc) for details). the following script will set IUNLdm with unit No. 8 and NSAVL equal to 10: LDWRite IUNL 8 NSAVL 10 (e) RMBOnd and RMANgle Since each energy term is scaled by lambda, RMBOnd and RMANgle can prevent bond breaking caused by such scaling during dynamic simulations. Alternatively one can fix bonds (and angles) using SHAKE. But is is not always possible. END  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/cadpac.doc $Revision: 1.2 $  File: Cadpac, Node: Top, Up: (chmdoc/commands.doc), Next: Description Combined Quantum Mechanical and Molecular Mechanics Method Based on CADPAC in CHARMM by Paul Lyne paul@tammy.harvard.edu * Menu: * Description:: Description of the CADPAC commands * Using:: How to run CADPAC in CHARMM * Installation:: How to install CADPAC in CHARMM environment * Status:: Status of the interface code  File: Cadpac, Node: Description, Up: Top, Next: Usage, Previous: Top The CADPAC QM potential is initialized with the CADPac command. [SYNTAX CADPac] CADPac [REMOve] [EXGRoup] (atom selection) REMOve: Classical energies within QM atoms are removed. EXGRoup: QM/MM Electrostatics for link host groups removed. The syntax of the CADPAC command in CHARMM follows closely that of the GAMESS command.  File: Cadpac, Node: Usage, Up: Top, Next: Status, Previous: Description For complete information about CADPAC input see Chapter 1 in the CADPAC distribution. A QM-MM job using CADPAC needs four input files. The first is the normal CHARMM input file containing the CADPac command. The second file is the CADPAC input file specifying the basis set to be used and the Hamiltonian that is needed. The third and fourth files are libfil.dat and modpot.dat respectively. These are the library and model potential files that are supplied with CADPAC. Cadpac Input File ----------------- For the CADPAC input file the following cards must be present: TITLE, BASIS, ATOMS, RUNTYP, START, FINISH. TITLE: The keyword is always at the start of the input file and is followed by a one-line title on the next line of the input. BASIS: This descirbes the basis set to be used for the QM region if a generic basis set is required. Examples include STO3G,321G,631G,321G*,631G*. These are the most common. Other basis sets are descibed in the CADPAC documentation. It is also possible to run a calculation using specific basis sets for individual atoms. If this feature is required then the BASIS keyword should be ommitted and the LIBRARY keyword is used for each atom in the QM region. For a more detailed description of the library command please refer to the official CADAPC documentation. All the basis sets that are supported by CADPAC are found in the files libfil.dat and modpot.dat. ATOMS: This keyword is always required. RUNTYP: For the purposes of QM-MM calculations this will either be ENERGY for a single point calculation or GRADIENT if the forces are also required. For any minimization or dynamics calculations the GRADIENT keyword should be used. START: This keyword is always required. FINISH: This keyword is always required. Hamiltonians ------------ The Hamiltonian is HF unless otherwise specified. The Hamiltonian can be changed by inseerting the appropriate keyword after the RUNTYP key. For example MP2 Performs an MP2 calculation MP3 Performs an MP3 calcualtion CI Performs a Configuration Interaction calculatio (please refer to the official CADAPC manual) For DFT calculations use the KOHNSHAM keyword: KOHNSHAM LDA MEDIUM GRDWT Performs an LDA calculation with a medium sized grid for numerical quadrature. KOHNSHAM BLYP LARGE GRDWT Performs a non-local BLYP calculation with a large sized grid For other functionals see the official CADPAC manual. CADPAC I/O ---------- CADPAC has hard wired units 1,2 and 3 for the libfil.dat, modpot.dat and cadpac input file so avoid using these elsewhere in the CHARMM stream. Other units that CADPAC commonly uses for the grid, integrals etc are 13,14,18,35,53,and 54. Examples -------- An example of a CADPAC input file to run with CHARMM: TITLE ! Required this is a test ! Put whatever you like on one line BASIS STO3G ! Generic basis set to be used ATOMS ! Required GRADIENT ! Run type. Use this for optimizations START ! Required FINISH ! Required The above input file tells CADPAC to use an STO-3G basis for the atoms in the QM region. CADPAC will perform a gradient evaluation each time that it is called by CHARMM. If you require just a single point calculation without gradients just use ENERGY instead of GRADIENT. The input file above will perform a HF calculation. A DFT calculation is invoked as follows: TITLE ! Required this is a test ! Put whatever you like on one line BASIS STO3G ! Generic basis set to be used ATOMS ! Required GRADIENT ! Run type. Use this for optimizations KOHNSHAM LDA MEDIUM GRDWT START ! Required FINISH ! Required DF jobs are invoked by the KOHNSHAM card which takes the type of functional and grid to be used as arguments. In this case an LDA functional is used. Alternatives include BLYP, B3LYP. For details see the CADPAC distribution. A sample shell script to run CHARMM with CADPAC is: #!/bin/tcsh -f # parameters: # 1 data file name # echo starting date echo $1 set HOME= {where CADPAC data files are} # data set and output in home directory # set data=$HOME/$1.inp set output=$HOME/$1.out2 # make a temporary directory to hold the workfiles cd /tmp mkdir $1 cd $1 # basis set library file assigned to fort.1 # pseudopotential library on fort.2 # the CADPAC input file is copied to UNIT 3 cp $HOME/$1.str fort.3 cp $HOME/$1.par . cp $HOME/libfil.dat fort.1 #cp $HOME/modpot.dat fort.2 # # run the program charmm.exe < $data rm -r ../$1 An example file can be found in test/c25test/cwat.inp. This input file also uses cwat.str and the sample run script runcwat.  File: Cadpac, Node: Status, Up: Top, Next: Top, Previous: Usage CADPAC/CHARMM interface status (February 1997) - CADPAC, GAMESS and QUANTUM keywords cannot coexist in pref.dat - CADPAC recognizes atoms by their masses as specified in the RTF file - The program runs on ALPHA, SGI, C90, IBMRS, HPUX platforms. - There are references to a parallel version in the code. This has not been fully tested yet and so won't be included until a future release. CHARMM Element doc/cff.doc 1.1 ^_ File: CFF, Node: Top, Up: (chmdoc/commands.doc), Next: Usage Consistent Force Field (CFF) * Menu: * Usage:: How to use CFF with CHARMM standalone * Status:: Current status of CFF implementation in CHARMM * Theory:: Basis for, parameterization and performance of CFF * Funcform:: Functional form of the CFF energy expression * Refs:: References to papers describing CFF ^_ File: CFF, Node: Usage, Up: Top, Next: Status, Previous: Top In order to use CFF in CHARMM, the user has to issue the following commands: 1. use cff 2. read cff parameter file 3. (a) read rtf name , or (b) read psf name 4. read sequence ! if input is via the rtf route (step 3 (a)) 5. generate 6. read coord, or ic build ! if input is via the read rtf/sequence route. When using CFF95 or later Step 3a requires a CFF-capable rtf file. This means a file in which BOND records have been replaced by analogous DOUBLE records for cases in which the chemical structure has a double bond. Note that CFF-capable rtf files are *back compatible*. That is, such rtf files can equally well be used for calculations that utilize the CHARMM force field. Thus, it is *not* necessary to maintain two versions of the rtf files. NOTE: 1. no binary parameter files are supported for CFF. 2. CFF is an all hydrogen force field -- i.e., extended atoms are not supported Examples of CFF usage in CHARMM are given in the ccfftest directory. ^_ File: CFF, Node: Status, Up: Top, Next: Theory, Previous: Usage Status of CFF implementation into CHARMM (May 1999) ============================================================= This implementation of CFF in CHARMM is principally due to Rick Lapp (MSI) and William Young (MSI). Features currently supported in CHARMM/CFF (1) energy and first derivatives (2) minimization (3) dynamics (4) most ATOM based cutoff options Major features NOT currently implemented in CHARMM/CFF: (1) second derivatives (2) bonds between primary atoms and image atoms. (3) Cutoff options currently not supported are group-based cutoffs, distance shifting and force-based switching. (4) Fast multipoles. Other known limitations: (1) correlation analysis tools have not been implemented for CFF specific energy terms -- e.g. it is not possible to calculate the correlation function for an out-of-plane bending angle, etc ... (2) only all-atom models (no extended atoms) There are probably other problems/limitations/bugs. Your comments about limitations of the current CFF implementation in CHARMM (and bugs) will be very valuable. Please direct comments to: William Young, MSI e-mail: wyoung@msi.com phone: (619)799-5348 KNOWN BUGS: ^_ File: CFF, Node: Theory, Up: Top, Next: Refs, Previous: Status The aim of the CFF development is a force field that is: * broad, covering a relatively large number of differing functional groups, * accurate, achieved via accurate reproduction of the quantum mechanical energy surfaces, * consistent between differing phases and molecular environments, * applicable to a wide range of molecular properties, * consistent between differing types of molecules, such as interaction of protein active sites with ligands, or assemblies of proteins with nucleic acids or with solvent. Quantum mechanical forcefields The intramolecular parameters constituting the current generation of forcefields are based on the energies and energy derivatives computed by ab initio quantum mechanical procedures for a series of model compounds. CFF uses quantum computations in the Hartree-Fock approximation with the 6-31G* basis set to expand the wavefunctions [1][2]. The quantum mechanical energies and the energy first derivatives (gradients) and second derivatives (Hessians) were computed for the equilibrium molecular structures, at conformational energy barriers, and for a set of distorted structures. The distorted molecular structures were generated by randomly deforming all the internal coordinates, as well as by systematically rotating about individual bonds. These quantum observables were fit to the energy expression to obtain the Class II parameters [3][4]. Many of the atomic partial charges were also determined quantum mechanically. The intermolecular parameters of the forcefield may also, in principle, be computed quantum mechanically [5]. The remaining CFF forcefield intermolecular or nonbond parameters were computed by fitting to experimental crystal lattice constants and sublimation energies of crystals [6][7][8]. Internal energy terms The energy of the molecule or assembly is expressed in terms of internal coordinates such as bond lengths, bond angles, and dihedral angles. For Class II forcefields this set of descriptors is greatly expanded by including cross terms, that is, the interactions between bond lengths and angles, between pairs of angles, etc. CFF contains, in all, twelve types of energy terms: bond stretching, valence angle bending, valence dihedral angles, out-of-plane deformation, and eight cross terms. The cross terms extend the accuracy and range of application of the forcefield by including the effect of neighboring atomic positions on each of the bond lengths, valence angles, and dihedral angles. ^_ File: MMFF, Node: Funcform, Up: Top, Next: Refs, Previous: Theory Energy functional forms The energy expression may be decomposed into diagonal terms that depend on a single molecular internal coordinate such as a bond length, coupling terms between internal coordinates, and nonbond internuclear distances. This energy is fit to the quantum mechanical energy. 1. Bond stretching. Ebond = K2 * (b - b0)^2 + K3 * (b - b0)^3 + K4 * (b - b0)^4 (1) where K2, K3 and K4 are the quadratic, cubic and quartic forcefield parameters or force constants, b is the bond length, and b0 is the reference value of the bond length. 2. Angle bending. Eangle = K2 * Delta^2 + K3 * Delta^3 + K4 * Delta^4 (2) where Delta = Theta - Theta0 is the difference between the actual and reference bond angles. 3. Out-of-plane bending. Eoop = K * (Chi - Chi0)^2 (3) where chi is an out-of-plane coordinate as defined by Wilson et al.[9] 4. Torsion energy, in order to reflect differing hybridizations about the bonded atoms, must contain one-, two-, and threefold periodic terms: Etorsion = SUM(n=1,3) { V(n) * [ 1 - cos(n*Phi - Phi0(n)) ] } (4) where phi is a dihedral angle. 5. Stretch-Stretch interaction between two bonds in a valence angle. Ebond-bond = K(b,b') * (b - b0) * (b' - b0') (5) 6. Stretch-Bend interaction between an angle and its bonds. Ebond-angle = K * (b - b0) * (Theta - Theta0) (6) 7. Bend-Bend-Twist interaction between a dihedral angle and its two valence angles. Eangle-angle-torsion = K * (Theta - Theta0) * (Theta' - Theta0') * (Phi - Phi1(0)) (7) 8. Stretch-Twist interaction between a dihedral angle and its end bonds. Eend_bond-torsion = (b - b0) * SUM { V(n) * cos[n*phi] } (8) 9. Stretch-Twist interaction between a dihedral angle and its middle bond. Emiddle_bond-torsion = (b - b0) * { F(1) * cos(phi) + F(2) * cos(2 * phi) + F(3) * cos(3 * phi) } (9) 10. Bend-Twist interaction between a dihedral angle and its valence angles. Eangle-torsion = (Theta - Theta0) * { F(1) * cos(phi) + F(2) * cos(2 * phi) + F(3) * cos(3 * phi) } (10) 11. Bend-Bend interaction between two valence angles with a common vertex atom. Eangle-angle = K * (Theta - Theta0) * (Theta' - Theta0') (11) 12. Stretch-Stretch interaction between the two end bonds in a dihedral angle. Ebond-bond_1_3 = K(b,b') * (b - b0) * (b' - b0') (12) Finally, the nonbond energy between atoms in different molecules or between atoms separated by three or more bonded atoms is given by the sum of the Coulombic electrostatic interaction and a van der Waals energy of the 9-6 form: 13. Coulombic electrostatic interaction. Ecoul = 332.0716*qi*qj/(D*Rij) (13) where qi and qj are the atomic partial charges on atoms i and j, Rij is the distance between them and D is the dielectric constant. 14. Van der Waals interaction. Evdw = eps(ij) [2*r*(ij)/r(ij)**9 - 3*r*(ij)/r(ij)**6] (14) where r*(ij) = [(r(i)**6 + r(j)**6))/2]**(1/6) (15) eps(ij) = 2 sqrt(eps(i) * eps(j)) * r(i)^3 * r(j)^3/[r(i)^6 + r(j)^6] (16) where eps(ij) and r*(ij) are the negative of the minimum van der Waals energy and that distance between atoms i and j where the minimum occurs, respectively. Eps(ij) and r*ij are computed from the individual atomic parameters eps(i), eps(j), r*i, and r*j by the Waldman-Hagler combination rules [10]. The Hartree-Fock method, and to a lesser extent other quantum mechanical methods, results in systematic deviations from experiment. For example, bond lengths tend to be too short and bond-stretching vibrational frequencies too high [11]. However, by comparison with experimental gas-phase molecular structures and vibrational frequencies, these deviations may be compensated for. In general, the energy expression may be scaled using five constant factors, one for each of the classes of energy terms: bonds, angles, torsion angles, out-of-planes and all coupling terms [12]. The scaled energy is then: Ediagonal = Sb * SUM{Ebond} + Stheta * SUM{Eangle} + Sphi * SUM{Etorsion} + Schi * SUM{Eoop} (17) Ecross = Sc * SUM{eight cross terms} (18) The reference values b0 and q0 are also adjusted to fit experimental data. All these values may differ among different types of bonds, bond angles, and torsion angles. For the special case of hydrocarbons, the corrections are especially well determined by gas-phase measurements. For hydrocarbons, the best values of the scale factors are: Sb(C-C) 0.88 Sb(C-H) 0.83 Stheta (all angles) 0.81 Sphi (all torsions) 0.84 Schi (all out-of-planes) 1.00 Sc (all cross terms) 0.87 The reference bond lengths for hydrocarbons were also adjusted. Although the use of the quantum calculation greatly amplifies the available data so that only a few such corrections are necessary for the complete Class II forcefield, for the majority of functional groups (molecular types) no accurate gas-phase data are available. However, the Sb, Stheta, Sphi, and Sc constants are transferrable among different types of bonds, bond angles, and torsion angles. Therefore, the same scale factors are used in Eq. 17 and Eq. 18 in the final empirically scaled forcefield. In general, the reference values b0 and theta0 are determined from high-level quantum mechanical calculations on the model compounds. Validation of the CFF forcefield Table 1 shows the accuracy of the CFF forcefield for several common classes of molecules, compared with experimental gas-phase results. Table 1. Summary of rms deviations between experimental and CFF-calculated structural parameters, vibrational frequencies, and energy differences. bond valence torsion freq. energy length angle angle diff. (Ang) (deg) (deg) (cm-1) (kcal mol-1) hydrocarbons 0.02 0.9 1.2 40 0.93 alcohols 0.02 1.7 1.7 37 0.71 aldehydes & ketones 0.01 1.1 2.3 32 0.62 amines 0.00 0.9 -- 18 0.62 carboxylic acids 0.02 1.6 1.0 34 0.78 esters 0.02 1.7 0.5 42 1.88 ethers 0.01 0.9 1.1 41 0.40 heterocycles 0.01 1.0 0.0 35 --- sulfides 0.01 1.4 2.5 45 --- disulfides 0.01 0.9 2.0 43 --- thiols 0.01 1.6 1.0 -- 0.21 average 0.01 1.2 1.3 37 0.77 The frequencies are harmonic vibrational frequencies and the energy differences include conformational energy differences and energy barriers to internal rotation between stable conformers. ^_ File: MMFF, Node: Refs, Up: Top, Previous: Funcform, Next: Top References [1] Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 28, 213-222 (1973). [2] Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 77, 3654-3665 (1982). [3] Dinur, U.; Hagler, A. T. In Reviews in Computational Chemistry, Vol. 2, K. B. Lipkowitz; D. B. Boyd, Eds., VCH Publishers: New York, 99-164 (1991). [4] Maple, J. R.; Hwang, M.-J.; Stockfisch, T. P.; Dinur, U.; Waldman, M.; Ewig, C. S.; Hagler, A. T. J. Comp. Chem. 15, 162-182 (1994). [5] Dinur, U.; Hagler, A. T. J. Amer. Chem. Soc. 111, 5149-5151 (1989). [6] Hagler, A. T.; Huler, E.; Lifson, S. J. Amer. Chem. Soc. 96, 5319-5327 (1974). [7] Hagler, A. T.; Lifson, S.; Dauber, P. J. Amer. Chem. Soc. 101, 5122-5130 (1979a). [8] Hagler, A. T.; Dauber, P.; Lifson, S. J. Amer. Chem. Soc. 101, 5131-5141 (1979b). [9] Wilson, E. B., Jr; Decius, J. C.; Cross, P. C., Molecular Vibrations; Dover: New York, 1955, Chapter 4. [10] Waldman, M.; Hagler, A. T. J. Comp. Chem. 14, 1077-1084 (1993). [11] Michalska, D.; Schaad, L. J.; Carsky, P.; Hess, Jr., B. A.; Ewig, C. S. J. Comp. Chem., 9, 495 (1988). [12] Hwang, M.-J.; Stockfisch, T. P.; Hagler, A. T. J. Amer. Chem. Soc. 116, 2515-2525 (1994). CHARMM Element doc/cfti.doc $Revision: 1.1 $  File: CFTI, Node: Top, Up: (chmdoc/perturb.doc), Next: Constraints CFTI: conformational energy/free energy calculations * Menu: * Constraints:: Note on constrained optimization implementation * CFTINT:: Description and syntax of standard conformational free energy thermodynamic integration * CFTIM:: Description and syntax of multidimensional onformational free energy thermodynamic integration  File: CFTI, Node: Constraints, Up: Top, Previous: Top, Next: CFTINT Constraints: Energy minimization with holonomic constraints has been implemented. There are no special commands for this option. Charlie Brook's TSM module allows for MD simulations with constrained values of selected conformational coordinates - distances, atoms, dihedrals. This has been expanded to also allow energy minimization using several algorithms. The method is an alternative to using harmonic restraints in generating structures of flexible molecules with desired properties, or generating adiabatic profiles. To use this option, simply enter the 'TSM' module and give set of 'FIX' commands to define set of fixed internal coordinates (see perturb.doc for details). Next specify an energy minimization (see minmiz.doc). Algorithms that work: SD, CONJ, POWE (ABNR works also, for reasons unclear to me, KK)  File: CFTI, Node: CFTINT, Up: Top, Previous: Constraints, Next: CFTIM CFTI: standard (one-dimensional) conformational thermodynamic integration Description of method Method expands the capabilities of the TSM module. The TSM module employs the Thermodynamic Perturbation (TP) approach to conformational free energy simulations. The basis of the calculation is a MD simulation with a constrained value of a conformational coordinate. With minimal modifications, the alternative Thermodynamic Integration (TI) method is added on. In the modified code the user has the option of using TP only (as previously) or activating TI, in which case the same simulation and data files are used to give both TP and TI results. [SYNTAX CFTI] All commands are parsed by the TSM command parser, so should be within a 'TSM ... END' block. CFTI command activates the thermodynamic integration calculation the context of use should be the same as for a thermodynamic perturbation run, i.e. some coordinates should be fixed by 'FIX ,,,', saving data to a disk file should be specified by 'SAVI ...', and one perturbation should be defined by 'MOVE ...' The derivative dA/dx is calculated for the coordinate deifned in the 'MOVE ...' statement. This coordinate has to be also fixed with 'FIX ...'; other coordinates may also be fixed if desired. See test case cftigas.inp for details Notes: 1) The formatted data file generated by 'SAVI ...' may be read by both TI postprocessing command (CFTJ) and TP postprocessing (POST). The SAVI 'NWIN' keyword has meaning only for TP, it can be set to an arbitrary value if TI only is to be used. 2) For consistency with TP, the 'BY ' part of the 'MOVE' command was retained. The value has meaning for TP only, it can be set to an arbitrary number if TI only is to be used. CFTJ [TEMP ] [UICP ] [CONT ] Command to calculate the conformational free energy derivative dA/dx = as well as the energy-entropy components: d/dx, -TdS/dx Data is read in from the formatted file generated by the 'SAVI ...' command TEMP - specifies temperature, needed for energy-entropy components UICP - specifies unit with data CONT - defines length of data block for error analysis e.g. if data file has 1000 entries, 'CONT 100' will divide data into 10 blocks and calculate the standard deviation of the mean of the block averages CFTA [FIRSt ] [NUNIt ] [BEGIn ] [STOP ] [SKIP ] [CONT ] [TEMP ] Command activates analysis of CFTI-generated trajectory. Trajectory coordinate file(s) should be in consecutive units FIRST, NUNI, BEGIN, STOP, SKIP - define trajectory reading CONT, TEMP - as in CFTJ Examples of usage : see test case cftigas.inp  File: CFTI, Node: CFTIM, Up: Top, Previous: CFTINT, Next: Top CFTM: multidimensional conformational thermodynamic integration Description of method This is a new approach. MD simulations are performed with several conformational coordinates simultaneously constrained to fixed values. The partial derivatives of the conformational free energy with respect to all the coordinates in the fixed set are calculated from this one simulation. The free energy gradient may be used in different ways to explore conformational free energy surfaces of flexible molecules. Method expands the capabilities of the TSM module. Only TI calculations possible, no corresponding TP analysis possible. [SYNTAX CFTM] All commands are parsed by the TSM command parser, so should be within a 'TSM ... END' block. CFTM command activates the multidimensional TI method the context of use should be the same as for a thermodynamic perturbation run, i.e. several coordinates should be fixed by 'FIX ,,,', saving data to a disk file should be specified by 'SAVI ...'. and a perturbation should be defined by a 'MOVE ...' statement for each of the fixed coordinates. See test case cftmgas.inp for details Notes: 1) The formatted data file generated by 'SAVI ...' may be read by both TI postprocessing command (CFTJ) and TP postprocessing (POST). The SAVI 'NWIN' keyword has no meaning in CFTM. The CFTM data file format is different than that for CFTI. 2) For consistency with TP, the 'BY ' part of the 'MOVE' command was retained. The value has no meaning in CFTM. 'INTE' keyword has to be specified within the 'MOVE' command. CFTC [TEMP ] [UICP ] [CONT ] Command to calculate the conformational free energy derivatives dA/dx_i = as well as the energy-entropy components: d/dx_i, -TdS/dx_i Data is read in from the formatted file generated by the 'SAVI ...' command TEMP - specifies temperature, needed for energy-entropy components UICP - specifies unit with data CONT - defines length of data block for error analysis e.g. if data file has 1000 entries, 'CONT 100' will divide data into 10 blocks and calculate the standard deviation of the mean of the block averages Output includes all individual partial derivatives, and optionally their analysis into groups. The derivative with respect to a path direction is also calculated. CFTB [FIRSt ] [NUNIt ] [BEGIn ] [STOP ] [SKIP ] [CONT ] [TEMP ] Command activates analysis of CFTM-generated trajectory. Trajectory coordinate file(s) should be in consecutive units FIRST, NUNI, BEGIN, STOP, SKIP - define trajectory reading CONT, TEMP - as in CFTJ Output is the free energy gradient with respect to the set of fixed coordinates, the derivative along a specified direction (see DIRE) and optionally a group contribution analysis. CFTS [FIRSt ] [NUNIt ] [BEGIn ] [STOP ] [SKIP ] [CONT ] [TEMP ] [DUNI ] Analogous to CFTB, additionally writes out potential energy and dU/dx_i to a disk file specified by DUNI. NCOR NUMB NUMB specifies the number of internal coordinates involved (=NICP). Used in calculating the path derivative. DIRE LAMB The LAMB value specifies number of step (progress along reaction path). The following line(s) contain NICP real numbers defining a path vector. The vector will be normalized in side the program. The unit vector will be used to calculate derivatives of dA/dl, d/dl, -TdS/dl along the path from the gradients. Note: the vector components are read in free format CFTG NGRUp Define groups for group contribution analysis to free energy NGRUP is the number of groups. The following line(s) contain the integer group numbers of the coordinates (LGRUP(J),J=1,NICP) in free format After that follow line(s) with group symbols (i.e. tags that will be used to denote the groups) in (20A4) format (GSYM(J),J=1,NGRUP) Example of usage: The system is a decapeptide, we calculate derivatives with respect to all phi and psi backbone dihedrals (NICP=18). In the 18 'MOVE ...' commands we specify the 9 phi first and the 9 psi at the end. The following will calculate and print out an aggregate of all phi and all psi contributions labelled by the tags 'PHI' and 'PSI': cftg ngrup 2 1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2 PHI PSI cfts Note on sign of derivatives: In both CFTI and CFTM it is possible to obtain a derivative value with incorrect sign by cleverly manipulating the atom selections in the 'MOVE ...' command. A simple way of checking the sign is to run a 1-D test case using both TI and TP postprocessing (see test case cftigas.inp). A general rule is to think about how the coordinate is defined and how motions of fragments influence it. E.g. for a distance between atoms A and B, the coordinate is the length of the vector from A to B. Perturbations (TP) involve actual displacements of A and B along the =vector= from A to B; Derivative calculations (TI) do not involve actual motions of atoms, but rather predictions of how atomic positions will vary with infinitesimal coordinate changes. Moving B along this coordinate by delta > 0 will increase the coordinate, while moving A by delta will decrease the coordinate. I.e. to get correct sign of derivative you have to use the following scheme: either FIX DIST MOVE DIST BY 1.0 INTE - sele end or FIX DIST MOVE DIST BY 1.0 INTE - sele end sele end see test cases cftifas.inp and cftmgas.inp for examples. 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 each node are recorded by CHARMM developers to indicate new and modified features of CHARMM during the development cycle, i.e., the alpha version period. ------------------------------------------------------ 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 c24b2 Release February 15, 1996 c24g1 Release August 15, 1996 c24g2 Release February 15, 1997 CHARMM 25.0 c25a0 Developmental August 15, 1995 c25a1 Developmental February 15, 1996 c25a2 Developmental August 15, 1996 c25a3 Developmental February 15, 1997 c25b1 Release August 15, 1997 ------------------------------------------------------ * 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 * C24-C25:: Major enhancements and developments in CHARMM25  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 w