# Matthias 343
### __defs.top
# Define additional TOP variables needed for MODELLER:

DEFINE_INTEGER VARIABLES = STARTING_MODEL ENDING_MODEL RSTRS_REFINED
DEFINE_INTEGER VARIABLES = LOOP_STARTING_MODEL LOOP_ENDING_MODEL
DEFINE_INTEGER VARIABLES = MAX_ITERATIONS_STORE WRITE_INTERMEDIATES
DEFINE_INTEGER VARIABLES = IREPEAT REPEAT_OPTIMIZATION EXIT_STAGE
DEFINE_INTEGER VARIABLES = CREATE_RESTRAINTS REFINE_HOT_ONLY
DEFINE_INTEGER VARIABLES = MAX_VAR_ITERATIONS FINAL_MALIGN3D INITIAL_MALIGN3D
DEFINE_INTEGER VARIABLES = DO_LOOPS ID1B
DEFINE_REAL    VARIABLES = VIOL_REP_STORE MAX_MOLPDF
DEFINE_REAL    VARIABLES = MAX_CA-CA_DISTANCE MAX_N-O_DISTANCE
DEFINE_REAL    VARIABLES = MAX_SC-SC_DISTANCE MAX_SC-MC_DISTANCE
DEFINE_STRING  VARIABLES = MODEL MODEL2 CODE CODE2 ALNFILE MODEL2_FIT
DEFINE_STRING  VARIABLES = CSRFILE KNOWNS SCHFILE FINAL_MODEL
DEFINE_STRING  VARIABLES = GENERATE_METHOD RAND_METHOD MD_LEVEL
DEFINE_STRING  VARIABLES = SEGFILE PDB_EXT TOPLIB PARLIB FAMILY FIT_IN_REFINE
DEFINE_STRING  VARIABLES = ATOM_IDS1 ATOM_IDS2 OUTPUT2
DEFINE_STRING  VARIABLES = LOOP_CSRFILE LOOP_INI_MODEL
DEFINE_STRING  VARIABLES = LOOP_MD_LEVEL LOOP_INI_MODEL
DEFINE_STRING  VARIABLES = LOOP_MODEL
DEFINE_STRING VARIABLES = SEG1 SEG2


# default values for options in comparative modeling by MODELLER:

SET STARTING_MODEL= 1    # the index of the first model; 
                         # determines how many models are calculated;

SET ENDING_MODEL  = 13   # the index of the last model; 
                         # determines how many models are calculated;

SET DEVIATION     = 10.0  # the amount of randomization of the initial model
                         # must be > 0 if different final models are wanted;

SET DO_LOOPS = '1'       # whether or not to do automatic loop refinement
                         # for each model *.B???????? (0 for no, 1 for yes)
                         # can rely on automatic loop definition or 
                         # re-define select_loop_atoms routine.
# SET LOOP_STARTING_MODEL = 1 # how many loop models to generate for
# SET LOOP_ENDING_MODEL = 25  # each *.B????????
# SET LOOP_MD_LEVEL = 'refine_3' # the same as for MD_LEVEL, but for loops

# Options that are not changed frequently:

SET LIBRARY_SCHEDULE = 1      # 1 ... thorough var target func schedule
                              # 4 ... faster var target func schedule
SET MAX_VAR_ITERATIONS = 200  # maximal numb of iterations for the cycles
                              # of the variable target function method
SET MD_LEVEL = 'refine_3'     # what kind of optimization is done after
                              # the variable target function method:
                              # 'nothing'            ... nothing;
                              # 'refine_1'            ... very fast MD annealing;
                              # 'refine_2'            ... fast MD annealing;
                              # 'refine_3'            ... slow MD annealing;
                              # 'refine_4'            ... very slow MD annealing;
                              # 'refine_5'            ... very slow/large dt MD annealing;
SET REFINE_HOT_ONLY = 0       # 1 ... select and optimize only HOT atoms in refine;
                              # 0 ... select and optimize all atoms in refine;
                              # usually about half of the atoms are hot; in such cases,
                              # 0 is faster for sequences longer than about 100 aa
                              # because a faster non-bonded pairs algorithm can be used.
SET RSTRS_REFINED = 2         # the types of restraints used to define
                              # hot spots when MD_LEVEL <> 'nothing':
                              # 0 ... stereochemistry only;
                              # 1 ... stereochemistry and dihedral;
                              # 2 ... all restraints;
SET EXIT_STAGE    = 0         # 0 ... no effect;
                              # 1 ... exit without any optimization after
                              #       restraints and an initial model are 
                              #       calculated (more efficient than
                              #       REPEAT_OPTIMIZATION=0);
                              # 2 ... exit after the initial model is calculated
                              #       (restraints are not calculated)
SET REPEAT_OPTIMIZATION = 1   # how many times the whole optimization schedule (variable
                              # target function method and refinement) is repeated 
                              # for each initial model;
SET TRACE_OUTPUT =  15        # every which CG or MD cycle is reported;

SET MAX_MOLPDF   = 100E3      # abort optimization of the current model if 
                              # the molecular pdf is larger than this and
                              # continue with the next model;
SET TOPLIB = '${LIB}/top_heav.lib'  # topology library (all non-hydrogen atoms);
SET TOPOLOGY_MODEL = 3              # corresponding topology model;
SET PARLIB = '${LIB}/par.lib'       # parameters library;
SET WRITE_INTERMEDIATES = 0         # 0 ... do not write out intermediate 
                                    #       atom files during optimization;
                                    # 1 ... write out intermediate atom files;
SET INITIAL_MALIGN3D = 0            # 0 ... do not do MALIGN3D before
                                    #       TRANSFER_XYZ
                                    # 1 ... do that.
SET FINAL_MALIGN3D = 1              # 0 ... do not do MALIGN3D and write 
                                    #       superposed templates & models
                                    #       at the end of 'model'
                                    # 1 ... do that.
SET GENERATE_METHOD= 'transfer_xyz' # how to build the initial model: 
                                    # 'generate_xyz' from internal coordinates 
                                    #                and write them to a file;
                                    # 'transfer_xyz' from template coordinates
                                    #                and write them to a file;
                                    # 'read_xyz'     read coordinates from 
SET RAND_METHOD   = 'randomize_xyz' # a method to perturb the initial model:
                                    # 'randomize_dihedrals' ... uses DEVIATION in degrees;
                                    # 'randomize_xyz'       ... uses DEVIATION in angstroms;
                                    # 'nothing'
SET CREATE_RESTRAINTS =  1          # 0 ... read the restraints from a file;
                                    # 1 ... make the restraints and write them to a file 
                                    #       before reading them for the optimization; in 
                                    #       addition to the default restraints, the TOP
                                    #       routine 'special_restraints', which may be 
                                    #       re-defined in the user TOP file, is called for
                                    #       any user defined restraints that are then also 
                                    #       written to the same file.
SET SPLINE_ON_SITE = on             # on  ... convert some restraints into splines
                                    # off ... no conversion
SET OUTPUT_CONTROL = 1 1 1 1 1      # write real_output, notes, warnings, errors, dynmem
# SET OUTPUT_CONTROL = 1 0 0 1 0      # write real_output, notes, warnings, errors, dynmem

# Set maximal values for various distance restraints:
SET MAX_CA-CA_DISTANCE = 14.0
SET MAX_N-O_DISTANCE   = 11.0
SET MAX_SC-MC_DISTANCE =  5.5
SET MAX_SC-SC_DISTANCE =  5.0

