ICM Manual v.3.9
by Ruben Abagyan,Eugene Raush and Max Totrov
Jun 5 2024

 Prev ICM Language Reference How to plot Next

# How to make a simple plot y=F( x)

If your function is determined on an uneven set of abscissa values, stored in array x, then use:
```
plot display x y
```
If your function is determined on a regular 1D-grid of values varying from x0 to x1, and is stored in the array y, use:
```
plot display Rarray(Nof(y), x0, x1) y
```
Use on-line help:
```
help plot
```
or to see the concise list of arguments
```
help command plot
```

Another way to make a file with a plot is to use the make plot command, e.g.

```
add column t Random(100,0.,1.) Random(100,0.,1.) Random(100,0.,1.)
make plot t "x=A;y=B;color=C;size=6;style=dots;;" output="aaa.png" size={500,500} # or pdf,eps,jpg
```
The type of file will be determined by the extension.

# How to plot a histogram

Function Histogram returns the 2*N matrix which may be used on the fly with a command like
```
plot display BAR red Histogram(a 50)
```
More examples:
```
a=Random(0. 100. 10000)
b=.04*(Count(1 50)*Count(1 50))
annot = {"Histogram with quadratic ruler" "Random value" "N"}
plot display BAR red Histogram(a b) annot

b=Sqrt(100.*Count(1 100))
annot = {"Histogram with square root ruler" "Random value" "N"}
plot display green BAR Histogram(a b) annot
```

# How to make a 3D-surface plot of a 2D-function

In the most common case you may have a matrix N*M my_matrix representing a two-dimensional function of x and y real arrays. First create a graphics object, then display it and save as a picture if you wish:
```
# let's try a "wire" representation first
x=Rarray(Nof(ram),-180.,180.)  #this scale may be nonlinear
y=x
make grob "matwire" ram x y # scales for x,y. Def. z-scale=1.
display matwire
# now let's have a look at a solid surface
make grob solid "matsolid" ram x y 0.9 # scales for x,y and z
display matsolid solid
# adjust the view by rotating/translating
#   connect matsolid         # you may move them with respect to each other
# and save image
write image rgb "img_matsolid"
```

# How to create a new graphics object of a specific shape

First, use one of the provided graphics objects, stored in *.gro files to try these kind of ICM non-molecular objects. You can scale graphics objects by multiplying them by a real number and shifting (e.g. g1=10.*g + {1.,0.,0.}) and/or adjust their positions either by translate or rotate commands or by connecting to the selected grob and dragging/rotating them with the mouse in the graphics window, for example:
```
my_icos = g_icos*0.2
display my_icos g_icos
translate add my_icos {5. 15. 0.}
```
You may also create a new or edit an existing graphics objects (see provided *.gro files to learn about format).

# Flexible peptide docking

Let us consider a problem of flexible peptide docking to a receptor. The peptide may be completely free or constrained (e.g. have a helical fragment in the middle), and the receptor will be assumed to have a fixed conformation. We will be going through a series of interactive or non-interactive steps. Feel free to put the commands below into several scripts and run them sequentially.
We will follow a particular example through all the steps. In this example we will be docking a constrained 11-residue alpha-helical peptide to a Bcl-XL molecule. It is a good idea to create a separate directory for the project, say, projectX , to keep all the files in one place.
Preparing the receptor binding surface for peptide docking.
First, we need to get the receptor model as an ICM object. It may come from several sources. The most common one is a pdb-structure converted with the convertObject macro, or the build model command which builds a model by homology, followed by regularization and refinement with the refineModel macro.
We will create two objects:
 a_1.|the original pdb object with receptor with some other peptide a_2.|converted receptor (without the peptide)

Example:
```
% cd projectX
% icmgl -G
nice "1g5j"  # an open form of the receptor
delete a_w*  # delete waters but keep bound peptide
copy a_ "rec"
delete a_rec.2  # keep only the main chain
undisplay window
convertObject a_rec. yes yes yes no # also optimizes hydrogens
#
display box Box( a_1.2 4. )  # display a box around the peptide
```
Now you may rotate the scene and adjust the box size by pulling a box corner with the leftMouseClick. The box needs to cover the receptor surface you need to include in the simulation, but needs to be minimized to reduce the computational burden and memory use.
Sometimes the orientation of the box with respect to the receptor surface is not favorable. In that case, use connect a_*. and rotate the objects with the respect to the origin ( use axisLength = 15. and display origin to to see the axis in addition to seeing the box ). If you change the absolute position of the receptor model, you may want to save the new orientation of the pdb object (e.g. write pdb a_1. "template" )
The result Now we have two objects: (i) a pdb object in which we deleted the unnecessary molecules, but kept the peptide for comparison and convenience in the box definition; and (ii) an ICM-object with a receptor for which we will be calculating the grid maps.
Create grid potential maps
Now calculate six grid potentials with the make map potential command:
```
make map potential "gc,gh,ge,gs,gb" Box() 0.5 # additional 'gl' grid is made for 'gc'
```
This command uses all molecules in the current object as a source for the map calculations. The Box() function will return the box you see on your screen. Alternatively you can just provide the 6 numbers defining that box.
0.5 is the grid size. Even though the ICM grid potential is extremely fast compared to other programs, for large boxes you may need to increase the grid cell size. We do not advise to go beyond 0.7A, however.
The result of this step is five grid maps around the docking surface of interest.
Building the peptide
The peptide of interest can be build with the following simple command:
```
build string "ECLKRIGDELD"   # peptide sequence from Bax
rename a_ "pep"
set a_/* "HHHHHHHHHHH"
superimpose a_1.2/310:320/ca a_/1:11/ca align
```

