Sent from Mail for Windows 10
From: Sehrish Naz
Sent: 19 March 2016 13:16
To: sehrish.naz.outlook.com
Subject: problem in plotting total energy graph by using perl script
Dear all,I had run 10 ns simulation of large compounds in TIP3P box to check their stable conformation with time. the input file I used for production is NPT production with no restrains &cntrl imin=0, ntx=7, irest=1, ntrx=1, ntxo=1, ntpr=500, ntwx=500, ntwv=200, ntwe=200, ntf=2, ntb=2, cut=10.0, nsnb=100, igb=0, nstlim=2000000, t=0.0, dt=0.001, ntt=3, gamma_ln=1.0, tempi=300.0, temp0=300.0, vlimit=20, ntp=1, taup=1.0, pres0=1.0, comp=44.6, ntc=2, tol=0.00000001, /after getting the production files. I tried to calculate the total energy graphs of the system by using perl script as follows . #!/usr/bin/perlif ($#ARGV 0) { print " Incorrect usage...\n"; exit;}foreach $i ( 0..$#ARGV ) { $filein = $ARGV[$i]; $checkfile = $filein; $checkfile =~ s/\.Z//; if ( $filein ne $checkfile ) { open(INPUT, "zcat $filein |") || die "Cannot open compressed $filein -- $!\n"; } else { open(INPUT, $filein) || die "Cannot open $filein -- $!\n"; } print "Processing sander output file ($filein)...\n"; &process_input; close(INPUT);}print "Starting output...\n";@sortedkeys = sort by_number keys(%TIME);@sortedavgkeys = sort by_number keys(%AVG_TIME);foreach $i ( TEMP, TSOLUTE, TSOLVENT, PRES, EKCMT, ETOT, EKTOT, EPTOT, DENSITY, VOLUME, ESCF ) { print "Outputing summary.$i\n"; open(OUTPUT, " summary.$i"); %outarray = eval "\%$i"; foreach $j ( @sortedkeys ) { print OUTPUT "$j ", $outarray{$j}, "\n"; } close (OUTPUT); print "Outputing summary_avg.$i\n"; open(OUTPUT, " summary_avg.$i"); %outarray = eval "\%AVG_$i"; foreach $j ( @sortedavgkeys ) { print OUTPUT "$j ", $outarray{$j}, "\n"; } close (OUTPUT); print "Outputing summary_rms.$i\n"; open(OUTPUT, " summary_rms.$i"); %outarray = eval "\%RMS_$i"; foreach $j ( @sortedavgkeys ) { print OUTPUT "$j ", $outarray{$j}, "\n"; } close (OUTPUT);}sub by_number { if ($a $b) { -1; } elsif ($a == $b) { 0; } elsif ($a $b) { 1; }}sub process_input { $status = 0; $debug = 0; while ( INPUT ) { $string = $_; print $_ if ( ! /NB-upda/ && $debug ); if (/A V E R A G E S/) { $averages = 1; ($averages_over) = /.*O V E R.*(\d*).*S T E P S/; } $rms = 1 if (/R M S/); if (/NSTEP/) { ($time, $temp, $pres) = /NSTEP =.*TIME.* =(.*\d*\.\d*).*TEMP.* =(.*\d*\.\d*).*PRESS = (.*\d*\.\d*)/; if ( $debug ) { print $_; print "time is $time, temp is $temp, pres is $pres\n"; } $_ = INPUT; if (/Etot/) { ($etot, $ektot, $eptot) = /Etot.*=(.*\d*\.\d*).*EKtot.*=(.*\d*\.\d*).*EPtot.*=(.*\d*\.\d*)/; if ( $debug ) { print $_; print "Etot is $etot, ektot is $ektot, eptot is $eptot\n"; } $_ = INPUT; } if (/BOND.*ANGLE.*DIHED/) { ($bond, $angle, $dihedral) = /BOND.*=(.*\d*\.\d*).*ANGLE.*=(.*\d*\.\d*).*DIHED.*=(.*\d*\.\d*)/; if ( $debug ) { print $_; print "bond is $bond, angle is $angle, dihedral is $dihedral\n"; } $_ = INPUT; } if (/1-4 NB/) { ($nb14, $eel14, $nb) = /1-4 NB.*=(.*\d*\.\d*).*1-4 EEL.*=(.*\d*\.\d*).*VDWAALS.*=(.*\d*\.\d*)/; if ( $debug ) { print $_; print "nb14 is $nb14, eel14 is $eel14, vdwaals is $nb\n"; } $_ = INPUT; } if (/EELEC/) { ($eel, $ehbond, $constraint) = /EELEC.*=(.*\d*\.\d*).*EHBOND.*=(.*\d*\.\d*).*CONSTRAINT.*=(.*\d*\.\d*)/; if ( $debug ) { print $_; print "eel is $eel, ehbond is $ehbond, constraint is $constraint\n"; } $_ = INPUT;## check to see if EAMBER is in the mdout file (present when# NTR=1)# if ( /EAMBER/ ) { $_ = INPUT; } } if (/EKCMT/) { ($ekcmt, $virial, $volume) = /EKCMT.*=(.*\d*\.\d*).*VIRIAL.*=(.*\d*\.\d*).*VOLUME.*=(.*\d*\.\d*)/; if ( $debug ) { print $_; print "Ekcmt is $ekcmt, virial is $virial, volume is $volume\n"; } $_ = INPUT; } if (/T_SOLUTE/) { ($tsolute, $tsolvent) = /T_SOLUTE =(.*\d*\.\d*).*T_SOLVENT =(.*\d*\.\d*)/; if ( $debug ) { print $_; print "Temp solute is $tsolute, temp solvent is $tsolvent\n"; } $_ = INPUT; } if (/Density/) { ($density) = /.*Density.*=(.*\d*\.\d*)/; if ( $debug ) { print $_; print "Density is $density\n"; } $_ = INPUT; } if (/Etot/) { ($etot, $ektot, $eptot) = /Etot.*=(.*\d*\.\d*).*EKtot.*=(.*\d*\.\d*).*EPtot.*=(.*\d*\.\d*)/; if ( $debug ) { print $_; print "Etot is $etot, ektot is $ektot, eptot is $eptot\n"; } $_ = INPUT; } if (/ESCF/) { ($escf) = /.*ESCF.*=(.*\d*\.\d*)/; if ( $debug ) { print $_; print "ESCF is $escf\n"; } $_ = INPUT; }# update arrays if ( $averages == 1 ) { $AVG_TIME{$time} = $time; $AVG_TEMP{$time} = $temp; $AVG_PRES{$time} = $pres; $AVG_ETOT{$time} = $etot; $AVG_EKTOT{$time} = $ektot; $AVG_EPTOT{$time} = $eptot; $AVG_BOND{$time} = $bond; $AVG_ANGLE{$time} = $angle; $AVG_DIHEDRAL{$time} = $dihedral; $AVG_NB14{$time} = $nb14; $AVG_EEL14{$time} = $eel14; $AVG_NB{$time} = $nb; $AVG_EEL{$time} = $eel; $AVG_EHBOND{$time} = $ehbond; $AVG_CONSTRAINT{$time} = $constraint; $AVG_EKCMT{$time} = $ekcmt; $AVG_VIRIAL{$time} = $virial; $AVG_VOLUME{$time} = $volume; $AVG_TSOLUTE{$time} = $tsolute; $AVG_TSOLVENT{$time} = $tsolvent; $AVG_DENSITY{$time} = $density; $AVG_ESCF{$time} = $escf; $averages = 0; } elsif ( $rms == 1 ) { $RMS_TIME{$time} = $time; $RMS_TEMP{$time} = $temp; $RMS_PRES{$time} = $pres; $RMS_ETOT{$time} = $etot; $RMS_EKTOT{$time} = $ektot; $RMS_EPTOT{$time} = $eptot; $RMS_BOND{$time} = $bond; $RMS_ANGLE{$time} = $angle; $RMS_DIHEDRAL{$time} = $dihedral; $RMS_NB14{$time} = $nb14; $RMS_EEL14{$time} = $eel14; $RMS_NB{$time} = $nb; $RMS_EEL{$time} = $eel; $RMS_EHBOND{$time} = $ehbond; $RMS_CONSTRAINT{$time} = $constraint; $RMS_EKCMT{$time} = $ekcmt; $RMS_VIRIAL{$time} = $virial; $RMS_VOLUME{$time} = $volume; $RMS_TSOLUTE{$time} = $tsolute; $RMS_TSOLVENT{$time} = $tsolvent; $RMS_DENSITY{$time} = $density; $RMS_ESCF{$time} = $escf; $rms = 0; } else { $TIME{$time} = $time; $TEMP{$time} = $temp; $PRES{$time} = $pres; $ETOT{$time} = $etot; $EKTOT{$time} = $ektot; $EPTOT{$time} = $eptot; $BOND{$time} = $bond; $ANGLE{$time} = $angle; $DIHEDRAL{$time} = $dihedral; $NB14{$time} = $nb14; $EEL14{$time} = $eel14; $NB{$time} = $nb; $EEL{$time} = $eel; $EHBOND{$time} = $ehbond; $CONSTRAINT{$time} = $constraint; $EKCMT{$time} = $ekcmt; $VIRIAL{$time} = $virial; $VOLUME{$time} = $volume; $TSOLUTE{$time} = $tsolute; $TSOLVENT{$time} = $tsolvent; $DENSITY{$time} = $density; $ESCF{$time} = $escf; } } }}But after running this perl script on the output files I got graph displaying half values of output files as attached here. Can anybody guide me whats wrong with the perl graph or my input file as my output files are printing all values for energy but perl script plot half values??? Regards,
Sehrish Naz
Jr. Research Fellow,Computational Chemistry Unit.
Dr. Panjwani Center for Molecular Medicine and Drug Research,
International Center for Chemical and Biological Sciences,
University of Karachi, Karachi-75270.
E-mail: mjazse.gmail.com
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Mar 19 2016 - 01:30:03 PDT
