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()

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