Chimie Théorique » scripts_chimie4psmn

######################################################################################## 

#       Python script to calculate the mean density of a                                #

#               H2O along z.                                                                #

#                                                                                       #

#       Written by Rodrigo F. de Morais                                                #

#       Date: 08/04/2016                                                               #

#       Updated: 25/04/2016                                                             #

########################################################################################

#!/usr/bin/python

import sys

from numpy import sqrt

from netCDF4 import Dataset

###################### Functions #############################################

#----------------------------------------------------------------------------

#

########################### Main code #######################################

#---Definition of variables---

# ---Definition of code inicialiazation---

try:

   type=sys.argv[1] # Name of the PDB file to read

   type=sys.argv[2] # Traj file

except:

    print '----------------------------------------------------------------------'

    print 'The correct use of this script is:    '

    print 'python MeanRhoH2O_z.py prmtop mdcrd' 

    print '----------------------------------------------------------------------'

    sys.exit(2)

# ----------- Read Topology file

file=open(sys.argv[1],"r")

linhas=file.readlines()

file.close()

# ----------- Get total atoms

for i in range(0,len(linhas)):

        if len(linhas[i].split())>0 and linhas[i].split()[0]=='%FLAG':

                if linhas[i].split()[1]=='POINTERS':

                        n_atoms = linhas[i+2].split()[0] 

                        break

for i in range(0,len(linhas)):

        if len(linhas[i].split())>0 and linhas[i].split()[0]=='%FLAG':

                if linhas[i].split()[1]=='ATOM_NAME':

                        inicio=i+1

        if len(linhas[i].split())>0 and linhas[i].split()[0]=='%FLAG':

                if linhas[i].split()[1]=='CHARGE':

                        fim=i-1

                        break

AtomNames=[]

PosOAtoms=[]

PosPtAtoms=[]

# ------------ Get the position of O and Pt atoms

for i in range(inicio+1,fim):

        for j in range(0,20):

                AtomNames.append(linhas[i].split()[j])

for j in range(0,len(linhas[fim].split())):

                AtomNames.append(linhas[fim].split()[j])

for j in range(0,len(AtomNames)):

        if AtomNames[j]=='O':

                PosOAtoms.append(j)

        if AtomNames[j]=='Pt' or AtomNames[j]=='PU':

                PosPtAtoms.append(j)

#---------- Open the traj file in netcdf format

file=Dataset(sys.argv[2],"r",format='NETCDF4')

faixa=18.0

step=0.25

volume=24.031*19.424*step

goz=[0]*int(faixa/step)

zmax=0.0

for j in range(0,len(PosPtAtoms)):

                if zmax<float(file.variables['coordinates'][0][PosPtAtoms[j]][2]):

                       zmax=float(file.variables['coordinates'][0][PosPtAtoms[j]][2])

for i in range(0,len(file.variables['time'])/100): 

  CoordOtempz=[]

  for k in range(0,len(PosOAtoms)):

        CoordOtempz.append(float(file.variables['coordinates'][i*100][PosOAtoms[k]][2]))

  for j in range(0,int(faixa/step)):

        temp=0.0

        pos=0.0

        for k in range(0,len(CoordOtempz)):

                if CoordOtempz[k]>=zmax+step*(j-1) and CoordOtempz[k] <zmax+step*(j):

                                temp=temp+1

        if temp>0:

                goz[j]=goz[j]+temp/volume

FFile=open('goz',"w\n")

for i in range(0,len(goz)):

        FFile.write("%6.3f %6.3f %6.3f %6.3f\n" %(step*i,goz[i]/(len(file.variables['time'])/100)/0.033))

FFile.close()

file.close()

#----------------------------------------------------------------------------