Code
-- ------------------------------------
-- Ligand Docker v2.5
-- Credits: Raven_pl,Susume,LociOiling,lynnai,
-- robgee,jeff101
-- last updated April 9 2024 428am
-- ------------------------------------
-- 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.
-- 8/10/2020 for v2.3, textboxes not sliders now input
-- minDist & maxDist. Let WigBandCI vary.
-- 10/17/2022 for v2.4, change default
-- inputs to jeff101's favorites for 2.3a3
-- and then test changes made between
-- *2.3a3.txt 8/10/2020 657am and
-- *2.3a6.txt 8/19/2020 648am versions.
-- 2/17/2024 goals for v2.5:
-- At end, copy key quicksaves into undo graph & ir_solution files.
-- Let user choose if CI values are biased toward lower values.
-- Include atom #'s in hbond lists.
-- List bonus/objective/filters along with scores.
-- Sometimes flip the ligand 0 to 180 degrees.
-- Sometimes revert to lower-scoring quicksaves.
-- Revise sorting of quicksaves.
-- Change Verbosity for some messages.
-- ------------------------------------
-- Ideas for later:
-- Break up wiggling without bands into shorter intervals.
-- If gain doesn't pass a threshold, quit wiggling.
-- If gain passes a threshold, continue wiggling.
-- Have LevyCycles<0 subtract Levy from minDist & maxDist
-- and perhaps reverse Levy if minDist & maxDist
-- go out of bounds.
-- ------------------------------------
-- ---------- GLOBAL VARS -------------
-- ------------------------------------
Recipe = "Ligand Docker"
Version = "v2.5"
ReVersion = Recipe .. " " .. Version
OrigBands = band.GetCount() -- # of user-supplied bands at start
OrigBandsOn = 0 -- # of enabled user-supplied bands at start
OrigBandsOff = 0 -- # of disabled user-supplied bands at start
for jj=1,OrigBands do
if band.IsEnabled(jj) then
OrigBandsOn = OrigBandsOn +1
else
OrigBandsOff = OrigBandsOff + 1
end -- if band
end -- for jj
-- -
-- ------------- OPTIONS --------------
-- -
Cycles = 10000 --was 30 before 8/7/2020.
RevertCycles = 1 --# of cycles between revert to best score.
--0 means never revert to best score.
--was 5 before 8/9/2020.
--was 2 in v2.3a3 & v2.3a6.
--changed to 1 in v2.4a1 10/17/2022.
RevertFac = 1500 -- 0 means revert to best score only.
-- > 0 means favor best score, but sometimes
-- revert to lower scores from quicksaves
-- instead. Larger RevertFac lets lower
-- scoring quicksaves be used. Like kT
-- in a partition function, but for Foldit
-- scores instead of energy and temperature.
FlipProb = 0.6666 -- Probability of flipping the ligand.
-- 0 means never flip. 1 means always flip.
FlipAngMin = 30 -- lowest allowed flip angle in degrees.
FlipAngMax = 180 -- highest allowed flip angle in degrees.
FlipHubProb = 0.5 -- Probability to use a midpoint (hub) when flip ligand.
-- 0 means never use midpoint (hub) when flip ligand.
-- 1 means always use midpoint (hub) when flip ligand.
LevyCycles = 0 --# of cycles between uses of Levy.
--0 means never use Levy.
--was 8 before 8/9/2020.
--was 3 before 8/10/2020.
--was 3 in v2.3a3 & 0 in v2.3a6.
minLevy = 500
maxLevy = 800
WigWithBands = 10 --how long to wiggle with bands.
--was 1 before 8/7/2020.
--was 30 before 7/17/2023.
Wiggles = 80 --how long to wiggle with no bands.
--was 3 before 8/7/2020.
--was 200 before 7/17/2023.
BandedWigOpt = 2 --1 wiggles ligand only, 2 wiggles ligand & protein.
--default was 1 before 7/12/2023.
UnbandedWigOpt = 2 --1 wiggles ligand only, 2 wiggles ligand & protein.
--default was 2 before 7/12/2023.
AtomsInLigand = 0 --# of atoms in the ligand.
BandsToLigand = 5 --# of bands from ligand to space. was 2 before 8/7/2020.
minDist = 5 --min starting length for bands in 0.001 A
--was 10 before 8/7/2020.
maxDist = 6299 --max starting length for bands in 0.001 A
--was 50 before 8/10/2020.
--was 50 in v2.3a3 & 2000 in v2.3a6.
--changed to 6299 in v2.4a1 10/17/2022.
minGoal = 0 --min target length for bands in 0.001 A
maxGoal = 0 --max target length for bands in 0.001 A
Levy = 0
BandStrength = 1.50
UseRecentBest = 3 -- Used in random_bands() function below:
-- 0 means don't use recent best,
-- 1 means use recent best at end of each cycle,
-- 2 means use recent best when banded wiggle ends & again when unbanded wiggle ends,
-- 3 means use recent best just when banded wiggle ends.
LigName = '' -- A brief name for the ligand. It will appear in titles
-- for files where final quicksaves will be copied.
Verbosity = 2 -- 0 = the least output messages, 3 = the most output messages
-- was 0 before 8/10/2020.
-- v2.3a3 had 0, v2.3a6 had 2.
-- are below variables still needed? 7/12/2023
CapCI = false
userSetMaxCI = behavior.GetClashImportance()
WigBandCI = 0.88 -- the CI to use for wiggling with bands
-- are above variables still needed? 7/12/2023
atoms4orig_str = ""
atoms4orig_list = {} -- will be a list of 1's & 0's for each atom in the ligand.
-- 1 means can band to it in original banded wiggle way. 0 means don't band to it.
atoms4flip_str = ""
atoms4flip_list = {} -- will be a list of 1's & 0's for each atom in the ligand.
-- 1 means can band to it in flipping banded wiggle way. 0 means don't band to it.
-- Below are for the Behavior Sliders and Importances
cinams={'Clash','SidechainHBond','BackboneHBond','Packing','Hiding','Pairwise','Density'}
-- cinams are the official internal LUA names for each slider
ciabbr={'Clashing','Sidechain H-Bonds','Backbone H-Bonds','Packing','Hiding','Pairwise','Density'}
ciabbr2={'Clashing','Sidechain H-B','Backbone H-B','Packing','Hiding','Pairwise','Density'}
ciabbr3={'Clas','Side','Back','Pack','Hidi','Pair','Dens'}
-- ciabbr are user-friendly versions of the names in cinams
-- ciabbr2 & ciabbr3 are even shorter user-friendly names
cimins={0,0,0,0,0,0,0} -- min for each Importance
cimaxs={1,3,3,3,3,3,3} -- max for each Importance
cidefs={1,1,1,1,1,1,1} -- default values for each Importance
ciolds={}
cilos={}
cihis={}
civals={}
for jj=1,7 do
ciolds[jj]=behavior['Get'..cinams[jj]..'Importance']()
cilos[jj]=cimins[jj]
cihis[jj]=cimaxs[jj]
civals[jj]=ciolds[jj]
end -- for jj
cibias=1 -- 0 means no bias, so ci values (0 to 1/2) & (1/2 to 1) are equally likely,
-- and ci values (0 to 1) & (1 to 2) & (2 to 3) are all equally likely.
-- 1 means use bias so that lower ci values are more likely:
-- for example, ci values (0 to 1/3) & (1/3 to 1) are equally likely,
-- and ci values (0 to 1) & (1 to 3) are equally likely.
-- In effect, cibias was 0 before v2.5 in 2024.
--
-- To do cibias==1 right, let x be 0 to 1 or 0 to 2 with all x values equally likely.
-- When x=0 1/2 1 or 2, we want ci = y to have the values
-- y=0 1/3 1 or 3 respectively.
-- Using y=x^n, we get 0=0^n, (1/3)=(1/2)^n, 1=1^n, and 3=2^n.
-- Using logarithms, we get ln(3) = ln(2^n) = n*ln(2)
-- so that ln(3)/ln(2) = n ~ 1.5849625.
cibiasn1=math.log(3.0)/math.log(2.0) -- for natural logarithms ln(3)/ln(2)
cibiasn2=math.log10(3.0)/math.log10(2.0) -- for base-10 logarithms
print('cibias = '..cibiasn1..' = '..cibiasn2..' ~ 1.5849625') -- all outputs should be the same
cibiasn=cibiasn1
-- END OPTIONS
segCnt = structure.GetCount()
numSegs =structure.GetCount()
BASESCORE = current.GetEnergyScore()
highscore = BASESCORE
--
-- This recipe uses Quicksave (QS) slots 1-4 as below:
-- QS1 stores OG pose for backup ('Start pose')
-- QS2 stores best solution for the run ('Best pose')
-- QS3 stores VeryBest solution for the run ('Top score')
-- QS4 is a temporary storage slot
--
QS3score = BASESCORE
-- setnoteseg() -- not yet a defined function, so run it in Main() instead
-- setendnote() -- not yet a defined function, so run it in Main() instead
save.Quicksave(3) -- store current solution in QS3
--
-- QS slots QSlo to QShi will hold best results in the run.
QSmin=5
QSmax=100
QSlo=QSmin
QShi=QSmax
QSInUse=QSlo-1 -- start with slot QSlo-1,
-- then increment up to slot QShi.
-- once reach slot QShi, start re-using slots QSlo to QShi.
QSdat={} -- QSdat will hold lots of data for slots 1 to 3 and QSlo to QShi
-- -------------------------------------------
-- List all the recipe's functions below:
-- -------------------------------------------
function setnoteseg()
-- uses numSegs to find SegForNote.
-- SegForNote is the lowest segment #
-- with no Tab Notes. There should be no Tab Notes
-- for segments SegForNote to numSegs.
local seg
SegForNote=numSegs
for seg=numSegs,1,-1 do
if structure.GetNote(seg)~="" then break end
SegForNote=seg
end -- for seg
end -- setnoteseg()
-- puts final Note in segment SegForNote.
-- seems should do this right before every
-- save of a new result during a run.
function setendnote()
local mssg
mssg=string.format("(%s) %.3f + %s -> %.3f",
user.GetPlayerName(),BASESCORE,
ReVersion,current.GetEnergyScore())
structure.SetNote(SegForNote,mssg)
end -- setendnote()
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()
-- below is for printing scores the way
-- the Foldit gui seems to do so 2/24/2024
function prtsc(x)
return string.format('%.3f',x)
end -- prtsc()
function round ( x )--to 3 places
return x - x % 0.001
end -- round()
-- below sets ci between cilo and cihi
-- and biases ci to lower values if cibias = 1
function getci(cilo,cihi)
local n,invn,xlo,xhi,xval,cival
if cibias==0 then
cival=rndreal(cilo,cihi)
else -- bias ci to lower values
-- so that cival=xval^n and xval=cival^(1/n)=cival^invn
-- n below should be ln(3)/ln(2) = log(3)/log(2) = 1.5849625...
n=cibiasn
invn=1/n
xlo=cilo^invn
xhi=cihi^invn
xval=rndreal(xlo,xhi)
cival=xval^n
end -- if cibias
return cival
end -- getci()
-- below finds a real number [vlo,vhi]
-- ranging from vlo to vhi
function rndreal(vlo,vhi)
local tlo,thi,val,isok,trng
-- below sets tlo to the min of vlo,vhi
-- and thi to the max of vlo,vhi
tlo=math.min(vlo,vhi)
thi=math.max(vlo,vhi)
--
if tlo==thi then
val=tlo
else
trng=(thi-tlo)
isok=0 -- 0 means val is not correct yet.
-- isok will be 1 when val is correct.
while isok==0 do
-- math.random() gives random float [0,1)
-- from 0 to just below 1.
-- Below is to compensate for this, so it is
-- possible to get val equal to vlo vhi
-- or a random # in between.
val=tlo-0.0001*trng+math.random()*1.0002*trng
-- above should give some val values
-- outside tlo to thi range.
-- below will reject the ones outside
-- the tlo to thi range by
-- keeping isok set to 0.
if val>=tlo then
if val<=thi then
isok=1 -- val is correct now!
end -- if val
end -- if val
end -- while isok
end -- if tlo
--
return val
end -- rndreal()
-- listhbonds() below was inspired by
-- https://foldit.fandom.com/wiki/Foldit_Lua_Function_structure.GetHBonds
--
-- howgot tells how the current score was obtained:
-- 3 means at the end of an unbanded wiggle
-- 2 means as a recent best from an unbanded wiggle
-- 1 means anything from a banded wiggle (when in doubt, use this setting)
-- 0 means no need to store it in a quicksave
--
-- some messages will only print if Verbosity >= checkverb
--
-- In below, qs123 is an integer from 1 to 3 for quicksaves 1 to 3.
-- qs123 is 0 for quicksaves 5 to QSInUse.
-- quicksaves 1 to 3 need special treatment here,
-- so when the recipe ends, can sort & list qs's 1-3
-- & put them in undo graph & save them
-- as I do for qs's 5 to QSInUse.
function listhbonds(howgot,checkverb,qs123)
local bonds,nbonds,bonds2list,ngot
local nsegs,ii,jj,res,listit,resnum,atoms
local tstr,tstr2,tstr3,tmpscore,tmpbonus
local iires,gotit,slotused,worstslot
local numremoved,totbonus,bonustxt
--
-- next line lists the hbonds selected for
-- viewing, like the non-protein ones
bonds=structure.GetHBonds()
--
nbonds=(#bonds) -- # of hbonds found
--
nsegs=structure.GetCount() -- # of segments in protein
-- Usually the ligand is at segment nsegs
--
-- In bonds2list below want hbonds going
-- between segment nsegs (the ligand) and some other segment.
-- This will need changes if you ever get 2 ligands in same puzzle,
-- hbonds ligand makes with itself, or disulfide bonds.
bonds2list={}
ngot=0
for ii=1,nbonds do
res={}
atoms={}
-- print('ii='..ii)
res[1]=bonds[ii].res1
-- print('res[1]='..res[1])
res[2]=bonds[ii].res2
-- print('res[2]='..res[2])
atoms[1]=bonds[ii].atom1
atoms[2]=bonds[ii].atom2
listit=0 -- 0 means don't list hbond ii
-- 1 means list res[1], which ~= nsegs
-- 2 means list res[2], which ~= nsegs
if bonds[ii].bond_type == 0 then -- it is an hbond
-- bond_type is 1 for disulfides, which I don't want to list here
if res[1] == nsegs then
if res[2] ~= nsegs then
listit=2
end -- if res[2]
else -- res[1] ~= nsegs
if res[2] == nsegs then
listit=1
end -- if res[2]
end -- if res[1]
end -- if bonds[ii].bond_type
if listit > 0 then
ngot=ngot+1
bonds2list[ngot]={}
bonds2list[ngot][0]=res[listit]
bonds2list[ngot][1]=atoms[listit]
bonds2list[ngot][2]=atoms[3-listit] -- 3-listit gives 1 for 2 or 2 for 1
end -- if listit
end -- for ii
--
-- below outputs the results
tstr3='' -- Will be like ' w148-4-7 r193-2-9 r198-2-12 r198-4-13' at the end.
-- In w148-4-7 above, w is AA's 1-letter-code for tryptophan,
-- 148 is AA's segment/residue #, -4 is atom #4 on the AA w148,
-- and -7 is atom #7 on the ligand.
-- Should you worry about equivalent atoms on the same segment,
-- like O's in -COO groups or H's in -NH2 or -NH3 groups?
if qs123==0 then
if Verbosity >= checkverb then
totbonus=filter.GetBonusTotal()
if totbonus>=0 then
bonustxt=string.format('+%.1f',totbonus)
else
bonustxt=string.format('%.1f',totbonus)
end -- if totbonus
-- see https://foldit.fandom.com/wiki/Foldit_Lua_Function_filter.GetBonus
print(
string.format('For %s %s w/bonus %s and torsion %.1f, ',
LigName,prtsc(current.GetEnergyScore()),
bonustxt,filter.GetBonus('Torsion Quality'))..
string.format('%s%d has score=%.1f\n (clashing=%.1f pairwise=%.1f ',
structure.GetAminoAcid(nsegs),nsegs,
current.GetSegmentEnergyScore(nsegs)*20.0, -- 20 seems the ratio needed 4/9/2024
current.GetSegmentEnergySubscore(nsegs,'Clashing'),
current.GetSegmentEnergySubscore(nsegs,'Pairwise'))..
string.format('packing=%.1f hiding=%.1f bonding=%.1f ideality=%.1f)\n and ',
current.GetSegmentEnergySubscore(nsegs,'Packing'),
current.GetSegmentEnergySubscore(nsegs,'Hiding'),
current.GetSegmentEnergySubscore(nsegs,'Bonding'),
current.GetSegmentEnergySubscore(nsegs,'Ideality'))..
string.format('makes %d hbonds to other segs w/the bonding subscores below:',
ngot))
-- someday might list s-s bonds, which use 'Disulfides' instead of 'Bonding' above.
-- https://foldit.fandom.com/wiki/Foldit_Lua_Function_current.GetSegmentEnergySubscore
-- gives more general code than above, if other subscores become more common.
end -- if Verbosity
end -- if qs123
tstr=' '
if ngot <= 0 then
tstr=(tstr..' none')
else
for ii=1,ngot do
resnum=bonds2list[ii][0]
tstr2=(structure.GetAminoAcid(resnum)..resnum)
tstr2=string.format('%s-%d-%d',tstr2,bonds2list[ii][1],bonds2list[ii][2])
tstr=string.format('%s %s=%.1f',
tstr,tstr2,current.GetSegmentEnergySubscore(resnum,'Bonding'))
-- someday might list s-s bonds, which use 'Disulfides' instead of 'Bonding' above.
tstr3=(tstr3..' '..tstr2)
end -- for ii
end -- if ngot
if qs123==0 then
if Verbosity >= checkverb then
print(tstr..'.')
end -- if Verbosity
end -- if qs123
--
if howgot > 0 then
-- below will store key ones in QSdat and a QS slot if needed
tmpscore=current.GetEnergyScore()
tmpbonus=filter.GetBonusTotal()
if qs123>0 then
storeinslot(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,qs123)
elseif QSInUse==QSlo-1 then
-- store it all in slot QSInUse+1 (slot QSlo here)
QSInUse=QSInUse+1
storeinslot(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,QSInUse)
else -- QSInUse must be QSlo to QShi here
--- check if in a previous slot and store it there
gotit=0 -- 0 means not yet found in previous slot
for ii=QSlo,QSInUse do
if gotit==0 then
iires=isnewinold(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,ii)
if iires >= 1 then
-- new is in old
gotit=iires -- gotit should be 1 or 2 now
if iires == 2 then
-- new is in old and outdoes old
storeinslot(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,ii)
slotused=ii
end -- if iires
end -- if iires
end -- if gotit
end -- for ii
-- gotit should be 0 1 or 2 now
if gotit==2 then -- new result was put in slot w/# slotused.
-- should see if new result overrides other slots and
-- remove those slots.
numremoved=0 -- # slots removed so far
for ii=QSInUse,slotused+1,-1 do -- for ii=QSInUse to slotused+1 step -1
iires=isnewinold(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,ii)
if iires == 2 then -- new solution overrides one in slot ii
numremoved=numremoved+1
print('Removing Quicksave slot '..ii..
' and shifting larger-numbered slots down.')
for jj=ii,QSInUse-numremoved do
copyslot1to2(jj+1,jj)
end -- for jj
QSdat[QSInUse-numremoved+1]=nil
-- jeff101 doesn't think Foldit has a LUA command
-- to remove Quicksaves. If it did, here would be
-- a good place to remove Quicksave slot #
-- QSInUse-numremoved+1 (8/2/2023).
end -- if iires
end -- for ii
if numremoved > 0 then
QSInUse=QSInUse-numremoved
print('After removing '..numremoved..
' slots, only Quicksave slots up to '..QSInUse..
' are still in use.')
end -- if numremoved
elseif gotit==0 then -- new result wasn't in a previous slot
if QSInUse<QShi then
-- store new result in slot QSInUse+1
QSInUse=QSInUse+1
storeinslot(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,QSInUse)
else -- There are no more QS slots to use,
-- so find the worst slot, and if the new result
-- outdoes the worst slot, put the new result there.
-- jeff101 doubts this case will ever happen 8/2/2023.
worstslot=findworst()
iires=isnewbetter(tmpscore,tmpbonus,ngot,howgot,worstslot)
if iires==1 then
storeinslot(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,worstslot)
end -- if iires
end -- if QSInUse
end -- if gotit
end -- if qs123
end -- if howgot
--
end -- listhbonds()
function copyslot1to2(n1,n2)
local jj,kk
QSdat[n2]={}
QSdat[n2].score=QSdat[n1].score
QSdat[n2].bonus=QSdat[n1].bonus
QSdat[n2].howgot=QSdat[n1].howgot
QSdat[n2].ngot=QSdat[n1].ngot
QSdat[n2].hbstr=QSdat[n1].hbstr
QSdat[n2].hblist={}
for jj=1,QSdat[n1].ngot do
QSdat[n2].hblist[jj]={}
for kk=0,2 do
QSdat[n2].hblist[jj][kk]=QSdat[n1].hblist[jj][kk]
end -- for kk
end -- for jj
setendnote()
save.Quicksave(4) -- store present pose into slot 4
save.Quickload(n1)
setendnote()
save.Quicksave(n2)
save.Quickload(4) -- restore present pose from slot 4
end -- copyslot1to2()
-- 3/8/2024 Should the function below save QS 1-3 to files ? yes
-- 3/8/2024 Should the function below list QS 1-3 ? yes
function listallqs()
local sortedqs,ii,jj,kk,gotit,tstr,t1,t2,wigp,ftitles
local scores,qsnums,ngot,nuscores,indx,qsgot,qs2do,qii
qs2do={}
for ii=1,3 do
qs2do[ii]=ii
end -- for ii
qsgot=3
if QSInUse==QSlo-1 then
print('No quicksaves from '..QSlo..' to '..QShi..' are in use,')
print(' so sorting quicksaves 1 to 3 only.')
else
for ii=QSlo,QSInUse do
qsgot=qsgot+1
qs2do[qsgot]=ii
end -- for ii
print('Sorting '..qsgot..' quicksaves (1 to 3 and '..QSlo..' to '..QSInUse..').')
end -- if QSInUse
sortedqs={}
sortedqs[1]=qs2do[1]
for ii=2,qsgot do
qii=qs2do[ii]
gotit=0
for jj=1,ii-1 do -- jj < ii holds
if gotit==0 then
if isnewbetter(QSdat[qii].score,QSdat[qii].bonus,
QSdat[qii].ngot,QSdat[qii].howgot,sortedqs[jj])==1 then
-- above means ii is better than jj & so should go before jj
-- so all from jj to ii-1 will shift up to later in the list
-- and jj will be replaced by ii
for kk=ii-1,jj,-1 do -- for kk=ii-1 to jj step -1
sortedqs[kk+1]=sortedqs[kk]
end -- for kk
sortedqs[jj]=qii
gotit=1
end -- if isnew
end -- if gotit
end -- for jj
if gotit==0 then
sortedqs[ii]=qii
end -- if gotit
end -- for ii
print(' ')
--
wigp=behavior.GetWigglePower()
if wigp=='a' then
wigp='auto'
elseif wigp=='l' then
wigp='lo'
elseif wigp=='m' then
wigp='med'
elseif wigp=='h' then
wigp='hi'
end -- if wigp
--
scores={}
qsnums={}
ftitles={}
ngot=0
print('QS# how #hb score bonus hbonds:')
for ii=1,qsgot do
jj=sortedqs[ii]
ngot=ngot+1
scores[ngot]=QSdat[jj].score
qsnums[ngot]=jj
if QSdat[jj].bonus>=0 then
tstr=string.format('+%.1f',QSdat[jj].bonus)
else
tstr=string.format('%.1f',QSdat[jj].bonus)
end -- if QSdat
-- %8s below right justifies tstr
print(string.format('%3d %3d %3d %12.3f %8s %s.',
jj,QSdat[jj].howgot,QSdat[jj].ngot,
QSdat[jj].score,tstr,QSdat[jj].hbstr))
--
-- Below loads each quicksave into undo graph
-- and saves it to a puzzle*ir_solution file.
ftitles[jj]=string.format('%.3f %s %s',QSdat[jj].score,tstr,wigp)
-- Next, if LigName is 1 or more characters, then include it in ftitle.
if string.len(LigName)> 0 then
ftitles[jj]=(LigName..' '..ftitles[jj])
end -- if LigName
end -- for ii
print(' ')
--
-- Below will print the list of scores et al again
-- but in order from smallest to largest score.
nuscores,indx=sortit(scores)
-- sortit() above returns nuscores which are the
-- values in scores from highest to lowest.
setendnote()
save.Quicksave(4) -- store pre-loop pose in slot 4
print('QS# how #hb score bonus hbonds:')
for ii=ngot,1,-1 do -- for ii=ngot to 1 step -1
jj=qsnums[indx[ii]]
save.Quickload(jj) -- load pose from quicksave slot jj
if QSdat[jj].bonus>=0 then
tstr=string.format('+%.1f',QSdat[jj].bonus)
else
tstr=string.format('%.1f',QSdat[jj].bonus)
end -- if QSdat
-- %8s below right justifies tstr
print(string.format('%3d %3d %3d %12.3f %8s %s.',
jj,QSdat[jj].howgot,QSdat[jj].ngot,
QSdat[jj].score,tstr,QSdat[jj].hbstr))
-- print('Next line tries save.SaveSolution('..ftitles[jj]..').')
setendnote()
save.SaveSolution(ftitles[jj]) -- save present pose into a puzzle*ir_solution file.
-- date & time are automatically put in the comment field.
-- Seems like can only save one file per second,
-- otherwise their filenames will be the same
-- and they will overwrite each other.
t1=os.time()
repeat
t2=os.time()
until os.difftime(t2,t1)>1 -- want t1 and t2 to be at least 1 second apart
end -- for ii
save.Quickload(4) -- restore pre-loop pose from slot 4
end -- listallqs()
-- Below was adapted from DistMap1q.txt 1/24/16 312am code
-- at https://fold.it/recipes/101679 lines 813-1004.
--
-- sortit(x) should do like Matlab's sort for a vector x.
-- It should end with x[indx[i]]=s[i] and s[1]>=s[2]>=s[3] etc.
-- so s ends with x's elements in descending order.
function sortit(x)
local s={}
local sx={}
local indx={}
local i,dim
dim=#x
for i=1,dim do
sx[i]={}
sx[i][1]=x[i] -- copy x into sx
sx[i][2]=i -- store original indices
end -- for i
table.sort(sx,comp)
-- Above sorts rows of matrix sx in place from max to min
-- using all cols of sx if needed. Having orig index in the
-- last col of sx above makes equal items appear in reverse
-- their original order.
for i=1,dim do
s[i]=sx[i][1] -- retrieve sorted values
indx[i]=sx[i][2] -- retrieve sorted indices
end -- for i
return s,indx
end -- sortit()
--
-- comp() below is used with table.sort()
-- to sort things from highest to lowest.
-- Using comp(si,sj) below put things lowest to highest 1/24/16.
function comp(sj,si)
local qflag,swapem,k,ncols
-- here si=s[i] & sj=s[j]
ncols=#si
qflag=0
swapem=false -- won't swap i,j below
k=1 -- 1st try to sort using col 1
while qflag==0 do
if si[k]<sj[k] then -- use < so get largest element (here j) first
swapem=true -- will swap i,j below
qflag=1 -- will quit while loop
elseif si[k]==sj[k] then -- if can't sort w/col k, try w/col k+1
k=k+1
if k>ncols then -- have tried all cols, so i,j are same, let em be
qflag=1 -- will quit while loop
end -- if k
else -- si[k]>sj[k] already holds, so no need to swap i,j
qflag=1 -- will quit while loop
end -- if s
end -- while qflag
return swapem
end -- comp()
--
-- Above was adapted from DistMap1q.txt 1/24/16 312am code
-- at https://fold.it/recipes/101679 lines 813-1004.
function findworst()
local wslot=QSlo -- default value for worst slot
local wscore,wbonus,wgot,whow,ii,res
wscore=QSdat[wslot].score
wbonus=QSdat[wslot].bonus
wgot=QSdat[wslot].ngot
whow=QSdat[wslot].howgot
for ii=QSlo+1,QSInUse do
res=isnewbetter(wscore,wbonus,wgot,whow,ii)
-- above asks if wscore,wbonus,wgot,whow is better than ii's results
-- this is like asking if ii's results are worse than wscore,wbonus,wgot,whow
-- should get 1 if ii's result are worse
if res == 1 then -- ii is worse than wslot, so replace wslot with ii
wslot=ii
wscore=QSdat[ii].score
wbonus=QSdat[ii].bonus
wgot=QSdat[ii].ngot
whow=QSdat[ii].howgot
end -- if res
end -- for ii
print(string.format('findworst: slot %3d (%d %2d %11.3f) is the worst.',
wslot,whow,wgot,wscore))
return wslot -- return the # of the worst slot
end -- findworst()
function isnewbetter(tmpscore,tmpbonus,ngot,howgot,snum)
local res -- Will be 0 if new equals old (snum's results),
-- 1 if new is better than old (snum's results), and
-- -1 if new is worse than old (snum's results).
local oscore,obonus,ogot,ohow,symb
oscore=QSdat[snum].score
obonus=QSdat[snum].bonus
ogot=QSdat[snum].ngot
ohow=QSdat[snum].howgot
if howgot > ohow then
res = 1
elseif howgot == ohow then
if ngot > ogot then
res = 1
elseif ngot == ogot then
if tmpscore > oscore then
res = 1
elseif tmpscore == oscore then
res = 0
else -- tmpscore < oscore
res = -1
end -- if tmpscore
else -- ngot < ogot
res = -1
end -- if ngot
else -- howgot < ohow
res = -1
end -- if howgot
--[[ 2/24/2024 The above way of sorting should rank things as below:
howgot ngot score
best 3 4 4500
3 3 6000
3 3 5900
3 2 7000
1 7 8000
1 7 7500
worst 1 7 7200
]]--
if res==1 then
symb='beats'
elseif res==0 then
symb='equals'
else -- res == -1
symb='loses to'
end -- if res
if Verbosity>2 then
print(string.format('isnewbetter: new (%d %2d %12.3f)',
howgot,ngot,tmpscore))
print(string.format(' vs old=%3d (%d %2d %12.3f) gives res = %d (new %s old).',
snum,ohow,ogot,oscore,res,symb))
end -- if Verbosity
return res
end -- isnewbetter()
function storeinslot(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,snum)
local ii,kk,tstr
setendnote()
save.Quicksave(snum) -- save current solution in slot snum
tstr='Replacing Quicksave slot '..snum..' that contains'
if QSdat[snum] == nil then
print(tstr..' nothing with')
else
print(string.format('%s\n %12.3f %d %2d:%s with',
tstr,QSdat[snum].score,QSdat[snum].howgot,
QSdat[snum].ngot,QSdat[snum].hbstr))
end -- if QSdat[snum]
print(string.format(' %12.3f %d %2d:%s.',
tmpscore,howgot,ngot,tstr3))
QSdat[snum]={}
QSdat[snum].score=tmpscore
QSdat[snum].bonus=tmpbonus
QSdat[snum].ngot=ngot
QSdat[snum].hbstr=tstr3
QSdat[snum].howgot=howgot
QSdat[snum].hblist={}
for ii=1,ngot do
QSdat[snum].hblist[ii]={}
for kk=0,2 do
QSdat[snum].hblist[ii][kk]=bonds2list[ii][kk]
end -- for kk
end -- for ii
end -- storeinslot()
function isnewinold(tmpscore,tmpbonus,ngot,bonds2list,tstr3,howgot,snum)
local res = 0
-- res will end as 0 if new is not in old one (may save new one in a diff snum)
-- 1 if new is in old but is same or worse than old (will skip the new one)
-- 2 if new is in old and is better than old (will save new one into snum)
local oscore,obonus,ogot,ostr,ohow,olist,ii,kk
oscore=QSdat[snum].score
obonus=QSdat[snum].bonus
ogot=QSdat[snum].ngot
ostr=QSdat[snum].hbstr
ohow=QSdat[snum].howgot
olist={}
for ii=1,ogot do
olist[ii]={}
for kk=0,2 do
olist[ii][kk]=QSdat[snum].hblist[ii][kk]
end -- for kk
end -- for ii
--
if ngot==ogot then -- both have same # of hbonds
if tstr3 == ostr then -- both have exact same hbonds
if tmpscore > oscore then
if howgot >= ohow then
res=2 -- will save new one in snum
else
res=0 -- may save new one in another snum
end -- if howgot
elseif tmpscore == oscore then
if howgot > ohow then
res=2 -- will store in snum
else
res=1 -- will skip it
end -- if howgot
else -- tmpscore < oscore
if howgot > ohow then
res=0 -- may save new in another snum
else
res=1 -- will skip new one
end -- if howgot
end -- if tmpscore
else -- same # hbonds but diff set of hbonds
res=0
end -- if tstr3
-- 2/24/2024 jeff101 doesn't want to sort like below any more.
-- The way below pools solutions with different #'s of hbonds.
--[[
elseif ngot > ogot then
-- is old list in new list?
if isL1inL2(ogot,olist,ostr,ngot,bonds2list,tstr3)==1 then
-- yes
if tmpscore >= oscore then
if howgot >= ohow then
res=2 -- will store new in snum
else
res=0
end -- if howgot
else -- tmpscore < oscore
res=0
end -- if tmpscore
else
res=0
end -- if isL1inL2
else -- ngot < ogot
-- is new list in old list?
if isL1inL2(ngot,bonds2list,tstr3,ogot,olist,ostr)==1 then
-- yes
if tmpscore <= oscore then
if howgot <= ohow then
res=1 -- will skip new list
else
res=0
end -- if howgot
else -- tmpscore > oscore
res=0
end -- if tmpscore
else
res=0
end -- if isL1inL2
]]--
-- The way above pools solutions with different #'s of hbonds.
-- 2/24/2024 jeff101 doesn't want to sort like above any more.
end -- if ngot
if Verbosity>2 then
print(string.format(
'isnewinold: new (%d %2d %12.3f%s)\n vs old=%3d (%d %2d %12.3f%s) gives res=%d.',
howgot,ngot,tmpscore,tstr3,snum,ohow,ogot,oscore,ostr,res))
end -- if Verbosity
return res
end -- isnewinold()
-- In the function below,
-- n1 is # of items in list L1 (each item has 3 elements #'d 0-2),
-- str1 is a string listing all the elements in L1,
-- n2 is # of items in list L2 (each item has 3 elements #'d 0-2), and
-- str2 is a string listing all the elements in L2.
function isL1inL2(n1,L1,str1,n2,L2,str2)
local res, ii, jj, ngot, L1used, L2used
-- When done, this function will return res which will be
-- 2 if lists L1 & L2 are identical,
-- 1 if all elements of L1 are also in L2,
-- -1 if all elements of L2 are also in L1,
-- 0 if neither list is entirely contained in the other.
-- L1used & L2used will be lists of n1 or n2 0's.
-- These 0's will change to 1's when an item
-- is found in both L1 and L2.
if n1 > n2 then
res = -isL1inL2(n2,L2,str2,n1,L1,str1)
elseif n1 == n2 then
if str1 == str2 then
res = 2 -- lists L1 and L2 are identical
else
res = 0 -- neither list is entirely in the other list
end -- if str1
else -- n1 < n2 must hold here
res = 1 -- if below can't find all elements
-- of L1 in L2, it will change res to 0
L1used={}
for ii=1,n1 do
L1used[ii]=0
end -- for ii
L2used={}
for ii=1,n2 do
L2used[ii]=0
end -- for ii
ngot = 0 -- the # of items found so far in both lists
for ii=1,n1 do
if ii>ngot+1 then
-- this means at least one element of L1
-- has not been found in L2, so it is time to quit
res=0
else
if L1used[ii]==0 then
for jj=1,n2 do
if L1used[ii]==0 then -- as jj rises
-- this can change from 0 to 1
if L2used[jj]==0 then
if L1[ii][0]==L2[jj][0] then -- protein segment #
-- The above checks that L1 & L2 use the same segment #'s.
-- Below checks if L1 & L2 use the same atom #'s (new to v2.5).
if L1[ii][1]==L2[jj][1] then -- protein atom #
if L1[ii][2]==L2[jj][2] then -- ligand atom #
ngot=ngot+1
L1used[ii]=1
L2used[jj]=1
end -- if L1[ii][2]
end -- if L1[ii][1]
end -- if L1[ii][0]
end -- if L2used
end -- if L1used
end -- for jj
end -- if L1used
end -- if ii
end -- for ii
if ngot<n1 then
res=0
end -- if ngot
end -- if n1
if Verbosity>2 then
print(string.format(
'isL1inL2: Are all elements of L1 (%2d:%s)',n1,str1))
print(string.format(
' inside L2 (%2d:%s)? res=%d.',n2,str2,res))
end -- if Verbosity
return res
end -- isL1inL2()
-- -----------------------------------------------------
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 (by LociOiling?)
-- above keeps seed between 10000000
-- and 4294967296=2^32
-- below rounds seed to an integer
seed = seed - seed % 1
return seed
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()
-- is below line still needed? 7/12/2023
behavior.SetClashImportance( userSetMaxCI )
-- is above line still needed? 7/12/2023
local kk
for kk=1,7 do
behavior['Set'..cinams[kk]..'Importance'](ciolds[kk])
end -- for kk
print( "Final score: " ..prtsc( highscore ) )
print( "Total gain: " ..prtsc( highscore - BASESCORE ) )
listhbonds(0,0,0)
setendnote() -- puts final Note in segment SegForNote
listallqs()
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 prep_to_flip_ligand()
local atom1,atom2,ligseg,bos,b1f,b2f,b1m,b2m,b1,b2
local d1t2,d2t1,angdeg,angrad,olen,flen,lenb1m,lenb2m
local hd1t2,hd2t1,hpi,ds,dslen,dstr,th,ph,jj
local Xatom,Yatom,Oatom,Xaxis,Yaxis
ligseg=LIGAND -- segment # for ligand
--
-- Below picks 2 different atoms (atom1 & atom2) on the ligand
-- that will try to swap positions. Make sure they are both
-- allowed atoms (ones with 1's in atoms4flip_list).
repeat
atom1=math.random(1,AtomsInLigand)
until atoms4flip_list[atom1]==1
repeat
repeat
atom2=math.random(1,AtomsInLigand)
until atoms4flip_list[atom2]==1
until atom1~=atom2
--
-- Pick options for band.Add so bos[1] & bos[2] point in random directions
-- to 0.0001 in space away from positions of atom1 & atom2.
bos={}
for jj=1,2 do
th = rndreal(0.0,math.pi) -- gives real # from 0 to pi, including the end points
ph = rndreal(0.0,2*math.pi) -- gives real # from 0 to 2*pi, including the end points
Xaxis=math.random( 1, ligseg-1 ) -- gives integers 1 to ligseg-1
repeat
Yaxis=math.random( 1, ligseg-1 ) -- gives integers 1 to ligseg-1
until Yaxis ~= Xaxis
-- here should have ligseg, Xaxis, and Yaxis be 3 different integers
Xatom = 0
Yatom = 0
if jj==1 then
Oatom = atom1
else
Oatom = atom2
end -- if jj
bos[jj]=band.Add(ligseg, Xaxis, Yaxis, 0.0001, th, ph, Oatom, Xatom, Yatom )
end -- for jj
b1f=band.AddToBandEndpoint(ligseg,bos[2],atom1) -- band from atom1 to space near atom2
b2f=band.AddToBandEndpoint(ligseg,bos[1],atom2) -- band from atom2 to space near atom1
d1t2=band.GetLength(b1f)
d2t1=band.GetLength(b2f)
angdeg=rndreal(FlipAngMin,FlipAngMax) -- angdeg should be an angle between 0 and 180 degrees
angrad=math.rad(angdeg) -- angrad is angdeg converted to radians, so 0 to pi radians
olen=d1t2*math.sqrt(0.5*(1.0-math.cos(angrad)))
flen=d1t2*math.sqrt(0.5*(1.0+math.cos(angrad)))
-- angdeg= 0 degrees should give olen=0 & flen=d1t2
-- angdeg= 90 degrees should give olen=flen=d1t2/sqrt(2)
-- angdeg=180 degrees should give olen=d1t2 & flen=0
band.SetGoalLength(bos[1],olen) -- can set length to 0.0 to 10000.0
band.SetGoalLength(bos[2],olen)
band.SetGoalLength(b1f,flen)
band.SetGoalLength(b2f,flen)
band.SetStrength(bos[1],BandStrength) -- can set strength to 0.0 to 10.0
band.SetStrength(bos[2],BandStrength)
band.SetStrength(b1f,BandStrength)
band.SetStrength(b2f,BandStrength)
print(string.format(
' Trying to flip ligand atoms %d and %d by %.2f degrees.',atom1,atom2,angdeg))
print(' 180 degrees should swap them completely, but smaller angles move them less.')
print(string.format(
' Since d1t2=%.3f and d2t1=%.3f here, we get olen=%.3f and flen=%.3f.',
d1t2,d2t1,olen,flen))
--
-- FlipHubProb should be a real # from 0 to 1, including the endpoints.
-- math.random() gives a real # from 0 to 1 that can equal 0 but not 1.
if math.random() < FlipHubProb then -- use the midpoint below as a hub.
-- Below finds the midpoint in space between atom1 and atom2.
-- Want to keep atom1 & atom2 a distance hd1t2=d1t2/2 from the midpoint.
hpi=math.pi/2.0
b1={}
b1[1]=band.Add(ligseg,Xaxis,Yaxis,1.0,hpi,0.0,atom1,Xatom,Yatom)
b1[2]=band.Add(ligseg,Xaxis,Yaxis,1.0,hpi,hpi,atom1,Xatom,Yatom)
b1[3]=band.Add(ligseg,Xaxis,Yaxis,1.0,0.0,0.0,atom1,Xatom,Yatom)
ds={}
b2={}
for jj=1,3 do
b2[jj]=band.AddToBandEndpoint(ligseg,b1[jj],atom2)
ds[jj]=band.GetLength(b2[jj])
ds[jj]=(d1t2^2 + 1.0 - ds[jj]^2)/2.0
band.Delete(b2[jj])
end -- for jj
dslen=0
for jj=3,1,-1 do -- remove the bands in reverse the order you made them
band.Delete(b1[jj])
dslen=dslen+ds[jj]^2
end -- for jj
dslen=math.sqrt(dslen)
-- th=math.acos(ds[3]/d1t2) -- 0 to pi in radians -- code before 4/6/2024
-- how can d1t2 differ from dslen? shouldn't they be very equal?
th=math.acos(ds[3]/dslen) -- 0 to pi in radians --
ph=math.atan2(ds[2],ds[1]) -- want 0 to 2*pi in radians but atan2 gives -pi to +pi
if ph<0.0 then
ph=ph+2.0*math.pi
end -- if ph
if math.abs(dslen-d1t2)>0.001 then
dstr=' ***'
else
dstr=''
end -- if math.abs
print(' Using a hub band to the midpoint in space between ligand atoms '
..atom1..' and '..atom2..'.')
print(string.format(
' For ds[1-3]=%.3f,%.3f,%.3f, dslen,d1t2,d2t1=%.3f,%.3f,%.3f, and dslen-d1t2=%.3f%s,',
ds[1],ds[2],ds[3],dslen,d1t2,d2t1,dslen-d1t2,dstr))
print(string.format(' we get th,ph=%6.3f,%6.3f in radians,',th,ph))
print(string.format(
' and th,ph=%6.2f,%6.2f in degrees.',math.deg(th),math.deg(ph)))
hd1t2=0.5*d1t2
hd2t1=0.5*d2t1
b1m=band.Add(ligseg,Xaxis,Yaxis,hd1t2,th,ph,atom1,Xatom,Yatom)
-- Above makes b1m be a band from atom1 to midpoint.
-- It wants th = 0 to pi radians and ph = 0 to 2*pi radians.
b2m=band.AddToBandEndpoint(ligseg,b1m,atom2) -- b2m is band from atom2 to midpoint
band.SetGoalLength(b1m,hd1t2)
band.SetGoalLength(b2m,hd1t2)
band.SetStrength(b1m,BandStrength)
band.SetStrength(b2m,BandStrength)
lenb1m=band.GetLength(b1m)
lenb2m=band.GetLength(b2m)
print(string.format(' hd1t2=%.3f, hd2t1=%.3f, lenb1m=%.3f, and lenb2m=%.3f also hold.',
hd1t2,hd2t1,lenb1m,lenb2m))
end -- if math.random()
end -- prep_to_flip_ligand()
function random_bands()
local tstr, kk, anum, howgot
local midscore -- will be the score between banded & unbanded wiggle stages
local hscore = current.GetScore() -- should this be current.GetEnergyScore() instead ?
print("Cycle start score: "..prtsc( hscore ) )
listhbonds(0,0,0)
setendnote()
recentbest.Save()
local OG = LIGAND
if BandStrength >10.0 then BandStrength = 10.0 end
------------ below preps for banded wiggle -------------
if math.random() < FlipProb then
-- FlipProb should be a real # from 0 to 1, including the endpoints.
-- math.random() gives a real # from 0 to 1 that can equal 0 but not 1.
print('Flipping the ligand -- this overrides Levy and other settings:')
prep_to_flip_ligand()
else
print('Doing banded wiggle the original way:')
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 AtomsInLigand
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>2 then
if origval==0 then
print(' Banding the '..ntoget..' ligand atoms below:')
else -- origval is 1
print(' Banding all but the '..ntoget..' ligand atoms below:')
end
end -- if Verbosity
while ngot<ntoget do
repeat
anum=math.random(1,AtomsInLigand)
until atoms4orig_list[anum]==1-origval
-- math.random above returns an integer for anum
-- [1,AtomsInLigand] from 1 to AtomsInLigand.
-- repeat ... until runs until atoms4orig_list[anum] equals 1-origval.
-- origval=0 should seek anum's with atoms4orig_list[anum]==1
-- (atoms allowed to band).
-- origval=1 should seek anum's with atoms4orig_list[anum]==0
-- (atoms not allowed to band).
if HasBand[anum]==origval then
ngot=ngot+1
HasBand[anum]=1-origval
if Verbosity>2 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 -- this means atom # anum has a band to it
ngot=ngot+1
local RHO = math.random( minDist, maxDist ) / 1000 --RHO is a distance in Angstroms
-- while Levy, minDist, maxDist are distances in 0.001 Angstroms
local theta = rndreal(0.0,math.pi) -- gives 0 to pi
local phi = rndreal(0.0,2*math.pi) -- gives 0 to 2*pi
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
local dst2 = math.random( minGoal, maxGoal ) / 1000
band.SetGoalLength(qed,dst2)
-- before 8/19/2020 set band qed's goal length to 0.
-- after 8/19/2020 set it to random value in Angstroms
-- from minGoal/1000 to maxGoal/1000.
band.SetStrength(qed,BandStrength)
if Verbosity>0 then
print(string.format(
' (%3d) Band ligand atom %3d with start and goal lengths=%7.3f and %7.3f Angstroms.',
ngot,anum,RHO,dst2))
end -- if Verbosity
end -- if qed
end -- if HasBand
end -- for anum
-- ngot should equal BandsToLigand here
-- make all BandsToLigand bands above
end -- if FlipProb
------------ above preps for banded wiggle -------------
selection.DeselectAll ()
selection.Select ( OG ) -- selects the ligand
CI(WigBandCI) -- sets clashing importance for wiggling with bands
-- for banded wiggles do
print('Setting Importances as below:')
tstr=' '
for kk=1,7 do
civals[kk]=getci(cilos[kk],cihis[kk],cibias)
tstr=string.format('%s %s=%5.3f',tstr,ciabbr3[kk],civals[kk])
behavior['Set'..cinams[kk]..'Importance'](civals[kk])
end -- for kk
print(tstr..'.')
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.
if BandedWigOpt == 1 then -- this was default before 7/12/2023
structure.WiggleSelected(witers, true, true) -- moves ligand through space and bends it.
else
structure.WiggleAll(witers , true, true) -- why not use this line? jeff101 7/9/2023
end -- if BandedWigOpt
midscore = current.GetEnergyScore()
delete_temp_bands()
if Verbosity>=0 then
print('After wiggling with bands: '..prtsc(midscore))
end -- if Verbosity
listhbonds(1,0,0) -- 2nd # here means some messages print only if Verbosity >= 0
-- 1 above sets howgot=1 inside listhbonds()
-- I think here recentbest.GetEneryScore() >= midscore & hscore
if recentbest.GetEnergyScore() > QS3score then
QS3score = recentbest.GetEnergyScore()
setendnote()
save.Quicksave(4) -- temporarily store present structure in QS4
recentbest.Restore() -- load recent best's structure
delete_temp_bands()
setendnote()
listhbonds(1,0,3) -- use howgot=1 and save recent best's structure into QS3
save.Quickload(4) -- restore present structure from QS4
-- present structure should have score midscore
end -- if recentbest
if UseRecentBest==2 or UseRecentBest==3 then
-- midscore <= recentbest.GetEnergyScore() always holds,
-- but below uses recent best only when it beats hscore.
if hscore < recentbest.GetEnergyScore() then
recentbest.Restore()
delete_temp_bands()
midscore = current.GetEnergyScore()
if Verbosity>=0 then
print('Using recent best score: '..prtsc(midscore))
end -- if Verbosity
listhbonds(1,0,0) -- 2nd # here means some messages print only if Verbosity >= 0
-- 1 above sets howgot=1 inside listhbonds()
end -- if hscore
end -- if UseRecentBest
CI(1.0) -- sets clashing importance for wiggling without bands
-- for unbanded wiggles do
print('Setting all Importances to 1.')
for kk=1,7 do
behavior['Set'..cinams[kk]..'Importance'](1)
end -- for kk
--structure.ShakeSidechainsSelected( 2 )
if UnbandedWigOpt == 2 then -- this was default before 7/12/2023
structure.WiggleAll( Wiggles , true, true )
-- above moves ligand thru space and lets it bend.
-- above also moves protein backbone and sidechains.
else
structure.WiggleSelected( Wiggles, true, true )
-- above just moves ligand thru space and lets it bend.
end -- if UnbandedWigOpt
newscore = current.GetEnergyScore()
print('After wiggling without bands: '..prtsc(newscore))
listhbonds(3,0,0)
-- 3 above sets howgot=3 inside listhbonds()
howgot=3
-- I think here recentbest.GetEneryScore() >= newscore
if recentbest.GetEnergyScore() > QS3score then
QS3score = recentbest.GetEnergyScore()
setendnote()
save.Quicksave(4) -- temporarily store present structure in QS4
recentbest.Restore() -- load recent best's structure
delete_temp_bands()
setendnote()
print(string.format('Got to A w/howgot=%d and QS3score=%s.',
howgot,prtsc(QS3score)))
listhbonds(2,0,3) -- use howgot=2 and save recent best's structure into QS3
save.Quickload(4) -- restore present structure from QS4
-- present structure should have score newscore
end -- if recentbest
if UseRecentBest==1 or UseRecentBest==2 then
if newscore < recentbest.GetEnergyScore() then
recentbest.Restore()
delete_temp_bands()
newscore = current.GetEnergyScore()
print(string.format('Got to B w/howgot=%d.',howgot))
print('Using recent best score: '..prtsc(newscore))
listhbonds(2,0,0)
-- 2 above sets howgot=2 inside listhbonds()
howgot=2
end -- if newscore
end -- if UseRecentBest
if newscore > highscore then
setendnote()
print(string.format('Got to C w/howgot=%d, QS3score=%s, and newscore=%s.',
howgot,prtsc(QS3score),prtsc(newscore)))
listhbonds(howgot,0,2) -- send howgot into listhbonds
-- and save result into quicksave slot 2
-- recentbest.Save( ) -- not needed? 7/9/2023 jeff101
local gain = newscore - highscore
--
-- 2/18/2024 can you pool below few messages into just one line?
print('gain: '..prtsc(gain))
highscore = newscore
print("++ High score = "..prtsc(highscore))
listhbonds(0,0,0)
print('Got High score above using banded wiggles with:')
print(tstr..'.')
-- Above line lists the CI values used for banded wiggle stage.
-- 2/18/2024 can you pool above few messages into just one line?
--
setendnote() -- puts final Note in segment SegForNote
end -- if newscore
selection.DeselectAll()
print('')
end -- random_bands()
-- ----------------------------------------------------------
function LigandDock()
highscore = current.GetScore( )
local revertscore = current.GetScore( )
local i
for i = 1, Cycles do
Levy = 0
revertscore = current.GetScore( )
if revertscore > highscore then
-- 2/17/2024 how is this even possible?
-- if it happens, let it be, and skip the following.
else -- revertscore <= highscore
-- RevertCycles == 0 should skip the following every cycle i
if RevertCycles ~= 0 then
if i % RevertCycles == 0 then
if RevertFac <= 0 then
print(' -- reverting to best pose (quicksave slot 2 w/howgot='
..QSdat[2].howgot..') --')
save.Quickload(2)
else -- RevertFac > 0
local qs2use=pickqs()
print(' -- reverting to quicksave slot '..qs2use..
' w/howgot='..QSdat[qs2use].howgot..' --')
save.Quickload(qs2use)
-- 2/17/2024 pickqs() above lets RevertFac > 0
-- load quicksaves besides quicksave #2.
end -- if RevertFac
end -- if i
end -- if RevertCycles
-- LevyCycles == 0 should skip the following every cycle i
if LevyCycles ~= 0 then
if i % LevyCycles == 0 then
print(' -- Implementing Levy flight --')
Levy = math.random( minLevy, maxLevy )
-- Levy should be an integer from minLevy=500 to maxLevy=800 now
end -- if i
end -- if LevyCycles
end -- if revertscore
print('Cycle '..i..' of '..Cycles..' at '..os.date()..':')
random_bands()
end -- for i
end -- LigandDock()
-- ----------------------------------------------------------
function pickqs()
local scores,howgots,ngot,qsnums,terms,sums,j,num,pickj,k,isrepeat
-- Below loads quicksave scores & their slot #'s into scores & qsnums.
-- Make sure quicksave slot 2 is included. Include slots 1 3 as well? ok.
-- How about slot 4? no.
-- Should it prefer quicksaves made at the end of an unbanded wiggle stage
-- over ones from the end of a banded wiggle or the middle of an unbanded wiggle?
-- This means to check howgot for each quicksave, but as of 3/8/2024,
-- howgot is not recorded for quicksaves 1-3.
-- Should it only have one quicksave per score?
-- Some quicksaves like 1-3 often have the same score as other quicksaves.
-- Should the repeated scores be excluded? Keep them listed, but
-- set their terms[j] values to 0. Assume that identical scores
-- means the exact same protein + ligand pose.
scores={}
howgots={}
qsnums={}
setendnote()
save.Quicksave(4) -- put present pose in storage
for j=1,3 do
save.Quickload(j) -- load quicksave slot j
scores[j]=current.GetScore() -- get score for quicksave slot j
howgots[j]=QSdat[j].howgot
qsnums[j]=j
end -- for j
save.Quickload(4) -- get present pose back from storage
ngot=3
for j=QSlo,QSInUse do
ngot=ngot+1
scores[ngot]=QSdat[j].score
howgots[ngot]=QSdat[j].howgot
qsnums[ngot]=j
end -- for j
--
-- Below finds a Boltzmann-like term terms[j] for each scores[j].
-- scores[j] are like -1*energy = -E[j], so hi scores are for lowest energy.
-- Using scores[j] instead of the energy E[j] lets you skip the - sign
-- in the terms[j]=exp(-E[j]/kT) expressions below.
-- If scores[j] repeats scores[k], where k<j, set terms[j] to 0.
-- 3/8/2024 Would it be better to set terms[j]=0 for the lower of j,k in a pair?
-- 3/8/2024 Right now, it sets terms[j]=0 for the higher of j,k in a pair.
-- 3/8/2024 Should it set terms[j]=0 for ones w/howgots[j] < 3 ? yes, if j>1.
terms={}
for j=1,ngot do
isrepeat=0
if j > 1 then
if howgots[j] < 3 then
isrepeat=1
end -- if howgots
end -- if j
for k=1,j-1 do -- keeps k<j
if isrepeat==0 then
if scores[j]==scores[k] then
isrepeat=1
end -- if scores
end -- if isrepeat
end -- for k
if isrepeat==1 then
terms[j]=0.0
else
terms[j]=math.exp(scores[j]/RevertFac) -- RevertFac > 0 should hold
end -- if isrepeat
end -- for j
--
-- Below totals all values in terms to get the values in sums.
sums={}
sums[0]=0
for j=1,ngot do
sums[j]=sums[j-1]+terms[j]
end -- for j
-- I think sums[ngot] is like Z in the partition function
-- so the probability of getting quicksave j is
-- terms[j]/Z = terms[j]/sums[ngot].
-- Factors common to all terms[j] (like constants in the
-- conversion from energy E[j] to scores[j]) cancel out when
-- divided by sums[ngot].
--
-- Below num & sums are used to pick among the quicksaves,
-- making higher-scoring quicksaves more likely than
-- lower-scoring ones.
-- Will below work if some terms[j]==0 so that
-- some sums[j] are identical to each other?
-- I think it will work because pickj gets
-- changed from 0 to a positive integer
-- only for the first time a particular
-- value of sums[j] is used.
num=sums[ngot]*math.random()
-- math.random() gives real # form [0,1) or >=0 but < 1.
pickj=0
for j=1,ngot do
if pickj==0 then
if num>=sums[j-1] then
if num<sums[j] then
pickj=j
end -- if num
end -- if num
end -- if pickj
end -- for j
if pickj==0 then
print('pickj==0 error in pickqs()')
end -- if pickj
return qsnums[pickj] -- Returns the # of the quicksave to load.
-- If RevertFac is small, will it usually pick quicksave #2?
-- Will it pick other quicksaves with the same pose and score
-- as quicksave #2 instead?
end -- pickqs()
-------------------------------------------------------------------------
-- Below were adapted from Loop_rebuild_9.0f1.txt Feb 12 2023 1116pm code
-- to list atoms allowed for banded wiggles the orig way & the flip way.
-- below assumes templist is AtomsInLigand long
-- and contains just 0's & 1's. It totals the
-- # of 1's in templist.
function countones(templist)
local i
local tot=0
for i=1,AtomsInLigand do
tot=tot+templist[i]
end -- for i
return tot
end -- countones()
-- http://www.lua.org/manual/5.2/manual.html#6.4
-- helped make this function 11/16/17
function getlist(liststr)
local newlist={}
local i,ilo,ihi,idir,substr
for i=1,AtomsInLigand do
newlist[i]=0
end -- for i
-- below reads from liststr a series of substr's
-- where each substr contains an integer
-- followed by one or more '-' signs
-- followed by an integer
for substr in string.gmatch(liststr,"(%d+%-+%d+)") do
-- substr includes one or more '-' characters in a row
-- below gets ilo & ihi from substr
ilo=string.gsub(substr,"(%d+)%-+(%d+)","%1")+0
ihi=string.gsub(substr,"(%d+)%-+(%d+)","%2")+0
-- above gets ilo & ihi from substr
idir=1 -- the increment to use from ilo to ihi
if ilo>ihi then
idir= -1
end -- if ilo
for i=ilo,ihi,idir do -- for i=ilo to ihi step idir
if i>=1 and i<=AtomsInLigand then
newlist[i]=1
end -- if i
end -- for i
end -- for substr
-- below reads from liststr a series of substr's
-- where each substr contains an integer
for substr in string.gmatch(liststr,"(%d+)") do
i=substr+0 -- converts substr into the number i
if i>=1 and i<=AtomsInLigand then
newlist[i]=1
end -- if i
end -- for substr
return newlist
end -- getlist()
-- http://www.lua.org/manual/5.2/manual.html#6.4
-- helped make this function 11/16/17
function printlist(templist)
local outstr=''
local i,gotone,strlen
gotone=0
for i=1,AtomsInLigand do
if templist[i]==1 then
gotone=gotone+1
if gotone==1 then
outstr=(outstr..i)
end -- if gotone
else -- templist[i]==0
if gotone>1 then
outstr=(outstr..'-'..(i-1)..' ')
elseif gotone==1 then
outstr=(outstr..' ')
end -- if gotone
gotone=0
end -- if templist
end -- for i
if gotone>1 then
outstr=(outstr..'-'..AtomsInLigand)
end -- if gotone
strlen=string.len(outstr)
if string.sub(outstr,strlen,strlen)==' ' then
-- if outstr ends with a blank space,
-- then remove that blank space
outstr=string.sub(outstr,1,strlen-1)
end
return outstr
end -- printlist()
-- Above were adapted from Loop_rebuild_9.0f1.txt Feb 12 2023 1116pm code
-- to list atoms allowed for banded wiggles the orig way & the flip way.
-- ----------------------------------------------------------
function Prelim()
local kk,tstr
print( ReVersion ..' on '..ui.GetPlatform() )
print("Levy flight version, with random distance")
print('Starting pose for run in Quicksave 1')
print('Most-wanted pose for run in Quicksave 2')
print('Top-scoring pose for run in Quicksave 3')
print('Temporary storage for run in Quicksave 4')
print('Key poses for run in Quicksaves '..QSlo..' to '..QShi)
QSInUSe = QSlo-1
print('')
print(string.format('Puzzle %s (%s)',
puzzle.GetPuzzleID(),puzzle.GetName() ))
print("Protein length: "..structure.GetCount() )
print('# of Ligand Atoms: '..AtomsInLigand)
print(string.format('# of User-Supplied Bands: %d (%d enabled, %d disabled)',
OrigBands,OrigBandsOn,OrigBandsOff))
print('')
print('Inputs from Menu 1:')
print('Total # of Cycles to do: '..Cycles)
print('# of Cycles per Revert: '..RevertCycles..' (0 means never Revert)')
print('Revert Factor: '..RevertFac..' (0 always reverts to top score, ')
print(' more positive values can revert to lower scores from quicksaves)')
print('Flip Probability: '..FlipProb..' (0 means never Flip)')
print('Flip Angle Range: '..getmssg(FlipAngMin,FlipAngMax)..' in degrees')
print('Flip Hub Probability: '..FlipHubProb..
' (0 never uses a Flip Hub, 1 always uses a Flip Hub)')
print('# of Cycles per Levy: '..LevyCycles..' (0 means never use Levy)')
print('Levy Step: '..getmssg(minLevy,maxLevy)..' ('..getmssg(minLevy/1000,maxLevy/1000)..' Angstroms)')
print('# of Wiggles per Cycle with Bands: '..WigWithBands)
print('# of Wiggles per Cycle without Bands: '..Wiggles)
tstr={'wiggle ligand only','wiggle ligand and protein'}
print('Banded Wiggle Option: '..BandedWigOpt..' ('..tstr[BandedWigOpt]..')')
print('Unbanded Wiggle Option: '..UnbandedWigOpt..' ('..tstr[UnbandedWigOpt]..')')
print('# of Bands from Ligand to Space: '..BandsToLigand)
print('Band Start Lengths: '..getmssg(minDist,maxDist)..' ('..getmssg(minDist/1000,maxDist/1000)..' Angstroms)')
print('Band Goal Lengths: '..getmssg(minGoal,maxGoal)..' ('..getmssg(minGoal/1000,maxGoal/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('UseRecentBest: '..UseRecentBest)
print('Ligand Name for output files: ('..LigName..') without the outermost ()')
print('Random Seed: '..RNDSD)
print('Verbosity: '..Verbosity)
print('')
print('Inputs from Menu 2:')
--
-- are next 2 inputs used any more? 7/12/2023
print('Maximum CI recipe will use: '..userSetMaxCI )
print('CI for Banded Wiggles: '..WigBandCI )
-- are prev 2 inputs used any more? 7/12/2023
--
print('atoms4orig_list = \"'..atoms4orig_str..'\" = \"'..printlist(atoms4orig_list)..'\"')
print('atoms4flip_list = \"'..atoms4flip_str..'\" = \"'..printlist(atoms4flip_list)..'\"')
print('Below are Importances for banded wiggles:')
for kk=1,7 do
print(ciabbr[kk]..': '..getmssg(cilos[kk],cihis[kk]))
end -- for kk
print('CI bias: '..cibias..' (0 is off, 1 favors low CI values)')
print('Quicksave slots for key results: '..getmssg(QSlo,QShi))
print('')
print("Recipe start score: "..prtsc(current.GetEnergyScore()))
listhbonds(0,0,0)
print('')
math.randomseed(RNDSD)
math.random( ); math.random( ); math.random( ) --See stackoverflow
end -- Prelim()
function DialogOptions()
local dlog, dialogresult, mssg, kk, tstr, tlen
local gotinput=0 -- 0 means input isn't right
-- 1 would mean input is fine
while gotinput==0 do
dlog = dialog.CreateDialog(ReVersion..' Menu 1')
dlog.iters = dialog.AddTextbox("Total Cycles",Cycles)
dlog.label0 = dialog.AddLabel('Total Cycles = 1 integer 1 to 10000')
dlog.revcyc = dialog.AddSlider("Cycles/Revert",RevertCycles,0,10,0)
dlog.revfac = dialog.AddTextbox('Revert Factor',RevertFac)
dlog.label1 = dialog.AddLabel('Revert Factor = real # 0 to 20000')
dlog.flipprob = dialog.AddTextbox('Flip Probability',FlipProb)
dlog.label2 = dialog.AddLabel('Flip Probability = real # 0 to 1')
mssg=getmssg(FlipAngMin,FlipAngMax)
dlog.flipangs = dialog.AddTextbox('Flip Angle Range',mssg)
dlog.label3 = dialog.AddLabel('Flip Angle Range = 1-2 reals 0 to 180 degrees')
dlog.fliphubprob = dialog.AddTextbox('Flip Hub Prob',FlipHubProb)
dlog.label4 = dialog.AddLabel('Flip Hub Prob: 0 no hub, 1 flip uses a hub')
dlog.levycyc = dialog.AddSlider("Cycles/Levy",LevyCycles,-10,10,0)
mssg=getmssg(minLevy,maxLevy)
dlog.levystep = dialog.AddTextbox("Levy Step Range",mssg)
dlog.spacer1 = dialog.AddLabel('Levy Step Range = 1-2 integers 0 to 1000')
dlog.spacer1x = dialog.AddLabel('Levy Step & Band Lengths are in 0.001 Angstroms')
dlog.wwb = dialog.AddSlider("Banded Wiggles",WigWithBands,1,500,0)
dlog.wigls = dialog.AddSlider("Free Wiggles",Wiggles,1,500,0)
dlog.bandwigopt = dialog.AddSlider('Banded Option',BandedWigOpt,1,2,0)
dlog.freewigopt = dialog.AddSlider('Free Option',UnbandedWigOpt,1,2,0)
dlog.spacer1a = dialog.AddLabel('Option 1 wiggles ligand only, 2 wiggles protein too')
dlog.bnd2lig = dialog.AddSlider('# Ligand Bands',BandsToLigand,1,AtomsInLigand,0)
mssg=getmssg(minDist,maxDist)
dlog.startlens = dialog.AddTextbox('Band Start Lengths',mssg)
dlog.spacer2 = dialog.AddLabel('Band Start Lengths = 1-2 integers 5 to 150000')
mssg=getmssg(minGoal,maxGoal)
dlog.goallens = dialog.AddTextbox('Band Goal Lengths',mssg)
dlog.spacer3 = dialog.AddLabel('Band Goal Lengths = 1-2 integers 0 to 3500')
dlog.bandstren = dialog.AddSlider('Band Strength',BandStrength,0.1,5,1)
dlog.rcntbst = dialog.AddSlider('UseRecentBest',UseRecentBest,0,3,0)
dlog.LigName = dialog.AddTextbox('Ligand Name',LigName)
dlog.randseed = dialog.AddTextbox('Random Seed',RNDSD)
dlog.spacer4 = dialog.AddLabel('Random Seed = 1 integer 10000000 to 4294967296')
dlog.vrb = dialog.AddSlider('Verbosity',Verbosity,0,3,0)
dlog.ok = dialog.AddButton("OK",1)
dlog.cancel = dialog.AddButton("Cancel",0)
dialogresult = dialog.Show(dlog)
gotinput=1 -- 1 means input is fine
-- 0 would mean input isn't fine
--
if dialogresult > 0 then
-- below checks the inputs
Cycles,gotinput = readtextbox(dlog.iters.value,1,10000,1,gotinput)
-- in the above, the 1 before gotinput means to expect an integer as input
RevertCycles = dlog.revcyc.value
RevertFac,gotinput = readtextbox(dlog.revfac.value,0,20000,0,gotinput)
FlipProb,gotinput = readtextbox(dlog.flipprob.value,0,1,0,gotinput)
FlipAngMin,FlipAngMax,gotinput = readtextbox2(dlog.flipangs.value,0,180,0,gotinput)
FlipHubProb = readtextbox(dlog.fliphubprob.value,0,1,0,gotinput)
-- in the above, the 0 before gotinput means to expect real #'s as input
LevyCycles = dlog.levycyc.value
minLevy,maxLevy,gotinput = readtextbox2(dlog.levystep.value,0,1000,1,gotinput)
-- in the above, the 1 before gotinput means to expect integers as input
WigWithBands = dlog.wwb.value
Wiggles = dlog.wigls.value
BandedWigOpt = dlog.bandwigopt.value
UnbandedWigOpt = dlog.freewigopt.value
BandsToLigand = dlog.bnd2lig.value
minDist,maxDist,gotinput = readtextbox2(dlog.startlens.value,5,150000,1,gotinput)
minGoal,maxGoal,gotinput = readtextbox2(dlog.goallens.value,0,3500,1,gotinput)
-- in the above, the 1's before gotinput mean to expect integers as input
BandStrength = dlog.bandstren.value
UseRecentBest = dlog.rcntbst.value
LigName = purgeblanks(dlog.LigName.value)
RNDSD,gotinput = readtextbox(dlog.randseed.value,10000000,4294967296,1,gotinput)
Verbosity = dlog.vrb.value
if gotinput ~= 1 then
print('Something was input wrong. Please try again.')
end -- if gotinput
end -- if dialogresult
--
end -- while gotinput
--
if dialogresult > 0 then
gotinput=0 -- 0 means input isn't right
-- 1 would mean input is fine
while gotinput==0 do
dlog = dialog.CreateDialog(ReVersion..' Menu 2')
--
-- are next 2 inputs used any more? 7/12/2023
dlog.usermaxci = dialog.AddSlider("Maximum CI",userSetMaxCI,0.00,1,2)
dlog.wbci = dialog.AddSlider("Wiggle Bands CI",WigBandCI,0.00,1,2)
-- are prev 2 inputs used any more? 7/12/2023
--
dlog.alabel1=dialog.AddLabel('Below are allowed atoms for banded wiggles:')
dlog.atoms4orig = dialog.AddTextbox ( "Orig Way" , atoms4orig_str )
dlog.atoms4flip = dialog.AddTextbox ( "Flip Way" , atoms4flip_str )
dlog.alabel2 = dialog.AddLabel('List the allowed atoms like below:')
dlog.alabel3 = dialog.AddLabel('1-3,6,19-30,45-62 or 1-3 6 19-30 45-62')
dlog.label0=dialog.AddLabel('Below are Importances for banded wiggles:')
for kk=1,7 do
dlog['getci'..kk]=dialog.AddTextbox(ciabbr2[kk],getmssg(cilos[kk],cihis[kk]))
tstr=string.format('%s = 1-2 decimals %5.3f to %5.3f',ciabbr[kk],cimins[kk],cimaxs[kk])
dlog['label'..kk]=dialog.AddLabel(tstr)
end -- for kk
dlog.cibias = dialog.AddSlider("CI bias",cibias,0,1,0)
dlog.label8 = dialog.AddLabel('CI bias = 0 makes all CI equally likely')
dlog.label9 = dialog.AddLabel('CI bias = 1 favors low CI values')
dlog.getqs = dialog.AddTextbox('QS slots',getmssg(QSmin,QSmax))
dlog.qslabel = dialog.AddLabel('QS slots will store key results')
dlog.ok = dialog.AddButton("OK",1)
dlog.cancel = dialog.AddButton("Cancel",0)
dialogresult = dialog.Show(dlog)
gotinput=1 -- 1 means input is fine
-- 0 would mean input isn't fine
--
if dialogresult > 0 then
-- below checks the inputs
userSetMaxCI = dlog.usermaxci.value
WigBandCI = dlog.wbci.value
atoms4orig_list = getlist(dlog.atoms4orig.value)
atoms4orig_str = printlist(atoms4orig_list)
tlen=countones(atoms4orig_list)
if tlen<BandsToLigand then
gotinput=0
print('atoms4orig has only '..tlen..' of '..BandsToLigand..' elements.')
end -- if tlen
atoms4flip_list = getlist(dlog.atoms4flip.value)
atoms4flip_str = printlist(atoms4flip_list)
tlen=countones(atoms4flip_list)
if tlen<2 then
gotinput=0
print('atoms4flip has only '..tlen..' of 2 elements.')
end -- if tlen
for kk=1,7 do
cilos[kk],cihis[kk],gotinput =
readtextbox2(dlog['getci'..kk].value,cimins[kk],cimaxs[kk],0,gotinput)
-- 0 before gotinput means it expects decimal numbers like 0.832
end -- for kk
cibias = dlog.cibias.value
QSlo,QShi,gotinput=readtextbox2(dlog.getqs.value,QSmin,QSmax,1,gotinput)
-- 1 before gotinput means to expect integers
if gotinput ~= 1 then
print('Something was input wrong. Please try again.')
end -- if gotinput
end -- if dialogresult
--
end -- while gotinput
end -- if dialogresult
--
return dialogresult
end -- DialogOptions()
-- below makes a string mssg listing vlo to vhi.
-- if vlo==vhi, mssg lists vlo only.
function getmssg(vlo,vhi)
local mssg
if vlo<vhi then
mssg=(vlo..'-'..vhi)
elseif vhi<vlo then
mssg=(vhi..'-'..vlo)
else
mssg=vlo
end
return mssg
end -- getmssg()
-- below converts the string vtxt to a number val from vmin to vmax.
-- if intflag==1, val is an integer.
-- if intflag==0, val is a real #.
-- gotinput1 is input status before this function (0 or 1).
-- gotinput2 is input status after this function (0 or 1).
function readtextbox(vtxt,vmin,vmax,intflag,gotinput1)
local valo,val,gotinput2,ttxt
--
gotinput2=0 -- 0 means input isn't right
-- 1 would mean input is fine
--
-- below converts string vtxt to the number val
valo=tonumber(vtxt)
if valo==nil then -- vtxt not a valid number
val=vmin
else -- vtxt gave a valid number valo
if intflag==1 then
-- below rounds valo to the integer val
val = valo - valo % 1
else -- let val be a real #
val = valo
end -- if intflag
-- below puts val in range vmin to vmax
if val<vmin then
val=vmin
elseif val>vmax then
val=vmax
end -- if val
if val==valo then -- vtxt,valo were fine
gotinput2=1
end -- if val
end -- if valo
if gotinput2==0 then -- some input is bad
ttxt='real'
if intflag==1 then
ttxt='integer'
end -- if intflag
print('Error: Got '..vtxt..' when wanted an '..ttxt..' from '..vmin..' to '..vmax..'.')
end -- if gotinput2
--
-- below combines latest input status
-- with previous input status
gotinput2=gotinput2*gotinput1
--
return val,gotinput2
end -- readtextbox()
-- if intflag==1, below converts the string vtxt
-- to 2 integers vals[1] <= vals[2] from vmin to vmax.
-- vtxt should accept the formats
-- '432-546', '546-432', '489-489', or '489'.
--
-- if intflag==0, below accepts positive decimal
-- real #'s in vtxt instead of integers.
--
-- gotinput1 is input status before this function (0 or 1).
-- gotinput2 is input status after this function (0 or 1).
function readtextbox2(vtxt,vmin,vmax,intflag,gotinput1)
local ngot,vals,vtmp,kk,gotinput2,ttxt
--
gotinput2=1 -- 1 means input is fine
-- 0 would mean input isn't right
ngot=0
vals={}
--
-- print('vtxt='..vtxt)
if intflag==1 then
-- next line sends each integer in vtxt to vtmp
-- and ignores all other characters in vtxt.
for vtmp in string.gmatch(vtxt,"%d+") do
-- in %d+, %d looks for the digits 0-9, and
-- + looks for the longest possible group
-- of 1 or more of them in a row.
ngot=ngot+1
-- print('vals['..ngot..']='..vtmp)
vals[ngot]=tonumber(vtmp)
end -- for vtmp
else -- do same as above but for positive decimal real #'s like 2.999
for vtmp in string.gmatch(vtxt,"%d*%.?%d+") do
-- %d* looks for the longest possible group
-- of 0 or more digits 0-9 in a row.
-- %.? looks for 0 or 1 decimal point.
-- %d+ looks for the longest possible group
-- of 1 or more digits 0-9 in a row.
-- I think the above will skip '.' that are not followed by
-- 1 or more digits 0-9.
-- I think it will read '3', '.3', '3.3', '333.333' as desired.
-- Will it read '3.' as '3' ?
-- See http://www.lua.org/manual/5.1/manual.html#5.4.1 about patterns.
ngot=ngot+1
-- print('vals['..ngot..']='..vtmp)
vals[ngot]=tonumber(vtmp)
end -- for vtmp
end -- if intflag
--
if ngot==1 then
vals[2]=vals[1]
ngot=2
elseif ngot==2 then
if vals[1]>vals[2] then
vals[1],vals[2]=vals[2],vals[1]
end -- if vals[1]
else
gotinput2=0 -- the input is bad
vals[1]=vmin
vals[2]=vmax
end -- if ngot
if gotinput2==1 then
for kk=1,ngot do
if vals[kk]<vmin then
vals[kk]=vmin
elseif vals[kk]>vmax then
vals[kk]=vmax
end -- if vals[kk]
end -- for kk
else -- some input is bad
ttxt='reals'
if intflag==1 then
ttxt='integers'
end -- if intflag
print('Error: Got '..vtxt..' when wanted 1-2 '..ttxt..' from '..vmin..' to '..vmax..'.')
end -- if gotinput2
-- print('gotinput2='..gotinput2..' vals[1]='..vals[1]..' vals[2]='..vals[2])
--
-- below combines latest input status
-- with previous input status
gotinput2=gotinput2*gotinput1
--
return vals[1],vals[2],gotinput2
end -- readtextbox2()
-- Below takes the string txtin and removes blanks at start and end of it.
-- Then it returns the result as txtout.
function purgeblanks(txtin)
local txtout,tstr
tstr='purgeblanks('..txtin..') gives ('
txtout=string.gsub(txtin,'%s+','',1)
-- above replaces the first group of 1 or more
-- successive space characters in txtin with ''
txtout=string.reverse(string.gsub(string.reverse(txtout),'%s+','',1))
-- above replaces the last group of 1 or more
-- successive space characters in txtout with ''
print(tstr..txtout..') without the outermost ().')
return txtout
end -- purgeblanks()
function Main()
-- Below tests sortit():
local x={70,-50,55,40,90,-30,0}
local sx,indx,dim,ii
sx,indx=sortit(x)
dim=#x
print('Below tests sx,indx=sortit(x):')
print('ii x[ii] indx[ii] sx[ii]:')
for ii=1,dim do
print(string.format('%d %5d %8d %6d',ii,x[ii],indx[ii],sx[ii]))
end -- for ii
print('Above tests sortit().')
print('Ideally sx[ii]=x[indx[ii]] above for ii=1 to '..dim..'.')
print(' ')
-- Above tests sortit().
--
setnoteseg() -- finds SegForNote
setendnote() -- puts final Note in segment SegForNote
listhbonds(1,0,1) -- save OG pose for backup in qs slot 1 with howgot=1
-- when in doubt, use howgot=1
for kk=2,3 do
copyslot1to2(1,kk) -- copy slot 1 into slots 2 and 3
end -- for kk
-- qs slot 2 will be most wanted pose in run
-- qs slot 3 will be top-scoring pose in run
RNDSD=SeedRandom()
HasLigandCheck( ) -- this finds AtomsInLigand
if HASLIGAND == false then
print('This is NOT a LIGAND puzzle ... exiting')
return
end
--
local templist
templist = string.format('1-%d',AtomsInLigand)
atoms4orig_str = templist
atoms4orig_list = getlist(templist)
atoms4flip_str = templist
atoms4flip_list = getlist(templist)
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 )
highscore=current.GetScore()
setendnote() -- puts final Note in segment SegForNote
print('Restoring best score: '..prtsc(highscore))
listhbonds(0,0,0)
undo.SetUndo( true ) -- what does this do?
-- 2/18/2024 It lets changes be written to the undo graph.
-- https://foldit.fandom.com/wiki/Foldit_Lua_Function_undo.SetUndo
-- says using undo.SetUndo(false) blocks changes from writing
-- to the undo graph. This would let a single undo undo all
-- changes that the recipe makes at once.
-- is below line still needed? 7/12/2023
behavior.SetClashImportance( userSetMaxCI ) -- return to user setting
-- is above line still needed? 7/12/2023
local kk
for kk=1,7 do
behavior['Set'..cinams[kk]..'Importance'](ciolds[kk])
end -- for kk
listallqs()
end -- Main()
-- Main()
-- comment above line out once recipe ~behaves
-- comment below line out during testing
xpcall( Main, Cleanup )
--Ligand Docker v2.5