# Routine 'user_after_single_model' can be redefined to do whatever at the end
# of each model calculation (e.g. comparison with X-ray structure).

# To write out reports on individual optimizations:
SET OUTPUT = 'NO_REPORT SHORT'

# The alignment file format (I/O):
SET ALIGNMENT_FORMAT = 'PIR'

# The extension added to all *.Bxxxxnn filenames:
# SET PDB_EXT = ' '
SET PDB_EXT = '.modlr.pdb'

# to prevent SUPERPOSE in refine() if molecules are too small:
SET FIT_IN_REFINE = 'NO_FIT'

# To enable default filename generation if not explicitly defined:
#SET MODEL   = 'undefined'
SET MODEL = 'hERG_ali4.ini'
SET CSRFILE = 'undefined'

# Call this routine before calling 'model' if you want real fast optimization
SUBROUTINE ROUTINE = 'very_fast'

  SET STARTING_MODEL = 1
  SET ENDING_MODEL = 10
# SET STARTING_MODEL = 1
# SET ENDING_MODEL = 1
  SET MAX_CA-CA_DISTANCE = 10.0
  SET MAX_N-O_DISTANCE   =  6.0
  SET MAX_SC-MC_DISTANCE =  5.0
  SET MAX_SC-SC_DISTANCE =  4.5
  # Note that all models will be the same if you do not change RAND_METHOD
  SET RAND_METHOD = 'randomize_xyz'
  SET MAX_VAR_ITERATIONS = 50
  SET LIBRARY_SCHEDULE = 7
  SET MD_LEVEL = 'refine_1'

  RETURN
END_SUBROUTINE

### get-model.top
# PRIMER: STEP 5

SET ALNFILE  = '../../ALIGNMENT/Thias1.pir'           # alignment filename
SET KNOWNS   = 'KCSA_CUT'       # codes of the templates
SET SEQUENCE = 'HERG'               # code of the target
SET ATOM_FILES_DIRECTORY = './:../../PDB/' # directories for input atom files
SET STARTING_MODEL= 1               # index of the first model 
SET ENDING_MODEL  = 13              # index of the last model
                                    # (determines how many models to calculate)
SET DEVIATION     = 10.0          # e' gia' definita all'inizio  # have to be >0 if more than 1 model
SET RAND_SEED     = -11542          # to have different models from another TOP file

CALL ROUTINE = 'model'              # do homology modelling

### __model.top
# Purpose: calculate model(s) using some restraints, possibly obtained
#          from the alignment of a sequence with related structure(s).
# Input  : ALNFILE, KNOWNS, SEQUENCE, DATABASE_DIRECTORY,
#           DEVIATION, TOPLIB, PARLIB, GENERATE_METHOD, RAND_METHOD, 
#           MD_LEVEL
# Output : stereochemical and homology restraints in file CSRFILE
#          schedule in file SCHFILE, final model in file FINAL_MODEL,
#          debugging output in .D
#          MODEL, CSRFILE, SCHFILE, FINAL_MODEL, REFINED_MODEL
SUBROUTINE ROUTINE = 'model'
  # get the filenames from SEQUENCE:
  CALL ROUTINE = 'getnames'
  # getting restraints and the initial model from the alignment:
  CALL ROUTINE = 'homcsr'
  # exit early?
  IF ARGUMENTS = EXIT_STAGE 1, OPERATION = 'EQ', THEN = 'RETURN'
  # Read all restraints once for the whole job (except when loops are
  # done when restraints are read for each *.B9999???? model):
  CALL ROUTINE = 'rd_restraints'
  # getting model(s) (topology library must be in memory; ensured
  # now by one of the three GENERATE_METHOD routines):
  CALL ROUTINE = 'multiple_models'

  IF ARGUMENTS=FINAL_MALIGN3D 1, OPERATION='NE', THEN='GO_TO NO_MALIGN3D'
      CALL ROUTINE = 'fit_models_on_template'
     
      IF ARGUMENTS=DO_LOOPS 1, OPERATION='NE', THEN='GO_TO NO_MALIGN3D'
        DO ID2 = STARTING_MODEL ENDING_MODEL 1
          CALL ROUTINE = 'fit_loop_models'
        END_DO
      
  LABEL NO_MALIGN3D

RETURN
END_SUBROUTINE

# _fit.pdb extension for the output whole models, superposed on the template
# Input: ALNFILE, KNOWNS, SEQUENCE, ROOT_NAME, FILE_EXT
#
SUBROUTINE ROUTINE = 'fit_models_on_template'
      READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS SEQUENCE
      SET FILE_ID = '.B'
      EXPAND_ALIGNMENT EXPAND_CONTROL = 9999 9999 STARTING_MODEL ENDING_MODEL 0
      # To take care of the '.' in segment specs:
      WRITE_ALIGNMENT FILE = '.tmp.ali', ALIGNMENT_FORMAT = 'PIR'
      READ_ALIGNMENT FILE = '.tmp.ali', ADD_SEQUENCE  = OFF
      SYSTEM COMMAND = 'rm .tmp.ali'
      # superpose all on the first structure in the alignment:
      COMPARE FIT = on
      MALIGN3D GAP_PENALTIES_3D = 0 3.0, WRITE_WHOLE_PDB = off, WRITE_FIT = on, ;
               FIT = off, CURRENT_DIRECTORY = on
     
RETURN
END_SUBROUTINE

# _fit.pdb extension for the output loop models, superposed on their whole .B*_fit.pdb structure
# Input: ID2
#
SUBROUTINE ROUTINE = 'fit_loop_models'
    READ_MODEL FILE = 'default', FILE_EXT = '_fit.pdb', FILE_ID = '.B', ROOT_NAME = SEQUENCE, ID1 = 9999
    SEQUENCE_TO_ALI ALIGN_CODES = SEQUENCE, ATOM_FILES = LOOP_MODEL, ADD_SEQUENCE = off
    SEQUENCE_TO_ALI ADD_SEQUENCE = on, ALIGN_CODES = ALIGN_CODES SEQUENCE, ATOM_FILES = ATOM_FILES LOOP_MODEL
    DO ID1 = LOOP_STARTING_MODEL,LOOP_ENDING_MODEL,1
      READ_MODEL2 FILE = 'default' , FILE_EXT = '', FILE_ID = '.BL'
      PICK_ATOMS PICK_ATOMS_SET = 1, ATOM_TYPES = 'CA'
      SUPERPOSE FIT = 'ON'
      WRITE_MODEL2 FILE = 'default', FILE_EXT = '_fit.pdb', FILE_ID = '.BL'
    END_DO
    RETURN
END_SUBROUTINE

