User:Tlinnet

From PyMOLWiki
Revision as of 17:13, 19 November 2015 by Tlinnet (talk | contribs)
Jump to navigation Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

Author of following PyMOL scripts

PROPKA for PyMOL 2015/11/19 42k
Displacement map for conformational changes 2015/11/19 43k
F?orster distance calculator 2015/11/19 22k
Rotation matrix calculator 015/11/19 17k
Color by displacement 2015/11/19 13k
Surface Cysteine pKa predictor 2015/11/19 12k

149 k

Tutorial

Biochemistry student intro 2015/11/19 60k

Test

  1. from __future__ import print_function

from pymol import cmd import os import sys import math from time import localtime, strftime

  1. Thx for inspiration from Simple scriptin PymMOl http://www.pymolwiki.org/index.php/Simple_Scripting
  2. Made by Ma.Sc student. Troels Linnet, 2011-08. troels.linnet@bbz.uni-leipzig.de
  3. Based solely on the work by:
  4. Maik H. Jacob, Dan Amir, Vladimir Ratner, Eugene Gussakowsky, and Elisha Haas
  5. Predicting Reactivities of Protein Surface Cysteines as Part of a Strategy for Selective Multiple Labeling. (Biochemistry 2005, 44, 13664-13672)
  1. Example of pymol script: Directory "predict_reactivity" has script file cyspka.py and cysteine residue pdb file: cys.pdb
  2. import cyspka
  3. fetch 4AKE, async=0
  4. create 4AKE-A, /4AKE//A and not resn HOH
  5. delete 4AKE
  6. hide everything
  7. show cartoon, 4AKE-A
  8. cyspka 4AKE-A, A, 18


