Code
--iwdn hbnw
--version 1.3 28.02.2022
--H-Bond NetWorker: Script to semi-automatic find H-Bond networks.
--What it does:
--This recipe tries to extend the excellent work of spvincent's H-Network Probe.
--In hbnw not the H-Bond network-bonus is used as criterion to valuate if a net has formed but the function "structure.GetHBonds()" is used instead.
--See the post "https://fold.it/portal/node/2007744" for more details.
--This recipe will try to bond sidechains together as much as possible.
--It uses direct mutation and then shaking and wiggling for this.
--In order to use this you should first freeze all sidechains (you can use for example the public recipe "iwdn ufreeze SC 1.0" for this) and then de-freeze the ones that you want this recipe to act on.
--Basically, the unfrozen sidechains should be mostly core-residues. Especially in symmetry-puzzles, these should be the core-residues at the interface between monomer and symmetric copies (plus maybe some surrounding additional sidechains).
--Since version 1.2 there are different modes which can be used. The previous mode in older versions is now called single-mode and the new mode is called Auto-Grow (AG).
--I recommend using AG in the future. Especially for symmetry-puzzles with HBond-bonus there is a better chance to find HBbonus-nets than for single-mode. Got astonishing results with AG.
--In AG-mode it is tried in multiple iterations (iter) to first find a basic net and then to extend it in further iterations. Bonds are only accepted and incorporated into the net when they are connecting to the net from the previous iteration.
--Like this, the net will grow in AG-mode all over the structure.
--Be aware though that there is no valuation of the net done by hbnw. So it may find you nets with open ends. You will have to manually optimize your net by removing unwanted elements from the net. This is possible in most cases and leads to higher-quality nets which then give good HBbonus.
--But in any case the AG-mode will give you a good impression of where it is promising to build up your net.
--In single-mode it is gone randomly through all possible sidechain-mutation-poses and detected (after shake&wiggle) with "structure.GetHBonds()" if one or more new HBonds have formed. Duplicate poses are avoided.
--The nice thing is then that it is checked if such a net was already found. This is done by comparing the residues on which it has formed and also the atoms of the sidechains on which it has formed.
--This will reduce the number of found nets but it increases dramatically the quality of the found nets.
--Be aware though that some nets may not be kept. There are sidechains which have the same atoms at the same sidechain-position and only ONE of these nets will be kept.
--So it might be worthwile to play around with the found nets and maybe change one or more sidechains in it to an "equivalent" mutations.
--The stored nets are quick-saved and after the recipe has finished they are loaded into undo-stack so that you can flip through them via ctrl+z/y (thanks again spvincent for this nice function!).
--The ultimate goal of this recipe for me is to find AND EXTEND nets automatically. This is now implemented in a first step with the AG-mode.
--In single-mode, you should run the recipe on a few residues. Then check the found nets and re-run this recipe on non-bonded residues in the vicinity of the found net. Like this your net will grow.
--In AG-mode, this is done automatically for you. You just need to specify an initial center-segment around which your net shall grow (or you can leave the option "AG center seg(s)" set to 0 and hbnw will randomly choose a center-segment for you. You can watch the net grow then in the further iterations.
--I got very good results in AG-mode on puzzle 1980. But you will probably have to hand-optimize the found nets and your protein in general to get rid of bad-scoring configurations in the found net.
--It should be noted though that it may be better to fist find some initial nets and then try to extend them by re-running this recipe. In AG-mode you have the option "AG init base seg(s)" for this. In this you can specify the segments of an already present net and this will be used in further iterations to connect to. So this initial net will be extended. So if you have a net with e.g. HBbonus in a symmetry puzzle, you can use this to extend this net and get better HBbonus. In this case, the center segment should be chosen to one of the initial net segments.
--This recipe is general so it can be used on any kind of puzzle and not only ones with H-Bond bonus! This is very nice I think and I thank the developers very much for implementing "structure.GetHBonds()"!
--Be aware that some parts of this recipe still need speed-optimization. Sometimes it looks like it may have stalled but it should go on after a while. In single-mode, the duplicate-check may take quite a bit of time especially if already many nets were found.
--The function structure.getHBonds() itself is not the fastest but that is nothing I can do something about.
--Please be aware that there may be bugs in this recipe since due to the increased number of options now it was not possible to debug and test everything in detail. Getting this core-functionality going was already quite painsome due to its complexity and the plan is to improve this recipe over time.
--Nevertheless I got very good and promising results especially on symmetry-puzzles with HBbonus.
--version history
--1.0:
--Initial release. Base-functionality and GUI. Nets are not extended yet automatically. You have to run the recipe multiple times to extend nets.
--1.1:
--Smarter routine to set mutation-poses. This does not require any pre-stored large matrix anymore.
--Therefore you can now select almost an arbitrary number of unfrozen segments.
--Better ending-sequence. Even when the program is cancelled all the good found nets are loaded to undo-stack.
--Fixed a bug where the recipe would halt on puzzles which do not have HB-bonus
--Fixed a bug where wrong quicksaves were loaded upon end of the recipe
--Improved speed of duplicate-net detection. Only nets with the same order (bonds count) are compared in-depth.
--Fixed numerous bugs and did some further improvements.
--1.2:
--Added settings for behavior.SetSidechainHBondImportance(3) and behavior.SetBackboneHBondImportance(3) since we want to create as many bonds as possible and loose as little as possible.
--Added setting AAs to test via GUI
--Added option to keep only nets with added HBbonus
--Added option for linker puzzles: Sidechains on locked segments are mutable.
--Implemented "Auto-grow" modus. Found nets are stored and reloaded.
--Auto-grow: found nets are extended in multiple turns until stall. Like this large nets can quickly be found which often have good HBbonus.
--Fixed a bug in the setting of poses (data-format-limit). There is now no limit on the number of sidechains anymore. For large number of unfrozen sidechain direct mutation is used in AG with no stored poses.
--Added settings to limit the time to work on a single net in AG-mode (give-up-criterions: maximum number of poses to try or/and maximum number of iterations to try)
--Added option to set random center-segment for AG. Set "0" in center segment list!
--Added auto-grow chaos mode. All found bonds are kept. The found nets will in most cases be far less nice than without this option but maybe it can help in some instances.
--The HBfilter is now always active if it is available.
--Fixed a serious bug in single-mode where only few or no new nets were found when there were already bonds between the sidechains to mutate. Now uses netDiff instead of simple bond-count check
--Added verbosity and timestamps for a bit more detailed status-output.
--Modified distance-check to use bands in every case because bands need anyway to be used for symmetric puzzles.
--Exchanged "structure.LocalWiggleSelected()" with "structure.WiggleSelected()" because the Local function did not behave properly with the sidechain/backbone only option.
--Now in every case the initial solution before any modification is stored in the first quickstart-slot. This is necessary for AG and also helpful to identify the start when flipping through the results.
--Added mode to extend a given net with additional auto-grow. The given net must be specified in in the GUI in "AG init base seg(s)".
--Heavy bugfixing and performance improvement.
--30.04.2021:
--Fixed a small bug where the last iteration was not fully completed in AG-mode when the iteration-limit was reached.
--Added a text-output when all poses are tried and the recipe ends.
--For AG-mode the initial, unused AAs are restored before the solution is quicksaved. This is helpful if you want to directly use the result to go on.
--Mutable Cysteines are now detected as segments to work on. Since they have no sidechain, they cannot be flagged to NOT being used. If you don't want to use them, then mutate them to something else and freeze the sidechain before running this recipe.
--09.05.2021:
--Added second dialog for more AG options:
--It can now be chosen if already existing bonds of a secondary net (which gets connected via a bond to the main net) should get incorporated into the net or not. This was not done previously and it lead to not incorporating bonds that never change (e.g. on binder)
--Special mode for ligands etc: the center never changes and connections are only accepted to it. Like this you can quickly bind (as much as possible) to a ligand when you specify it as center. But the net will not grow further beyound segments that are connected to the center directly.
--It can now be chosen if for AG the unused segments (the ones not connecting to the net) should be reverted to their inital AA before storing. This gives much cleaner/nicer looking saves/results.
--Added option to input the segments to work on as space-separated text-list instead of unfrozen elements.
--1.3:
--06.09.2021:
--Prepared for public release
--Implemented the possibility to specify CI for shake and wiggle. Initial CI before running this recipe is restored when the recipe is done. This is also done for Sidechain- and BackboneImportance.
--28.04.2022:
--Added the option to leave all filters completely untouched. Like this, you can choose for yourself which filters you want to have active while running this recipe by en- or disabling them before executing this recipe.
--Variables and Constants
--************* Change if desired *******************
local qs_idx_start = 3 --startindex for quick-saves. Is overwritten from value in GUI
--*********** End: Change if desired ****************
--Internal global variables
local version = "1.3"
local dia = {} --Global GUI
local dia2 = {}
local extNets = {} --Stores all found nets
local qs_idx = qs_idx_start
local posesDone = {}
local guardInt01 = 1 --in secs
local turnFiltersOn = false --Can be used to keep filters off. May make thinks more stable?
local grow_iters = 1 --needs quasi-global access!
local CI_initial = 1 --CI before running this recipe. This is overwritten at start of the recipe
local CI_SC_initial = 1
local CI_BB_initial = 1
--Function definitions
function wait(seconds)
local start = os.time()
repeat until os.time() > start + seconds
end
function mysplit(inputstr, sep)
if sep == nil then
sep = "%s"
end
local t={}
for str in string.gmatch(inputstr, "([^"..sep.."]+)") do
table.insert(t, str)
end
return t
end
function round(inNum)
return inNum-inNum%0.001
end
function getScore()
return current.GetEnergyScore() --This function is newer and should be used instead of GetScore()
end
function checkFilters()
--Check the objectives of interest and store for later reference. Currently this is only HBbonus.
local currFilters = {}
local filterNames = filter.GetNames()
local filterBonus_tmp = nil
if (filterNames[1] ~= nil) then
for m=1,#filterNames do
local filterName = filterNames[m]
if (filter.HasBonus(filterName)) then
--Store here filters that may not degrade in score. These will be protected with the save/restore-scheme.
if (string.find(string.lower(filterName), "hydrogen") ~= nil and string.find(string.lower(filterName), "network") ~= nil) then
filterBonus_tmp = filter.GetBonus(filterName)
table.insert(currFilters, {filterBonus_tmp, filterName})
end
end
end
end
return currFilters
end
function checkMutables(cm_toWorkOn)
local outArr = {}
for m=1,#cm_toWorkOn do
outArr[m] = {}
end
for m=1,#cm_toWorkOn do
if (structure.IsLocked(cm_toWorkOn[m]) == true) then --Strange: On binder-puzzles the sidechains on the target cannot mutate but CanMutate states that it can?!
if (dia.cbx00.value == false) then --This is not very nice. The dialog is accessed from within this function. Clean it up someday.
outArr[m][1] = structure.GetAminoAcid(cm_toWorkOn[m])
else
--This means that non-locked sidechains on locked segments are mutable. This is the case for linker-puzzles.
for n=1,#aasAll do
if (structure.CanMutate(cm_toWorkOn[m], aasAll[n]) == true) then
table.insert(outArr[m], aasAll[n])
else
end
end
end
else
for n=1,#aasAll do
if (structure.CanMutate(cm_toWorkOn[m], aasAll[n]) == true) then
table.insert(outArr[m], aasAll[n])
end
end
end
end
return outArr
end
function getDistances(gd_testArr, gd_center)
--Determines reliably the minimal distances between the sidechain of the center-element and ALL others
--Distance refers to the beta-carbon. Only bands can be used in symmetry puzzles for this since structure.GetDistance() cannot be used for symmetric copies.
--Therefore I use bands in every case.
--For symmetry puzzles, only the minimum distance from center to any of the copies is output. So only one value for each residue in gd_testArr.
local gd_symCount = structure.GetSymCount() --0...
--local tmpBdIdxs = {}
--local tmpBdsIdxs = {}
local tmpDists = {}
--local tmpDist2 = {}
local outArr = {}
--local tmpBdIdx_last = {}
for m=1,#gd_testArr do
if (gd_testArr[m] == gd_center) then
outArr[m] = -1 --Just set an illegal value if center is in gd_testArr so that it will be ignored.
else
local tmpBdIdxs = {}
local tmpDists = {}
for n=1,gd_symCount+1 do
local tmpBdIdx = band.AddBetweenSegments(gd_center, gd_testArr[m], 5, 5, n-1) --Distance-check to beta-carbon. Should work also for Glycine which has no beta-carbon.
table.insert(tmpBdIdxs, tmpBdIdx)
local tmpDist = band.GetLength(tmpBdIdx)
table.insert(tmpDists, tmpDist)
end
--Now all bands from center to the symmetric copies of the current segment are drawn.
--Determine the shortest one and output this as distance.
local tmpMinDist = tmpDists[1]
local tmpMinDistIdx = 1
for n=2,#tmpDists do
if (tmpDists[n] < tmpMinDist) then
tmpMinDist = tmpDists[n]
tmpMinDistIdx = n
end
end
--local tmpString = ""
--for n=1,#tmpDists do
-- tmpString = tmpString .. round(tmpDists[n]) .. " "
--end
--wait(5) --for debug
--Delete all bands. Somehow this needs to be done in reverse since the indices need to be consecutive?
for n=#tmpDists,1,-1 do
band.Delete(tmpBdIdxs[n])
end
--Debug-check:
--local tmpBdIdx = band.AddBetweenSegments(gd_center, gd_testArr[m], 5, 5, tmpMinDistIdx-1) --Distance-check to beta-carbon. Attention! Glycine will be ignored here!
--print("Seg " .. gd_testArr[m] .. " dists: " .. tmpString .. " taken: " .. round(band.GetLength(tmpBdIdx)))
--wait(5) --for debug
--band.Delete(tmpBdIdx)
outArr[m] = tmpDists[tmpMinDistIdx]
end
end
return outArr
end
function rand()
return math.random() --0...1 (boarder-values are INCLUSIVE!)
end
function netDiff(nd_old, nd_new)
--Compares 2 results of structure.GetHBonds() for differences.
--ATTENTION! Only the TYPE of net is detected (same atom-positions).
--But basically there are different types of mutations possible.
local oldDoneVec = {}
local newDoneVec = {}
for m=1,#nd_old do
oldDoneVec[m] = 0 --Done indicator
end
for m=1,#nd_new do
newDoneVec[m] = 0
end
--Go through the nets and see where there is a difference.
local tmp_nd_old = {}
local tmp_nd_new = {}
for m=1,#nd_old do
tmp_nd_old = nd_old[m] --prestore this vector for speed.
for n=1,#nd_new do
--Check if HBond of nd_old is present also in nd_new.
--This is the case if res1, res2, atom1 and atom2 are the same!
--The subsequent "interlaced" scheme is WAY faster than testing each condition in a single if!!!
tmp_nd_new = nd_new[n] --prestore this vector for speed.
if (newDoneVec[n] == 0) then
if (tmp_nd_old["res1"] == tmp_nd_new["res1"]) then
if (tmp_nd_old["res2"] == tmp_nd_new["res2"]) then
if (tmp_nd_old["atom1"] == tmp_nd_new["atom1"]) then
if (tmp_nd_old["atom2"] == tmp_nd_new["atom2"]) then
oldDoneVec[m] = 1 --flags that this item was done/found.
newDoneVec[n] = 1 --flags that this item was done/found.
break
end
end
end
end
end
end
end
--Go through the DoneVecs and build return-Array
local outArr = {}
outArr[1] = {} --This would only be needed if we want to detect if a bond out of the old net was broken. This is not used at the moment.
outArr[2] = {}
for m=1,#oldDoneVec do
if (oldDoneVec[m] == 0) then
--outArr[1][#(outArr[1])+1] = nd_old[m] --We lost this bond in the old net!
table.insert(outArr[1], nd_old[m])
end
end
for m=1,#newDoneVec do
if (newDoneVec[m] == 0) then
--outArr[2][#(outArr[2])+1] = nd_new[m] --We gained this bond in the new net!
table.insert(outArr[2], nd_new[m])
end
end
--print ("netDiff size: " .. " " .. #(outArr[1]) .. " " .. #(outArr[2]))
return outArr
end
function filtDisAllExcl(fdae_excludedNames)
--Disables all filters excluding the ones with names in table fdae_excludedNames.
--Useful if for example HBnet bonus shall be always on.
if (#fdae_excludedNames == 0) then
filter.DisableAll()
else
local filterNames = filter.GetNames()
local found = false
for m=1,#filterNames do
found = false
for n=1,#fdae_excludedNames do
if (filterNames[m] == fdae_excludedNames[n]) then
found = true
break
end
end
if (found == false) then
filter.Disable(filterNames[m])
end
end
end
end
function getTimestamp()
return os.date("%a %X") --Abbreviated weekday (e.g. Wed) and time in format 23:12:01
end
function cleanup(err)
if (string.find(string.lower(err), "cancelled") ~= nil) then
print("User Cancelled. Cleaning up...")
--Debug: print posesDone. There may be NO jumps in it!
local debug = 0
if (debug == 1) then
print("posesDone: ")
for m=1,#posesDone do
print(posesDone[m])
end
end
print("Resetting importances and filters...")
if (dia2.cbx04.value == false) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
behavior.SetSidechainHBondImportance(CI_SC_initial)
behavior.SetBackboneHBondImportance(CI_BB_initial)
behavior.SetClashImportance(CI_initial)
if (grow_iters >= 1) then
print("Loading all up to now found Nets...")
undo.SetUndo(true)
for n=1,qs_idx-qs_idx_start do
save.Quickload(qs_idx_start+n-1) --Puts found nets into restore. Use ctrl+z/y to flip through.
if (dia2.cbx04.value == false) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
print("Net in qs" .. qs_idx_start+n-1 .. " (" .. n .. "/" .. qs_idx-qs_idx_start .. ") loaded.")
end
end
else
print("ERROR: " .. err)
end
print("iwdn HBNW " .. version .. " has finished!")
end
--***** MAIN MAIN MAIN *****
function main()
print("Running iwdn HBNW " .. version .. "...")
--Store initial CIs before running this recipe
CI_initial = behavior.GetClashImportance()
CI_SC_initial = behavior.GetSidechainHBondImportance()
CI_BB_initial = behavior.GetBackboneHBondImportance()
--Show GUI
dia = dialog.CreateDialog("HBNW " .. version)
dia.lbl01 = dialog.AddLabel("Settings for global operations:")
dia.txb02 = dialog.AddTextbox("AAs to try", "wdetsyhqn") --Def.: wdetsyhqn. Only promising ones: hdeqn
dia.cbx00a = dialog.AddCheckbox("Use text-input instead of unfrozen segs:", false)
dia.txb01 = dialog.AddTextbox("ToWorkOn segs", "0")
dia.cbx00 = dialog.AddCheckbox("SCs on locked segs are mutable", false)
dia.sli01 = dialog.AddSlider("Shake Iters", 1, 0, 20, 0)
dia.sli02 = dialog.AddSlider("Wiggle Iters", 2, 0, 20, 0)
dia.cbx00b = dialog.AddCheckbox("Wiggle backbone", false)
dia.sli04 = dialog.AddSlider("Max Nets to find", 15, 1, 50, 0)
dia.sli04b = dialog.AddSlider("Qs startslot", 3, 1, 50, 0)
dia.cbx00c = dialog.AddCheckbox("Turn filters on", true)
dia.sli04c = dialog.AddSlider("Verbosity", 0, 0, 2, 0)
dia.lbl01b = dialog.AddLabel("")
dia.lbl02 = dialog.AddLabel("Settings for Single:")
dia.sli03 = dialog.AddSlider("Min Order", 1, 1, 10, 0)
dia.cbx01 = dialog.AddCheckbox("Store only HBbonus nets", false)
dia.lbl02b = dialog.AddLabel("")
dia.lbl03 = dialog.AddLabel("Settings for Auto-Grow (AG):")
dia.cbx02 = dialog.AddCheckbox("Auto-grow nets", true)
dia.sli05 = dialog.AddSlider("AG radius (segs)", 30, 1, 100, 0)
dia.sli05b = dialog.AddSlider("AG max iters", 15, 2, 50, 0) --can be used to limit the max. complexity of the found net and to save time.
dia.sli06 = dialog.AddSlider("AG maxPoses", 20, 5, 200, 0)
dia.txb03 = dialog.AddTextbox("AG center seg(s)", "0")
dia.txb03a = dialog.AddTextbox("AG init base seg(s)", "0")
dia.cbx04 = dialog.AddCheckbox("AG store additional HB bonus nets", true)
dia.but01 = dialog.AddButton("OK", 1) --In Foldit-lua a pressed button will always hide the current window!
dia.but02 = dialog.AddButton("Cancel", 0)
dia.but03 = dialog.AddButton("More", 2)
--Secondary window for further AG-options:
dia2 = dialog.CreateDialog("HBNW " .. version .. " further options")
dia2.cbx00 = dialog.AddCheckbox("AG chaos mode", false)
dia2.cbx01 = dialog.AddCheckbox("AG revert unused segments", true)
dia2.cbx02 = dialog.AddCheckbox("AG add already connected segs", true)
dia2.cbx03 = dialog.AddCheckbox("AG connect only to fixed center", false)
dia2.sli07 = dialog.AddSlider("Shake CI", CI_initial, 0, 1, 2)
dia2.sli08 = dialog.AddSlider("Wiggle CI", CI_initial, 0, 1, 2)
dia2.cbx04 = dialog.AddCheckbox("Leave all filters untouched", false)
dia2.but01 = dialog.AddButton("OK", 0)
--Window management
local tmpBack = dialog.Show(dia)
while true do
local tmpBack2 = {}
if (tmpBack == 2) then --Secondary options dialog
tmpBack2 = dialog.Show(dia2)
tmpBack = dialog.Show(dia)
end
if (tmpBack == 0) then
print("iwdn HBNW " .. version .. " has finished.")
return
elseif (tmpBack == 1) then
break --continue with normal operation.
end
end
--set verbosity from GUI
local verbosity = dia.sli04c.value --0: minimal output
--Handling of filters
turnFiltersOn = dia.cbx00c.value --If true then only for score-readout all filters are turned on. This will get more accurate score because filter bonuses are included then.
--Set random seed
math.randomseed(os.time()-1000*os.clock()) --Disable this line if you want to have always the same result/sequence for each run.
--Set AAs from GUI
aasAll = {}
local aasAllString = "ivlfcmagtswyphdeqnkr"
local aasAccString = "" --stores the accepted identifiers
for m=1,string.len(dia.txb02.value) do
if (string.find(aasAllString, string.lower(string.sub(dia.txb02.value, m, m))) ~= nil) then --Ignore erroneous strings
aasAll[m] = string.sub(dia.txb02.value, m, m)
--print(aasAll[m])
aasAccString = aasAccString .. aasAll[m]
end
end
--Convert initial centers to num array
local ag_tmpCenter = mysplit(dia.txb03.value, " ")
for m=1,#ag_tmpCenter do
ag_tmpCenter[m] = tonumber(ag_tmpCenter[m])
-- print("ag_tmpCenter: " .. ag_tmpCenter[m])
end
--Convert initial base-segs to num array
local ag_initialBase = mysplit(dia.txb03a.value, " ")
for m=1,#ag_initialBase do
ag_initialBase[m] = tonumber(ag_initialBase[m])
-- print("ag_initialBase: " .. ag_initialBase[m])
--Since the user may have forgotten to freeze these base-residues we simply do this so that these residues are neglected in the determination of toWorkOn
if ag_initialBase[m] ~= 0 then
freeze.Freeze(ag_initialBase[m], false, true) --backbone, sidechain
end
end
--Text-input of segments toWorkOn
local txtInp_toWorkOn = {}
if (dia.cbx00a.value == true) then
--Convert initial base-segs to num array
txtInp_toWorkOn = mysplit(dia.txb01.value, " ")
for m=1,#txtInp_toWorkOn do
txtInp_toWorkOn[m] = tonumber(txtInp_toWorkOn[m])
end
end
--Check if HBond-Bonus is available.
local tmpFilters = checkFilters()
local filterHBname = {}
local hasHBbonus = false
local filt2exclArr = {}
local HBbonus_start = 0
if (#tmpFilters > 0) then
if (#(tmpFilters[1]) > 1) then
hasHBbonus = true
HBbonus_start = tmpFilters[1][1]
filterHBname = tmpFilters[1][2]
print("Puzzle has H-Bond Bonus. Current HB bonus: " .. round(HBbonus_start) .. " Pts.")
filt2exclArr[1] = filterHBname
end
end
--Check for symmetric copies
local symCount = structure.GetSymCount()
local residues = structure.GetCount()
print("Puzzle has " .. residues .. " residues and " .. symCount .. " symmetric copies.")
--List all settings from the GUI here at the start.
print("")
print("Selected global settings:")
print("Mutations to try: " .. aasAccString)
print("Use text-input instead of unfrozen segs:" .. tostring(dia.cbx00a.value))
print("ToWorkOn segments (space-delimited):" .. dia.txb01.value)
print("Sidechains on locked segments are mutable (linker/interface): " .. tostring(dia.cbx00.value))
print("Shake iterations: " .. dia.sli01.value)
print("Wiggle iters: " .. dia.sli02.value)
print("Wiggle backbone(s) in addition to sidechain(s): " .. tostring(dia.cbx00b.value))
print("Maximum number of nets to find: " .. dia.sli04.value)
print("Quicksave slot to start from: " .. dia.sli04b.value)
print("Turn filters on (more exact score-readout. Leaving it off may lead to better stability (?) and a bit of time-saving.): " .. tostring(dia.cbx00c.value))
print("Verbosity-level for output information: " .. dia.sli04c.value)
print("")
if (dia.cbx02.value == false) then --Non auto-grow
print("Selected mode: Single. Every pose with more new bonds than specified in 'Min order' will be stored as new Net.")
print("Nets to keep must have this minimum order: " .. dia.sli03.value)
print("Store only HB-Bonus-nets: " .. tostring(dia.cbx01.value))
else
--Auto-Grow
print("Selected mode: Auto-Grow (AG). Net will be grown in multiple iters around center.")
print("AG radius (segments): " .. dia.sli05.value)
print("AG maximum interations (give-up criterion): " .. dia.sli05b.value)
print("AG max poses per iteration (give-up criterion): " .. dia.sli06.value)
print("AG initial center segments to cycle (set 0 for random or e.g. 11 31 45 for specific starters.): " .. dia.txb03.value)
print("AG initial base segments (use these to try to extend to an already existing net): " .. dia.txb03a.value)
print("AG store additional HB-Bonus-nets (requires more qs-slots than 'max number of nets'): " .. tostring(dia.cbx04.value))
--Second dialog
print("AG chaos mode (keep all bonds even if they are not connected to center): " .. tostring(dia2.cbx00.value))
print("AG revert unused segments (which do not connect to the net) to original AA: " .. tostring(dia2.cbx01.value))
print("AG add already connected segs to the net: " .. tostring(dia2.cbx02.value))
print("AG connect only to fixed center: " .. tostring(dia2.cbx03.value))
print("Shake CI: " .. tostring(dia2.sli07.value))
print("Wiggle CI: " .. tostring(dia2.sli08.value))
print("Leave all filters untouched: " .. tostring(dia2.cbx04.value))
end
print("")
--Find unfrozen sidechains. User has to freeze all other sidechains beforehand!
local toWorkOn = {}
local toWorkOnString = ""
if (dia.cbx00a.value == false) then --Normal input: detection of unfrozen segments
for m=1,residues do
--IsLocked cannot really be used because this will be true for EVERY locked backbone.
--Therefore a segment will be detected as locked even when it's sidechain is mutable.
--This can be circumvented by checking available rotamers. Locked segments have only one. But cysteine and alanine and maybe others under certain circumstances have only one as well!
local frzback, frzside = freeze.IsFrozen(m)
if (frzside == false) then
--We only need to consider unfrozen sidechains. All others are irrelevant in this case.
if (structure.IsMutable(m) == true) then
--We need to consider this segment in every case because it can mutate to different sidechains. So it can potentially be part of a net when it gets mutated and bonded to.
--This will now include freely mutable cysteines as well even though they have no sidechain!
table.insert(toWorkOn, m)
toWorkOnString = toWorkOnString .. m .. " "
else
--In this case the sidechain is non-frozen and non-mutable. So it can be a completely locked segment. This needs to be detected.
if (rotamer.GetCount(m) > 1) then
--This means we have a non-locked sidechain. The bb may or may not be locked - we don't have to care about it. This segment we also want to use because it can shake and wiggle and therefore potentially connect to the net (if it has open atoms).
table.insert(toWorkOn, m)
toWorkOnString = toWorkOnString .. m .. " "
else
--This means we have a non-frozen, non-mutable AA with single rotamer.
--This can be a fully locked segment (bb and sc).
--But this can also be a cysteine, an alanine or another AA with only one sidechain.
if (structure.IsLocked(m) == false) then
--If it is also non-locked, then this could potentially be a usable sidechain. But only if it has at least one open atom
--Therefore check for open atom(s) here.
local tmpAA = structure.GetAminoAcid(m)
if (string.find("wdetsyhqnkr", string.lower(tmpAA)) ~= nil) then
table.insert(toWorkOn, m)
toWorkOnString = toWorkOnString .. m .. " "
end
end
end
end
end
end
print(#toWorkOn .. " unfrozen segments toWorkOn: " .. toWorkOnString)
else
--Text-input of segments toWorkOn
--Convert initial base-segs to num array
local warn = false
for m=1,#txtInp_toWorkOn do
if (txtInp_toWorkOn[m] >= 1 and txtInp_toWorkOn[m] <= residues and math.floor(txtInp_toWorkOn[m]) == txtInp_toWorkOn[m]) then --includes check for integer
table.insert(toWorkOn, txtInp_toWorkOn[m])
toWorkOnString = toWorkOnString .. txtInp_toWorkOn[m] .. " "
else
warn = true
end
end
if (warn == true) then
print("Warning: Please check text-input for segments in toWorkOn because at least one index could not be read!")
end
print(#toWorkOn .. " text-input segments in toWorkOn: " .. toWorkOnString)
if (#toWorkOn == 0) then
print("Error: no segments toWorkOn were found! This recipe is terminating.")
cleanup("No segments in toWorkOn found Error")
end
--freeze ALL sidechains not specified here (this may not be necessary?!)
local tmpFound = {}
for m=1,residues do
tmpFound = false
for n=1,#toWorkOn do
if (m == toWorkOn[n]) then
tmpFound = true
break
end
end
if (tmpFound == false) then
freeze.Freeze(m, false, true) --backbone, sidechain
end
end
end
--tmpBase0 = getDistances(toWorkOn, 60) --Attention! The result of getDistances() is unsorted.
--wait(1000)
--Determine current AAs of toWorkOn in order to restore the ones where no net was found before storing it to file.
--This may be helpful because you may not have to run mutate after finishing this recipe. Like this the result will be as close as possible to the initial configuration.
--There may still be differences because it is not possible to EXACTLY reset the initial structure unfortunately.
local toWorkOnAAs = {}
for m=1,#toWorkOn do
toWorkOnAAs[m] = structure.GetAminoAcid(toWorkOn[m])
end
--Auto-grow(AG): the subsequent steps need to be run in a loop.
local ag_looper = 1 --Is used for signalling that for autogrow multiple nets should be found and stored.
if (dia.cbx02.value == true) then
grow_iters = dia.sli05b.value --Set this to a very high value. Basically this should be infinite. But >10 should be ok.
ag_looper = dia.sli04.value --Set this here to max number of nets to find.
end
--Store initial configuration for restore in AG-loops
qs_idx_start = dia.sli04b.value
save.Quicksave(qs_idx_start)
qs_idx = qs_idx_start+1
print("Storing initial in qs" .. qs_idx_start .. " for later reference.")
--Disable all filters. Otherwise this is too slow.
behavior.SetSidechainHBondImportance(3) --Set this as high as possible since we want to find sidechain-bonds at any cost.
behavior.SetBackboneHBondImportance(3) --We don't want to loose backbone-H-Bonds
if (dia2.cbx04.value == false) then
filtDisAllExcl(filt2exclArr) --if filterHBname is empty then it will not lead to an error.
end
for aa=1,ag_looper do
if (aa > 1) then --Means we are in auto-grow
--reload initial net because we always want to start from the same start-condition.
--otherwise there would be unwanted clashes etc right from the start.
save.Quickload(qs_idx_start)
end
local ag_net_breaker = false --used to flag that a net shall be finished.
local foundNetRes = {} --included residues in the found net as indices.
local HBbonusArr = {} --stores the diverse HB bonus values found in all iterations for AG.
local tmpCenter0 = {}
--local ag_tmpCenterMutate = true
for a=1,grow_iters do
--Out of all segments to work on select a region in which to look for HB-nets
local tmpCenter = {}
local tmpBase = {} --All residues in this array are mutated randomly in the current iteration. All others are not.
local tmpBase0 = {} --Temporarily required for distance-check.
if (a == 1) then
if (grow_iters == 1) then --No auto-grow
--Simply choose ALL unfrozen segments. No distance-check is required.
tmpBase = toWorkOn
print("Non-auto-grow mode (Single) started...")
else
--Auto-grow is chosen. Select a center and find closest residues.
if (ag_tmpCenter[1] ~= 0) then
tmpCenter = ag_tmpCenter[(aa-1) % #ag_tmpCenter + 1] --Get residue-number of center from GUI. Needs to be converted from string to number!
else
--This means we want to start with a random segment out of all toWorkOn
local tmpVal = rand() * #toWorkOn
if (tmpVal == #toWorkOn) then
tmpVal = #toWorkOn-1 --Special treatment for rand()==1
end
tmpCenter = toWorkOn[math.floor(tmpVal)+1]
end
if (ag_initialBase[1] == 0) then
tmpBase[1] = tmpCenter --In the first turn the center needs to be in the base. Because we want to mutate it as well to get a HBond.
else
--Special mode: connect to a given net.
--check if current tmpCenter is in ag_initialBase. If this is the case then it may not be mutated and therefore not be put into tmpBase.
local tmpFound = false
for m=1,#ag_initialBase do
if (tmpCenter == ag_initialBase[m]) then --In this case the tmpCenter shall NOT be mutated because it is part of the net already.
tmpFound = true
break
end
end
if (tmpFound == false) then
--tmpCenter is not in initBase. It may be mutated.
tmpBase[1] = tmpCenter
end
end
tmpCenter0 = tmpCenter --Needed for summary in AG
--Find closest residues and put them in base.
tmpBase0 = getDistances(toWorkOn, tmpCenter) --Attention! The result of getDistances() is unsorted.
local tmpBase_start = #tmpBase
for m=1,dia.sli05.value-tmpBase_start do --There may already be segments assigned to tmpBase. Therefore -tmpBase_start here.
--Go through the distances and find smallest ones
local tmpMinDistIdx = {}
local tmpMinDist = -2
for n=1,#tmpBase0 do
if (tmpBase0[n] > 0) then --Protection against invalid distances
if (tmpMinDist == -2) then --Initial
tmpMinDist = tmpBase0[n]
tmpMinDistIdx = n
else
if (tmpBase0[n] < tmpMinDist) then
tmpMinDist = tmpBase0[n]
tmpMinDistIdx = n
end
end
end
end
--It can happen that more residues than there are available shall be worked on.
--Avoid this here by breaking this loop.
if (tmpMinDist == -2) then --no more distance-segment was found!
break
end
--Store found minimum
table.insert(tmpBase, toWorkOn[tmpMinDistIdx]) --The smallest distance corresponds to a segment stored in toWorkOn!
tmpBase0[tmpMinDistIdx] = -1 --set to a useless value so that this value is not considered anymore in the next turn.
end
print("Auto-growing Net " .. aa .. "/" .. dia.sli04.value .. " center: " .. tmpCenter .. ". Please wait...")
if (verbosity >= 1) then
print(getTimestamp() .. ": Auto-grow iter 1 started on " .. #tmpBase .. " segments...")
end
if (verbosity >= 2) then
local tmpBase_string = ""
for m=1,#tmpBase do
tmpBase_string = tmpBase_string .. tmpBase[m] .. " "
end
tmpBase_string = string.sub(tmpBase_string, 1, #tmpBase_string-1)
print("Center (init/curr): " .. tmpCenter0 .. "/" .. tmpCenter .. " tmpBase: " .. tmpBase_string)
end
end
else
--For further loops the center needs to be chosen out of the already found net. Prefereably the last found segment?!
--TODO: Avoid taking a center residue twice.
--Check the residues of the net for a new center for this turn. Simply pick the last entry. This may not be smart.
--It should be detected which residue is actually the open end. But that will take more doing.
if (dia2.cbx03.value == true) then --In this case/mode the center shall always stay the same! This is always the first element in foundNetRes!
tmpCenter = foundNetRes[1]
else
tmpCenter = foundNetRes[#foundNetRes] --Standard mode: the center will vary
end
--Find closest residues and put them in base.
local toWorkOn_tmp = {} --subtract the already found net and store it in this variable.
local foundEntry = 0
for m=1,#toWorkOn do
foundEntry = 0
for n=1,#foundNetRes do
if (toWorkOn[m] == foundNetRes[n]) then
foundEntry = 1
end
end
if (foundEntry == 0) then
table.insert(toWorkOn_tmp, toWorkOn[m])
end
end
foundEntry = nil --invalidate
local debug = 0
if (debug == 1) then
local toWorkOn_tmpString = ""
for m=1,#toWorkOn_tmp do
toWorkOn_tmpString = toWorkOn_tmpString .. toWorkOn_tmp[m] .. " "
end
--print("Unfrozen segments toWorkOn: " .. toWorkOn_tmpString)
end
tmpBase0 = getDistances(toWorkOn_tmp, tmpCenter) --Attention! The result of getDistances() is unsorted.
for m=1,dia.sli05.value do --There cannot be segments assigned to tmpBase already in this case!
--Go through the distances and find smallest ones
local tmpMinDistIdx = {}
local tmpMinDist = -2
for n=1,#tmpBase0 do
if (tmpBase0[n] > 0) then --Protection against invalid distances
if (tmpMinDist == -2) then --Initial
tmpMinDist = tmpBase0[n]
tmpMinDistIdx = n
else
if (tmpBase0[n] < tmpMinDist) then
tmpMinDist = tmpBase0[n]
tmpMinDistIdx = n
end
end
end
end
--It can happen that more residues than there are available shall be worked on.
--Avoid this here by breaking this loop.
if (tmpMinDist == -2) then --no more distance-segment was found!
break
end
--Store found minimum
table.insert(tmpBase, toWorkOn_tmp[tmpMinDistIdx]) --The smallest distance corresponds to a segment stored in toWorkOn_tmp!
tmpBase0[tmpMinDistIdx] = -1 --set to a useless value so that this value is not considered anymore in the next turn.
end
if (verbosity >= 1) then
print(getTimestamp() .. ": Auto-grow iter " .. a .. " started on " .. #tmpBase .. " segments...")
end
if (verbosity >= 2) then
local tmpBase_string = ""
for m=1,#tmpBase do
tmpBase_string = tmpBase_string .. tmpBase[m] .. " "
end
tmpBase_string = string.sub(tmpBase_string, 1, #tmpBase_string-1)
print("Center (init/curr): " .. tmpCenter0 .. "/" .. tmpCenter .. " tmpBase: " .. tmpBase_string)
local tmpFoundNetRes_string = ""
for m=1,#foundNetRes do
tmpFoundNetRes_string = tmpFoundNetRes_string .. foundNetRes[m] .. " "
end
tmpFoundNetRes_string = string.sub(tmpFoundNetRes_string, 1, #tmpFoundNetRes_string-1)
print("Found net so far (" .. #foundNetRes .. " segs): " .. tmpFoundNetRes_string)
end
end --End: if (a == 1) then
--Now the base-residues are determined. Find the mutation-pose for them.
--For auto-grow it is not necessary to checkMutables again. Do this only in the first turn for all residues in toWorkOn.
--This saves some time.
local tmpBase_mutables = {}
if (dia.cbx02.value == false) then --non-auto-grow
tmpBase_mutables = checkMutables(tmpBase) --tmpBase is ALWAYS toWorkOn in this case!
else
--For auto-grow use this...
if (aa == 1 and a == 1) then
toWorkOn_mutables = checkMutables(toWorkOn) --find mutables for ALL possible residues to be used in the future.
end
--Go through toWorkOn_mutables and find the corresponding mutables for the current tmpBase.
for m=1,#tmpBase do
--tmpBase_mutables[m] = {} --init 2-dimensional array. Is this necessary?
for n=1,#toWorkOn do
if (tmpBase[m] == toWorkOn[n]) then
table.insert(tmpBase_mutables, toWorkOn_mutables[n])
break
end
end
end
end
--calculate the maximum number of poses and display it
local maxPoses = 1 --Initialize to 1 because we will do multiplication
--local mutables_string = "" --For debug only
for m=1,#tmpBase_mutables do
maxPoses = maxPoses * #(tmpBase_mutables[m]) --calculate max number of combinations
--mutables_string = mutables_string .. #(tmpBase_mutables[m]) .. " "
end
--print("maxPoses= " .. maxPoses)
--print("mutables_string: " .. mutables_string)
--Determine the number of HBonds BEFORE the search begins
local startBonds = structure.GetHBonds()
selection.DeselectAll()
for m=1,#tmpBase do
selection.Select(tmpBase[m]) --select residues in the vicinity
end
local extNetsMinOrder = dia.sli03.value
local extNetsBondsNr = {} --Stores the number of bonds for each found net in extNets.
--Basic approach: Don't prestore anything. Just store the poses that were done and randomly try new ones.
posesDone = {} --Stores the already done poses. This holds just a decimal value for each pose. Sorted with rising values.
local ag_breaker = false --necessary for auto-grow-mode to end a single grow-iteration.
for m=1,maxPoses do
--In each turn of m, another pose is tried with random indexing.
--Also it is checked if some new network has formed.
local tmpRand = math.floor(rand() * maxPoses)+1 --Random integer index from 1...maxPoses
if (tmpRand >= maxPoses) then --clipping to avoid value > maxPoses
tmpRand = maxPoses
end
--print("maxPoses: " .. maxPoses)
--print("currPose: " .. tmpRand)
local tmpComb = {}
if (maxPoses < 10^10 and dia.cbx02.value == true) then --Treatment for huge number of residues to try. Since we would reach the data-format-limit. For single-mode we NEED to use the first scheme because we need to fill the posesDone array.
--The subsequent if-condition and internals are a BEAST! Gave me some headache to figure it out.
--I would be thankful if someone could hint me to a simpler solution.
if (#posesDone == 0) then
posesDone[1] = tmpRand --First initialization.
else
--Check if this random pose was already used. If not then find free slot.
for n=1,#posesDone do --Run once from start to end over posesDone-Array
--Special case: first element
if (n == 1) then
if (tmpRand < posesDone[1]) then --Insert at front. Other cases will be treated in following steps.
table.insert(posesDone, 1, tmpRand)
break
end
end
--Special case: we hit the last element. No element follows in posesDone!
if (n == #posesDone) then
if (posesDone[n] == maxPoses) then --In this case we NEED to wrap! This case is: tmpRand = posesDone[#posesDone]!
--wrap around and find free slot!
for o=1,#posesDone-1 do --minus 1 because we can omit the current (last) slot.
if (posesDone[o+1] > posesDone[o]+1) then
table.insert(posesDone, o+1, posesDone[o]+1)
tmpRand = posesDone[o]+1
break
end
end
break --Break the loop over n since we placed the item.
else --this means: last element in posesDone is < maxPoses
if (tmpRand == posesDone[n]) then
posesDone[n+1] = posesDone[n]+1 --This slot MUST be free because we are smaller than maxPoses.
tmpRand = posesDone[n]+1
break
else --tmpRand is larger than the last element. Smaller is not possible since this case is treated in the next steps.
posesDone[n+1] = tmpRand
break
end
end
end
--Now the probably most common and difficult case:
--tmpRand lies either on a value (and another value follows in posesDone) or tmpRand lies between two values.
if (tmpRand >= posesDone[n] and tmpRand < posesDone[n+1]) then
if (tmpRand > posesDone[n] and tmpRand < posesDone[n+1]) then
--we have found a free slot. Store the new element in between these values in posesDone.
table.insert(posesDone, n+1, tmpRand)
break
else
--This means that tmpRand lies ON posesDone[n]. We have to increment to find a free slot.
--This is the most difficult case because it may happen that we will have to wrap around the end of posesDone.
local wrapOffset = 0
for o=1,#posesDone-1 do --minus 1 because we can omit the current (last) slot.
local tmpIdx = n+o-1+wrapOffset --first index is posesDone[n]!
--Treat the case where tmpIdx points to the FIRST element in posesDone. This can happen due to wrapping!
if (tmpIdx == 1) then
if (posesDone[1] > 1) then --Insert at front
table.insert(posesDone, 1, posesDone[1]-1)
tmpRand = posesDone[1]-1 --Step back
break
else
--Nothing needs to be done here? This is covered by subsequent cases?
end
end
if (tmpIdx == #posesDone) then
--We have reached the end of posesDone. Either there is a free slot or we will have to wrap.
if (posesDone[tmpIdx] == maxPoses) then --In this case we NEED to wrap!
--wrap around and find free slot!
wrapOffset = -#posesDone --All subsequent actions need to be executed shifted!
--Is this action sufficient or is there something else to be done?
else
--We have found a free slot at the end of posesDone
table.insert(posesDone, posesDone[tmpIdx]+1)
tmpRand = posesDone[tmpIdx]+1
break
end
else
--If this point is reached then there is still an element following the current one in posesDone.
--We can safely use the next element.
if (posesDone[tmpIdx+1] > posesDone[tmpIdx]+1) then
--We found a free slot.
table.insert(posesDone, tmpIdx+1, posesDone[tmpIdx]+1)
tmpRand = posesDone[tmpIdx]+1
break
else
--Do nothing. Wait for next turn.
end
end
end --End: internal for-loop (o)
break --Break loop over n
end
end --End: interval-case
--There is no other case that needs to be treated. If the above conditions are not satisfied then we just have to
--do another turn for n.
end --End: for n=1,#posesDone do
end --End: if (#posesDone == 0)
--Decode the current pose and apply
--Strange: this routine cannot take more than about 15 segments. Then some strange clipping occurs. Maybe the data-format?
--Yes, seems like it. So I should use for >15 a different scheme. I should not store anything but set the residues directly random.
--This is because it is VERY unlikely that with so many residues I will get the same pose twice.
tmpRand = tmpRand-1 --0...maxPoses-1
--tmpRand = maxPoses-1-5*9^8-2*9^3 --debug
local tmpIdx = 0
local tmpVal = tmpRand
local tmpVal2 = 0
for n=1,#tmpBase_mutables do
tmpVal2 = math.floor(tmpVal/#(tmpBase_mutables[n]))
tmpIdx = tmpVal-tmpVal2*#(tmpBase_mutables[n]) --Integer-index to the current slot in tmpBase_mutables.
--print(#(tmpBase_mutables[n]) .. " " .. tmpVal2 .. " " .. tmpIdx)
tmpVal = tmpVal2 --Prepare for next turn
tmpComb[n] = tmpBase_mutables[n][tmpIdx+1]
end
else
--Huge number of residues to try in auto-grow. Workaround to avoid data-format limitation
--Simply set each residue directly and don't store anything.
for n=1,#tmpBase_mutables do
if (#(tmpBase_mutables[n]) == 1) then
tmpComb[n] = tmpBase_mutables[n][1]
else
tmpRand = math.floor(rand() * #(tmpBase_mutables[n]))
if ( tmpRand > #(tmpBase_mutables[n])-1 ) then --special treatment of "1" in rand()
tmpRand = #(tmpBase_mutables[n])-1
end
tmpComb[n] = tmpBase_mutables[n][tmpRand+1]
end
end
--For compatibility: fill posesDone. This is used for debug output
table.insert(posesDone, 0)
end
--local tmpString = ""
--for n=1,#tmpComb do
-- print(tmpComb[n])
-- tmpString = tmpString .. tmpComb[n] .. " "
--end
--print ("tmpComb: " .. tmpString)
for n=1,#tmpComb do
structure.SetAminoAcid(tmpBase[n], tmpComb[n])
end
if (dia.sli01.value > 0) then
behavior.SetClashImportance(dia2.sli07.value) --Set shake-CI according to what is chosen/set in the GUI
structure.ShakeSidechainsSelected(dia.sli01.value) --Shake really helps to pre-order things. I think it should be used.
end
if (dia.sli02.value > 0) then
--Very strange! LocalWiggleSelected does not seem to work right with the backbone and sidechain options. If backbone is false than it does not execute.
--Submit this as bug or investigate further later. Workaround: WiggleSelected. This does work with these options. Strange!
--structure.LocalWiggleSelected(dia.sli02.value, dia.cbx00b.value, true) --backbone, sidechains <== DO NOT USE THIS!!!
behavior.SetClashImportance(dia2.sli08.value) --Set wiggle-CI according to what is chosen/set in the GUI
structure.WiggleSelected(dia.sli02.value, dia.cbx00b.value, true)
end
--If the puzzle has HBbonus and it shall be incorporated then we can just check for this here.
--This jumper probably only makes sense in single-mode since for auto-grow in many cases you can get HBbonus by modifiying the net manually. So better store them for auto-grow anyway.
--This should save quite some time in single-mode since GetHBonds does not need to be called.
local jumper = false --lua has no continue so this workaround is necessary.
if (hasHBbonus == true and dia.cbx02.value == false) then --Only for single-mode!
if (dia.cbx01.value == true) then
local tmpFilters = checkFilters()
if (tmpFilters[1][1] <= HBbonus_start) then
jumper = true --skip the rest of this turn since no gain in HBnet-bonus was achieved.
end
end
end
if (jumper == false) then
--check if some HBond has formed
tmpBonds = structure.GetHBonds()
--In non-autogrow-mode we can NOT simply check the difference in bond-count.
--This is because we may have bonds in the initial pose before the first mutation. These bonds will rip off due to pose-mutations and we would then have to find MORE new bonds than the ones that ripped off.
--This is not desired and therefore we will have to use netDiff here instead as it is done in the auto-grow-routine.
local diffBack = netDiff(startBonds, tmpBonds) --We need this for single- AND auto-grow mode!
if (dia.cbx02.value == false) then --This means we are in non-auto-grow mode!
if (#(diffBack[2]) >= extNetsMinOrder) then
if (#extNets == 0) then
extNets[1] = tmpBonds --We need to store each new found net for later reference to compare for duplicate. This is only necessary in single-mode.
extNetsBondsNr[1] = #tmpBonds
if (hasHBbonus == false) then
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
end
if (verbosity == 0) then
print("NET " .. #extNets .. ": +" .. #(diffBack[2]) .. " bds in qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts.")
else
print(getTimestamp() .. ": NET " .. #extNets .. ": +" .. #(diffBack[2]) .. " bds in qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts.")
end
save.Quicksave(qs_idx)
qs_idx = qs_idx+1
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filtDisAllExcl(filt2exclArr)
end
end
else
local tmpFilters = checkFilters()
if (dia.cbx01.value == false or (dia.cbx01.value == true and tmpFilters[1][1] > HBbonus_start)) then
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
end
if (verbosity == 0) then
print("NET " .. #extNets .. ": +" .. #(diffBack[2]) .. " bds in qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts HB " .. round(tmpFilters[1][1]) .. " Pts.")
else
print(getTimestamp() .. ": NET " .. #extNets .. ": +" .. #(diffBack[2]) .. " bds in qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts HB " .. round(tmpFilters[1][1]) .. " Pts.")
end
save.Quicksave(qs_idx)
qs_idx = qs_idx+1
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filtDisAllExcl(filt2exclArr)
end
end
end
end
else
--Non-auto-grow (single): further turns
for n=1,#extNets do
--Only nets with the same number of bonds need to be compared more in-depth.
--This will save a large amount of time.
if (#tmpBonds == extNetsBondsNr[n]) then
local diffBack = netDiff(extNets[n], tmpBonds)
if (#(diffBack[1]) == 0 and #(diffBack[2]) == 0) then
--This means we have found a duplicate net with the same bonds as it is stored in one of extNets already. Ignore.
--print("Found duplicate network to Nr " .. n .. ". Ignoring it.")
break
end
end
if (n == #extNets) then --Only evaluate in the last turn!
table.insert(extNets, tmpBonds)
table.insert(extNetsBondsNr, #tmpBonds)
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
end
if (hasHBbonus == false) then
if (verbosity == 0) then
print("NET " .. #extNets .. ": +" .. #(diffBack[2]) .. " bds in qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts.")
else
print(getTimestamp() .. ": NET " .. #extNets .. ": +" .. #(diffBack[2]) .. " bds in qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts.")
end
save.Quicksave(qs_idx)
qs_idx = qs_idx+1
else
local tmpFilters = checkFilters()
if (dia.cbx01.value == false or (dia.cbx01.value == true and tmpFilters[1][1] > HBbonus_start)) then
if (verbosity == 0) then
print("NET " .. #extNets .. ": +" .. #(diffBack[2]) .. " bds in qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts HB " .. round(tmpFilters[1][1]) .. " Pts.")
else
print(getTimestamp() .. ": NET " .. #extNets .. ": +" .. #(diffBack[2]) .. " bds in qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts HB " .. round(tmpFilters[1][1]) .. " Pts.")
end
save.Quicksave(qs_idx)
qs_idx = qs_idx+1
end
end
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filtDisAllExcl(filt2exclArr)
end
end
end
end --End: for n
end
end
if (#extNets >= dia.sli04.value) then --stop criterion
break
end
else --Auto grow mode!!!
--The order may not simply be derived by subtracting #startBonds from #tmpBonds.
--This is because in previous stages some bonds probably have formed on the other, irrelevant sidechains.
--These will break down and form new ones in further turns. But the broken ones will lead to an error in this equation.
--This will probably lead to ignoring promising new nets!
--It would be better to use netDiff for this because that function "splits" broken and new bonds.
--Simply accept ANY order in this mode. Anyway only bonds connected to the center/base-net will be kept!
if (#(diffBack[2]) >= 1) then --This means at least one new bond has formed.
--local foundNet_tmp = diffBack[2]
--check for connections of diffBack to foundNet. Since multi-connections via multiple segments in diffBack are possible, we need an infinite loop here.
if (dia2.cbx00.value == true) then --this means we are in chaos-mode.
--Keep external nets. This should be very fast but chaotic. We will see.
--Simply store ALL found new bonds and freeze the connected sidechains
--The center is not important in this mode. ONLY for determining the base. So this does not matter here.
--We need to check for duplicate residues in foundNetRes especially for symm puzzles.
for c=1,#foundNet_tmp do
if (#foundNetRes == 0) then
--Special treatment for first element. Simply put in the first found one since we want to accept ALL bonds
foundNetRes[1] = foundNet_tmp[c]["res1"]%(residues+1)
end
local tmpFound = false
for d=1,#foundNetRes do
--Go through this array and see if we already found this residue. Store only new ones.
if (foundNet_tmp[c]["res1"]%(residues+1) == foundNetRes[d]) then
tmpFound = true
break
end
end
if (tmpFound == false) then
table.insert(foundNetRes, foundNet_tmp[c]["res1"]%(residues+1))
freeze.Freeze(foundNet_tmp[c]["res1"]%(residues+1), false, true) --backbone, sidechain
end
tmpFound = false
for d=1,#foundNetRes do
--Go through this array and see if we already found this residue. Store only new ones.
if (foundNet_tmp[c]["res2"]%(residues+1) == foundNetRes[d]) then
tmpFound = true
break
end
end
if (tmpFound == false) then
table.insert(foundNetRes, foundNet_tmp[c]["res2"]%(residues+1))
freeze.Freeze(foundNet_tmp[c]["res2"]%(residues+1), false, true) --backbone, sidechain
end
end
ag_breaker = true --this flags that the current auto-grow iteration shall be ended because a new bond has been found.
else
--This means we are in AG-non-chaos mode
--Unfortunately it is possible that some bonds may not be considered in this simple scheme. Because as it is, only new-formed bonds which are connected somehow to the net will be kept.
--But especially for binder-puzzles it can happen that there are bonds that never change and still get connected to the net via a new-formed bond. These will also have to be incorporated!
--This can be done by checking ALL current bonds for connections to the net. But we have to exclude BB to BB bonds.
if (dia2.cbx02.value == true) then
foundNet_tmp = tmpBonds --Add already connected bonds on a secondary net to the net when it gets connected via a bond to the main net.
else
foundNet_tmp = diffBack[2] --Only newly formed bonds can get added to the net. This may leave out some.
end
if (a == 1) then
--special treatment in first turn.
--If it has not been done already then the first item in foundNetRes must be set to the initial center.
--Like this the same routine can be used as in further loops.
if (#foundNetRes == 0) then --Isn't this always the case for a==1?
if ag_initialBase[1] == 0 then
foundNetRes[1] = tmpCenter
--print("Setting tmpCenter seg " .. tmpCenter .. " as first element in foundNetRes.")
else
--Special mode: connect to given net.
for d=1,#ag_initialBase do
table.insert(foundNetRes, ag_initialBase[d])
end
end
end
end
--Only new bonds that are CONNECTED to the old net may be kept!
--This is the most prominent and difficult routine for the auto-grow modus
while true do
--210331: This is strange! The symmetric copies seem to start at residues+2! I have read that this is some Foldit-internal feature connected with Rosetta?!
--I will have to investigate this further!
--This is gonna get ugly. It can happen that a sidechain bonds to the backbone. This will also be detected as net-extension.
--Avoid this in the first step by detecting that the bond goes at least to the beta-carbon. ==> Removed this to see if it is problematic or helps. Uncomment below line for re-activation.
local c_ender = false
local foundNetRes_maxIdx = {}
if (dia2.cbx03.value == true) then
foundNetRes_maxIdx = 1 --Spechial mode: the center is fixed and we may ONLY connect to IT.
else
foundNetRes_maxIdx = #foundNetRes --Normal mode: We may connect to ANY already found net segment.
end
for c=1,#foundNet_tmp do
for d=1,foundNetRes_maxIdx do
--if (foundNet_tmp[c]["res1"]%(residues+1) == foundNetRes[d] and (foundNet_tmp[c]["atom1"] >= 5 and foundNet_tmp[c]["atom2"] >= 5)) then --Means that we found an extended connection. Between sidechains only!
if (foundNet_tmp[c]["res1"]%(residues+1) == foundNetRes[d] and (foundNet_tmp[c]["atom1"] >= 5 or foundNet_tmp[c]["atom2"] >= 5)) then --Allows for connections from BB to SC and vice versa. Butt BB to BB is excluded.
--if (foundNet_tmp[c]["res1"]%(residues+1) == foundNetRes[d]) then --This accepts bonds to the backbone as well. Also BB to BB!
--Check if second residue is also in foundNetRes. This can happen for symm-puzzles and it must be avoided that they are stored twice.
local tmpFound = false
for e=1,#foundNetRes do
if (foundNetRes[e] == foundNet_tmp[c]["res2"]%(residues+1)) then
tmpFound = true
break
end
end
if (tmpFound == true) then
--This means that we have found a duplicate. Invalidate it an do NOT store it in foundNet
table.remove(foundNet_tmp, c)
--print("found duplicate1")
c_ender = true
break --break loop over foundNetRes. But also the loop over foundNet_tmp must be broken! Simply go through foundNet_tmp again in the while-loop
else
if (a == 1 and #foundNetRes == 1) then
freeze.Freeze(tmpCenter, false, true) --Freeze the initial center.
end
table.insert(foundNetRes, foundNet_tmp[c]["res2"]%(residues+1)) --The second atom is the new one. The first is already in the foundNet.
freeze.Freeze(foundNet_tmp[c]["res2"]%(residues+1), false, true) --backbone, sidechain
table.remove(foundNet_tmp, c) --invalidate
ag_breaker = true --this flags that the current auto-grow iteration shall be ended because a new bond has been found.
c_ender = true
break
end
--elseif (foundNet_tmp[c]["res2"]%(residues+1) == foundNetRes[d] and (foundNet_tmp[c]["atom1"] >= 5 and foundNet_tmp[c]["atom2"] >= 5)) then
elseif (foundNet_tmp[c]["res2"]%(residues+1) == foundNetRes[d] and (foundNet_tmp[c]["atom1"] >= 5 or foundNet_tmp[c]["atom2"] >= 5)) then
--elseif (foundNet_tmp[c]["res2"]%(residues+1) == foundNetRes[d]) then --Accepts bonds to BB as well.
--Check if first residue is also in foundNetRes. This can happen for symm-puzzles and it must be avoided that they are stored twice.
local tmpFound = false
for e=1,#foundNetRes do
if (foundNetRes[e] == foundNet_tmp[c]["res1"]%(residues+1)) then
tmpFound = true
break
end
end
if (tmpFound == true) then
--This means that we have found a duplicate. Invalidate it an do NOT store it foundNet
table.remove(foundNet_tmp, c)
--print("found duplicate2")
c_ender = true
break --break loop over foundNetRes. But also the loop over foundNet_tmp must be broken! Simply go through foundNet_tmp again in the while-loop
else
if (a == 1 and #foundNetRes == 1) then
freeze.Freeze(tmpCenter, false, true) --Freeze the initial center.
end
table.insert(foundNetRes, foundNet_tmp[c]["res1"]%(residues+1)) --The second atom is the new one. The first is already in the foundNet.
freeze.Freeze(foundNet_tmp[c]["res1"]%(residues+1), false, true) --backbone, sidechain
table.remove(foundNet_tmp, c) --invalidate
ag_breaker = true --this flags that the current auto-grow iteration shall be ended because a new bond has been found.
c_ender = true
break
end
end
end --End: d
if (c_ender == true) then
break --Break loop over c
end
end --End: c
if (c_ender == false) then
break --break infinite while loop
end
end --End: Infinite while loop
end --End: keep external nets
end --End: if (#(diffBack[2]) >= 1) then
end --End: auto-grow mode!!!
end --End: if (jumper == false) then
--print out the number of completed poses for auto-grow-loop
if (verbosity >=1) then
if (ag_looper > 1) then
if (#posesDone % 5 == 0) then
print(getTimestamp() .. ": AG posesDone in this iter: " .. #posesDone .. "/" .. dia.sli06.value)
end
end
end
if (ag_looper > 1 and #posesDone >= dia.sli06.value) then --AG Giveup-criterion: Max number of poses to try reached.
--This means that the net is not growing anymore.
--So we can store the net and start a new net in the next turn about (aa)
ag_breaker = true --break the loop about grow_iters
ag_net_breaker = true
end
if (a >= grow_iters and grow_iters > 1 and ag_breaker == true) then --AG Giveup-criterion: Max number of iters to try reached; grow_iters > 1: this is an override for single-mode. Leave it in!
ag_net_breaker = true
end
if (ag_net_breaker == true) then
--Restore of AAs
--Only do this for non-HBbonus-puzzles since residues outside of the found net may get this bonus!
--As of 210505 this can be set via GUI. The user should be aware that this should be deactivated probably on puzzles with HB-bonus!
if (dia2.cbx01.value == true and dia.cbx02.value == true) then --Only do this in AG-mode at the moment.
if (verbosity > 0) then
print("Restoring AAs of non-used segments in toWorkOn for qs...")
end
selection.DeselectAll()
--Only act on the residues which are NOT in foundNetRes but are in toWorkOn.
--The corresponding AAs to restore were stored in toWorkOnAAs.
local foundEntry = 0
for n=1,#toWorkOn do
foundEntry = 0
for o=1,#foundNetRes do
if (toWorkOn[n] == foundNetRes[o]) then
foundEntry = 1
break
end
end
if (foundEntry == 0) then
--Reset the initially saved AA for this segment
structure.SetAminoAcid(toWorkOn[n], toWorkOnAAs[n])
selection.Select(toWorkOn[n])
end
end
foundEntry = nil --invalidate
--Optional: Shake/wiggle the restored sidechains to get back to how they were placed before?!
structure.ShakeSidechainsSelected(1)
structure.WiggleSelected(1, false, true) --backbone, sidechains
selection.DeselectAll()
end
--store the net in extNets.
tmpBonds = structure.GetHBonds()
table.insert(extNets, tmpBonds)
table.insert(extNetsBondsNr, #tmpBonds)
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
end
if (hasHBbonus == false) then
print("NET " .. aa .. ": qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts center " .. tmpCenter0 .. " net " .. #foundNetRes .. " segs " .. a .. " iters.")
save.Quicksave(qs_idx)
qs_idx = qs_idx+1
else
local tmpFilters = checkFilters()
if (dia.cbx01.value == false or (dia.cbx01.value == true and tmpFilters[1][1] > HBbonus_start)) then
print("NET " .. aa .. ": qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts HB " .. round(tmpFilters[1][1]) .. " Pts center " .. tmpCenter0 .. " net " .. #foundNetRes .. " segs " .. a .. " iters.")
save.Quicksave(qs_idx)
qs_idx = qs_idx+1
end
end
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filtDisAllExcl(filt2exclArr)
end
end
else
--If the puzzle has HBbonus and it shall be incorporated then we can just check for this here.
--As it is, EVERY pose that has a HBbonus is stored as long as its bonus differs from the already found bonusses for this net. So as little as possible bonus-nets are probably lost.
if (hasHBbonus == true and ag_net_breaker == false and dia.cbx02.value == true) then --Only execute in intermediate iters and NOT at end. Only needs to be executed in auto-grow mode!
local tmpFilters = checkFilters()
if (tmpFilters[1][1] > 0 and tmpFilters[1][1] ~= HBbonus_start) then
--Search in HBbonusArr if this bonus was already achieved in previous iterations
if (#HBbonusArr == 0) then
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
end
HBbonusArr[1] = tmpFilters[1][1]
print("NET " .. aa .. "_" .. #HBbonusArr .. ": qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts HB " .. round(tmpFilters[1][1]) .. " Pts center " .. tmpCenter0 .. " net " .. #foundNetRes .. " segs " .. a .. " iters.")
save.Quicksave(qs_idx)
qs_idx = qs_idx+1
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filtDisAllExcl(filt2exclArr)
end
end
else
for o=1,#HBbonusArr do
if (tmpFilters[1][1] == HBbonusArr[o]) then
break
end
if (o == #HBbonusArr) then
--We found a new HB bonus. Store this as new quicksave!
table.insert(HBbonusArr, tmpFilters[1][1])
if (dia.cbx04.value == true) then
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
end
print("NET " .. aa .. "_" .. #HBbonusArr .. ": qs" .. qs_idx .. " Score " .. round(getScore()) .. " Pts HB " .. round(tmpFilters[1][1]) .. " Pts center " .. tmpCenter0 .. " net " .. #foundNetRes .. " segs " .. a .. " iters.")
save.Quicksave(qs_idx)
qs_idx = qs_idx+1
if (dia2.cbx04.value == false) then
if (turnFiltersOn == true) then
filtDisAllExcl(filt2exclArr)
end
end
end
end
end
end
end
end
end
if (m == maxPoses) then
print("Warning: All possible poses were tried. Provide more unfrozen sidechains to avoid this.")
end
if (ag_breaker == true) then
break --Break loop over maxPoses (m)
end
end --End: for m=1,maxPoses do
--Debug: print posesDone. There may be NO jumps in it!
local debug = 0
if (debug == 1) then
print("posesDone: ")
for m=1,#posesDone do
print(posesDone[m])
end
end
if (ag_net_breaker == true) then
break
end
end --END: for a=1,grow_iters do
if (dia.cbx02.value == true and #extNets >= dia.sli04.value) then --stop criterion
break
end
end --END: for aa=1,ag_looper do
print("Resetting importances and filters...")
behavior.SetSidechainHBondImportance(CI_SC_initial)
behavior.SetBackboneHBondImportance(CI_BB_initial)
behavior.SetClashImportance(CI_initial)
if (dia2.cbx04.value == false) then
filter.EnableAll()
end
print("Restoring saves to undo...")
undo.SetUndo(true)
for n=1,qs_idx-qs_idx_start do
save.Quickload(qs_idx_start+n-1) --Puts found nets into restore. Use ctrl+z/y to flip through.
if (dia2.cbx04.value == false) then
filter.EnableAll()
wait(guardInt01) --This seems to be necessary in order for things to settle.
end
print("Net in qs" .. qs_idx_start+n-1 .. " (" .. n .. "/" .. qs_idx-qs_idx_start .. ") loaded.")
end
print("iwdn HBNW " .. version .. " has finished!")
end --main
xpcall(main, cleanup) --For debug: disable this line and instead enable the line below!
--main()