### __getnames.top
# Purpose: To create the filenames for various file types
# Input  : SEQUENCE (1 entry)
#
SUBROUTINE ROUTINE = 'getnames'

  # Redefine only if undefined to allow definition by a user:
  STRING_IF STRING_ARGUMENTS = MODEL 'undefined', OPERATION = 'EQ', THEN = ;
    'STRING_OPERATE OPERATION = CONCATENATE, STRING_ARGUMENTS = SEQUENCE .ini, RESULT = MODEL'

  STRING_IF STRING_ARGUMENTS = CSRFILE 'undefined', OPERATION = 'EQ', THEN = ;
    'STRING_OPERATE OPERATION = CONCATENATE, STRING_ARGUMENTS = SEQUENCE .rsr,  RESULT = CSRFILE'

  STRING_OPERATE OPERATION = 'CONCATENATE', ;
                  STRING_ARGUMENTS = SEQUENCE '.sch',  RESULT = SCHFILE

  STRING_OPERATE OPERATION = 'CONCATENATE', ;
                 STRING_ARGUMENTS = SEQUENCE '.mat',  RESULT = MATRIX_FILE

  SET ROOT_NAME = SEQUENCE

  RETURN
END_SUBROUTINE

### __homcsr.top
# Purpose: to construct typical comparative modelling restraints
# Input  : ALNFILE, SEQUENCE, MODEL, CSRFILE (necessary)
#          KNOWNS, DATABASE_DIRECTORY (for homology restraints)
#          TOPLIB, PARLIB
#
# Output : stereochemical and homology restraints in file CSRFILE

SUBROUTINE ROUTINE = 'homcsr'
  # Check the alignment
  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS SEQUENCE

  CHECK_ALIGNMENT
  # make topology and build/read the atom coordinates:
  CALL ROUTINE = 'GENERATE_METHOD'
  # exit early?
  IF ARGUMENTS = EXIT_STAGE 2, OPERATION = 'EQ', THEN = 'RETURN'
  # make and write the stereochemical, homology, and special restraints?
  IF ARGUMENTS = CREATE_RESTRAINTS 0, OPERATION = 'EQ', THEN='GO_TO __SKIP_RSRS'
    CALL ROUTINE = 'mkhomcsr'
  LABEL __SKIP_RSRS

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'mkhomcsr'

  MAKE_RESTRAINTS RESTRAINT_TYPE = 'stereo', ADD_RESTRAINTS = OFF
  # go_to 'next'

  # SET ALIGNMENT_FORMAT = 'PIR'
  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS SEQUENCE

  MAKE_RESTRAINTS RESTRAINT_TYPE = 'phi-psi_binormal', ADD_RESTRAINTS = ON

  SET SPLINE_RANGE = 4.0, SPLINE_DX = 0.3, SPLINE_MIN_POINTS = 5

  MAKE_RESTRAINTS RESTRAINT_TYPE = 'omega_dihedral', ADD_RESTRAINTS = ON
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'chi1_dihedral', ADD_RESTRAINTS = ON
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'chi2_dihedral', ADD_RESTRAINTS = ON
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'chi3_dihedral', ADD_RESTRAINTS = ON
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'chi4_dihedral', ADD_RESTRAINTS = ON

  SET SPLINE_RANGE = 4.0, SPLINE_DX = 0.7, SPLINE_MIN_POINTS = 5

  # Only do the standard residue types for CA, N, O, MNCH, SDCH dst rsrs
  # (no HET or BLK residue types):
  SET RES_TYPES = 'STD'

  SET DISTANCE_RSR_MODEL = 5, MAXIMAL_DISTANCE = MAX_CA-CA_DISTANCE
  SET RESIDUE_SPAN_RANGE = 2 99999, RESIDUE_SPAN_SIGN = on
  SET RESTRAINT_GROUP = 9
  PICK_ATOMS PICK_ATOMS_SET = 2, ATOM_TYPES = 'CA'
  PICK_ATOMS PICK_ATOMS_SET = 3, ATOM_TYPES = 'CA'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'distance', ADD_RESTRAINTS = 'ON'
  
  SET DISTANCE_RSR_MODEL = 6, MAXIMAL_DISTANCE = MAX_N-O_DISTANCE
  SET RESIDUE_SPAN_RANGE = 2 99999, RESIDUE_SPAN_SIGN = off
  SET RESTRAINT_GROUP = 10
  PICK_ATOMS PICK_ATOMS_SET = 2, ATOM_TYPES = 'N'
  PICK_ATOMS PICK_ATOMS_SET = 3, ATOM_TYPES = 'O'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'distance', ADD_RESTRAINTS = 'ON'

  SET DISTANCE_RSR_MODEL = 6, MAXIMAL_DISTANCE = MAX_SC-MC_DISTANCE
  SET RESIDUE_SPAN_RANGE = 1 2, RESIDUE_SPAN_SIGN = off
  SET RESTRAINT_GROUP = 23, RESTRAINT_STDEV = 0.5 1.5
  PICK_ATOMS PICK_ATOMS_SET = 2, ATOM_TYPES = 'SDCH'
  PICK_ATOMS PICK_ATOMS_SET = 3, ATOM_TYPES = 'MNCH'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'distance', ADD_RESTRAINTS = 'ON'

  SET DISTANCE_RSR_MODEL = 6, MAXIMAL_DISTANCE = MAX_SC-SC_DISTANCE
  SET RESIDUE_SPAN_RANGE = 2 99999, RESIDUE_SPAN_SIGN = on
  SET RESTRAINT_GROUP = 26, RESTRAINT_STDEV = 0.5 2.0
  PICK_ATOMS PICK_ATOMS_SET = 2, ATOM_TYPES = 'SDCH'
  PICK_ATOMS PICK_ATOMS_SET = 3, ATOM_TYPES = 'SDCH'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'distance', ADD_RESTRAINTS = 'ON'

  # Generate intra-HETATM and HETATM-protein restraints:
  CALL ROUTINE = 'hetatm_restraints'
  # Generate intra-BLK and BLK-protein restraints:
  CALL ROUTINE = 'blk_restraints'
  # Special restraints have to be called last so that possible cis-proline
  # changes are reflected in the current restraints:
  CALL ROUTINE = 'special_restraints'

  # label 'next'
  CONDENSE_RESTRAINTS

  WRITE_RESTRAINTS FILE = CSRFILE

  SET RESIDUE_SPAN_RANGE = -999 -999, RESIDUE_SPAN_SIGN = on

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'hetatm_restraints'

  # Select MODEL=7 where you can specify stand. dev. explicitly
  SET RESTRAINT_TYPE = 'distance'
  SET DISTANCE_RSR_MODEL = 7
  SET MAXIMAL_DISTANCE = 7.0
  SET ADD_RESTRAINTS = ON
  SET RESTRAINT_GROUP = 27

  # Inter- and intra- residue:
  SET RESTRAINT_STDEV = 0.2 0.0
  # Note: there are going to be duplicated HETATM-HETATM restraints:
  SET RESIDUE_SPAN_RANGE = 0 99999, RESIDUE_SPAN_SIGN = off
  PICK_ATOMS PICK_ATOMS_SET = 2, ATOM_TYPES = 'ALL', RES_TYPES = 'ALL'
  PICK_ATOMS PICK_ATOMS_SET = 3, ATOM_TYPES = 'ALL', RES_TYPES = 'HET'
  MAKE_RESTRAINTS

  RETURN
END_SUBROUTINE


# Re-define the auxillary restraints routine to include the restraints between
# the protein and the BLK residues. The BLK atoms will be restrained by their 
# distances to the protein CA atoms that are within MAXIMAL_DISTANCE angstroms 
# of the selected BLK atoms in the templates. Note: this only works because 
# the BLK atoms have unique atom names.
SUBROUTINE ROUTINE = 'blk_restraints'

  # To derive restraints from all (one) templates, comment out (uncomment):
