ICM Manual v.3.8
by Ruben Abagyan,Eugene Raush and Max Totrov
Jun 20 2021

 Prev ICM Language Reference Energetics and electrostatics Next

# How to plot the distance dependence of a van der Waals interaction

The following script will plot three energy-interatomic distance plots for three possible van der Waals terms defined by the vwMethod preference ("exact"-black,"soft"-blue and "old soft"-red). Simply mark and paste the following lines into your ICM session:
```
build string "o;o"  # two oxygens
set term "vw" only
set a_2//o Sum(Xyz(a_1//o ))+{.0 .0001 .0}
n=200
a=Rarray(n)
b=a
c=a
r=Rarray( n .03 3.)
vwSoftMaxEnergy = 14.5
for i=1,n
translate a_2 add {0.03 0. 0.}
#
vwMethod="exact"
r[i]=Distance( a_2//o  a_1//o ) # use Sum(Distance(..)) for the old version
show ey mute
a[i]=Energy("vw")
#
vwMethod="old soft"
show ey mute
b[i]=Energy("vw")

vwMethod="soft"
show ey mute
c[i]=Energy("vw")
endfor
s=Sarray(3*n)
s[n+1]="_red line"
s[2*n+1]="_blue line"
plot ds r//r//r a//b//c s {0. 6.01 1. 1. ,-1. 17. 1. 1.}
```

# How to calculate the electrostatic free energy by the REBEL-method

This short script solves the Poisson equation by the "Rapid-Exact- Boundary ELement (REBEL) method for crambin.
Examples:
```
electroMethod="boundary element"
delete a_w*       # get rid of water molecules
show energy "el"
show Energy("el")-r_out, r_out # Coulomb and solvation components
```

To extract the surface polarization charge per atom use the Rarray (a_//*) function, e.g.

```
electroMethod = 4  # REBEL
show energy "el"
Rarray( a_//* )    # returns polarization charges
# if you want to correct the partial charges induced polarization
set charge a_//* Charge(a_//*) - Rarray( a_//* )
```

# How to evaluate the pK shift

Suppose we mutate a surface Asp into Ala and want to evaluate how the pK of the neighboring His is changed. The pK shift may be evaluated as the difference of potential at Nd1 His and Ne2 His nitrogens due to the mutation. Since Ala may be considered as uncharged, the shift is simply the potential at the nitrogens due to the Asp charge.
Example (pKshift of His34 from Asp22):
```
# we assume that the positive charge is
# equally distributed between the two nitrogens
make boundary
pKshift=Sum({0.5, 0.5}*Potential(a_/his/nd1,ne2 a_/22/*))/ \
(Boltzmann*300.*Log(10.))
```

# How to evaluate the binding energy

There are many different approaches to the evaluation of binding energy. One of the reasonable approximations has the following features:
• van der Waals/hydrogen bonding interaction is excluded since it has close magnitudes for protein-protein and for protein-solvent interactions;
• electrostatic free energy change is calculated by the REBEL method (see also the section "How to calculate the electrostatic free energy ... ") above);
• side-chain entropy change is calculated by standard ICM entropic term based on exposed surface area of flexible side-chains;
• hydrophobic energy change is calculated using surface term with constant surface tension of 20. cal/Angstrom.

Example:
```
electroMethod="boundary element"
surfaceMethod="constant tension"
surfaceTension=0.020
dielConst = 12.7
set terms "sf,el,en"
show energy a_1 a_1 mute
e1  =Energy("el,sf,en")
show energy a_2 a_2 mute
e2  =Energy("el,sf,en")
show energy mute
e12 =Energy("el,sf,en")
print "Binding energy = ", e12 - e1 - e2
```

# How to calculate an ensemble average

The following is an example of calculating the average of an inter-atomic distance over a set of conformations collected in the conformational stack. This calculation is written as a macro. Feel free to change it. You may also use trajectory file and load frame instead of stack and load conf, respectively.
```
# first, define the macro
macro ensembleAverage r_temperature
l_commands = no
l_info = no

load conf 0                      # extract the lowest energy
e0 = Energy("func")

ansAver = 0.                     # the statistical sum initialization
statsum = 0.
r_temperature = r_temperature * Boltzmann

for i = 1,Nof(conf)              # loop through all the stack
# conformations
prob = Exp((e0-Energy("func"))/r_temperature)

# averaging distance between two ca
# is just an example
ansAver = ansAver + Distance(a_/2/ca a_/4/ca)*prob
statsum = statsum+ prob
endfor
r_out = ansAver/statsum
print " Ensemble average is: ", r_out
endmacro
# Now you can calculate your average

read stack                    # stack file is assumed to have
# the same name
ensembleAverage 600.          # sometimes you may use the elevated
# temperature to account for relaxation
```

# How to evaluate helicity of a peptide from the BPMC simulation

1. Run the _folding script first. Make sure the procedure converges by running several simulations (say _f1 _f2 _f3) from different random starting conformations. E.g.
```
cp \$ICMHOME/_folding _f1      # adjust the script
icm _f1 > f1.ou &
cp _f1 _f2
icm _f2 > f2.ou &
cp _f2 _f3
icm _f3 > f3.ou &
```
2. You can evaluate helicity for each simulation. If they converge the result will be about the same.
3. Helicity is just the ensemble average of the parameter which can be calculated as the relative number of the helical residues. Therefore you need to assign secondary structure for a particular trajectory frame or stack conformation and count number of helical residues. See macro _helicity averaging helicity over the trajectory frames.
```
macro helicity s_trjFileName r_temperature
# attention: 'temperature' is extremely important.
# You may use elevated temperature to account for relaxation.
l_commands=no
l_info=no
e0=Energy("func") # the lowest energy
av=0.
ssum = 0.
r_temperature = r_temperature * Boltzmann
res = Real(Nof(a_/*))
for i=1,Nof(frame)
assign sstructure
prob = Exp((e0-Energy("func"))/r_temperature)
av = av + prob*Nof(Sstructure(a_1/*),"H")/res
ssum= ssum+prob
endfor
print " The best E=", e0, "  Helicity= " av*100./ssum
endmacro
```

# How to merge and compress several conformational stacks

You may run several montecarlo simulations and accumulate several conformational stacks ( *.cnf files). To unite them it is essential that they have been created with the same energy function, because the compression algorithm takes the energy into account to decide which structure is more valuable. If it is not the case, you can always recalculate energies for the stack conformations by the following procedure:
```
read stack "f2" append # the second stack will be appended
for i=1,Nof(conf)
show energy s_correctTerms  # say, "vw,14,to,el,sf,en"
store conf i
endfor
```

Now, to unite all the stacks and compress them you may do the following (just the idea):

```
delete stack
read stack "f1"        # first simulation
read stack "f2" append # second simulation
read stack "f3" append # third simulation
show stack             # look at what you have now
compare v_//phi,psi  # use the comparison criterion from the simulation script