def cyspka(molecule, chain, residue, SeeProgress='yes', pH=7.2, MoveSGatom='no', SGatom=str((0, 0, 0))): 

   # If SeeProgress='yes', computation time will take 10-20% extra, but nice to follow.
   cmd.refresh()
   RotationRange = 360
   RotationDegree = 1
   # For error checking, the energies can be printed out
   printMC = 'no'
   printSC = 'no'
   # Parameters
   DieElecSpheDist = 7.0
   DieElecWaterDist = 1.4
   DieElecWater = 78.5
   DieElecCore = 4.0
   BornPenaltyB = 1.0
   AvogadroR = 8.31446216
   Temp = 298
   DeltapKMCSC = 0
   pK1 = 9.25
   pK2 = 8.0
   NotPopuDist = 2.4
   PopEnergyPenalty = 10000000
   # Side chain discrete charges
   DieElecSC = 40.0
   SCchargeASP = -1
   SCchargeGLU = -1
   SCchargeOXT = -1
   SCchargeARG = +1
   SCchargeHIS = +1
   SCchargeLYS = +1
   SCchargeMET1 = +1
   # Main chain partial charges
   NrMainchainNeighBours = 5
   DieElecMC = 22.0
   MCchargeC = +0.55
   MCchargeO = -0.55
   MCchargeN = -0.35
   MCchargeH = +0.35
   MCchargeProCA = +0.1
   MCchargeProCD = +0.1
   MCchargeProN = -0.2

   # Loading an Cys residue, give it a logic name, and aligning it. The oxygen atom can not be aligned in many cases, and are skipped.
   # We use only this molecule, to find the initial position of the SG atom, and to rotate the SG atom around the CA-CB bond. The molecule atom positions are not used for electric potential calculatons.
   Cysmolecule = str(molecule) + str(residue) + "Cys"
   cmd.fragment("cys")
   cmd.set_name('cys', Cysmolecule)
   # We use pair_fir, since align and super gets unstable with so few atoms
   pairfitCys(Cysmolecule, molecule, chain, residue)
   # Give nice representations quickly
   cmd.show("sticks", Cysmolecule)
   cmd.select(str(molecule) + str(residue) + "Res", "/" + molecule + "//" + chain + "/" + residue)
   print("/" + molecule + "//" + chain + "/" + residue)
   cmd.show("sticks", str(molecule) + str(residue) + "Res")
   cmd.disable(str(molecule) + str(residue) + "Res")
   # Find out what is the residuename we are investigating for
   Respdbstr = cmd.get_pdbstr(str(molecule) + str(residue) + "Res")
   Ressplit = Respdbstr.split()
   residueName = Ressplit[3]

   print("")
   print("# Hello, PyMOLers. It should take around 1 minute per residue.")
   print("# molecule: %s , chain: %s, residue: %s %s, pH: %s " % (molecule, chain, residueName, residue, pH))

   # Determine the range of neighbour residues possible.
   Maxresidues = cmd.count_atoms("/" + molecule + "//" + chain + " and name CA")
   for i in range(NrMainchainNeighBours + 1):
       if int(residue) - i >= 1:
           Minresidue = int(residue) - i
       else:
           break
   for i in range(NrMainchainNeighBours + 1):
       if int(residue) + i <= Maxresidues:
           Maxresidue = int(residue) + i
       else:
           break

   # Get the position and the vector for the CA->CB bond.
   dihedN = "/" + Cysmolecule + "//" + "/" + "/N"
   dihedCA = "/" + Cysmolecule + "//" + "/" + "/CA"
   dihedCB = "/" + Cysmolecule + "//" + "/" + "/CB"
   dihedSG = "/" + Cysmolecule + "//" + "/" + "/SG"
   dihedralPosCA = cmd.get_atom_coords(dihedCA)
   dihedralPosSG = cmd.get_atom_coords(dihedSG)
   dihedralVector = AtomVector(dihedCA, dihedCB)

   # To compare with article, we can move the SGatom to a starting position. The rotation is still determined around the CA-CB bond.
   if MoveSGatom == 'yes':
       SGatom = [float(SGatom[1:-1].split(",")[0]), float(SGatom[1:-1].split(",")[1]), float(SGatom[1:-1].split(",")[2])]
       Translate = [(SGatom[0] - dihedralPosSG[0]), (SGatom[1] - dihedralPosSG[1]), (SGatom[2] - dihedralPosSG[2])]
       cmd.translate(Translate, dihedSG, state=0, camera=0)
       dihedralPosSG = cmd.get_atom_coords(dihedSG)
   # Create a pymol molecule, that in the end will hold and show all SG atoms. Gives the representation of the rotameric states.
   SGName = str(molecule) + str(residue) + "SG"
   cmd.create(SGName, "None")
   # Create a pymol molecule, that in the end will hold and show all Amide protons. Gives a nice representation, and easy to delete.
   AmideName = str(molecule) + str(residue) + "NH"
   cmd.create(AmideName, "None")
   # Check if there are any nearby SG atoms, which could make a SG-SG dimer formation. The
   breakDimer = "no"
   breakDimer = CheckDimer(dihedSG, molecule, chain, residue)
   # Create a list for appending the calculated energies.
   ListofEnergies = []
   ListofRotamerDiscarded = []
   # print "Angle before rotation", cmd.get_dihedral(dihedN,dihedCA,dihedCB,dihedSG)
   # Enter into the loop of rotameric states
   for i in range(int(math.floor(RotationRange / RotationDegree))):
       Angle = i * RotationDegree
       # Create pymol molecule/SG atom for which we will calculate for.
       SGNameAngle = str(residue) + "SG" + str(Angle)
       cmd.create(SGNameAngle, dihedSG)
       # Calculate new coordinates for rotation around CA->CB bond. Then translate the created SG atom.
       SGNewPos = fRotateAroundLine(dihedralPosSG, dihedralPosCA, dihedralVector, Angle)
       Translate = [(SGNewPos[0] - dihedralPosSG[0]), (SGNewPos[1] - dihedralPosSG[1]), (SGNewPos[2] - dihedralPosSG[2])]
       cmd.translate(Translate, SGNameAngle, state=0, camera=0)
       # If one wants to "see it happen" while its making the states. But it will take extra computation time.
       if SeeProgress == 'yes':
           cmd.refresh()
       # Calculate number of neighbours within 2.4 Angstrom. Amide hydrogens are not considered, and are actually not build yet.
       nameselect = "(((/" + molecule + "//" + chain + " and not /" + molecule + "//" + chain + "/" + residue + ") or /" + molecule + "//" + chain + "/" + residue + "/N+CA+C+O)  within " + str(NotPopuDist) + " of /" + SGNameAngle + "//" + "/" + "/SG) and not resn HOH"
       # print nameselect
       cmd.select("NotPop", nameselect)
       NotPopNr = cmd.count_atoms("NotPop")
       # print Angle, NotPopNr, cmd.get_dihedral(dihedN,dihedCA,dihedCB,SGNameAngle)
       # If no neighbours, then proceed calculating
       if NotPopNr == 0:
           SumAllWMC = 0.0
           # Now calculate the electric potential due to the side chains.
           SumWSC = fSumWSC(molecule, SGNameAngle, chain, residue, DieElecSC, SCchargeASP, SCchargeGLU, SCchargeOXT, SCchargeARG, SCchargeHIS, SCchargeLYS, SCchargeMET1, printSC)
           # Now we calculate for the flanking 5 peptide groups on each side of the Cysteine CA atom.
           # For the first residue, only calculate for the tailing C,O atom in the peptide bond. No test for Proline.
           SumWMCFirst = fSumWMCFirst(molecule, SGNameAngle, chain, residue, Minresidue, DieElecMC, MCchargeC, MCchargeO, printMC)
           # For the residue itself, we dont test for PRO, since it should be a Cysteine.
           SumWMCresidue = fSumWMCresidue(molecule, SGNameAngle, chain, residue, int(residue), DieElecMC, MCchargeC, MCchargeO, MCchargeN, MCchargeH, AmideName, printMC)
           # For the last residue, we test for Proline. We only calculate for the N,H atom, or if Proline, N,CA and CD atom.
           SumWMCLast = fSumWMCLast(molecule, SGNameAngle, chain, residue, Maxresidue, DieElecMC, MCchargeN, MCchargeH, MCchargeProCA, MCchargeProCD, MCchargeProN, AmideName, printMC)
           # Then loop over rest of the residues in the chain.
           for j in (list(range(Minresidue + 1, int(residue))) + list(range(int(residue) + 1, Maxresidue))):
               MCNeighbour = j
               # print "Looking at neighbour", j
               SumWMC = fSumWMC(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeC, MCchargeO, MCchargeN, MCchargeH, MCchargeProCA, MCchargeProCD, MCchargeProN, AmideName, printMC)
               SumAllWMC = SumAllWMC + SumWMC
               # print "Rotation: %s Neighbour: %s " % (Angle, j)
           # Since the SG atom is negative, we multiply with -1.
           SumMCSC = -1 * (SumWSC + SumWMCFirst + SumWMCresidue + SumWMCLast + SumAllWMC)
           # Makes the neighbour count. Everything in 'molecule" within 7 ang of aligned SG atom. Not counting 'residue'. Adding 5 for 'residue' N,CA,C,O,CB
           ListNeighbourCount = fNeighbourCount(molecule, SGNameAngle, chain, residue, DieElecSpheDist)
           # Calculate the weighted electric potential and alter the b factor for coloring. Then add the rotated SG into bucket of SG atoms.
           SG_MCSC_Weight = fBoltzSingleState(SumMCSC, AvogadroR, Temp) * SumMCSC
           cmd.alter(SGNameAngle, 'b="%s"' % SG_MCSC_Weight)
           cmd.alter(SGNameAngle, 'name="S%s"' % Angle)
           cmd.create(SGName, SGName + " + " + SGNameAngle)
           # Then save the calculated values
           ListofEnergies.append([Angle, SumMCSC, ListNeighbourCount, NotPopNr, SG_MCSC_Weight, cmd.get_atom_coords(SGNameAngle)])
           cmd.delete(SGNameAngle)
       else:
           SumMCSCPenalty = PopEnergyPenalty
           ListNeighbourCount = fNeighbourCount(molecule, SGNameAngle, chain, residue, DieElecSpheDist)
           ListofRotamerDiscarded.append([Angle, SumMCSCPenalty, ListNeighbourCount, NotPopNr, 0, cmd.get_atom_coords(SGNameAngle)])
           cmd.delete(SGNameAngle)
   # Now show all the SG atoms as the available rotameric states.
   cmd.show("nb_spheres", SGName)
   cmd.delete("NotPop")
   cmd.spectrum("b", selection=SGName)
   AvailRotStates = len(ListofEnergies)
   # print "Available Rotational States: ", AvailRotStates

   # Do the calculations according to eq 5.
   # Find the partition function
   BoltzPartition = 0.0
   for i in range(len(ListofEnergies)):
       Boltz = fBoltzSingleState(ListofEnergies[i][1], AvogadroR, Temp)
       BoltzPartition = BoltzPartition + Boltz
   # Find the summed function
   BoltzSumNi = 0.0
   for i in range(len(ListofEnergies)):
       BoltzNi = fBoltzSingleState(ListofEnergies[i][1], AvogadroR, Temp) * ListofEnergies[i][1]
       BoltzSumNi = BoltzSumNi + BoltzNi

   # Check if there was any possible rotamers

   nostates = "no"
   if len(ListofEnergies) == 0:
       print("####################################################")
       print("########### WARNING: No states available ###########")
       print("########### Did you mutate a Glycine?    ###########")
       print("####################################################")
       BoltzSumNi = 0
       BoltzPartition = 0
       BoltzMCSC = 0
       DeltapKMCSC = 99
       NeighbourCount = 0
       nostates = "yes"
   else:
       # Final calculation
       BoltzMCSC = (BoltzSumNi) / (BoltzPartition)
       DeltapKMCSC = fDeltapK(BoltzMCSC, AvogadroR, Temp)

   # Find average number of neighbours
   NCSum = 0.0
   NCWeightedSum = 0.0
   for i in range(len(ListofEnergies)):
       NCi = ListofEnergies[i][2]
       NCSum = NCSum + NCi
       NCWeightedi = fBoltzSingleState(ListofEnergies[i][1], AvogadroR, Temp) * ListofEnergies[i][2] / BoltzPartition
       NCWeightedSum = NCWeightedSum + NCWeightedi
   # print "Weighted neighbour", int(round(NCWeightedSum))
   #NeighbourCount = int(round(NCSum/len(ListofEnergies)))
       NeighbourCount = round(NCWeightedSum, 1)
   # If we found dimers
   if breakDimer == "yes":
       print("####################################################")
       print("########### WARNING: Dimer formation?    ###########")
       print("####################################################")
       BoltzSumNi = 0
       BoltzPartition = 0
       BoltzMCSC = 0
       DeltapKMCSC = 99
       NeighbourCount = 0

   # Calculate the BornPenalty based on the neighbour count. It's a wrapper script for equation 13, 12, 11.
   EnergyBornPenalty = fEnergyBornPenalty(DieElecSpheDist, DieElecWaterDist, NeighbourCount, DieElecWater, DieElecCore, BornPenaltyB)
   DeltapKB = fDeltapK(EnergyBornPenalty, AvogadroR, Temp)

   # Do the calculations according to eq 3 and 9.
   pKm1 = fpKm1(DeltapKMCSC, pK1)
   pKm2 = fpKm2(DeltapKMCSC, DeltapKB, pK2)
   FracCysm1 = fFracCys(pKm1, pH)
   FracCysm2 = fFracCys(pKm2, pH)

   # Lets make a result file, and write out the angle, the SumMCSC, and the number of neighbours for this state.
   Currentdir = os.getcwd()
   Newdir = os.path.join(os.getcwd(), "Results")
   if not os.path.exists(Newdir):
       os.makedirs(Newdir)
   filename = os.path.join(".", "Results", "Result_" + molecule + "_" + chain + "_" + residue + ".txt")
   filenamelog = os.path.join(".", "Results", "Result_log.log")
   logfile = open(filenamelog, "a")
   outfile = open(filename, "w")
   timeforlog = strftime("%Y %b %d %a %H:%M:%S", localtime())
   logfile.write("# " + timeforlog + "\n")
   logfile.write("# molecule: %s , chain: %s, residue: %s %s, pH: %s " % (molecule, chain, residueName, residue, pH) + "\n")
   logfile.write("# BoltzSumNi:  BoltzPartition:  BoltzMCSC" + "\n")
   logfile.write(("# %.4f  %.4f  %.4f" + '\n') % (BoltzSumNi, BoltzPartition, BoltzMCSC))
   logfile.write("#    Res  NC    States  pKmcsc  pK1   pKB     pK2  pKm1     pKm2    f(C-)m1   f(C-)m2" + "\n")
   logfile.write(("; %s %s   %s  %s     %.4f  %s  %.4f  %s  %.4f  %.4f  %.6f  %.6f" + '\n') % (residueName, residue, NeighbourCount, AvailRotStates, DeltapKMCSC, pK1, DeltapKB, pK2, pKm1, pKm2, FracCysm1, FracCysm2))
   if nostates == "yes":
       logfile.write("##### ERROR; No states available ###" + "\n")
   if breakDimer == "yes":
       logfile.write("##### ERROR; Dimer formation ###" + "\n")
   logfile.write('\n')
   outfile.write("# molecule: %s , chain: %s, residue: %s %s, pH: %s " % (molecule, chain, residueName, residue, pH) + "\n")
   outfile.write("# BoltzSumNi:  BoltzPartition:  BoltzMCSC" + "\n")
   outfile.write(("# %.4f  %.4f  %.4f" + '\n') % (BoltzSumNi, BoltzPartition, BoltzMCSC))
   outfile.write("#    Res  NC    States  pKmcsc  pK1   pKB     pK2  pKm1     pKm2    f(C-)m1   f(C-)m2" + "\n")
   outfile.write(("; %s %s   %s  %s     %.4f  %s  %.4f  %s  %.4f  %.4f  %.6f  %.6f" + '\n') % (residueName, residue, NeighbourCount, AvailRotStates, DeltapKMCSC, pK1, DeltapKB, pK2, pKm1, pKm2, FracCysm1, FracCysm2))
   if nostates == "yes":
       outfile.write("##### ERROR; No states available ###" + "\n")
   if breakDimer == "yes":
       outfile.write("##### ERROR; Dimer formation ###" + "\n")
   outfile.write('\n')
   outfile.write("#Ang  SumMCSC   NC rNC MCSC_Weight       SG[X,Y,Z]" + "\n")
   for i in range(len(ListofEnergies)):
       outfile.write("%4.1d %10.3f %2.1d %1.1d %10.3f [%8.3f, %8.3f, %8.3f]" % (ListofEnergies[i][0], ListofEnergies[i][1], ListofEnergies[i][2], ListofEnergies[i][3], ListofEnergies[i][4], ListofEnergies[i][5][0], ListofEnergies[i][5][1], ListofEnergies[i][5][2]) + '\n')
   for i in range(len(ListofRotamerDiscarded)):
       outfile.write("%4.1d %10.3f %2.1d %1.1d %10.3f [%8.3f, %8.3f, %8.3f]" % (ListofRotamerDiscarded[i][0], ListofRotamerDiscarded[i][1], ListofRotamerDiscarded[i][2], ListofRotamerDiscarded[i][3], ListofRotamerDiscarded[i][4], ListofRotamerDiscarded[i][5][0], ListofRotamerDiscarded[i][5][1], ListofRotamerDiscarded[i][5][2]) + '\n')
   outfile.close()

   # Now, we are done. Just print out. The ; is for a grep command to select these lines in the output.
   print("# residue: %s %s. Average NeighbourCount NC= %s " % (residueName, residue, NeighbourCount))
   print("# From residue %s to residue %s" % (Minresidue, Maxresidue))
   print("# BoltzSumNi:  BoltzPartition:  BoltzMCSC")
   print("# %.4f  %.4f  %.4f" % (BoltzSumNi, BoltzPartition, BoltzMCSC))
   print("# Result written in file: %s" % (filename))
   print("#    Res  NC    States  pKmcsc  pK1   pKB     pK2  pKm1     pKm2    f(C-)m1   f(C-)m2")
   print("; %s %s   %s  %s     %.4f  %s  %.4f  %s  %.4f  %.4f  %.6f  %.6f" % (residueName, residue, NeighbourCount, AvailRotStates, DeltapKMCSC, pK1, DeltapKB, pK2, pKm1, pKm2, FracCysm1, FracCysm2))
   if nostates == "yes":
       print("##### ERROR; No states available ###")
   if breakDimer == "yes":
       print("##### ERROR; Dimer formation ###")