#  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = '3b5c' SEQUENCE
#  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS SEQUENCE

  # Select MODEL=7 where you can specify stand. dev. explicitly
  SET RESTRAINT_TYPE = 'distance'
  SET DISTANCE_RSR_MODEL = 7
  SET MAXIMAL_DISTANCE = 10.0
  SET ADD_RESTRAINTS = ON
  SET RESTRAINT_GROUP = 27

  # Intra-residue:
  SET RESTRAINT_STDEV = 0.05 0.0
  SET RESIDUE_SPAN_RANGE = 0 0, RESIDUE_SPAN_SIGN = on
  PICK_ATOMS PICK_ATOMS_SET = 2, ATOM_TYPES = 'ALL',  RES_TYPES = 'BLK'
  PICK_ATOMS PICK_ATOMS_SET = 3, ATOM_TYPES = 'ALL',  RES_TYPES = 'BLK'
  MAKE_RESTRAINTS

  # Inter-residue:
  SET RESTRAINT_STDEV = 0.2 0.0
  SET RESIDUE_SPAN_RANGE = 1 99999, RESIDUE_SPAN_SIGN = off
  # There may be some duplicated CA BLK - CA BLK restraints:
  PICK_ATOMS PICK_ATOMS_SET = 2, ATOM_TYPES = 'CA',  RES_TYPES = 'ALL'
  PICK_ATOMS PICK_ATOMS_SET = 3, ATOM_TYPES = 'ALL', RES_TYPES = 'BLK'
  MAKE_RESTRAINTS

  # Clean after yourself:
  DELETE_ALIGNMENT

  RETURN
END_SUBROUTINE

### __multiple_models.top
# Purpose : calculate several models from the initial structure and restraints
# Input   : MODEL, CSRFILE, SCHFILE, ROOT_NAME, DEVIATION, MD_LEVEL,
#           RAND_METHOD
#
# Output  : schedule in file SCHFILE, final models in files ROOT_NAME.B*,
#           debugging output in files ROOT_NAME.D*
#

SUBROUTINE ROUTINE = 'multiple_models'

   DO ID2 = STARTING_MODEL, ENDING_MODEL, 1
     SET  FINAL_MODEL = 'default'
     CALL ROUTINE     = 'single_model'
   END_DO

   IF ARGUMENTS = DO_LOOPS 1, OPERATION = 'NE', THEN = 'GO_TO NO_LOOPS'

   # Do thorough loop modeling if so selected, for each *.B9999???? model:
   DO ID2 = STARTING_MODEL, ENDING_MODEL, 1
     SET ROOT_NAME = SEQUENCE, FILE_ID = '.B', FILE_EXT = PDB_EXT, ID1 = 9999
     CALL ROUTINE = 'loop', LOOP_MODEL = '${DEFAULT}'
   END_DO

   LABEL NO_LOOPS

RETURN
END_SUBROUTINE

### __single_model.top
# Purpose : calculate a single model from the initial structure and 
#           restraints. Possibly randomize it before optimization
#           (RAND_METHOD and DEVIATION <> 0) and possibly refine it
#           after variable target function CG optimization (MD_LEVEL).
#
# Input   : MODEL, CSRFILE, SCHFILE, FINAL_MODEL, DEVIATION,
#           RAND_METHOD, MD_LEVEL, REPEAT_OPTIMIZATION
#
# Output  : schedule in file SCHFILE, final model in file FINAL_MODEL,
#           debugging output in .D, energy profile in .E, violations in .V
#
SUBROUTINE ROUTINE = 'single_model'
   SET MAX_ITERATIONS = MAX_VAR_ITERATIONS

   SET ID1 = 0
   SWITCH_TRACE FILE = 'default', FILE_EXT = '', FILE_ID = '.D'

   # Vary the initial structure
   # Note that you are counting on some MODEL arrays not being deleted by 
   # READ_MODEL (ie the charge(1:natm) array generated by GENERATE_TOPOLOGY)
   READ_MODEL FILE = MODEL
   CALL ROUTINE = 'select_atoms'
   CALL ROUTINE = RAND_METHOD

   # Generate the var targ func schedule, write it to a file
   MAKE_SCHEDULE
   WRITE_SCHEDULE FILE = SCHFILE

##   IF ARGUMENTS = WRITE_INTERMEDIATES 1, OPERATION = 'EQ', ;
##      THEN = 'WRITE_MODEL FILE = default, FILE_EXT = PDB_EXT', FILE_ID = '.B'

   # repeat the basic optimization cycle REPEAT_OPTIMIZATION times:
   DO  IREPEAT = 1, REPEAT_OPTIMIZATION, 1
     CALL ROUTINE = 'single_model_pass'
   END_DO

   IUPAC_MODEL

   SET ID1 = 9999
   SET OUTPUT2 = OUTPUT
#   ENERGY OUTPUT = 'LONG'
   ENERGY OUTPUT = 'LONG ENERGY_PROFILE', FILE = 'default', FILE_ID = '.E', FILE_EXT = ''
   ENERGY OUTPUT = 'LONG VIOLATIONS_PROFILE', FILE = 'default', FILE_ID = '.V', FILE_EXT = ''
   SET OUTPUT = OUTPUT2

   # write the final model:
   WRITE_MODEL FILE = FINAL_MODEL, FILE_ID = '.B', FILE_EXT = PDB_EXT
   SET FINAL_MODEL = FILE

RETURN
END_SUBROUTINE


# The basic optimization cycle:
SUBROUTINE ROUTINE = 'single_model_pass'
 # do the variable target function optimization:

 READ_MODEL FILE = MODEL, MODEL_SEGMENT = '539:A' '672:D'
## # Define the two segments (chains in this case) to be identical:
 CALL ROUTINE = 'defsym', SEG1 = '539:A' '672:A', SEG2 = '539:B' '672:B'
 CALL ROUTINE = 'defsym', SEG1 = '539:A' '672:A', SEG2 = '539:C' '672:C'
 CALL ROUTINE = 'defsym', SEG1 = '539:A' '672:A', SEG2 = '539:D' '672:D'

   DO SCHEDULE_STEP = 1, N_SCHEDULE, 1
       OPERATE ARGUMENTS = ID1 1, OPERATION = 'SUM', RESULT = ID1
       # CALL ROUTINE = 'rd_restraints' # Not needed each time anymore
       PICK_RESTRAINTS ADD_RESTRAINTS = off

       OPTIMIZE 
###
## PROVA:
#       OPTIMIZE OPTIMIZATION_METHOD = 3,TEMPERATURE = 320,MD_RETURN = 'Minimal',INIT_VELOCITIES = on, MAX_ITERATIONS = 70
###
       IF ARGUMENTS = WRITE_INTERMEDIATES 1, OPERATION = 'EQ', ;
          THEN = 'WRITE_MODEL FILE = default, FILE_EXT = PDB_EXT, FILE_ID = .B'
       IF ARGUMENTS = ERROR_STATUS 0, OPERATION = 'GT', ;
         THEN = 'GO_TO __ABORT'
       IF ARGUMENTS = MOLPDF MAX_MOLPDF, OPERATION = 'GT', ;
         THEN = 'GO_TO __ABORT'
   END_DO

   CALL ROUTINE = 'refine', SCHEDULE_STEP = N_SCHEDULE

   LABEL __ABORT
   # Recover the status
   SET ERROR_STATUS = 0

RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'randomize_dihedrals'

  # Just in case, generate topology again (ROTATE needs bonds)
  # (could replace with GENERATE_TOPOLOGY if no special patches)
  CALL ROUTINE = 'create_topology'
  # Optimize all dihedral angles:
  SET DIHEDRALS = 'phi psi omega chi1 chi2 chi3 chi4'
###
## PROVA:
  # Uncomment the following block for optimization instead of randomization:
#  CALL ROUTINE = 'rd_restraints'
#  PICK_RESTRAINTS ADD_RESTRAINTS = off, SCHEDULE_STEP = N_SCHEDULE
#  CONDENSE_RESTRAINTS
#  ROTATE_DIHEDRALS CHANGE = 'optimize'
###
  ROTATE_DIHEDRALS CHANGE = 'randomize', DEVIATION = 3.0 # uses DEVIATION
###
 RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'randomize_xyz'

  RANDOMIZE_XYZ DEVIATION = 3.0, RAND_SEED = -11237      # uses DEVIATION

RETURN
END_SUBROUTINE



SUBROUTINE ROUTINE = 'rd_restraints'

   READ_RESTRAINTS FILE = CSRFILE, ADD_RESTRAINTS = 'off'

RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'local_optimization'

   # locally optimize sidechains
   ROTATE_DIHEDRALS CHANGE = 'optimize', DIHEDRALS = 'chi1 chi2 chi3 chi4' 
   # remove bad clashes one final time, with all restraints
#  OPTIMIZE
   OPTIMIZE OPTIMIZATION_METHOD = 1, MAX_ITERATIONS = 100

RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'nothing'
RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'select_atoms'
  # Select all atoms in the system:
   PICK_ATOMS SELECTION_SEGMENT='@:@' 'X:X', SELECTION_SEARCH='segment',;
             PICK_ATOMS_SET=1, RES_TYPES='all', ATOM_TYPES='all', ;
             SELECTION_FROM='all', SELECTION_STATUS='initialize'
RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'defsym'

   SET ATOM_TYPES = 'ALL'
   SET SELECTION_STATUS = 'INITIALIZE'
   SET SELECTION_SEARCH = 'SEGMENT'
##
   SET SYMMETRY_WEIGHT = 1.0
   PICK_ATOMS PICK_ATOMS_SET = 2, SELECTION_SEGMENT = SEG1
   PICK_ATOMS PICK_ATOMS_SET = 3, SELECTION_SEGMENT = SEG2
   DEFINE_SYMMETRY ADD_SYMMETRY = on off
##
  RETURN
END_SUBROUTINE

### __refine.top
# Purpose: refine hot spots in a model using Molecular Dynamics and 
#          Conjugate Gradients. If RSTRS_REFINED = 0, hot spots are
#          those atoms that violate stereochemical restraints (excluding
#          mainchain and sidechain dihedral angles phi,psi,omega,chi_i, but
#          including hard sphere overlap) than 4 standard deviations. These are
#          then refined. If 1, atoms that violate important homology-derived
#          restraints for more than 4 standard deviations are also hot spots.
#          If 2, any sufficiently violated restraints are hotspots.
#
# Input  : CSRFILE, RSTRS_REFINED
#
# Assuming that a topology library and the model are already in memory.

SUBROUTINE ROUTINE = 'refine'
   STRING_IF STRING_ARGUMENTS = MD_LEVEL 'none', OPERATION = 'EQ', THEN = ;
     RETURN

   # Get precise radii:
#    SET TOPOLOGY_MODEL = 3
#    READ_TOPOLOGY FILE = '${LIB}/top_heav.lib'
   # Get model with hot spots:
   READ_MODEL FILE = MODEL
   # obtain bond connectivities for local optimization (sequence from MODEL)
   # (this is only needed when RSTRS_REFINED >= 1)
   DELETE_ALIGNMENT
   CALL ROUTINE = 'create_topology'

   # Save the current model:
   STRING_IF STRING_ARGUMENTS = FIT_IN_REFINE 'NO_FIT', OPERATION = 'NE', ;
     THEN = 'WRITE_MODEL FILE = TO_BE_REFINED.TMP'

   # Reset schedule to all restraints with factors of 1:
#   MAKE_SCHEDULE
#   SET SCHEDULE_STEP = N_SCHEDULE

   # Get all static restraints:
   # (already there by default)
   CALL ROUTINE = 'rd_restraints' 
   CALL ROUTINE = 'select_atoms'
   PICK_RESTRAINTS ADD_RESTRAINTS = off

   SET MAX_ITERATIONS_STORE = MAX_ITERATIONS

   # Possibly skip selecting hot atoms only and optimize all atoms:
   IF ARGUMENTS = REFINE_HOT_ONLY 1, OPERATION = 'NE', THEN ;
   'GO_TO __BEGIN_MD'

   SET VIOL_REP_STORE  = VIOL_REPORT_CUT

   # Refine only hotspots that have badly violated stereochemical restraints:
   IF ARGUMENTS = RSTRS_REFINED 0, OPERATION = 'EQ', THEN ;
   'SET VIOL_REPORT_CUT = 4 4 4 4 4 999 999 999 999 999 999 999 999 999 999 999 999 4 4 4 999 999 999 999 999 999 999 999 999 999 999'

   # Refine hotspots that have badly violated stereochemical
   # restraints and the important homology-derived restraints:
   IF ARGUMENTS = RSTRS_REFINED 1, OPERATION = 'EQ', THEN ;
   'SET VIOL_REPORT_CUT = 4 4 4 4 4 999 999 999 4 4 999 999 4 999 999 999 999 4 4 4 999 999 4 999 4 4 999 999 999 999'

   # Refine hotspots that have badly violated any kind of restraints
   IF ARGUMENTS = RSTRS_REFINED 2, OPERATION = 'EQ', THEN ;
   'SET VIOL_REPORT_CUT = 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4'
   
   # Pick hot atoms (must use RESIDUE mode because of sidechains):
   SET PICK_HOT_CUTOFF = 4.5, SELECTION_MODE  = 'RESIDUE', EXTEND_HOT_SPOT = 0
   PICK_HOT_ATOMS 

   SET VIOL_REPORT_CUT = VIOL_REP_STORE

   # Pick all corresponding (violated and others) restraints:
   PICK_RESTRAINTS ADD_RESTRAINTS = off

   # Local optimization to prevent MD explosions:
   OPTIMIZE OPTIMIZATION_METHOD = 1, MAX_ITERATIONS = 150

   LABEL __BEGIN_MD

   # Do simulated annealing MD:
   CALL ROUTINE = MD_LEVEL

   # Possibly skip 'HOT CG' after MD:
   IF ARGUMENTS = REFINE_HOT_ONLY 1, OPERATION = 'NE', THEN ;
   'GO_TO __AFTER_MD'

   # Get conjugate gradients refined hot spots:
   OPTIMIZE OPTIMIZATION_METHOD = 1, MAX_ITERATIONS = 200

   # Get all static restraints again and select all atoms
   CALL ROUTINE = 'rd_restraints'
   CALL ROUTINE = 'select_atoms'
   PICK_RESTRAINTS ADD_RESTRAINTS = off

   LABEL __AFTER_MD

   # Get a final conjugate gradients refined structure:
   OPTIMIZE OPTIMIZATION_METHOD = 1, MAX_ITERATIONS = 200

   # Evaluate gross changes between the initial and final refined model:
   STRING_IF STRING_ARGUMENTS = FIT_IN_REFINE 'NO_FIT', OPERATION = 'INDEX', ;
     THEN = 'GO_TO __NO_FIT'
   DELETE_ALIGNMENT
   READ_MODEL2 FILE = 'TO_BE_REFINED.TMP'
   PICK_ATOMS ATOM_TYPES = 'CA', PICK_ATOMS_SET = 1
   SUPERPOSE
   PICK_ATOMS ATOM_TYPES = 'ALL'
   SUPERPOSE
   CALL ROUTINE = 'select_atoms'
   # Delete TO_BE_REFINED.TMP
   SET OPERATION = 'CONCATENATE'
   STRING_OPERATE STRING_ARGUMENTS = ' ./' OUTPUT_DIRECTORY, RESULT = COMMAND
   STRING_OPERATE STRING_ARGUMENTS = 'rm' COMMAND,  RESULT = COMMAND
   STRING_OPERATE STRING_ARGUMENTS = COMMAND '/TO_BE_REFINED.TMP'
   SYSTEM

   LABEL __NO_FIT

   # Restore some TOP variables which are likely to be used afterwards
   SET MAX_ITERATIONS  = MAX_ITERATIONS_STORE

