CHARMm Principles



2       Preparing Models for Energy Calculations

 A molecular model must be specified before energy calculations can be performed. Model information, ultimately stored in the principal structure file, comes or is generated from a variety of sources including residue topology and sequences files. Cartesian coordinates must be defined and a set of parameters, including force constants and equilibrium geometries, is also needed. This chapter describes the information that must be available for each molecular model and the sources of that information.

The chapter explains:

Script examples are included throughout the chapter.


Constructing and using residue topology files

Information about residues, the basic chemical units that compose all models, is stored in residue topology files (RTFs). The atoms, atomic properties, bonds, bond angles, torsion angles, improper torsion angles, hydrogen bond donors, acceptors, antecedents, and non-bonded exclusions are all specified on a per residue basis.

The purpose of residue topology files is to store the information for generating a representation of a molecule from its sequence. RTFs can be created by editing a file and inserting the appropriate commands and keywords to describe the molecule. Typically, these files are automatically created by QUANTA or Insight II when you construct or modify a molecule. Libraries of topology files for amino acids, nucleic acids, saccharides, and polymers are also supplied with CHARMm.

As with all text files read into CHARMm, the first section of an RTF is a title in the usual format delimited by a line containing only an asterisk in the first column.

The remaining information is read in free-field format as commands to define the RTF. The ordering of the commands is important. For example, atoms of a residue must be defined before the bonds between the atoms can be defined.

File structure

If you are constructing an RTF, the recommended structure for the file is as follows:

Title lines

Version numbers

Mass specification for each atom type

Declarations of out-of-residue definitions

Defaults for patching on the first and last residues

Autogeneration of angles or dihedrals

Name and total charge specification

Atom definitions within the residue

Group definitions of atoms

Bond specifications

Angle specifications

Dihedral angle specifications

Improper dihedral angle specifications

Donor specifications

Acceptor specifications

Internal coordinate information

Patching residues to use if defaults are not desired

Example: Creating an RTF containing a water residue

An RTF with one residue for water is written as follows:


* Topology definitions for HOH
*
18 1
MASS 1 H 1.00800
MASS 56 OH2 15.99940
AUTOGENERATE ANGLES DIHEDRALS
RESI WAT 0.00000
GROUP
ATOM OH2 OH2 -0.40000
ATOM H1 H 0.20000
ATOM H2 H 0.20000
BOND OH2 H1
BOND OH2 H2
DONO H1 OH2 -O -O
DONO H2 OH2 -O -O
ACCE OH2
PATCH FIRST NONE LAST NONE
END
This is a simple example. The mass of the atom types is set first. The residue name, WAT, is defined and its charge listed.

For the residue to be recognized by CHARMm, the residue name used in the sequence must match that in the RTF. The same definition for WAT might also be used for H2O, HOH, and OH2, residue names commonly used in PDB files for water molecules.

The ATOM keyword matches the atom name to the atom type and defines the partial charge of that atom. The BOND keyword defines the atoms that are bonded. The DONOr and ACCEptor keywords affect the hydrogen bonding behavior of the molecule.

Example: Creating an RTF containing an alanine residue

An RTF containing a single alanine residue is written as follows:


* Topology definitions for ALA
*
18 1
MASS 1 H 1.00800
MASS 2 HA 1.00800
MASS 11 C 12.01100
MASS 15 CT 12.01100
MASS 31 NP 14.00670
MASS 51 O 15.99940
MASS 56 OH2 15.99940
DECL -C
DECL -O
DECL +N
DECL +H
DECL +CA
DEFA FIRST NTERM LAST CTERM
AUTOGENERATE ANGLES DIHEDRALS
RESI ALA 0.00000
GROUP
ATOM N NP  -0.35
ATOM H H    0.25
ATOM CA CT  0.00
ATOM HA HA  0.10
GROUP
ATOM CB CT -0.30
ATOM HB1 HA 0.10
ATOM HB2 HA 0.10
ATOM HB3 HA 0.10
GROUP
ATOM C C    0.45
ATOM O O   -0.45
BOND N CA N H CA C  CA HA CA CB C O C +N CB ND1 CB NB2- CB NB3
IMPH N -C CA H 
IMPH C CA +N O
DONO H N -C CA
ACCE O C
IC -C CA *N H  0.0000 0.00 180.00 0.00 0.0000
IC -C N CA C   0.0000 0.00 180.00 0.00 0.0000
IC N CA C +N   0.0000 0.00 180.00 0.00 0.0000
IC +N CA *C O  0.0000 0.00 180.00 0.00 0.0000
IC CA C +N +CA 0.0000 0.00 180.00 0.00 0.0000
IC N C *CA CB  0.0000 0.00 120.00 0.00 0.0000
IC N CA CB HB1 0.0000 0.00 180.00 0.00 0.0000
Two major differences exist between the RTFs for water and for alanine. The first is the statement: 


DEFA FIRST NTERM LAST CTERM
This specifies that the first amino acid residue in a sequence will be converted to the N-terminus and the last will be converted to the C-terminus automatically with a special patch for the N- and C-termini. For more information about using patches, see Modifying a PSF.

The second difference is the addition of special atoms in the alanine RFT declared as -C, -O, +N, +H and +CA. These refer to out-of-residue atoms belonging to the preceding (-) and following (+) residue in the sequence. These atoms are then used in the BOND, ANGLE, DIHEDRAL, and IMPROPER statements. In this case, the angles and dihedral angles are generated automatically. Each statement adds to the list of bonds and angles that are included in an energy calculation.

Example: Creating an RTF definition for glycerol

This CHARMm input file uses RTF definitions to create a model of glycerol. This example uses or generates the following files:


*...
* Copyright (c) 1994
* Molecular Simulations Inc. 
* All Rights Reserved
*
! Read the RTF from the CHARMm command file
READ RTF CARD
* RTF FOR GLYCEROL
*
 20 1! Version number
MASS 1 H 1.00800! Hydrogen bonding hydrogen (neutral group)
MASS 3 HA 1.00800! Aliphatic or aromatic hydrogen
MASS 10 CT 12.01100! Aliphatic carbon (tetrahedral)
MASS 45 OT 15.99940! Hydroxyl oxygen (tetra.)/Ionizable acid -			 oxygen
! Autogenerate all angles in the topology definitions
AUTOGEN ANGLES
! During generation, provide no patches on the terminal
! residues
DEFAULT FIRST NONE LAST NONE
! Topology definition for glycerol
RESIDUE GLYC 0.0 
GROUP             !     H11 O1 -- H1
 ATOM C1 CT 0.05  !      \ /
 ATOM H11 HA 0.10 !       C1
 ATOM H12 HA 0.10 !      / \
 ATOM O1 OT -0.65 !    H12  \
 ATOM H1 H 0.40   !         |
GROUP             !         |
 ATOM C2 CT 0.15  !         |
 ATOM H21 HA 0.10 !  H21 -- C2 -- O2 -- H2
 ATOM O2 OT -0.65 !         |
 ATOM H2 H 0.40   !         |
GROUP             !         |
 ATOM C3 CT 0.05  !     H31 /
 ATOM H31 HA 0.10 !      \ /
 ATOM H32 HA 0.10 !       C3
 ATOM O3 OT -0.65 !      / \
 ATOM H3 H 0.40   !     H32 O3 -- H3
