Code
max_rebuild_length = 5 -- Maximum rebuild length.
min_rebuild_length = 5 -- - Must be less than or equal to max_rebuild_length. Suggest not going below 4.
number_of_rebuilds = 8 -- It'll look for the best score from this number of rebuilds. Change to taste.
start_idx = 0 -- if this is a residue id, use it as the start point. If not
-- use a semi-random number computed from the score
ss_control = 3
-- ss_control == 1 Only rebuilds loops
-- ss_control == 2 Rebuilds everything. Helices and sheets will be temporarily converted to loops
-- ss_control == 3 Rebuild loops but allows a single non-loop residue in the rebuild segment.
-- they'll be temporarily converted to loops
loop_forever = true -- just keep going. Good for overnight running.
min_residue = 1 -- Change these if you want to rebuild part of a protein only. For example, to
-- rebuild 53-64, set min_residue to 53 and max_residue to 64. To restore, set to
-- 1 and 9999 (or any large number) respectively
max_residue = 999
kClashImportanceDuringRebuild = 1 -- Might try reducing this to get more a different rebuild than the highest scoring -- one, which might overstate
-- the importance of clashing which can be remedied by a quick shake.
-- If, after rebuild, shake, wiggle, etc. , the score is close enough to the starting score, attempt some quakes.
-- close enough is defined by the entry in kQuakingThresholds indexed by the size of segment being rebuilt
kQuakingThresholds = { 8 , 8 , 8 , 15 , 25 , 30 , 35 }
nQuakingThresholds = 7
-- Use of quicksave slots
kOriginalStructureOrNewBest = 1
kBestPostRebuild = 2
-- Stores the best non-improved structure in Quicksave slot 10 (but can change this below) . The structure that's saved
-- is the best before WS etc, so if this script fails to produce (or even if it does) you can load the structure in
-- slot 10 and try WS in different ways.
kBestNonImprovedStructure = 10
original_secondary_structure = {}
n_residues = 0
curr_score = 0
multiplier = 1
-----
-- The Local Quake bit. Copied from LQ1.2 for expediency
-- Bands radiate from a segment to spatially close neighbours
-- Foregoes Fuses and generally tries to optimize for speed
use_ranked_score = 1 -- if 1, use the score being used for ranking (best normally). If 0, use the
-- old style scoring system
max_bands_per_residue = 6 -- Increase this for more bands at each segment
start_idx = 0 -- Set to a valid residue to start with that residue: otherwise
-- start with a semi-random one,
second_offset = 3 -- Set to 0 for a single centre. This is the most likely one to change.
opposite_centre = false -- overrides second_offset if true
target_wiggle_reduction = 120 -- On wiggle, try and reduce the score by this amount after creating the bands.
-- If it's wildly out, reduce/increase band strength
-- for the next try (don't try to adjust for the current scramble).
max_segment_length = 20 -- Don't create bands greater than this length. Set to 999 or something if you're
-- not worried about the potential for long bands creeping in,
min_exclusion = 999 -- These can be set to establish an "exclusion zone", within which no bands will
-- be created. This might be used to prevent banding to a wayward loop.
-- If min_exclusion > max_exclusion no exclusion will be set.
max_exclusion = 1
n_quakes = 10 -- A reasonable number
band_to_minima = true -- Band to minima or maxima in the distance. See bands_from_segment for meaning.
do_second_shake = false -- After creating the bands and pulling with them do a Shake/Wiggle/Shake
-- if do_second_shake is false (best I think) don't do the last shake.
-- Speeds things up quite a bit
scale_band_strength = true -- If true, band strength will be scaled inversely with length
monit = false
---------- A few more globals
candidate_ids_start = {}
candidate_ids_end = {}
candidate_distances = {}
n_candidates = 0
best_score = 0
function GetScore ()
if ( use_ranked_score == 1 ) then
return current.GetScore ()
else
return current.GetEnergyScore ()
end
end
function Coprime ( n )
if ( n < 31 ) then
return 1
elseif ( n >= 31*37 ) then
return 1
elseif ( n%31 ~= 0 ) then
return 31
else
return 37
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 roundto3 ( i )
-- printing convenience
return i - i % 0.001
end
function delete_all_bands ()
band.DeleteAll ()
end
function ShellSort ( ids , distances , n )
-- Adapted from Numerical Recipes in C
local inc = 1
repeat
inc = inc * 3 + 1
until inc > n
repeat
inc = inc / 3
inc = inc - inc % 1
for i = inc + 1 , n do
v = distances [ i ]
w = ids [ i ]
j = i
flag = false
while ( flag == false and distances [ j - inc ] > v ) do
distances [ j ] = distances [ j - inc ]
ids [ j ] = ids [ j - inc ]
j = j - inc
if ( j <= inc ) then
flag = true
end
end
distances [ j ] = v
ids [ j ] = w
end
until inc <= 1
end
function test_for_minimum ( i , seg , seg_distances )
if ( ( i ~= seg ) and
( seg_distances [ i ] <= seg_distances [ i - 1 ] ) and
( seg_distances [ i ] <= seg_distances [ i - 2 ] ) and
( seg_distances [ i ] <= seg_distances [ i + 1 ] ) and
( seg_distances [ i ] <= seg_distances [ i + 2 ] ) ) then
return true
else
return false
end
end
function test_for_maximum ( i , seg , seg_distances )
if ( ( seg_distances [ i ] >= seg_distances [ i - 1 ] ) and
( seg_distances [ i ] >= seg_distances [ i - 2 ] ) and
( seg_distances [ i ] >= seg_distances [ i + 1 ] ) and
( seg_distances [ i ] >= seg_distances [ i + 2 ] ) ) then
return true
else
return false
end
end
function get_band_segment_data ( seg )
seg_distances = {}
-- Get the function that looks like this:
-- x = residue number
-- y = distance of that segment from seg
-- ignore the degenerate minimum
n_segments = 0
segment_ids = {}
segment_distances = {}
for i = 1, n_residues do
segment_distances [ i ] = structure.GetDistance ( seg , i )
end
for i = 3 , n_residues - 3 do
k = segment_distances [ i ]
if ( ( band_to_minima == true ) and
( test_for_minimum ( i , seg , segment_distances ) == true ) ) then
n_segments = n_segments + 1
segment_ids [ n_segments ] = i
segment_distances [ n_segments ] = segment_distances [ i ]
elseif ( ( band_to_minima == false ) and
( test_for_maximum ( i , seg , segment_distances ) == true ) ) then
n_segments = n_segments + 1
segment_ids [ n_segments ] = i
segment_distances [ n_segments ] = segment_distances [ i ]
end
end
ShellSort ( segment_ids , segment_distances , n_segments )
if ( n_segments >= max_bands_per_residue ) then
n_c = max_bands_per_residue
else
n_c = n_segments
end
for i = 1 , n_c do
if ( ( segment_distances [ i ] <= max_segment_length ) and
( seg < min_exclusion or seg > max_exclusion ) and
( segment_ids [ i ] < min_exclusion or segment_ids [ i ] > max_exclusion ) ) then
n_candidates = n_candidates + 1
candidate_ids_start [ n_candidates ] = seg
candidate_ids_end [ n_candidates ] = segment_ids [ i ]
candidate_distances [ n_candidates ] = segment_distances [ i ]
end
end
most_distant_minimum = segment_ids [ n_segments ]
end
function create_bands ()
if ( scale_band_strength == true ) then
total_band_length = 0
for i = 1 , n_candidates do
total_band_length = total_band_length + candidate_distances [ i ]
end
average_band_length = total_band_length / n_candidates
end
for i = 1 , n_candidates do
band.AddBetweenSegments ( candidate_ids_start [ i ] , candidate_ids_end [ i ] )
if ( scale_band_strength == true ) then
band.SetStrength ( i , band_strength * average_band_length / candidate_distances [ i ] )
else
band.SetStrength ( i , band_strength )
end
end
end
function Quake ( seg )
selection.SelectAll ()
init_score = GetScore ()
structure.WiggleSelected ( 1 )
score_after_wiggle = GetScore ()
score_reduction_after_wiggle = init_score - score_after_wiggle
-- Adaptively change band strength based on how much the previous wiggle overshot/undershot the target. Could
-- be made much more elaborate.
if ( score_reduction_after_wiggle > target_wiggle_reduction * 1.2 ) then
band_strength = band_strength * ( 0.85 + 0.1 *pseudorandom ( score_after_wiggle ) )
elseif ( score_reduction_after_wiggle < target_wiggle_reduction * 0.8 ) then
band_strength = band_strength + ( 1 - band_strength ) * ( 0.05 + 0.1 *pseudorandom ( score_after_wiggle ) )
end
delete_all_bands ()
structure.MutateSidechainsSelected ( 1 )
gas = GetScore () - score_after_wiggle
if ( monit ) then
print ( "law " .. roundto3 ( init_score - score_after_wiggle ) .. " gas " .. roundto3 ( gas ) .. " band_strength " .. roundto3 ( band_strength ) )
end
structure.WiggleSelected ( 8 )
if ( do_second_shake == true ) then
score_a2w = GetScore ()
structure.MutateSidechainsSelected ( 1 )
score_a2s = GetScore ()
if ( score_a2s - score_a2w > 1 ) then
structure.WiggleSelected ( 4 )
end
end
end
function local_quake ( start_idx )
band_strength = 0.5
behavior.SetClashImportance ( 1 )
recentbest.Save ()
best_score_lq = GetScore ()
inc = Coprime ( n_residues )
idx_lq = start_idx
for i = 1, n_quakes do
recentbest.Restore ()
-- Need to record any improvements from the previous iteration,.
score_lq = GetScore()
if ( score_lq > best_score_lq ) then
best_score_lq = score_lq
print ( " Gain to " .. roundto3 ( score_lq - score_lq % 0.01 ) )
end
delete_all_bands ()
n_candidates = 0
get_band_segment_data ( idx_lq )
if ( opposite_centre == true ) then
second_idx = most_distant_minimum
get_band_segment_data ( most_distant_minimum )
elseif ( second_offset ~= 0 ) then
second_idx = idx_lq + second_offset
while ( second_idx > n_residues ) do
second_idx = second_idx - n_residues
end
get_band_segment_data ( second_idx )
end
if ( n_candidates > 0 ) then
create_bands ()
if ( second_offset ~= 0 ) then
print ( " i: " .. i .. " seg ".. idx_lq .. "," .. second_idx )
else
print ( " i: " .. i .. " seg ".. idx_lq )
end
Quake ( idx_lq )
end
idx_lq = idx_lq + inc
if ( idx_lq > n_residues ) then
idx_lq = idx_lq - n_residues
end
end -- for
recentbest.Restore ()
delete_all_bands ()
end
------
function IsSegmentALoop ( first , last )
-- Returns number of non-loop residues (helices or sheets)
n_non_loops = 0
for i = first , last do
if ( structure.GetSecondaryStructure ( i ) ~= "L" ) then
n_non_loops = n_non_loops + 1
end
end
return n_non_loops
end
function ConvertToLoop ( first , last )
for k = first , last do
if ( original_secondary_structure [ k ] ~= "L" ) then
selection.DeselectAll ( )
selection.Select ( k )
structure.SetSecondaryStructureSelected ( "L" )
end
end
end
function RestoreStructure ( first , last )
for k = first, last do
selection.DeselectAll ()
selection.Select ( k )
structure.SetSecondaryStructureSelected ( original_secondary_structure [ k ] )
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.WiggleSelected ( 1 )
behavior.SetClashImportance ( 1 )
structure.WiggleSelected ( 8 )
recentbest.Restore ()
end
function WiggleMutateWiggle ( start_idx , end_idx )
-- returns the score after doing a WMW: shortcircuiting if things don't look promising
score_a = GetScore ()
structure.WiggleSelected ( 15 )
recentbest.Restore ()
score_after_wiggle = GetScore ()
if ( score_after_wiggle - score_a < 10*multiplier and best_score - score_after_wiggle > 30*multiplier ) then
-- The wiggle got stuck and didn't achieve anything
nudge ( score_after_wiggle )
score_after_wiggle = GetScore ()
end
if ( monit ) then
--print ( "After first wiggle (bb) : ".. roundto3 ( best_score - score_after_wiggle ) )
end
if ( best_score - score_after_wiggle > 100*multiplier ) then
-- Not worth continuing
return score_after_wiggle
end
structure.MutateSidechainsSelected ( 1 )
score_after_second_shake = GetScore ()
if ( true ) then --score_after_second_shake - score_after_wiggle > 0.5 ) then
structure.WiggleSelected ( 10 )
recentbest.Restore ()
score_after_second_wiggle = GetScore ()
if ( true ) then -- score_after_second_wiggle - score_after_second_shake < 1 and best_score - score_after_second_wiggle < 50*multiplier ) then
-- The second shake did something but the subsequent wiggle didn't
nudge ( score_after_second_wiggle )
score = GetScore ()
return score
else
return score_after_second_wiggle
end
else
return score_after_second_shake
end
end
function post_rebuild ( score_after_rebuild , start_idx , end_idx )
-- Tries to recover after a rebuild.
if ( best_score - score_after_rebuild > 3000*multiplier ) then
structure.MutateSidechainsSelected ( 2 )
elseif ( best_score - score_after_rebuild > 1500*multiplier ) then
structure.MutateSidechainsSelected ( 1 )
end
-- First try WMW
curr_score = WiggleMutateWiggle ( start_idx , end_idx )
end
function rebuild_seg ( start_idx , end_idx )
-- return codes
--
-- 0 rebuild had no effect: seg probably locked
-- 1 no improvement
-- 2 improvement
selection.DeselectAll ()
selection.SelectRange ( start_idx , end_idx )
-- for some reason, do_local_rebuild ( 1) often leaves the protein unchanged
k = 0
repeat
k = k + 1
behavior.SetClashImportance ( kClashImportanceDuringRebuild )
structure.RebuildSelected ( k )
behavior.SetClashImportance ( 1 )
score_after_rebuild = GetScore ()
until score_after_rebuild ~= best_score or k > 2
if ( score_after_rebuild == best_score ) then
-- If the scores are still equal we probably have a locked segment
return 0
end
recentbest.Save ()
if ( number_of_rebuilds > 1 ) then
behavior.SetClashImportance ( kClashImportanceDuringRebuild )
structure.RebuildSelected ( number_of_rebuilds - 1 )
behavior.SetClashImportance ( 1 )
recentbest.Restore ()
end
-- Now do the post rebuild wiggle and shakes
score_after_rebuild = GetScore ()
selection.SelectAll ()
-- Compute a value delta which will be used later to determine if we are
-- "sufficiently" close to our original best score. If it turns out later that
-- we are, then nudge it around
rebuild_length = end_idx - start_idx + 1
if ( rebuild_length > nQuakingThresholds ) then
delta = kQuakingThresholds [ nQuakingThresholds ] * multiplier -- sanity check for lengthy rebuild segs
else
delta = kQuakingThresholds [ rebuild_length ] * multiplier
end
save.Quicksave ( kBestPostRebuild )
post_rebuild ( score_after_rebuild , start_idx , end_idx )
if ( best_score - curr_score < delta ) then
print ( " cur score " .. roundto3 ( curr_score ) )
local_quake ( start_idx )
end
if ( ss_control > 1 ) then
save.LoadSecondaryStructure ()
end
curr_score = GetScore ()
if ( curr_score > best_score ) then
best_score = curr_score
improvement_made = true
save.Quicksave ( kOriginalStructureOrNewBest )
print ( "Improvement to ".. roundto3 ( best_score ) )
end
if ( curr_score > qs10_best_score and improvement_made == false ) then
recentbest.Save ()
qs10_best_score = curr_score
save.Quickload ( kBestPostRebuild )
-- Bracket the area with frozen segments for easy recognition
selection.DeselectAll ()
if ( start_idx > 1 ) then
selection.Select ( start_idx - 1 )
end
if ( end_idx < n_residues ) then
selection.Select ( end_idx + 1 )
end
freeze.FreezeSelected ( true , false )
save.Quicksave ( kBestNonImprovedStructure )
recentbest.Restore ()
end
if ( improvement_made == true ) then
return 2
else
return 1
end
end
function rebuild ( rebuild_length )
n_possible_segs = n_residues - rebuild_length + 1
inc = Coprime ( n_possible_segs )
-- The point of this Coprime, start_idx business is so we go through the protein in a
-- nonuniform way. Rather than rebuilding 1-4, 2-5 , 3-6 for example, we might do 21 4, 53-57, 3-6
-- or something like that
for i = 1 , n_possible_segs do
improvement_made = false
end_idx = start_idx + rebuild_length - 1
if ( start_idx >= min_residue and end_idx <= max_residue and
( ss_control == 2 or
( ss_control == 1 and IsSegmentALoop ( start_idx , end_idx ) == 0 ) or
( ss_control == 3 and IsSegmentALoop ( start_idx , end_idx ) < 2 ) ) ) then
rb = rb + 1
print ( "rb ".. rb .. " " .. start_idx .."-" .. end_idx .. " (" ..i .. "/" .. n_possible_segs .. ")" )
save.Quickload ( kOriginalStructureOrNewBest )
if ( ss_control > 1 ) then
ConvertToLoop ( start_idx , end_idx )
end
best_score = GetScore ()
k = rebuild_seg ( start_idx, end_idx )
if ( ss_control > 1 ) then
save.LoadSecondaryStructure ()
end
end
start_idx = start_idx + inc
if ( start_idx > n_possible_segs ) then
start_idx = start_idx - n_possible_segs
end
end -- for i
end
-- Tidy up if necessary
band.DeleteAll ()
n_residues = structure.GetCount ()
if ( ss_control > 1 ) then
for i = 1, n_residues do
original_secondary_structure [ i ] = structure.GetSecondaryStructure ( i )
end
end
save.SaveSecondaryStructure ()
rb = 0
save.Quicksave ( kOriginalStructureOrNewBest )
behavior.SetClashImportance ( 1 )
best_score = GetScore ()
print ( "Start score : " .. roundto3 ( best_score ) )
if ( start_score == 0 ) then
print ( "Suggest quitting as starting score < = 0" )
end
if ( ss_control == 1 ) then
print ( "Loops only" )
elseif ( ss_control == 2 ) then
print ( "Loops, sheets, and helices" )
else
print ( "Loops with one non-loop allowed" )
end
if ( max_residue > n_residues ) then
max_residue = n_residues
end
if ( min_residue < 1 ) then
min_residue = 1
end
if ( max_residue <= min_residue ) then
print ( "Rebuild range error" )
loop_forever = false
end
multiplier = current.GetExplorationMultiplier ()
if ( multiplier <1 ) then
multiplier = 1
end
if ( true ) then -- min_residue > 1 or max_residue < n_residues ) then
print ( "Rebuild range " .. min_residue .. " to " .. max_residue )
else
print ( "Rebuild all" )
end
qs10_best_score = 0
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
repeat
for rebuild_length = max_rebuild_length , min_rebuild_length , -1 do
rebuild ( rebuild_length )
end
until loop_forever == false
save.Quickload ( kOriginalStructureOrNewBest )
