[ Abs | Acc | Acos | Acosh | _Align | Angle | Area | Areacontactmatrix | Asin | Asinh | Ask | Atan | Atan2 | Atanh | Atom | Augment | Axis | Bfactor | Boltzmann | Box | Bracket | Cad | Ceil | Cell | Charge | Cluster | _Color | Consensus | Corr | Cos | Cosh | Count | Date | Deletion | Det | Disgeo | Distance | Eigen | Energy | Error | Exist | Existenv | Extension | Exp | Field | Fielduser | File | _Find | Floor | Getenv | Gradient | _Grob | Group | Histogram | Iarray | Iarrayatoms | nvis | IcmSequence | Index | _Indexx | Insertion | Info | Integer | Integral | Interrupt | Label | Length | LinearFit | Link | Log | Map | Mass | Matrix | Max | MaxHKL | Mean | Min | Money | Mod | Mol | Name | Namex | Next | Nextcovalentneighbors | Nof | Norm | Obj | Occupancy | Path | Parray | Pattern | Pi | Potential | Power | Probability | Profile | Putenv | Radius | Random | Rarray | Real | Remainder | Replace | Res | Resali | Resolution | Rfactor | Rfree | Rmsd | Rot | Sarray | Score | Select | Sequence | reversecompdnaseq | Sign | Sin | Sinh | Site | Smiles | Smooth | Split | Sql | Sqrt | Sphere | _Split | Srmsd | String | Sstructure | Sum | Symgroup | Table | alitable | Tablestack | Tan | Tanh | Tensor | Temperature | Time | Tolower | Torsion | Toupper | Tr123 | Tr321 | Trace | Trans | Transpose | Trim | Turn | Type | Unix | Value | Vector | Version | Volume | View | Warning | Xyz ]

 
 show Mean(Random(1.,3.,10))  

 
 a=Abs(-5.)                 # a=5.  
 print Abs({-2.,0.1,-3.})   # prints rarray {2., 0.1, 3.}  
 if (Abs({-3, 1}=={3 1}) print "ok" 

 
                                   # let us select interface residues  
 read object s_icmhome+"complex" 
                                   # display all surface residues 
 show surface area 
 display Acc( a_/* ) 
                                   # now let us show the interface residues 
 display a_1,2 
 color a_1 yellow 
 color a_2 blue 
 show surface area a_1 a_1         # calculate surface of  
                                   # the first molecule only  
 
                                   # select interface residues  
                                   # of the first molecule  
 color red Sphere(a_2/* a_1/* 4.) & Acc(a_1/*)      
 
 read object "crn" 
 show energy "sf" 
 display 
 display cpk Acc(a_//* 0.1)   # display accessible atoms  
 
 show surface area            # prior to invoking Acc function  
                              # surface area should be calculated  
 color Acc(a_/*) red          # color residues with relative  
                              # accessibility > 25% red  

 
 print Acos(1.)               # equal to 0.  
 print Acos(1)                # the same  
 
 print Acos({-1., 0., 1.})    # returns {180. 90. 0.}  

 
 print Acosh(1.)               # returns 0  
 print Acosh(1)                # the same  
 
 print Acosh({1., 10., 100.})  # returns {0., 2.993223, 5.298292}  

• i_out - the number of identical residues in the alingment
• r_out - contains Log( Probability_of_structural_dissimilarity ) only for pairwise alignments
• r_2out - percent identity of the alignment.

 
   read sequences s_icmhome+"sh3.seq" # read 3 sequences 
   print Align(Fyn,Spec)              # align two of them  
   Align( )                           # the first two 
   a=Align( sequence[1] sequence[3] ) # 1st and 3rd 
   if(r_out > 5.) print "Sequences are struct. related"  

 
 read pdb "1lbd" 
 show surface area  
 make sequence 
  Info> sequence  1lbd_m  extracted  
 1lbd_m                       # see the relative areas 
 read pdb sequence "1fm6.a/"  # does not have areas 
 Info> 1 sequence 1fm6_a  read from /data/pdb/fm/pdb1fm6.ent.Z    
 ali3d = Align( 1lbd_m 1fm6_a area )  

 
 set area seq1 R_weights  # must be > 0. and less than 2.37 
 Align( seq1 seq2 area ) 

 
 seq1 = Sequence("WEARSLTTGETGYIPSA") 
 seq2 = Sequence("WKVEVNDRQGFVPAAY") 
 Align() 
 # Consensus   W.#.  .~~.~G%#P^ 
           WEARSLTTGETGYIPS-- 
           WKVE--VNDRQGFVPAAY 
 m = Matrix(17,16,0.) 
 m[10,4] = 3. # reward alignment of E in seq1[10]  and E in seq2[4] 
 Align(seq1 seq2 m ) 
 # Consensus   W.#      E    ~G%#P^ 
           WEARSLTTGE----TGYIPS-- 
           WKV------EVNDRQGFVPAAY 

• sequences must be linked to 3D molecules to access the coordinate information;
• two 3D structures must have a superimposable subsets

 
             # 14 sequences  
 read alignment msf s_icmhome + "azurins"  
             # extract a pairwise alignment by names  
 aa = Align(azurins,Azu2_Metj,Azur_Alcde)  
             # reordered sub-alignment extracted by numbers  
 bb = Align(azurins,{2 5 3 4 10 11 12})    

 
 macro reorderAlignmentSeq( ali_ ) 
  nn=Name(ali_)  # names in the alignment order 
  ii=Iarray(Nof(nn)) 
  j=0 
  for i=1,Nof(sequence)  # the original order 
   ipos = Index( nn, Name(sequence[i] ) ) 
   if ipos >0 then 
     j=j+1 
     ii[j] = ipos 
   endif 
  endfor 
  ali_new = Align( ali_ ii ) 
  keep ali_new 
 endmacro 

 
 build string "se ala his leu gly trp ala" "a"  # obj. a  
 build string "se his val gly trp gly ala" "b"  # obj. b  
 set tether a_2./1:3 a_1./2:4 align                 # impose tethers  
 show Align(a_2.1)         # derive alignment from tethers  
 write Align(a_2.1) "aa"   # save it to a file  

 
 print Angle(  # and then click the atom of interest. 

 
 d=Angle( a_/4/c )                         # d equals N-Ca-C angle  
 print Angle( a_/4/ca a_/5/ca a_/6/ca )    # virtual Ca-Ca-Ca planar angle  

 
read object "1crn" 
show surface area 
a=Area(a_/* )      # absolute conformation dependent residue accessilities  
b=Area(a_/* type ) # maximal residue accessilities in the extended conformation 
c = a/b            # relative (normalized) accessibilities 

 
 read object "crn"  # good old crambin  
 s=String(Sequence(a_/A)) 
 PLOT.rainbowStyle="blue/rainbow/red" 
 plot area Area(a_/A, a_/A) comment=s//s color={-50.,50.} \ 
    link transparent={0., 2.} ds 
 
 read object "complex" 
 plot area Area(a_1/A, a_2/A) grid color={-50.,50.} \ 
    link transparent={0., 2.} ds 

 
 build                   # build a molecule according to the sequence 
                         # from file def.se (default) 
 show area surface       # calculate surface area  
 a = Area(a_//o*)        # individual accessibilities of oxygens  
 
 stdarea = Area("lys")   # standard accessibility of lysine  
 
# More curious example 
 read object "crn"           
 show energy "sf"        # calculate the surface energy contribution 
                         # (hence, the accessibilities are 
                         # also calculated) 
 
 assign sstructure a_/* "_"     
                         # remove current secondary structure assignment 
                         # for tube representation 
 display ribbon 
                         # calculate smoothed relative accessibilities  
                         # and color tube representation accordingly  
 color ribbon a_/* Smooth(Area(a_/*)/Area(a_/* type) 5) 
                         # plot residue accessibility profile 
 plot Count(1 Nof(a_/*)) Smooth(Area(a_/*)/Area(a_/* type) 5) display 

 
 print Asin(1.)               # equal to 90 degrees  
 print Asin(1)                # the same  
 
 print Asin({-1., 0., 1.})    # returns {-90., 0., 90.}  

 
 print Asinh(1.)               # returns 0.881374  
 print Asinh(1)                # the same 
 
 print Asinh({-1., 0., 1.})    # returns {-0.881374, 0., 0.881374} 

 
 windowSize=Ask("Enter window size",windowSize) 
 s_mask=Ask("Enter alignment mask","xxx----xxx") 
 
 grobName=Ask("Enter grob name","xxx") 
 display $grobName 
 
 show Ask("Enter string, it will be interpreted by ICM:", "")  
          #e.g. Consensus( myAlignm ) 
 
 show Ask("Enter string:", "As Is",simple)     
          #your input taken directly as a string 

 
 print Atan(1.)                    # equal to 45.  
 print Atan(1)                     # the same.  
 
 print Atan({-1., 0., 1.})         # returns {-45., 0., 45.}  

 
 print Atan2(1.,-1.)                        # equal to 135.  
 print Atan2({-1., 0., 1.},{-0.3, 1., 0.3}) # returns phases {-106.7 0. 73.3}  

 
 print Atanh(0.)                    # returns 0.  
 print Atanh(1.)                    # returns error  
 
 print Atanh({-0.9999, 0., .9999})  # returns { -4.951719, 0., -4.951719 }  

 
 asel=Acc(a_2/his)             # select accessible His residues of  
                               # the second molecule  
 show Atom(asel)               # show atoms of these residues 
 show Atom( v_//phi )          # carbonyl Cs 

 
 a1  a2  a3  | a4 
 a5  a6  a7  | a8 
 a9  a10 a11 | a12 
 ------------+---- 
 0.  0.  0.  | 1. 

 
 M_inv = Power(Augment(R_12direct),-1)) 
 R_12inv = Vector(M_inv) 

 
 display a__crn.        # load and display crambin: P21 group  
 obl = Augment(Cell( ))  # extract oblique matrix  
 A = obl[1:3,1]         # vectors A, B, C 
 B = obl[1:3,2] 
 C = obl[1:3,3] 
 g1=Grob("cell",Cell( )) # first cell 
 g2=g1+ (-A)            # second cell 
 display g1 g2          

• r_out rotation angle (degrees);
• r_2out helix rise;
• R_out 3-rarray with a point on the axis.

 
 color ribbon a_/ Trim(Bfactor( a_/ simple ),-0.5,3.)//-0.5//3. # or 
 color a_// Trim(Bfactor( a_// simple ),-0.5,3.)//-0.5//3.  # for atoms 

 
 avB=Min(Bfactor(a_//ca))     # minimal B-factor of Ca-atoms  
 show Bfactor(a_//!h*)        # array of B-factors of heavy atoms  
 color a_//* Bfactor(a_//*)   # color previously displayed atoms  
                              # according to their B-factor  
 color ribbon a_/A Bfactor(a_/A) # color the whole residue by mean B-fac. 

 
 deltaE = Boltzmann*temperature  # energy  

 
 build string "se ala his"  # a peptide  
 display box Box(a_/2 1.2)  # surround the a_/2 by a box with 1.2A margin 
 color a_//* & Box( )  

 
 make map potential "gh,gc,gb,ge,gs" a_1 Box() 
 m_ge = Bracket(m_ge, Box( a_1/15:18,33:47 )) # redefine m_ge 

[ Cad1 | Cadalign ]

• between two conformations of the same protein based on full atom residue-residue contact area calculation, Cad(..)
• between two conformations of the same protein based on Cbeta-Cbeta distance evaluation (`Cad1{Cad}(.. distance ) .ICM uses an empirically derived ContactStrength( Cb-distance ) function.
• between two homologous structures based preservation of the residue contacts through the alignment ( Cad (.. alignment )) . The contact strength in this case is also derived from the interresidue distances.

• by a single selection rs_A1 : all pairs between selected residues (is equivalent to rs_A1 rs_A1 )
• by two residue selections rs_A1 rs_A2: cross pairs between two sets of selected residues (e.g. the contacts between two subunits)

• distance option allows to compare approximations of the inter-residue contact areas by the Ca and Cbeta positions. This allows to calculated deformations between two homologous proteins which is not possible in the default mode in which two chemically identical molecules are compared. The residue pairs in two homologues are equivalenced according to the alignments linked to the molecules. Residues deleted in a homologue are considered to have zero contact.

 
# Ab initio structure prediction, Overall models by homology  
   read pdb "cnf1"   # one conformation of a protein  
   read pdb "cnf2"   # another conformation of the same protein 
   show 1.8*Cad(a_1. a_2.)   # CAD=0. - identical; =100. different 
   show 1.8*Cad(a_1.1 a_2.1) # CAD between the 1st molecules (domains) 
   show 1.8*Cad(a_1.1/2:10 a_2.1/2:10) # CAD in a window 
   PLOT.rainbowStyle = 2 
   plot area grid M_out comment=String(Sequence(a_1,2.1)) link display 
 
# Loop prediction: 0% - identical; ~100% totally different 
#    CAD for loop 10:20 and its interactions with the environment 
   show 1.8*Cad(a_1.1/10:20 a_1.1/* a_2.1/10:20 a_2.1/*)   
#    CAD for loop 10:20 itself 
   show 1.8*Cad(a_1.1/10:20 a_1.1/10:20 a_2.1/10:20 a_2.1/10:20)   
 
# Evaluation of docking solutions: 0% - identical; 100% totally different 
   read pdb "expr"   # one conformation of a complex  
   read pdb "pred"   # another conformation of the same complex 
   show Cad(a_1.1 a_1.2 a_2.1 a_2.2)  # CAD between two docking solutions  
# 
# ANOTHER EXAMPLE: the most changed contacts  
   read object "crn" 
   copy a_ "crn2" 
   randomize v_ 5. 
   Cad(a_1. a_2.) 
   show s_out 
   read column group input= s_out name="cont" 
   sort cont.1 
   show cont 
# the table looks like this (the diffs can be both + and -): 
#>T cont 
#>-1-----------2-----------3---------- 
   -39.        a_crn.m/38  a_crn.m/1 
   -36.4       a_crn.m/46  a_crn.m/4 
   -32.1       a_crn.m/46  a_crn.m/5 
   -29.8       a_crn.m/30  a_crn.m/9 
   -25.2       a_crn.m/37  a_crn.m/1 
... 
   42.5        a_crn.m/43  a_crn.m/5 
   45.1        a_crn.m/44  a_crn.m/6 
   45.2        a_crn.m/43  a_crn.m/6 
   55.3        a_crn.m/46  a_crn.m/7 
   56.         a_crn.m/45  a_crn.m/7 

• formal : return formal charges
• mmff : return formal charges calculated according to mmff atom types and rules. Note: do not confuse this option with a function to return the mmff charges.

 
 buildpep "ala his glu lys arg asp" 
 show Charge(a_1)      # charge per molecule 
 show Charge( a_1/* )  # charge per residue 
 show Charge( a_1//* ) # charge per atom  
 
 avC=Charge(a_/15)        # total electric charge of 15th residue  
 avC=Sum(Charge(a_/15/*)) # another way to calculate it 
 show Charge(a_//o*)      # array of oxygen charges  
 
# to return mmff charges: 
 set type mmff 
 set charge mmff 
 Charge( a_//* ) 
 
# to return total charges per molecular object: 
 read mol s_icmhome+"ex_mol.mol" 
 set type mmff 
 set charge mmff 
 Charge( a_*. ) 

 
# let us make a distance matrix D 
# we will cook it from 5 vectors {0. 0. 0.} 
 m=Matrix(5,3)  # initialize 5 vectors 
 m[2,1:3]={1. 0. 0.}   # v2 
 m[3,1:3]={1. 1. 0.}   # v3 
 m[4,1:3]={1. 1. 1.}   # v4 
 m[5,1:3]={1. 0.1 0.1} # v5 close to v2 
 
 D = Distance( m )  # 5x5 distance matrix created 
 
 Cluster( D , 0.2 ) # v2 and v5 are assigned to cluster 1 
 Cluster( D , 0.1 ) # radius too small. All items are singlets 
 Cluster( D , 2.  ) # radius too large. All items are in cluster 1 

 
 build string "se his" 
 display xstick 
 make grob image name="g_" 
 display g_ only smooth 
 M_clr = Color( g_ ) 
 for i=1,20    # shineStyle = "color" makes it disappear completely 
   color g_ (1.-i/20.)*M_clr 
 endfor 
 color g_ M_clr  

 
 read alignment "sx"  # load alignment  
 read pdb "x"         # structure  
 display ribbon 
        # multiply rs_ by a mask like "  A C   N  .." 
 cnrv = a_/A & Replace(Consensus(sx cd59),"[.^~#]"," ")  
 display cnrv red  
 display residue label cnrv 

 
 r=Corr(a,b)                   # two vectors a and b  
 if (Abs(r_out) < 0.3) print "it is actually as good as no correlation" 

 
 show Cos(60.)                           # returns 0.5  
 show Cos(60)                            # the same  
 
 rho={3.2 1.4 2.3}                       # structure factors  
 phi={60. 30. 180.}                      # phases  
 show rho phi rho*Cos(phi) rho*Sin(phi)  # show in columns rho, phi,  
                                         # Re, Im  

 
 show Cosh(1.)                           # 1.543081  
 show Cosh(1)                            # the same  
 
 show Cosh({-1., 0., 1.})                # returns {1.543081, 1., 1.543081}  

 
 show Count(-2,1)         # returns {-2,-1,0,1}  
 show Count(4)            # returns {1,2,3,4}  

 
Date()  # uses the default format 
Date("%m %DD,%Y %h%pm") # today returns 
 Oct 08,2002 5pm 

• "all" (default: no string option) select deletions of all types
• "nter" select only N-terminal fragments
• "cter" select only C-terminal fragments
• "loop" select only the internal zones of deleted loops

 
 a=Rot({0. 0. 1.}, 30.)     # Z-rotation matrix by 30 degrees  
 print Det(a)               # naturally, it is equal to 1.  

 
 group sequence se1 se2 se2 se4 mySeqs 
 align mySeqs 
 distMatr=Distances(mySeqs) 

 
 read sequences "zincFinger" # read sequences from the file, 
 list sequences              # see them, then ...  
 group sequence alZnFing     # group them, then ...  
 align alZnFing              # align them, then ...  
 a=Distance(alZnFing)        # a matrix of pairwise distances  
 n=Nof(a)                    # number of points  
 b=Disgeo(a)                 # calculate principal components  
 corMat=b[1:n,1:n-1]         # coordinate matrix [n,n-1] of n points  
 eigenV=b[1:n,n]             # vector with n sorted eigenvalues  
 xplot= corMat[1:n,1] 
 yplot= corMat[1:n,2] 
 plot xplot yplot CIRCLE display # call plot a 2D distribution  

[ Distanceiarray | Distancerarray | Distanceas_ | Distanceas_rarray | Distancematrix | DistanceTanimoto | Distancematrixtwosets | Distancetether | Distancedayhofffunction | Distancealignment | Distancetwoalignments ]

• Display your molecule
• type Dist , press TAB
• Ctrl-RightMB click on the atom you want (or double click for a residue) and press RETURN

 
  distWeighted = 1. - Sum( Wi_of_common_On_Bits ) / Sum( Wi_of_On_Bits ) 

 
Distance({3} {1} 32 simple ) # returns 0. 
Distance({1} {3} 32 simple ) # returns 0.5 

 
Distance({1 2 3},{1 2 3},32) 
 #>M 
 0. 1. 0.5 
 1. 0. 0.5 
 0.5 0.5 0. 

read pdb "1crn" convert tether # keeps tethers to the pdb original deviations = Distance( a_//!h*,vt* , 9.9) perResDevs = Group( deviations, a_//!h*,vt* ,"max") # find max.devs per residue display ribbon color ribbon a_/* perResDevs # Another example Distance( a_//T ) # selects only tethered atoms #>R 1.677 1.493 1.386 1.435 1.645 1.570 2.165 1.399

1. d1 = 1.0 - (nResidueIdentities/Min(Length(Seq1), Length(Seq2)) (d1 belongs to [0.,1.] range)
2. if there is no identity option the distance is corrected: Distance(Seq1,Seq2) = DayhoffTransformation( d1 )

 
 1.-(nResidueIdentities-gapPenalty)/Min(Length(Seq1), Length(Seq2)) 

 
 read alignment msf "azurins"            # read azurins.msf  
 NormCoord = Disgeo(Distance(azurins))   # 2D sequence diversity in  

 
 show a1   # La1 = 3 
 ABC---XYZ 
 ABCDEF--- 
 show a2   # La2 = 6 
 ABCXYZ 
 ABCDEF 
 Distance(a1,a2)   # a1 is a subalignment of a2, distance is 0. 
 0. 

 
 Distance(a1,a2,exact)  # returns 0.5 for the above a1 and a2 

 
 read sequence msf "azurins.msf"  
 gapOpen =2.2 
 a=Align(Azu2_Metj  Azup_Alcfa)   # the first alignment 
 gapOpen =1.9                     # smaller gap penalty and .. 
 b=Align(Azu2_Metj  Azup_Alcfa)   # the alignment changes 
 show 100*Distance(a b )          # 20% difference 
 show 100*Distance(a b exact )    # 21.7% difference 
 show a b   

 
 A = Matrix(3, 3, 0.)      # create a zero square matrix... 
 A[1:3,1] = {1.,-2.,-1.}   # and set its elements 
 A[2,2] =  4. 
 for i = 1, 3-1            # the matrix must be symmetric 
   for j = i+1, 3 
     A[i,j] = A[j,i] 
   endfor 
 endfor 
 X = Eigen(A)              # calculate eigenvectors... 
 V = R_out                 # and save eigenvalues in rarray V 
 printf "eigenvalue 1 eigenvalue 2 eigenvalue 3\n" 
 printf "%12.3f %12.3f %12.3f\n", V[1], V[2], V[3] 
 printf "eigenvector1 eigenvector2 eigenvector3\n" 
 for i = 1, 3 
   printf "%12.3f %12.3f %12.3f\n", X[i,1], X[i,2], X[i,3] 
 endfor 

• Allowed terms in the string are "vw,14,hb,el,to,af,bb,bs,cn,tz,rs,xr,sf";
• "func" stands for the total of all the terms, both energy and penalty;
• "ener" is only the energy part (i.e. "vw,14,hb,3l,to,af,bb,bs,sf" );
• "pnlt" is only the penalty part (i.e. "cn,tz,rs,xr" ).
• load conf and load frame commands fill out all the energy/penalty terms, which are stored in both stacks and movies (of course the values also depend on a set of free variables). You can get the energy/penalty terms of the loaded conformation without explicitly recalculating them.

 
 build 
 show energy 
 print Energy("vw,14,hb,el,to")  # ECEPP energy  
 
 read stack "f1" 
 load conf 0 
 print Energy("func")  # extract the best energy without recalculating it 

• electrostatic ( "el" ) term and electroMethod = "boundary element", "MIMEL", or "generalized Born"

 
  read pdb "1crn" 
  delete a_W 
  convert 
  set terms "vw,14,hb,el,to,en,sf" 
  group table t Energy( a_/A ) "energy" Label(a_/A ) "res" 
  show t 
  unfix V_//* 
  group table tBondsAngles Energy( a_/A "bs,bb" ) "covalent" Label(a_/A ) "res" 
  show tBondsAngles 

• sarray of the energy term names ( .hd ) and
• rarray of energy values for each energy term ( .ey ) and

 
 set terms only "vw,14,hb,el,to"  # set energy terms  
 show energy v_//xi*              # calculate energy with only  
                                  # side chain torsions unfixed  
           # energy depends on what variables are fixed since  
           # interactions inside rigid bodies are not calculated,  
           # and rigid body structure depends on variables  
 
 a = Energy("vw,14")      # a is equal to the sum of two terms  
 
 electroMethod="MIMEL"    # MIMEL electrostatics  
 set terms only "el,sf"   # set energy terms  
 show energy 
 print Energy("ener")     # total energy  
 print Energy("sf")       # only the surface part of the solvation energy  
 print Energy("el")       # electrostatic energy  
 print r_out              # electrostatic part of the solvation energy  

 
 read pdb "1mng"  # this file contains strange 28-th residue  
 if (Error) print "These alternative positions will kill me" 
 
 read pdb "1abcd"  # file does not exist 
 errorMessage = Error(string) 
 ier = i_out 
 if ier == 3102  print "What kind of name is that?" 
 
 if=2 
 if Error(2073) print "syntax error" # detect error [2073] 
 read pdb "1mok" 
 clear error 

 
 show 1.-Error(Sqrt(0.5)) # P of being inside +-sigma (about 68%) 
 show Error(2.*Sqrt(0.5)) # P of being outside +- 2 sigma 

 
 x=Rarray(1000 0. 5. ) 
 plot display x Error(x ) {0. 5. 1. 1. 0. 1. 0.1 0.2 } 
 plot display x Log(Error(x ),10.) {0. 5. 1. 1.}   
       #NB: can be approximated by a parabola 
       #to deduce the appr. inverse function. 
       #Used for the Seq.ID probabilities. 

 
 if (!Exist("/data/pdb/") then 
   unix mkdir /data/pdb 
 endif 
 
 if(!Exist(key,"Ctrl-B")) set key "Ctrl-B" "l_easyRotate=!l_easyRotate" 
 
 if !Exist(gui)  gui simple 

 
 if(Existenv("ICMPDB")) s_pdb=Getenv("ICMPDB") 

 
 print Extension("aaa.bbb.dd.eee")   # returns ".eee"  
 show Extension({"aa.bb","122.22"} dot)         # returns {"bb","22"} 
 read sarray "filelist"                                     
 if (Extension(filelist[4])==".pdb") read pdb filelist[4]   

 
 print Exp(deltaE/(Boltzmann*temperature))   # probability  
 print Exp({1. 2.})                          # returns { E, E squared }  

 
 Field("a b c",3," ")  # space 
 Field("a:b:c",3,":")  # colon 

 
 Field("a%b:c",3,"%:")  # percent OR colon 

 
 Field("a  b   c",3,"  ")    # two==multiple spaces in field delim 
 Field("a%b::::c",3,"%::")   # a single percent or multiple colons 

 
 s=Field("1 ener glu 1.5.",3)    # returns "glu"  
 show Field("aaa:bbb",2,":")     # returns "bbb"  
 show Field("aaa 12\nbbb 13","bbb",1) # returns "13"  
 show Field("aaa 12\nbbb 13 14","bbb",2,"  \n\n") # two spaces and two \n .  
# another example 
 read object s_icmhome+"all" 
       # energies from the object comments, the 1st field after 'vacuum' 
 show Rarray(Field(Namex(a_*.),"vacuum",1))  

 
 show Field({"a:b","d:e"},2,":") # returns {"b","e"}  
 s=Field({"aa 2 3.3", "bb 4 1.3", "cc 31a 1.1 3"},2) 
         # returns {"2","4","31a"}  
 s=Field({"aa 2 3.3", "bb 4 1.3", "cc 31a 1.1 3"},4)   
         # returns {"","","3"}  

 
set field a_//* Random(0.,1.,Nof(a_//*)) 
show Field( a_//* ) 

 
set field 2 a_/A Random(0.,1.,Nof(a_/A)) # set the 2nd field to random values 
color a_/* Field( a_/A  2 )              # color by it 

 
read pdb "/home/nerd/secret/hiv.ob" 
File( a_ ) 
 /home/nerd/secret/hiv.ob 

1. type character:
   • 'f' - regular file
   • 'd' - directory
   • 'l' - symbolic link
   • 'c' - character special file
   • 'p' - pipe
2. 'r' if you can read the file (or from the directory)
3. 'w' if you can write to this file (or directory)
4. 'x' if you can execute this file (or cd to this directory)
5. file size in bytes

 
 if File("/opt/icm/icm.rst")=="- - - - 0" print "No such file"   
 
 if Field(File("PDB.tab"),2)!= "w"  print "can not write"   
 
 if ( Indexx( File("/home/bob/icm/") , "d ? w x *" ) ) then 
    print "It is indeed a directory to which I can write" 
 endif 
             # Here the Indexx function matched the pattern.  
 
 if ( Integer(Field(File(s_name),5)) < 10 ) return error "File is too small"  

 
 (t.a=="word1" | t.b == "word1") & (t.a=="word2" | t.b == "word2"). 

 
 read database "ref.db"  # database of references  
 group table ref $s_out  # group created arrays into a table  
 show Find(ref,"energy profile") & ref.authors == "frishman" 

 
 user = Getenv("USER")      # extract user's name from the environment  
 if (user=="vogt") print "Hi, Gerhard" 

 
 show energy                # to calculate the gradient and its components  
 if (Gradient( ) > 10.) minimize 
 show Max(Gradient(a_//c*)  # show maximum "force" applied to the carbon atoms  

[ Grobselectbycolor ]

 
 GROB.relArrowSize = 0.1 
 g_arr = Grob("arrow",Box( ))  # return arrow between corners of displayed box 
 display g_arr red            # display the arrow  
 
 g_arr1 = Grob("ARROW100",{1. 1. 1.}) 
 display g_arr1 

 
 read csd "qfuran" 
 gcell = Grob("CELL",Cell( ) )      # solid cell  
 display a_//* gcell transparent   # fancy stuff  

 
  build string "se ala his trp" 
  g = Grob( "distance", a_/1/ca a_/2/ca ) 
  display g 
  GRAPHICS.displayLineLabels = no 
  display new 

 
 display a__crn.//ca,c,n 
 g = Grob("line",{0.,0.,0.,5.,5.,5.})   # a simple line (just as an example) 
 display g yellow 
 gCa = Grob("line",Rarray(Xyz(a_//ca))) # connect Cas with lines  
 display gCa pink                       # display the grobs  

 
 display a__crn.//ca,c,n 
               # make grob and translate to a_/5/ca  
               # Sum converts Matrix 1x3 into a vector 
 g=Grob("SPHERE",5.,2)+Sum(Xyz(a_/5/ca))   
               # mark it with dblLeftClick and 
               # play with Alt-X, Alt-Q and Alt-W 
 display g red   

 
buildpep "ADERD" # a peptide  
dsRebel a_  no no 
g=Grob(g_skin {0.9,1.,-0.1,0.1,-0.1,0.1} )  # red color 
display g_skin transparent 
display g 
show Res(Sphere( g, a_//* 1.5)) 

• "min" - stores the minimal atomic property for each selected residue
• "max" - stores the maximal atomic property for each selected residue
• "avg" - (syn. "mean") stores the mean of properties for each selected residue
• "rms" - stores the root-mean-square deviation of properties for each selected residue
• "sum" - stores the sum of properties for each selected residue
• "first" - stores the property of the first atom in selected residue
• "count" - stores the number of selected atoms in selected residue

 
read pdb "1crn" 
show Group( a_A//* "count"  ) # numbers of atoms in residues 
show Group( Mass( a_A//* ) , a_A//* "sum"  ) # residue masses 
show Group( Mass( a_A//* ) , a_A//* "rms"  ) # residue mass rmsd 

 
 plot display Histogram({ -2, -2, 3, 10, 3, 4, -2, 7, 5, 7, 5}) BAR 
 
 a=Random(0. 100. 10000)                   
 u=Histogram(a 50)                         
 s_legend={"Histogram at linear sampling curve" "Random value" "N"} 
 plot display regression BAR u s_legend 
 
 a=Random(0. 100. 10000)                   
 b=.04*(Count(1 50)*Count(1 50))           
 u=Histogram(a b)                          
 s_legend={"Histogram at square sampling curve" "Random value" "N"} 
 plot display BAR u s_legend 
 b=Sqrt(100.*Count(1 100))                 
 s_legend={"Histogram at square root sampling curve" "Random value" "N"} 
 plot display green BAR Histogram(a b) s_legend 

[ Iarrayinverse ]

 
 a=Iarray(5)                        # returns {0 0 0 0 0}  
 a=Iarray(5,3)                      # returns {3 3 3 3 3}  
 b=Iarray({2.1, -4.3, 3.6})         # returns {2, -4, 4}  
 c=Iarray({"2", "-4.3", "3.6"})     # returns {2, -5, 3}  

 
Iarray({1 2 3} reverse)  # returns {3 2 1} 

 
build string "se ala" 
ii = Iarray( a_//c* ) # returns {6,8,12} 
Select( a_  ii )      # returns three carbons 

 
 show stack 
 iconf>       1       2       3       4       5 
 ener>    -15.3   -15.1   -14.9   -14.8   -13.3 
 rmsd>     84.5    75.3     6.4    37.2   120.8 
 naft>        3       0       4       0       2 
 nvis>       10       9       8       1       4 
Integer(stack)  # returns { 10 9 8 1 4 } 

• a sequence, e.g. IcmSequence(1crn_m)
• a string, e.g. "ASDGFRE", or "SfGDA;WER" .
• or residue selection, rs_ , (e.g. a_2,3/* ).

• standard L amino-acids: upper case one-letter code (B,J,X,Z are illegal), e.g. ACD
• D-amino acids: lower case for a corresponding amino acid (e.g. AaA for ala Dala ala )
• new molecule: use semicolon or dot as a chain separator ( ; ) ( e.g. AAA;WWW )

 
IcmSequence( a_/* ) # C-terminal residue "cooh" will be added if oxt is found 
IcmSequence( a_/* "","" )	 # no terminal groups will be added 
IcmSequence( a_/* "","@coo-" )    # "coo-" will be added only if oxt is found 
IcmSequence( a_/* "nh3+","coo-" ) # "nh3+" and "coo-" will always be added  

 
 write IcmSequence(seq1) "seq1.se"      # create a sequence  
                                        #  file for build command  
 
 show IcmSequence("FAaSVMRES","nh3+","coo-")  # one peptide with Dala 
 
 show IcmSequence("FAAS.VMRES","nter","cooh") # two peptides  
 show IcmSequence("AA;MRES","nter","cooh") # two peptides  
 
 read pdb "2ins" 
 write IcmSequence(a_b,c/* ,"nter","@cooh") "b.se" # .se file for b  
                                                  # and c chains  

 
 buildpep "SDSRAARESW;KPLKPHYATV"  # two 10-res. peptides 

 
 show Index("asdf","sd")  # returns 2  
 show Index("asdf" "wer") # return 0  
 a=Sequence("AGCTTAGACCGCGGAATAAGCCTA") 
 show Index(a "AATAAA")   # polyadenylation signal 
 show Index(a "CT" last)  # returns 22 
 show Index(a "CT" 10)    # starts from position 10. returns 22  
 show Index(a "CT", -10)  # search only the last 10 positions  

 
EX = "xxxx xxxxxx xxxxxxxxxxxxxx xxxxxx  xxx " 
sp=0 
while(yes) 
  x=Index(EX "x" sp) 
  if(x==0) break 
  sp=Index(EX " " x) 
  print x sp-1 
endwhile 

 
   group table t {33 22 11} "A" {"a","b","c"} "B" 
   Index(t.A==22)  # returns 2 for 2nd row 
 #>I 
  2 
   t.B[ Index(t.A== 22 )[1] ]   # returns B according to A value 
 b 

 
 show Index({"Red Dog","Amstal","Jever"}, "Jever")      # returns 3  
 show Index({"Red Dog","Amstal","Jever"}, "Bitburger")  # returns 0  

 
 l_commands = no 
 read pdb "1crn" 
 read object "crn" 
 printf "The object a_crn. is the %d-nd, while ...\n", Index(object) 
 set object a_1. 
 printf "the object a_1crn. is the %d-st.\n", Index(object) 

• * any string including an empty string;
• ? any single character;
• [ string ] any of the enclosed characters;
• [! string ] any but the enclosed characters.
• ^ beginning of a string
• $ string end

 
 show Indexx("asdf","s[ed.]")       # returns 2  
 show Indexx("asdfff","ff$")       # returns 5 (not 4) 
 show Indexx("asdf" "w?r")          # return 0  

• "all" (or no string option) select insertions of all types
• "nter" select only N-terminal fragments
• "cter" select only C-terminal fragments
• "loop" select only the internal loops

 
 read pdb "1phc.m/"   # read the first molecule form this pdb-file 
 read pdb "2hpd.a/"   # do the same for the second molecule 
 make sequences a_*.  # you may also read the sequence and 
                      # the alignment from a file 
 aaa=Align( )         # on-line seq. alignment. 
                      # You may read the edited alignment  
                      # worm representation 
 assign sstructure a_*. "_"  
 display ribbon 
 link a_*. aaa        # establish connection between sequences and 3D obj. 
 superimpose a_1. a_2. aaa 
 display ribbon a_*. 
 color a_1. ribbon green 
 color ribbon Insertion(a_1.1 aaa) magenta 
 color ribbon Insertion(a_2.1 aaa) red 
 show aaa 

 
 show Integer(2.2), Integer(-3.1)      # 2 and -3  
 jj=Integer("256aaa")                  # jj will be equal to 256  

 
# Let us integrate sqrt(x)  
 
 x=Rarray( 1000 0. 10. ) 
 plot x Integral( Sqrt(x) 10./1000. ) grid {0.,10.,1.,5.,0.,25.,1.,5.} display 
  
# Let us integrate x*sin(x). Note that Sin expects the argument in degrees 
    
 x=Rarray( 1000 0. 4.*Pi )              
# 1000 points in the [0.,4*Pi] interval 
 plot x Integral( x*Sin(x*180./Pi) x[2]-x[1] ) \ 
    {0., 15., 1., 5., -15., 10., 1., 5. } grid display  
# x[2]-x[1] is just the increment  

    
 x=Rarray(100 ,-.9999, .9999 ) 
 x=x*x*x 
 plot display x Integral((3*x*x-1.) x) cross 

 
 if (Error | Interrupt) return 

 
 build 
 resLabelStyle = "Ala 5"  # other styles also available  
 aa = Label(a_/2:5)       # extract residue name and/or residue number info  
 show aa                  # show the created string array  

 
 len=Length("asdfg")     # len is equal to 5  
 
 a=Matrix(2,4)            # two rows, four columns  
 nCol=Length(a)           # nCol is 4  
 
 read profile "prof"      # read sequence profile  
 show Length(prof)        # number of residue positions in the profile  
 
 vlen=Length({1 1 1})     # returns 1.732051  

 
 X = Random(1., 10., 10)  
 Y = 2.*X + 3. + Random(-0.1, 0.1, 10)  
 lfit = LinearFit(X Y) 
 printf "Y = %.2f*X + %.2f\n",     lfit[1],lfit[2]  
 printf "s.d. = %.2f; r = %.3f\n", lfit[3],lfit[4]  
 show column X, Y, X*lfit[1]+lfit[2], R_out  

 
 M1=Transpose(X)*X   
 M2=Power(M1,-1)  
 W =(M2*Transpose(X))*Y  

 
read pdb "1avx" 
Link(a_1) 
 #>S string_array 
 TRYP_PIG 
rename a_1 Link(a_1)[1] 

 
 print Log(2.)          # prints 0.693147  
 print Log(10000, 10)    # decimal logarithm 
 print Log({1.,3.,9.}, Sqrt(3.)) # {0. 2. 4.}  

 
 read object "crn" 
 read map "crn"   
 display a_//ca,c,n m_crn 
 m1 = Map(m_crn, {0 0 0 22 38 38})  # half of the m_crn  
 m2 = Map(m_crn, {0 0 0 88 38 38})  # double of the m_crn  
 display m1 
 display m2 

 
 buildpep "ala his trp glu"  
 objmasses = Mass( a_*. ) 
 molmasses = Mass( a_* ) 
 resmasses = Mass( a_/* ) 
 
 masses=Mass( a_//!?vt* )   # array of masses of nonvirtual atoms  
 molweight = Mass( a_1  )   # mol.weight of the 1st molecule  
 molweight = Sum(Mass( a_1//* )) # another way to calculate 1st mol. weight 

 
Matrix(3,{1. 2. 3.}) 
 #>M 
 1. 0. 0. 
 0. 2. 0. 
 0. 0. 3. 

 
 Matrix({1. 2. 3. 4. 5. 6.},3) 
 #>M 
 1. 2. 3. 
 4. 5. 6. 
 Matrix({1. 2. 3. 4. 5. 6.},-3) 
 #>M 
 1. 4. 
 2. 5. 
 3. 6. 
 Matrix({1. 2. 3. 4. 5. 6.},4) 
 Error>  non-matching dimension [4] and vector size [6] 

 
 Matrix( Matrix(3) 0,0,1, 2)  # first two columns 

 
icm/def> m 
#>M m 
1. 0. 0. 
0. 1. 0. 
7. 0. 7. 
icm/def> Matrix( m right ) 
#>M 
1. 0. 0. 
0. 1. 0. 
0. 0. 7. 
icm/def> Matrix( m left ) 
#>M 
1. 0. 7. 
0. 1. 0. 
7. 0. 7. 

 
icm/def> Matrix(comp_matrix "CAR") 
#>M 
2.552272 0.110968 -0.488261 
0.110968 0.532648 -0.133162 
-0.488261 -0.133162 1.043102 

 
 mm=Matrix(2,4)          # create empty matrix with 2 rows and 4 columns 
 
 mm=Matrix(2,4,-5.)      # as above but all elements are set to -5. 
 
 show Matrix(3)          # a unit matrix [1:3,1:3] with diagonal 
                         # elements equal to 1. 
 a=Matrix({1. 3. 5. 6.}) # create one row matrix [1:1,1:4 ] 
 
 Matrix({1.,0.},{0.,1.}) # tensor product 
 #>M 
 0. 1. 
 0. 0. 

 
icm/def> Matrix(Random(0. 5. 20) Random(0. 5. 20) {0. 5. 0. 5. 1. 1.}) 
#>M 
1. 0. 2. 1. 0. 
1. 1. 1. 1. 2. 
0. 2. 1. 1. 0. 
0. 1. 0. 0. 1. 
1. 0. 0. 1. 2. 

 
  g_skin_1 
  g_skin_2 

 
 show Max({2. 4. 7. 4.})          # 7. will be shown  

 
 print Mean({1,2,3})          # returns 2.  
 
 show Mean(Xyz(a_2/2:8))      # shows {x y z} vector of geometric  
                              # center of the selected atoms  
 Mean({1. 2. 3.} {2. 3. 4.}) 
 #>R 
   1.5 
   2.5 
   3.5 

 
 show Min({2. 4. 7. 4.})      # 2. will be shown  
 show Min(2., 4., 7., 4.)     # 2. will be shown  

• %m specification for the rounded integral amount;
• %.m specification to add cents after dot. The default is "$%.m", i.e. Money(1222.33) returns $1,222.33.
• %M the same as %m but with dot instead of comma in the European style
• %.M the same as %.m dot and comma are inverted in the European style

 
 Money(1452.39)  # returns "$1,452.39"  
 Money(1452.39,"DM %m")  # returns "DM 1,452"  
 Money(1452.39,"%.M FF")  # inverts comma and dot "1.452,39 FF"  

 
 phi = Mod(phi)    # transform angle to [0., 360.] range 
 a   = Mod(17,10)  # returns 7 

 
 show Mol( Sphere(a_1//* 4.) )      
           # molecules within a 4 A vicinity of the first one  
           # Sphere function Sphere(as_atoms) selects atoms.  

• simple : removes non-alpha-numeric simbols from a string and replaces them by underscore.
• unique : checks if the name s_ exists in the ICM-shell. If the name does not exist, it is returned without changes, otherwise a number is appended to the name to guarantee its uniqueness.

 
  Name( " %^23 a 2,3 xreno-77-butadien" simple)     
 23_a_2_3_xreno_77_butadien 
 
  a=1 
  Name("a",unique) 
 a1 

 
 read alignment msf "azurins"        # load alignment  
 seqnames = Name(azurins)            # extract sequence names  
 
 show Name( Acc( a_/* ) ) # array of names of exposed residues  

 
 remarks = Namex( s_icmhome+"log3.ob") # get remarks directly without reading 
 group table t Rarray(Field(remarks,"vacuum",1,"\n")) "vacuum" 
 group table t append Rarray(Field(remarks,"hexadecane",1,"\n")) "hex" 
 show t 
 
 read object s_icmhome + "log3"  #  read multiple object file  
# extract numbers following the 'LogP' word in the object comments 
 logPs = Rarray(Namex(a_*.),"LogP",1,"  \n")  

 
 read index s_inxDir+"/SWISS.inx"  
 read sequence SWISS[2]  # read the 2nd sequence from Swissprot 
 show Namex( sequence[0] ) 

 
 show Next( a_/4 ) )      # show residues number 5  
 phi5=Torsion( a_/5/c ) 
                          # psi5 formally belongs to the next residue  
 psi5=Torsion( a_//n & Next(a_/5) ) 

 
 BS his   
 display 
 display a_/his/he2 ball red 
 display Next( a_/his/he2 bond ) ball magenta # show atom preceding he2 
 
 cd2_neigh = Next( a_//cd2 bond ) 
 for i=1,Nof(cd2_neigh) 
   nei = cd2_neigh[i] 
   print "  Distance between a_//cd2 and ",Sum(Name(nei)), " = ", Distance( a_//cd2 nei) 
 endfor 

 
 read mol s_icmhome+"ex_mol.mol" 
 Nof( a_*. atom ) 

 
 for i=1,Nof("def.cnf",conf)  # stack is NOT loaded 
   read conf i 
 endfor 

 
 if(Nof(String(dn1),"[!ACGT]" pattern) > 0.5*Length(dn1)) print " Warning> Bad DNA sequence" 

 
 nseq = Nof(sequences)         # number of sequences currently loaded  
 if(Nof(object)==0) return error "No objects loaded" 
 
 if ( Nof( sequence select ) == 2 ) a = Align( select ) 

 
 show Obj( a_*./dod )             # show objects containing heavy water  

 
 avO=Min(Occupancy(a_//ca))     # minimal occupancy of Ca-atoms  
 show Occupancy(a_//!h*)        # array of occupancy of heavy atoms  
 color a_//* Occupancy(a_//*)   # color previously displayed atoms  
                                # according to their occupancy 
 color ribbon a_/A Occupancy(a_/A) # color residues by mean occupancy 

 
 sPath=Path("/usr/mitnick/hacker.loot") # returns "/usr/mitnick/"  

 
read all Path(preference)+"icm.ini" # restore the settings 

 
read table mol "ex_mol.mol" name="t" 
sss  = String(t.mol[1]) # sss contains mol/sdf text 
t.mol[1] = Parray(sss mol)  # sss is parced and converted 

 
 read sequence s_icmhome + "zincFing" 
 group sequence aaa 
 align aaa 
 show Pattern("#~???A%  ?P")  # symbols from consensus string 
 show Pattern(aaa) 
 show Pattern(aaa exact) 

 
 print Pi/2. </tt> 

 
 read object "crn" 
            # prepare electrostatic boundary descriptions 
 make boundary  
            # potential from oe*, od* at cz of two args 
 show Potential(a_/arg/cz a_/glu,asp/o?* ) 
 print 0.5*Charge(a_//*)*Potential(a_//* a_//* )  
            # the total electrostatic energy which is 
            # actually calculated directly by show energy "el"  

 
 Power(2.,{1. 2. 3.}) # returns {2.,4.,8.} 

 
 size=Power(tot_volume,1./3.)        # cubic root  
 
 read matrix "LinearEquationsMatrix" # read matrix [1:n,1:n] 
 read rarray b                       # read the right-hand column [1:n] 
 x=Power(LinEquationsMat,-1) * b     # solve system of linear equations 
 
 a=Rot({0. 1. 0.}, 90.0)                 
       # create rotation matrix around Y axis by 90 degrees  
 if (Power(a,-1) != Transpose(a)) print "Wrong rotation matrix" 
                            # the inverse should be  
                            # equal to the transposed  
 rotate a_1 Power(a,3)      # a-matrix to the third is  
                            # three consecutive rotations  

 
                   # chance to find residues RGD at a given position 
 show Probability("RGD")        
                   # a more tricky pattern  
 show Probability("[!P]?[AG]")  

• no argument: alignment score
• identity: the number of identical residues * 100 and divided by the minimal sequence length
• similarity: alignment score without gap penalties *100. and divided by the minimal sequence length
• comp_matrix: alignment score without gap penalties (unnormalized)
• sort: an additional score for ranking

 
 read alignment s_icmhome+ "sx"  # 2 seq. alignment 
 read pdb  s_icmhome+"x" 
 p= -Log(Probability(sx)) 
 display ribbon a_ 
 color ribbon a_/A -Rarray(p,sx,cd59)  
 # Rarray projects the alignment array to the sequence 

 
 read pdb sequence "2mhb" 
 read pdb sequence "4mbn" 
 m=Probability(2mhb_a 4mbn_m 7 30) 
 print " pProbability: Min=" Min(Min(m)) "Max=" Max(Max(m)) 
 PLOT.rainbowStyle="white/rainbow/red" 
# show probability of the chance matching (comparable to the BLAST P-value) 
 plot area m display color={.2 0.001} transparent={0.2 1.} link grid 
# OR show -Log10(Probability)  
 plot area -Log(m,10) display color={0.7 3.} transparent={0. 0.7} link grid 

 
 show Putenv("aaa=bbb")   
      # change/add variable 'aaa' with value 'bbb' to environment 
 show Getenv("aaa")        
      # check if it has been successful 

 
 print Random(5)            # one of the following: 1 2 3 4 or 5  
 print Random(2,5)          # one of the following: 2 3 4 or 5  
 print Random(2.,5.)        # random real in [2.,5.]  
 randVec=Random(-1.,1.,3)   # random 3-vector with components in [-1. 1.]  
 randVec=Random(3,-1.,1.)   # the same as the previous command  
 randMat=Random(-1.,1.,3,3) # random 3x3 matrix with components in [-1. 1.]  
 randMat=Random(3,3,-1.,1.) # the same as the previous command   
 Random(0., 1., 10, "gauss" # normal distribution 
 a=Random(seq_crn)          # a contains a randomized crambin sequence 

[ rarraysequenceprojection | rarrayalignmentprojection | rarrayproperties | RarrayAlignment ]

1. all elements by rows (the default)
2. all elements by columns
3. the upper triangle,
4. the lower triangle,
5. the diagonal elements
6. the upper triangle without diagonal elements
7. the lower triangle without diagonal elements

 
 a=Rarray(54)                     # create 54-th dimensional vector of zeros 
 a=Rarray(3,-1.)                  # create vector {-1.,-1.,-1.} 
 a=Rarray(5,1.,3.)                # create vector {1., 1.5, 2., 2.5, 3.} 
 a=Rarray({1 2 3})                # create vector {1. 2. 3.} 
 a=Rarray({"1.5" "2" "-3.91"})    # create vector {1.5, 2., -3.91} 
# 
 M=Matrix(2) 
 M[2,2]=2 
 M[1,2]=3 
 Rarray( M ) 
 Rarray( M 1 ) 
 Rarray( M 2 ) 
 Rarray( M 3 ) 
 Rarray( M 4 ) 
 Rarray( M 5 ) 

[ Rarrayinverse ]

1. seq1 to ali : RA = Rarray( R1 seq1 ali r_gapValues )
2. ali to seq2 : R2 = Rarray( RA ali seq2 )

 
Rarray({1. 2. 3.} reverse)  # returns {3. 2. 1.} 

 
read alignment s_icmhome+"sh3" 
t = Table(sh3) 
group table t Count(Nof(t)) "n" append  # add a column with 1,2,3,.. 
show t  # t looks like this: 
 #>T t 
 #>-cons--------Fyn---------Spec--------Eps8-------n--- 
    " "         0           1           1          1 
    " "         0           2           2          2 
    " "         0           3           3          3 
    " "         0           4           4          4 
    .           1           5           5          5 
 ... 
t2forFyn =  t.Fyn == 2   # table row for position 2 in seq. Fyn 
t2forFyn.n               # corresponding alignment position  

 
 s= Sequence("TTCCPSIVARSNFNVCRLPGTPEAICATYTGCIIIPGATCPGDYAN") # crambin sequence 
# 26-dim. hydrophobicity vector for A,B,C,D,E,F,.. 
 h={1.8,0.,2.5,-3.5,-3.5,2.8,-.4,-3.2,4.5,0.,-3.9,3.8,1.9,-3.5,.0,-1.6,-3.5,-4.5,-.8,-.7,0.,4.2,-.9,0.,-1.3,0.} 
 hs=Rarray(s,h)              # h-array for each sequence position 
 hh = Smooth(Rarray(s,h), 5) # window average 

 
read alignment "sh3" 
show Rarray(sh3) 
# 
a=Rarray(sh3 simple)  # a number for each alignment position 
# to project a to a particular sequence, do the following 
b=Rarray(a,sh3,Spec)  # a number for each Spec residue 
String(Rarray(a, sh3, Spec ))//String(Spec)   # example  

 
 s = "5.3"         
 a = Real(s)      # a = 5.3  
 s = "5.3abc"     # will ignore 'abc'    
 a = Real(s)      # the same, a = 5.3  

 
 a[2,3]=Real(Value(v_/2/phi))  # a is a matrix, a[2,3] expects a real 

 
              # transform angle to the standard  
              # [-180., 180.] range. (Period=360 is implied)  
 phi=Remainder(phi)  
 
              # we assume that you have two objects 
              # with different conf. of the same molecule 
 phiPsiVec1 = Value(v_1.//phi,psi)    
 phiPsiVec2 = Value(v_2.//phi,psi)    
              # average angular  
              # deviation  
 angDev=Mean(Abs(Remainder(phiPsiVec1-phiPsiVec2)))   
              # cut and paste these examples into the ICM-shell 
 print Remainder(13,10 )  Mod(13,10 ) 
 print Remainder(17,10 )  Mod(17,10 ) 
 print Remainder(-13,10 ) Mod(-13,10 ) 
 print Remainder(-17,10 ) Mod(-17,10 ) 

 
 a=Replace(" 1crn "," ","")  # remove empty space 

 
 invertedSeq = String(0,1,"GTAAAGGGGTTTTCC")  # result: CCTTT.. 
 complSeq=Replace(invertedSeq,{"A","C","G","T"},{"T","G","C","A"}) 
# result: GGAAA... 

 
 cleanStr=Replace("XXTEXTYYTEXT",{"XX","YY"},"") 

 
 aa={"Terra" "Tera" "Teera" "Ttera"} 
 show column aa Replace(aa "er?" "ERR") Replace(aa "*[tT]" "Shm")  

 
 show Res( Sphere(a_1/1/* 4.) )   # show residues within 4 A  
                                  # vicinity from the firsts one  

 
 read pdb "1crn" 
 read alignment "bb" # a short sublalignment of 1crn_m 
 link a_*. bb        # link objects and sequences 
 show bb 
 #>ali bb 
 # Consensus   .C CP~.#A.^.# 
 1crn_m        TC-CPSIVARSNF------ 
 a             PCGCPDGIAIARIYPFAVG 
 #1crn_m       EE E__HHHHHHH 
   
 # 1crn_m a nID 4 Lmin 46 ID 8.7 % Score 4.75 Sim 15.54 Gap 2.40 nOverlap 12 pP 0.16 
 #MATGAP gonnet 2.4 0.15 
 
 Res(  bb 1  ) 
 - Num  Res. Type ----{SS Molecule}-- Object - sf - sfRatio 
     2  thr   Amino    T E  m            1crn     0.0  0.00  . 
     3  cys   Amino    C E  m            1crn     0.0  0.00  . 
     4  cys   Amino    C E  m            1crn     0.0  0.00  . 
     5  pro   Amino    P _  m            1crn     0.0  0.00  . 
     6  ser   Amino    S _  m            1crn     0.0  0.00  . 
     7  ile   Amino    I H  m            1crn     0.0  0.00  . 
     8  val   Amino    V H  m            1crn     0.0  0.00  . 
     9  ala   Amino    A H  m            1crn     0.0  0.00  . 
    10  arg   Amino    R H  m            1crn     0.0  0.00  . 
    11  ser   Amino    S H  m            1crn     0.0  0.00  . 
    12  asn   Amino    N H  m            1crn     0.0  0.00  . 
    13  phe   Amino    F H  m            1crn     0.0  0.00  . 

 
 res=Resolution(a_1crn.) </tt>  
 print "PDB structure 1crn: resolution = ", res, " A" 

 
 read factor "igd"   # read h,k,l,fo table from a file 
 read pdb "1igd"     # cell is defined there 
 defCell = Cell(a_)  # extract the cell parameters from the object 
 group table append igd Resolution(igd) "res" 
 show igd 

 
 read pdb "1mbn"                          # load myoglobin  
 read pdb "1pbx"                          # load alpha and beta 
                                          # subunits of hemoglobin  
 print Rmsd(a_1.1 a_2.1 align)            # myo- versus alpha subunit 
                                          # of hemo- all atoms 
 print Rmsd(a_1.1//ca a_2.1//ca align)    # myo- versus alpha subunit 
                                          # of hemo- Ca-atoms 
 
 print Rmsd(a_1./4,29/ca a_2.1/2,102/cb exact) # exact match  

 
         #  rotate molecule by 30 deg. around z-axis  
 rotate a_* Rot({0. 0. 1.},30.) 

 
# rotate by 30 deg. around {0.,1.,0.} axis through the center of mass 
 nice "1crn" 
 R_pivot = Mean(Xyz(a_//*)) 
 transform a_* Rot(R_pivot,{0. 1. 0.}, 30.) 

[ Sarrayinverse ]

• append : will merge residue ranges
• name : will return a string array of residueName residueNumber records
• residue : will return an array of selection strings a_obj.mol/residueNumber records

 
Sarray(a_/2,4:5 name) 
 #>S string_array 
 def.a1/ala2 
 def.a1/trp4 
 def.a1/glu5 
Sarray(a_/2,4:5 residue) 
 #>S string_array 
 def.a1/2 
 def.a1/4 
 def.a1/5 
Field(Sarray(a_/2,4:5 name),2,"/") # extract residues 
#>S string_array 
 ala2 
 trp4 
 glu5 

• String( rs_ ) which returns one string;
• Label( rs_ | as_ ) which will format the output string according to the resLabelStyle or atomLabelStyle preference.
• l_showResCodeInSelection

• 'B': -200 < phi < -80 , 140 < psi < 200
• 'A': -101 < phi < -24 , -81 < psi < 4
• 'g': -169 < phi < -15 , -64 < psi < 54
• 'd': -211 < phi < -5 , 8 < psi < 136
• 'L': 24 < phi < 101 , -4 < psi < 81
• '_': the rest

• 'M': -120 <= xi1 < 0
• 'P': 0 <= xi1 < 120
• 'T': 120 <= xi1 < 240

 
 show Sarray(stack,v_/2:10/x*)      # coding of side-chain conformations 
 show Sarray(stack,v_//phi,psi)     # backbone conformation character coding 
 show Sarray(stack,v_/2:10/phi,PSI) # character coding of a chain fragment 

 
 ss=Sarray(5)       # create empty sarray of 5 elements  
 ss[2]="thoughts"   # assign string to the second element of the sarray  
 
 sa=Sarray("the first element") 
 
 show Sarray(Count(1 100)) # string array of numbers from 1 to 100 

 
Sarray({"one","two"} reverse)  # returns {"two","one"} 

 
a={"123","12345"} 
Sarray(a,2,3) 
 {"23","23"} 
Sarray(a,5,2) 
 {"32","5432"} 

[ Scoreali_ ]

 
 show Score({1. 2. 5. 3.} {3. 1.5 1.5 5.})  # 0. perfectly overlapping arrays 
 show Score({2. 5. 3.} {1. 1.5 0.5})        # 1. no overlap R_1 > R_2 
 show Score({1. 1.5 0.5} {2. 5. 3.})        # -1. no overlap R_2 > R_1 
 show 1.-Abs(Score({1. 3. 2.5} {2. 5. 3.})) # relative overlap between R1 and R2 

• starting from distances (e.g. RMSD) from the correct answer Di (Di >=0.)
• calculating their well-behaved and inverted version, exp( -w*D ) [0., 1.]
• calculating the normalized Boltzmann average if the previous similarity measure
• taking -Log of the previous average

 
 read sequence "seqs" 
 a = Score( Azur_Alcde Azur_Alcfa ) # it is around 90. 

• a column is extracted from a linked N-sequence alignment ( see the linkcommand )
• Si = Sum( Cij )/N where j=1,..N and Cij is the residue comparison value
• simple mode: Cij = 1. for two identical residues and 0. otherwise
• the default mode: Cij is taken from a normalized comp_matrix . Its elements are calculated as Cij_norm = Cij/Sqrt(Cii*Cjj) .

 
  read alignment s_icmhome+"sh3.ali" 
  read pdb "1fyn" 
  make sequence a_a 
  group sequence sh3 
  align sh3 
  display ribbon 
  color ribbon a_a/A Score( a_a/A simple ) 
  show surface area 
  show Mean( Score( Acc(a_a/*) ) )    # conservation score for the surface 
  show Mean( Score( a_a & !Acc(a_/*)))# conservation score for the buried 

• no second argument : the straight Needleman and Wunsch score: aligned residues score according to the residue comparison table, gaps according to the gapOpen and gapExtension parameters.
• identity the number of identical residues in the alignment divided by the smallest sequence length and multiplied by 100 %.
• comp_matrix the alignment score without the gap component. It contains only the total score of the aligned residues calculated from the residue comparison table and does not include penalty term.
• similarity the alignment score without the gap component multiplied by 100. and divided by the smallest sequence length.
• sort a score occasionally used for ranking/sorting the alignments in fold recognition. Currently it is equal to the comp_matrix_score - 1.3*totalGapPenalty

 
 Score( 150, 1./55000.,identity)  

 
  display skin Select(residue) 

• x,y,z atomic coordinates ("x","y","z")
• bfactor ("bfactor","b","B")
• occupancy ("o")
• charge ("charge","c","q")
• accessible surface area ("area","a")
• user-field ("u") which can be set with set field and extracted with the Field function.
• residue user-fields: "u","v","w" for the 1st, 2nd and 3rd field, respectively.

 
 build string "se glu arg" 
 show Select(a_//* "charge < 0.")|Select(a_//c* "x> -2.4") 
 show Select(a_//c* , "x>", -2.4) 
 
 show Select(a_/* , "w>3.") # 3rd res. user field greater than 3. 

 
 buildpep "ASD" 
 aa = a_/2/c*          # selection in the current obj a_ 
 copy a_ "b"           # a copy of the source object 
 bb = Select(aa,a_b.)  # selection aa moved to a_b. 

 
 seqA = Sequence( a_1./15:89 )     # create sequence object  
                                   # with fragment 15:89  
 
 show Align(seq1, Sequence("HFGD--KLS AREWDDIPYQ")  
                                   # non-characters will be squeezed out 
 a=Sequence("ACTGGGA", nucleotide)  
 Type(a , 2)  # returns the type-string : 
 nucleotide 

 
     nucleotide      |complement 
_____________________|__________ 
                     |   
 A = Adenosine       | T (replace by U for RNA) 
 C = Cytidine        | G 
 G = Guanosine       | C 
 T = Thymidine       | A 
 U = Uridine         | A 
 R = puRine    (G A) | Y 
 Y = pYrimidine(T C) | R 
 K = Keto      (G T) | M 
 M = aMino     (A C) | K 
 S = Strong    (G C) | S 
 W = Weak      (A T) | W 
 B = !A      (G T C) | V 
 D = !C      (G A T) | H 
 H = !G      (A C T) | D 
 V = !T      (G C A) | B 
 N = aNy             | N 

 
 Sign(-23)  
 -1 
 Sign(-23.3)  
 -1. 
 Sign({-23,13})  
 {-1,1} 
 Sign({-23.0,13.1})  
 {-1.,1.} 

 
 print Sin(90.)                  # equal to 1  
 print Sin(90)                   # the same  
 
 print Sin({-90., 0., 90.})      # returns {-1., 0., 1.}  

 
 print Sinh(1.)                  # equal to 1.175201  
 print Sinh(1)                   # the same  
 print Sinh({-1., 0., 1.})       # returns {-1.175201, 0., 1.175201}  

 
 nice "1est"  # contains some sites  
 delete site a_1 Site("CONFLICT",a_1)  

[ smoothr | Smoothrs | smoothali | smoothmap ]

 
 {-1.,1.}/Xstep             # for the first derivative 
 {1.,-2.,1.}/(Xstep*Xstep)  # for the second derivative 
 {-1.,3.,-3.,1.}/(Xstep*Xstep*Xstep)  # for the third derivative 
 #  ... etc. 

 
 gauss=Exp( -Power(Rarray(31,-1.,1.) , 2) )   # N(0.,1.) distribution on a grid 
 x = Rarray(361,-180.,180.)                   # x-array grows from 0. to 180.  
 a = Sin(x) + Random(-0.1,0.1,361)            # noisy sine  
 
 b = Smooth(a,gauss)          # gauss averaging   
                              # see how noise and smooth signals look  
 plot x//x a//b display {-180.,180.,30.,10.} 
                              # take the first derivative of Sin(x)  
 c = Smooth(Sin(x),{-1., 1.}) * 180.0 / Pi  
                              # plot the derivative  
 plot x c display {"X","d(Sin(X))/dX","Derivative"} 

 
 nice "1tet"   # it is a macro displaying ribbon++ 
 R = Bfactor(a_/A ) # an array we will be 3D-averaging 
 color ribbon a_/A Smooth(a_/A R 1.)//5.//30.  # averaging with 1A radius 
 color ribbon a_/A Smooth(a_/A R 5.)//5.//30.  # with 5A radius 
 color ribbon a_/A Smooth(a_/A R 10.)//5.//30. # with 10A radius 
  # 5.//30. are appended for color scaling from 5. (blue) to 30.(red) 
  # rather than automated rescaling to the current range 
 
 set field a_/A Smooth(a_/A R 5.) 
 show Select( a_/A "u>30." ) # select residues with 1st field > 30. 

 
 m_gc = Smooth(Smooth(m_gc "expand"), "expand" )  

 
 p = Parray( "O=O.NN", smiles ) 
 Nof(p) # 1 element 
 1 
 Split( p, mol ) 
 O=O 
 NN 
 

• it supports one connection at a time
• requires a running MySql server (see www.mysql.com) with a database
• the record columns from the database are converted into the ICM sarrays, rarrays and iarrays with one exception. The mysql BIGINT values can not be converted into iarray and are stored as a string array (sarray).

• s_host (default "localhost") - the host name
• s_loginName (default "root")
• s_password - the database password
• s_dbName - the database name

 
 if !Sql( connect "localhost","john","secret","swiss") print "Error"  
 id=24 
 T =Sql( "SELECT * FROM swissprothits  WHERE featureid="+id ) 
 sort T.featureid 
 web T 
# Another example 
 tusers = Sql("select * from user where User=" + s_usrName ) 
 Sql(off) 

 
 show Sqrt(4.)                # 2.  
 show Sqrt({4. 6.25})         # {2. 2.5}  

• a group of atoms ( as_)
• any vertex point of a grob ( g_)
• a point in space ( R_xyz)
• a group of points in space ( M_xyz) (e.g. see the Xyz function)

 
 show Sphere( a_subA/14:15/ca,c,n,o , 5.2) 
 Res(Sphere( a_1.2 a_2.)) # residues of a_2. around ligand a_1.2 

 
 lines=Split("a 1 \n 2","\n")  # returns 2-array of {"a 1" " 2 "}   
 flds =Split("a b c")          # returns 3-array of {"a" "b" "c"}   
 flds =Split("a b:::c",":")    # returns 4-array of {"a b","","","c"}  
 resi =Split("ACDFTYRWAS","")  # splits into individual characters  
                               # {"A","C","D","F",...}  

 
 superimpose a_1.1 a_2.1            # two similar objects, each  
                                    # containing two molecules  
 print Srmsd(a_1.2//ca a_2.2//ca)   # compare how second molecule  
                                    # deviates if first superimposed  

[ strinv | stringali | ali_seq_project | seq_ali_project | Stringas | chemformula ]

 
 file=s_tempDir//String(Energy("ener"))   # tricky file name  
 show Index(String(seq),"AGST")              # use Index to find seq. pattern  
 tenX = String("X",10)                       # generate "XXXXXXXXXX"  
 
 read matrix  
 show String(def," ..:*#") 

 
 String(1,3,"12345")  # returns substring "123" 
 String(4,2,"12345")  # returns substring "432" 
 String(1,0,"12345")  # returns "12345" 
 String(0,1,"12345")  # returns INVERTED string "54321" 
 String(-1,1,"12345") # returns "4321" 

 
 read alignment s_icmhome+"sh3" 
 offs=Mod(Indexx(String(sh3),"--NSGDG"),Length(sh3)+1) 
         # extract alignment into a string, (+1 to account for '\n') 
 iSeq = 1 + Indexx(String(sh3),"--NSGDG")/(Length(sh3)+1) 
         # identify which sequence contains the pattern 

 
 read alignment s_icmhome+"sh3" 
 show String(sh3 tree) 

 
 read alignment s_icmhome+"sh3"  # 3 seq.  
 cc = Consensus(sh3) 
 show String(Spec)//String(cc,sh3,Spec) 

 
 read alignment s_icmhome+"sh3"  # 3 seq.  
 ssFyn = Sstructure(Fyn) 
 set sstructure Spec String(String(ssFyn,Fyn,sh3,"_"),sh3,Spec) 
 show Spec 

 
 nice "1crn" 
 l_showResCodeInSelection = no 
 nei = String( Res(Sphere( a_/leu a_/!leu , 4.)) )   
 show nei 
 a_1crn.m/14:17,19:20 
 display xstick $nei 

 
 read pdb "2ins" 
 for i=1, Nof( a_//c* ) 
    print String( a_//c* i ) 
 endfor 

 
 build string "se ala"     # alanine  
 show String(a_//!h* all ) # returns no hydrogen chemical formula: C3NO   
 show String(a_//* all )   # returns chemical formula: C3H5NO   
 show String(a_//* sln )   # returns SLN notation: NHCH(CH3)C=O  
 show String(a_//* smiles) # returns SMILES string: [NH][CH]([CH3])C=O  

 
 show Sstructure("HHHHHHH_____EEEEE",compress)  # returns string "HE"  
#    
 read object "crn" 
 show Sstructure( a_/A , compress)  # returns string "EHHEB"  

 
 show Index(Sstructure(a_1crn.,"HHHHHH"))  # first occurrence of  
                                           # helix in crambin  
 
 read sequence "sh3"    # load 3 sequences (the full name is s_icmhome+"sh3")  
 show Sstructure(Spec)  # secondary structure prediction for one of them 
 show Sstructure("AAAAAAAAAAAAA")  # sec. structure prediction for polyAla 
 
 read sequence "fasta_results.seq" 
 group sequences a unique 0.05     # remove redundant sequences 
 show Sstructure(my_seq_name)      # the actual prediction, be patient 
 plot number M_out display         # plot 3 propensities and reliability 

 
 show Sum({4 1 3})                     # 8  
 show Sum(Mass(a_1//*))                # mass of the first molecule  
 show Sum({"bla" "blu" "bli"})         # "bla blu bli" string  
 show Sum({"bla" "blu" "bli"},"\t")    # separate words by TAB  
 show Sum({"bla" "blu" "bli"},"\n")    # create a multiple line string  

 
 iGroup = Symgroup("P212121")  # find the group number=19 
 print "N_mol. in the cell =", Nof(Symgroup(iGroup))/12 

 
 read string       # read from stdin in to the ICM s_out string 
 a=Table(s_out)    # create table a with arrays a.name and a.value 
 show a            # show the table 
 for i=1,Nof(a)    # just a loop accessing the array elements 
   print a.name[i] a.value[i] 
 endfor 

 
#>T pos 
#>-cons----seq1-----seq2------- 
   " "       0        1           
   " "       0        2           
   C         1        3           
   " "       2        0           
   ~         3        4           
   C         4        5           
   " "       0        6           
# for the following alignment: 
# Consensus    C ~C   
seq1         --CYQC-  
seq2         LQC-NCP  

 
 read alignment "sh3"    
 t = Table(sh3 number)  # arrays t.1 t.2 t.3  
 t = Table(sh3)         # arrays t.cons t.Fyn t.Spec t.Eps8  
# 
 cc = t.cons ~ "[A-Z]"  # all the conserved positions  
 show cc                # show aa numbers at all conserved positions 
 show t.Fyn>=10 & t.Fyn<=20 # numbers of other sequences in this range 

 
% icm 
 buildpep "ala his trp" 
 montecarlo 
 show stack 
 iconf>       1       2       3       4       5       6       7 
 ener>    -15.1   -14.6   -14.6   -14.2   -13.9   -11.4    -1.7 
 rmsd>      0.3    39.2    48.0    44.1    27.4    56.6    39.3 
 naft>        1       0       0       1       1       1       0 
 nvis>        4       1       1       4       4       4       1 
t= Table(stack) 
show t 
 #>T t 
 #>-i--ener--------rmsd--------naft--------nvis------- 
    1  -15.126552  0.295555    1           4 
    2  -14.639667  39.197378   0           1 
    3  -14.572973  47.996203   0           1 
    4  -14.220515  44.058755   1           4 
    5  -13.879041  27.435388   1           4 
    6  -11.438268  56.636246   1           4 
    7  -1.654792   39.265912   0           1 
t1= Table(stack v_//phi,psi)  # show also five phi-psi angles 
 #>T 
 #>-ener----rmsd--naft-nvis------v1------v2------v3------v4-------v5----- 
    1 -15.12   0.29   1   4     -79.10  155.59  -75.30  146.99  -141.13 
    2 -14.63  39.19   0   1    -157.22  163.56  -78.25  139.51  -137.30 
    3 -14.57  47.99   0   1    -157.26  166.87  -85.08   92.55   -84.74 
    4 -14.22  44.05   1   4     -67.65   80.43  -76.67  103.05   -81.85 
    5 -13.87  27.43   1   4     -82.72  155.86  -85.02   93.11   -81.46 
    6 -11.43  56.63   1   4     -78.28  152.80 -154.79   66.26   -77.61 
    7  -1.65  39.26   0   1     -78.17  169.41 -133.89   96.39   -96.03 

 
 show Tan(45.)               # 1.  
 show Tan(45)                # the same  
 
 show Tan({-30., 0. 60.})    # returns {-0.57735, 0., 1.732051}  

 
 show Tanh(1)                 # returns 0.761594  
 show Tanh({-2., 0., 2.})     # returns -0.964028, 0., 0.964028  

 
 Transpose( xyz ) * xyz / Nof(xyz) 

 
 build string "se ala ala ala ala"  # let is define the ellipsoid 
 display virtual 
 a = Tensor(Xyz(a_//!h*))  # Xyz returns matrix K by 3 
 b=Eigen(a)                # 3x3 matrix of 3 eigenvectors 
 b[?,1] = b[?,1] / Length(  b[?,1] ) # normalize V1 in place 
 b[?,2] = b[?,2] / Length(  b[?,2] ) # normalize V2 
 b[?,3] = Vector( b[?,1], b[?,2] )   # V3 is a vector product V1 x V2 
 rotate a_ Transpose( b )  # b is the rotation matrix now 
  # Transpose(b) is the inverse rotation 
 set view  # set default X Y Z view 

 
 if (Time( ) > 3660.) print "Tired after " Time(string) " of work?" 

 
 show Tolower("HUMILIATION") 
 
 read sarray "text.tx"  #create sarray 'text' (file extension is ignored) 
 text1 = Tolower(text) 

 
 d=Torsion( a_/4/c )                   # d equals C-Ca-N-C angle  
 
 print Torsion(a_/4/ca a_/5/ca a_/6/ca a_/7/ca)  # virtual Ca-Ca-Ca-Ca  
                                                 # torsion angle  

 
 show Toupper("promotion") 
 
 read sarray "text.tx" 
 text1 = Toupper(text)