BOND C1 H11 C1 H12 C1 O1 O1 H1
BOND C1 C2
BOND C2 H21 C2 O2 O2 H2
BOND C2 C3
BOND C3 H31 C3 H32 C3 O3 O3 H3
DIHE H1 O1 C1 C2
DIHE O1 C1 C2 C3
DIHE H2 O2 C2 C1
DIHE O3 C3 C2 C1
DIHE H3 O3 C3 C2
! Internal coordinate definitions
IC H1 O1 C1 C2 0.0000 0.00 180.00 0.00 0.0000
IC O1 C1 C2 C3 0.0000 0.00 180.00 0.00 0.0000
IC H2 O2 C2 C1 0.0000 0.00 180.00 0.00 0.0000
IC O3 C3 C2 C1 0.0000 0.00 180.00 0.00 0.0000
IC H3 O3 C3 C2 0.0000 0.00 180.00 0.00 0.0000
IC O1 C2 *C1 H11 0.0000 0.00 120.00 0.00 0.0000
IC O1 C2 *C1 H12 0.0000 0.00 240.00 0.00 0.0000
IC C3 C1 *C2 O2 0.0000 0.00 120.00 0.00 0.0000
IC C3 C1 *C2 H21 0.0000 0.00 240.00 0.00 0.0000
IC O3 C2 *C3 H31 0.0000 0.00 120.00 0.00 0.0000
IC O3 C2 *C3 H32 0.0000 0.00 240.00 0.00 0.0000
END
! Open and read the parameter file (binary)
OPEN READ UNIT 11 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARAMETERS FILE UNIT 11
CLOSE UNIT 11
! Read sequence information
READ SEQUENCE CARDS
* Sequence for glycerol
*
1
GLYC
! Generate the PSF and set up the internal coordinate tables
GENERATE GLYC SETUP
! Generate coordinates from parameters
IC PARAMETERS
IC SEED 1 O1 1 C1 1 C2
IC BUILD
PRINT IC
PRINT COORDINATES
! Write coordinates to disk file
OPEN WRITE UNIT 08 CARD NAME GLYCINI.CRD
WRITE COORDINATES CARD UNIT 08
* Initial coordinates for glycerol generated
* from parameter vales
*
STOP


Using sequence information

 After the topology information has been read into CHARMm, you must specify the sequence of residues that will form the molecular model you are building. Sequence information can be included directly in the input file or in a separate sequence file (.seq) that is read into CHARMm. Sequence information can also be read from CRD or PDB files.

Multiple sets of residue sequences can be specified. Each set of residues is included in a segment definition as part of the PSF. Each molecule in CHARMm is composed of one or more segments, with each segment formed from a sequence of residues. A segment can include a single macromolecular chain, a collection of water molecules, a prosthetic group, or a collection of small organic molecules. Special commands are available to generate sequences of water molecules in segments that allow you to create bulk solvent without specifying each individual water residue.

For information on the maximum number of residues and segments that CHARMm can handle, see CHARMm array dimensions.

Sequence file format

The structure of a sequence file is:

For example, a sequence file to read in the residues for the enkephalin peptide would be written as follows:


* Enkephalin sequence
*
5
TYR GLY GLY PHE MET


Generating and using a protein structure file

 Before any energy calculations are computed for any structure in CHARMm, a molecular model is specified in a protein structure file (PSF). The PSF is the central molecular model data structure in CHARMm. It is the concatenation of information in an RTF, specifying information for the entire structure.

The PSF maintains every bond, bond angle, torsion angle, and improper torsion angle, as well as information needed to generate the hydrogen bonds and the nonbonded list essential for the calculation of the energy of the molecular system. Separate data structures deal with symmetric images of the atoms, parameters, and Cartesian coordinates.

PSFs are created within CHARMm after the sequence information is read and before Cartesian coordinates are generated or read. After a PSF is created, it can be saved as an external file and used repetitively. A PSF can be created automatically in QUANTA or Insight II and read as a single file into CHARMm. Generally, it is unnecessary to make manual modifications to a PSF after it is created.

Organization and structure

PSFs have a hierarchical organization with atoms collected into groups, groups into residues, and residues into segments that comprise the structure. Each atom is uniquely identified within a residue by its IUPAC name, residue identifier, and segment identifier.

The structure of a PSF is as follows:

Atom number

Segment name

Residue identifier

Residue name

Atom name

Atom type

Atomic charge

Atomic mass

A flag to indicate whether the atom is constrained

Creating a PSF

 Two principal methods are available for creating a PSF:

This method translates an existing model into atomic information (charges and types) and internal coordinates. The information is then grouped and written to an external file. This method is the one used by QUANTA or Insight II when it transfers molecular models to CHARMm.

Modifying a PSF

 After a PSF is created from a sequence of residues, it can be modified by using patch residues. Modifications include the addition of a blocking group to the terminal ends of a polymer, changing a neutral functional group to a charged one, adding disulfide bridges, adding links within a single segment or between segments, deleting atoms, and so on.

The format of a patch residue is similar to that of a standard residue. However, patch residues can also add or delete atoms and energy terms from the PSF, and can modify the atom type or charge of a currently defined atom.

Patches differ from residues in another important respect. Because patch residues are applied after a PSF is generated, it is not possible to autogenerate angle or dihedral energy terms within the patch. Therefore, all patch residues that are used after the PSF is generated must contain all necessary angle and dihedral energy terms. The exception is for the patches that are routinely included in the PSF before it is generated. These patches are defined within the RTF as the default first and last patches of a sequence of residues.

Patches can be applied to single or multiple residues within one or more segments. Because unique atom names exist within each residue, no ambiguity arises in applying patches to a single residue. For example, a patch that would protonate an amine group would only need to specify an additional hydrogen atom and the nitrogen to which it is to be bonded. Even though more than one nitrogen can exist in the residue, each has a unique name.

However, with multiple residues, ambiguity often arises. For example, when a disulfide bridge is formed between two cysteine residues, each side chain contains the same atom names. To properly place patches between multiple residues, CHARMm allows you to add an index number before an atom name in a patch residue definition to specify the residue that a particular atom is referencing. This index, using numbers one through nine, corresponds to the order in which the residues are given when the patches are applied.

Example: Adding a patch residue to create a sulfur bridge

A patch residue to add a sulfur bridge between two isoprene units is written as follows:


