In [ ]:

    

import sasmol.sasmol as sasmol
import numpy as np


    

In [ ]:

    

def create_pdb(coor, pdb_name=None, **kwargs):
    """
    coor is an n_atoms x 3 array
    create a pdb using coors
    """
    coor = np.asanyarray(coor)
    n_atoms = len(coor[0])
    index = np.arange(n_atoms) + 1

    create = sasmol.SasMol(0)
    create.setAtom(['ATOM'] * n_atoms)
    create.setIndex(index)
    create.setName(['Ar'] * n_atoms)
    create.setLoc([' '] * n_atoms)
    create.setResname(['VDW'] * n_atoms)
    create.setChain(['X'] * n_atoms)
    create.setResid(np.ones(n_atoms, dtype=np.int))
    create.setRescode([' '] * n_atoms)
    create.setCoor(coor.reshape(1, n_atoms, 3))
    create.setOccupancy(['1.00'] * n_atoms)
    create.setBeta(['0.00'] * n_atoms)
    create.setSegname(['{:04d}'.format(i) for i in index])
    create.setElement(['Ar'] * n_atoms)
    create.setCharge(['  '] * n_atoms)

    if not pdb_name:
        pdb_name = 'default_name.pdb'
    create.write_pdb(pdb_name, 0, 'w')
def position_molecule(coor, xyz_position):
    """
    this function is designed to generate a modified version of the input
    coordinates (coor3)
    1) center the lysozyme pdb
    2) store the centered lysozyme coordinates in a 4 vector, coor4 (1's in col_4)
    3) generate a random rotation matrix, rotate
    4) gererate translation matrix, translate (4x4 identity matrix, [3, :3] = [dx, dy, dz])
    5) take the dot product of the rt=rotate.dot(translate)
    6) transform the coordinates, coor4.dot(rt)
    """
    import numpy
    # 1)
    coor = np.asanyarray(coor)
    xyz_position = np.asanyarray(xyz_position)
    n_atoms = len(coor)
    # 2)
    # initialize vector arrays for coordinates and orientation vectors
    # changing them from 3 component vectors into 4 component vectors
    # this incorporates transaltions into the matrix math

    coor4 = numpy.ones((n_atoms, 4), dtype=numpy.float)
    coor4[:, 0:3] = coor
#     print(post_coor4)
#     post_coor4[:, 0:3] = post_coor3

    # 3)
    # create the translation-rotation matrix
    # This is intended to be multiplied from the right (unlike standard matrix
    # multiplication) so as not to require transposing the coordinate vectors.
    theta_xyz = numpy.random.random(3) * numpy.pi * 2 # radians

    [cx, cy, cz] = numpy.cos(theta_xyz)
    [sx, sy, sz] = numpy.sin(theta_xyz)

    # initialize the rotation
    # consolidated method of defining the rotation matrices
    rotate = numpy.eye(4, dtype=numpy.float)
    rotate[0][0:3] = [ cy*cz,          cy*sz,          -sy   ]
    rotate[1][0:3] = [ sx*sy*cz-cx*sz, sx*sy*sz+cx*cz, sx*cy ]
    rotate[2][0:3] = [ sx*sz+cx*sy*cz, cx*sy*sz-sx*cz, cx*cy ]

    # 4)
    translate = numpy.eye(4, dtype=numpy.float)
    translate[3, :3] = xyz_position

    # 5)
    rt = rotate.dot(translate) # check that the resulting matrix should rotate, then translate coordinates

    # 6)
    coor4 = coor4.dot(rt) # check that this does as expected
    return coor4[:, :3]


    

In [ ]:

    

#Load the lysozyme + center it
mol = sasmol.SasMol(0)
mol.read_pdb('/home/data/sascalc_pbc/pdbs/lysozyme_ca.pdb')
mol.center(0)
nLys = len(mol.coor()[0])

#Load the lennard Jones pdb
LJ_mol = sasmol.SasMol(0)
LJ_mol.read_pdb('/home/data/sascalc_pbc/pdbs/run_0.pdb')
LJ_mol.read_dcd('/home/data/sascalc_pbc/pdbs/run_1.dcd')
LJ_coors = LJ_mol.coor()[-2]

#Determine scaling factor so no overlap
from scipy.spatial.distance import cdist
fact = np.max(cdist(mol.coor()[0],mol.coor()[0]))/3.76 #3.76 Ar radius

#Create and populate the output coords
out_coors = np.zeros((1,len(LJ_coors)*len(mol.coor()[0]),3))
for i,c in enumerate(LJ_coors):
    out_coors[0][i*nLys: (i+1)*nLys] = position_molecule(mol.coor()[0],c*fact)
    
#create new pdb
create_pdb(out_coors,pdb_name='/home/data/sascalc_pbc/pdbs/argLys.pdb')
print("DONE")


    

In [ ]: