# Example for: model.patch(), topology(), parameters.read()

# This will define a CYS-CYS disulfide bond between residues 3 and 22.

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

# Read the sequence: 
code = '1fas'
mdl = model(env, file=code)
# have two copies of the sequence in the alignment, for TRANSFER_XYZ later:
aln = alignment(env)
aln.append_model(mdl, atom_files=code, align_codes=code)
aln.append_model(mdl, atom_files=code+'.ini', align_codes=code+'-ini')
mdl.generate_topology(aln, sequence=code+'-ini')

# Create the disulfide bond:
mdl.patch(residue_type='DISU', residue_ids=('3', '22'))

# Get MODEL's coordinates from the template, using the alignment (1:1 here):
mdl.transfer_xyz(aln)
# Calculate missing coordinates using internal coordinates:
mdl.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')

# Create the stereochemical restraints
mdl.restraints.make(aln, restraint_type='stereo', spline_on_site=False)

# Calculate the energy to test the disulfide:
mdl.energy()