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

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.}

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_//* )

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.))

• 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.

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

# 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

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

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