Code
----
---- For each rotamer position of each sidechain, does a flip, freeze, shake, unfreeze, wiggle, shake, wiggle.
---- Skips the rest of the operation sequence and continues with the next rotamer position if the first shake does not improve the score.
----
-- if this is a residue id, use it as the start point. If not
-- use a semi-random number computed from the score
start_idx = 0
function Dialog_ChooseLevel()
ask = dialog.CreateDialog("Level?")
--ask.level = dialog.AddSlider( "Level", 2, 1, 4, 0 )
--ask.ok = dialog.AddButton( "OK", 0 )
ask.level1 = dialog.AddButton( "Level 1", 1 )
ask.level2 = dialog.AddButton( "Level 2", 2 )
ask.level3 = dialog.AddButton( "Level 3", 3 )
ask.level4 = dialog.AddButton( "Level 4", 4 )
return dialog.Show( ask )
end
-- level 1: Flip Leu and Ile only. Flipping these tend to be the most productive, although Glu is a close third.
-- level 2: Default. Best trade-off between potential gain and speed.
-- level 3: Tries more positions. Unlike the lower 2 levels, does not require that on the first shake an additional side chain be flipped.
-- level 4: Foregoes the freeze-shake-unfreeze bit of the algorithm entirely. This can result in severe disruption to the protein structure.
-- Not terribly useful except perhaps in the opening.
level = 2
manual_level_set = false
-- post_scramble = 1: no sidechain scramble. Won't leak walk juice.
-- post_scramble = 2. does a brief scramble iff a flip gets within 2 (see kScrambleThreshold) points of the original.
post_scramble = 2
kScrambleThreshold = 2
-- verbosity level:
-- 0: prints no progress output (not recommended)
-- 1: prints the residue the recipe is working on
-- 2: also prints the rotamer snap indices
verbosity = 2
-- Amino acids. Hydrophobics first.
single_letter_codes = { "g","a","v","l","i","m","f","w","p","s","t","c","y","n","q","d","e","k","r","h"}
three_letter_codes = { "Gly","Ala","Val","Leu","Ile","Met","Phe","Trp","Pro","Ser","Thr","Cys","Tyr","Asn","Gln","Asp","Glu","Lys","Arg","His"}
-- Use of save.Quicksave slots: change if desired
kBest = 1 -- used to hold the initial structure and any improved structures
kClosest = 10 -- holds the best scoring structure that does not exceed that in kBest
function get_aa3 ( i )
k = structure.GetAminoAcid ( i )
for j = 1, 20 do
if ( k == single_letter_codes [ j ] ) then
return three_letter_codes [ j ]
end
end
return "Unk"
end
--
-- returns a positive integer which has no common prime factor with the given integer
-- more mathematically: for each integer n: Coprime(n) > 0 such that gcd(n,Coprime(n)) == 1
-- something of a hack, but it seems unlikely we'll ever have puzzles with more than 1146 residues
--
function Coprime( n )
if( n == 0 ) then
-- every integer is a factor of 0
-- but we only need gcd(0,Coprime(0)) == 1
return 1
elseif( n < 0 ) then
n = -n
end
if( n%31 ~= 0 ) then
return 31
elseif( n%37 ~= 0 ) then
return 37
else
return 1
end
end
function pseudorandom ( score )
-- Returns a pseudorandom number >= 0 and < 1.
-- Based on the fractional part of the score
if ( score >= 0 ) then
return score % 1
else
return (-score ) % 1
end
end
function nudge( score )
selection.SelectAll ()
recentbest.Save ()
behavior.SetClashImportance ( 0.2 + 0.2 * pseudorandom ( score ) )
structure.WiggleAll ( 1 )
behavior.SetClashImportance ( 1 )
structure.WiggleAll ( 8 )
recentbest.Restore ()
end
function quick_nudge( score )
selection.SelectAll ()
recentbest.Save ()
behavior.SetClashImportance ( 0.1 + 0.1 * pseudorandom ( score ) )
structure.ShakeSidechainsAll ( 1 )
behavior.SetClashImportance ( 1 )
structure.WiggleAll ( 8 )
structure.ShakeSidechainsAll ( 1 )
structure.WiggleAll (1,false,true)
recentbest.Restore ()
end
function FlipAllSidechains()
-- go through the protein in a less than regular way
segment_count = structure.GetCount()
increment = Coprime( segment_count )
while( increment >= segment_count ) do
increment = increment - segment_count
end
segment = start_idx
for i = 1,segment_count do
-- flip the sidechain
amino_acid = get_aa3( segment )
if( verbosity >= 1 ) then
print( i, "/", segment_count .. ":", amino_acid, segment )
end
if( leu_ile_only == false or amino_acid == "Leu" or amino_acid == "Ile" ) then
FlipSidechain( segment )
end
-- move on to next segment
segment = segment + increment
if( segment > segment_count ) then
segment = segment - segment_count
end
end
end
function FlipSidechain( segment )
r = 1
improvement_made = false
-- Lysines and arginines can have 40 or rotamers and its wasteful to go through them all: flipping the
-- sidechains of those amino acids rarely works anyway
-- Till there's something in the API that allows identification of the residue, restrict the number of attempts
n_attempts = 0
while ( r <= rotamer.GetCount(segment) and improvement_made == false and n_attempts < max_attempts ) do
rotamer.SetRotamer( segment, r )
score_after_snap = current.GetEnergyScore()
snap_loss = best_score - score_after_snap
if ( snap_loss < 0 ) then
-- Can conceivably happen
save.Quicksave( kBest )
best_score = score_after_snap
improvement_made = true
print ( "Immediate improvement at segment", segment )
print ( "best score:", best_score )
n_attempts = n_attempts + 1
elseif ( snap_loss > min_snap_loss ) then
-- If snapping changes it by less than 100, this technique is unlikely to produce a gain
n_attempts = n_attempts + 1
if( verbosity >= 2 ) then
print( "", "attempt " .. n_attempts .. ":", "rotamer", r, "/", rotamer.GetCount( segment ) )
end
gain_from_shake = 0
if( level < 4 ) then
selection.DeselectAll()
selection.Select( segment )
freeze.FreezeSelected ( false , true )
selection.SelectAll()
structure.ShakeSidechainsAll ( 1 )
score_after_shake = current.GetEnergyScore()
gain_from_shake = score_after_shake - score_after_snap
freeze.UnfreezeAll()
else
score_after_shake = score_after_snap
gain_from_shake = 0
end
if( gain_from_shake > lowest_allowable_shake_gain and best_score - score_after_shake < max_shake_loss ) then
--print ( "score after shake " , best_score - score_after_shake )
structure.WiggleAll ( 15 )
score_after_wiggle = current.GetEnergyScore ()
structure.ShakeSidechainsAll ( 1 )
score_after_second_shake = current.GetEnergyScore ()
if ( score_after_second_shake - score_after_wiggle > 0.5 ) then
structure.WiggleAll ( 10 )
score_after_second_wiggle = current.GetEnergyScore ()
if ( score_after_second_wiggle - score_after_second_shake < 1 and
best_score - score_after_second_wiggle < 30 ) then
nudge ( score_after_second_wiggle )
score = current.GetEnergyScore ()
else
score = score_after_second_wiggle
end
else
score = score_after_second_shake
end
if( best_score - score < kScrambleThreshold and post_scramble == 2 ) then
quick_nudge ( score )
score = current.GetEnergyScore()
end
if ( score > best_score ) then
structure.ShakeSidechainsAll( 1 )
save.Quicksave ( kBest )
best_score = current.GetEnergyScore()
--print( "best score ", best_score , " snap loss " , snap_loss , " shake gain " , gain_from_shake )
print( "new best score:", best_score )
improvement_made = true;
end
if ( improvement_made == false ) then
if ( score > score_in_qs10 ) then
save.Quicksave ( kClosest )
score_in_qs10 = score
--print ( "Save.Quicksave 10 " , score )
end
end
end -- if shake improved things
end -- if not current rotamer position
r = r + 1
save.Quickload( kBest )
end -- loop over all rotamers
end
---- main script ----
if( manual_level_set ) then
level = Dialog_ChooseLevel()
end
n_residues = structure.GetCount()
behavior.SetClashImportance( 1 )
-- tidy up the protein
freeze.UnfreezeAll()
band.DisableAll()
--?? was there any reason for using this loop instead of band.DisableAll() -- Wipf, 08. Mar 2019
--band_count = band.GetCount()
--for i = band_count, 1, -1 do
-- band.Disable ( i )
--end
save.Quicksave( kBest )
best_score = current.GetEnergyScore()
score_in_qs10 = 0
print( "best score ", best_score )
if( level == 1 ) then
leu_ile_only = true
else
leu_ile_only = false
end
if( level == 1 or level == 2 ) then
lowest_allowable_shake_gain = 0.5
else
lowest_allowable_shake_gain = -0.5
end
if( level == 1 or level == 2 or level == 4 ) then
max_attempts = 3
else
max_attempts = 6
end
min_snap_loss = 60
if( level == 1 or level == 2 ) then
max_shake_loss = 3000
else
max_shake_loss = 100000
end
if( start_idx < 1 or start_idx > n_residues ) then
r = pseudorandom ( best_score )
start_idx = 1 + n_residues * r
start_idx = start_idx - start_idx % 1
end
FlipAllSidechains()
