Code
-- Ligand docker v2.3
-- Credits: Raven_pl,Susume,LociOiling,lynnai,
-- robgee,jeff101
-- dev version : Ligand docker BiS RW v1.5
-- It makes BiS bands to move the Ligand thru space,
-- hopefully finding a better position,
-- every so often it moves the Ligand extra far
-- to try and optimize its search strategy for finding a better position.
-- This type of search is based on a Levy flight (e should be accented on Levy),
-- Sharks and Albatrosses use it.
-- It is not a proper random walk/Levy flight because
-- frankly i dont know how to do that with spacebands.
-- To Do:
-- Band X,Y atoms ????
-- toggle for wiggle(backbone),wiggle(all)
-- ------------------------------------
-- 7/01/2020 for v2.1b, jeff101 changed input upper limits
-- for Wiggles to 500 and for minDist & maxDist to 50000.
-- Setting minDist to 3000 & maxDist to 25000
-- lets the actual band lengths used range from 3-25 Angstroms.
-- This recipe seems to assume the protein's final segment
-- is the one and only ligand.
-- 7/16/2020 for v2.2, jeff101 allowed more than 2 bands
-- to the ligand. Unclear if Levy random walk will still work
-- as intended, even using just 2 bands.
-- 8/07/2020 for v2.3, jeff101 let # cycles between reverts
-- and Levy use vary. Changed minDist & maxDist input upper
-- limits to 150000. Changed many default values.
-- 8/09/2020 for v2.3, jeff101 changed goal length for bands
-- from 3.5 (the default value) to 0.
-- Ideas for later:
-- Use textboxes not sliders to input minDist and maxDist.
-- Let user pick default random seed or input his own value.
-- Break up wiggling without bands into shorter intervals.
-- If gain doesn't pass a threshold, quit wiggling.
-- If gain passes a threshold, continue wiggling.
-- last updated 8/10/2020 657am
-- ------------------------------------
-- ---------- GLOBAL VARS -------------
-- ------------------------------------
Recipe = "Ligand docker"
Version = "v2.3"
ReVersion = Recipe .. " " .. Version
OrigBands = band.GetCount() -- # of user-supplied bands at start
-- -
-- ------------- OPTIONS --------------
-- -
Cycles = 1000 --was 30 before 8/7/2020.
RevertCycles = 2 --# of cycles between revert to best score.
--0 means never revert to best score.
--was 5 before 8/9/2020.
LevyCycles = 3 --# of cycles between uses of Levy.
--0 means never use Levy.
--was 8 before 8/9/2020.
WigWithBands = 30 --how much to wiggle with bands. was 1 before 8/7/2020.
Wiggles = 200 --how much to wiggle with no bands. was 3 before 8/7/2020.
AtomsInLigand = 0 --# of atoms in the ligand.
BandsToLigand = 5 --# of bands from ligand to space. was 2 before 8/7/2020.
minDist = 5 --was 10 before 8/7/2020.
maxDist = 50
Levy = 0
BandStrength = 1.50
CapCI = false
userSetMaxCI = behavior.GetClashImportance()
Verbosity = 0 -- 0 = the least output messages, 2 = the most output messages
-- END OPTIONS
segCnt = structure.GetCount()
numSegs =structure.GetCount()
highscore = 0
BASESCORE = current.GetEnergyScore()
-- -------------------------------------------
-- List all the recipe's functions below:
-- -------------------------------------------
function HasLigandCheck( )
alwaysLigand = false
segCnt2 = segCnt
while structure.GetSecondaryStructure(segCnt2)=="M" do
segCnt2 = segCnt2-1
end -- while
-- line below looks odd to jeff101 8/7/2020
HASLIGAND = segCnt2 < segCnt
-- is line above the same as the line below?
-- HASLIGAND = (segCnt2 < segCnt)
if HASLIGAND then
LIGAND = segCnt
AtomsInLigand = structure.GetAtomCount(LIGAND)
alwaysLigand = true
end -- if HASLIGAND
end -- HasLigandCheck()
function round ( x )--to 3 places
return x - x % 0.001
end -- round()
-- -----------------------------------------------------
function CI( RecipeCI ) -- Cap CI else relative CI
local Using_CI = RecipeCI
if CapCI == true then
if RecipeCI > userSetMaxCI then
Using_CI = userSetMaxCI
end
else
Using_CI = RecipeCI * userSetMaxCI -- relative CI
end
if Using_CI < 0 or Using_CI > 1 then --sanity check
print( 'CI ERROR, CI now set to 1' )
Using_CI = 1
end
--print('Recipe set CI '..RecipeCI..', Using_CI '..Using_CI )
behavior.SetClashImportance( Using_CI )
end-- Thnx to TvdL
function SeedRandom( )
local seed = os.time ( ) / math.abs ( current.GetEnergyScore ( ) )
seed = seed % 0.001
seed = 1 / seed
while seed < 10000000 do
seed = seed * 1000
end
while seed > 2 ^ 32 do
seed = seed / 10
end--LO
seed = seed - seed % 1
math.randomseed( seed )
math.random( ); math.random( ); math.random( ) --See stackoverflow
end-- Rav3n_pl, LociOiling
-- ------------------------------------------------------
function Cleanup(err) -- cancel recipe function
if CLEANUPENTRY ~= nil then
return
end
CLEANUPENTRY = true --LO
print( '========================' )
print( ReVersion..': Cleaning Up.' )
save.Quickload( 2 )
selection.DeselectAll()
delete_temp_bands()
behavior.SetClashImportance( userSetMaxCI )
print( "Final score: "..round( highscore ) )
print( "Total gain: " ..round( highscore - BASESCORE ) )
print( "Error "..err) -- ..err prints no space before err
-- print( "Error",err) -- ,err prints tab before err
end -- Cleanup()
-- -----------------------------------------
function delete_temp_bands()
local TempBands=band.GetCount()
for bnum=TempBands, OrigBands+1, -1 do
band.Delete(bnum)
end -- for bnum
end -- delete_temp_bands()
-- -----------------------------------------
function random_bands()
local hscore = current.GetScore()
print("Cycle start score: "..round( hscore ) )
local OG = LIGAND
if BandStrength >10.0 then BandStrength = 9.99 end
minDist = minDist + Levy
maxDist = maxDist + Levy
-- minDist, maxDist, Levy are all in 0.001 Angstroms above
local HasBand={} -- an array of 0's or 1's with length BandsToLigand
local origval=0
local ngot=0
local ntoget=BandsToLigand
if 2*BandsToLigand>AtomsInLigand then
origval=1
ntoget=AtomsInLigand-BandsToLigand
end
for anum=1, AtomsInLigand do
HasBand[anum]=origval
end -- for anum
if Verbosity>1 then
if origval==0 then
print('Banding the '..ntoget..' ligand atoms below:')
else
print('Banding all but the '..ntoget..' ligand atoms below:')
end
end -- if Verbosity
while ngot<ntoget do
anum=math.random(1,AtomsInLigand)
if HasBand[anum]==origval then
ngot=ngot+1
HasBand[anum]=1-origval
if Verbosity>1 then
print(string.format('(%3d) %3d',ngot,anum))
end -- if Verbosity
end -- if HasBand
end -- while ngot
-- ngot should equal ntoget here
-- make all BandsToLigand bands below
if Verbosity>=0 then
print('Levy,minDist,maxDist='..Levy..','..minDist..','..maxDist..' in 0.001 Angstroms')
end -- if Verbosity
ngot=0
for anum=1, AtomsInLigand do
if HasBand[anum]==1 then
ngot=ngot+1
local dst = math.random( minDist, maxDist )
local RHO=(dst / 1000) --RHO is a distance in Angstroms
-- while Levy, dst, minDist, maxDist are distances in 0.001 Angstroms
--print('Distance '..RHO)
local theta = math.acos(math.random( ) )
local phi = 2 * math.pi * math.random( )
repeat
Xaxis=math.random( 1, numSegs-1 )
Yaxis=math.random( 1, numSegs-1 )
until OG ~= Xaxis and Yaxis ~= Xaxis and OG ~= Yaxis
--print(OG,Xaxis,Yaxis) all 3 must be different segment #'s
-- {segmentOrigin, segmentXAxis, segmentYAxis, rho, theta, phi, Oatom, Xatom, Yatom}
-- Xatom & Yatom below are for alpha-carbons on segments Xaxis & Yaxis
local Xatom = 0
local Yatom = 0
local OriginAtom = anum
local qed = band.Add(OG, Xaxis, Yaxis, RHO, theta, phi, OriginAtom, Xatom, Yatom )
-- above gives 3.5 as the default goal length for band qed
-- above gives 1.0 as the default strength for band qed
-- RHO is the initial actual length for band qed
if qed >= 1 then
band.SetGoalLength(qed,0) -- set band qed's goal length to 0
band.SetStrength(qed,BandStrength)
if Verbosity>0 then
print(string.format('(%3d) Band ligand atom %3d with %7.3f Angstroms.',
ngot,anum,RHO))
end -- if Verbosity
end -- if qed
end -- if HasBand
end -- for anum
-- ngot should equal BandsToLigands here
-- make all BandsToLigand bands above
selection.DeselectAll ()
selection.Select ( OG ) -- selects the ligand
CI(0.88)
local witers = WigWithBands
--structure.ShakeSidechainsSelected(witers)
--structure.WiggleSelected(witers, true, true) -- moves ligand through space and bends it.
--structure.WiggleSelected(witers, true, false) -- moves ligand through space but does not bend it.
structure.WiggleSelected(witers, true, true) -- moves ligand through space and bends it.
-- -----------------------------------------------------------------------------------------------------
delete_temp_bands()
if Verbosity>=0 then
print('After banded wiggle: '..round(current.GetEnergyScore()))
end
CI(1.0)
--structure.ShakeSidechainsSelected( 2 )
structure.WiggleAll( Wiggles , true, true )
-- above moves ligand thru space and lets it bend.
-- above also moves protein backbone and sidechains.
newscore = current.GetEnergyScore()
print('After unbanded wiggle: '..round(newscore))
if newscore > highscore then
save.Quicksave( 2 )
recentbest.Save( )
local gain = newscore - highscore
print('gain: '..round(gain))
highscore = newscore
print("++ High score = "..round(highscore))
end
selection.DeselectAll()
print('')
end -- random_bands()
function LigandDock()
highscore = current.GetScore( )
local revertscore = current.GetScore( )
for i = 1, Cycles do
Levy = 0
revertscore = current.GetScore( )
-- RevertCycles == 0 should skip the following every cycle i
if RevertCycles ~= 0 then
if i % RevertCycles == 0 and revertscore <= highscore then
print(' -- reverting to best position --')
save.Quickload( 2 )
end
end
-- LevyCycles == 0 should skip the following every cycle i
if LevyCycles ~= 0 then
if i % LevyCycles == 0 and revertscore <= highscore then
print(' -- Implementing Levy flight --')
Levy = math.random( 500, 800 )
-- Levy should be an integer from 500 to 800 now
end
end
print('Cycle '..i..' of '..Cycles..' at '..os.date()..':')
random_bands()
end -- for i
end -- LigandDock()
--
--- ----------------------------------------------------------
--
function Prelim()
print( ReVersion )
save.Quicksave( 1 ) --OG pose for backup
save.Quicksave( 2 ) --for high score
print('Start pose in Quicksave 1')
print('Best pose in Quicksave 2')
print("Levy flight version, with random distance")
print("Protein length: "..structure.GetCount() )
print('# of Ligand Atoms: '..AtomsInLigand)
print('# of User-Supplied Bands: '..OrigBands)
print('')
print('Total # of Cycles to do: '..Cycles)
print('# of Cycles per Revert: '..RevertCycles..' (0 means never Revert)')
print('# of Cycles per Levy: '..LevyCycles..' (0 means never use Levy)')
print('# of Wiggles with Bands per Cycle: '..WigWithBands)
print('# of Wiggles without Bands per Cycle: '..Wiggles)
print('# of Bands from Space to Ligand: '..BandsToLigand)
print('Minimum Band Length: '..minDist..' ('..(minDist/1000)..' Angstroms)')
print('Maximum Band Length: '..maxDist..' ('..(maxDist/1000)..' Angstroms)')
-- 8/10/2020 jeff101 thinks minDist & maxDist are given in 0.001 Angstroms
-- so that math.random(minDist,maxDist) gives an integer from minDist to maxDist
print('Band Strength: '..BandStrength )
print('Maximum CI recipe will use: '..userSetMaxCI )
print('Verbosity: '..Verbosity)
print("Recipe start score: "..round(current.GetEnergyScore()))
print('')
end -- Prelim()
function DialogOptions()
local dlog = dialog.CreateDialog(ReVersion)
dlog.iters = dialog.AddSlider("Cycles",Cycles,1,1000,0)
dlog.revcyc = dialog.AddSlider("Revert Cycles",RevertCycles,0,10,0)
dlog.levycyc = dialog.AddSlider("Levy Cycles",LevyCycles,0,10,0)
dlog.wwb = dialog.AddSlider("Banded Wiggles",WigWithBands,1,500,0)
dlog.wigls = dialog.AddSlider("Free Wiggles",Wiggles,1,500,0)
dlog.spacer1 = dialog.AddLabel('')
dlog.bnd2lig = dialog.AddSlider('# Ligand Bands',BandsToLigand,1,AtomsInLigand,0)
dlog.minidist = dialog.AddSlider('Min Band Length',minDist, 5,150000,0)
dlog.maxidist = dialog.AddSlider('Max Band Length',maxDist,20,150000,0)
dlog.bandstren = dialog.AddSlider('Band Strength',BandStrength,0.1,5,1)
dlog.usermaxci = dialog.AddSlider("Maximum CI",userSetMaxCI,0.00,1,2)
dlog.spacer3 = dialog.AddLabel('')
dlog.vrb = dialog.AddSlider('Verbosity',Verbosity,0,2,0)
dlog.ok = dialog.AddButton("OK",1)
dlog.cancel = dialog.AddButton("Cancel",0)
dialogresult = dialog.Show(dlog)
if dialogresult > 0 then
Cycles = dlog.iters.value
RevertCycles = dlog.revcyc.value
LevyCycles = dlog.levycyc.value
WigWithBands = dlog.wwb.value
Wiggles = dlog.wigls.value
BandsToLigand = dlog.bnd2lig.value
minDist = dlog.minidist.value
maxDist = dlog.maxidist.value
BandStrength = dlog.bandstren.value
userSetMaxCI = dlog.usermaxci.value
Verbosity = dlog.vrb.value
end
return dialogresult
end -- DialogOptions()
function Main()
SeedRandom()
HasLigandCheck( )
if HASLIGAND == false then
print('This is NOT a LIGAND puzzle ... exiting')
return
end
local go = DialogOptions( )
if go == 1 then
Prelim( ) --print the settings
LigandDock( )
else
print( 'Cancelled' )
return
end
-- restore best score and print it below
save.Quickload( 2 )
local tmpscore=current.GetScore()
print('Restoring best score: '..round(tmpscore))
undo.SetUndo( true )
behavior.SetClashImportance( userSetMaxCI ) -- return to user setting
end -- Main()
xpcall( Main, Cleanup )
--Ligand docker v2.3
