Code
-- global variables
-- don't let any depend on functions defined below
recipename='RaiseRigidProteinScore0'
rmin=0.001 -- minimum length for a band
pi=math.pi
dpi=2*pi
hpi=pi/2
subscores=puzzle.GetPuzzleSubscoreNames()
tot=structure.GetCount() -- # of AA's in protein
anum=0 -- atom # on each residue to band to/from
-- 0=default
-- 0,2=alpha-carbon 1=N 2=Ca 3=C 4=O for most segments
-- even for prolines
-- but the segments 1 and tot they could be different
-- RaiseRigidProteinScore0 begun 319pm 12/11/15
-- last updated 12/29/15 1128pm
function Main()
local bandsi,bandsf,initial_score
local note_segment,username,htot
local res,rand_seed,ask,runit
local score,score1,score2,score3
local nwigs,loopno,nloops,smin,smax
local gotbetter,s1,s2,s3,b1,b2,b2tmp
local str1,str2,strlo,strhi,qflag
local r1,r1f,r2,r2e,r2i,r2f,r12,r12i,r12f
local th1,th2,th2e,ph1,ph2,ph2e
local slot1,slot2,bway,wway
local bbeg,bmid,bfin,rbeg,rmid,rfin
local b1o,b2o,b3o,r1o,r2o,r3o
bandsi=band.GetCount() -- # of user-supplied bands
initial_score=current.GetScore()
slot1=21 -- quicksave slot for structure 1 below
slot2=22 -- quicksave slot for structure 2 below
note_segment=noteSegment() -- returns 1-tot if a Note is blank, -1 otherwise
username=user.GetPlayerName()
htot=round0(tot/2)
print(recipename..' begins at '..os.date()..' on puzzle ID '..puzzle.GetPuzzleID())
print(' with name '..puzzle.GetName()..',')
print(' score '..trunc3(initial_score)..', '..tot..' residues, and '..bandsi..' bands.')
print(' ')
-- default values
rand_seed=os.time() -- seed for random # generator
nloops=50 -- # of loops
nwigs=30 -- # times to wiggle each loop
wway=1 -- 1=global wiggle, 2=local wiggle, 3=does both
bway=1 -- makes r1,r2,r3 be fixed (0) or random (1)
-- r12 is starting distance between end pts s1 and s2
smin=0.1 -- smin*r12 sets smallest random band length
smax=1 -- smax*r12 sets largest random band length
-- strlo,strhi should both be from 0-10
strlo=0.2 -- min random band strength
strhi=1 -- max random band strength
-- read inputs
ask=dialog.CreateDialog(recipename)
ask.Label0=dialog.AddLabel(' \n#loops: use 0 to run indefinitely\n ')
ask.nloops=dialog.AddSlider("#loops: ",nloops,0,100,0)
ask.nwigs=dialog.AddSlider("#wiggles/loop: ",nwigs,1,200,0)
ask.Label1=dialog.AddLabel(' \nwway=1 wiggles globally, wway=2 wiggles locally,\nwway=3 does both\n ')
ask.wway=dialog.AddSlider('wway: ',wway,1,3,0)
ask.Label2=dialog.AddLabel(' \nsegments s1,s2,s3 are fixed for bway=0\nand random for bway=1\n ')
ask.bway=dialog.AddSlider('bway: ',bway,0,1,0)
ask.Label3=dialog.AddLabel(' \nsmin & smax: min & max scale factor:\n(band length)/(s1 to s2 distance)\n ')
ask.smin=dialog.AddSlider("smin: ",smin,0.1,10,1)
ask.smax=dialog.AddSlider("smax: ",smax,0.1,10,1)
ask.Label4=dialog.AddLabel(' \nstrlo & strhi: min & max band strength\n ')
ask.strlo=dialog.AddSlider("strlo: ",strlo,0,10,1)
ask.strhi=dialog.AddSlider("strhi: ",strhi,0,10,1)
ask.rand_seed = dialog.AddTextbox("rand_seed: ",rand_seed)
ask.OK = dialog.AddButton("OK", 1)
ask.Cancel = dialog.AddButton("Cancel", 0)
runit = dialog.Show(ask)
if runit==1 then
nloops=ask.nloops.value+0
nwigs=ask.nwigs.value+0
wway=ask.wway.value+0
bway=ask.bway.value+0
smin=ask.smin.value+0
smax=ask.smax.value+0
strlo=ask.strlo.value+0
strhi=ask.strhi.value+0
for res in string.gmatch(ask.rand_seed.value,'[%s%a]*([%d%-%+]+)') do
rand_seed=res
end -- for res
rand_seed=rand_seed+0 -- force to be a number
math.randomseed(rand_seed)
-- want smin<=smax
if smin>smax then
res=smin
smin=smax
smax=res
end -- if smin
-- want strlo<=strhi
if strlo>strhi then
res=strlo
strlo=strhi
strhi=res
end -- if strlo
-- nloops==0 means to run forever
if nloops==0 then
res='(runs forever) '
else
res=''
end -- if nloops
-- print the inputs
print('Got the following inputs by '..os.date()..':')
print(' '..nloops..' loops '..res..'with '..nwigs..' wiggles each (wway='..wway..').')
print(' bway='..bway..' smin='..smin..' smax='..smax..' strlo='..strlo..' strhi='..strhi..' rand_seed='..rand_seed..'.')
print(' ')
print('Clashing importance='..behavior.GetClashImportance()..'.\n ')
bbeg=makezlb(1,tot,htot)
bfin=makezlb(tot,1,htot)
bmid=makezlb(htot,1,tot)
print(' ')
gotbetter=0
qflag=0 -- want qflag==1 to quit main loop
loopno=0
while qflag==0 do -- start main loop
loopno=loopno+1
if nloops>0 then
if loopno>=nloops then
qflag=1 -- will quit after present loop
end -- if loopno
end -- if nloops
if nloops>0 then
print('Loop '..loopno..' of '..nloops..' ('..os.date()..'):')
else -- nloops == 0 means to go forever
print('Loop '..loopno..' ('..os.date()..'):')
end -- if nloops
if gotbetter==1 then
print('Reusing same band parameters.')
else -- pick new random parameters for bands b1 & b2
print('Picking new random band parameters.')
if bway==1 then
-- s1,s2,s3 below are random residue or segment #'s
-- s1,s2 are for 2 diff alpha-carbons for the end points of b1 & b2
-- direction from s1 to s2 is the +x axis
-- s3 is for a 3rd alpha-carbon used to set the xyz axis system for both bands
s1=randi(1,tot)
repeat
s2=randi(1,tot)
until math.abs(s1-s2)>0
-- want s1,s2 > 0 residues apart
repeat
s3=randi(1,tot)
until s3~=s1 and s3~=s2
-- s1,s2,s3 should all be different now
else -- bway==0
s1=1
s3=htot
s2=tot
-- s1,s2,s3 should all be different now
end -- if bway
print(string.format('s1=%d s2=%d s3=%d',s1,s2,s3))
getscores(s1,s2,s3)
-- below are random strengths for bands b1 & b2
str1=randf(strlo,strhi)
str2=randf(strlo,strhi)
-- below are random spherical coordinates for bands b1 & b2, all in radians
th1=randf(0,pi) -- want 0 to pi
th2=randf(0,pi)
ph1=randf(0,dpi) -- want 0 to 2*pi
ph2=randf(0,dpi)
freeze.FreezeAll() -- should freeze both bb & sc
r12=structure.GetDistance(s1,s2) -- get starting distance between end pts
-- below are lengths of bands b1 & b2 in Angstroms
r1=r12*randf(smin,smax) -- can be length 0.0 to 10000.0
r2=r12*randf(smin,smax)
-- below makes sure r1,r2 are both >= rmin (tiny & >0)
if r1<rmin then
r1=rmin
end -- if r1
if r2<rmin then
r2=rmin
end -- if r2
-- get spherical coordinates for band b2tmp
-- also gets forces and torques about center of r12 segment
r2e,th2e,ph2e=getcoords(r12,str1,r1,th1,ph1,str2,r2,th2,ph2)
end -- if gotbetter
--
-- http://foldit.wikia.com/wiki/Foldit_Lua_Functions gives:
--
-- Function: integer band.Add(integer segmentOrigin, integer segmentXAxis, integer segmentYAxis,
-- number rho, number theta, number phi,
-- [integer atomIndexOrigin], [integer atomIndexXAxis], [integer atomIndexYAxis])
-- Description: Add a band to empty space. Returns band number.
--
-- want rho>0 and <10000, th=theta=0 to pi, ph=phi=0 to 2*pi
--
-- for bands b1 and b2tmp, s1 is the origin, and the +x axis goes from s1 to s2.
-- s3 is used to orient the y and z axes.
-- b1 and b2tmp share the same xyz coordinate system.
b1=band.Add(s1,s2,s3,r1,th1,ph1,anum,anum,anum) -- place band b1 from end pt s1
b2tmp=band.Add(s1,s2,s3,r2e,th2e,ph2e,anum,anum,anum) -- place band b2tmp from end pt s1
-- want same xyz coordinate system for bands b1 & b2
-- so next line will use band b2tmp to place band b2
b2=band.AddToBandEndpoint(s2,b2tmp,anum)
-- above places band b2 from end pt s2
-- to where b2tmp ends
--
-- http://foldit.wikia.com/wiki/Foldit_Lua_Functions gives:
--
-- Function: integer band.AddToBandEndpoint(integer segmentIndex, integer bandIndex,
-- [integer atomIndex])
-- Description: Add a band to the endpoint of an existing band. Returns band number.
--
band.Delete(b2tmp) -- remove band b2tmp, should shift all higher band #'s (like b2) down by 1
b2=b2tmp -- band # for band b2 has changed to b2tmp, so account for this here
r2i=band.GetLength(b2) -- get starting length of band b2, should match r2
if r2~=r2i then
print('ERROR: Desired r2='..round2(r2)..' and actual r2='..round2(r2i)..' differ by '..(r2i-r2)..'.\n ')
end -- if r2
band.SetGoalLength(b1,0) -- can be length 0.0 to 10000.0
band.SetGoalLength(b2,0)
band.SetStrength(b1,str1) -- can be strength 0.0 to 10.0
band.SetStrength(b2,str2)
b1o=makezlb(s1,s2,s3)
b2o=makezlb(s2,s1,s3)
b3o=makezlb(s3,s1,s2)
print(' ')
score1=current.GetScore()
print('Score is '..trunc3(score1)..' at '..os.date()..' in loop '..loopno..'.')
save.Quicksave(slot1) -- save structure1
recentbest.Save() -- set recentbest
if wway==1 or wway==3 then
print('Wiggling Globally all backbone up to '..nwigs..' steps.')
structure.WiggleAll(nwigs,true,false)
-- nwigs steps wiggling backbone not sidechains
end -- if wway
if wway==3 then
print('Score is '..trunc3(current.GetScore())..' at '..os.date()..'.')
end -- if wway
if wway==2 or wway==3 then
-- local wiggle below was very slow vs global wiggle 12/15/15
print('Wiggling Locally all backbone up to '..nwigs..' steps.')
structure.LocalWiggleAll(nwigs,true,false)
-- nwigs steps wiggling backbone not sidechains
end -- if wway
score2=current.GetScore()
print('Score is '..trunc3(score2)..' at '..os.date()..'.')
save.Quicksave(slot2) -- save structure2
recentbest.Restore() -- get recentbest
score3=current.GetScore() -- score for recentbest
print('Score changed from '..trunc3(score1)..' to '..trunc3(score2)..' w/recent best '..trunc3(score3)..'.')
if score3>score1 then
if score3>=score2 then -- score3 (recentbest) wins
gotbetter=1 -- score improved
-- recentbest is already in memory, so no need to re-load it
else -- score2 wins
gotbetter=1
save.Quickload(slot2) -- load structure2
end -- if score3
else -- score1>=score3
if score1>=score2 then -- score1 wins
gotbetter=0 -- score stayed same
save.Quickload(slot1) -- load structure1
else -- score2 wins
gotbetter=1
save.Quickload(slot2) -- load structure2
end -- if score1
end -- if score3
score=current.GetScore()
if gotbetter==0 then
print('Keeping original structure w/score '..trunc3(score)..'.')
else -- gotbetter==1
print('Continuing w/new structure w/score '..trunc3(score)..'.')
getscores(s1,s2,s3)
r1f=band.GetLength(b1) -- final length of band b1
r2f=band.GetLength(b2) -- final length of band b2
r12f=structure.GetDistance(s1,s2) -- get final distance between end pts s1 & s2
print('In loop '..loopno..':')
print(string.format(' r1 went from %6.2f to %6.2f (changed by ',round2(r1),round2(r1f))..(r1f-r1)..').')
print(string.format(' r2 went from %6.2f to %6.2f (changed by ',round2(r2),round2(r2f))..(r2f-r2)..').')
print(string.format(' r12 went from %6.2f to %6.2f (changed by ',round2(r12),round2(r12f))..(r12f-r12)..').')
r1o=band.GetLength(b1o)
r2o=band.GetLength(b2o)
r3o=band.GetLength(b3o)
print(' C-alpha s1='..s1..' moved '..r1o..'.')
print(' C-alpha s2='..s2..' moved '..r2o..'.')
print(' C-alpha s3='..s3..' moved '..r3o..'.')
if note_segment>0 then
structure.SetNote(note_segment,string.format("(%s) %.3f + %s (%i) %.3f",
username,trunc3(initial_score),recipename,loopno,trunc3(score)))
end -- if note_segment
end -- if gotbetter
rbeg=band.GetLength(bbeg)
rmid=band.GetLength(bmid)
rfin=band.GetLength(bfin)
print('In all '..loopno..' loops so far:')
print(' C-alpha 1 has moved '..rbeg..'.')
print(' C-alpha '..htot..' has moved '..rmid..'.')
print(' C-alpha '..tot..' has moved '..rfin..'.')
print(string.format('b3o=%d b2o=%d b1o=%d b2=%d b1=%d bmid=%d bfin=%d bbeg=%d',
b3o,b2o,b1o,b2,b1,bmid,bfin,bbeg))
band.Delete(b3o) -- remove band b3o, b2o, b1o, b2, and then b1
band.Delete(b2o)
band.Delete(b1o)
band.Delete(b2)
band.Delete(b1)
print('Just deleted bands b3o b2o b1o b2 and b1.\n ')
end -- while qflag==0
band.Delete(bmid) -- remove band bmid,bfin, and then bbeg
band.Delete(bfin)
band.Delete(bbeg)
print('Just deleted bands bmid bfin and bbeg.\n ')
freeze.UnfreezeAll() -- should unfreeze both bb & sc
end -- if runit
bandsf=band.GetCount() -- final # of bands
score=current.GetScore()
print('Clashing importance='..behavior.GetClashImportance()..'.')
print('Ending with '..bandsf..' bands and score '..trunc3(score)..' at '..os.date()..'.')
end -- Main()
-- make a disabled 0-length and 0-strength band
function makezlb(s1,s2,s3)
local bnum
bnum=band.Add(s1,s2,s3,rmin,hpi,0,anum,anum,anum)
-- make a tiny length band along the +x axis (th=hpi=pi/2=90 degrees, ph=0)
-- where +x is the direction from residue s1 to s2
print(string.format('Making band %d (with %d,%d,%d) be disabled with zero length and strength.',bnum,s1,s2,s3))
band.SetStrength(bnum,0)
band.SetGoalLength(bnum,0)
band.Disable(bnum)
return bnum
end -- makezlb()
-- below finds the effective spherical coordinates
-- r2e,th2e,ph2e for the end pt s2 of band b2
--
-- also finds forces and torques about center of r12 segment between s1 & s2
function getcoords(r12,str1,r1,th1,ph1,str2,r2,th2,ph2)
-- here r1 & r2 are lengths of bands b1 & b2 from residues s1 & s2
local x1,y1,z1,x2,y2,z2 -- coordinates of band end pts s1 & s2
local r2e,th2e,ph2e -- effective spherical coords for band b2
local k1,k2 --- spring constants for bands b1 & b2
local rx -- lever arm for torques
local fx1,fy1,fz1,fx2,fy2,fz2,fx,fy,fz -- forces
local tx1,ty1,tz1,tx2,ty2,tz2,tx,ty,tz -- torques
-- below gets x1,y1,z1 coordinates of band b1's end pt s1
x1=r1*sin(th1)*cos(ph1)
y1=r1*sin(th1)*sin(ph1)
z1=r1*cos(th1)
print(string.format(' r1=%9.2f th1=%9.2f ph1=%9.2f',round2(r1),round2(math.deg(th1)),round2(math.deg(ph1))))
print(string.format(' x1=%9.2f y1=%9.2f z1=%9.2f',round2(x1),round2(y1),round2(z1)))
-- below gets x2,y2,z2 coordinates of band b2's end pt s2
x2=r2*sin(th2)*cos(ph2)+r12
y2=r2*sin(th2)*sin(ph2)
z2=r2*cos(th2)
print(string.format(' r2=%9.2f th2=%9.2f ph2=%9.2f',round2(r2),round2(math.deg(th2)),round2(math.deg(ph2))))
print(string.format(' x2=%9.2f y2=%9.2f z2=%9.2f x2+r12=%9.2f',round2(x2-r12),round2(y2),round2(z2),round2(x2)))
-- shifting x2 by r12 above will
-- change the coordinates r2,th2,ph2 to r2e,th2e,ph2e
-- http://foldit.wikia.com/wiki/Foldit_Lua_Function_band.Add
-- s1=origin, s2=x-axis, s3=defines y and z axes
-- so that the xyz system obeys a right-hand rule
-- next convert x2,y2,z2 to r2e,th2e,ph2e for band b2
r2e=math.sqrt(x2^2+y2^2+z2^2) -- want r2e>0 and tiny
if r2e<rmin then
r2e=rmin -- make r2e>0 and tiny
th2e=0 -- put it along z axis
ph2e=0
else -- r2e is not tiny
th2e=math.acos(z2/r2e) -- should give 0 to pi
if x2==0 and y2==0 then
ph2e=0
else -- x2 or y2 is nonzero
ph2e=math.atan2(y2,x2) -- should give 0 to 2*pi
if ph2e<0 then
ph2e=ph2e+dpi
end -- if ph2e
end -- if x2
end -- if r2e
print(string.format(' r2e=%9.2f th2e=%9.2f ph2e=%9.2f',round2(r2e),round2(math.deg(th2e)),round2(math.deg(ph2e))))
print(' ')
-- shift x2 back by r12
x2=x2-r12
-- below gets forces
k1=str1^2 -- a band's strength is the square
k2=str2^2 -- of its Hooke's Law spring constant
fx1=k1*x1
fy1=k1*y1
fz1=k1*z1
fx2=k2*x2
fy2=k2*y2
fz2=k2*z2
fx=fx1+fx2
fy=fy1+fy2
fz=fz1+fz2
print(string.format('str1=%9.2f str2=%9.2f k1=%9.2f k2=%9.2f',round2(str1),round2(str2),round2(k1),round2(k2)))
print(string.format(' fx1=%9.2f fy1=%9.2f fz1=%9.2f',round2(fx1),round2(fy1),round2(fz1)))
print(string.format(' fx2=%9.2f fy2=%9.2f fz2=%9.2f',round2(fx2),round2(fy2),round2(fz2)))
print(string.format(' fx=%9.2f fy=%9.2f fz=%9.2f',round2(fx),round2(fy),round2(fz)))
print(' ')
-- below gets torques about center of r12 segment,
-- as if s1 is centered at x= -rx= -r12/2
-- and s2 is centered at x= rx= r12/2
-- forces in x direction give no torques
rx=r12/2
tx1= 0
ty1= rx*fz1 -- CCW>0 since s2's x > s1's x (s2's x = s1's x + r12)
tz1= -rx*fy1 -- CW<0
tx2= 0
ty2= -rx*fz2 -- CW<0
tz2= rx*fy2 -- CCW>0
tx=tx1+tx2
ty=ty1+ty2
tz=tz1+tz2
print(string.format(' rx=%9.2f r12=%9.2f',round2(rx),round2(r12)))
print(string.format(' tx1=%9.2f ty1=%9.2f tz1=%9.2f',round2(tx1),round2(ty1),round2(tz1)))
print(string.format(' tx2=%9.2f ty2=%9.2f tz2=%9.2f',round2(tx2),round2(ty2),round2(tz2)))
print(string.format(' tx=%9.2f ty=%9.2f tz=%9.2f',round2(tx),round2(ty),round2(tz)))
print(' ')
return r2e,th2e,ph2e
end -- getcoords()
-- list subscores for whole protein
-- and then for each of residues s1,s2,s3
function getscores(s1,s2,s3)
local fn,n,nn
local ngot,aa,ss,sstr
fn={}
fn.score1=current.GetScore()
fn.segscore1={}
fn.subscore={}
fn.score2=0 -- total of fn.segscore1 for all segments
fn.segscore2={} -- total of fn.subscore for all subscores
fn.score3=0 -- total of fn.subscore for all segments & subscores
-- same as total of fn.segscore2 for all segments
fn.subscoretot={} -- total of fn.subscore for all segments
for nn=1,#subscores do
fn.subscoretot[nn]=0
end -- for nn
for n=1,tot do
fn.segscore1[n]=current.GetSegmentEnergyScore(n)
fn.segscore2[n]=0
fn.score2=fn.score2+fn.segscore1[n]
fn.subscore[n]={}
for nn=1,#subscores do
fn.subscore[n][nn]=current.GetSegmentEnergySubscore(n,subscores[nn])
fn.subscoretot[nn]=fn.subscoretot[nn]+fn.subscore[n][nn]
fn.segscore2[n]=fn.segscore2[n]+fn.subscore[n][nn]
end -- for nn
fn.score3=fn.score3+fn.segscore2[n]
end -- for n
print(' ')
print('Total score 1,2 = '..fn.score1..', '..fn.score2..' (differ by '..(fn.score2-fn.score1)..').')
print('Total score 1,3 = '..fn.score1..', '..fn.score3..' (differ by '..(fn.score3-fn.score1)..').')
print('Total score 2,3 = '..fn.score2..', '..fn.score3..' (differ by '..(fn.score3-fn.score2)..').')
for nn=1,#subscores do
if fn.subscoretot[nn]~=0 then
print('Total for '..subscores[nn]..' = '..fn.subscoretot[nn]..'.')
end -- if fn
end -- for nn
print(' ')
ngot=0
for n=1,tot do
if n==s1 or n==s2 or n==s3 then
if n==s1 then
sstr='s1'
elseif n==s2 then
sstr='s2'
elseif n==s3 then
sstr='s3'
else
sstr='s?'
end
ngot=ngot+1
aa=structure.GetAminoAcid(n)
ss=structure.GetSecondaryStructure(n)
print('('..sstr..') Segment '..aa..n..' has ss '..ss..' and the subscores below:')
for nn=1,#subscores do
if fn.subscore[n][nn]~=0 then
print(' '..subscores[nn]..' = '..fn.subscore[n][nn]..'.')
end -- fn
end -- for nn
print('('..sstr..') Segment '..aa..n..' has score '..fn.segscore1[n]..' = '..fn.segscore2[n]..' (differ by '..(fn.segscore2[n]-fn.segscore1[n])..').')
if ngot<=3 then
print(' ')
end -- if ngot
end -- if n
end -- for n
end -- getscores()
function cos(val)
return math.cos(val)
end -- cos()
function sin(val)
return math.sin(val)
end -- sin()
-- below returns a random integer from xlo to xhi
function randi(xlo,xhi)
return math.random(xlo,xhi)
end -- randi()
-- below returns a random real number from xlo to xhi
function randf(xlo,xhi)
local res
res=xlo+(xhi-xlo)*math.random() -- gives [xlo,xhi) since math.random() gives [0,1)
return res
end -- randf()
-- below rounds down to 3 decimal places
function trunc3(numi)
local numo=trunc0(numi*1000)/1000
return numo
end -- trunc3()
-- below rounds val down to nearest integer
function trunc0(val)
return math.floor(val)
end -- trunc0()
-- below rounds to 3 decimal places
function round3(numi)
local numo=round0(numi*1000)/1000
return numo
end -- round3()
-- below rounds to 2 decimal places
function round2(numi)
local numo=round0(numi*100)/100
return numo
end -- round2()
-- below rounds to 1 decimal places
function round1(numi)
local numo=round0(numi*10)/10
return numo
end -- round1()
-- below rounds val to the nearest integer to get res
function round0(val)
local res=math.floor(val)
if val-res>=0.5 then
res=res+1
end -- if val
return res
end -- round0()
-- finds 1st residue with a blank note
-- assuming all notes after that one are blank
function noteSegment()
local i
for i=tot,1,-1 do
if structure.GetNote(i) ~= "" then
if i == tot then return(-1) else return(i+1) end
end -- if
end -- for
return(1)
end -- function
Main()
