# Example for: model.energy()

# This will calculate the stereochemical energy (bonds,
# angles, dihedrals, impropers) for a given model.

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

aln = alignment(env)
code = "1fas"
mdl = model(env, file=code)
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')
# Must patch disulfides here to calculate the non-bonded
# energy properly. Also, when you use hydrogens, disulfides
# must always be patched so that sulfhydril hydrogens are
# removed from the model.
for ids in [ ('17', '39'),
             ( '3', '22'),
             ('53', '59'),
             ('41', '52') ]:
    mdl.patch(residue_type='DISU', residue_ids=ids)

mdl.transfer_xyz(aln)
mdl.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')

mdl.restraints.make(aln, restraint_type='stereo', spline_on_site=False)
molpdf = mdl.energy(edat=energy_data(dynamic_sphere=True))