# Example for: model.generate_topology(), model.build()

# This will read a specified atom file, generate all hydrogen atoms,
# add atomic radii and charges, and write the model to a PDB file in
# the GRASP format. This can be used with GRASP to display electrostatic
# properties without assigning charges and radii in GRASP.

log.verbose()
env = environ()

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

code = '1fas'
mdl = model(env, file=code)

aln = alignment(env)
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)

# Have to patch the topology here to remove sulfhydril hydrogens:
mdl.patch(residue_type='DISU', residue_ids=('17', '39'))
mdl.patch(residue_type='DISU', residue_ids=( '3', '22'))
mdl.patch(residue_type='DISU', residue_ids=('53', '59'))
mdl.patch(residue_type='DISU', residue_ids=('41', '52'))

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

mdl.write(file='1fas.ini1', model_format='GRASP')
mdl.write(file='1fas.ini2', model_format='PDB')