cmd.extend("cyspka", cyspka)


def loopcyspka(molecule, chain, residue, SeeProgress='no', pH=7.2, MoveSGatom='no', SGatom=str((0, 0, 0))): 

   residue = residue.split('.')
   residueList = []
   for i in residue:
       if '-' in i:
           tmp = i.split('-')
           residueList.extend(list(range(int(tmp[0]), int(tmp[-1]) + 1)))
       if '-' not in i:
           residueList.append(int(i))
   print("Looping over residues")
   print(residueList)
   for i in residueList:
       cyspka(molecule, chain, str(i), SeeProgress, pH, MoveSGatom, SGatom)

cmd.extend("loopcyspka", loopcyspka)


def fNeighbourCount(molecule, Cysmolecule, chain, residue, DieElecSpheDist): 

   nameselect = "(((/" + molecule + "//" + chain + " and not /" + molecule + "//" + chain + "/" + residue + ") or /" + molecule + "//" + chain + "/" + residue + "/N+CA+C+O)  within " + str(DieElecSpheDist) + " of /" + Cysmolecule + "//" + "/" + "/SG) and not resn HOH"
   # print nameselect
   cmd.select(residue + "NC", nameselect)
   # Adding 1 for CB
   Neighbours = cmd.count_atoms(residue + "NC") + 1
   cmd.delete(residue + "NC")
   return Neighbours