RETURN
END_SUBROUTINE

## nuovo:
SUBROUTINE ROUTINE = 'refine_1'
   # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
   SET UPDATE_DYNAMIC = 0.39, CONTACT_SHELL = 4.0
   SET CAP_ATOM_SHIFT = UPDATE_DYNAMIC, MD_TIME_STEP = 4.0

   # Equilibration:
   SET MAX_ITERATIONS = 50, EQUILIBRATE = 10, OPTIMIZATION_METHOD = 3

   OPTIMIZE TEMPERATURE = 150.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = on
   OPTIMIZE TEMPERATURE = 400.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = off
   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'

   # Sampling:
   SET EQUILIBRATE = 100, MAX_ITERATIONS = 300

   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 800.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 500.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 300.0, MD_RETURN = 'FINAL'

   RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'refine_2'
   # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
   SET UPDATE_DYNAMIC = 0.39, CONTACT_SHELL = 4.0
   SET CAP_ATOM_SHIFT = UPDATE_DYNAMIC, MD_TIME_STEP = 4.0

   # Equilibration:
   SET MAX_ITERATIONS = 100, EQUILIBRATE = 20, OPTIMIZATION_METHOD = 3

   OPTIMIZE TEMPERATURE = 150.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = on
   OPTIMIZE TEMPERATURE = 250.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = off
   OPTIMIZE TEMPERATURE = 500.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'

   # Sampling:
   SET EQUILIBRATE = 100, MAX_ITERATIONS = 400

   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 800.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 500.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 300.0, MD_RETURN = 'FINAL'

   RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'refine_3'
   # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
   SET UPDATE_DYNAMIC = 0.39, CONTACT_SHELL = 4.0
   SET CAP_ATOM_SHIFT = UPDATE_DYNAMIC, MD_TIME_STEP = 4.0

   # Equilibration:
   SET MAX_ITERATIONS = 200, EQUILIBRATE = 20, OPTIMIZATION_METHOD = 3

   OPTIMIZE TEMPERATURE = 150.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = on
   OPTIMIZE TEMPERATURE = 250.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = off
   OPTIMIZE TEMPERATURE = 400.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 700.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'

   # Sampling:
   SET EQUILIBRATE = 200, MAX_ITERATIONS = 600

   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 800.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 600.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 500.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 400.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 300.0, MD_RETURN = 'FINAL'

   RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'refine_4'
   # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
   SET UPDATE_DYNAMIC = 0.39, CONTACT_SHELL = 4.0
   SET CAP_ATOM_SHIFT = UPDATE_DYNAMIC, MD_TIME_STEP = 4.0

   # Equilibration:
   SET MAX_ITERATIONS = 300, EQUILIBRATE = 20, OPTIMIZATION_METHOD = 3

   OPTIMIZE TEMPERATURE = 150.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = on
   OPTIMIZE TEMPERATURE = 250.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = off
   OPTIMIZE TEMPERATURE = 400.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 700.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE =1300.0, MD_RETURN = 'FINAL'

   # Sampling:
   SET EQUILIBRATE = 200, MAX_ITERATIONS = 1000

   OPTIMIZE TEMPERATURE =1300.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 800.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 600.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 500.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 430.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 370.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 320.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 300.0, MD_RETURN = 'FINAL'

   RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'refine_5'
   # at T=1000, max_atom_shift for 4fs is cca 0.15 A.
   SET UPDATE_DYNAMIC = 0.39, CONTACT_SHELL = 4.0
   SET CAP_ATOM_SHIFT = UPDATE_DYNAMIC, MD_TIME_STEP = 10.0

   # Equilibration:
   SET MAX_ITERATIONS = 200, EQUILIBRATE = 20, OPTIMIZATION_METHOD = 3

   OPTIMIZE TEMPERATURE = 150.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = on
   OPTIMIZE TEMPERATURE = 250.0, MD_RETURN = 'FINAL', INIT_VELOCITIES = off
   OPTIMIZE TEMPERATURE = 400.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 700.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE =1500.0, MD_RETURN = 'FINAL'

   # Sampling:
   SET EQUILIBRATE = 200, MAX_ITERATIONS = 2000

   OPTIMIZE TEMPERATURE =1500.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE =1000.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 800.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 600.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 500.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 400.0, MD_RETURN = 'FINAL'
   OPTIMIZE TEMPERATURE = 300.0, MD_RETURN = 'FINAL'

   RETURN
END_SUBROUTINE


### __generate_model.top
# Purpose: build a structure from a sequence in the alignment file:
#          either build from scratch (GENERATE_METHOD='generate_xyz') or
#          copy the equivalent coordinates from the homologue(s)
#          (GENERATE_METHOD='transfer_xyz').
#
# Input  : ALNFILE, SEQUENCE, MODEL
#
# System : TOPLIB, PARLIB
#
# Output : file with a model, and a file with stereochemical restraints

SUBROUTINE ROUTINE = 'generate_xyz'

  # SET ALIGNMENT_FORMAT = 'PIR'
  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = SEQUENCE

  READ_TOPOLOGY   FILE = TOPLIB
  READ_PARAMETERS FILE = PARLIB

  CALL ROUTINE = 'create_topology'

  BUILD_MODEL INITIALIZE_XYZ = ON, BUILD_METHOD = '3D_INTERPOLATION'
  WRITE_MODEL FILE = MODEL

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'transfer_xyz'

  # SET ALIGNMENT_FORMAT = 'PIR'
  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS
# READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS

  IF ARGUMENTS = INITIAL_MALIGN3D 0, OPERATION = 'EQ', THEN = 'GO_TO NO_INITIAL_MALIGN3D'

  MALIGN3D FIT = off, GAP_PENALTIES_3D = 0 4

  LABEL NO_INITIAL_MALIGN3D

  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS SEQUENCE

  READ_TOPOLOGY   FILE = TOPLIB
  READ_PARAMETERS FILE = PARLIB

  CALL ROUTINE = 'create_topology'

  TRANSFER_XYZ 
  BUILD_MODEL INITIALIZE_XYZ = OFF
  RENAME_SEGMENTS SEGMENT_IDS = 'A' 'B' 'C' 'D', RENUMBER_RESIDUES = 539 539 539 539

  WRITE_MODEL FILE = MODEL

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'read_xyz'

  # SET ALIGNMENT_FORMAT = 'PIR'
  # READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = SEQUENCE

  READ_TOPOLOGY   FILE = TOPLIB
  READ_PARAMETERS FILE = PARLIB

  DELETE_ALIGNMENT
  READ_MODEL FILE = MODEL
  CALL ROUTINE = 'create_topology'
  
  # READ_MODEL FILE = MODEL

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'create_topology'

  GENERATE_TOPOLOGY ADD_SEQUENCE = OFF 
  CALL ROUTINE = 'default_patches'
  CALL ROUTINE = 'special_patches'

  RETURN
END_SUBROUTINE

### __default_patches.top
# Purpose: generate default patches
# Input  : MODEL (after GENERATE_MODEL)
# Output : patched MODEL

# Derive the presence of disulfides from template structures
# (you can still define the additional disulfide in the
#  special_patches TOP routine)

SUBROUTINE ROUTINE = 'default_patches'
  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS SEQUENCE
  PATCH_SS_TEMPLATES
  RETURN
END_SUBROUTINE

### __loop.top
# This file defines the current loop modeling routines. An example of
# how to use CALL ROUTINE = 'loop' is given in the comments following
# this. 
#
# Input: none, edit this file
# Output: generated models in *.M* files, calculated energies in *.E* files
#
SET SEQUENCE = 'HERG'             # the sequence code
SET LOOP_MODEL = 'HERG.B9999????.modlr.pdb'       # the structure into which the loop is built
# SET LOOP_MODEL = 'HERG.B????????.modlr.pdb'
# 
SET LOOP_STARTING_MODEL = 1
SET LOOP_ENDING_MODEL = 25  # how many loop models are calculated
SET RAND_SEED = -34871     # to get different starting models each time
SET ID2 = 0                
# 
CALL ROUTINE = 'loop'                      # do the loop models
# Try to get an even better model by clustering (optional):
CALL ROUTINE = 'cluster', ID1 = LOOP_STARTING_MODEL, ID2 = LOOP_ENDING_MODEL
#
# This routine picks model residues that need to be refined (necessary):
#
SUBROUTINE ROUTINE = 'select_loop_atoms'
#   # Make sure your atom selection is NOT a function of coordinates!!!!
#   # (because atom selections for restraints construction and optimization
#   #  will be different!) This means SELECTION_SEGMENT = 'SPHERE' is bad!
# 
  PICK_ATOMS SELECTION_SEGMENT =   '569:A'  '584:A', SELECTION_STATUS = 'initialize'
  PICK_ATOMS SELECTION_SEGMENT =   '569:B'  '584:B', SELECTION_STATUS = 'add'
  PICK_ATOMS SELECTION_SEGMENT =   '569:C'  '584:C', SELECTION_STATUS = 'add'
  PICK_ATOMS SELECTION_SEGMENT =   '569:D'  '584:D', SELECTION_STATUS = 'add'

  PICK_ATOMS SELECTION_SEGMENT =   '594:A'  '603:A', SELECTION_STATUS = 'add'
  PICK_ATOMS SELECTION_SEGMENT =   '594:B'  '603:B', SELECTION_STATUS = 'add'
  PICK_ATOMS SELECTION_SEGMENT =   '594:C'  '603:C', SELECTION_STATUS = 'add'
  PICK_ATOMS SELECTION_SEGMENT =   '594:D'  '603:D', SELECTION_STATUS = 'add'

  RETURN
END_SUBROUTINE

# This routine adds any special restraints (optional):
#
SUBROUTINE ROUTINE = 'special_restraints'
  SET ADD_RESTRAINTS = on
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'ALPHA', RESIDUE_IDS = '585:A' '593:A'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'ALPHA', RESIDUE_IDS = '604:A' '611:A'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'ALPHA', RESIDUE_IDS = '585:B' '593:B'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'ALPHA', RESIDUE_IDS = '604:B' '611:B'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'ALPHA', RESIDUE_IDS = '585:C' '593:C'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'ALPHA', RESIDUE_IDS = '604:C' '611:C'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'ALPHA', RESIDUE_IDS = '585:D' '593:D'
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'ALPHA', RESIDUE_IDS = '604:D' '611:D'

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'loop'
  CALL ROUTINE = 'getnames-loop'

  READ_TOPOLOGY FILE = TOPLIB
  READ_TOPOLOGY FILE = '$(LIB)/top_heav.lib'

  READ_ATOM_CLASSES ATOM_CLASSES_FILE = '$(LIB)/atmcls-melo.lib'

  READ_PARAMETERS FILE = PARLIB
  READ_PARAMETERS FILE = '$(LIB)/par.lib', ADD_PARAMETERS = off
  READ_PARAMETERS FILE = '$(LIB)/melo1-dist.lib', ADD_PARAMETERS = on

  # Best is 7 A, but have to comment it out to be able to set it in the
  # the main .top file:
  SET  CONTACT_SHELL   = 7.00
  
  SET OUTPUT = 'LONG'

  # Get the initial coordinates from the input atom file, make sure
  # that the disulfides occuring in the input atom file are patched
  # in the model and that the residue numbering is also inherited.
  # Also, it is not necessary that the input atom file has either
  # all atoms or atoms in the standard RTF order.

  READ_MODEL FILE = LOOP_MODEL           # MODEL_SEGMENT !
  SEQUENCE_TO_ALI ALIGN_CODES = SEQUENCE, ATOM_FILES = LOOP_MODEL, ADD_SEQUENCE = off
  SEQUENCE_TO_ALI ADD_SEQUENCE = on, ;
     ALIGN_CODES = ALIGN_CODES SEQUENCE, ATOM_FILES = ATOM_FILES LOOP_MODEL
  GENERATE_TOPOLOGY
  # Create the disulfide bonds using equivalent disulfide bonds in templates
  # (in case there are disulfides bonds in loops):
  PATCH_SS_TEMPLATES
  CALL ROUTINE = 'special_patches'
  TRANSFER_XYZ
  BUILD_MODEL INITIALIZE_XYZ = off
  READ_MODEL2 FILE = LOOP_MODEL
  TRANSFER_RES_NUMB 

  CALL ROUTINE = 'select_loop_atoms'

  # Select all atoms to be modeled 'ab initio':
  CALL ROUTINE = 'loop_restraints'

  # Select corresponding restraints only:
  # only necessary to eliminate inefficiencies in 'special_restraints' 
  # because MAKE_RSRS works with selected atoms now:
##  PICK_RESTRAINTS ADD_RESTRAINTS = off 
  CONDENSE_RESTRAINTS
