import os
import random
import time
import math
import matplotlib.pyplot as plt
import numpy as np





def truncate(a):
    j=0
    b=""
    for i in range(len(a)):
        if a[i]==" ":
            pass
        else:
            b+=a[i]
            j+=1
    return b

def av(x,y):
    a=(x[0]+y[0])*0.5
    b=(x[1]+y[1])*0.5
    c=(x[2]+y[2])*0.5
    return [a,b,c]

def meanstdv(x):
	sum,mean,std,stdv=0,0,0,0
	for a in x:
		sum=sum+a
	mean=sum/len(x)
	for a in x:
		std=std+(a-mean)**2
	stdv=np.sqrt(std/(len(x)-1))
	return mean, stdv

def dis(x,y):
    return ((x[0]-y[0])**2+(x[1]-y[1])**2+(x[2]-y[2])**2)**0.5



def distance(coor1,coor2):
	t=[]
	for c1 in coor1:
		for c2 in coor2:
			t.append(dis(c1,c2))

	return t


def sorting(l1, l2):
    # l1 and l2 has to be numpy arrays
	newL1=np.array(l1)
	newL2=np.array(l2)
	idx = np.argsort(l1)
	return l1[idx], l2[idx]
	
	
	
	

directory='.'

extension = ".pdb"
files = [file for file in sorted(os.listdir(directory)) if file.lower().endswith(extension)]



yjuli=[]
for f in files:
	#print f
	res=[]
	w=open(f)
	lines=w.readlines()
	w.close()
	coor=[]
	startsentence='ATOM      1  N   PHE E  27'
	endsentence='ATOM  16311  HG  CYS E1059'
	for line in lines:
		if startsentence in line:
			startidx=lines.index(line)
		if endsentence in line:
			endidx=lines.index(line)+1
	for i in range(startidx,endidx):
		line=lines[i]
		newline=line[0:21] + line[22:]
		line=newline
		if line[:4]=='ATOM':
			if 'NH1 ARG  241' in line:
				x1=float(line.split()[5])
				y1=float(line.split()[6])
				z1=float(line.split()[7])
				c1=[x1,y1,z1]
				
			elif 'NH2 ARG  241' in line:
				x1=float(line.split()[5])
				y1=float(line.split()[6])
				z1=float(line.split()[7])
				c2=[x1,y1,z1]
	#print c1,c2
	
	
	
	
	gaba=[]
	for i in range(endidx,len(lines)):
		line=lines[i]
		if 'ATOM' in line[:4]:
			if line.split()[2] in ['O1','O2']:
				x1=float(line.split()[4])
				y1=float(line.split()[5])
				z1=float(line.split()[6])
				gaba.append([x1, y1, z1])
	#print len(gaba)
	#print distance([c1,c2],gaba)
	if min(distance([c1,c2],gaba))<6:
		print f,distance([c1,c2],gaba)
		
	
	
	