Recipe: Ligand Docker v2.5b

created by jeff101

Profile

• Name: Ligand Docker v2.5b
• ID: 108720
• Shared with: Public
• Parent: None
• Children: None
• Created on: May 29, 2024 at 04:37 AM UTC
• Updated on: May 30, 2024 at 00:29 AM UTC
• Description:

  May 27 2024 912pm *v2.5b0e.txt code

Code

-- ------------------------------------
-- Ligand Docker v2.5b
-- Credits:  Raven_pl,Susume,LociOiling,lynnai,
--           robgee,jeff101
-- last updated May 27 2024 912pm
-- ------------------------------------
-- 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.

-- 5/11/2024 for v2.5b, give user options to 
-- wiggle ligand only, ligand + sidechains, 
-- or ligand + sidechains + backbone.
-- Let there be a different chance in each
-- wiggle stage that it does one of the 
-- options above.

-- ------------------------------------
-- 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.5b"
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.
BandedWigOptVals={9,3,1} -- {9,3,1} here means for every 13 wiggles,
                    -- 9 will do WigOpt=1, 3 will do WigOpt=2,
                    -- & 1 will do WigOpt=3.
                    --
                    -- {8,0,2} here would mean for every 10 wiggles, 
                    -- 8 will do WigOpt=1, 0 will do WigOpt=2, 
                    -- & 2 will do WigOpt=3.
                    --
                    --1 wiggles ligand only, 2 wiggles ligand & sidechains,
                    --3 wiggles ligand & sidechains & backbone.
                    --
                    --default was like {1,0,0} (always WigOpt=1) before 7/12/2023.
                    --default was like {0,0,1} (always WigOpt=3) before 5/11/2024.
UnbandedWigOptVals={7,7,7} -- Put {7,7,7} here if want unbanded to  
                    -- use the same WigOpt's as banded in each cycle.
                    --
                    -- {8,2,1} here means for every 11 wiggles, 
                    -- 8 will do WigOpt=1, 2 will do WigOpt=2,
                    -- & 1 will do WigOpt=3.   
                    --
                    --1 wiggles ligand only, 2 wiggles ligand & sidechains,
                    --3 wiggles ligand & sidechains & backbone.
                    --
                    --default was like {0,0,1} (always WigOpt=3) before 7/12/2023.
                    --default was like {0,0,1} (always WigOpt=3) before 5/11/2024.
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 [0,1), 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()

-- Below preps for a 4-band flip. Its inputs are len, the distance between
-- the 2 atoms to flip (atoms 1 and 2, which start at points A and B), 
-- th,ph, the angles in radians for B when A is at the origin, 
-- and alp, the angle in radians between A and A' as well as B and B'.
-- Overall, A A' B B' lie in the same plane and are the corners of a rectangle.
function prep4bandflip(len,th,ph,alp)
local x1y1z1s={} -- x1 y1 z1 coordinates for A B A' B' 
local x2y2z2s -- x2 y2 z2 coordinates for A B A' B' 
local xyzs -- x y z coordinates for A B A' B'
local rhothphs -- rho values and theta, phi values (angles in radians) for A B A' B'
local hlen=len/2.0 -- half the length len
local gam=math.random()*2.0*math.pi -- a random angle in radians = [0,1)*2*pi 
   -- so gam >= 0 and gam < 2*pi holds.
local ca,sa,cg,sg,ct,st,cp,sp -- cosines and sines for alp, gam, th, and ph
local rot1t2={} -- a matrix to rotate coordinates from x1 y1 z1 to x2 y2 z2
local rot2txyz={} -- a matrix to rotate coordinates from x2 y2 z2 to x y z
ca=math.cos(alp)
sa=math.sin(alp)
cg=math.cos(gam)
sg=math.sin(gam)
ct=math.cos(th)
st=math.sin(th)
cp=math.cos(ph)
sp=math.sin(ph)
x1y1z1s[1]={       0, 0,    0       } -- A  = starting position for atom 1
x1y1z1s[2]={       0, 0,  len       } -- B  = starting position for atom 2
x1y1z1s[3]={ hlen*sa, 0, hlen*(1-ca)} -- A' = target/goal position for atom 1
x1y1z1s[4]={-hlen*sa, 0, hlen*(1+ca)} -- B' = target/goal position for atom 2
rot1t2[1]={cg,-sg,  0}
rot1t2[2]={sg, cg,  0} -- a matrix to rotate coordinates from x1 y1 z1 to x2 y2 z2
rot1t2[3]={ 0,  0,  1}
rot2txyz[1]={cp*ct, -sp, cp*st}
rot2txyz[2]={sp*ct,  cp, sp*st} -- a matrix to rotate coordinates from x2 y2 z2 to x y z
rot2txyz[3]={  -st,   0,    ct}
x2y2z2s=matbyvecs(rot1t2,x1y1z1s)
xyzs=matbyvecs(rot2txyz,x2y2z2s)
rhothphs=xyztspher(xyzs)
return rhothphs
end -- prep4bandflip()

