import chimera
from chimera import runCommand

def findMAV():
	# locate and return the (presumably only) running instance
	# of MultAlign Viewer
	from MultAlignViewer.MAViewer import MAViewer
	from chimera.extension import manager
	mavs = [inst for inst in manager.instances
					if isinstance(inst, MAViewer)]
	if len(mavs) != 1:
		raise AssertionError("not exactly one MAV")
	return mavs[0]


# read a file named "structList" whose lines have 3 fields:  two PDB files
# and an output file name for the alignment
for line in open("structList", "r"):
	pdb1, pdb2, output = line.strip().split()

	# open the pdb files
	mol1 = chimera.openModels.open(pdb1, type="PDB")[0]
	mol2 = chimera.openModels.open(pdb2, type="PDB")[0]

	# use MatchMaker to get them superimposed
	from MatchMaker import match, CP_BEST, defaultMatrix, \
				defaultAlignAlgorithm, defaultGapOpen, \
				defaultGapExtend, defaultIterateCutoff
	# match the best-aligning chain of pdb1 onto pdb2, specifying
	# the similarity matrix, alignment algorithm, and the gap
	# open/extend penalties
	atoms1, atoms2, rmsd = match(CP_BEST, (mol1, [mol2]), defaultMatrix,
		defaultAlignAlgorithm, defaultGapOpen, defaultGapExtend,
		iterate=defaultIterateCutoff, showAlignment=True)[0]

	# find the Multalign Viewer that is now showing
	mav = findMAV()

	# find the original (ungapped) sequences corresponding to the
	# gapped sequences being displayed
	seqs = []
	for gapped in mav.seqs:
		for ungapped in gapped.molecule.sequences():
			if gapped.residues[0] in ungapped.residues:
				seqs.append(ungapped)
				break
		else:
			raise AssertionError("Could not find ungapped sequence"
					" for %s, %s" % (gapped.molecule.name,
					gapped.name))
	# dismiss this mav instance
	mav.Quit()

	# use Match->Align to create the superposition-based alignment
	import StructSeqAlign
	from StructSeqAlign.align import match2align
	# use 4.0 angstroms as the distance cutoff for matching, period (".")
	# as the gap character, and don't look for circular permutations
	match2align(seqs, 4.0, None, ".", False)

	# again find the new MAV instance
	mav = findMAV()

	# write out the alignment in MSF format
	# other format choices:  Clustal ALN, Aligned FASTA, Aligned NBRF/PIR,
	#	Pfam, GCG RSF, and Stockholm
	from MultAlignViewer.output import saverFunc
	out = open(output, "w")
	saverFunc["MSF"](out, mav, mav.seqs, mav.fileMarkups)
	out.close()

	# write out the RMSD
	import os.path
	base = os.path.splitext(output)[0]
	out = open(base + ".rmsd", "w")
	print>>out, "%.3f" % rmsd, "for", len(atoms1), "final residue pairs"
	# write out matched residues
	for r1, r2 in zip([a.residue for a in atoms1],
						[a.residue for a in atoms2]):
		print>>out, r1, r2
	out.close()

	# again quit the MAV
	mav.Quit()

	# write transformed PDB
	runCommand("write relative 0 1 " + base + ".pdb")

	# close the PDBs
	chimera.openModels.close([mol1, mol2])
raise chimera.ChimeraSystemExit("script finished\n")