def fNeighbourWater(DieElecSpheDist, DieElecWaterDist, NeighbourCount): 

   Waters = 0.74 * math.pow(DieElecSpheDist, 3) / math.pow(DieElecWaterDist, 3) - NeighbourCount
   return Waters


def fDieElecEF(NeighbourWater, DieElecWater, NeighbourCount, DieElecCore): 

   DieElecEF = (NeighbourWater * DieElecWater + NeighbourCount * DieElecCore) / (NeighbourWater + NeighbourCount)
   return DieElecEF


def fBornPenalty(BornPenaltyB, DieElecEF, DieElecWater): 

   BornPenalty = (1.39 * math.pow(10, 6)) / (2 * BornPenaltyB) * (1.0 / DieElecEF - 1.0 / DieElecWater)
   return BornPenalty


def fEnergyBornPenalty(DieElecSpheDist, DieElecWaterDist, NeighbourCount, DieElecWater, DieElecCore, BornPenaltyB): 

   NeighbourWater = fNeighbourWater(DieElecSpheDist, DieElecWaterDist, NeighbourCount)
   DieElecEF = fDieElecEF(NeighbourWater, DieElecWater, NeighbourCount, DieElecCore)
   BornPenalty = fBornPenalty(BornPenaltyB, DieElecEF, DieElecWater)
   return BornPenalty


