# 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'
     WRITE_MODEL FILE = default, FILE_EXT = PDB_EXT, FILE_ID = '.B'
   END_IF

   # 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
   STRING_IF STRING_ARGUMENTS = TOP_VERSION 'accelrys', OPERATION = 'EQ'
     # Accelrys code here:
     ENERGY OUTPUT = 'LONG ENERGY_PROFILE', NORMALIZE_PROFILE = OFF, ;
          FILE = 'default', FILE_ID = '.E', FILE_EXT = ''
     ENERGY OUTPUT = 'LONG ENERGY_PROFILE', NORMALIZE_PROFILE = ON, ;
          FILE = 'default', FILE_ID = '.NE', FILE_EXT = ''

     # The new request from Lisa/Azat to print out only stereochemical
     # restraint violations (6/24/03):
     # see doc_INT/energy-cmds.tex for information about other options
     # SET SMOOTHING_WINDOW = 1
     # SET NORMALIZE_PROFILE = off
     # select only stereochemical restraints (maybe add dihedral angles?):
     SET VIOL_REP_STORE  = SCHEDULE_SCALE
     SET SCHEDULE_SCALE = 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
     ENERGY OUTPUT = 'ENERGY_PROFILE', FILE = 'default', FILE_ID = '.ES', FILE_EXT = '', NORMALIZE_PROFILE = OFF
     ENERGY OUTPUT = 'ENERGY_PROFILE', NORMALIZE_PROFILE = ON, ;
            FILE = 'default', FILE_ID = '.NES', FILE_EXT = ''
     # ENERGY OUTPUT = 'VIOLATIONS_PROFILE'
     SET SCHEDULE_SCALE = VIOL_REP_STORE
   ELSE
     # Academic code here:
     ENERGY OUTPUT = 'LONG VIOLATIONS_PROFILE', FILE = 'default', FILE_ID = '.V', FILE_EXT = ''
   END_IF

   SET OUTPUT = OUTPUT2

   CALL ROUTINE = 'user_after_single_model'

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

   # Assess the model if desired
   ASSESS_MODEL

RETURN
END_SUBROUTINE


SUBROUTINE ROUTINE = 'user_after_single_model'
RETURN
END_SUBROUTINE



# The basic optimization cycle:

SUBROUTINE ROUTINE = 'single_model_pass'
   SET ABORT_OPT = 0
   # do the variable target function optimization:
   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
       IF ARGUMENTS = WRITE_INTERMEDIATES 1, OPERATION = 'EQ'
         WRITE_MODEL FILE = default, FILE_EXT = PDB_EXT, FILE_ID = .B
       END_IF
       IF ARGUMENTS = ERROR_STATUS 0, OPERATION = 'GT'
         SET ABORT_OPT = 1
         EXIT
       END_IF
       IF ARGUMENTS = MOLPDF MAX_MOLPDF, OPERATION = 'GT'
         SET ABORT_OPT = 1
         EXIT
       END_IF
   END_DO
   SET SCHEDULE_STEP = N_SCHEDULE

   IF ARGUMENTS = ABORT_OPT 0, OPERATION = 'EQ'
     CALL ROUTINE = 'refine', SCHEDULE_STEP = N_SCHEDULE
   END_IF

   # 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

  # 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 # uses DEVIATION

RETURN
END_SUBROUTINE




SUBROUTINE ROUTINE = 'randomize_xyz'

  # ROTATE_DIHEDRALS CHANGE = optimize, DIHEDRALS = chi1 chi2 chi3 chi4 
  RANDOMIZE_XYZ # 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

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