Saving the results
Let us now save both the objects and the maps for future use.
```
write m_gc m_gh m_ge m_gs m_gb # gc.map .. files are created
write object a_*. "all"  # write 3 objects in one multi-object file
```

Running a simple simulation from one starting conformation
The actual simulation does not need to be run interactively. An ICM script file _dockpep is described below. Run the script by
```
% _dockpep >! f1.out  # the output file
```
The file contains the following commands:
```
#!/usr/local/bin/icmng
read object "all"  # contains a_
set object a_pep.  # the peptide is the current object
fix v_/3:10/phi,psi,omg       # fix some backbone angles
nvar = Nof( v_//phi,psi,H,P ) # number of essential variables
# let us choose iteration limits depending on nvar
mncallsMC    = 10000 + Integer(0.008*nvar*nvar*nvar*nvar*nvar)
mncalls      = 170+nvar*3
print mncallsMC mncalls
# make sure that these two numbers are not insane
# or set smaller mncallsMC and mncalls (e.g. 10000 and 100) as a test
temperature  = 600   # optimal temperature for the simulation
tolGrad      = 0.001 # exit minimization when gradient is < 0.001
mnhighEnergy = 15 + nvar/2   #
set vrestraints a_/* # load preferred backbone and side-chain angle zones
#
vicinity = 2.5
compare a_pep. static  # use sRmsd of the peptide
set terms "vw,14,hb,el,to,en,gc,gh,ge,gb,gs"
minimizeMethod = "conjugate"
#
# maps
GRID.margin = Distance( a_pep./1/ca a_pep./11/ca )
#
# run montecarlo from a single start
# the reverse option will make more reasonables moves
#
# impose tethers ?
montecarlo reverse trajectory

# run montecarlo from multiple starts and merge the stacks
for i=1,Nof(<starting conditions>)
montecarlo reverse trajectory append
endfor
```

This simulation created several files:
 f1.cnf a stack file with mnconf (or less) best non-redundant conformations f1.trj a trajectory file with all accepted conformations (if you used the trajectory option) f1.ou|the file with the text output of the script
The output file f1.ou contains information about every accepted conformation of the simulation, e.g.
```
...
__ __      300   5   asp BPMC   db   da   da    -32.42    -23.66  100   1.29 97810
DY Visi    300   1   phe BPMC   fb   fb   fb    -32.42    -32.42  101   0.57 97911
__ __      300   4   ser BPMC   sa   sg   sd    -32.42    -24.45   88   0.06 97999
DY Visi    300   5   asp BPMC  dmn  dmn  dmn    -32.42    -32.43   49   0.08 98048
...
```

Analyzing the results of a single simulation
Energy profile: Is the simulation long enough
The energies can be plotted to analyze the progression of the global energy optimization.
```
plotEnergy "f1.ou" 50.
plotBestEnergy "f1.ou" 50. "display append"
```
The first macro plots energies of each conformation generated and accepted in the course of the ICM stochastic global energy optimization. Due to the nature of the procedure, the energy may go up and down. The second macro plots (and appends the plot to the previously generated def.eps file) only the lowest energy achieved by a given iteration. This plot shows if and when a better conformation is found. Naturally, the energy can not be improved after the global minimum is found.
However, having a long flat plateau does not guarantee that the global minimum is found since a new drop of energy may come after a very long simulation. The time required for finding the global minimum in general depends on the number of variables, number of atoms (the cost of a single energy calculation) and complexity of the potential energy surface.
The best empirical way to make sure that the minimum is found is to compare the lowest energies found from completely different starting conformations (read about multiple start simulations below).
We use a window of 50 kcals/mole from the lowest energy. It helps to avoid the dependence of the scale on the initial energy which may be quite high.
To analyze the improvements only use plotBestEnergy .
Graphical analysis of the results
You can generate all stack conformations as independent ICM objects with the mkStackConf macro,e.g.
```
read object "all"  # three objects, the peptide is the current object
mkStackConf 1 10   # make ICM objects from 10 best energy structures
```
The objects will be called a_pep1. , a_pep2. , etc. Now these objects can be displayed simultaneously or separately and further analyzed, e.g.
```
display a_pep*.
color a_pep?*. Energy( stack ) # in case you extracted all confs
```
You can loop through the conformations interactively either with the original stack, e.g.
```
# create some view you like
for i=1,10 # or Nof(conf)
pause   # rotate and zoom, hit Return
endfor
```
, or use the objects you created with the mkStackConf macro.
Running a multiple start simulation
You can figure it out yourself :-) .
Convergence analysis
The key question we would like to ask is if at least two independent simulations from random starting conformations identified a nearly identical conformation with a similar energy as the lowest energy conformation. A low-tech example with just two simulations f1 and f2 :
```