Sb intro1

Starting up CHARMM
While you can run CHARMM interactively, one usually tells the program what to do by means of a script. Under Unix (at least for non-parallel versions of the program), this means that in order to execute a (short) CHARMM calculation, one runs from the command line (Unix Shell prompt)

charmm_executable < charmm_input_script.inp

exploiting input redirection available under all Unix shells. (TODO: what about other OS (Windows, OSX ?), explain non-parallel remark. Possibly should state somewhere that tutorial assumes Unix/Linux??). Since as we shall see shortly CHARMM output tends to be verbose, one normally also redirects the output to a file, thus ending up with

charmm_executable < charmm_input_script.inp > charmm_output_file.out

Of course, instead of charmm_executable use the filename starting up CHARMM at your location and replace charmm_input_script.inp and charmm_output_file.out by the names of the actual script/output file which you want to run (to which you want to save your output).

Titles
First, let's do something really silly and start up charmm reading from an empty file; which can be easily accomplished by executing

charmm_executable < /dev/null

CHARMM prints a header telling you copyright info, version and some more stuff, followed by a warning

RDTITL> No title read. ***** LEVEL 1 WARNING FROM  ***** ***** Title expected. ******************************************     BOMLEV (  0) IS NOT REACHED. WRNLEV IS 5

The job finishes by printing some status info. The interesting part is the warning from which we learn that CHARMM expected a "title". Indeed, each CHARMM script should start with a title, and if the main script tells CHARMM to read from another file, the program also expects to find a title at the beginning of that file.

A title should not be confused with comments. E.g., it can only occur at the beginning of a file (we'll explain the apparent exceptions when we encounter them). Title lines start with a star *; to indicate the end of the title give a line containing only a star. (A title can consist of up to 32 consecutive lines) Thus,

* This would be a short title *

If you start CHARMM with a short file containing the above snippet (=title), you get the title echoed in uppercase letters

RDTITL> * THIS WOULD BE A SHORT TITLE RDTITL> *

instead of the warning when using the empty file.

Comments
Having blabbered so much about titles, what are comments: A comment in a CHARMM script is everything following an exclamation mark ! I.e.,

! this is a comment on a line by itself

and this would be a line containing a CHARMM command, followed by a comment

ENERgy ! as you might expect, this command calculates an energy

Ending a CHARMM script
So far, CHARMM finished when it reached the end of the script file (as the line

NORMAL TERMINATION BY END OF FILE

in the output informs you. It's OK to end a CHARMM script in this manner, but the preferred way of stopping CHARMM is by issuing the

stop

command. We, thus, can create a first, completely useless CHARMM script, looking like

* A CHARMM script doing nothing * ! we really should be doing something, but in the meantime all we know is stop

In addition to the title, the comment is echoed as well. Note that CHARMM now prints upon finishing

NORMAL TERMINATION BY NORMAL STOP

A first full example
One of the major hurdles in getting started with CHARMM is the initial setting up of one's simulation system (usually starting from a PDB file of the system (protein) of interest). This is where sites such as http://www.charmming.org prove so helpful, and we shall take a detailed look at what is involved in the initial setup shortly. For the time being, however, let's assume that we have accomplished this first step and that we want to carry out a short calculation with our system, e.g., a minimization. The model system is crambin (PDB entry: 1CRN), and you can download the necessary files here (Tim: how do I upload to the main server?. Create a directory and unpack the tar.gz archive there. You end up with five files The first four are data files needed to set up the system, the fifth file is the CHARMM script itself. Before we can dig in, we need to say a few words about CHARMM data structures
 * 1) top_all27_prot_na.rtf
 * 2) par_all27_prot_na.prm
 * 3) 1crn.psf
 * 4) 1crn_init.crd
 * 5) 1crn_1.inp