PRES VULC -0.20	! Name of the patch
ATOM S ST -0.12	! These atoms are either new (S)
ATOM 1C1 CT -0.04	!  or are modified 
ATOM 2Cl CT -0.04
DELE ATOM 1H2C1	! These atoms are deleted from
DELE ATOM 2H2C1	!  the PSF.
BOND S 1C1	! New Bonds, angles and 
BOND S 2C1	!  dihedrals must be defined.
ANGL S 1C1 1C2
ANGL S 1C1 1H1C1
ANGL S 1C1 1-C4
ANGL S 2C1 2C2
DIHE 1H1C1 1C1 S 2C1
DIHE 1C1 S 2C1 2H1C1
IC 1H1C1 1C1 S 2C1   0.00 0.00 180.00 0.00 0.00
IC 1C1 S 2C1 2H1C1   0.00 0.00 180.00 0.00 0.00
IC 1-C4 1C2 *1C1 S   0.00 0.00 240.00 0.00 0.00
IC 2-C4 2C2 *2C1 S   0.00 0.00 240.00 0.00 0.00
IC 1H1C1 S *1C1 1C2  0.00 0.00 120.00 0.00 0.00
IC 1H1C1 S *1C1 1-C4 0.00 0.00 240.00 0.00 0.00
IC 2H1C1 S *2C1 2C2  0.00 0.00 120.00 0.00 0.00
IC 2H1C1 S *2C1 S-C4 0.00 0.00 240.00 0.00 0.00
This script includes the following information:

S is a new atom and will be added.

C1 is a currently existing atom in the isoprene residue. Its charge is to be changed to balance the extra negative charge of the sulfur atom.

Two hydrogens must be deleted so that proper valence is maintained about the C1 atom. All energy terms that include an H2C1 will also be automatically deleted. The 1 and 2 that precede the atom name refer to the first and second residues specified by the PATCH command.

New bonds, angles, and dihedrals that result from the addition of the sulfur atom are defined. Out-of-residue references are indicated by the - or + signs that immediately precede the atom name. New internal coordinates table entries are provided. Many of these new entries seem redundant, but are needed to construct coordinates of the polymer segment across the bridge given the coordinate on one side.


Modeling hydrogens

 CHARMm supports two types of atom models for constructing molecules: the extended atom model and the all-hydrogen (or all-atom) model.

Extended-atom model

For large molecular systems, you can increase modeling efficiency by combining hydrogens with the neighboring heavy atoms to which they are bound. This combined atom is referred to as an extended atom.

For peptide groups and amino acid side chains, the extended-atom dihedral angle potentials can essentially reproduce those obtained from an all-atom representation The extended-atom representation also has been shown to provide a satisfactory representation of the internal vibrations and bulk properties of small molecules including simple peptides.

Although large molecular systems are most efficiently modeled in the extended-atom mode, keep in mind the advantages and disadvantages of this representation. These are outlined below.

Advantages

Disadvantages

In a modified form of the extended-atom model called the polar hydrogen model, only aliphatic hydrogens that are not significantly charged and cannot participate in hydrogen bonds are presented as extended atoms. The polar hydrogen model allows hydrogen bonding and dipole-dipole interactions to be effectively simulated, eliminating several of the disadvantages of the extended-atom model.

All-hydrogen model

The all-hydrogen (or all-atom) model contains all atoms, including both polar and non-polar hydrogens. For small molecules, this is the model of choice.


Generating and using coordinates

A complete set of Cartesian coordinates must be generated prior to any energy computation. Situations often exist where all coordinates are not available or where some or all of a structure must be modified or built.

CHARMm has generalized facilities for the construction, manipulation, transformation, and interconversion of internal coordinates (ICs) and Cartesian coordinates.

Although crystallographic and simulation results are usually obtained and stored in Cartesian coordinates, it is frequently desirable to transform them either to re-oriented Cartesian frames or to internal coordinates. All of the reorientations and transformations can be applied to an easily specified set of atoms. Because a flexible atom selection facility is provided in CHARMm, internal coordinates can be specific without regard to the bonding between atoms.

CHARMm internal coordinates

Given the positions of any three atoms, the position of a fourth atom can be defined. In CHARMm, this is how internal coordinates are developed. The position of the fourth atom is determined relative to the positions of the three defined atoms using a bond distance, bond angle, and dihedral specification.

In CHARMm, the internal coordinates definition for a chain of atoms takes the form:

IC I J K L Rij ijk ijkl jkl Rkl

The data structure uses five values to define four atoms. In this way, if either end-point is unknown, but three atoms are determined, the position of the fourth atom can be found.

The improper type of internal coordinates data structure is used for branching structures. The definition for a branched chain of atoms takes the form:

IC I J *K L Rij ijk ijkl jkl Rkl

The asterisk preceding the Kth atom indicates that this is an improper dihedral, and K is the central atom.

Because five values exist in the data structure for every atom, the positions are over-specified. When internal coordinates are used to determine atomic positions, CHARMm uses the first acceptable value to build a structure and ignores any redundancies.

The internal coordinates used to build a model can be specified arbitrarily, taken from an existing structure, or chosen from minimum energy values in parameter tables. For example, specified internal coordinates could be used to construct coordinates for a peptide backbone when only values of the phi and psi dihedral angles are available.

The internal coordinate set can be manipulated so that you can build or modify a coordinate set in terms of bond angles and distances. The ability to take a Cartesian coordinate structure, fill and edit the internal coordinates table to specify modifications in bond distances, angles, or dihedral values, and then construct a new structure is extremely useful. With this feature, entire sections of a molecular model can be moved relative to others in an arbitrary but easily specified manner.

Building Cartesian coordinates from internal coordinates

X-ray crystal coordinates are sometimes ill-defined or incomplete. In these cases, a method to compute the Cartesian coordinates of any atom in the molecular system is a requirement. Additionally, coordinates need to be easily built for some modeling applications that require modifications to be made in structures.

In CHARMm, internal coordinates of any structure with some known coordinate positions can be used to construct Cartesian coordinates for atoms whose positions are unknown. Given the positions of any three atoms, the position of a fourth atom can be defined in relative terms using three values: a distance, an angle, and a dihedral specification.

The following commands can be used to generate missing coordinates if the internal coordinates tables have been set up during PSF generation:

Example: Constructing a model using generated internal coordinates

 This example illustrates the standard steps for constructing a model. The first part of the script generates a PSF. An RTF file is read, the sequence of enkephalin is entered, and a segment named ENKP is defined. Generating the segment creates a PSF within CHARMm.

The second part of the script generates a set of Cartesian coordinates for enkephalin using internal coordinates. When three atoms have coordinates, the remaining atoms can be placed by using the relationships defined by the PSF.

This example uses or generates the following files:

AMINO.RTF

PARM.PRM 


*...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
* 
! Open and read an all hydrogen amino acid topology file
OPEN READ UNIT 11 FILE NAME "$CHM_DATA/AMINOH.BIN"
READ RTF UNIT 11 FILE
CLOSE UNIT 11
! Read the sequence for met-enkephalin. Following the
! title, the number of residues in the sequence is 
! specified. The next line(s) contain the residue
! sequence.
READ SEQUENCE CARD
* Pentapeptide sequence for met-enkephalin
*
 5
TYR GLY GLY PHE MET
! Open and read the parameter file
OPEN READ UNIT 12 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARAMETERS UNIT 12 FILE
CLOSE UNIT 12
! Generate the PSF with a segment identifier of ENKP.
! Set up the internal coordinate tables.
GENERATE ENKP SETUP
! Construct initial Cartesian coordinates from the
! parameter values. Initially place (seed) the first
! three atoms of the peptide sequence.
IC PARAMETERS
IC SEED 1 N 1 CA 1 C
IC BUILD
! Print the coordinates in the output file.
PRINT COORDINATES
! Use the WRITE command to output coordinates to an external 
! file. All WRITE commands must be preceded with an OPEN 
! command.
OPEN WRITE UNIT 08 CARD NAME ENKPINI.CRD
WRITE COORDINATES UNIT 08 CARD
* Initial coordinates for met-enkephalin created by
* using parameter data and the internal coordinate
* tables.
*
! Exit CHARMm
STOP
This example only saves the coordinates, but it can be modified to also save the PSF by adding the following lines to the script before exiting CHARMm:


WRITE PSF CARD NAME ENKP.PSF
*PSF for enkephalin
*
The saved PSF can be used in other CHARMm scripts.

Example: Converting a PDB file to a PSF

This example presents a situation encountered frequently by CHARMm users. A set of coordinates is available for a structure from the Brookhaven Protein Database (PDB). The data must be converted and additional data generated before any CHARMm calculations can be performed.

A PDB file contains sequence information, atom names, and Cartesian coordinates for a molecule. Although proper PDB files contain connect records that describe the connectivity, CHARMm does not use these. An RTF for amino acids must be supplied by CHARMm. The connectivity, atom types, and charges are defined by the RTF.

PDB files do not contain hydrogens. Hydrogens are created in CHARMm and stored in the PSF. Generally, CHARMm uses polar hydrogens defined in the RTF file for this task. Cartesian coordinates for these atoms are also generated by CHARMm.

In the following script, the RTF and parameters are first read into CHARMm. A PDB file is opened and the sequence is read directly from the file. This provides sufficient information to generate the PSF. CHARMm then returns to the beginning of the PDB file (REWIND) and reads the Cartesian coordinates from it.

This example uses or generates the following files:

AMINO.BIN

PARM.BIN

1gcn.pdb 


*...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
*
! Open and read polar hydrogen amino acid topology file(binary)
OPEN READ UNIT 11 FILE NAME "$CHM_DATA/AMINO.BIN"
READ RTF UNIT 11 FILE
CLOSE UNIT 11
! Open and read parameter file (binary)
OPEN READ UNIT 12 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARAMETERS UNIT 12 FILE
CLOSE UNIT 12
! Read the sequence from the protein data bank file
OPEN READ UNIT 13 CARD NAME 1GCN.PDB
READ SEQUENCE PDB UNIT 13
! Generate segment 1GCN and set up internal coordinate table
GENERATE 1GCN SETUP
! Rewind the protein data bank file for reading of coordinates REWIND UNIT 13
! Read glucagon PDB coordinates
READ COORDINATES PDB UNIT 13
CLOSE UNIT 13
! Print coordinates to output file. Missing coordinates
! (hydrogens) are assigned values of 9999.0
PRINT COORDINATES 
! Construct hydrogen positions
HBUILD
! Print coordinates to output file. All hydrogen
! positions constructed
PRINT COORDINATES
! If necessary, any additional missing coordinates
! can be placed using the following commands
IC FILL PRESERVE
IC PARAMETERS
IC BUILD
! Write the coordinates to a disk file
OPEN WRITE UNIT 14 CARD NAME GLUCINI.CRD
WRITE COORDINATES CARD UNIT 14
* Glucagon coordinates
* Heavy atom positions taken from PDB file
* Hydrogen positions generated
*
STOP

Coordinate arrays

 CHARMm maintains two parallel internal arrays for coordinates: the main coordinate array and the comparison coordinate array.

In the name of the second array, the word comparison is historical. This array was used almost exclusively for comparison purposes in early versions of CHARMm. CHARMm now utilizes this array for a variety of functions related to coordinate operations. For example, atomic forces (direction and magnitude) can be stored in the comparison coordinate array. Examples in this book use the comparison coordinate array for temporary storage or for analysis of coordinates.

Most CHARMm commands operating on the main coordinate array can also operate on the comparison coordinate array by simply adding the keyword COMP to the command line. This allows great flexibility in studying multiple conformations of structures or in comparing structures at different points of a trajectory. Most often, coordinates are copied to the comparison array with the line: 


COOR COPY COMP

Format of the CHARMm coordinate file

 CHARMm coordinate files contain information about the location of each atom in Cartesian (3D) space. The format of the ASCII (CARD) CHARMm coordinate files is:

The coordinate lines contain specific information about each atom in the model and consist of the following fields:

The FORTRAN FORMAT statement for the coordinate lines is:


I5, I5, 1X, A4, 1X, A4, 3(F10.5), 1X, A4, 1X, A4, F10.5

Example: Coordinate file for a two-segment model

 The following example illustrates a coordinate file for a multi-segment model with two segments (SEG1 and PEP2):


* Two dipeptides, each in a separate segment
*
36
 1    1 ALA  N      0.00000   0.00000   0.00000 SEG1 1 0.00000
 2    1 ALA  HT1   -0.34665  -0.98053   0.00000 SEG1 1 0.00000
 3    1 ALA  HT2   -0.34665   0.49026   0.84916 SEG1 1 0.00000
 4    1 ALA  HT3   -0.34665   0.49026  -0.84916 SEG1 1 0.00000
 5    1 ALA  CA     1.47000   0.00000   0.00000 SEG1 1 0.00000
 6    1 ALA  CB     1.95113  -0.74250  -1.24824 SEG1 1 0.00000
 7    1 ALA  C      2.02771   1.40861   0.00000 SEG1 1 0.00000
 8    1 ALA  O      1.31114   2.40216   0.00000 SEG1 1 0.00000
 9    2 ALA  N      3.37163   1.46259   0.00000 SEG1 2 0.00000
10    2 ALA  H      3.92937   0.63258   0.00000 SEG1 2 0.00000
11    2 ALA  CB     3.51425   3.52492   1.24824 SEG1 2 0.00000
12    2 ALA  C      5.50913   2.71647   0.00000 SEG1 2 0.00000
13    2 ALA  OCT1   6.04832   1.61095   0.00000 SEG1 2 0.00000
14    2 ALA  OCT2   6.14262   3.77078   0.00000 SEG1 2 0.00000
15    2 ALA  CA     3.99557   2.78255   0.00000 SEG1 2 0.00000
16    3 THR  N      0.00000   0.00000   0.00000 PEP2 1 0.00000
17    3 THR  HT1   -0.34665  -0.98053   0.00000 PEP2 1 0.00000
18    3 THR  HT2   -0.34665   0.49026   0.84916 PEP2 1 0.00000
19    3 THR  HT3   -0.34665   0.49026  -0.84916 PEP2 1 0.00000
20    3 THR  CA     1.47000   0.00000   0.00000 PEP2 1 0.00000
21    3 THR  CB     1.95287  -0.73967  -1.24925 PEP2 1 0.00000
22    3 THR  OG1    3.38262  -0.75331  -1.27229 PEP2 1 0.00000
23    3 THR  HG1    3.63799  -1.21844  -2.05787 PEP2 1 0.00000
24    3 THR  CG2    1.43958  -0.04232  -2.51065 PEP2 1 0.00000
25    3 THR  C      2.02771   1.40861   0.00000 PEP2 1 0.00000
26    3 THR  O      1.31114   2.40216   0.00000 PEP2 1 0.00000
27    4 THR  N      3.37163   1.46259   0.00000 PEP2 2 0.00000
28    4 THR  H      3.92937   0.63258   0.00000 PEP2 2 0.00000
29    4 THR  CB     3.51754   3.52536   1.24925 PEP2 2 0.00000
30    4 THR  OG1    4.10882   4.82690   1.28535 PEP2 2 0.00000
31    4 THR  HG1    3.78782   5.25025   2.07049 PEP2 2 0.00000
32    4 THR  CG2    3.92118   2.75819   2.50997 PEP2 2 0.00000
33    4 THR  C      5.50913   2.71647   0.00000 PEP2 2 0.00000
34    4 THR  OCT1   6.04832   1.61095   0.00000 PEP2 2 0.00000
35    4 THR  OCT2   6.14262   3.77078   0.00000 PEP2 2 0.00000 
36    4 THR  CA     3.99557   2.78255   0.00000 PEP2 2 0.00000

