How to use "gold list - better sampling around active sites"?

To make extra conformers around sites of special interest

1. Identify residues of interest

2. Save the corresponding grepped lines from a PDB file into a file that you will name
(Any line from a residue defines the whole residue.)

3. Amend run.prm

3.1 Turn rotamer making flag on:
t Use control file "head1.lst" for rotamer making (ROT_SPECIF), and
t Do rotate? (PACK)

3.2 If you want to limit conformer making for other residues, then amend the ROTATIONS line to:
1 number of rotamers in a bond rotation (ROTATIONS)

4. Run step 1 and 2 Step 1 will create head1.lst based on residues within 4Å to those in the will
have 12 rotation steps in head1.lst. Step 2 will then create rotamers based on head1.lst.

How to keep conformers as fixed?

If you have a group you want to fix always (eg Sidechains around an ion)

A. You can edit head1.lst for those groups so that they don't make rotate (R), swing (S) or H optimization (H) conformers. This fixes your choices latter and if you want to let these move (eg if you are looking at the apo protein) you have to run again from the beginning.


B. Keep all conformers. In step 4 you can control the run with head3.lst.
These are the first columns to the left of head3.lst
00001 SER01_0002_001 f 0.00

FL f means you do sampling; FL t means you fix the conformer (sorry)
To fix a particular conformer (eg the correct neutral His ligand to be the ligand to a metal)
set all conformers to
t 0.0 (these will be fixed off)
the conformer you want to
t 1.0

To fix a residue in it's ionized state.
set all neutral conformers to
t 0.0
keep all occupied conformers as
f 0.0

The most direct check of the resultant occupancy of each conformer is in fort.38
This is the raw conformer occupancy for each conformer at each pH or Eh
sum_crg.out and pK.out are derived from fort.38

How to use Hic-up website to make tpl files

1.Download the pdbdict file which contains connectivity of your cofactor.
  • If your cofactor is on the hic-up site, you can download the pdb-dict file for the cofactor in the form of (external link)
XXX and xxx are the residue name in capital case and little case respectively.
  • Example: If the cofactor is ADP. In unix you can run
curl (external link) > PDB_dict.txt
to download this file and save it as PDB_dict.txt

2.Convert the pdbdict file to tpl format
There is a tool in the bin directory of MCCE package named which does this job. Run PDB_dict.txt > adp.tpl
This converts the file into MCCE format and save it as adp.tpl
  • the first two lines of this output tells you the format of the CONNECT parameter
  • orbital type is guessed based on number of connected atoms, so you need to double check that.
  • all the atoms are placed in ADP01 conformer.
If you want to move them to another conformer (eg. BK), you need to do that before going forward.

3.Make other parameters using mk_iatom tool
Parameters including CONFLIST, NATOM, IATOM and ATOMNAME can be made automatically with mk_iatom tool based on CONNECT parameters. You need to avoid manually changing them. Run
mk_iatom adp.tpl adp.tpl
first adp.tpl is the input file name, second is the output name. You can name them the same, the output wouldn't be messed up.

4.Add other parameters
So far only the parameters related to the structure are taken care of. You still need to make parameters for CHARGE, RADIUS, PKA, EM etc.
Check (external link) for the rest of the making process.

How to make additional conformers for non-covalently bound ligands?

Some cofactors in the protein have more degrees of freedom than rotation around single bonds. Three more options are available to add rotamers.


1. Make sure the following two parameter are in the run.prm file:
1 Number of translation steps (N_TRANS)
0.5 Translation distance for each step (TRANS_DIST)

2. Add translation parameters into the tpl file
#ParaNam|Res |Atom|Param/toggle

Rotating the whole cofactor around an axis defined by two atoms

1. This is similar as the regular rotamer making step, using ROTATE and SWING subroutines. So make sure the parameters in run.prm are correctly configured. You may want to turn on SWING for the small movement of the cofactor.
2. In the tpl file, ROTAMER parameter is used for controling this rotamer making step. However, instead of list all the atoms being rotated