-- Below finds the product of the matrix mat with the vectors in the rows of vecs.
-- It returns the results as vectors in the rows of res.
function matbyvecs(mat,vecs)
local matnr,matnc -- # of rows, columns for the matrix mat
local vecsnr,vecsnc -- # of rows, columns for the matrix vecs
local ii,jj,kk
local qflag=0 -- will quit if this gets reset to 1
local res={} -- matrix of results, should be vecsnr rows and vecsnc columns
matnr=#mat
matnc=#(mat[1])
vecsnr=#vecs
vecsnc=#(vecs[1])
qflag=0
if matnc~=vecsnc then
   print('Error in matbyvecs: matnc='..matnc..' differs from vecsnc='..vecsnc..'.')
   qflag=1
end -- if matnc
if matnr~=matnc then
   print('Error in matbyvecs: matnr='..matnr..' differs from matnc='..matnc..'.')
   qflag=1
end -- if matnr
if qflag==0 then
   for ii=1,vecsnr do
      res[ii]={}
      for jj=1,matnr do
         res[ii][jj]=0.0
         for kk=1,matnr do
            res[ii][jj]=res[ii][jj]+mat[jj][kk]*vecs[ii][kk]
         end -- for kk
      end -- for jj
   end -- for ii
   return res
else
   -- there's an error ... what to do?
end -- if qflag
end -- matbyvecs()

-- Below converts the x,y,z coordinates in each row of xyzs
-- into spherical rho,th,ph coordinates (th,ph are angles in radians) 
-- and returns them in the rows of res.
function xyztspher(xyzs)
local nrows,k,x,y,z,pi,rho,th,ph
local res={}
pi=math.pi
nrows=#xyzs
for k=1,nrows do
   x=xyzs[k][1]
   y=xyzs[k][2]
   z=xyzs[k][3]
   rho=math.sqrt(x^2 + y^2 + z^2)
   th=math.acos(z/rho) -- 0 to pi in radians
   ph=math.atan2(y,x) -- 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   
   res[k]={rho,th,ph}
end -- for k
return res
end -- xyztspher()

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 -------------

    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
    selection.DeselectAll () -- selects nothing
    selection.Select ( OG ) -- selects the ligand
    --structure.ShakeSidechainsSelected(witers)
    BandedWigOpt=pickopt(BandedWigOptVals)
    --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.
    elseif BandedWigOpt == 3 then -- wiggle ligand + sidechains + backbone
       structure.WiggleAll(witers, true, true) -- why not use this line? jeff101 7/9/2023
    else -- below should wiggle ligand + sidechains only -- unclear how to fix 5/11/2024.
       structure.WiggleAll(witers, false, true) 
       -- Above keeps protein bb fixed & blocks ligand xlation
       -- but lets protein sc move & ligand bend.
       -- Maybe can do better than this later.
    end -- if BandedWigOpt
    selection.DeselectAll()

    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

    selection.DeselectAll () -- selects nothing
    selection.Select ( OG ) -- selects the ligand
    --structure.ShakeSidechainsSelected( 2 )
    if all7s(UnbandedWigOptVals)==1 then -- if all values are 7, then reuse BandedWigOpt below
       UnbandedWigOpt=BandedWigOpt
       print('Reusing Banded Wiggle\'s option '..BandedWigOpt..' for Unbanded Wiggle.')
    else
       UnbandedWigOpt=pickopt(UnbandedWigOptVals)
    end -- if all7s()
    if UnbandedWigOpt == 3 then -- this was default before 7/12/2023
       structure.WiggleAll( Wiggles, true, true )
        -- above moves ligand thru space (1st true) and lets it bend (2nd true).
        -- above also moves protein backbone (1st true) and sidechains (2nd true).
    elseif UnbandedWigOpt == 1 then
       structure.WiggleSelected( Wiggles, true, true )
       -- above just moves ligand thru space (1st true) and lets it bend (2nd true).
    else -- wiggle ligand + sidechains only. -- unclear how to fix 5/11/2024.
       structure.WiggleAll( Wiggles, false, true ) 
       -- Above keeps protein bb fixed & blocks ligand xlation
       -- but lets protein sc move & ligand bend.
       -- Maybe can do better than this later.
    end -- if UnbandedWigOpt
    selection.DeselectAll()

    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()

