CHARMM: A user

1
CHARMM A users perspective
V. M. Dadarlat, March 2007

• CHARMM - a research program developed at Harvard
  University for energy minimization and dynamics
  simulation of proteins, nucleic acids and lipids
  in vacuum, solution or crystal environments
  (Harvard CHARMM Web Page http//yuri.harvard.edu/)
  .
• - Models the dynamics and mechanics of
  macromolecular systems using empirical and mixed
  empirical/quantum mechanical force fields.
• - Designed to investigate the structure and
  dynamics of large molecules. It performs free
  energy calculations of mutations and drug binding
  as well as conformational folding of peptides.
• - Uses classical mechanical methods to
  investigate potential energy surfaces derived
  from experimental and "ab initio" quantum
  chemical calculations.
• - Mixed quantum mechanical/classical systems can
  be defined to investigate chemical processes such
  as enzyme catalysis.

2

• CHARMM - a command line program, i.e. a command
  is read from the input stream (typed, or from a
  file) and acted upon.
• References
• CHARMM A Program for Macromolecular Energy,
  Minimization, and Dynamics Calculations, J. Comp.
  Chem. 4, 187-217 (1983), by B. R. Brooks, R. E.
  Bruccoleri, B. D. Olafson, D. J. States, S.
  Swaminathan, and M. Karplus.
• CHARMM The Energy Function and Its
  Parameterization with an Overview of the Program,
  in The Encyclopedia of Computational Chemistry,
  1, 271-277, P. v. R. Schleyer et al., editors
  (John Wiley Sons Chichester, 1998), by A. D.
  MacKerell, Jr., B. Brooks, C. L. Brooks, III, L.
  Nilsson, B. Roux, Y. Won, and M. Karplus.
• http//www.psc.edu/general/software/package
  s/charmm/charmm.html
• http//www.psc.edu/general/software/package
  s/charmm/tutorial/index.html
• http//www.ch.embnet.org/MD_tutorial/

3
Potential energy functions in CHARMM - CHARMM22
- all atom potential function for -
proteins (MacKerell et al 1998), -
nucleic acids (MacKerell et al. 1995), -
lipids (Schlenkrich et al. 1996) and
carbohydrates (Ha et al. 1988). - CHARMM27 -
CHARMM19
4
What CHARMM can do for you

• Main features of the CHARMM program
• Molecular Mechanics (Energy Minimization)
• Molecular Dynamics
• Classical Simulations
• Crystal Simulations
• Specialized techniques (Umbrella Sampling,
  Constant Pressure)
• Miscellaneous Mean Field Potentials (MMFP)
• Langevin/Implicit Euler Dynamics
• Monte Carlo simulation package
• Free Energy Perturbation Calculations
• Block program
• Perturb program
• Thermodynamic Simulation Method (TSM)

5
- Normal Mode and Molecular vibrational analysis
facility - Reaction Path Determination - Genetic
Algorithm - Advanced implicit solvent models -
Analytical Continuum Electrostatics (ACE) model
- Effective Energy Function 1 (EEF1) -
Generalized Born Solvation Energy (GBorn)
Interfaces - Multi-body dynamics
(MBOND) - Merck Molecular Force Field
(MMFF) - Quantum and Molecular
Mechanical Force Field (QM/MM) - Analysis
facility - Time series and correlation function
calculations - NMR analysis facility -
Poisson-Boltzmann Equation Solver - Reference
Interaction Site Model (RISM)
6
Controlling a CHARMM Run Commands

• IF command-parameter test-spec
  comparison-string command-spec
• GOTO label-string
• LABEL label-string STREAM UNIT
  integer file-specification
• RETURN
• SET command-parameter string
• INCRement command-parameter BY real
• DECRement command-parameter BY real
• IF tests to conditionally execute a single
  command
• GOTO and LABEL transfers within a file
• STREAM and RETURN transfers to different command
  files.
• The logical unit in OPEN, CLOSE, and REWIND
  commands are useful in working with streams.

7
CHARMM units for quantities The AKMA system
(Angstroms, Kilocalories/Mole, Atomic mass
units). - Distances in Angstroms
- Energies in kcal/mole - mass in
atomic mass units - charge is in units
of electron charge. - angles in
degrees (analysis and constraint sections)
the force constants for angles, dihedrals,
and dihedral constraints are specified in
kcal/mole/radian/radian.
8
Running MD with CHARMM

• A. BUILD STRUCTURE
• B. EQUILIBRATION
• C. PRODUCTION RUN
• EQUILIBRIUM DYNAMICS

9
Input files

• Structure files
• Topology file (top.inp) Parameter file
  (parainp)
• Data structures, information about the
  molecule
• composition, chemical
  connectivity,
• atomic properties, internal
  coordinates and
• parameters for the energy
  function.
• For a specific molecule, the data is extracted
• from these files and stored in the Protein
  Structure File (PSF).
• The Coordinates (COOR)
• Coordinate file from .pdb (from the Protein
  Data Bank) to a
• CHARMM specific format .crd
• The coordinates are Cartesian for all the atoms
  in the PSF
• 2 coordinate files main and comparison sets
• CHARMM script file kidkixdynamics.inp

10
Output Files

• Standard output .out (time, energy,
  temperature, statistics, etc.)
• Trajectory file (binary) .dcd (or .trj) a
  time series, collection of snapshots of the
  positions of each atom in the system as a
  function of time
• A restart file, .rst contains positions and
  velocities of all atoms at the end of the current
  dynamics run
• A time series of velocities can also be saved

11
A. BUILD STRUCTURE
1fsc.pdb structure representation showing
disulfide bonds.

1. Download 1fsc.pdb from PDB
2. Import pdb file in QUANTA gui
3. Pdb format to crd format
4. Can build hydrogens
5. Find and show disulfide bonds, patches
6. Structure comparisons

12

• Initial decisions and analysis read the pdb file
  carefully
• - Titratable groups what to do with Histidines,
  Glutamic and Aspartic acids?
• - Xray waters
• - Missing residues
• - bound ions and coordination complexes
• Disulfide bonds
• Solution conditions if NMR
• 20 NMR lowest energy structures which one?
• Minimization and solvation
• Minimize structure with CHARMM force field,
  center the molecule, get structural info Xmin,
  Xmax Ymin, Ymax Zmin, Zmax
• - Patches for protein and ligands and any
  additional rtf and parameter information all into
  the solvation script.
• - Decide on periodic boundary conditions image
  or crystal, cubic or octahedral (box)
• Build adequate water box, with sizes Xmin-12 A,
  Xmax12 A
• Equilibrate water box
• Overlay protein in water box by removing all
  water molecules that overlap on the protein

13
Histidine HSD resid 6
Histidine HSE resid 29
14
Figuring out Titratable groups for
Histidine Charmm lt histidines.inp gt
histidines.out CHARMM command COOR MINDist will
give the closest atom to each ring nitrogen in
each histidine important output as follows. 
Residue 6 has an electron donor near ND1
indicating it should be protonated
NE2 has an electron acceptor near it 
indicating it should not be protonated.
We should decide to make residue 6 a HSD.
Residue 29 a HSE since there is nothing near
ND1 but an acceptor near
NE1. MINIMUM DISTANCE FOR ALL
SELECTED ATOMS 83 1FSC HSD 6 ND1 - 106
1FSC THR 7 O 2.8497
MINIMUM DISTANCE FOR ALL SELECTED ATOMS
88 1FSC HSD 6 NE2 - 583 1FSC ARG 37 HH12
1.9869 MINIMUM
DISTANCE FOR ALL SELECTED ATOMS 446 1FSC HSD
29 ND1 - 420 1FSC ARG 28 HB1 3.6827
MINIMUM DISTANCE FOR ALL
SELECTED ATOMS 451 1FSC HSD 29 NE2 - 690
1FSC ASN 47 OD1 2.6855
15

• Equilibrate
• minimize solvent with solute fixed to relieve bad
  solvent vdW contacts
• Run several ps CPT MD at 300K with weak
  constraints on the solute
• Run several hundred ps  CPT MD to equilibrate
  all.
• Analyze energy
• Analyze pressure and temperature
• Production Run
• Depending on the system and investigated
  properties you have to decide how long your
  simulation should run. For some systems, e. g.
  folded proteins, few ns would be enough to assess
  stability, volume fluctuations, etc. For PMF
  calculations from MD simulations, hundreds of ns
  are needed for sampling the conformational space
• Analyze
• RMSD with respect to the x-Ray (NMR) structure
• Atomic Fluctuations
• rmsd time series
• Energy, temperature, density time series
• P2 correlation functions and order parameters

16
EXAMPLE MD SIMULATIONS of 2 ALANINE DIPEPTIDE
MOLECULES IN WATER
17

• CHARMM Input file - script generate a potential
  energy surface for alanine dipeptide(example)
• This is a sample command file for
  CHARMM which calls a stream file
• to build a structure and then maps
  out an adiabatic potential
• surface defined by a pair of
  dihedrals
• OPEN UNIT 10 READ FORM NAME
  makestruc.inp
• STREAM UNIT 10
• SET 1 -180.
• SET 2 -180.
• LABEL LOOP
• CONS CLDH
• CONS DIHE first-dihedral-angle-spec
  FORCE 100.0 MIN _at_1
• CONS DIHE second-dihedral-angle-spec
  FORCE 100.0 MIN _at_2
• MINI minimization-spec
• INCR 1 BY 30.0
• IF 1 LT 170. GOTO LOOP
• SET 1 -180.
• INCR 2 BY 30.0

18
file adp2.dynam.inp EQUILIBRIUM DYNAMICS
this program starts 5000 CPT dynamics steps with
a timestep of 0.002 bomlev -2 set 1
/bio/terry3/charmm/c30a2x1/toppar open unit 1
read form name _at_1/top_all22_model.inp read rtf
card unit 1 close unit 1 open unit 2 read form
name _at_1/par_all22_prot.inp read param card unit
2 close unit 2 read sequence card Sequence for
generating adp 1 ALAD generate ADP1 read
sequence card Sequence for generating adp
1 ALAD generate ADP2 read sequence tip3 990
19
generate BULK noangle nodihe open unit 1 form
read name 2adp990wateq.crd read coor card unit 1
close unit 1 ! Define the box for PBC set
theta 90.00 set alongx 30.0 set alongy
30.0 set alongz 30.0 CRYSTAL DEFINE CUBIC
_at_alongx _at_alongy _at_alongz _at_theta _at_theta
_at_theta CRYSTAL BUILD CUTOFF 15.0 NOPERATIONS
0 !coor dist sele type CA .and. segid ADP1 end
sele type CA .and. segid ADP2 end !coor copy
comp Image byres select resname TIP3 end Image
byseg select segid ADP1 end Image byseg select
segid ADP2 end open unit 20 read card name
2adp5free.rst.in open unit 30 write file name
2adp5free.trj open unit 40 write card name
2adp5free.rst
20
update inbfrq -1 imgfrq 20 - cutim 14.0
ctonnb 10.0 wmin 1.0 - CTOFNB 12.0 CUTNB
14.0 - ! atom vatom vswitch - ewald
pmewald order 6 kappa 0.340 fftx 32 ffty 32 fftz
32 !mini abnr nstep 200 LATTICE shake bonh param
tol 1.0e-09 !!!!!!Dynamics Specifications DYNAMI
CS nstep 20000 timestep 0.002 nprint 50
iprfrq 500 - restart
firstt 278.0 finalt 278.0 - !
start firstt 278.0 finalt 278.0
- - CPT hoover -
Pcons Pref 1.0 Pmass
500.0 Pgamma 25.0 -
refT 278.0 Tmass 1000.0 - -
inbfrq -1 imgfrq 50 -
ntrfrq 100 - iunrea 20 -
iunwri 40 - iuncrd
30 nsavc 50 - iunvel -1
nsavv 0
21
open write card unit 1 name 2adp990wat.crd write
coor card unit 1 After Nps coor dist sele
segid ADP1 .and. type CA end - sele
segid ADP2 .and. type CA end Stop _______________
______________________________________________
22
(No Transcript)
23
(No Transcript)
24
ANALYZE the TRAJECTORY file analyze.inp open
unit 34 read unfo name adp5freep.trj traj iread
34 nread 1 skip 50 set 5 1 label loop traj read
coor dipo sele segid ADP1 end coor rgyr sele
segid ADP1 end inte sele segid ADP1 end sele
segid BULK end incr 5 by 1 if 5 lt 401 goto
loop stop close unit 34
25
File analyze !/bin/bash STYPE2adp278nowat DPAT
H/caracas/voichi/thermprot/2adp278ana FPATHcarac
as-scratch/voichi/2adp278snowat max201 icnt100
max2 max1501 icnt1001 icnt1376 while
"icnt" ! "max" do cp FPATH/STYPE.icnt.dc
d 2adp278nowat.dcd /bio/terry3/charmm/c30a2x1/exec
/gnu.serial/charmm lt analyzenowat.inp gt
analyze.out /bio/terry3/charmm/c30a2x1/exec/gnu.s
erial/charmm lt phipsicalc.inp gt analyze.out grep
"" fort.11 gtgt psi1.check.out grep "" fort.10 gtgt
phi1.check.out grep -A 3 -f rotmat.grep
analyze.out gtgt rotmat1.out grep -f
surf1apol.grep analyze.out awk -f surfapol.awk
gtgt partsurfapol.ADP1.out grep -f surf2apol.grep
analyze.out awk -f surfapol.awk gtgt
partsurfapol.ADP2.out grep -f surfpol.grep
analyze.out awk -f surfpol.awk gtgt
surfpol.ADP2.out grep -f distcm.grep analyze.out
awk -f distcm.awk gtgt distcm.out grep -f
rmsd.grep analyze.out awk -f rmsd.awk gtgt
rmsd1.out grep -f rotangle.grep analyze.out
awk -f rotangle.awk gtgt rotangle1.out grep -f
dipo.grep analyze.out awk -f dipo.awk gtgt
dipo22.all.out grep -f rgyr.grep analyze.out
awk -f rgyr.awk gtgt rgyrADP1ADP2.out grep -f
rotangle.grep analyze2.out awk -f rotaxis.awk
gtgt rotaxis.out grep -f ener.grep analyze.out
awk -f ener.awk gtgt tote.alawat.out grep -f
ener.grep analyze.out awk -f ener.awk gtgt
tote12.out grep -f vdw.grep analyze.out awk -f
vdw.awk gtgt vdw12.out
26
CONTINUE ANALYZE (note useful to know grep and
awk) grep -f elec.grep analyze.out awk -f
elec.awk gtgt vdwelec.alawat.out grep -f elec.grep
analyze.out awk -f elec.awk gtgt elec12.out grep
-f cm.grep analyze.out awk -f cm.awk gtgt
cm1.out grep -f end2enddist.grep analyze.out
awk -f end2enddist.awk gtgt endtoenddist2.out grep
-f x.grep analyze.out awk -f file.awk gtgt
xcl1.out grep -f y.grep analyze.out awk -f
file.awk gtgt ycl1.out grep -f z.grep analyze.out
awk -f file.awk gtgt zcl1.out grep -f x.grep
analyze.out awk -f file.awk gtgt x2ca.out grep
-f y.grep analyze.out awk -f file.awk gtgt
y2ca.out grep -f z.grep analyze.out awk -f
file.awk gtgt z2ca.out grep -f x1.grep analyze.out
awk -f file.awk gtgt x1ca.out grep -f y1.grep
analyze.out awk -f file.awk gtgt y1ca.out grep
-f z1.grep analyze.out awk -f file.awk gtgt
z1ca.out icnt((icnt 1)) done
27
Some Results from Analysis of Trajectories -
Potentials of Mean Force - Conformational Space
of alad-alad system
28
Running CHARMM Under Unix General syntax
(assuming /bin/csh) gt charmm lt filename.inp gt
filename.out paramvalue ...
filename.inp - a text file containing CHARMM
input commands. filename.out - the log file
for the CHARMM run, containing echoed
commands, and various amounts of command
output. The output print level may be increased
or decreased, and procedures such as
minimization and dynamics have printout frequency
specifications. - the optional ampersand
will place the program in the background
under most Unix shells.
29
Run Multiple Dynamics Trajectories
runDyna !/bin/bash STYPE2adp5free DPATH/carac
as/voichi/thermprot/2adp278 FPATH/pucc/fortress/v
oichi/2adp278sml max2001 icnt1 while "icnt"
! "max" do /bio/terry3/charmm/c30a2x1/exec/gnu
.serial/charmm lt STYPE.dynam.inp gt
STYPE.dynam.out cp DPATH/STYPE.rst
DPATH/trj/STYPE.rst.icnt mv DPATH/STYPE.rst
DPATH/STYPE.rst.in cp STYPE.dynam.out
DPATH/trj/dynam.STYPE.out.icnt cp
DPATH/STYPE.trj DPATH/trj/STYPE.trj.icnt gzip
DPATH/trj/dynam.STYPE.out.icnt gzip
DPATH/trj/STYPE.rst.icnt mkdir
STYPEtrjicnt mv DPATH/trj/dynam.STYPE.out.i
cnt.gz STYPEtrjicnt mv DPATH/trj/STYPE.rst.
icnt.gz STYPEtrjicnt mv DPATH/trj/STYPE.trj.
icnt STYPEtrjicnt tar -cf STYPEtrjicnt.ta
r STYPEtrjicnt rm -r STYPEtrjicnt mv
STYPEtrjicnt.tar FPATH/. icnt((icnt
1)) done
30
Submit as batch jobs the rcac ITAP computers
have a batch job queueing system which should be
used for running long MD simulations with CHARMM.
CHARMM calculations should be limited to around
12 to 16 hours, to both promote resource sharing
and to minimize data loss due to machine or
network failures. To accomplish this, the user
has to carefully choose the number of time steps
for the calculations of dynamics trajectories.
Trajectories are stored in the .dcd (or .trj)
binary files and typically contain 20 to 40 ps
each. Example qsub q sp_huge nodes8,
walltime160000 rundyna
31
Interfacing to CHARMM A mechanism is provided
to allow users of CHARMM to write their own
special purpose subroutines which can be
incorporated into the system without threatening
its integrity. There are six "hooks" into CHARMM
which have been specially provided for casual
modifiers. For detailed descriptions of each of
these hooks, consult the routine in
/charmm/source/main/usersb.src on UNIX machines
or ...CHARMM.SOURCE.MAINUSERSB.SRC under
VAX/VMS.
32
RUNNING CHARMM at NIH CHARMM at NIH (Bernie
Brooks) http//www.lobos.nih.gov/Charmm/
Available on several computer systems at
NIH CHARMM at Pittsburg Supercomputing Center
? http//www.lobos.nih.gov/Charmm/chmdoc.html
NHLBI LBC Computational Biophysics SectionCHARMM
Documentation Index                             
                                                  
                                                  
                             
33
Use with caution! A tale from the literature J.
Chem. Phys. 2002, 117 8892-8897. An examination
of the five-site potential (TIP5P) for
water Martin Lísal, Ji í Kolafa, Ivo
Nezbeda Parameterization of the five-site model
(TIP5P) for water M. W. Mahoney and W.
L. Jorgensen, J. Chem. Phys. 112, 8910 (2000)
has been examined by several computer simulation
methods accounting properly for long-range
forces. ... It is shown that the simple spherical
cutoff method used in the original simulations to
find optimized parameters of this five-site model
yields results that differ from those obtained by
both the Ewald summation and reaction field
methods. Consequently, the pivot property to
which the parameters were adjusted, the location
of the density maximum at 1 atm, does not agree
with experimental values. The equilibrium
properties then show only a fair agreement with
experimental data and are uniformly inferior to
those of the four-site TIP4P water over the
entire coexistence range.
34
Use with caution J. Phys. Chem. B, 104 (15),
3668 -3675, 2000. Molecular Dynamics Simulations
of a Polyalanine Octapeptide under Ewald
Boundary Conditions Influence of Artificial
Periodicity on Peptide Conformation Wolfgang
Weber, Philippe H. Hünenberger, and J. Andrew
McCammon Ewald and related mesh methods are
nowadays routinely used in explicit-solvent
simulations of solvated biomolecules, although
they impose an artificial periodicity in systems
which are inherently nonperiodic. In the present
study, we investigate the consequences of
this approximation for the conformational
equilibrium of a polyalanine octapeptide... We
report three explicit-solvent molecular dynamics
simulations of this peptide in cubic unit cells
of edges L 2, 3, and 4 nm, using the P3M
method... The initial configuration of the
peptide is helical. In the largest unit cell (L
4 nm), the helix unfolds quickly... By contrast,
in the two smaller unit cells (L 2 and 3 nm),
the helix remains stable during 2 ns. ... The
helical conformation is stabilized by artificial
periodicity relative to any other configuration
sampled during the trajectories. This artificial
stabilization is larger for smaller unit cells,
and is responsible for the absence of unfolding
in the two smaller unit cells... These results
suggest that artificial periodicity imposed by
the use of infinite periodic (Ewald) boundary
conditions in explicit-solvent simulations
of biomolecules may significantly perturb the
potentials of mean force for conformational equili
bria, and even in some cases invert the relative
stabilities of the folded and unfolded states.