def fDeltapK(Energy, AvogadroR, Temp): 

   DeltapK = -1 * math.log10(math.exp(-Energy / (AvogadroR * Temp)))
   return DeltapK


def fRotateAroundLine(OriPoint, ThroughLinePoint, LineVector, AngleDeg): 

   # See http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/. Section 6.1
   AngleRad = math.radians(AngleDeg)
   x = OriPoint[0]
   y = OriPoint[1]
   z = OriPoint[2]
   a = ThroughLinePoint[0]
   b = ThroughLinePoint[1]
   c = ThroughLinePoint[2]
   u = LineVector[0]
   v = LineVector[1]
   w = LineVector[2]
   L = math.pow(u, 2) + math.pow(v, 2) + math.pow(w, 2)
   xPos = ((a * (math.pow(v, 2) + math.pow(w, 2)) - u * (b * v + c * w - u * x - v * y - w * z)) * (1 - math.cos(AngleRad)) + L * x * math.cos(AngleRad) + math.sqrt(L) * (-c * v + b * w - w * y + v * z) * math.sin(AngleRad)) / L
   yPos = ((b * (math.pow(u, 2) + math.pow(w, 2)) - v * (a * u + c * w - u * x - v * y - w * z)) * (1 - math.cos(AngleRad)) + L * y * math.cos(AngleRad) + math.sqrt(L) * (c * u - a * w + w * x - u * z) * math.sin(AngleRad)) / L
   zPos = ((c * (math.pow(u, 2) + math.pow(v, 2)) - w * (a * u + b * v - u * x - v * y - w * z)) * (1 - math.cos(AngleRad)) + L * z * math.cos(AngleRad) + math.sqrt(L) * (-b * u + a * v - v * x + u * y) * math.sin(AngleRad)) / L
   NewPos = [xPos, yPos, zPos]
   return NewPos


def fWSC(charge, DieElecSC, DistR): 

   # print charge, DistR
   WSC = 1.39 * math.pow(10, 6) * charge / (DieElecSC * DistR)
   return WSC


def fSumWSC(molecule, SGNameAngle, chain, residue, DieElecSC, SCchargeASP, SCchargeGLU, SCchargeOXT, SCchargeARG, SCchargeHIS, SCchargeLYS, SCchargeMET1, printSC): 

   SumWSC = 0.0
   SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
   # Sidechain ASP
   nameselect = "/" + molecule + " and resn ASP and name CG and not resi " + residue
   cmd.select("SC", nameselect)
   SClist = cmd.identify("SC")
   for i in range(len(SClist)):
       ResDist = cmd.dist(residue + 'distASP', SGnameselect, molecule + " and id " + str(SClist[i]))
       WSC = fWSC(SCchargeASP, DieElecSC, ResDist)
       SumWSC = SumWSC + WSC
       if printSC == 'yes':
           print("SC ASP ", str(SClist[i]), " ", SCchargeASP, " ", DieElecSC, " ", ResDist, " ", WSC)
   cmd.delete(residue + 'distASP')
   # Sidechain GLU
   nameselect = "/" + molecule + " and resn GLU and name CD and not resi " + residue
   cmd.select("SC", nameselect)
   SClist = cmd.identify("SC")
   for i in range(len(SClist)):
       ResDist = cmd.dist(residue + 'distGLU', SGnameselect, molecule + " and id " + str(SClist[i]))
       WSC = fWSC(SCchargeGLU, DieElecSC, ResDist)
       SumWSC = SumWSC + WSC
       if printSC == 'yes':
           print("SC GLU ", str(SClist[i]), " ", SCchargeGLU, " ", DieElecSC, " ", ResDist, " ", WSC)
   cmd.delete(residue + 'distGLU')
   # print "GLU", cmd.count_atoms("SC"), SumWSC
   # Sidechain OXT
   nameselect = "/" + molecule + " and byres name OXT and not resi " + residue
   cmd.select("SC", nameselect)
   cmd.select("SC", "SC and name C")
   SClist = cmd.identify("SC")
   for i in range(len(SClist)):
       ResDist = cmd.dist(residue + 'distOXT', SGnameselect, molecule + " and id " + str(SClist[i]))
       WSC = fWSC(SCchargeOXT, DieElecSC, ResDist)
       SumWSC = SumWSC + WSC
       if printSC == 'yes':
           print("SC OXT ", str(SClist[i]), " ", SCchargeOXT, " ", DieElecSC, " ", ResDist, " ", WSC)
   cmd.delete(residue + 'distOXT')
   # print "OXT", cmd.count_atoms("SC"), SumWSC
   # Sidechain ARG
   nameselect = "/" + molecule + " and resn ARG and name CZ and not resi " + residue
   cmd.select("SC", nameselect)
   SClist = cmd.identify("SC")
   for i in range(len(SClist)):
       ResDist = cmd.dist(residue + 'distARG', SGnameselect, molecule + " and id " + str(SClist[i]))
       WSC = fWSC(SCchargeARG, DieElecSC, ResDist)
       SumWSC = SumWSC + WSC
       if printSC == 'yes':
           print("SC ARG ", str(SClist[i]), " ", SCchargeARG, " ", DieElecSC, " ", ResDist, " ", WSC)
   cmd.delete(residue + 'distARG')
   # print "ARG", cmd.count_atoms("SC"), SumWSC
   # Sidechain HIS
   nameselect = "/" + molecule + " and resn HIS and name CD2 and not resi " + residue
   cmd.select("SC", nameselect)
   SClist = cmd.identify("SC")
   for i in range(len(SClist)):
       ResDist = cmd.dist(residue + 'distHIS', SGnameselect, molecule + " and id " + str(SClist[i]))
       WSC = fWSC(SCchargeHIS, DieElecSC, ResDist)
       SumWSC = SumWSC + WSC
       if printSC == 'yes':
           print("SC HIS ", str(SClist[i]), " ", SCchargeHIS, " ", DieElecSC, " ", ResDist, " ", WSC)
   cmd.delete(residue + 'distHIS')
   # print "HIS", cmd.count_atoms("SC"), SumWSC
   # Sidechain LYS
   nameselect = "/" + molecule + " and resn LYS and name NZ and not resi " + residue
   cmd.select("SC", nameselect)
   SClist = cmd.identify("SC")
   for i in range(len(SClist)):
       ResDist = cmd.dist(residue + 'distLYS', SGnameselect, molecule + " and id " + str(SClist[i]))
       WSC = fWSC(SCchargeLYS, DieElecSC, ResDist)
       SumWSC = SumWSC + WSC
       if printSC == 'yes':
           print("SC LYS ", str(SClist[i]), " ", SCchargeLYS, " ", DieElecSC, " ", ResDist, " ", WSC)
   cmd.delete(residue + 'distLYS')
   # print "LYS", cmd.count_atoms("SC"), SumWSC
   # Sidechain MET1
   nameselect = "/" + molecule + " and resn MET and res 1 and not resi " + residue
   cmd.select("SC", nameselect)
   cmd.select("SC", "SC and name N")
   SClist = cmd.identify("SC")
   for i in range(len(SClist)):
       ResDist = cmd.dist(residue + 'distMET1', SGnameselect, molecule + " and id " + str(SClist[i]))
       WSC = fWSC(SCchargeMET1, DieElecSC, ResDist)
       SumWSC = SumWSC + WSC
       if printSC == 'yes':
           print("SC MET1 ", str(SClist[i]), " ", SCchargeMET1, " ", DieElecSC, " ", ResDist, " ", WSC)
   cmd.delete(residue + 'distMET1')
   # print "MET1", cmd.count_atoms("SC"), SumWSC
   cmd.delete("SC")
   return SumWSC


