Recipe: Ligand Docker v2.5d.short

created by NinjaGreg

Profile

• Name: Ligand Docker v2.5d.short
• ID: 108948
• Shared with: Public
• Parent: Ligand Docker v2.5d
• Children:
  • Ligand Docker v2.5d.short2
• Created on: December 23, 2024 at 01:15 AM UTC
• Updated on: December 23, 2024 at 01:16 AM UTC
• Description:

  Dec 1 2024 147am *v2.5d0b.txt code

Code

-- ------------------------------------
-- Ligand Docker v2.5d
-- Credits:  Raven_pl,Susume,LociOiling,lynnai,
--           robgee,jeff101
-- last updated Dec 1 2024 147am
-- ------------------------------------
-- 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.

-- 5/31/2024 for v2.5c, make menus less tall.

-- 7/11/2024 for v2.5d, make 2 new ways to flip ligand atoms.
-- 12/1/2024 for v2.5d, let it run for non-ligand puzzles too.
--           It will assume final segment is the only ligand
--           in the puzzle, even if final segment is not a ligand.
-- ------------------------------------
-- 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.5d"
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
--  comment out next line 12/1/2024
--  while structure.GetSecondaryStructure(segCnt2)=="M" do 
       segCnt2 = segCnt2-1 
--  end -- while
--  comment out prev line 12/1/2024 
    -- 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,torstxt
--
-- 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
   -- print('ta')
      totbonus=filter.GetBonusTotal()
   -- print('totbonus='..totbonus)
      if totbonus>=0 then
         bonustxt=string.format('+%.1f',totbonus)
      else
         bonustxt=string.format('%.1f',totbonus)
      end -- if totbonus
   -- print('bonustxt='..bonustxt)
      -- see https://foldit.fandom.com/wiki/Foldit_Lua_Function_filter.GetBonus
      if structure.GetSecondaryStructure(nsegs)=="M" then -- check if ligand puzzle
         -- Standard ligand puzzles have just one ligand at segment # nsegs
         -- and have a Torsion Quality bonus (as of 12/1/2024).
         torstxt=string.format('torsion %.1f',filter.GetBonus('Torsion Quality'))
      else -- non-ligand puzzles have no torsion bonus
         torstxt='no torsion'
      end -- if ligand puzzle
   -- print('torstxt='..torstxt)
      print(
         string.format('For %s %s w/bonus %s and %s, ',
            LigName,prtsc(current.GetEnergyScore()),
            bonustxt,torstxt)..
         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:')
shortliststart = ngot
if ngot > 5 then
    ngot = 5
end
    
for ii=shortliststart,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()
-- -----------------------------------------
--
-- For flipping below, it would help a lot if Foldit would allow
-- bands to directly connect ligand's atom1 to ligand's atom2.
-- Then it would be easy to maintain the distance between atom1
-- and atom2 as you flip them (swap their positions or rotate them
-- around the midpoint between their starting positions).
--
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 preps for an 8-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.
-- 3 other useful points are M M1 & M2. M is the midpoint between points A & B.
-- M is also the midpoint between points A' & B'. M lies in the same plane as
-- A A' B B' (the x1 z1 plane in the x1 y1 z1 coordinate system).
-- M1, M, & M2  are on the y1 axis at y1 = -hlen, 0, & hlen respectively (hlen=len/2).
-- A, M, & B are along the z1 axis at z1 = -hlen, 0, & hlen respectively.
-- A, M, & B are along the z2 axis at z2 = -hlen, 0, & hlen respectively.
-- M1 M M2 all line in the x2 y2 plane.
function prep8bandflip(len,th,ph,alp)
local x1y1z1s={} -- x1 y1 z1 coordinates for M M1 M2 A B A' B' 
local x2y2z2s -- x2 y2 z2 coordinates for M M1 M2 A B A' B' 
local x3y3z3s -- x3 y3 z3 coordinates for M M1 M2 A B A' B' 
local xyzs -- x y z coordinates for M M1 M2 A B A' B'
local rhothphs -- rho values and theta, phi values (angles in radians) for M M1 M2 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 rot2t3={} -- a matrix to rotate coordinates from x2 y2 z2 to x3 y3 z3
local shift3xyz={} -- a vector to shift coordinates from x3 y3 z3 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   } -- M  = midpt between A & B = midpt between A' & B'
x1y1z1s[2]={       0, -hlen,     0   } -- M1 = endpt of vector perp to plane of M A B A' B'
x1y1z1s[3]={       0,  hlen,     0   } -- M2 = endpt of vector perp to plane of M A B A' B' 
x1y1z1s[4]={       0,     0, -hlen   } -- A  = starting position for atom 1
x1y1z1s[5]={       0,     0,  hlen   } -- B  = starting position for atom 2
x1y1z1s[6]={ hlen*sa,     0, -hlen*ca} -- A' = target/goal position for atom 1
x1y1z1s[7]={-hlen*sa,     0,  hlen*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} -- it keeps M at the origin.
rot2t3[1]={cp*ct, -sp, cp*st}
rot2t3[2]={sp*ct,  cp, sp*st} -- a matrix to rotate coordinates from x2 y2 z2 to x3 y3 z3.
rot2t3[3]={  -st,   0,    ct} -- it keeps M at the origin.
shift3txyz={hlen*cp*st, hlen*sp*st, hlen*ct} -- a vector to shift coordinates from
   -- x3 y3 z3 to x y z. it moves point A to the origin.
