import numpy as np
import numpy
import scipy
import os


def rmsd(crds1,crds2):
	crds1,crds2=np.array(crds1), np.array(crds2)
	"""Returns RMSD between 2 sets of [nx3] numpy array"""
	n_vec = numpy.shape(crds1)[0]
	correlation_matrix = numpy.dot(numpy.transpose(crds1), crds2)
	v, s, w = numpy.linalg.svd(correlation_matrix)
	is_reflection = (numpy.linalg.det(v) * numpy.linalg.det(w)) < 0.0
	if is_reflection:
		s[-1] = - s[-1]
	E0 = sum(sum(crds1 * crds1)) + \
		sum(sum(crds2 * crds2))
	#print E0
	rmsd_sq = (E0 - 2.0*sum(s)) / float(n_vec)
	rmsd_sq = max([rmsd_sq, 0.0])
	return numpy.sqrt(rmsd_sq)

  

'''
a=[0,0,0]
b=[5,5,5]
c=[1]

c1=np.array([a,a,a,b,a,a])
c2=np.array([b,b,b,a,a,b])
'''



#print rmsd(np.array(c1),np.array(c2))


coor1=[]
w=open('receptor.pdb').readlines()

for line in w:
	if 'ATOM' in line[:4]:
		if 'C' == line.split()[-1]:
			x=line.split()[-6]
			y=line.split()[-5]
			z=line.split()[-4]
			c=[float(x),float(y),float(z)]
			coor1.append(c)
		
print len(coor1)		



coor2=[]
w=open('complex.pdb').readlines()

for line in w:
	if 'ATOM' in line[:4]:
		if 'C' == line.split()[-1]:
			x=line.split()[-6]
			y=line.split()[-5]
			z=line.split()[-4]
			c=[float(x),float(y),float(z)]
			coor2.append(c)
		
print len(coor2)		

c1=coor1
c2=coor2

print rmsd(np.array(c1),np.array(c2))