Constructing hydrogen positions

 Coordinates that are provided from experimental data normally give positions only for heavy, non-hydrogen atoms. However, CHARMm requires the coordinates for all atoms given in the topological definitions of the individual residues in order to compute the empirical energy functions.

If necessary, CHARMm builds and optimizes all undefined hydrogen coordinates for the molecule as well as for the surrounding water. In the first part of hydrogen building, all proton positions of the molecule are constructed. The protons are initially classified according to their respective environments. These classifications are defined in the following list. Note that the donor atom is the heavy atom to which the proton is connected.

Class 1
A proton bound to a donor with at least two heavy atom donor antecedents

Class 2
A proton bound to a donor with one heavy atom donor antecedent and no other proton

Class 3
A proton bound to a donor with one heavy atom antecedent and one other proton

Class 4
A proton bound to a donor with one heavy atom antecedent and two other protons

Proton positions are generated using the following steps:

1.   All protons of class 1 are placed using equilibrium bond lengths, angles, and dihedrals. These protons are over-determined if there is more than one heavy atom donor antecedent. In these cases, all possible ways to place the proton are averaged to give one position.

2.   Protons of all other classes (2-4) are constructed. For these classes, there is an additional degree of freedom for placing the protons. To find an optimum position, the dihedral angle with the symmetry axis antecedent-donor is modified in small steps (using the keyword PHIStp) over a certain range determined by the symmetry of the donor group:

    a.   For each dihedral angle, the protons of the donor are placed according to their equilibrium geometry.

    b.   Relative energy of the corresponding configuration is evaluated. The energy is determined by using the hydrogen bond potential, the van der Waals term, the electrostatic term, and the dihedral term derived from the full energy expression of CHARMm.

    c.   The dihedral angle that yields the lowest energy is chosen and the protons of the donor group are placed using this optimum dihedral angle.

In this step, this procedure is performed in the order given by the residue sequence of the molecule. Protons that have not yet been constructed have no influence on the current energy evaluations.

3.   After construction of all explicit protons, water protons are constructed:

    a.   A sequence of water molecules is determined independent of any other sequence.

    b.   Waters are ordered with respect to the minimum distance of the water oxygen to any atoms in the molecule.

    c.   Protons of waters near the molecule are constructed first. Three classes of water molecules are handled by this method:

    Water able to form two different hydrogen bonds to an acceptor atom -- Protons are placed by performing a rotation of the water molecule in the plane defined by the two best hydrogen bonds, taking the minimum energy configuration.

    Water able to form only one hydrogen bond to acceptor atoms -- One proton is placed on the (linear) hydrogen bond and the other proton is placed using the equilibrium geometry. The minimum energy configuration is taken. The evaluated relative energy is the sum of the van der Waals, the electrostatic, and the hydrogen bonding terms.

    Water forming no hydrogen bonds -- Protons are placed in a standard manner (that is, H1 on the X axis, and H2 in the XY plane) after all other protons have been placed.


Using parameter files

 Parameter files contain values for force constants, equilibrium geometries, van der Waals radii, stretching constants, and other constants necessary to calculate the energy of a molecule. The parameter files allow CHARMm to impart actual physical constraints upon any structure.

File structure

A parameter file must start with a valid CHARMm title. The remaining information is read in free-field format as commands to read in parameter data. Each section of parameters is initiated by a single CHARMm parameter keyword. The structure of a parameter file is:

No wildcard usage is allowed for bonds or angles.

For dihedrals, two parameter specifications can be used:

Or

For impropers, five different parameter specifications can be used:

When an improper dihedral is classified, the first acceptable match (from the above order) is chosen. The match can be made in either direction.

Wildcard atom-types can be used in specifying NBOND, NBFIX, or HBOND interaction parameters.

The periodicity value for dihedrals and improper dihedral terms must be an integer. If it is positive, then a cosine functional form is used. Only positive values 1, 2, 3, 4, 5, and 6 are allowed. Phase is either 0.0 or 180.0 for dihedrals with the minimum staggered or eclipsed, respectively.

When periodicity is given as 0.0, a harmonic restoring potential phi-phi_min is used. The phase value gives phi_min for this option. This functional form can be used for either dihedrals or impropers. For dihedrals, the selection is usually based on the middle two atoms. For impropers, the selection is usually based on the outer two atoms. When periodicity is given as 0.0 for other than the first dihedral in a multiple dihedral set, amplitude is a constant added to the energy. This is needed for the Ryckaert-Bellemans potential of hydrocarbons.

Nonbonded parameters are determined based on how values are listed:

Nonbonded parameters for 1-4 interactions can be specified by placing the extra set of parameters after the first. By default, the same parameters are used for 1-4 interactions as for all other interactions.

The NBFIX set of parameters allows individual atom type VDW pair interactions to be specified. These parameters are processed in order. In the case of duplicate specifications, the last specification is used.

Ryckaert-Bellemans torsional potential

To calculate the Ryckaert-Bellemans torsional potential (see Ryckaert and Bellemans 1975 and 1978) for butane and other extended-atom hydrocarbons, the following terms should be included in the parameter file:


PHI
CH3E CH2E CH2E CH3E 0.470467 5 0.0
CH3E CH2E CH2E CH3E 0.783947 4 0.0
CH3E CH2E CH2E CH3E 2.53516 3 0.0
CH3E CH2E CH2E CH3E 1.56789 2 0.0
CH3E CH2E CH2E CH3E 2.34787 1 0.0
CH3E CH2E CH2E CH3E -4.70368 0 0.0
This potential should be used with SHAKE (described in Chapter 5, Setting Constraints and Periodic Boundaries). A zero periodicity term should not be the first in the set or the term is treated as an improper torsion.

Example: Parameter file PARM.PRM

For an example of a parameter file, see the file $CHM DATA/ PARM.PRM in your CHARMm software.


Summary examples

The following scripts contain typical examples of model preparation. These examples cover generation and use of residue topology files, PSFs, internal and Cartesian coordinates, parameter files, and other information that can be used for preparing models for energy calculations and analysis.

Example: Constructing polymer segments joined by a sulfur bridge

This script illustrates how to use a patch residue to modify a PSF. Two polymer segments are created using the ISOP residue, then a sulfur bridge is created to connect them. The sulfur bridge uses the patch residue VULC.