Some comments about CHARMM data structures
Incidentally, the 1-4 you just downloaded correspond to the four most important data structures required in almost any CHARMM run:


 * RTF (residue topology file)
 * The so-called residue topology file contains the information how the building blocks of larger molecules (e.g., amino acids in the case of proteins) look like. In the context of CHARMM the term residue is used not only for amino acids, but also, e.g., for nucleic acids or sugars, as well as small molecules, e.g., water. For each residue, the RTF provides information about the atoms it contains, the atom type, and the atomic charges. Further, it specifies the connectivity of the atoms (which atom is (chemically) bonded to what other atoms). From these data, CHARMM can usually derive information about angles and dihedral angles itself, although one can/could also specify them by hand. CHARMM has RTFs for nucleic acids, lipids, proteins and carbohydrates. E.g., the file top_all27_prot_na.rtf we use here is part of the CHARMM distribution and contains residue information for amino acids and nucleic acids parametrized with the CHARMM27 force field (TODO:  should probably give some instructive refs). Without going into details, a short look at the file should give you an idea of the type of information it contains.


 * Parameters
 * The force field parameters in (more or less) human readable form; e.g., force constants of the bond stretching terms, Lennard-Jones radii and  well-depths for all atom types etc. (The notable exception are partial atomic charges, which are listed in the RTF --- can you think of the reason why?!) The file  par_all27_prot_na.prm  used in our example is the companion of top_all27_prot_na.rtf and contains parameters for all compounds described by this RTF file (amino acids and nucleic acids). Again, a short glimpse at the file should give you an idea of the info it contains.


 * PSF (protein structure file)
 * The RTF and parameter files (and once processed by CHARMM the corresponding data structures) contain information that applies to, e.g., all proteins by specifying how amino acids look like. Obviously, you want to simulate a specific protein, characterized by its unique primary sequence. The PSF datastructure contains exactly this system specific information, such as the exact list of all bond stretching terms in your protein, all angle bending terms etc. and much more. Before you can do anything useful with CHARMM, this information must have been created internally (we'll see later how this is done). Once created, the PSF data structure can be saved to a file for later reuse; this is what we shall be doing in our crambin example.
 * Note that the PSF does not contain the force field parameters. As just mentioned, all parameters are read from a parameter file and stored in (an) appropriate data structure(s). The parameters that are needed for the specific system of interest (e.g. crambin) are internally stored in something called the CODES data structure. Normally, you don't have to be concerned with this; however, if you get a CODES error, this normally indicates that your parameter file was incomplete or inconsistent.


 * Coordinates
 * Finally, you need coordinates for each atom in the system. The most frequent source of coordinates in real work (with proteins) are PDB (protein database) files, but they need some "massaging" before they can be used by CHARMM. Therefore, we use preprepared coordinates for our example.

Normally, one RTF and parameter file suffices for a complete class of systems (e.g., proteins), and these file are provided together with CHARMM; they constitute what is referred to as "the CHARMM force field". PSF and coordinates, on the other hand, are specific to the system being studied. At some point, they have to be created by the user of the program (again, we'll see shortly how this is done!)

Back to the crambin minimization example
Provided one has a PSF and coordinates for the system of interest, many typical applications of CHARMM follow a rather simple pattern:
 * (A) Read RTF and parameter file (i.e., load force field information into CHARMM)
 * (B) Read PSF and coordinate file (i.e., load information about specific system)
 * (C) If necessary, do some more preparatory stuff
 * (D) Do some work
 * (E) Possibly do "on the fly" analysis of what was done in (D)
 * (F) Frequently, (D) changed the coordinates; so save the coordinates and whatever else may have changed.

Reading files
Please do take a look at the downloaded example CHARMM script 1crn_1.inp now. As you may expect, it contains a title; followed by a lengthy comment paraphrasing what was just said. And, then it starts off by indeed reading the RTF file, which is accomplished by the following three commands

open unit 1 form read name top_all27_prot_na.rtf read rtf card unit 1 clos unit 1

An analogous triplet of commands is carried out for reading the parameter file, the PSF and the coordinate file, so it is worthwhile to spend a few words on what's going on here: I/O (Input/Output) in CHARMM follows the model of most programming languages; since CHARMM is programmed in FORTRAN, there are significant similarities to the FORTRAN I/O model. First, a file has to be opened ("connected" to a "unit"); the remaining arguments to the open command specify the mode of operation. Note: units 5 and 6 are reserved; on most systems you can have unit numbers ranging from 1 to 99. Obviously, you have to specify the file name (name ). The keyword read indicates that the file is opened for read access (meaning that it is an error if the file does not exist!). Further, we specify that the file is "formatted" (keyword form (alternatively one could use card)), i.e., it is an ASCII text file as opposed to a binary file (in which case the keyword would be file</tt>; file</tt> is actually the default, so form</tt> must be specified!).

The next command ( read rtf card 1</tt>) does the actual reading. The first keyword after read</tt> specifies what type of file is going to be read; not surprisingly, the keyword is rtf</tt> for RTFs, para</tt> for parameter files etc. (see 1crn_1.inp). The second keyword after the read</tt> command has to do with the fact that many types of files can be read by CHARMM in more than one format (we'll see examples later!). In all read</tt> statements in 1crn_1.inp we specify the respective native CHARMM format by the card</tt> keyword. Finally, we have to tell CHARMM from which unit</tt> we want to read; typically (as is the case here) this is the unit that was opened just before the read</tt> command ( read ... unit 1</tt>).

Finally, files opened for reading, should be explicitly closed when one is done. (This is not necessary, but good practice).

A brief comment on the minimization
The coordinates we read are basically the coordinates from the 1CRN pdb file of crambin, with missing atoms (e.g., hydrogens added). The simple command

energy

computes the energy of the system, using a multitude of default settings. This is one reason why this example is somewhat artificial! Many default settings in CHARMM are, unfortunately, not optimal, and normally one would carefully tune options here (more on this, in particular the necessary background, later). Further, we have crambin without surrounding solvent, so the subsequent minimization of the protein in "vacuum" is not very realistic. For the moment, we can ignore these caveats and take a look at the output of the energy</tt> commmand instead.

CHARMM>   energy NONBOND OPTION FLAGS: ELEC    VDW      ATOMs    CDIElec  SHIFt    VATOm    VSWItch BYGRoup NOEXtnd  NOEWald CUTNB = 14.000 CTEXNB =999.000 CTONNB = 10.000 CTOFNB = 12.000 WMIN  =  1.500 WRNMXD =  0.500 E14FAC =  1.000 EPS    =  1.000 NBXMOD =     5 There are       0 atom  pairs and        0 atom  exclusions. There are       0 group pairs and        0 group exclusions. <MAKINB> with mode  5 found   1835 exclusions and   1713 interactions(1-4) <MAKGRP> found   578 group exclusions. Generating nonbond list with Exclusion mode = 5 == PRIMARY == SPACE FOR  184479 ATOM PAIRS AND        0 GROUP PAIRS General atom nonbond list generation found: 119614 ATOM PAIRS WERE FOUND FOR ATOM LIST 4229 GROUP PAIRS REQUIRED ATOM SEARCHES ENER ENR: Eval#     ENERgy      Delta-E         GRMS ENER INTERN:         BONDs       ANGLes       UREY-b    DIHEdrals     IMPRopers ENER CROSS:          CMAPs ENER EXTERN:       VDWaals         ELEC       HBONds          ASP          USER --      -    -    -    -    - ENER>        0   -296.51905      0.00000     15.23577 ENER INTERN>     117.51821    131.42640      9.84634    215.59668       2.29425 ENER CROSS>      -64.40694 ENER EXTERN>     -91.86576   -616.92824      0.00000      0.00000       0.00000 --      -    -    -    -    -

The lines immediately following the energy command report (some of) the settings used to compute the energy (e.g., the use of a 12 A cut-off radius (CTOFNB), of a shifting function for the electrostatic (SHIFt) and a switching function (VSWItch) for Lennard Jones interactions). The following lines contain some information about the generation about the generation of the nonbonded lists. (TODO: this should point towards my material (possibly updated) about energy calculation) The actual output of the energy calculation can be found in the lines starting with ENER</tt>. The header lines should be mostly self-explanatory, in particular the total energy can be found as entry ENERgy</tt>. UREY-b</tt> denotes the Urey-Bradley like bonded terms in the 1-3 distance, <tt>CMAP</tt> is the new dihedral cross term (TODO: some ref?), and <tt>VDWaals</tt> in CHARMM jargon refers to the Lennard-Jones interactions.

The next two commands do the real work (minimization). Obviously, <tt>mini</tt> initiates the minimization. The keyword following <tt>mini</tt> indicates the minimizer to be used; here we start with steepest descent (SD), followed by adopted basis set Newton Raphson (ABNR) [give reference]. Details of these minimizers will be discussed later. The option <tt>nstep 100</tt> tells the program to carry out 100 minimization steps. The behavior of the minimizers could be finetuned by method specific additional options, but this is beyond the scope of this introduction. The output from the minimization resembles that of the energy calculation

MINI MIN: Cycle     ENERgy      Delta-E         GRMS    Step-size MINI INTERN:         BONDs       ANGLes       UREY-b    DIHEdrals     IMPRopers MINI CROSS:          CMAPs MINI EXTERN:       VDWaals         ELEC       HBONds          ASP          USER --      -    -    -    -    - [... snip ...] MINI>      10   -492.83381    196.31476      1.39253      0.00156 MINI INTERN>      31.06297     99.50559      7.06564    212.92448       4.12544 MINI CROSS>      -68.59167 MINI EXTERN>    -129.29653   -649.62974      0.00000      0.00000       0.00000

More about minimization later on. At present, we just want to point your attention to the fields <tt>Delta-E</tt> and <tt>GRMS</tt>. The first reports the change in energy to the previous printout, i.e., if you add the 196.31 kcal/mole of the <tt>MINI> 10</tt> line to the energy at this step (-492.83 kcal/mole), you obtain the initial energy of the system (-296.52 kcal/mole, see the <tt>MINI> 0</tt> line, or the output of the energy command). The gradient root mean square (GRMS) field reports the accompanying root mean square deviation in the forces. If you skim through the output, you can see how the energy goes down; if you set <tt>nstep</tt> of the ABNR minimization to something very high (e.g., 10000), you would most likely see that both <tt>Delta-E</tt> and <tt>GRMS</tt> eventually become zero; possibly, the minimization would even stop if some internal convergence criteria are fulfilled.

It should be pointed out that the <tt>ener</tt> command is redundant; we could have just carried out the minimization. However, you'll often find an energy command after everything is set up, and we maintain this tradition. If, for some reason something went wrong up to this point, the energy calculation will most likely fail and the CHARMM run will abort at this (early) point. Also, one would usually add to the <tt>ener</tt> command all the various options pertaining to how energy and forces are calculated exactly; these are remembered throughout the remainder of the CHARMM script and commands like <tt>mini</tt> or <tt>dyna</tt> (molecular dynamics) need to have only options controlling details of minimization or MD, rather than energy options as well. We shall see this in later examples, where attention is paid to these details.

Writing files
At the end of the minimization, our system does not only have a much lower energy, but its coordinates have changed as well. Most likely, one wants to save these coordinates for later use, and so we have to learn how to write files. Writing files in CHARMM is very similar to reading files. First, the file is opened (with option <tt>writ</tt>

open unit 11 form writ name 1crn_somemini.crd

instead of <tt>read</tt>; in this case, the file need not exist), followed by the

writ coor card unit 11

statement, which writes to the just opened unit 11. As in the <tt>read</tt> command we have to tell CHARMM what type of file (<tt>coor</tt> for coordinates) we plan to write, and in what format (<tt>card</tt> for native CHARMM format in the above snippet). If you look at <tt>1crn_1.inp</tt>, you'll see that we actually write the coordinates twice, once in CHARMM format (keyword <tt>card</tt>, but also in PDB format (keyword <tt>pdb</tt>). The native CHARMM format is intended for reuse by CHARMM, whereas the PDB file can be used for data exchange with other programs, in particular (but not limited to) graphics programs. The novel part (compared to reading) is that the <tt>writ</tt> command is followed by a title

writ coor pdb unit 11 * CRAMBIN, from pdb entry 1crn, * with CHARMM27 parameter set * after 100 steps SD + 100 steps ABNR minimization *

This provides the possibility to pass some information into the respective output file (take a look at <tt>1crn_somemini.crd</tt>, the title is there, augmented by some "accounting information"). The title can be ommitted, but providing an informative(!) title is good practice. After most type of writes, files are closed automatically; so there is usually no need for an explicit <tt>clos</tt> command.

A closing comment
Congratulations, you have just worked your way through a complete (and almost meaningful) CHARMM run/script. Lets take one more look at <tt>1crn_1.inp</tt> from a different point of view. Strip the file of all comments and titles, and then count the number of commands. Given that the <tt>ener</tt> command is actually superfluous, the real work is initiated by just two lines (the two <tt>mini</tt> commands). All other commands (I count 18!) do routine stuff like reading in the force field (RTF and parameters) or write out the final coordinates. This ratio is not too untypical; while the "hard work" commands of CHARMM do very much controlled by a single command line, they need to be surrounded by multiple "lesser" commands taking care of the small details. There are two other characteristics of (most) CHARMM scripts that you'll notice soon. (i) Many of them are very similar. Had we decided to run a short MD simulation as a first example, our script would have looked very similar. (ii) If you want to apply our minimization script to a different protein (provided you have a PSF and coordinates), all you would have to do is to change the filenames. In the following, we want to show you capabilities of CHARMM scripting that help you make your scripts more general (e.g., to write one script that could be used to minimize any system for which you have a PSF and coordinate file), as well as modular (i.e., to split the script into recurring parts, such as reading of RTF and parameters, as opposed to the work part, which changes depending whether you want to do MD or minimization or normal coordinate analysis or ...

Command line parameters, @variables
So far, the CHARMM scripting language seems to be a concatenation of individual commands. It does contain, however, most (if not all) elements of a (imperative) programming language (even though it is not an extremely comfortable one). One key element of a programming language is the ability to set and manipulate variables. In the context of running CHARMM, it is extremely useful to pass certain variables into the program from the outside. Run the following miniscript (we name it title2.inp)

* Example of passing a variable from the command line * The variable myvalue was initialized to @myvalue
 * stop

by executing charmm_executable myvalue=500 < title2.inp

and study the output (you could also state <tt>myvalue:500</tt>): Before the title is echoed, you see that an argument was passed

Processing passed argument "myvalue:500" Parameter: MYVALUE <- "500"

and in the echoing of the title, @myvalue is replaced by the value assigned to the variable myvalue, i.e., 500 in our case. The little example highlights something to keep in mind. The first thing that is done when a line in a script (containing a title or a command) is parsed, is to scan for @variables, which are then replaced by their value (no variable replacement is done in comments!).

Thus, we can make an easy generalization to our minimization script <tt>1crn_1.inp</tt>. Rather than hardcoding the number of minimization steps, we could pass the number from the command line. Replace, e.g., the line mini sd nstep 100

by mini sd nstep @nstep

and execute the script as

charmm_executable nstep:10 < 1crn_1.inp > my_output_name

and you will carry out 10 steps of SD minimization. If, however, you specified <tt>nstep:1000</tt>, you would carry out 1000 SD minimization steps. (Obviously, you could do the same for the ABNR minimization command).

If you follow some naming conventions, you can generalize you script even further. Note that all files pertaining to our model protein (crambin) that are read or written by <tt>1crn_1.inp</tt> start with "1crn" (the PDB code of the crambin structure used). If you replace "1crn" by "@system" and pass system:1crn on the command line, you have exactly the functionality of the original script. However, suppose that you have a psf and coordinate file for some other protein, then you can use the very same script for this protein as well. Let's assume that your modified <tt>1crn_1.inp</tt> script is named <tt>simplemini_1.inp</tt> to distinguish it from the original downloaded version, then executing

charmm_executable nstep:100 system:1crn < simplemini_1.inp > some_output_name

does exactly what the original <tt>1crn_1.inp</tt> script did, where everything was hardcoded. If you specify

charmm_executable nstep:500 system:1ccn < simplemini_1.inp > some_other_output_name

then you carry out 500 steps of (SD) minimization for a protein to which you have given the label 1ccn (1CCN would be the PDB code for an NMR structure of crambin).

Our new script is not free from problems. First, it contains hardcoded titles/comments referring to crambin. Second, if you forget to pass one or more variables on the command line, you CHARMM job will "bomb" (we'll show how to do some error handling for this a bit later). Finally, if you use the @system trick, you have the variable as part of a longer string. E.g, you want that @system_init.crd gets replaced by 1crn_init.crd. In this case this works fine, but suppose you wanted to substitute @system in @systemproperty. This will fail, since CHARMM searches for a variable "systemproperty". You can avoid this by writing @{system}property, which (given that @system expands to 1crn) expands to 1crnproperty. Whenever in doubt, use @{varname} rather than just @varname!

Needless to say, one can also set variables in a script itself, manipulate them, and calculate with them. Before we show some examples, let us show you two ways how one can see/check the result of the calculations/manipulations.

The <tt>echo</tt> command: does more or less the same as its Unix counterpart; everything following <tt>echo</tt> is written to standard output (one can also force it to write to a different file unit, see miscom.doc). Variables are expanded, so the output will contain the value of the variables rather than the @variable string. Suppose <tt>@test</tt> contains the value 625. Then

echo value of test is @test

will result in the output line

VALUE OF TEST IS 625

in the CHARMM output

The <tt>write title</tt> command: We have seen that @variables get expanded in titles. Just as one can write coordinates (<tt>writ coor</tt>) or PSFs etc., one can also write an arbitrary title. This is the most flexible way to get "generic" output through CHARMM (i.e, to write an almost arbitrary piece of text). The trick here consists in creating a title on the fly and writing it to standard output (or, alternatively, a separate file). In Unix, standard output is associated with (file) unit 6, and so we are going to write our title to unit 6. Look at the following code snippet writ title unit 6 * some title on the fly * which may contain @variable values which will be substituded *

The <tt>writ title</tt> command (as any <tt>writ</tt> command) expects that a title follows, on which normal command line substitution is carried out. This title will be written to the specified unit, in the above example standard output, i.e., it gets merged into the normal CHARMM output. Consider the following example (which we name title3.inp) * Example of passing a variable from the command line * The variable myvalue was initialized to @myvalue * ! we now set a second variable ourselves set 2value 100 ! and check that we can use it write title unit 6 * checking values of variables * myvalue : @myvalue * 2value : @2value stop

Run it as charmm_executable myvalue:500 < title3.inp and observe the output. You see that the title lines following the <tt>write title</tt> appear twice; once, as in the beginning of the script following the <tt>RDTITL> *</tt> string, and a second time all by themselves, without leading <tt>*</tt>. (BTW: try replacing <tt>unit 6</tt> by, e.g., <tt>unit 51</tt> and precede the command by  open writ form name some_file_name unit 51

If you rerun this modified title3.inp script, you will notice that the lines CHECKING VALUES OF VARIABLES MYVALUE : 500 2VALUE : 100

are now missing from the output, but you will find them in the new file named <tt>some_file_name</tt> instead.

Manipulating variables
Based on the framework of title3.inp we can now illustrate how variables can be manipulated. Please try the following examples and make sure that you understand what's going on:

varmanip1.inp:

* First example of manipulating variables * ! set a variable set value 100 ! and check its value write title unit 6 *checking values of variables *value : @value ! now increment it by one incr value by 1 ! and check its value write title unit 6 *checking values of variables *value : @value stop

Of particular interest is the <tt>calc</tt> facility of CHARMM, which permits fairly complex calculations. Please take a look at calc.doc (TODO: get these refs right); the following just skims the surface

varmanip2.inp:

* Example of calculating with calc * ! set a variable set value 100 ! and check its value write title unit 6 *checking values of variables *value : @value ! now increment it by one calc value = 2 * @value calc sqvalue = sqrt ( @value ) ! and check the results write title unit 6 *checking values of variables *value  : @value *sqvalue : @sqvalue stop

The example should be self-explanatory. Some care is required when and where to use @ (typically on the right hand side of an expression), and to leave spaces between and mathematical operators and variable names, e.g.

calc sqvalue = sqrt ( @value )

in the example (see the documentation for further details!)

A note on ?variables
For completeness sake, a preliminary comment on a second type of variable. Many CHARMM commands, aside from producing more or less direct output, place some key values in "internal" variables. E.g., the energy of a system calculated with the most recent <tt>ENER</tt> command is put into a variable named <tt>ener</tt>. To avoid name clashes with variables set by users, these value of these variables can be accessed by preceding the variable name by a question mark ?. I.e., if CHARMM encounters (in a title or in a command) the string <tt>?ener</tt>, it will attept to replace it with the energy value from the most recent <tt>ENERgy</tt> command. This is quite handy, and we'll see many examples later on, once we reached the stage where we can use CHARMM to work with biomolecular systems.

A first revised version of our full example
As an example of the things we have learned, let's study a revised version of our first full example (<tt>1crn_1.inp</tt>). We call the file <tt>simplemini.inp</tt>, and you can download it here. It incorporates the ideas outlined earlier and uses ?variables to report the energy change upon minimization in the coordinate files saved.

Some shortcomings remain: (1) If you forget to pass the necessary command line variables, your job dies somewhere in the middle of the script. It would be nice if some error checking were done at the beginning of the script and / or default values were provided were possible. We'll show shortly how to address this. (2) Regardless of the number of minimization steps requested, coordinates are saved to the same files, potentially overwriting earlier outputs. It really depends on you (what you want to accomplish) whether you consider this a problem or not, so there is no generic solution to offer.

Branching
An essential element of any programming language is branching. In CHARMM scripting you can do things along the lines

IF comparison THEN ! do something ELSE ! do something else ENDIF

There is no explicit <tt>ELSE IF</tt>, but <tt>IF</tt>s can be nested.

Example 1: In our work, we occasionally run MD simulations of a system in the gas phase and with an implicit solvent model. Since in this case the same PSF and coordinate file can be used, w prefer to use one and the same script for both environments (gas phase and "solution"). Depending on the environment, the options for energy calculations have to be set differently; also, for the calculations with implicit solvent some additional commands are required. We pass from the command line a variable <tt>env</tt>, which should be set to "gasp" or "solv" (with obvious meaning). Our script(s) then contain(s) code like the following

if @env eq gasp then energy ! energy options for gas phase, e.g., infinite cut-off etc. endif if @env eq solv then energy ! energy options as required by implicit solvent model ! setup implicit solvent model energy ! final energy command to check that everything is OK endif

See miscom.doc for the full syntax; more examples will follow shortly.

Labels and goto
to be done

Loops
to be done

Stream files
to be done

Another closing comment
Hopefully we could convince you that you can do some nifty things with CHARMM scripts. We use many of these "tricks" to keep our own scripts managable and (reasonably) generic. On the other hand, there are limits what CHARMM scripting can do, and we discourage you "to stretch the limits". Fancy CHARMM scripts can rapidly become incomprehensible (no matter how many comments you write!); hence, if things go wrong (and we mean the type of wrong which you notice after days of simulations), it's often very tedious to spot the error. Despite all its features, the programming capabilities of CHARMM scripts are vastly inferior to those of Perl, Python or even the bash shell. Thus, if you find yourself constructing complex loops within CHARMM scripts or issueing tons of <tt>system</tt> commands, you should ask yourself whether it would not be better to construct the loops in Perl (Python, bash, ruby ...) and call simpler CHARMM jobs from this external script (passing job specific info via the command line). Yet another approach is taken, e.g., by. The site goes one step further and writes CHARMM scripts specific for your task and system on the fly.

As you journey into the depths of CHARMM progresses, keep these options in mind and do not forgot to weigh elegance and generality versus simplicity.