ICM Manual v.3.8
by Ruben Abagyan,Eugene Raush and Max Totrov
Copyright © 2018, Molsoft LLC
Oct 18 2018

Contents
 
Introduction
Reference Guide
Command Line User's Guide
 ICM-shell
 ICM graphics
 Str.Analysis
 Sequence
 Molcart
 Pharmacophores
 Energy
  Plotting vw faq
  REBEL faq
  PK shift
  Binding energy
  Ensemble average
  Helicity
  Stacks merge
 Molecules
 Animation
 Symmetry
 X-ray
 Plotting
 Docking/VLS
 Examples
 _chemSuper
 Chemical Conformation Generator
References
Glossary
 
Index
PrevICM Language Reference
Energetics and electrostatics
Next

[ Plotting vw faq | REBEL faq | PK shift | Binding energy | Ensemble average | Helicity | Stacks merge ]

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" 
 read object s_icmhome+"crn" 
 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):
 
 read object s_icmhome+"rinsr" # load insulin  
                               # 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" 
 read object s_icmhome+"2ptc"  
 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 
    load conf i 
    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 object "my_peptide" 
 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 
     read conf 0 s_trjFileName 
     e0=Energy("func") # the lowest energy 
     av=0. 
     ssum = 0. 
     r_temperature = r_temperature * Boltzmann 
     res = Real(Nof(a_/*)) 
     read trajectory s_trjFileName 
     for i=1,Nof(frame) 
        load frame i s_trjFileName 
        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 object "f" 
 read stack "f1" 
 read stack "f2" append # the second stack will be appended 
 for i=1,Nof(conf)  
   load conf i 
   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):

 
 read object "f"         
 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 
                      # compose a new one, see also the compare command 
 vicinity = 40.       # you may redefine the vicinity parameter 
 compress stack 
 show stack           # look at the compressed stack 
 write stack "f4"     # save the result 


Prev
Ph4 object faq
Home
Up
Next
Molecules

Copyright© 1989-2018, Molsoft,LLC - All Rights Reserved. Copyright© 1989-2018, Molsoft,LLC - All Rights Reserved. This document contains proprietary and confidential information of Molsoft, LLC. The content of this document may not be disclosed to third parties, copied or duplicated in any form, in whole or in part, without the prior written permission from Molsoft, LLC.