x2y2z2s=matbyvecs(rot1t2,x1y1z1s)
x3y3z3s=matbyvecs(rot2t3,x2y2z2s)
xyzs=vecplusvecs(shift3txyz,x3y3z3s)
rhothphs=xyztspher(xyzs)
return rhothphs
end -- prep8bandflip()

-- 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 finds the sum of the vector vec with the vectors in the rows of vecs.
-- It returns the results as vectors in the rows of res.
function vecplusvecs(vec,vecs)
local vecnc -- # of columns for the vector vec. vecnc should equal 3.
local vecsnr,vecsnc -- # of rows, columns for the matrix vecs.
   -- vecsnr is the # of points in vecs. vecsnc is the # of cols in each point's vector.
   -- vecsnc should equal 3.
local ii,jj
local qflag=0 -- will quit if this gets reset to 1
local res={} -- matrix of results, should be vecsnr rows and vecsnc columns
vecnc=#vec
vecsnr=#vecs
vecsnc=#(vecs[1])
qflag=0
if vecnc~=vecsnc then
   print('Error in vecplusvecs: vecnc='..vecnc..' differs from vecsnc='..vecsnc..'.')
   qflag=1
end -- if matnc
if qflag==0 then
   for ii=1,vecsnr do
      res[ii]={}
      for jj=1,vecsnc do
         res[ii][jj]=vec[jj]+vecs[ii][jj]
      end -- for jj
   end -- for ii
   return res
else
   -- there's an error ... what to do?
end -- if qflag
end -- vecplusvecs()

-- 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('')
    print('Inputs from Menu 2:')
    print('Band Start Lengths: '..getmssg(minDist,maxDist)..' ('..getmssg(minDist/1000,maxDist/1000)..' Angstroms)')
    print('Band Goal Lengths: '..getmssg(minGoal,maxGoal)..' ('..getmssg(minGoal/1000,maxGoal/1000)..' Angstroms)')
-- 8/10/2020 jeff101 thinks minDist & maxDist are given in 0.001 Angstroms
-- so that math.random(minDist,maxDist) gives an integer from minDist to maxDist
    print('Band Strength: '..BandStrength )
    print('UseRecentBest: '..UseRecentBest)
    print('Ligand Name for output files: ('..LigName..') without the outermost ()')
    print('Random Seed: '..RNDSD)
    print('Verbosity: '..Verbosity)
    --
    -- 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('')
    print('Inputs from Menu 3:')
    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')
    -- try to keep <= 25 dialog items between Menu line & OK line.
    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)
    -- try to keep <= 25 dialog items between Menu line & OK line.
    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
        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')
        -- try to keep <= 25 dialog items between Menu line & OK line.
        dlog.spacer2o = dialog.AddLabel('Band Lengths below are in 0.001 Angstroms')
        mssg=getmssg(minDist,maxDist)
        dlog.startlens = dialog.AddTextbox('Band Start Lengths',mssg)
        dlog.spacer2 = dialog.AddLabel('Band Start Lengths = 1-2 integers 5 to 150000')
        mssg=getmssg(minGoal,maxGoal)
        dlog.goallens = dialog.AddTextbox('Band Goal Lengths',mssg)
        dlog.spacer3 = dialog.AddLabel('Band Goal Lengths = 1-2 integers 0 to 3500')
        dlog.bandstren = dialog.AddSlider('Band Strength',BandStrength,0.1,5,1)
        dlog.rcntbst = dialog.AddSlider('UseRecentBest',UseRecentBest,0,3,0)
        dlog.LigName = dialog.AddTextbox('Ligand Name',LigName)
        dlog.randseed = dialog.AddTextbox('Random Seed',RNDSD)
        dlog.spacer4 = dialog.AddLabel('Random Seed = 1 integer 10000000 to 4294967296')
        dlog.vrb = dialog.AddSlider('Verbosity',Verbosity,0,3,0)
        --
        -- 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')
        -- try to keep <= 25 dialog items between Menu line & OK line.
        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
            minDist,maxDist,gotinput = readtextbox2(dlog.startlens.value,5,150000,1,gotinput)
            minGoal,maxGoal,gotinput = readtextbox2(dlog.goallens.value,0,3500,1,gotinput)
            -- in the above, the 1's before gotinput mean to expect integers as input
            BandStrength = dlog.bandstren.value
            UseRecentBest = dlog.rcntbst.value
            LigName = purgeblanks(dlog.LigName.value)
            RNDSD,gotinput = readtextbox(dlog.randseed.value,10000000,4294967296,1,gotinput)
            Verbosity = dlog.vrb.value
            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
            if gotinput ~= 1 then
               print('Something was input wrong. Please try again.')
            end -- if gotinput
        end -- if dialogresult
        --
    end -- while gotinput
end -- if dialogresult
--
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 3')
        -- try to keep <= 25 dialog items between Menu line & OK line.
        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')
        -- try to keep <= 25 dialog items between Menu line & OK line.
        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
            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.5d