-- below uses the 3 integers in vals
-- and a random # rnum to pick option 1 2 or 3
function pickopt(vals)
local opt,sums,j,rnum,tstr
sums={}
sums[0]=0
for j=1,3 do
   sums[j]=vals[j]+sums[j-1]
end -- for j
rnum=math.random()*sums[3] 
-- math.random() gives real # from [0,1) or >=0 but < 1
-- so rnum is [0,sums[3])    or >=0 but < sums[3]
opt=0
for j=1,3 do
   if opt==0 then
      if rnum>=sums[j-1] then
         if rnum<sums[j] then
            opt=j
         end -- if rnum
      end -- if rnum
   end -- if opt
end -- for j
if opt==0 then
   print('opt==0 error in pickopt()')
end -- if opt
tstr={'wiggle ligand only','wiggle ligand + sidechains','wiggle ligand + sidechains + backbone'}
print('rnum='..rnum..' so will do option '..opt..' ('..tstr[opt]..') next.')
return opt
end -- pickopt()
-- ----------------------------------------------------------
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 printvals(vals)
   local outstr=''
   local i
   for i=1,#vals do
      outstr=(outstr..vals[i])
      if i<#vals then
         outstr=(outstr..' ')
      end -- if i
   end -- for i
   return outstr
end -- printvals()

-- Below gives 1 if vals[1-3] are all 7's
-- and 0 otherwise.
function all7s(vals)
local res=0
if vals[1]==7 then
   if vals[2]==7 then
      if vals[3]==7 then
         res=1
      end -- if vals[3]
   end -- if vals[2]
end -- if vals[1]
return res
end -- all7s()

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 sidechains',
       'wiggle ligand, sidechains, and backbone'}
    print('Banded Wiggles ('..printvals(BandedWigOptVals)..'):')
    for kk=1,3 do
       print('   Do '..BandedWigOptVals[kk]..' with option '..kk..' ('..tstr[kk]..')')
    end -- for kk
    print('Unbanded Wiggles ('..printvals(UnbandedWigOptVals)..'):')
    if all7s(UnbandedWigOptVals)==1 then -- all 3 UnbandedWigOptVals are 7's
       print('   Will copy Unbanded Wiggle Options from Banded Wiggle Options')
    else
       for kk=1,3 do
          print('   Do '..UnbandedWigOptVals[kk]..' with option '..kk..' ('..tstr[kk]..')')
       end -- for kk
    end -- if all7s()
    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('')
    print('Inputs from Menu 2:')
    print('Ligand Name for output files: ('..LigName..') without the outermost ()')
    print('Random Seed: '..RNDSD)
    print('Verbosity: '..Verbosity)
    --
    -- 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('Above are Importances for banded wiggles')
    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, jj, 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)
    mssg=printvals(BandedWigOptVals)
    dlog.bandwigopt = dialog.AddTextbox('Banded Option',mssg)
    mssg=printvals(UnbandedWigOptVals)
    dlog.freewigopt = dialog.AddTextbox('Free Option',mssg)
    dlog.spacer1a = dialog.AddLabel('Above are likelihoods to pick Option 1, 2, or 3.')
    dlog.spacer1b = dialog.AddLabel('1 wiggles ligand only, 2 wiggles sidechains too,')
    dlog.spacer1c = dialog.AddLabel('and 3 wiggles ligand, sidechains, and backbone.')
    dlog.spacer1d = dialog.AddLabel('Use 7 7 7 to copy Banded Option to Free Option.')
    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.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
        BandedWigOptVals,gotinput = readtextbox3(dlog.bandwigopt.value,gotinput)
        UnbandedWigOptVals,gotinput = readtextbox3(dlog.freewigopt.value,gotinput)
        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
        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')
        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)
        --
        -- 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 jj=1,7 do
           dlog['getci'..jj]=dialog.AddTextbox(ciabbr2[jj],getmssg(cilos[jj],cihis[jj]))
           if jj==1 then
              for kk=1,2 do
                 if kk==1 then
                    mssg=ciabbr[kk]
                 else -- kk==2 here
                    mssg='The rest below'
                 end -- if kk
                 tstr=string.format('%s = 1-2 decimals %5.3f to %5.3f',
                    mssg,cimins[kk],cimaxs[kk])
                 dlog['label'..kk]=dialog.AddLabel(tstr) 
              end -- for kk
           end -- if jj
        end -- for jj
        dlog.label8 = dialog.AddLabel('Above are Importances for banded wiggles')
        dlog.cibias = dialog.AddSlider("CI bias",cibias,0,1,0)
        dlog.label9 = dialog.AddLabel('CI bias = 0 makes all CI equally likely')
        dlog.label10 = 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
            LigName = purgeblanks(dlog.LigName.value)
            RNDSD,gotinput = readtextbox(dlog.randseed.value,10000000,4294967296,1,gotinput)
            Verbosity = dlog.vrb.value
            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()

