Profile
- Name
- RaiseRigidProteinScore1
- ID
- 101743
- Shared with
- Public
- Parent
- RaiseRigidProteinScore0
- Children
- None
- Created on
- January 05, 2016 at 23:34 PM UTC
- Updated on
- January 05, 2016 at 23:34 PM UTC
- Description
*1b.txt 1/5/16 531pm code
Best for
Code
-- global variables
-- don't let any depend on functions defined below
recipename='RaiseRigidProteinScore1'
rmin=0.001 -- minimum length for a band
pi=math.pi
dpi=2*pi
hpi=pi/2
verbose=0
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 for segments 1 and tot, they could be different
--
-- RaiseRigidProteinScore0 (RRPS0) begun 319pm 12/11/15 by jeff101
--
-- RaiseRigidProteinScore1 (RRPS1) begun 949am 1/5/16 by jeff101
-- adapted from RRPS0i.txt 12/29/15 1128pm code
--
-- last updated 1/5/16 531pm
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,wflag
local cibandson,cibandsoff,ciold
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)
ciold=behavior.GetClashImportance()
print(recipename..' begins at '..os.date()..' on puzzle ID '..puzzle.GetPuzzleID())
print(' with name '..puzzle.GetName()..',')
print(' score '..trunc3(initial_score)..', '..tot..' residues, '..bandsi..' bands, and ci='..ciold..'.')
print(' ')
-- default values
nloops=50 -- # of loops
nwigs=30 -- # times to wiggle each loop
wway=1 -- 1=global wiggle, 2=local wiggle, 3=global then local, 4=local then global
wflag=true -- true wiggles w/bands then w/o bands, false just wiggles w/bands
cibandson=0.1 -- clashing importance for wiggling with bands enabled
cibandsoff=1 -- clashing importance for wiggling with bands disabled
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
verbose=0 -- 0=most quiet, 1=less quiet, 2=noisiest
rand_seed=os.time() -- seed for random # generator
-- 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, 2 wiggles locally,\n3 globally then locally, 4 locally then globally\n ')
ask.wway=dialog.AddSlider('wway: ',wway,1,4,0)
ask.wflag=dialog.AddCheckbox('Also wiggle w/disabled bands',wflag)
ask.Label2=dialog.AddLabel(' \ncibandson is clashing importance for enabled bands,\ncibandsoff is ci for disabled bands\n ')
ask.cibandson=dialog.AddSlider('cibandson: ',cibandson,0,1,2)
ask.cibandsoff=dialog.AddSlider('cibandsoff: ',cibandsoff,0,1,2)
ask.Label3=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.Label4=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.Label5=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.Label6=dialog.AddLabel(' \nverbose: 0=least, 1=more, 2=most output\n ')
ask.verbose=dialog.AddSlider('verbose: ',verbose,0,2,0)
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
wflag=ask.wflag.value
cibandson=ask.cibandson.value+0
cibandsoff=ask.cibandsoff.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
verbose=ask.verbose.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..' wflag='..tostring(wflag)..').')
print(' cibandson='..cibandson..' cibandsoff='..cibandsoff..' bway='..bway..' smin='..smin..' smax='..smax..'.')
print(' strlo='..strlo..' strhi='..strhi..' verbose='..verbose..' rand_seed='..rand_seed..'.')
print(' ')
bbeg=makezlb(1,tot,htot)
bfin=makezlb(tot,1,htot)
bmid=makezlb(htot,1,tot)
if verbose>=0 then
print(' ')
end -- if verbose
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:')
printpars(s1,s2,s3,r1,th1,ph1,str1,r2,th2,ph2,str2)
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
-- 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
printpars(s1,s2,s3,r1,th1,ph1,str1,r2,th2,ph2,str2)
if loopno==1 or verbose>1 then
getscores(s1,s2,s3,3) -- print both total & segment scores
else
getscores(s1,s2,s3,2) -- only print segment scores
end -- if loopno
-- 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 math.abs(r2-r2i)>=0.01 then
res='ERROR: '
else
res=''
end -- if math.abs
print(res..'Desired r2='..round2(r2)..' and actual r2='..round2(r2i)..' differ by '..(r2i-r2)..'.\n ')
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
-- wway: 1=global wiggle, 2=local wiggle, 3=global then local wiggle, 4=local then global wiggle
dowigs(wway,wflag,nwigs,cibandson,cibandsoff,ciold,b1,b2)
if wway==3 then
dowigs(2,wflag,nwigs,cibandson,cibandsoff,ciold,b1,b2) -- do local wiggle w/no message
elseif wway==4 then
dowigs(1,wflag,nwigs,cibandson,cibandsoff,ciold,b1,b2) -- do global wiggle w/no message
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)..'.')
-- do below 2 lines in case recentbest comes from wiggling w/bands disabled
band.Enable(b1)
band.Enable(b2)
-- do above 2 lines in case recentbest comes from wiggling w/bands disabled
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 > score3 > score1, so score2 wins
gotbetter=1
save.Quickload(slot2) -- load structure2
end -- if score3
else -- score1>=score3
if score1>=score2 then -- score1 >= score2,score3, so score1 wins
gotbetter=0 -- score stayed same
save.Quickload(slot1) -- load structure1
else -- score2 > score1 >= score3, so 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)..'.\n ')
else -- gotbetter==1
print('Continuing w/new structure w/score '..trunc3(score)..'.\n ')
getscores(s1,s2,s3,3) -- print both total & segment scores
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..'.\n ')
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.')
if qflag==0 then
print(' ')
end -- if qflag
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()
behavior.SetClashImportance(ciold)
print('Ending with '..bandsf..' bands, ci='..ciold..', and score '..trunc3(score)..' at '..os.date()..'.')
end -- Main()
function printpars(s1,s2,s3,r1,th1,ph1,str1,r2,th2,ph2,str2)
print(string.format(' s1=%d s2=%d s3=%d',s1,s2,s3))
print(string.format(' r1=%6.2f th1=%6.2f ph1=%6.2f str1=%5.2f',round2(r1),round2(math.deg(th1)),round2(math.deg(ph1)),round2(str1)))
print(string.format(' r2=%6.2f th2=%6.2f ph2=%6.2f str2=%5.2f',round2(r2),round2(math.deg(th2)),round2(math.deg(ph2)),round2(str2)))
print(' ')
end -- printpars()
-- 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
if verbose>=0 then
print(string.format('Making band %d (with %d,%d,%d) be disabled with zero length and strength.',bnum,s1,s2,s3))
end -- if verbose
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)
if verbose>0 then
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)))
end -- if verbose
-- 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)
if verbose>0 then
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)))
end -- if verbose
-- 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
if verbose>0 then
print(string.format(' r2e=%9.2f th2e=%9.2f ph2e=%9.2f',round2(r2e),round2(math.deg(th2e)),round2(math.deg(ph2e))))
print(' ')
end -- if verbose
-- 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
if verbose>1 then
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(' ')
end -- if verbose
-- 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
if verbose>1 then
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(' ')
end -- if verbose
return r2e,th2e,ph2e
end -- getcoords()
-- list subscores for whole protein
-- and then for each of residues s1,s2,s3.
-- sflag=1 just shows totals
-- sflag=2 just shows segment scores
-- sflag=3 shows both
function getscores(s1,s2,s3,sflag)
local fn,n,nn
local ngot,aa,ss,sstr
if verbose>=0 then
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
end -- if verbose
if verbose>=0 and (sflag==1 or sflag==3) then
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(' ')
end -- if verbose
if verbose>1 and (sflag==2 or sflag==3) then
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 -- if n
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 -- if 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 -- if verbose
end -- getscores()
-- in below, wway=1 for global wiggle, 2 for local wiggle, 3 for global wiggle then mssg, 4 for local wiggle then mssg
function dowigs(wway,wflag,nwigs,cibandson,cibandsoff,ciold,b1,b2)
behavior.SetClashImportance(cibandson)
if wway==1 or wway==3 then
print('Globally wiggling backbone '..nwigs..' steps w/bands b1='..b1..' b2='..b2..' enabled and ci='..cibandson..'.')
structure.WiggleAll(nwigs,true,false) -- nwigs steps wiggling backbone not sidechains
else -- if wway==2 or wway==4 then
print('Locally wiggling backbone '..nwigs..' steps w/bands b1='..b1..' b2='..b2..' enabled and ci='..cibandson..'.')
structure.LocalWiggleAll(nwigs,true,false) -- nwigs steps wiggling backbone not sidechains
end -- if wway
behavior.SetClashImportance(ciold)
if wflag then
print('Score is '..trunc3(current.GetScore())..' at '..os.date()..'.')
band.Disable(b1)
band.Disable(b2)
behavior.SetClashImportance(cibandsoff)
if wway==1 or wway==3 then
print('Globally wiggling backbone '..nwigs..' steps w/bands b1='..b1..' b2='..b2..' disabled and ci='..cibandsoff..'.')
structure.WiggleAll(nwigs,true,false) -- nwigs steps wiggling backbone not sidechains
else -- if wway==2 or wway==4 then
print('Locally wiggling backbone '..nwigs..' steps w/bands b1='..b1..' b2='..b2..' disabled and ci='..cibandsoff..'.')
structure.LocalWiggleAll(nwigs,true,false) -- nwigs steps wiggling backbone not sidechains
end -- if wway
behavior.SetClashImportance(ciold)
band.Enable(b1)
band.Enable(b2)
end -- if wflag
if wway==3 or wway==4 then
print('Score is '..trunc3(current.GetScore())..' at '..os.date()..'.')
end -- if wway
end -- dowigs()
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 -- if i
end -- if structure
end -- for i
return(1)
end -- noteSegment()
Main()