2.2. Submitting and monitoring mcce job,
a simple example
Here is an example of calculating pKa values for Hen egg lysozyme. In
this example, mcce was installed under /usr/local/mcce-1.0/, and the directory
/usr/local/mcce-1.0/bin/ is in my excutable searching path. First make a
directory for this set of run. You get the structure file 4LZT from the rscb
data site with command "getpdb" included in MCCE package:
jmao@metis jmao]$ mkdir 4lzt
[jmao@metis jmao]$ cd 4lzt
[jmao@metis 4lzt]$ getpdb 4lzt
Inquiring the remote file 4LZT.pdb ...
Saving as 4LZT.pdb ...
Then copy a "run.prm" file to this directory and edit it:
[jmao@metis 4lzt]$ cp /usr/local/mcce/run.prm.quick run.prm
You will need to edit this "run.prm" file. In the edtitor, make first
string of line (INPDB) "4LZT.pdb", and turn flags for all 4 steps to "t".
You may want to run one step at a time and look at the output and run.log
messages for each step.
The Most modified section will look like this:
Most modified entries
Input and Output:
t step 1: pre-run, pdb->mcce pdb (DO_PREMCCE)
t step 2: make rotatmers (DO_ROTAMERS)
t step 3: do energy calculations (DO_ENERGY)
t step 4: monte carlo sampling (DO_MONTE)
MCCE can be invoked simply by running the executable file "mcce". As the
program will take hours to days to finish, you may want to run the program
at the background and redirect the output message to a log file using command
"mcce > run.log &".
[jmao@metis 4lzt]$ mcce > run.log&
While running, MCCE writes out the progress, warning messages and error
messages directed here into file "run.log". Check mcce run time messages
and serious errors will appear in it. The "debug.log" and "progress.log" provide extra warning messages and progress report, which can be safely ignored
In running 4LZT I noticed a line in the run.log file:
Error! premcce_confname(): Can't get conformer list of this residue NO3
This means a cofactor NO3 was not recognized by mcce. The program continues
anyway and treats all 4 atoms of NO3 with 0 charge and radii of 1.7 Angstroms.
You may want to stop the program and make a TPL (residue paramter file) file,
but we will ignore the error. If NO3 is on the surface, it won't cause much
error in pKa becuase of 1) the screening effect of the solvent and 2) Poisson-Boltzmann
equation will take care of the salt effect if this surface NO3 is stripped
off like surface water by the mcce program.
After step 1, the file "step1_out.pdb" is a mcce format pbd file with
terminal residues split off to form an independently titratable NTR and CTR,
groups renamed based in file "name.txt", and surface water and salt molecules
with 10% solvent accessibility stripped off. Undefined cofactors will be
identified as made up of non-charged atoms. A temporary parameter file "new.tpl"
will be written out. This can help you star to make a tpl file for new groups.
The file "acc.res" lists the solvent accessibility of residues and the file
"acc.atm" lists the solvent accessibility of atoms. No conformers are made
for groups with more than 30% surface accessibility in the configuration
After step 2, you will get a multi-conformation pdb file "step2_out.pdb",
which is the only file used by step 3. The rotamer making statistics are
dynamically updated in file "rot_stat".
The DelPhi runs in step 3 are the most time consuming step in mcce. If
you look at the bottom of run.log when DelPhi is running it will look like
[jmao@metis 4lzt]$ tail run.log
Creating connectivity table...
Preparing DelPhi runs ...
3 focusing runs required for this protein.
Computing pairwise from conformer 1 to 281 of 281 total conformers
see progress.log for progress...
So, we take a look at the additional progress file "progress.log":
[jmao@metis 4lzt]$ tail progress.log
Doing pairwise 72 of 281 conformers. 5 seconds
Doing pairwise 73 of 281 conformers. 2 seconds
Doing pairwise 74 of 281 conformers. 4 seconds
Doing pairwise 75 of 281 conformers. 2 seconds
Doing pairwise 76 of 281 conformers. 3 seconds
Doing pairwise 77 of 281 conformers. 3 seconds
Doing pairwise 78 of 281 conformers. 4 seconds
Doing pairwise 79 of 281 conformers. 4 seconds
Doing pairwise 80 of 281 conformers. 5 seconds
Doing pairwise 81 of 281 conformers.[jmao@4lzt]$
This progress is reported in "progress.log" and it can be used to estimate
how long step 3 will take. A pairwise calculation and a reaction field energy
calculation will be performed for each conformer. So step 3 will take about
4 seconds*281 conformers*2 DelPhi's= 2248 seconds or 37 minutes. The output
directory of step 3 is "energies" and the output file is "head3.lst".Both
are needed as input of step 4.
Step 4 is Monte Carlo sampling, which simulates a pH or Eh titration.
It writes the progress to file "mc_out". The occupancy of each conformer
at each pH is in fort.38. The net charge of each ionizable residue is in
"sum_crg.out". The pKa or Ems are in "pK.out". These are obtained by non-linear
fitting of the titration curves in file "sum_crg.out".
One reason to run simulations is to be able to ask WHY. To understand
what contributes to the ionization energy of a residue at a specific pH,
use the tool "mfe.py". This program approximates the interactions between
residues by Mean Field approach. To look at the free energy of the reaction
GLU -> GLU- + H+ for Glu 58 at pH 5.078, the midpoint of the titration
curve of this residue run.
[jmao@metis 4lzt]$ mfe.py GLU-_0035_ 5.078 0.2
Residue GLU-_0035_ pKa/Em=5.078
Terms pH meV Kcal
vdw0 0.03 1.98 0.05
vdw1 0.08 4.91 0.12
tors -0.00 -0.00 -0.00
ebkb -1.39 -80.41 -1.89
dsol 2.72 157.67 3.71
offset -0.99 -57.45 -1.35
pH&pK0 -0.33 -19.04 -0.45
Eh&Em0 0.00 0.00 0.00
-TS 0.11 6.35 0.15
residues -0.15 -8.47 -0.20
TOTAL 0.10 5.54 0.13
ASP_0048_ 0.22 13.02 0.31
ASP_0052_ 0.95 54.90 1.29
ARG_0112_ -0.29 -16.64 -0.39
ARG_0114_ -0.24 -13.93 -0.33
The program mfe.py is bundled with the mcce package. The syntax of the
first argument about how to refer to residue should be taken from pK.out.
The second argument means the Boltzmann distribution of conformers for all
residues at 5.078 will be used. The Last argument says only specific interactions
with individual residues over 0.2 kcal/mol will be printed out individually.
The result explains the free energy of the reaction GLU -> GLU+
at pH 5.078, the midpoint of the titration curve of this residue, where reactant
and product are equally favored. So it is expected here that the ionization
energy in line TOTAL is close to 0. The small error is caused by mean field
approximation and this error may indicate the degree of coupled ionization
or conformational changes with other residues. The breakdown of ionization
energy helps identifying the factors controlling the pKa of the residue.
Negative number favors ionization, so ASP52 in this protein shifts up the
pKa of GLU35 by 0.95 pH unit.
Listed in the table are other energy terms. All these terms are the difference
between the contributions from the disignated term to the ionized conformers
and netural conformers. "vdw0" is van der Waals interaction of the side chain
of this residue. "vdw1" is the van der Waals interaction of the side chain
of this residue to the backbone atoms. tors is the torsion energy of this
residue. "ebkb" is the electrostatic interaction of the side chain to the
backbone. "dsolv" is the loss of solvation energy. "offset" is a correction
for each residue type defined in "extra.tpl". "pH&pK0" is the pH effect,
which equals solution pH minus solution pKa of that residue if the residue
is an acid or solution pKa of that residue minus solution pH if the residue
is a base. "Eh&Em0" is the redox potential effect, which is the current
solution potential minus the standard redox potential of that cofactor. "-TS"
is entropy effect. As the protonated and deprotonated residue have have different
constraints in a protein, the entropy of the two ionization states could
be different. "residues" refers to the mean field interaction with residues
given their conformations at this pH. This is the term that is pH or Eh
A mcce web interface analysis tool is provided separately. If it is installed,
just point your browser to "http://localhost/cgi-bin/mcce.py"
(replace localhost by the host name of mcce server if you run mcce on a remote
computer). For more information, refer to section 4, Understanding MCCE Results.
<Prev | Contents | Next >