PBC_Code.py

#Code written by Hanh Nguyen. hanhn3@uci.edu

import MDAnalysis
import numpy as np
from MDAnalysis import analysis
import matplotlib.pyplot as plt
import math 

traj = ['prod1.trr', 'prod2.trr', 'prod3.trr', 'prod4.trr', 'prod5.trr', 'prod6.trr', 'prod7.trr', 'prod8.trr', 'prod9.trr',
        'prod10.trr', 'prod25.trr', 'prod50.trr', 'prod75.trr', 'prod100.trr']

gro = ['mol1.gro', 'mol2.gro', 'mol3.gro', 'mol4.gro', 'mol5.gro', 'mol6.gro', 'mol7.gro', 'mol8.gro', 'mol9.gro', 
       'mol10.gro', 'mol25.gro', 'mol50.gro', 'mol75.gro', 'mol100.gro'] 

DMSO_Percentage = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 25, 50, 75, 100]

#Set periodic boundary conditions
MDAnalysis.core.flags['use_periodic_selections'] = True 
MDAnalysis.core.flags['use_KDTree_routines'] = False 

#Arrays to store number of Oxygen and Sulfur atoms 
Ave_S_dmso = []
Ave_O_wat = []
Ave_S_dmso_bulk = []
Ave_O_wat_bulk = []
expected_S_dmso = []
expected_O_wat = []

#Clear previous figures and data 
plt.cla()
plt.clf()
plt.gca().set_color_cycle(['red', 'green', 'blue',  'black', 'purple', 'pink', 'orange', 'maroon'])

# markers = ['o', 'v', '*','x', '+']

# Radius of the sphere 
rad = 8.0
for percentage,traj_file,gro_file  in zip(DMSO_percentage, traj, gro):
    print "Working on %s percentage of dmso "%(percentage)
    frames = MDAnalysis.Universe(gro_file,traj_file)    
    # the groups of atoms of O1 of water to search for 
    wat_group = frames.select_atoms('resname wat and name O1')
    # the groups of atoms of dmso to search for 
    dmso_group = frames.select_atoms('bynum 18:27 or resname DMS and name O1')
    # Select the solute 
    sol = frames.select_atoms('bynum 1:17')
    
    
    #Reset the counter for O and S atoms and box volume for each percentage of DMSO
    O_total = 0 
    S_total = 0
    O_total_bulk = 0 
    S_total_bulk = 0 
    box_volume = 0 
    Nexpected_wat = 0 
    Nexpected_DMSO = 0
    ##  NOW WORKING ON NUMBER OF S AND O ATOMS NEAR SOLUTE AND IN BULK##    
    # Loop through all the frames 
    for ts in frames.trajectory:
        #Find the center of mass of solute
        center_sol = sol.center_of_mass()
        xcoord = center_sol[0]
        ycoord = center_sol[1]
        zcoord = center_sol[2]
        ##           WATER PART          ##
        #Select all oxygen atoms of water within a sphere around the center of mass of that solute 
        sele_O_wat = frames.select_atoms('(point   %f   %f   %f   %f) and resname wat and name O1'%(xcoord,ycoord,zcoord,rad))
        #Sum up the number of selected atoms 
        O_total = len(sele_O_wat) + O_total
        ##           DMSO PART          ##
        #Select all oxygen atoms of water within a sphere around the center of mass of that solute 
        sele_S_dms = frames.select_atoms('(point   %f   %f   %f   %f) and name S1'%(xcoord,ycoord,zcoord,rad))
        #Sum up the number of selected atoms 
        S_total = len(sele_S_dms) + S_total 
        ##  NOW WORKING IN BULK  ##
        #Find oxygen atoms within ~16A radius of the center of mass of solute 
        O_22A_away_sol = frames.select_atoms('sphlayer  15.0  16.0 (bynum 1:17)') 
        #Take the first point from the different set and select that atom with that index 
        point = frames.select_atoms('bynum %s'%(O_22A_away_sol.indices[0]))
        pointx = point.positions[0][0]
        pointy = point.positions[0][1]
        pointz = point.positions[0][2]
        # Select all oxygen atoms of water around the center of mass of that solute  in bulk 
        sele_O_wat_bulk = frames.select_atoms('(point   %f   %f   %f    %f) and resname wat and name O1'%(pointx,pointy,pointz,rad))
        # Sum up the number of selected atoms  in bulk 
        O_total_bulk = len(sele_O_wat_bulk) + O_total_bulk 

        # Select all oxygen atoms of water around the center of mass of that solute in bulk 
        sele_S_dms_bulk = frames.select_atoms('(point   %f   %f   %f    %f) and name S1'%(pointx,pointy,pointz,rad))
        
        # Sum up the number of selected atoms in bulk 
        S_total_bulk = len(sele_S_dms_bulk) + S_total_bulk
        
        #Find the sum of volume of the box from each frame
        box_volume = box_volume + ts.volume 
        
    #Divide the total of Oxygens of water and Sulfurs of DMSO around the solute and in bulk over total frames (1001 frames)

    Ave_S_dmso_bulk.append(S_total_bulk / len(frames.trajectory) ) 
    
    Ave_O_wat_bulk.append(O_total_bulk / len(frames.trajectory) ) 
      
    Ave_S_dmso.append(S_total / len(frames.trajectory) ) 
    
    Ave_O_wat.append(O_total / len(frames.trajectory) )
