Copyright © 2020, Molsoft LLC Jun 5 2024
|
[ Plotting vw faq | REBEL faq | PK shift | Binding energy | Ensemble average | Helicity | Stacks merge ]
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.}
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_//* )
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.))
There are many different approaches to the evaluation of binding energy. One of the reasonable approximations has the following features:
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
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
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
|
Copyright© 1989-2024, Molsoft,LLC - All Rights Reserved. Copyright© 1989-2024, 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. |