In [ ]:

    

import numpy as np
import matplotlib.pyplot as plt


    

In [ ]:

    

# These are perfect crystal coordinates for the nodes

Q = np.array([
    [  9.985,   17.32050808,  16.58],
    [ -0.015,   34.6410162 ,  16.58],
    [ 29.970,   17.32050808,  16.58],
    [  9.985,   17.32050808,  33.16],
    [ -0.015,   34.6410162 ,  33.16],
    [ 29.970,   17.32050808,  33.16],
])


    

In [ ]:

    

def get_transmat(P, Q):
    
    
    def remove_com(P):
        old_com = np.mean(P, axis=0)
        newP = P - old_com
        return newP, old_com


    # get rotation vector from singular value decomposition (SVD)
    def kabsch(P0, Q0):
        A = np.dot(P0.transpose(), Q0)
        u,s,v = np.linalg.svd(A)  # a = U S V.H ... v==V.H, u==U
        
        U = u
        VT = v
        V = VT.transpose()
        UT = U.transpose()
        
        d = np.sign(np.linalg.det(np.dot(V, UT)))
        dmat = np.diag((1, 1, d))
        rotmat = np.dot(V, np.dot(dmat, UT))
        return rotmat

    
    # remove COM
    Q0, Qcom = remove_com(Q)
    P0, Pcom = remove_com(P)

    # get trans vector to go from P to Q
    trans = Qcom-Pcom
    rotmat = kabsch(P0, Q0)

    # convert to 4x4 matrix
    transmat = np.insert(rotmat, 3, trans, axis=1)
    transmat = np.insert(transmat, 3, [0,0,0,1], axis=0)
    
    return transmat


    

In [ ]:

    

def read_dumpfile(dumpfile):
    with open(dumpfile, 'r') as fin:
        data = {}
        for line in fin:
            if line[0]!='#':
                line_chunks=line.split()
                if len(line_chunks)==2:
                    timestep, N = map(int, line_chunks)
                    data[timestep] = []
                if len(line_chunks)==4:
                    serial, x,y,z = map(float, line_chunks)
                    ind = int(serial) - 1
                    data[timestep].append([x,y,z])
    return data


def rotate(xyz, transmat):
    xyz1 = np.append(xyz, np.ones((len(xyz),1)), axis=1).transpose()
    newxyz = np.dot(transmat, xyz1).transpose()
    return newxyz[:,:-1]


node_dumpfile = '/home/vargaslo/ws/_tmp_check_system/sarin/40/MDSim_298K/4_my_nodes.dump'
nodes = read_dumpfile(node_dumpfile)
mlc_dumpfile = '/home/vargaslo/ws/_tmp_check_system/sarin/40/MDSim_298K/4_my.dump'
mlc = read_dumpfile(mlc_dumpfile)


nodes_transmat = {}
for k,v in (nodes.iteritems()):
    P = np.array(v)
    transmat = get_transmat(P, Q)
    nodes_transmat[k] = transmat
    

dump = {}
for k,v in (mlc.iteritems()):
    xyz = np.array(v)
    transmat = nodes_transmat[k]
    xyz_ = rotate(xyz, transmat)
    dump[k] = xyz_


    

In [ ]:

    

def write_new_xyz(outfile):
    with open(outfile, 'w') as fout:
        for k in sorted(dump.iterkeys()):
            v = dump[k]
            N = len(v)
            fout.write('{} {}\n'.format(k, N))
            for ind, xyz in enumerate(v):
                serial = ind+1
                x,y,z = xyz
                fout.write('{} {:.6f} {:.6f} {:.6f}\n'.format(serial, x, y ,z))
    return

write_new_xyz(mlc_dumpfile+'.mod')