# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import numpy as np
import mdtraj as md
import sys
np.set_printoptions(threshold=sys.maxsize)


# +
t = md.load('OP19_7.pdb')
origin = 0 #index for first origin atom
atomspermol = 6 #how many atoms per mole
nmols = 3200 #total number of molecules
cutoff = 5 #cutoff in angstroms

atom_indices = np.zeros(nmols, dtype=int)
for i in range(len(atom_indices)):
    atom_indices[i] = int(i*atomspermol+origin)

t = t.atom_slice(atom_indices)
#get PBCs bounds.  This only works for orthohombic PBCs for the time being.  Might implement more general solution later.
pbc_vectors = np.diag(np.asarray((t.unitcell_vectors*10))[0])
coords = t.xyz*10 #defines coordinates.  multiply by 10 to get angstroms, since mdtraj converts to nm.

#computes distances from two sets of xyz coordinates
def distance(a,b,bounds):
    a = np.asarray((a))
    b = np.asarray((b))
    min_dists = np.min(np.dstack(((a - b) % bounds, (b - a) % bounds)), axis = 2)
    dist = np.sqrt(np.sum(min_dists ** 2, axis = 1))
    return dist



# +
neighbors = np.zeros((nmols,nmols))

for i in range(nmols):
    for j in range(i+1,nmols):
        if distance(coords[0][i],coords[0][j], pbc_vectors) <= cutoff:
            neighbors[i][j] = 1
            neighbors[j][i] = 1
nneighbors = np.sum(neighbors,axis=1)
# -



# +
f = open('OP19_7.cop','r')
g = open('normalOPs','w')

orderparameters = []
for x in f:
    orderparameters.append(x)
    
f.close()

bondops = np.array(orderparameters[(10+nmols):(10+2*nmols)],dtype=float)
normalops = bondops/(nneighbors+1)
np.savetxt('newOPs.txt',normalops)
# -



StarcoderdataJupyterScriptsDedupFiltered