# Res # Axis Rotated_Atoms?

simply use WHOLE_CONF to tell the program to rotate the whole cofactor

# Res # Axis Rotated_Atoms?

Rotating the whole cofactor along a plane defined by three atoms

1. Same as the last option, this step also use ROTATE and SWING subroutine to make rotamers, so make sure run.prm is configured correctly.

2. SLIDE parameter is used to tell the program which three atoms are used to define the plane. The axis of rotation goes through the first atom, perpendicular to the plane. All atoms in the conformer will be rotated.

# Res #
SLIDE UbQ 0 C9 - C3 - C6

How to compile new Delphi (version 0.4, release1.1 Jan,2005) program?

Needed because the binary files written by DelPhi and read by MCCE have different structure with different compilers

In the unzipped Delphi directory(hestia: /home/jun/Qui_oxi/Delphi_test/newdelphi/newdelphi/), the makefile is the instruction file to compile all the files. It determined which compiler it is using. The default compilers for Linux SGI computers are pgf77 and pgcc. We use f77 and gcc instead. Inside the /src directory, the source files are placed here. An important file listed in makefile /src/ classifies the source files into different compiling types. Basically, .f, .c. and .F files will be used by different compiler and options. In the original folder, qdiff4v.f, timef.f and qinttot.f files were replaced by qdiff4v.F, timef.F and qinttot.F files due to the compiling error. The major problem in the compiling procedure is the adjustable array is not supported by FORTRAN 77 Compiler 7.5a. The compromising way is use a constant value instead of the parameter called by adjustable array. The following parameters are replaced: resnummax in wrtsit4.f file is replaced by 10000. iisitsf in wrtsit4 file is replaced by 1. natm in wrtsit4.f and rforce.f files are replaced by 100000. After this process, the Delphi compilation could be achieved. Another parameter, nbra(1000) is adjusted into nbra(5000) in the vwtms2.f file according current Delphi modification. When we change the parameter ifrm from .true. to .false. in getatm2.f file, the new compiled executable Delphi will be changed from formatted pdb Delphi to unformatted pdb Delphi.

How to do energy correction using correction list after step3?

Use to modify pairwise interactions between specific sites using other information (such as interaction energies from a QM calculation: eg Song, Y., Mao, J., and Gunner, M. R. (2006) Electrostatic environment of hemes in proteins: pKas of hydroxyl ligands, Biochemistry 45, 7949-7958

1.Extract energy.opp to a energies directory: zopp -x energies
2.Copy the head3.lst and correct.lst into the energy directory. (An example of correct.lst is as folloing) correction program: hestia: /home/jun/bin/opp_correct
4.The opp files of the energies directory has been updated. If you need updated energy.opp, you can create energies.opp based on the updated energies directory: zopp -c energies
An example of correct.lst:

HA3+W CUB2W 0.52285
HA3+H CUB2W 0.78746
HA3+W CUB2H -0.57868

How to add a low dielectric slab to mimic membrane?

Amend the IPECE (Implement Protein Environment for Continuum Electrostatics) key/value pairs in run.prm.
The IPECE routine will add a 33Å thick slab of neutral atoms and save it as a new chain labeled M:

t    add neutral atoms to simulate a membrane slab (IPECE_ADD_MEM)
33.  the thichness of the membrane to be add (IPECE_MEM_THICKNESS)
M    the chain ID of the add atoms (IPECE_MEM_CHAINID)

Step 1 will determine the membrane position and save its coordinates into the file mem_pos.
Step 2 will add the neutral atoms after all rotamers are added.

For each position of the protein in the slab, the burial of ionizable group is determined by minimizing
Σ(giSAi), where SA is the exposed surface area without the slab.
  • g is 1 for each terminal O or N of Asp, Glu, and Arg
  • 2 for Lys in the slab and 0 if it is exposed.
Thus after a first trial, some residues can be mutated to Lys for proper membrane positioning whenever they
should not be buried by the slab.

For IPECE citation and further details, please reference the following paper:
Y. Song, J. Mao, and M.R. Gunner, Calculation of proton transfers in bacteriorhodopsin bR and M intermediates, Biochemistry 2003, 42, 9875-88

How to add waters or ions in a specified protein cavity using IPECE

Rationale for using IPECE to add either ions or water:

Typically, IPECE is used when one has identified a cavity around a cluster of atoms, or the portion of a channel delineated by a cluster of atoms. These local features are the target of IPECE: the program will pack as many ions or waters (inclusively if all species have an ION_NAME line setup) it can in the defined volume, and can then provide information about possible binding sites or solvent interactions.

Installation: IPECE is not yet part of the MCCE software package distribution, but it is available for download upon request.

How IPECE runs: ipece uses a fixed column format parameter file similar to run.prm (called ipece.prm.addwat). This file has five key/value parameters to set when adding ions or waters:

  • Key Value

    LATTICE_SCALE 0.4 # = Minimal distance of separation for the added species on the placement grid.

    # Consequently, ION_NAME value(3) >= LATTICE_SCALE value.

    # If a value lower than this is chosen for ION_NAME value(3), then the LATTICE_

    # _SCALE value needs reducing until the above relation is true again, else

    # the 'sieve' is useless! All distances in angstroms.
    ADD_ION t # t = ion/water added, f = none added. The species in ION_NAME lines will be

    # added sequentially
    CAV_POS 89.515 117.352 36.414 # = Anchor Position: coordinates from the pdb file.
    CAV_THR 5.0 # = Cavity Threshold: approximate radius of protein cavity.
    ION_NAME O HOH X 1.2 1.2 # (1) col 33-42: PDB(col 13-22);
    ION_NAME _CL Y 1.2 1.0 # (2) col 49-56: Ion radius (used for probing);

    # (3) col 57-64: Distance separation for ion placement.


  • I. Prerequisite MCCE files: acc.atm and step1_out.pdb from MCCE-Step1;

  • II. Edit ipece.prm.addwat:

    • 1) Identify one atom to be use as an anchor for the volume of waters/ions that will reside in the cavity.
    • 2) Pull this atom's xyz coordinates from prot.pdb (the pdb file in your working directory): these are the values needed for the CAV_POS key. IPECE matches CAV_POS in prot.pdb and outputs the coordinates of the waters/ions it can fit into the defined cavity.
    • 3) Save ipece.prm.addwat into ipece.prm.

  • III. Run ipece:

    • USAGE: ipece ipece.prm step1_out.pdb ipece_output.pdb
      • ipece.prm: parameters file for ipece routine;
      • input and output pdb files (step1_out.pdb and ipece_output.pdb, respectively) can have any name.

  • IV. Rerun MCCE steps 1 to 4 after setting ipece_output.pdb as the prot.pdb in run.prm (e.g.: after soft-linking the output of ipece as prot.pdb).

How to use make use of head3.lst?

The file "head3.lst" contains the self energy of each conformer and the control flag for Step 4.
The flag is: "t" for fixed occupancy 0 or 1, or "f" for free to sample. The energy unit is in kcal/mol.
For example, the three following lines from head3.lst pertain to a LYS residue:

iConf CONFORMER FL occ crg Em0 pKa0 ne nH vdw0 vdw1 tors epol dsolv extra history

02925 LYS01_0756_001 f 0.00 0.000 0 0.00 0 0 -5.962 -0.910 3.876 -0.755 0.574 0.000 01O000M000t
02926 LYS+1_0756_002 f 0.00 1.000 0 10.80 0 1 -6.045 -1.114 3.876 -3.510 2.888 0.000 +1O000M000t
02927 LYS+1_0756_003 f 0.00 1.000 0 10.80 0 1 -6.045 -1.114 3.876 -3.510 2.888 0.000 +1O000M000t