##  NOW WORKING ON EXPECTED NUMBER OF S AND O ATOMS ##    
    # Find average volume of the box over all frames: 
    average_box_volume = box_volume / len(frames.trajectory) 
    # solvent group 
    solv_group = len(dmso_group) + len(wat_group) 
    # calculate volume of sample: 
    vol_sample = (4/3) * math.pi * rad**3 
    # Number of expected DMSO 
    Nexpected_DMSO = len(dmso_group) * ((vol_sample ) /average_box_volume) 
    # Append this value to an array of expected number of s atoms of dmso 
    expected_S_dmso.append(Nexpected_DMSO)
    # Number of expected water
    Nexpected_wat = len(wat_group) * ( (vol_sample ) /average_box_volume) 
    # Append this value to an array of expected number of o atoms of water
    expected_O_wat.append(Nexpected_wat)

##  NOW WORKING ON GRAPHING ##  
    
#Plot the numbers of oxygen atoms near solute: and in bulk scripted values versus expected values 
plt.plot(DMSO_Percentage,expected_O_wat)
plt.plot(DMSO_Percentage,Ave_O_wat)
plt.plot(DMSO_Percentage,Ave_O_wat_bulk)
#Plot the numbers of S atoms near solute and in bulk: scripted values versus expected values 
plt.scatter(DMSO_Percentage,expected_S_dmso, color = 'purple' , marker = 'o' , linewidths='3')
plt.scatter(DMSO_Percentage,Ave_S_dmso, color = 'orange', marker = 'v', linewidths='4')
plt.scatter(DMSO_Percentage,Ave_S_dmso_bulk, color ='maroon', marker = 'D', linewidths='5')
plt.xlabel('Percentage of DMSO')
plt.ylabel('Average number of O and S atoms near solute and in bulk over 1001 frames')
plt.title('Average number of O-wat and S-dmso atoms around the solute and in bulk over total frames')
plt.legend(['Mathematically_expected_number_of_O_of_wat', 'MDAnalysis_calculated_number_of_O_of_wat', 'MDAnalysis_calculated_number_of_O_of_wat_16Aaway_solute', 'Mathematically_expected_number_of_S_of_dmso','MDAnalysis_calculated_number_of_S_of_dmso', 'MDAnalysis_calculated_number_of_S_of_dmso_16Aaway_solute'], loc='upper right')
plt.xlim(0,100)
plt.savefig('methylparabe_compute_neighbors_PBC')
plt.show()
print "Job Done."