getSplitWaters.py

# coding: utf-8

# In[1]:

import math
# gets the distance between two vec3
def get_distance(pos1, pos2):
    return (math.sqrt((pos1[0]-pos2[0])**2 +                       (pos1[1]-pos2[1])**2 +                       (pos1[2]-pos2[2])**2 ))

# checks if any OH distances are greater than the specified max_distance
def get_indices_of_split_waters(pdb_positions, max_distance):
    molecules_with_greater_distances = []
    for water_index in range(0,len(pdb_positions), 4):
        pos_O = pdb_positions[water_index]._value
        pos_H1 = pdb_positions[water_index+1]._value
        pos_H2 = pdb_positions[water_index+2]._value
        if ((get_distance(pos_O, pos_H1) > max_distance) or             (get_distance(pos_O, pos_H2) > max_distance)):
            molecules_with_greater_distances.append(water_index)
#            print(get_distance(pos_O, pos_H1))
#            print(get_distance(pos_O, pos_H2))
    
    print ("There are {} molecules with OH bonds that exceed the specified max distance".format(len(molecules_with_greater_distances)))
    print ("They are {}".format(molecules_with_greater_distances))
#    return molecules_with_greater_distances


# In[2]:

#run water256.py


# In[3]:

run water512.py


# In[4]:

get_indices_of_split_waters(pdb.positions, boxsize[2]/2)


# In[ ]: