TIMD extraction.py

From Rizzo_Lab
Jump to: navigation, search
import math, sys
import os.path

#to use the script, choose the right fn number and the right filename prefix
#################################################################################################################
## loops over "fn" of output files .
## fn = file number 
def read_TIMD_out(fileprefix,fn):

    ### get average and rms values of DV/DL, TEMP, PRESS, Density parameters

    average_temp=[];average_pres=[];average_dvdl=[];average_dens=[];rms_temp=[];rms_pres=[];rms_dvdl=[];rms_dens=[];
    for i in range(1, fn+1):
        file = fileprefix + str(i) +".out";
        file1 = open(file, 'r');

    ##flag for Average and RMS values
        flag_average = 0;
        flag_rms     = 0;

        lines  =  file1.readlines();
        num_lines = len(lines);

        for line in lines:
            linesplit = line.split() #split on white space
            if (len(linesplit) > 2):
                if (linesplit[0] == "A" and linesplit[1] == "V"):
                    flag_average = 1;
                    flag_rms     = 0;

                if (linesplit[0] == "R" and linesplit[1] == "M"):
                    flag_average = 0;
                    flag_rms     = 1;

                if (flag_average == 1 and linesplit[0] == "DV/DL"):
                    average_dvdl.append( float(linesplit[2]) );

                if (flag_average == 1 and linesplit[0] == "Density"):
                    average_dens.append( float(linesplit[2]) );

                if (flag_rms == 1 and linesplit[0] == "DV/DL"):
                    rms_dvdl.append( float(linesplit[2]) );

                if (flag_rms == 1 and linesplit[0] == "Density"):
                    rms_dens.append( float(linesplit[2]) );
                    flag_rms = 0; ###reading is finished, remove the possibility of other readings

            if (len(linesplit) > 6):
                if (flag_average == 1 and linesplit[6] == "TEMP(K)"):
                    average_temp.append( float(linesplit[8]) );
                    average_pres.append( float(linesplit[11]));

                if (flag_rms == 1 and linesplit[6] == "TEMP(K)"):
                    rms_temp.append( float(linesplit[8]) );
                    rms_pres.append( float(linesplit[11]));

        file1.close();
    return average_temp,average_pres,average_dvdl,average_dens,rms_temp,rms_pres,rms_dvdl,rms_dens;
#################################################################################################################
def write_TIMD_report(ksteps,average_temp,average_pres,average_dvdl,average_dens,rms_temp,rms_pres,rms_dvdl,rms_dens,file):

    fn = len(average_temp);
    file.write("Lambda  DVDL   rms   ksteps TEMP  rms    PRESS  rms   Density rms"+"\n");
    for i in range(1,fn+1):
        file.write( "%5.3f%9.3f%6.3f%4d%9.2f%5.2f%6.1f%9.1f%7.4f%8.4f" %(i/10.0,average_dvdl[i-1],rms_dvdl[i-1],ksteps,average_temp[i-1],rms_temp[i-1],average_pres[i-1],rms_pres[i-1],average_dens[i-1],rms_dens[i-1]) );
        file.write( "\n" );
    return;
#################################################################################################################
def main():
    
    fileprefix="bnz_prod_v0_l";
    num=9;##number of files to be read
    ksteps=100;

    average_temp,average_pres,average_dvdl,average_dens,rms_temp,rms_pres,rms_dvdl,rms_dens = read_TIMD_out( fileprefix, num );

    file_report = "report." + fileprefix;
    file = open( file_report, 'write' );
    write_TIMD_report(ksteps,average_temp,average_pres,average_dvdl,average_dens,rms_temp,rms_pres,rms_dvdl,rms_dens,file); 
    file.close()
#################################################################################################################
#################################################################################################################
main()