import sys
import math
import re

# Check usage
if (len(sys.argv) < 3):
  print "Error in usage!"
  sys.exit()

# read in PDB files
file1 = sys.argv[1];
file2 = sys.argv[2];
f = open(file1);
lines1 = f.readlines()
f.close()
f = open(file2);
lines2 = f.readlines()
f.close()


# Calculate and print RMSDs
i1 = -1
i2 = -1
pri = "" # previous residue index
prn = "" # previous residue name
sad = [] # array of squares of sidechain atom deviations of the current residue
while (1):
  f1 = 0
  for i1 in range(i1+1, len(lines1)):
    if not(re.match("^ATOM", lines1[i1])):
      continue
    rn1 = lines1[i1][17:20]  
    ri1 = lines1[i1][22:26]
    an1 = lines1[i1][12:16]
    if (re.match("(^\s*H\s*$|^\s*C\s*$|^\s*O\s*$|^\s*N\s*$|^\s*CA\s*$)", an1)):
      continue
    f1 = 1
    break

  f2 = 0
  for i2 in range(i2+1, len(lines2)):
    if not(re.match("^ATOM", lines2[i2])):
      continue
    rn2 = lines2[i2][17:20]  
    ri2 = lines2[i2][22:26]
    an2 = lines2[i2][12:16]
    if (re.match("(^\s*H|^\s*C\s*$|^\s*O\s*$|^\s*N\s*$|^\s*CA\s*$)", an2)):
      continue
    f2 = 1
    break

  if (f1 and f2 and not((rn1 == rn2) and (an1 == an2) and (ri1 == ri2))):
    print "Error - mismatched lines:\n%s%s" % (lines1[i1], lines2[i2])
    sys.exit()
  
  # have we started a new residue?
  if (not(ri1 == pri) or (not f1) or (not f2)):
    if (len(sad) > 0):
      rmsd = 0
      for i in range(len(sad)):
        rmsd += sad[i]
      rmsd = math.sqrt(rmsd/len(sad))
      print "%s %s %.2f" % (prn, pri, rmsd)
    # reset all arrays
    sad = []
    prn = rn1
    pri = ri1

  if (not(f1) or not(f2)):
    break

  # now process atom deviation
  x1 = float(lines1[i1][30:38])
  y1 = float(lines1[i1][38:46])
  z1 = float(lines1[i1][46:54])
  x2 = float(lines2[i2][30:38])
  y2 = float(lines2[i2][38:46])
  z2 = float(lines2[i2][46:54])
  sad.append((x1-x2)*(x1-x2) + (y1 - y2)*(y1 - y2) + (z1 - z2)*(z1 - z2))