def fWMC(charge, DieElecMC, DistR): 

   WMC = 1.39 * math.pow(10, 6) * charge / (DieElecMC * DistR)
   return WMC


def fSumWMCFirst(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeC, MCchargeO, printMC): 

   # print "First", MCNeighbour
   SumWMCFirst = 0.0
   SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
   NBnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour)
   cmd.select("MC", NBnameselect)
   MCpdbstr = cmd.get_pdbstr("MC")
   MCsplit = MCpdbstr.split()
   residueName = MCsplit[3]
   # print NBnameselect, residueName
   # Mainchain C
   Cnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/C"
   ResDist = cmd.dist(residue + 'distFirstC', SGnameselect, Cnameselect)
   WMC = fWMC(MCchargeC, DieElecMC, ResDist)
   SumWMCFirst = SumWMCFirst + WMC
   if printMC == 'yes':
       print("MC C ", MCNeighbour, " ", MCchargeC, " ", DieElecMC, " ", ResDist, " ", WMC)
   # Mainchain O
   Onameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/O"
   ResDist = cmd.dist(residue + 'distFirstO', SGnameselect, Onameselect)
   WMC = fWMC(MCchargeO, DieElecMC, ResDist)
   SumWMCFirst = SumWMCFirst + WMC
   if printMC == 'yes':
       print("MC O ", MCNeighbour, " ", MCchargeO, " ", DieElecMC, " ", ResDist, " ", WMC)
   cmd.delete(residue + 'distFirstC')
   cmd.delete(residue + 'distFirstO')
   cmd.delete("MC")
   return SumWMCFirst


def fSumWMCresidue(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeC, MCchargeO, MCchargeN, MCchargeH, AmideName, printMC): 

   # print "residue", MCNeighbour
   SumWMCresidue = 0.0
   SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
   NBnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour)
   cmd.select("MC", NBnameselect)
   MCpdbstr = cmd.get_pdbstr("MC")
   MCsplit = MCpdbstr.split()
   residueName = MCsplit[3]
   # print NBnameselect, residueName
   AmideProt = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/H01"
   Hnameselect = "/" + AmideName + "//" + chain + "/" + str(MCNeighbour) + "/H01"
   if cmd.count_atoms(AmideProt) == 0 and cmd.count_atoms(Hnameselect) == 0:
       HbuildSelect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
       cmd.h_add(HbuildSelect)
       cmd.create(AmideName, AmideName + " + " + AmideProt)
       cmd.remove(AmideProt)
   # Mainchain AmideH
   ResDist = cmd.dist(residue + 'distResH', SGnameselect, Hnameselect)
   WMC = fWMC(MCchargeH, DieElecMC, ResDist)
   SumWMCresidue = SumWMCresidue + WMC
   if printMC == 'yes':
       print("MC H ", MCNeighbour, " ", MCchargeH, " ", DieElecMC, " ", ResDist, " ", WMC)
   # Mainchain C
   Cnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/C"
   ResDist = cmd.dist(residue + 'distResC', SGnameselect, Cnameselect)
   WMC = fWMC(MCchargeC, DieElecMC, ResDist)
   SumWMCresidue = SumWMCresidue + WMC
   if printMC == 'yes':
       print("MC C ", MCNeighbour, " ", MCchargeC, " ", DieElecMC, " ", ResDist, " ", WMC)
   # Mainchain O
   Onameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/O"
   ResDist = cmd.dist(residue + 'distResO', SGnameselect, Onameselect)
   WMC = fWMC(MCchargeO, DieElecMC, ResDist)
   SumWMCresidue = SumWMCresidue + WMC
   if printMC == 'yes':
       print("MC O ", MCNeighbour, " ", MCchargeO, " ", DieElecMC, " ", ResDist, " ", WMC)
   # Mainchain N
   Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
   ResDist = cmd.dist(residue + 'distResN', SGnameselect, Nnameselect)
   WMC = fWMC(MCchargeN, DieElecMC, ResDist)
   SumWMCresidue = SumWMCresidue + WMC
   if printMC == 'yes':
       print("MC N ", MCNeighbour, " ", MCchargeN, " ", DieElecMC, " ", ResDist, " ", WMC)
   cmd.delete(residue + 'distResH')
   cmd.delete(residue + 'distResC')
   cmd.delete(residue + 'distResO')
   cmd.delete(residue + 'distResN')
   cmd.show("nb_spheres", AmideName)
   cmd.delete("MC")
   return SumWMCresidue


