disorder.py

import numpy as np

models = 20
res = 166
filename = "lis"

file = open(filename,'r')

ave = [[0 for x in xrange(3)] for x in xrange(res)]
stdev = [[0 for x in xrange(3)] for x in xrange(res)]
coord = [[[0 for x in range(3)] for x in xrange(res)] for x in xrange(models)]
#print coord

for i in range(models):
    for j in range(res):
        line = file.readline()
        words = string.split(line)
        coord[i][j][0] = float(words[6])
        coord[i][j][1] = float(words[7])
        coord[i][j][2] = float(words[8])
        for k in range(3):
             ave[j][k] += coord[i][j][k]/models

for i in range(models):
    for j in range(res):
        for k in range(3):
            stdev[j][k] += (coord[i][j][k]-ave[j][k])**2/models

ave_stdev = []
for j in range(res):
     temp = 0.0
     for k in range(3):
         stdev[j][k] = math.sqrt(stdev[j][k])
         temp += stdev[j][k]
#    print j, ",", temp/3.0
     ave_stdev.append(temp/3.0)

#tails are removed
start_res = 5 #actually 6 but we count from 0
end_res = 158 #actually 159 but we count from 0

v = np.asarray(ave_stdev[start_res:end_res])

for j in range(start_res,end_res):
    if ave_stdev[j] > np.mean(v)+np.std(v):
       print j+1