Can you try to use mdout_analyzer.py in AmberTools15? (It's pain to debug
Perl code)
Check
http://ambermd.org/doc12/Amber15.pdf (page 515)
Hai
On Sat, Mar 19, 2016 at 1:25 AM, Sehrish Naz Aijaz <sehrish.naz.outlook.com>
wrote:
> Dear all,
> I had run 10 ns simulation (AMBER12) 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/perl
>
>
> if ($#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???
>
>
>
> Sent from Mail for Windows 10
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Mar 19 2016 - 02:00:03 PDT