def fSumWMCLast(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeN, MCchargeH, MCchargeProCA, MCchargeProCD, MCchargeProN, AmideName, printMC): 

   # print "Last", MCNeighbour
   SumWMCLast = 0.0
   SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
   NBnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour)
   cmd.select("MC", NBnameselect)
   MCpdbstr = cmd.get_pdbstr("MC")
   MCsplit = MCpdbstr.split()
   residueName = MCsplit[3]
   # print NBnameselect, residueName
   if residueName == "PRO":
       # Proline CA
       CAnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/CA"
       ResDist = cmd.dist(residue + 'distLastProCA', SGnameselect, CAnameselect)
       WMC = fWMC(MCchargeProCA, DieElecMC, ResDist)
       SumWMCLast = SumWMCLast + WMC
       if printMC == 'yes':
           print("MC ProCA ", MCNeighbour, " ", MCchargeProCA, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Proline CD
       CDnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/CD"
       ResDist = cmd.dist(residue + 'distLastProCD', SGnameselect, CDnameselect)
       WMC = fWMC(MCchargeProCD, DieElecMC, ResDist)
       SumWMCLast = SumWMCLast + WMC
       if printMC == 'yes':
           print("MC ProCD ", MCNeighbour, " ", MCchargeProCD, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Proline N
       Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
       ResDist = cmd.dist(residue + 'distLastProN', SGnameselect, Nnameselect)
       WMC = fWMC(MCchargeProN, DieElecMC, ResDist)
       SumWMCLast = SumWMCLast + WMC
       if printMC == 'yes':
           print("MC ProN ", MCNeighbour, " ", MCchargeProN, " ", DieElecMC, " ", ResDist, " ", WMC)
   else:
       AmideProt = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/H01"
       Hnameselect = "/" + AmideName + "//" + chain + "/" + str(MCNeighbour) + "/H01"
       if cmd.count_atoms(AmideProt) == 0 and cmd.count_atoms(Hnameselect) == 0:
           HbuildSelect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
           cmd.h_add(HbuildSelect)
           cmd.create(AmideName, AmideName + " + " + AmideProt)
           cmd.remove(AmideProt)
       # Mainchain AmideH
       ResDist = cmd.dist(residue + 'distLastH', SGnameselect, Hnameselect)
       WMC = fWMC(MCchargeH, DieElecMC, ResDist)
       SumWMCLast = SumWMCLast + WMC
       if printMC == 'yes':
           print("MC H ", MCNeighbour, " ", MCchargeH, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Mainchain N
       Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
       ResDist = cmd.dist(residue + 'distLastN', SGnameselect, Nnameselect)
       WMC = fWMC(MCchargeN, DieElecMC, ResDist)
       SumWMCLast = SumWMCLast + WMC
       if printMC == 'yes':
           print("MC N ", MCNeighbour, " ", MCchargeN, " ", DieElecMC, " ", ResDist, " ", WMC)
   cmd.delete(residue + 'distLastProCA')
   cmd.delete(residue + 'distLastProCD')
   cmd.delete(residue + 'distLastProN')
   cmd.delete(residue + 'distLastH')
   cmd.delete(residue + 'distLastN')
   cmd.show("nb_spheres", AmideName)
   cmd.delete("MC")
   return SumWMCLast


def fSumWMC(molecule, SGNameAngle, chain, residue, MCNeighbour, DieElecMC, MCchargeC, MCchargeO, MCchargeN, MCchargeH, MCchargeProCA, MCchargeProCD, MCchargeProN, AmideName, printMC): 

   # print "chain", MCNeighbour
   SumWMC = 0.0
   SGnameselect = "/" + SGNameAngle + "//" + "/" + "/SG"
   NBnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour)
   cmd.select("MC", NBnameselect)
   MCpdbstr = cmd.get_pdbstr("MC")
   MCsplit = MCpdbstr.split()
   residueName = MCsplit[3]
   # print NBnameselect, residueName
   if residueName == "PRO":
       # Proline CA
       CAnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/CA"
       ResDist = cmd.dist(residue + 'distProCA', SGnameselect, CAnameselect)
       WMC = fWMC(MCchargeProCA, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC ProCA ", MCNeighbour, " ", MCchargeProCA, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Proline CD
       CDnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/CD"
       ResDist = cmd.dist(residue + 'distProCD', SGnameselect, CDnameselect)
       WMC = fWMC(MCchargeProCD, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC ProCD ", MCNeighbour, " ", MCchargeProCD, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Proline N
       Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
       ResDist = cmd.dist(residue + 'distProN', SGnameselect, Nnameselect)
       WMC = fWMC(MCchargeProN, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC ProN ", MCNeighbour, " ", MCchargeProN, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Proline C
       Cnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/C"
       ResDist = cmd.dist(residue + 'distProC', SGnameselect, Cnameselect)
       WMC = fWMC(MCchargeC, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC ProC ", MCNeighbour, " ", MCchargeC, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Proline O
       Onameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/O"
       ResDist = cmd.dist(residue + 'distProO', SGnameselect, Onameselect)
       WMC = fWMC(MCchargeO, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC ProO ", MCNeighbour, " ", MCchargeO, " ", DieElecMC, " ", ResDist, " ", WMC)
   else:
       AmideProt = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/H01"
       Hnameselect = "/" + AmideName + "//" + chain + "/" + str(MCNeighbour) + "/H01"
       if cmd.count_atoms(AmideProt) == 0 and cmd.count_atoms(Hnameselect) == 0:
           HbuildSelect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
           cmd.h_add(HbuildSelect)
           cmd.create(AmideName, AmideName + " + " + AmideProt)
           cmd.remove(AmideProt)
       # Mainchain AmideH
       ResDist = cmd.dist(residue + 'distH', SGnameselect, Hnameselect)
       WMC = fWMC(MCchargeH, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC H ", MCNeighbour, " ", MCchargeH, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Mainchain C
       Cnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/C"
       ResDist = cmd.dist(residue + 'distC', SGnameselect, Cnameselect)
       WMC = fWMC(MCchargeC, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC C ", MCNeighbour, " ", MCchargeC, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Mainchain O
       Onameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/O"
       ResDist = cmd.dist(residue + 'distO', SGnameselect, Onameselect)
       WMC = fWMC(MCchargeO, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC O ", MCNeighbour, " ", MCchargeO, " ", DieElecMC, " ", ResDist, " ", WMC)
       # Mainchain N
       Nnameselect = "/" + molecule + "//" + chain + "/" + str(MCNeighbour) + "/N"
       ResDist = cmd.dist(residue + 'distN', SGnameselect, Nnameselect)
       WMC = fWMC(MCchargeN, DieElecMC, ResDist)
       SumWMC = SumWMC + WMC
       if printMC == 'yes':
           print("MC N ", MCNeighbour, " ", MCchargeN, " ", DieElecMC, " ", ResDist, " ", WMC)
   cmd.delete(residue + 'distProCA')
   cmd.delete(residue + 'distProCD')
   cmd.delete(residue + 'distProN')
   cmd.delete(residue + 'distProC')
   cmd.delete(residue + 'distProO')
   cmd.delete(residue + 'distH')
   cmd.delete(residue + 'distC')
   cmd.delete(residue + 'distO')
   cmd.delete(residue + 'distN')
   cmd.show("nb_spheres", AmideName)
   cmd.delete("MC")
   return SumWMC


def fBoltzSingleState(SumMCSC, AvogadroR, Temp): 

   BoltzSingleState = math.exp(-SumMCSC / (AvogadroR * Temp))
   return BoltzSingleState


def fpKm1(DeltapKMCSC, pK1): 

   pKm1 = DeltapKMCSC + pK1
   return pKm1


def fpKm2(DeltapKMCSC, DeltapKB, pK2): 

   pKm2 = DeltapKMCSC + DeltapKB + pK2
   return pKm2


def fFracCys(pKm, pH): 

   FracCys = 1.0 / (math.pow(10, (pKm - pH)) + 1)
   return FracCys


def AtomVector(AtomStart, AtomEnd): 

   PosStart = cmd.get_atom_coords(AtomStart)
   PosEnd = cmd.get_atom_coords(AtomEnd)
   VectorDiff = [(PosEnd[0] - PosStart[0]), (PosEnd[1] - PosStart[1]), (PosEnd[2] - PosStart[2])]
   return VectorDiff


def pairfitCys(Cysmolecule, molecule, chain, residue): 

   RN = "/" + Cysmolecule + "//" + "/" + "/N"
   PN = "/" + molecule + "//" + chain + "/" + residue + "/N"
   RCA = "/" + Cysmolecule + "//" + "/" + "/CA"
   PCA = "/" + molecule + "//" + chain + "/" + residue + "/CA"
   RC = "/" + Cysmolecule + "//" + "/" + "/C"
   PC = "/" + molecule + "//" + chain + "/" + residue + "/C"
   RCB = "/" + Cysmolecule + "//" + "/" + "/CB"
   PCB = "/" + molecule + "//" + chain + "/" + residue + "/CB"
   cmd.select("CBatom", PCB)
   CBatomNr = cmd.count_atoms("CBatom")
   # If PRO or GLY, then only fit N, CA, C atoms
   if CBatomNr == 0:
       cmd.pair_fit(RN, PN, RCA, PCA, RC, PC)
   else:
       # cmd.pair_fit(RN,PN,RCA,PCA,RC,PC,RCB,PCB)
       cmd.pair_fit(RN, PN, RCA, PCA, RC, PC)
   cmd.delete("CBatom")


def CheckDimer(dihedSG, molecule, chain, residue): 

   breakDimer = "no"
   nameselect = "(/" + molecule + "//" + chain + " and name SG and not /" + molecule + "//" + chain + "/" + residue + ") within 5 of " + dihedSG
   cmd.select(str(molecule) + str(residue) + "Dimer", nameselect)
   DimerSG = cmd.count_atoms(str(molecule) + str(residue) + "Dimer")
   if DimerSG > 0:
       print("####################################################")
       print("########### WARNING: SG in near detected ###########")
       print("########### Is this a dimer?             ###########")
       print("####################################################")
       cmd.select(str(molecule) + str(residue) + "Dimer", "byres " + str(molecule) + str(residue) + "Dimer")
       cmd.show("sticks", str(molecule) + str(residue) + "Dimer")
       breakDimer = "yes"
   else:
       cmd.delete(str(molecule) + str(residue) + "Dimer")
   return breakDimer
Retrieved from "http://pymolwiki.org/index.php?title=User:Tlinnet&oldid=12233"