I use Modeller to induce point mutations on a PDB. This part works very well (see code below). 

My problem is the protein has several disulfides defined by CONECT records in the template PDB. The output pdb is missing these records and they are required later. Since the atom count changes after several point mutations, it can not be done as easily, and doing it by hand is not an option.

Is there a way to force Modeller to do so? or it lacks the code to recompute the sulfurs?. I read  it does if the output file is a PSF, but not with PDB. 

Thanks in advance.

Pedro 

       mdl = model(env, file=pdb_file)
       aln = alignment(env)
       aln.append_model(mdl, atom_files=pdb_file, align_codes=pdb_file)

       for mutation in mutations:
           for pos, new_residue in mutation.items():
               chain_id, res_id = pos.split(":")
               pointmut = f'{res_id}:{chain_id}'
               sel = selection(mdl.residues[pointmut])
               sel.mutate(residue_type=self.res1to3[new_residue]) #mutate position
           
       aln.append_model(mdl, align_codes='mut')
       mdl.clear_topology()
       mdl.generate_topology(aln['mut'])
       mdl.transfer_xyz(aln)
       mdl.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')
       mdl2 = model(env, file=pdb_file)
       mdl.res_num_from(mdl2,aln)
       mdl.write(file=output_file, model_format='PDB')