This example uses or generates the following files:


*...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
* ...
* This CHARMm command file builds two short polymer segments
* of natural rubber and joins them with a sulphur bridge.
* ...
* 
READ RTF CARD

* Topology file containing the isoprene residue, terminal 
* patches, and the sulphur bridge patch
*
   22    2		   ! CHARMm Version #
! MASS statements deleted for brevity
! Out of residue declarations
DECL -C4
DECL +C1
DECL +C2
! Autogenerate all angles during generation of the PSF AUTOGENERATE ANGLES DIHEDRALS
! Apply terminal patches TER1 and TER2 by default during
! generation of the PSF
DEFA FIRS TER1 LAST TER2
! The isoprene residue
RESI ISOP        0.00
ATOM C1   CT    -0.20
ATOM H1C1 HA     0.10
ATOM H2C1 HA     0.10
ATOM C2   CUA1   0.00
ATOM C3   CUA1  -0.10
ATOM H1C3 HA     0.10
ATOM C4   CT    -0.20
ATOM H1C4 HA     0.10
ATOM H2C4 HA     0.10
ATOM C5   CT    -0.30
ATOM H1C5 HA     0.10
ATOM H2C5 HA     0.10
ATOM H3C5 HA     0.10
BOND C1 C2   C2 C3   C3 C4   C2 C5   C4 +C1
BOND C1 H1C1  C1 H2C1
BOND C3 H1C3
BOND C4 H1C4   C4 H2C4
BOND C5 H1C5   C5 H2C5   C5 H3C5
IC -C4 C1 C2 C3        0.00  0.00  180.00  0.00  0.00
IC  C1 C2 C3 C4        0.00  0.00    0.00  0.00  0.00
IC  C2 C3 C4 +C1       0.00  0.00  180.00  0.00  0.00
IC  C3 C4 +C1 +C2      0.00  0.00  180.00  0.00  0.00
IC -C4 C2 *C1 H1C1     0.00  0.00  120.00  0.00  0.00
IC -C4 C2 *C1 H2C1     0.00  0.00  240.00  0.00  0.00
IC  C1 C3 *C2 C5       0.00  0.00  180.00  0.00  0.00
IC  C2 C4 *C3 H1C3     0.00  0.00  180.00  0.00  0.00
IC  C3 +C1 *C4 H1C4    0.00  0.00  120.00  0.00  0.00
IC  C3 +C1 *C4 H2C4    0.00  0.00  240.00  0.00  0.00
IC  C1 C2 C5 H1C5      0.00  0.00  180.00  0.00  0.00
IC  C2 H1C5 *C5 H2C5   0.00  0.00  120.00  0.00  0.00
IC  C2 H1C5 *C5 H3C5   0.00  0.00  240.00  0.00  0.00
! The terminal patches
PRES TER1        -0.20
ATOM H3C1 HA      0.10
ATOM C1   CT     -0.30
BOND H3C1 C1
DIHE H3C1 C1 C2 C3
IC  H3C1 C1 C2 C3       0.00  0.00  180.00  0.00  0.00
IC  H3C1 C2 *C1 H1C1    0.00  0.00  120.00  0.00  0.00
IC  H3C1 C2 *C1 H2C1    0.00  0.00  240.00  0.00  0.00
PRES TER2        -0.20
ATOM H3C4 HA      0.10
ATOM C4   CT     -0.30
BOND H3C4 C4
DIHE H3C4 C4 C3 C2
IC  H3C4 C4 C3 C2       0.00  0.00  180.00  0.00  0.00
IC  H3C4 C3 *C4 H1C4    0.00  0.00  120.00  0.00  0.00
IC  H3C4 C3 *C4 H2C4    0.00  0.00  240.00  0.00  0.00
! The sulphur bridge patch
PRES VULC     -0.20
ATOM S    ST  -0.12
ATOM 1C1  CT  -0.04
ATOM 2C1  CT  -0.04
DELE ATOM 1H2C1
DELE ATOM 2H2C1
BOND S 1C1   S 2C1
ANGL S 1C1 1C2   S 1C1 1H1C1  S 1C1 1-C4   S 2C1 2C2
ANGL S 2C1 2H1C1   S 2C1 2-C4
DIHE 1H1C1 1C1 S 2C1   1C1 S 2C1 2H1C1
IC 1H1C1 1C1 S 2C1    0.00  0.00  180.00  0.00  0.00
IC 1C1 S 2C1 2H1C1    0.00  0.00  180.00  0.00  0.00
IC 1-C4 1C2 *1C1 S    0.00  0.00  240.00  0.00  0.00
IC 2-C4 2C2 *2C1 S    0.00  0.00  240.00  0.00  0.00
IC 1H1C1 S *1C1 1C2   0.00  0.00  120.00  0.00  0.00
IC 1H1C1 S *1C1 1-C4  0.00  0.00  240.00  0.00  0.00
IC 2H1C1 S *2C1 2C2   0.00  0.00  120.00  0.00  0.00
IC 2H1C1 S *2C1 2-C4  0.00  0.00  240.00  0.00  0.00
END
! Read in the parameter file
OPEN READ UNIT 12 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARA UNIT 12 FILE
CLOSE UNIT 12
! Generate the first segment (3 isoprene residues)
READ SEQU CARD
* TEST
*
3
ISOP ISOP ISOP
GENE RBR1 SETUP
! Generate the second segment (3 isoprene residues)
READ SEQUENCE CARD
* TEST
*
3
ISOP ISOP ISOP
GENERATE RBR2 SETUP
! Apply the patch to bridge the two segments with a sulphur 
! atom
PATCH VULC RBR1 2 RBR2 2 SETUP
! Initialize the internal coordinate arrays and build
! the coordinates
IC PURGE
IC PARAMETERS
IC SEED 1 C1 1 C2 1 C3
IC BUILD
IC PRINT
! Merge the segments into one segment
JOIN RBR1 RBR2 RENU
! Rename the segment newly merged segment RBBR
RENAME SEGI RBBR SELE SEGI RBR1 END
! Write the initial coordinates to disk
OPEN WRITE UNIT 11 CARD NAME RUBBER.CRD
WRITE COORDINATES UNIT 11 CARD
* TEST
*
! The following commands are not referenced in this
! example
UPDATE INBFRQ -1 IHBFRQ 0 RDIE
MINI ABNR NSTEP 1000 NPRINT 50 TOLENR 0.000001 TOLGRD 0.000001
STOP

Example: Constructing an alpha helix of polyalanine

 This script builds an alpha helix of polyalanine. The script is quite simple, creating a polyalanine sequence and PSF, then building coordinates. The key lines in this script are in the four lines:


SET 1 1
SET 2 11
OPEN READ UNIT 18 CARD NAME HELIX.STR
STREAM UNIT 18
The first two lines of this group assign values to the variables 1 and 2. Variable 1 is set to 1 and variable 2 is set to 11. The variables can be used as counters for loops in scripts, and they can also be substituted as values using the @ sign. So the command:


SET 3 @1
sets the variable 3 to the value of variable 1.

The next line opens an external file. The STREAM command instructs CHARMm to read every line in the file as a CHARMm command.

