poging1.1.py

import sys
import numpy as np
from itertools import permutations as perm

rmsd_list=[1000]
final_atom_names=[]

def rmsd(V, W):
    """
    Calculate Root-mean-square deviation from two sets of vectors V and W.

    Parameters
    ----------
    V : array
        (N,D) matrix, where N is points and D is dimension.
    W : array
        (N,D) matrix, where N is points and D is dimension.

    Returns
    -------
    rmsd : float
        Root-mean-square deviation

    """
    D = len(V[0])
    N = len(V)
    rmsd = 0.0
    for v, w in zip(V, W):
        rmsd += sum([(v[i] - w[i])**2.0 for i in range(D)])
    return np.sqrt(rmsd/N)

def calculator(initial_file, final_file):

    p_atoms, p_all = "xyz", initial_file
    q_atoms, q_all = "xyz", final_file 
    #de ingelezen data

    # Calculate 'dumb' RMSD
    normal_rmsd =rmsd(p_all[:], q_all[:])
    if normal_rmsd < rmsd_list[0]:
        rmsd_list[0] = normal_rmsd 
        #de calculator
        return normal_rmsd

#--------------------------- inlezen --> xyz file 1 is altijd de initial file, xyz file 2 is altijd de final file

def inread(file_xyz, list_coordinates, list_names, initial):
    filename = file_xyz
    datafile = open(filename)

    lines= list(datafile.readlines())

    if initial == True:
        molecule_length.append(int(lines.pop(0).rstrip()))
        #hier isoleren we de molecuul lengte van het initial-molecule.
    elif initial == False:
        lines.pop(0) #haalt het eerste element uit de lijst, dit is de lengte van het molecuul
    else:
        print("Fourth statement should be True or False")

    for line in lines: 
    #lines bestaat uit een lijst met elk atoom en zijn coordinaten als apart element in deze lijst
        list_coordinates.append(map(str, line.split())) 
        #maakt van elk atoom en zijn coordinaten een aparte lijst en zet deze in een overkoepelende lijst (initial_molecule). 
    list_coordinates.pop(0) 
    #haalt het eerste element uit de lijst omdat dit een n/ was dus dit resulteert nu in een lege lijst

    for atoms in list_coordinates:
         list_names.append(atoms.pop(0).rstrip()) 
        #verwijdert uit alle lijsten in initial_molecule het eerste element (de atoomnaam) en plaatst deze in de lijst atom_names. 

    for atoms in list_coordinates: 
    #maakt van alle overgebleven coordinaten in initial_molecule een float
        for i, coordinates in enumerate(atoms):
            atoms[i]= float(coordinates)





#--------------------------- permuteren en vergelijken
#dit wordt de lijst waar elke keer het gepermuteerde molecuul in opgeslagen gaat worden
def permutate_and_calculate(molecule_length, initial_molecule_coordinates ,final_molecule_coordinates, final_atom_names):
    

    permutated_final_molecule=[] 
    permutated_atom_names=[]    
    counter = 0
    
    for i in perm(range(molecule_length)):
    #deze loop gaat elke keer 1 finial molecule permuteren en dan vergelijken met het initial molecule
    #i: lijst met permutatie nummers (bv: 1,2,3,4,5)
        for j in i:
        #j: is dit cijfer individueel
            permutated_final_molecule.append(final_molecule_coordinates[j])
            #haalt uit de final_molecule lijst de coordinaten met de indexen die hij krijgt via j
            permutated_atom_names.append(final_atom_names[j])   
            #haalt uit de atom_names_final lijst de atoom namen met de zelfde indexen die hij krijgt via j
        permutated_final_molecule= np.asarray(permutated_final_molecule) 
        #maakt van de final_molecule_lijst een array (is nodig als input)
        if calculator(initial_molecule_coordinates, permutated_final_molecule) != None: 
        #als de nieuwe berekende rmsd groter is dan de vorige returned het progamma een none, anders returned hij de nieuwe rmsd
            better_list[0] = permutated_atom_names
            better_list[1] = permutated_final_molecule
        permutated_final_molecule=[] 
        #leegt de lijsten zodat ze opnieuw gebruikt kunnen worden
        permutated_atom_names=[]
        counter += 1
    print ("Hij heeft {} loops gerunt".format(counter))


molecule_length=[]
initial_molecule = []
atom_names = []

inread(sys.argv[1], initial_molecule, atom_names, True)

final_molecule= []
atom_names_final = []

inread(sys.argv[2], final_molecule, atom_names_final, False)

initial_molecule= np.asarray(initial_molecule) 
#maakt van de lijst met lijsten met coordinaten een array.

better_list=[0,0]

permutate_and_calculate(molecule_length[0],initial_molecule, final_molecule, atom_names_final)
print (better_list[0])
print (better_list[1])