-- want below to convert vtxt into 3 positive integers.
function readtextbox3(vtxt,gotinput1)
local ngot,vals,vtmp,kk,gotinput2
--
gotinput2=1 -- 1 means input is fine
-- 0 would mean input isn't right
ngot=0
vals={}
--
print('vtxt='..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
   if ngot<=3 then
      print('vals['..ngot..']='..vtmp)
      vals[ngot]=tonumber(vtmp)
      if vals[ngot]<0 then
         gotinput2=0 -- the input is bad
         vals[ngot]=0
      end -- if vals
   end -- if ngot
end -- for vtmp
--
if ngot<3 then
   gotinput2=0 -- the input is bad
   for kk=ngot+1,3 do
      vals[kk]=0
   end -- for kk
end -- if ngot
--
if vals[1]==0 then 
   if vals[2]==0 then
      if vals[3]==0 then
         gotinput2=0 -- the input is bad
         vals[1]=1
      end -- if vals[3]
   end -- if vals[2]
end -- if vals[1]
--
print('gotinput2='..gotinput2..' vals[1]='..vals[1]..' vals[2]='..vals[2]..' vals[3]='..vals[3])
--
-- below combines latest input status 
-- with previous input status
gotinput2=gotinput2*gotinput1
--
return vals,gotinput2
end -- readtextbox3()

-- 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.5b

Comments

[jeff101]

For now, please run this recipe using its default values. It is set to run 10,000 cycles by default. I haven't determined yet when it is best to quit this recipe. Odds are you'll want to quit it before 10,000 cycles have finished. It should behave the same whether you wait for it to finish or quit it early.

By default, this recipe will use all 100 quicksave slots. You can change some input values to pick a smaller range of quicksave slots to use instead. When this recipe ends, it will output the contents of the quicksave slots into *ir_solution files and the undo graph. One input variable is the ligand name you want to appear in the *ir_solution files' titles. By default the ligand name is left blank.

This recipe varies all of the behavior sliders, so under "Behavior Options", check the "Recipe score modding" box before running it. This recipe also monitors hydrogen bonds between the protein and ligand, so under "View Options", check the "Show bonds (non-protein)" box.

The main new feature in v2.5b is to have 3 different ways to wiggle (1 = ligand only, 2 = ligand + sidechains, & 3 = ligand + sidechains + backbones). The default values are set to pick one of these ways at random each cycle. Right now, 9 of 13 cycles do way 1 (the fastest way), 3 of 13 do way 2, and 1 of 13 does way 3 (the slowest way). Right now, the 2nd wiggle option doesn't work as well as it could. It is limited by Foldit's LUA commands. It lets the sidechains wiggle and lets the ligand bend & twist, but it won't let the ligand translate. Meanwhile, wiggle options 1 and 3 let the ligand translate as well as bend & twist.

For puzzles like 2461, I run FreezeSelectSome1 (https://fold.it/recipes/102617) to freeze all the protein backbones (segments 1-512 in Puzzle 2461) before I run Ligand Docker v2.5b on the result.

Please Foldit e-mail me at jeff101 if you have any questions about or problems with this recipe. Thanks!

[jeff101]

In Puzzle 2464, it looks like freezing is not allowed, so to keep the backbones fixed, I switched from Ligand Docker's default values 9 3 1 in the previous comment to 9 3 0. Using 9 3 0 means that for every 12 cycles, 9 of 12 will use wiggle option 1 (ligand only, the fastest option), 3 of 12 will use wiggle option 2 (ligand + sidechains, the 2nd fastest option), and 0 of 12 will use wiggle option 3 (ligand + sidechains + backbones, the slowest option). Using 3 1 0 instead of 9 3 0 would act the same way: 9/12 (=3/4) of the cycles would do option 1, 3/12 (=1/4) of the cycles would do option 2, and no cycles would do option 3.