Right now, it has two ionization states, neutral (01) and protonated (+1), and a total of three conformers. All of them have
the possibility to be selected in the Monte Carlo sampling unless the flag tells Step 4 otherwise.
The control flag of the neutral state, for instance, can be switched so that this conformer will be fixed into the stated occupancy
(the number following the flag). In the example below, the neutral Lysine conformer is eliminated from MC sampling, while the
third conformer is forced into its conformation with 75% occupancy:

02925 LYS01_0756_001 t 0.00 0.000 0 0.00 0 0 -5.962 -0.910 3.876 -0.755 0.574 0.000 01O000M000t
02926 LYS+1_0756_002 f 0.00 1.000 0 10.80 0 1 -6.045 -1.114 3.876 -3.510 2.888 0.000 +1O000M000t
02927 LYS+1_0756_003 t 0.75 1.000 0 10.80 0 1 -6.045 -1.114 3.876 -3.510 2.888 0.000 +1O000M000t

Now only the two conformers in the protonated (+1) state can be occupied in the MC sampling of Step 4.
head3.lst thus allows for the full control of the conformers (position + ionization states) of each residue and cofactor.

How to visualize surface potential from MCCE output?

First use the col# in fort.38 script to extract the most occupied conformers at a pH.
Then, use script to convert the mcce formatted pdb to the pqr format that apbs can read.
Finally in pymol (with apbs plugin), use the pqr file (which contains charge and radii) to generate the surface potential.

  • Note the function call to retrieve the most occupied conformers at pH7: 8 (i.e. pHx + 1)

How does MCCE interact with Delphi?

Code in the MCCE source file energies.c calls the delphi executable through a system command. All input and output to and from Delphi is stored in a temporary folder, the location of which can be set with the (PBE_FOLDER) option in run.prm. Additionally, although the default behavior of MCCE is to delete the temporary Delphi folder when it is done with it, by setting the (CLEAN_DELPHI) option to 'f' in the run.prm file, you can make MCCE keep this folder around for you to examine later.

Important notes


-the potentials output by Delphi are in units of kT/mol
-the energies in MCCE's .opp files are in unit of kcal/mol
-the energies in MCCE's final output files are in units of pH, i.e. the energy required to shift the pK of an acidic residue up by one unit, or to shift the pK of a basic residue down by one unit

Unit conversions used by MCCE

-1 pH = 1.364 kcal/mol
- 1 kcal/mol = 1.688 kT/mol
- 1 kJ/mol = 0.239 kcal/mol

How to split step 3 into separate runs and combine them?

Step 3 in MCCE needs most time of the whole process. It is possible to separate this step into independent runs and runs off on multiple CPUs to save time. The mechanism is like the following example for illustrating purpose, and one can write script to automate it.

  • Inside a working directory, run step 1 and 2
  • Check the number of conformers in head2.lst, divide them into groups. For example, I have 2384 conformers, then I will divide them into 5 groups, with 500 conformers in each.
  • Make subdirectories subrun01, subrun02, subrun03, subrun04, subrun05
  • Inside each subrun directory, copy run.prm from its parent directory, soft link file step2_out.pdb, head2.lst. 
  • Inside each subrun directory, edit run.prm, to run step 3 only, and directive DELPHI_START and DELPHI_END to corresponding conformer numbers. For example, the first subrun directory should be:
    1         delphi start conformer number, 0 based            (DELPHI_START)
    500        delphi end conformer number, self included        (DELPHI_END)
  • Submit the jobs under each directory
  • Once finished, each subrun directory will have an energies.opp file. Go to working directory above the subrun directories, copy the first energies.opp file from subrun01, then run
    zopp -a subrun02 subrun03 subrun04 subrun05
    to merge all the runs.
  • To check if the results are merged successfully, view file head3.lst