The second part of this example is the file HELIX.STR. This file shows some of the flow control elements that can be used in a CHARMm script including labels, if statements, and goto statements. In this script, the dihedral angles of each successive residue are set to an alpha-helical conformation.

This example uses or generates the following files:

AMINO.BIN

PARM.BIN

HELIX.STR

Script part 1 


*...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
! Open and read the topology and parameter files (binary)
OPEN READ UNIT 11 FILE NAME "$CHM_DATA/AMINO.BIN"
READ RTF UNIT 11
CLOSE UNIT 11
OPEN READ UNIT 11 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARA UNIT 11 FILE
CLOSE UNIT 11
! Read the sequence for polyalanine
READ SEQU ALA 11
* Polyalanine (11 residues)
*
GENERATE HELX SETUP
! Generate coordinates from parameter data; first fill the
! internal coordinate
IC PARAMETERS
! Define the range of residues to be defined as an alpha
! helix; set CHARMm parameter 1 to be the first residue and
! CHARMm parameter 2 to be the last residue
SET 1 1
SET 2 11
OPEN READ UNIT 18 CARD NAME HELIX.STR
STREAM UNIT 18
! Initially place three coordinates, and build the rest
IC SEED 1 N 1 CA 1 C
IC BUILD
! Write coordinates to disk file
OPEN WRITE UNIT 04 CARD NAME ALPHAHLX.CRD
WRITE COORDINATES CARDS UNIT 04
* Initial helix coordinates constructed using IC EDIT
* and IC BUILD
*
STOP
Script part 2 


*...
* Copyright (c) 1994 
* Molecular Simulations Inc.
* All Rights Reserved
*...
* HELIX.STR
*...
* This stream file edits the internal coordinate table
* by defining the phi and psi angles to be an alpha helix
* for a range of residues (CHARMm parameters 1 and 2)
*
! Define variables for residues i+1 and i+2
SET 3 @1
SET 4 @1
INCR 3 BY 1
INCR 4 BY 2
LABEL START
! Invoke IC EDIT mode, and define phi and psi
! dihedral angles for an alpha helix
IC EDIT
DIHE @1 C @3 N @3 CA @3 C -57.0
DIHE @3 N @3 CA @3 C @4 N -47.0
END
! Increment counters and check for last specified residue
INCR 1 BY 1
INCR 3 BY 1
INCR 4 BY 1
IF 4 GT @2 GOTO STOP
GOTO START
LABEL STOP
RETURN

Example: Generating random conformations of enkephalin

 Another example of how to use stream files is provided below. Here, 100 random conformations are generated and their energies computed. The main script is used to set up the PSF while the generation of the conformations and energies is performed in the stream file. This generic stream file could easily be inserted into other input files or customized for other purposes.

This example uses or generates the following files:

ENKP.ICR -- Internal coordinates table

ENKPCONF000.CRD through ENKPCONFxxx.CRD -- Conformer Cartesian coordinate files

ENKPCONFENG.LST

AMINO.BIN

PARM.BIN

RANDCONF.STR

Script part 1 


*...
* Copyright (c) 1994 
* Molecular Simulations Inc.
* All Rights Reserved
! Open and read the topology and parameter files
OPEN READ UNIT 11 FILE NAME "$CHM_DATA/AMINOH.BIN"
READ RTF UNIT 11
CLOSE UNIT 11
OPEN READ UNIT 12 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARA UNIT 12
CLOSE UNIT 12
! Generate the PSF and initial coordinates
READ SEQU CARD
* Enkephalin
*
5
TYR GLY GLY PHE MET
GENE MAIN SETU
IC PARA
IC SEED 1 N 1 CA 1 C
IC BUIL
! Write the internal coordinate table to an external file
OPEN WRIT UNIT 15 CARD NAME ENKP.ICR
IC WRIT UNIT 15 CARD
* Enkephalin initial internal coordinate table
*
! Compute an initial energy; write out initial
! coordinates and open a log file for conformation
! energies
UPDATE
GETE
OPEN WRIT UNIT 24 CARD NAME ENKPCONF000.CRD
WRIT COOR UNIT 24 CARD
* Enkephalin initial conformation with energy?ener
*
OPEN WRIT UNIT 25 CARD NAME ENKPCONFENG.LST
WRIT TITL UNIT 25
* Enkephalin initial conformation energy = ?ener
*
! Stream random conformation generation file
OPEN READ UNIT 28 CARD NAME RANDCONF.STR
STREAM UNIT 28
CLOSE UNIT 24
CLOSE UNIT 28
STOP
Script part 2 


* ...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
* ...
* This stream file randomly generates values for the dihedral 
* angles phi and psi in all residues of a protein
* ...
*
! Set counter to 0
SET 1 0
! Main loop: increment and check counter against the maximum
! number of conformations desired
LABEL START
INCR 1 BY 1
IF 1 GT 100 GOTO STOP
! Temporarily save internal coordinate table; keep only those
! dihedrals that correspond to phi and psi and assign random
! values to these internal coordinates; print results to
! output file
IC SAVE
IC KEEP FIRST SELE ATOM * * N .OR. ATOM * * CA END
IC KEEP SECO SELE ATOM * * CA .OR. ATOM * * C END
IC KEEP THIRD SELE ATOM * * C .OR. ATOM * * N END
IC KEEP FOURTH SELE ATOM * * N .OR. ATOM * * CA END
IC RAND 
IC PRIN
! Restore full internal coordinate table preserving randomly
! assigned values of phi and psi; initialize and rebuild all
! Cartesian coordinates
IC REST PRES
COOR INIT
IC SEED 1 N 1 CA 1 C
IC BUIL
! Compute an energy and write to log file
GETE
WRITE TITLE UNIT 25
* Enkephalin conformation @1 energy = ?ener
*
! Check value of energy; if less than 0.0, write
! coordinates to unique file; if not, generate another
! conformation
SET 2 ?ENER
IF 2 GT 0.0 GOTO START
SET 3 ENKPCONF
IF 1 LT 100 SET 3 ENKPCONF0
IF 1 LT 10 SET 3 ENKPCONF00
OPEN WRITE UNIT 44 CARD NAME @3@1.CRD
WRIT COOR UNIT 44 CARD
* Enkephalin conformation @1 with energy ?ener
*
GOTO START
LABEL STOP
RETURN

Example: Constructing an N-methylacetamide dimer

 The following input file uses many of the coordinate manipulation commands to generate a copy of the N-methylacetamide molecule and position it so that the energetics of the dimer can be studied.

This example uses or generates the following files: 


*...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
* ...
* This CHARMm command file uses the coordinate manipulation
* tools of CHARMm to generate and model the N-methylacetamide
* dimer
* ...
* 
READ RTF CARDS
* RTF for N-methyl acetamide
*
   22    2
MASS  1 H   1.00800 ! Hydrogen bonding hydrogen (neutral group)
MASS 13 CH3E 15.03500 ! Extended atom C with three hydrogens
MASS 14 C  12.01100 ! Carbonyl or Guanidinium carbon
MASS 32 NP 14.00670 ! Peptide/amide nitrogen
MASS 40 O  15.99940 ! Carbonyl oxygen for acids/amides
AUTOGENERATE ANGLES DIHEDRALS
RESI NMA           0.00
GROU
 ATOM CL   CH3E    0.00