##  WRITE_RESTRAINTS FILE = LOOP_CSRFILE

  # Calculate energy for the original (raw) loop:
  ENERGY NONBONDED_SEL_ATOMS = 1

  # Prepare the starting structure (comment it out if
  # the input PDB file is a better initial structure):
  CALL ROUTINE = 'build_ini_loop'
  SET ROOT_NAME = SEQUENCE, FILE_ID = '.IL', FILE_EXT = PDB_EXT, ID1 = 0
  SET LOOP_INI_MODEL = '${DEFAULT}'

  WRITE_MODEL FILE = LOOP_INI_MODEL

  MAKE_SCHEDULE LIBRARY_SCHEDULE = 6
  SET WRITE_ALL_ATOMS = on
  DO  ID1 = LOOP_STARTING_MODEL, LOOP_ENDING_MODEL, 1
    READ_MODEL FILE = LOOP_INI_MODEL
    CALL ROUTINE = 'select_loop_atoms'
    RANDOMIZE_XYZ DEVIATION = 5.0
    SET ROOT_NAME = SEQUENCE, FILE_ID = '.DL'
    SWITCH_TRACE FILE = '${DEFAULT}', TRACE_OUTPUT = 15
  # Refine without the rest of the protein:
    CALL ROUTINE = 'optimize', NONBONDED_SEL_ATOMS = 2
    # Refine in the context of the rest of the protein:
    CALL ROUTINE = 'optimize', NONBONDED_SEL_ATOMS = 1
    ENERGY
    SET ROOT_NAME = SEQUENCE, FILE_ID = '.BL'
    IUPAC_MODEL
    WRITE_MODEL FILE = 'default'
  END_DO

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'loop_restraints'

# do not use template derived rsrs at all:
  SET DIH_LIB_ONLY = on 

  SET DYNAMIC_SPHERE = on
  SET DYNAMIC_MODELLER = on

  SET RESTRAINT_SEL_ATOMS = 1 # only restraints that have that many atoms in 
                              # the selected set are going to be calculated

  MAKE_RESTRAINTS RESTRAINT_TYPE = 'bond', ADD_RESTRAINTS = off
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'angle', ADD_RESTRAINTS = on
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'improper', ADD_RESTRAINTS = on
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'dihedral', ADD_RESTRAINTS = on
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'phi-psi_binormal', ADD_RESTRAINTS = ON

  SET SPLINE_RANGE = 4.0, SPLINE_DX = 0.3, SPLINE_MIN_POINTS = 5

  MAKE_RESTRAINTS RESTRAINT_TYPE = 'omega_dihedral', ADD_RESTRAINTS = ON

  MAKE_RESTRAINTS RESTRAINT_TYPE = 'chi1_dihedral', ADD_RESTRAINTS = ON
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'chi2_dihedral', ADD_RESTRAINTS = ON
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'chi3_dihedral', ADD_RESTRAINTS = ON
  MAKE_RESTRAINTS RESTRAINT_TYPE = 'chi4_dihedral', ADD_RESTRAINTS = ON

  CALL ROUTINE = 'special_restraints'

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'optimize'

  SET OPTIMIZATION_METHOD = 1, MAX_ITERATIONS = 200, MIN_ATOM_SHIFT = 0.001
  OPTIMIZE SCHEDULE_STEP = 1
  OPTIMIZE SCHEDULE_STEP = 2
  OPTIMIZE SCHEDULE_STEP = 3
  OPTIMIZE SCHEDULE_STEP = 4
  OPTIMIZE SCHEDULE_STEP = 5

  CALL ROUTINE = LOOP_MD_LEVEL

  SET OPTIMIZATION_METHOD = 1, MAX_ITERATIONS = 1000, MIN_ATOM_SHIFT = 0.00001
  OPTIMIZE

  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'getnames-loop'
  STRING_OPERATE OPERATION = 'CONCATENATE', ;
                 STRING_ARGUMENTS = SEQUENCE '.lrsr',  RESULT = LOOP_CSRFILE
  RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'build_ini_loop'

  UNBUILD_MODEL
  BUILD_MODEL BUILD_METHOD = '3D_INTERPOLATION', INITIALIZE_XYZ = off

  RETURN
END_SUBROUTINE


# This routine picks model residues that are modeled in __loops 
# (necessary, but can be re-defined just as other routines).
# This version of the atom selection routine automatically picks all
# model residues that occur in 'loops':
#
SUBROUTINE ROUTINE = 'select_loop_atoms'
  READ_ALIGNMENT FILE = ALNFILE, ALIGN_CODES = KNOWNS SEQUENCE, ADD_SEQUENCE = off
  SET SELECTION_STATUS = 'INITIALIZE'
  SET SELECTION_MODE = 'ATOM'
  SET SELECTION_SEARCH = 'SEGMENT'
  SET SELECTION_FROM = 'ALL'
# Uncomment, to allow user selection:
#   SET GAP_EXTENSION = 3 3 # gaps are extended for two residues in both
#                           # directions to get loops (insertions & deletions)
  SET MAX_LOOP_LENGTH = 25
  PICK_ATOMS ATOM_TYPES = 'ALL', SELECTION_SEGMENT = 'LOOPS' '', ;
             SELECTION_SEARCH = 'SEGMENT', GAP_EXTENSION = 3 3
# PICK_ATOMS ATOM_TYPES = 'ALL', SELECTION_SEGMENT = 'LOOPS' '', ;
#            SELECTION_SEARCH = 'SPHERE_SEGMENT', SPHERE_RADIUS = 6.0, ;
#            GAP_EXTENSION = 3 3
# SUPERPOSE
  RETURN
END_SUBROUTINE

### __principal.top
# Purpose: Do principal components on family FAMILY and produce an ASGL plot
#
# Input  : TOP variable FAMILY
#
# Output files: FAMILY.mat, FAMILY.dat.
#

SUBROUTINE ROUTINE = 'principal_components'

STRING_OPERATE OPERATION = CONCATENATE, STRING_ARGUMENTS = FAMILY '.ali', ;
               RESULT = FILE
READ_ALIGNMENT

STRING_OPERATE OPERATION = CONCATENATE, STRING_ARGUMENTS = FAMILY '.mat', ;
               RESULT = MATRIX_FILE
ID_TABLE

STRING_OPERATE OPERATION = CONCATENATE, STRING_ARGUMENTS = FAMILY '.dat', ;
               RESULT = FILE
PRINCIPAL_COMPONENTS

DENDROGRAM

RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'special_patches'
  RETURN
END_SUBROUTINE

### __complete.top
# To put any missing (heavy) atoms on the model and to reorder the atoms 
# according to the MODELLER/CHARMM order:
#
# Input : MODEL; filename of the input atom file
# Output: complete model in memory 

SUBROUTINE ROUTINE = 'complete'
# Use this to put on heavy atoms only:
# READ_TOPOLOGY FILE = '${LIB}/top_heav.lib'
# Use this to put on all atoms, including hydrogens:
  READ_TOPOLOGY FILE = '${LIB}/top.lib'
  READ_PARAMETERS FILE = '${LIB}/par.lib'

# read the sequence from the input atom file:
  READ_MODEL FILE = MODEL

# generate the internal coordinates
  GENERATE_TOPOLOGY
  CALL ROUTINE = 'special_patches'

# copy the known atom coordinates from MODEL to model in the memory
  DELETE_ALIGNMENT
  SEQUENCE_TO_ALI ATOM_FILES = MODEL, ALIGN_CODES = MODEL
  SEQUENCE_TO_ALI ADD_SEQUENCE = on, ATOM_FILES = ATOM_FILES MODEL, ;
                  ALIGN_CODES = ALIGN_CODES MODEL
  TRANSFER_XYZ

# Fill in the coordinates for the missing atoms:
  BUILD_MODEL INITIALIZE_XYZ = 'OFF'

 RETURN
END_SUBROUTINE

STOP