# Example for: model.optimize(), misc.switch_trace()

# This will optimize stereochemistry of a given model, including
# non-bonded contacts.

env = environ()
env.edat.dynamic_sphere = True

env.libs.topology.read(file='$(LIB)/top_heav.lib')
env.libs.parameters.read(file='$(LIB)/par.lib')

code = '1fas'
aln = alignment(env)
mdl = model(env, file=code)
aln.append_model(mdl, align_codes=code, atom_files=code)
aln.append_model(mdl, align_codes=code+'-ini', atom_files=code+'.ini')
mdl.generate_topology(aln, sequence=code+'-ini')
mdl.transfer_xyz(aln)
mdl.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')
mdl.write(file=code+'.ini')

# Generate the restraints:
mdl.restraints.make(aln, restraint_type='stereo', spline_on_site=False)
mdl.restraints.write(file=code+'.rsr')

mpdf = mdl.energy()
mdl.switch_trace(file=code+'.trc')
mdl.optimize(optimization_method=1, max_iterations=20, output='REPORT',
             trace_output=1)
mdl.optimize(optimization_method=3, temperature=300, max_iterations=50,
             output='REPORT', trace_output=1)
mdl.optimize(optimization_method=1, max_iterations=20, output='REPORT',
             trace_output=1)
mpdf = mdl.energy()

mdl.write(file=code+'.B')