GROU
 ATOM C    C       0.55
 ATOM O    O      -0.55
GROU
 ATOM N    NP     -0.35
 ATOM H    H       0.25
 ATOM CA   CH3E    0.10
BOND N CA     N H
BOND CL C     C N      C O
IMPH C CL N O
IMPH N C CA H
IC   N   CL   *C   O   0.0000   0.00  180.00    0.00   0.0000
IC   CL   C    N   CA  0.0000   0.00  180.00    0.00   0.0000
IC   C   CA   *N   H   0.0000   0.00  180.00    0.00   0.0000
PATCH FIRST NONE LAST NONE
END
! Open and read the parameter file (binary)
OPEN READ UNIT 11 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARA UNIT 11 FILE
CLOSE UNIT 11
! Read the sequence information
READ SEQUENCE CARDS
* SEQUENCE FOR N-METHYL ACETAMIDE
*
1
NMA
! Generate the first segment (A); generate the second
! segment (B) by duplicating the first
GENERATE A SETUP
GENERATE B SETUP DUPLICATE A
! Construct coordinates for the first segment
IC PARAMETERS
IC SEED  1 C  1 N  1 CA
IC BUILD
PRINT COORDINATES
! Orient the amide nitrogen and hydrogen of segment a
! along the x-axis
COOR ORIE SELE ATOM A 1 N .OR. ATOM A 1 H END
! Duplicate the coordinates of segment A to segment B
COOR DUPLICATE SELE ATOM A * * END SELE ATOM B * * END
! Temporarily save the coordinates in the comparison set
COOR COPY COMP
! Orient the coordinates such that the carbonyl oxygen
! and carbon of segment B is aligned with the x-axis;
! this will also reorient the coordinates of segment A
COOR ORIE SELE ATOM B 1 C .OR. ATOM B 1 O END
! Retrieve coordinates of segment A from the comparison set
COOR COPY SELE SEGI A END
! Define a vector between the amide hydrogen of
! segment A and the carbonyl oxygen of segment B;
! this command also determines the distance between
! these atoms
COOR AXIS SELE ATOM A 1 H END SELE ATOM B 1 O END
! Translate the coordinates of segment B along the x-axis
! by the distance determined from the COOR AXIS command;
! this will place the hydrogen of segment A and the oxygen
! of segment B in the same place
COOR TRANS SELE SEGI B END AXIS FACT -1.0
! Translate the coordinates of segment B another 1.8 A,
! the optimum hydrogen bond distance
COOR TRANS SELE SEGI B END XDIR 1.8
! Write the initial coordinates to disk
OPEN WRITE CARD UNIT 12 NAME NMAINI.CRD
WRITE COORDINATES CARD UNIT 12
* NMA INITIAL COORDINATES
*
STOP

Example: Constructing conformations of cyclohexane

 This example uses the CHARMm internal coordinates commands to construct the chair and boat conformations of cyclohexane. Note that Cartesian coordinates are generated directly from idealized torsion angles.

The following files are used or generated in this script:

CHAIR.CRD

CYCLO.TXT

BOAT.CRD 


*...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
* ...
* This CHARMm command file uses the internal coordinate (IC) 
* commands of CHARMm to generate the chair and boat 
* conformers of cyclohexane
* ...
*
! Read the RTF topology file for cyclohexane
READ RTF CARDS
* RTF for cyclohexane
*
   22    2                ! Version number
! MASS Statements removed for brevity
AUTOGENERATE ANGLES DIHEDRALS
RESIDUE CYCL     0.00000
GROUP
   ATOM C1   CH2E      0.00
   ATOM C2   CH2E      0.00
   ATOM C3   CH2E      0.00
   ATOM C4   CH2E      0.00
   ATOM C5   CH2E      0.00
   ATOM C6   CH2E      0.00
BOND   C1 C2   C2 C3   C3 C4   C4 C5   C5 C6   C6 C1
IC  C1 C2 C3 C4    0.0000    0.00   +60.0    0.00   0.0000
IC  C2 C3 C4 C5    0.0000    0.00   -60.0    0.00   0.0000
IC  C3 C4 C5 C6    0.0000    0.00   +60.0    0.00   0.0000
IC  C4 C5 C6 C1    0.0000    0.00   -60.0    0.00   0.0000
IC  C5 C6 C1 C2    0.0000    0.00   +60.0    0.00   0.0000
IC  C6 C1 C2 C3    0.0000    0.00   -60.0    0.00   0.0000
PATC FIRST NONE LAST NONE
END
! Open and read the parameter file (binary)
OPEN READ UNIT 11 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARAMETERS UNIT 11 FILE
CLOSE UNIT 11
! Read the sequence information and generate the PSF
READ SEQUENCE CARDS
* SEQUENCE FOR CYCLOHEXANE
*
1
CYCL
GENERATE A SETUP
! Generate coordinates for the chair conformer
! (defined by the IC commands in the RTF)
IC PARAMETERS
PRINT IC
IC SEED  1 C1  1 C2  1 C3
IC BUILD
OPEN WRITE CARD NAME CHAIR.CRD UNIT 13
WRITE COORDINATES CARD UNIT 13
* Chair conformation coordinates for cyclohexane
*
! Compute the static energy
NBONDS CUTNB 99.0 CDIE VSWITCH SWITCH
ENERGY
! Use the ?ener variable to output text to a file
OPEN WRITE UNIT 33 CARD NAME CYCLO.TXT
WRITE TITLE UNIT 33
* Energy for chair cyclohexane is ?ener kcal
*
! Close contact searches
COOR DIST ENERGY CLOSE 14EXCLUSIONS NONBONDS
! Generate the boat conformer by editing the internal
! coordinate tables
IC EDIT
DIHEDRAL  1 C1  1 C2  1 C3  1 C4      0.0
DIHEDRAL  1 C2  1 C3  1 C4  1 C5    -60.0
DIHEDRAL  1 C3  1 C4  1 C5  1 C6    +60.0
DIHEDRAL  1 C4  1 C5  1 C6  1 C1      0.0
DIHEDRAL  1 C5  1 C6  1 C1  1 C2    -60.0
DIHEDRAL  1 C6  1 C1  1 C2  1 C3    +60.0
END
! Initialize the coordinate values and build new
! coordinates
COOR INIT
IC SEED  1 C1  1 C2   1 C3
IC BUILD
WRITE COORDINATES CARD UNIT 45
* Boat conformation coordinates for cyclohexane
*
NBONDS CUTNB 99.0 CDIE VSWITCH SWITCH
ENERGY
COOR DIST ENERGY CLOSE 14EXCLUSIONS NONBONDS
STOP


References

 Slater, J. C. and Kirkwood, J. G., Phys. Rev., 37 682 (1931).

Ryckaert, J. P. and Bellemans, A., Chem. Phys. Lett., 30 123 (1975).

Ryckaert, J. P. and Bellemans, A., Disc. Farad. Soc., 66 95 (1978).



Last updated September 18, 1998 at 04:02PM PDT.
Copyright © 1997, Molecular Simulations, Inc. All rights reserved.