Code
-- shiftit0 begun 8/11/15 by jeff101
--
-- goal is to move each alpha-carbon to another
-- alpha-carbon's position using bands
--
-- last updated 8/13/15 238pm
function Main()
local rnam='shiftit0'
local qflag,ask,askres,verbose
local input_dir,dir,nstep,stepn,aaperstep,aaperband
local obands,keep_bands,rand_orient,rand_seed
local abandlist={1,2,3,4,6,8}
local abands,abandn,abandways,abandstr
local rtot,rnum,rnumo,rnumy,lastrnum
local bnum1,bnum2,blen,bstr,r,th,ph
local bxyz,tbxyz,rbxyz,trbxyz
local bdots,nbdots,angmat
local nwigs,wigway,newss
local sqrt3,sqrt2,th4,th8,th4d,th8d,pi
abandways=#abandlist
abandstr=''
for abandn=1,abandways do
abandstr=(abandstr..' '..abandlist[abandn])
end -- for abandn
sqrt3=math.sqrt(3)
sqrt2=math.sqrt(2)
pi=math.pi
th4=math.acos(-1/3) -- in radians
th8=math.acos(1/3) -- in radians
th4d=math.deg(th4) -- in degrees
th8d=math.deg(th8) -- in degrees
obands=band.GetCount() -- starting # of bands
rtot=structure.GetCount() -- # of residues
print('Running '..rnam..' at '..os.date()..' on puzzle ID '..puzzle.GetPuzzleID())
print(' '..puzzle.GetName())
print(' w/score='..trunc3(current.GetScore())..', '..rtot..' residues, and '..obands..' user-supplied bands.\n ')
print(string.format('sqrt(2)=%.3f sqrt(3)=%.3f pi=%.3f.\nth4d=%.2f and th8d=%.2f degrees.\n ',sqrt2,sqrt3,pi,th4d,th8d))
-- below are default values
nstep=1 -- # of moves up or down
aaperstep=10 -- # amino acids to move per step
aaperband=1 -- # amino acids between bands
input_dir=1 -- 0 moves down, 1 moves up
newss='H' -- 2nd str to put on new residues. want 'H' 'E' or 'L'.
nwigs=100 -- # of wiggles per move
wigway=3 -- # wiggling strategy
abands=1 -- # of bands per atom
rand_orient=0 -- 0 means not random, 1 means random
blen=10 -- length of bands to space, must be >0.0 and <10000.0
bstr=10 -- strength of bands, can be 0.0 to 10.0
rtiny=0.001 -- length to use instead of zero for bands to space (is 0.001 ok?)
keep_bands=1 -- 0 means remove all but obands at end, 1 means keep all bands at end
verbose=0 -- 0=quiet, 1=noisy, 2=more noisy
rand_seed = os.time() -- the default seed for random # generator
qflag=0
while qflag==0 do
-- get input below
ask=dialog.CreateDialog((rnam..' inputs: page 1 of 2'))
ask.Input1=dialog.AddSlider("nstep: ",nstep,1,rtot,0)
ask.Label1=dialog.AddLabel(" \nnstep: # of shifts to do.\n ")
ask.Input1a=dialog.AddSlider("aaperstep: ",aaperstep,1,rtot,0)
ask.Label1a=dialog.AddLabel(" \naaperstep: # residues moved by per shift.\n ")
ask.Input1b=dialog.AddSlider("aaperband: ",aaperband,1,rtot,0)
ask.Label1b=dialog.AddLabel(" \naaperband: # residues between sets of bands.\n ")
ask.Input2=dialog.AddSlider("input_dir: ",input_dir,0,1,0)
ask.Label2=dialog.AddLabel(" \ninput_dir: 0 shifts down, 1 shifts up.\n ")
ask.Input2a=dialog.AddTextbox("newss: ",newss)
ask.Label2a=dialog.AddLabel(" \nnewss: secondary structure for new residues: H E L.\n ")
ask.Input2b=dialog.AddSlider("nwigs: ",nwigs,1,500,0)
ask.Label2b=dialog.AddLabel(" \nnwigs: # of wiggles per shift.\n ")
ask.Input2c=dialog.AddSlider("wigway: ",wigway,1,4,0)
ask.Label2c=dialog.AddLabel(" \nwigway: 1=global only, 2=local only,\n 3=global then local, 4=local then global.\n ")
ask.Input3=dialog.AddSlider("abands: ",abands,1,abandways,0)
ask.Label3=dialog.AddLabel((" \nabands: # of bands per atom.\n 1-"..abandways.." here stand for"..abandstr..", respectively.\n "))
ask.OK = dialog.AddButton("OK", 1)
askres=dialog.Show(ask) -- should give a 1 here since 1 button
-- the +0's below make sure inputs become numbers
nstep=ask.Input1.value+0
aaperstep=ask.Input1a.value+0
aaperband=ask.Input1b.value+0
input_dir=ask.Input2.value+0
newss=ask.Input2a.value -- no +0 since want it to stay a string
nwigs=ask.Input2b.value+0
wigway=ask.Input2c.value+0
abands=ask.Input3.value+0
ask=dialog.CreateDialog((rnam..' inputs: page 2 of 2'))
ask.Input4=dialog.AddSlider("rand_orient: ",rand_orient,0,1,0)
ask.Label4=dialog.AddLabel(" \nrand_orient: 0=non-random, 1=random.\n ")
ask.Input5=dialog.AddTextbox("blen: ",blen)
ask.Label5=dialog.AddLabel(" \nblen: length of bands to space.\n Want 0.0 < blen < 10000.0.\n ")
ask.Input6=dialog.AddSlider("bstr: ",bstr,0,10,1) -- 1 is for increments of 0.1
ask.Label6=dialog.AddLabel(" \nbstr: strength of bands to space.\n ")
ask.Input7=dialog.AddTextbox("rtiny: ",rtiny)
ask.Label7=dialog.AddLabel(" \nrtiny: true length for zero-length bands.\n ")
ask.Input8=dialog.AddSlider("keep_bands: ",keep_bands,0,1,0)
ask.Label8=dialog.AddLabel(" \nkeep_bands: 0 only removes new bands when done,\n 1 keeps all bands when done.\n ")
ask.Input9=dialog.AddSlider("verbose: ",verbose,0,2,0)
ask.Label9=dialog.AddLabel(" \nverbose: 0=quiet, 1=noisy, 2=noisiest.\n ")
ask.Input9b=dialog.AddTextbox("rand_seed: ",rand_seed)
ask.Label9b=dialog.AddLabel(" \nrand_seed: sets sequence of random #'s.\n ")
ask.OK = dialog.AddButton("OK", 1)
askres=dialog.Show(ask) -- should give a 1 here since 1 button
rand_orient=ask.Input4.value+0
blen=ask.Input5.value+0
bstr=ask.Input6.value+0
rtiny=ask.Input7.value+0
keep_bands=ask.Input8.value+0
verbose=ask.Input9.value+0
rand_seed=ask.Input9b.value+0
-- input_dir: 0 moves down, 1 moves up
dir=2*input_dir-1
-- dir: -1 moves down, 1 moves up
qflag=1
if (newss~='H' and newss~='E' and newss~='L') then
print('newss='..newss..' should be H E or L.')
newss='H'
qflag=0
end -- if newss
if rtiny<=0 then
print('rtiny='..rtiny..'<=0 is not allowed.')
rtiny=0.001
qflag=0
end -- if rtiny
if blen<=0 or blen>=10000 then
print('blen='..blen..' but want 0<blen<10000.')
blen=10
qflag=0
end -- if blen
if rand_seed<0 then
print('rand_seed='..rand_seed..'<0 is not allowed.')
rand_seed = os.time()
qflag=0
elseif rand_seed~=math.floor(rand_seed) then
print('rand_seed='..rand_seed..' must be an integer.')
rand_seed = os.time()
qflag=0
end -- if rand_seed
end -- while qflag
math.randomseed(rand_seed) -- use rand_seed here
print('Since abands=1-'..abandways..' means'..abandstr..', abands='..abands..' becomes '..abandlist[abands]..'.\n ')
abands=abandlist[abands]
print('Got the inputs: nstep='..nstep..' aaperstep='..aaperstep..' aaperband='..aaperband)
print(' input_dir='..input_dir..' dir='..dir..' newss='..newss..' nwigs='..nwigs..' wigway='..wigway)
print(' abands='..abands..' rand_orient='..rand_orient..' blen='..blen..' bstr='..bstr..' rtiny='..rtiny)
print(' keep_bands='..keep_bands..' verbose='..verbose..' and rand_seed='..rand_seed..'.\n ')
-- In below try to have bands stick out in y and z dirs.
-- x axis will be along backbone of protein (unless rand_orient==1).
-- In bxyz, vector lengths don't matter but directions do matter.
if abands==1 then -- zero length band
bxyz={{0,0,0}} -- {{}} so bxyz=(1 row)x(3 col) matrix
elseif abands==2 then -- linear
bxyz={{0,1,0},{0,-1,0}} -- bxyz=(2 row)x(3 col) matrix
-- above 2 are 180 degrees apart and perp to x-axis
elseif abands==3 then -- trigonal planar
bxyz={{0,2,0},{0,-1,sqrt3},{0,-1,-sqrt3}}
-- each pair above is 120 degrees apart and perp to x-axis
elseif abands==4 then -- tetrahedral
bxyz={{1,1,1},{1,-1,-1},{-1,1,-1},{-1,-1,1}}
-- 4 non-adjacent corners of a cube
-- each pair above is th4=acos(-1/3) apart
elseif abands==6 then -- octahedral
bxyz={{sqrt2,-2,0},{sqrt2,1,sqrt3},{sqrt2,1,-sqrt3},{-sqrt2,2,0},{-sqrt2,-1,-sqrt3},{-sqrt2,-1,sqrt3}}
-- easy way is along +-x,y,z axes
-- but want em to stick out to the sides here
-- each pair above is either 90 or 180 degrees apart
elseif abands==8 then -- all 8 corners of a cube
bxyz={{-1,-1,-1},{-1,-1,1},{-1,1,1},{-1,1,-1},{1,1,-1},{1,1,1},{1,-1,1},{1,-1,-1}}
-- each pairs above is 180 degrees, th8=acos(1/3), or th4=acos(-1/3) apart
else
print('Error: abands='..abands..' is not allowed.')
end
listmat(bxyz,'bxyz=unnormalized band vectors',-3)
tbxyz=xpose(bxyz)
bdots=matxmat(bxyz,tbxyz) -- dot products of bxyz directions
nbdots=normmat(bdots) -- rescale bdots by bdots[1][1]
listmat(nbdots,'nbdots=dot products of normalized band vectors',-3) -- diagonal should be all 1's
angmat=getangs(nbdots) -- gets angles using acos(nbdots)
-- and lists them in degrees
listmat(angmat,'angmat=angles in degrees from dot products',-2)
print(' ')
for stepn=1,nstep do
print('Step '..stepn..' of '..nstep..' starts at '..os.date()..' w/score='..trunc3(current.GetScore())..':')
keepbands(obands)
obands=band.GetCount()
print(' Keeping just '..obands..' bands.')
if dir==1 then
rbeg=1
rend=rtot-aaperstep
else -- dir == -1
rbeg=rtot
rend=1+aaperstep
end -- if dir
lastrnum=rbeg -- lastrnum will keep track of last rnum done in loop below
-- if both aaperstep and aaperband are 1, loop should end with lastrnum=rend
for rnum=rbeg,rend,dir*aaperband do -- start,end,stepsize
-- aaperband tells how many residues to skip per set of bands.
-- aaperband=1 lets rnum go from rbeg to rend.
-- For general aaperband, rnum won't quite reach rend.
--
-- rnumy below is residue # used to set y direction for band-to-space coordinates.
if rnum~=rend then -- rnum<rend for dir=1. rnum>rend for dir= -1.
rnumy=rnum+dir*(1+aaperstep)
else -- rnum==rend, so doing last residue
rnumy=rnum-dir*aaperstep
end -- if rnum
--
-- in below we will add bands to move residue rnum to position of residue rnumo
rnumo=rnum+dir*aaperstep
--
-- want rnum,rnumo,rnumy to always be in range 1 to rtot
-- and different from each other
--
lastrnum=rnum -- record latest value of rnum in range rbeg to rend
-- if both aaperstep and aaperband are 1, loop should end with lastrnum=rend
--
if rand_orient==0 or abands==1 then
-- abands=1 is a zero-length band, so its orientation doesn't matter.
-- Randomly orienting it won't matter and wastes time.
rbxyz=copymat(bxyz)
trbxyz=copymat(tbxyz)
else -- rand_orient==1 (aims the bands at random)
rbxyz,trbxyz=randbxyz(bxyz,tbxyz)
end -- if rand_orient
--
-- just need rbxyz below
--
for abandn=1,abands do
r,th,ph=xyz2rtp(rbxyz[abandn],pi,verbose)
--
-- 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
if r==0 then -- use "zero" length band (actually length rtiny) to space
bnum1=band.Add(rnumo,rnum,rnumy,rtiny,th,ph) -- band target position rnumo to space
else -- use band of length blen to space
bnum1=band.Add(rnumo,rnum,rnumy,blen,th,ph) -- band target position rnumo to space
end -- if r==0
if bnum1<=0 then -- band bnum1 was not made fine
print(' Error making bnum1 for band '..abandn..' of '..abands..' to move residue '..rnum..' to '..rnumo..'.')
if r==0 then
print('rtiny='..rtiny..' may be too small.')
end -- if r==0
else -- band bnum1 was made fine
--
-- 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.
--
bnum2=band.AddToBandEndpoint(rnum,bnum1) -- band initial position rnum to endpoint in space of bnum1
if bnum2<=0 then -- band bnum2 was not made fine
print(' Error making bnum2 for band '..abandn..' of '..abands..' to move residue '..rnum..' to '..rnumo..'.')
else -- band bnum2 was made fine
-- below makes band bnum2 strength bstr with target length 0
band.SetStrength(bnum2,bstr) -- can be strength 0.0 to 10.0
band.SetGoalLength(bnum2,0) -- can be length 0.0 to 10000.0
end -- if bnum2
band.Delete(bnum1) -- It's safe to delete bnum1 now that bnum2 is done.
-- If you delete bnum1 before bnum2 is done, band #'s shift down,
-- and bnum2 will point nowhere, giving an error.
end -- if bnum1
end -- for abandn
structure.SetSecondaryStructure(rnum,structure.GetSecondaryStructure(rnumo))
end -- for rnum
--
-- At end of loop, final value of rnum is rend if both aaperstep and aaperband are 1,
-- but in general final value of rnum is lastrnum.
--
-- Below line worked fine when both aaperstep and aaperband were 1:
-- structure.SetSecondaryStructure(rend+dir,newss)
--
-- Below lines are for the more general case:
for rnum=(lastrnum+dir),(rend+dir*aaperstep),dir do
structure.SetSecondaryStructure(rnum,newss)
end -- for rnum
print(' Now have '..band.GetCount()..' bands on the protein.')
print(' Score='..trunc3(current.GetScore())..' at '..os.date()..'.')
-- Below lines wiggle the protein.
if wigway==1 or wigway==3 then
print(' Next global wiggle the backbone '..nwigs..' times.')
structure.WiggleAll(nwigs,true,false)
print(' Score='..trunc3(current.GetScore())..' at '..os.date()..'.')
end -- if wigway
if wigway==2 or wigway==3 or wigway==4 then
print(' Next local wiggle the backbone '..nwigs..' times.')
structure.LocalWiggleAll(nwigs,true,false)
print(' Score='..trunc3(current.GetScore())..' at '..os.date()..'.')
end -- if wigway
if wigway==4 then
print(' Next global wiggle the backbone '..nwigs..' times.')
structure.WiggleAll(nwigs,true,false)
print(' Score='..trunc3(current.GetScore())..' at '..os.date()..'.')
end -- if wigway
-- Above lines did WiggleAll (global) or LocalWiggleAll (local)
-- on backbone not sidechain for nwigs iterations each.
-- https://fold.it/portal/node/994671 discusses global vs local wiggling.
print(' Step '..stepn..' of '..nstep..' ends at '..os.date()..' w/score='..trunc3(current.GetScore())..'.\n ')
end -- for stepn
if keep_bands==0 then
keepbands(obands)
end
obands=band.GetCount()
print('End of '..rnam..' at '..os.date()..' w/score='..trunc3(current.GetScore())..' and '..obands..' bands.')
end -- Main()
-- below removes all bands with numbers > obands
function keepbands(obands)
local btot,bnum
btot=band.GetCount()
for bnum=btot,obands+1,-1 do -- start,end,stepsize
band.Delete(bnum)
end -- for bnum
end -- keepbands()
function randbxyz(bxyz,tbxyz)
local alpd,betd,gamd,alp,bet,gam
local sa,sb,sg,ca,cb,cg
local ralp,rbet,rgam,rba,rgba
local rbxyz,trbxyz
local bdots,nbdots,angmat
-- http://easyspin.org/documentation/eulerangles.html gives
alpd=360*math.random() -- returns [0,360) since math.random() gives [0,1)
betd=360*math.random() -- returns [0,360) (ideally want 0-180)
-- 0 and 360 give the same results, so [0,360) covering 0 not 360 is fine.
-- But 0 and 180 give different results.
-- [0,360) includes both 0 and 180.
-- [0,180) includes 0 not 180.
gamd=360*math.random() -- returns [0,360)
alp=math.rad(alpd) -- converts to radians so 0 to 2*pi
bet=math.rad(betd)
gam=math.rad(gamd)
print(string.format(" \n Using alp,bet,gam=%6.2f,%6.2f,%6.2f degrees\n or alp,bet,gam=%6.3f,%6.3f,%6.3f radians.\n ",alpd,betd,gamd,alp,bet,gam))
sa=math.sin(alp)
sb=math.sin(bet)
sg=math.sin(gam)
ca=math.cos(alp)
cb=math.cos(bet)
cg=math.cos(gam)
-- http://easyspin.org/documentation/eulerangles.html gives
ralp={{ca, sa, 0},{-sa, ca, 0},{ 0, 0, 1}} -- rot around z
rbet={{cb, 0, -sb},{ 0, 1, 0},{sb, 0, cb}} -- rot around y
rgam={{cg, sg, 0},{-sg, cg, 0},{ 0, 0, 1}} -- rot around z
rba=matxmat(rbet,ralp)
rgba=matxmat(rgam,rba)
-- rgba should now be rotation matrix
-- with a random orientation based on
-- random euler angles alp,bet,gam
listmat(rgba,'rgba=rotation matrix',-3)
trbxyz=matxmat(rgba,tbxyz)
rbxyz=xpose(trbxyz)
listmat(rbxyz,'rbxyz=rotated unnormalized band vectors',-3)
bdots=matxmat(rbxyz,trbxyz) -- dot products of rbxyz directions
nbdots=normmat(bdots) -- rescale bdots by bdots[1][1]
listmat(nbdots,'nbdots=dot products of normalized band vectors',-3) -- diagonal should be all 1's
angmat=getangs(nbdots) -- gets angles using acos(nbdots)
-- and lists them in degrees
listmat(angmat,'angmat=angles in degrees from dot products',-2)
print(' ')
return rbxyz,trbxyz
end -- randbxyz()
function xyz2rtp(xyz,pi,verbose)
-- converts vector xyz containing x,y,z coordinates
-- into r,theta,phi (spherical coordinates)
-- want th=theta=0 to pi, ph=phi=0 to 2*pi
-- want x=y=z=0 to give r=th=ph=0
local x,y,z,r,th,ph
x=xyz[1]
y=xyz[2]
z=xyz[3]
r=math.sqrt(x^2+y^2+z^2)
if r~=0 then
th=math.acos(z/r) -- gives 0 to pi
ph=math.atan2(y,x) -- want 0 to 2*pi
if ph<0 then
ph=ph+2*pi
end -- if ph
else -- r==0
th=0
ph=0
end -- if r~=0
if verbose>0 then
print(string.format('x,y,z=%.3f,%.3f,%.3f gives r=%.3f and th,ph=%.2f,%.2f degrees.',x,y,z,r,math.deg(th),math.deg(ph)))
end -- if verbose
return r,th,ph
end -- xyz2rtp()
function getangs(mat)
local rown,coln,nrow,ncol,omat
nrow=#mat
ncol=#(mat[1])
omat=mat
omat={}
for rown=1,nrow do
omat[rown]={}
for coln=1,ncol do
omat[rown][coln]=math.deg(math.acos(mat[rown][coln])) -- gets angle in degrees
end -- for coln
end -- for rown
return omat
end -- getangs()
function normmat(mat)
local rown,coln,nrow,ncol,cf,omat
cf=mat[1][1] -- Most bxyz have 2 nonzero-length vectors,
-- so each bxyz will give a dot product matrix
-- with all on diagonal elements equal and nonzero.
-- abands=1 has bxyz={{0,0,0}}, a single zero-length vector,
-- so its dot product matrix is {{0}}.
if cf==0 then -- this is for abands=1
cf=1 -- should act like copying mat to omat
end -- if cf
nrow=#mat
ncol=#(mat[1])
omat=mat
omat={}
for rown=1,nrow do
omat[rown]={}
for coln=1,ncol do
omat[rown][coln]=mat[rown][coln]/cf
end -- for coln
end -- for rown
return omat
end -- normmat()
function listmat(mat,str,dec)
local rown,coln,nrow,ncol,ostr
nrow=#mat
ncol=#(mat[1])
print(str..':')
for rown=1,nrow do
ostr=' '
for coln=1,ncol do
ostr=(ostr..roundn(mat[rown][coln],dec)..' ')
end -- for coln
print(ostr)
end -- for rown
end -- listmat()
-- Below finds matrix product pmat
-- of the matrices amat and bmat:
-- pmat = amat x bmat.
--
-- amat should be nrow rows x nmid cols.
-- bmat should be nmid rows x ncol cols.
-- pmat will be nrow rows x ncol cols.
function matxmat(amat,bmat)
local rown,midn,coln,nrow,nmid,nmid2,ncol,tot
local pmat={}
nrow=#amat
nmid=#(amat[1])
nmid2=#bmat
ncol=#(bmat[1])
if nmid~=nmid2 then
print(string.format(' Error in matxmat: amat=%dx%d bmat=%dx%d dimensions are wrong.',nrow,nmid,nmid2,ncol))
end -- if nmid
for rown=1,nrow do
pmat[rown]={}
for coln=1,ncol do
tot=0
for midn=1,nmid do
tot=tot+amat[rown][midn]*bmat[midn][coln]
end
pmat[rown][coln]=tot
end -- for coln
end -- for rown
return pmat
end -- matxmat()
-- Below converts matrix mat
-- into its transpose matrix tmat.
function xpose(mat)
local rown,coln,nrow,ncol
local tmat={}
nrow=#mat
ncol=#(mat[1])
-- mat is nrow rows x ncol cols
-- tmat will be ncol rows x nrow cols
for coln=1,ncol do
tmat[coln]={}
for rown=1,nrow do
tmat[coln][rown]=mat[rown][coln]
end -- for rown
end -- for coln
return tmat
end -- xpose()
function copymat(mat)
local rown,coln,nrow,ncol
local omat={}
nrow=#mat
ncol=#(mat[1])
for rown=1,nrow do
omat[rown]={}
for coln=1,ncol do
omat[rown][coln]=mat[rown][coln]
end -- for coln
end -- for rown
return omat
end -- copymat()
function trunc3(val)
return math.floor(val*1000)/1000
end -- trunc3()
function roundn(val,n)
local cf=10^n
-- n is power of 10 to round to:
-- -2 rounds to 0.01's, -1 to 0.1's, 0 to 1's, 1 to 10's, 2 to 100's
-- n rounds to 10^n's or cf's
return cf*round0(val/cf)
end -- roundn()
function round0(val)
local res=math.floor(val)
if val-res>=0.5 then
res=res+1
end -- if val
return res
end -- round0()
Main()
