Ticket #8587: dos.pl

File dos.pl, 25.5 KB (added by Samuel Jackson, 7 years ago)

Original Perl script

Line 
1#! /usr/bin/env perl
2#
3# Parse a "<>.castep" or "<>.phonon" or "<>.bands" output file
4
5use POSIX;
6use Getopt::Long;
7#use Math::Complex;
8use strict;
9#use warnings;
10
11sub usage {
12   printf STDERR "Usage: dos.pl [-xg] [-ps|-eps] [-np] <seed>.castep|<seed>.phonon ... \n";
13   printf STDERR "       Extract phonon or bandstructure data from .castep, .phonon or .bands files";
14   printf STDERR "       and optionally prepare a DOS plot using XMGRACE as a backend.";
15   printf STDERR "    -xg        Write a script and invoke GRACE to plot data.\n";
16   printf STDERR "    -gp        Write a script and invoke GNUPLOT to plot data.\n";
17   printf STDERR "    -ps        Invoke GRACE to plot data and write as a PostScript file.\n";
18   printf STDERR "    -eps       Invoke GRACE to plot data and write as an encapsulated.\n               PostScript (EPS) file.\n";
19   printf STDERR "    -np        Do not plot data, write a GRACE script.\n";
20   printf STDERR "    -bs        Read band-structure from <>.castep or <>.bands.\n";
21   printf STDERR "    -mirror    Plot spin-polarized electronic DOS using \"mirror\" plot.\n";
22   printf STDERR "    -b w       Set histogram resolution for binning (eV or cm**-1).\n";
23   printf STDERR "    -ir        Extract ir intensities and model (fundamentals-only) ir spectrum from .phonon.\n";
24   printf STDERR "    -raman     Extract raman intensities and model (fundamentals-only) raman spectrum from .phonon.\n";
25   printf STDERR "    -temp T    Temperature to use (in raman spectrum modelling).\n";
26   printf STDERR "    -expt FILE Read experimental data from EXPT and overplot.\n";
27   printf STDERR "    -dat       Reread standard output from previous run and plot.\n";
28   printf STDERR "    -w s       Set Gaussian/Lorentzian FWHM for broadening.\n";
29   printf STDERR "    -lorentz   Use Lorentzian broadening instead of Gaussian\n";
30   printf STDERR "    -v         Be verbose about progress\n";
31   printf STDERR "    -z         Print zero-point energy\n";
32die;
33}
34
35my $number  = qr/-?(?:\d+\.?\d*|\d*\.?\d+)(?:[Ee][+-]?\d{1-3})?/o;
36my $fnumber = qr/-?(?:\d+\.?\d*|\d*\.?\d+)/o;
37
38my $qtol = 5.0e-6;
39my $fzerotol = 3.0;
40my $bwidth = 1.0;
41my $gwidth = 10.0;
42my $pi=3.14159265358979;
43
44my $castep_phonon_re = qr/^ +\+ +q-pt= +\d+ \( *( *$fnumber) *($fnumber) *($fnumber)\) +($fnumber) +\+/o;
45my $dot_phonon_re    = '';
46
47my $castep_bs_re = qr/^  \+ +Spin=(\d) kpt= +(\d+) \( *( *$fnumber) *($fnumber) *($fnumber)\) +kpt-group= +\d +\+/o;
48my $dot_bands_re    = '';
49
50my $castep_re;
51
52my ($grflag, $gpflag, $plotflag, $datflag,  $irflag, $ramanflag, $psflag, $qlabels, $exptfile, $exptq, $abscissa,$verbose) = ("", "", "","","","","","","");
53my ($title, $fileroot, $xlabel, $ylabel, $fermi_u, $fermi_d, $i,$freq,$exptn, $wt, $b, $temperature, $lorentzian, $zpe, $nqlast,$spin_polarised);
54
55my (@freqs, @freqs_d, @qpts, @weights, @dos, @dos_u, @dos_d, @ir_intensities, @r_intensities, $base);
56
57my ($opt_xg, $opt_gp, $opt_np, $opt_ps, $opt_eps, $opt_bs, $opt_mirror, $opt_bw,$opt_ir, $opt_raman, $opt_temp, $opt_lorentz,  $opt_dat, $opt_expt, $opt_gwidth,$opt_help,$opt_v,$opt_z) = (0,0,0,0,0,0,0,"",0,0,0,0);
58
59&GetOptions("xg"   => \$opt_xg, 
60            "gp"   => \$opt_gp, 
61            "np"   => \$opt_np, 
62            "ps"   => \$opt_ps, 
63            "eps"  => \$opt_eps, 
64            "bs"   => \$opt_bs, 
65            "mirror"   => \$opt_mirror, 
66            "b=f"  => \$opt_bw,
67            "ir"   => \$opt_ir, 
68            "raman" => \$opt_raman, 
69            "temp=f" => \$opt_temp,
70            "dat"  => \$opt_dat, 
71            "expt=s" => \$opt_expt,
72            "lorentz" => \$opt_lorentz, 
73            "w=s"  => \$opt_gwidth, 
74            "h|help"    => \$opt_help,
75            "v"    => \$opt_v,
76            "z"    => \$opt_z) || &usage();
77
78if ($opt_help || $opt_ps && $opt_eps) {&usage()};
79
80$grflag= $opt_xg if $opt_xg;
81$gpflag= $opt_gp if $opt_gp;
82$plotflag = 1;
83$plotflag= 0 if $opt_np;
84$psflag = 0;
85$psflag= 1 if $opt_ps;
86$psflag= 2 if $opt_eps;
87$grflag = 1 if $psflag && ! ($grflag || $gpflag) ;
88$exptfile = "";
89$exptfile = $opt_expt if $opt_expt;
90$datflag = 1 if $opt_dat;
91$gwidth = $opt_gwidth if $opt_gwidth;
92$verbose = 0;
93$verbose = 1 if $opt_v;
94$opt_bs = 1 if( $#ARGV >= 0 && $ARGV[0]=~ /\.bands$/ );
95$opt_bs = 1 if $opt_mirror;
96$opt_bs = 0 if( $#ARGV >= 0 && $ARGV[0]=~ /\.phonon$/ );
97$irflag = $opt_ir;
98$ramanflag = $opt_raman;
99$temperature = 300;
100$temperature=$opt_temp if $opt_temp;
101$lorentzian = $opt_lorentz;
102
103if ($#ARGV >= 0) {
104  $title=$ARGV[0];
105  $title=~s/\.(phonon|castep|bands)//;
106} else {
107  $title = "";
108}
109$fileroot=$title;
110
111if( $opt_bs ) {
112  $castep_re = $castep_bs_re;
113  if( $grflag ) {
114    $xlabel = '\\f{Symbol}e\\f{} (eV)';
115    $ylabel = 'g(\\f{Symbol}e\\f{})';
116  } elsif ($gpflag ) {
117    $xlabel = '{/Symbol e} (eV)';
118    $ylabel = 'g({/Symbol e})';
119  }
120  $gwidth = 0.1;
121  $gwidth = $opt_gwidth if $opt_gwidth;
122  $bwidth = 0.05;
123  $fzerotol = 0.0;
124} else {
125  $castep_re = $castep_phonon_re;
126  if( $grflag ) {
127    $xlabel = '\\f{Symbol}w\\f{} (cm\\S-1\\N)';
128    $ylabel = 'g(\\f{Symbol}w\\f{})';
129  } elsif ($gpflag ) {
130    $xlabel = '{/Symbol w} (cm^{-1})';
131    $ylabel = 'g({/Symbol w})';
132  }
133}
134$bwidth = $opt_bw if( $opt_bw );
135
136
137if ($plotflag ) {
138  if( $grflag) {
139    if( $psflag ) {
140      open STDOUT, "| gracebat -pipe -nosafe" or die "Failed to open pipe to GRACEBAT process";
141    } else {
142      open STDOUT, "| xmgrace -pipe" or die "Failed to open pipe to GRACE process";
143    }
144  } elsif ($gpflag) {
145    if( $psflag ) {
146      open STDOUT, "| gnuplot -persist - 2>/dev/null 1>&2" or die "Failed to open pipe to GNUPLOT process";
147    } else {
148      open STDOUT, "| gnuplot -persist - 2>/dev/null 1>&2" or die "Failed to open pipe to GNUPLOT process";
149    }
150  }
151}
152
153
154#$pos = tell ARGV;
155$_ = <>;
156seek ARGV,0,0;
157if ($datflag) {
158  reread_data(\@freqs,\@qpts);
159} elsif ( /^ \+-+\+\s*$/) {
160  if(  $opt_bs ) {
161    read_dot_castep_bands(\@freqs,\@freqs_d,\@qpts, \$fermi_u, \$fermi_d, \@weights);
162  } else {
163    read_dot_castep_phonon(\@freqs,\@qpts, \@weights, \@ir_intensities, \@r_intensities);
164  }
165} else {
166  if(  $opt_bs ) {
167    read_dot_bands(\@freqs,\@freqs_d,\@qpts, \$fermi_u, \$fermi_u, \@weights);
168  } else {
169    read_dot_phonon(\@freqs,\@qpts, \@weights, \@ir_intensities, \@r_intensities);
170  }
171}
172
173$spin_polarised = $#{$freqs_d[0]} > 0;
174
175$base = 1.0e21;
176if( $irflag ) {
177  compute_dos(\@freqs, \@weights, \@dos_u, $bwidth, $gwidth, \$base, \@ir_intensities);
178} elsif( $ramanflag ) {
179  compute_raman(\@freqs, \@weights, \@dos_u, $bwidth, $gwidth, \$base, \@r_intensities);
180} else {
181  $zpe = compute_zpe(\@freqs, \@weights);
182  printf STDERR "Zero Point Energy = %12.6f eV\n", $zpe if $opt_z;
183  compute_dos(\@freqs, \@weights, \@dos_u, $bwidth, $gwidth, \$base);
184  compute_dos(\@freqs_d, \@weights, \@dos_d, $bwidth, $gwidth, \$base) if ($opt_bs && $spin_polarised) ;
185}
186
187@dos = (\@dos_u, \@dos_d);
188
189if( $grflag ) {
190  if(  $opt_bs ) {
191    graceout(\@dos, $base, $bwidth, $fileroot, $psflag, $xlabel, $ylabel, $title, $fermi_u);
192  } else {
193    graceout(\@dos, $base, $bwidth, $fileroot, $psflag, $xlabel, $ylabel, $title,"" ,$exptn);
194  }
195} elsif( $gpflag ) {
196  if(  $opt_bs ) {
197    gnuplotout(\@dos, $base, $bwidth, $fileroot, $psflag, $xlabel, $ylabel, $title, $fermi_u);
198  } else {
199    gnuplotout(\@dos, $base, $bwidth, $fileroot, $psflag, $xlabel, $ylabel, $title,"" ,$exptn);
200  }
201
202} else {
203  if ($opt_bs && $spin_polarised > 0) {
204    for $b (0..$#dos_u) {
205      printf "%12.6f %12.6f %12.6f\n", $b*$bwidth+$base, $dos_u[$b],$dos_d[$b];
206    }
207  } else {
208
209    for $b (0..$#dos_u) {
210      printf "%12.6f %12.6f\n", $b*$bwidth+$base, $dos_u[$b];
211    }
212  }
213}
214
215#
216# End of executable part -- subroutines follow.
217#
218
219sub reread_data {
220  my($freqs, $qpts) = @_;
221  my( $count, $pt);
222  my(@freqk,@q,@line, $nfreqs, $pltcoord );
223
224  $count = 0;
225
226  while ( <> ) {
227    #
228    # <>.txt
229    #
230    print STDERR $_ if $verbose;
231    @line = split;
232   
233    $pt = shift @line;
234    @q = splice @line, 0, 3;
235    $pltcoord = shift @line;
236    @freqk = @line;
237
238    die "Need sequentially numbered lines ($pt != $count)" unless $pt == $count;
239    if( $count > 1 ) {
240      die "Every line should have same number of values" unless $nfreqs == $#freqk;
241    }
242   
243    $count++;
244    $nfreqs = $#freqk;
245    push @$freqs, [@freqk];
246    push @$qpts, [@q];
247   
248  }
249}
250
251sub read_dot_castep_phonon {
252  my($freqs, $qpts, $weights, $intens, $rintens) = @_;
253  my(@freqk,@q,@intensk,@rintensk);
254  my ($weight);
255
256  while ( <>) {
257    #
258    # <>.castep
259    #
260    if (/$castep_phonon_re/o) {
261      @q=($1,$2,$3);
262      $weight = $4;
263      print STDERR "Reading $q[0] $q[1] $q[2] $weight\n" if $verbose;
264      #print "   ";
265      last if eof ARGV;
266      $_ = <>;
267      last if eof ARGV;
268      $_ = <> if /^ +\+ +Effective cut-off =/;
269      last if eof ARGV;
270      $_ = <> if /^ +\+ +q->0 along \( *( *$fnumber) *($fnumber) *($fnumber)\) +\+/;
271      last if eof ARGV;
272      $_ = <> if /^ +\+ -+ \+/;
273      @freqk=();
274      @intensk=();
275      @rintensk=();
276
277      while (<>) {
278        if (/^ +\+ +\d+ +($fnumber)( +\w)? *($fnumber)?( +[YN])? *($fnumber)?( +[YN])? *\+/) { 
279          push @freqk, $1;
280          push @intensk, $3;
281          push @rintensk, $5;
282#       } elsif (/^  \+ +q->0 along \( *( *$fnumber) *($fnumber) *($fnumber)\) +\+/ ) {
283        } elsif (/^ +\+ -+ \+/) {
284          last;
285        } elsif (/^ +\+ +.*\+/) {
286        } else {
287          last;
288        }
289      }
290      push @$freqs, [@freqk];
291      push @$intens, [@intensk];
292      push @$rintens, [@rintensk];
293      push @$qpts, [@q];
294      push @$weights, $weight;
295      last if eof ARGV;
296    }
297  }
298}
299
300sub read_dot_castep_bands {
301  my($freqs_u, $freqs_d, $qpts, $fermi_u, $fermi_d, $weights) = @_;
302  my(@freqk,@q,@kp_map,$kp_count,$spin,$kpt,$weight,$nk);
303 
304  $kp_count = 0;
305  while ( <>) {
306    #
307    # <>.castep
308    #
309    $$fermi_u = $1 if /^  \+  Fermi energy for spin  up  electrons is:\s+ ($number)[^+]+\+/;
310    $$fermi_d = $1 if /^  \+  Fermi energy for spin down electrons is:\s+ ($number)[^+]+\+/;
311    if (/$castep_bs_re/o) {
312      ($spin,$kpt,@q,$weight)=($1,$2,$3,$4,$5,$6);
313      print STDERR "Reading Spin=$spin;kpt=$kpt  ($q[0],$q[1],$q[2]) $weight\n" if $verbose;
314      #print "   ";
315      last if eof ARGV;
316      $_ = <>;
317      last if eof ARGV;
318      $_ = <> if /^  \+ -+ \+/;
319      @freqk=();
320
321      while (<>) {
322        if (/^  \+ +\d+ +($fnumber)( +\w)?( +$fnumber)* *\+/) { # This matches exptl IR -active format
323          push @freqk, $1;
324#       } elsif (/^  \+ +q->0 along \( *( *$fnumber) *($fnumber) *($fnumber)\) +\+/ ) {
325        } elsif (/^  \+ -+ \+/) {
326          last;
327        } elsif (/^  \+ +.*\+/) {
328        } else {
329          last;
330        }
331      }
332      print STDERR @freqk,"\n" if $verbose;
333      #
334      # Keep track of kpt ordering for spin=2.  Assumes all spin=1 precede spin=2.
335      #
336      if( $spin == 1 ) {
337        $kp_map[$kpt] = $kp_count++;
338        push @$freqs_u, [@freqk];
339        push @$qpts, [@q];
340        push @$weights, 1.0;
341      } else {
342        push @{$freqs_d[$kp_map[$kpt]]}, @freqk;
343      }
344      last if eof ARGV;
345    }
346  }
347  for $nk (0..$#weights) {
348    $weights[$nk] = 1.0/$kp_count;
349  }
350}
351
352sub min {
353  my ($a,$b) = @_;
354  return $a if $a < $b;
355  return $b;
356}
357
358sub max {
359  my ($a,$b) = @_;
360  return $a if $a > $b;
361  return $b;
362}
363
364sub read_dot_phonon {
365  my($freqs, $qpts, $weights, $intens, $rintens) = @_;
366  my(@freqk,@freqk_r,@freqk_o,@q,@eigvec, @oeigvec, @qmapt, @qmap, @intensk, @rintensk);
367  my(@coords,@xyz,@corrmat);
368  my($x,$y,$z,$nions,$nmodes, $ion,$aty, $mass, $count, $mode, $omode,$dot);
369  my($qmode,$qomode,$dot_same,$savedot,$maxdot,$save,$qmaxdot,$weight);
370
371  #
372  # Read frequencies from <>.phonon file
373  #
374  # Also reads eigenvectors and attempts to join up branches
375  # on the dispersion curve using a heuristic algorithm based
376  # assuming close to orthonormality of eigenvectors at
377  # adjacent q-points.
378  #
379  @$weights = ();
380  @$qpts    = ();
381  @$freqs   = ();
382  @$intens  = ();
383  while (<>) {
384    #
385    # <>.phonon
386    #
387    $nions = $1 if /^ Number of ions\s*(\d+)/;
388    $nmodes = $1 if /^ Number of branches\s*(\d+)/;
389    if( / Fractional Co-ordinates/../ end.*header/i) {
390      if( /^\s+\d+(\s+$fnumber){3}\s+[A-Za-z]{1,2}\s+$fnumber/ ) {
391        ($ion,$x,$y,$z,$aty,$mass) = split;
392        push @coords, [$x,$y,$z];
393      }
394    }
395    #
396    # Found start of block of frequencies
397    #
398    if (/^ +q-pt= +(\d+) +($fnumber) +($fnumber) +($fnumber) *($fnumber){0,4}/) {
399#    if (/^ +q-pt= +\d+ +($fnumber) *($fnumber) *($fnumber) ( *$fnumber){0,4}$/) {
400      @q=($2,$3,$4);
401      $weight = $5;
402      #
403      # Test for repeated q-point and assign zero weight to copies
404      #
405      $weight = 0.0 if( $1 == $nqlast );
406      $nqlast = $1;
407      print STDERR " ",$q[0]," ",$q[1]," ",$q[2]," ",$weight,"\n" if $verbose;
408      #print "   ";
409      @intensk=();
410      @rintensk=();
411      @freqk_r=();
412
413      while (<>) {
414        if (/^ +\d+ +($fnumber) +($fnumber) +($fnumber)/) {
415          push @freqk_r, $1;
416          push @intensk, $2;
417          push @rintensk, $3;
418        } elsif (/^ +\d+ +($fnumber) +($fnumber)/) {
419          push @freqk_r, $1;
420          push @intensk, $2;
421        } elsif (/^ +\d+ +($fnumber) */) {
422          push @freqk_r, $1;
423        } else {
424          last;
425        }
426      }
427#      print STDERR @freqk_r,"\n";
428      push @$qpts, [@q];
429      push @$freqs, [@freqk_r];
430      push @$intens, [@intensk];
431      push @$rintens, [@rintensk];
432      push @$weights, $weight;
433    }
434  }
435}
436
437sub read_dot_bands {
438  my($freqs_u, $freqs_d, $qpts, $fermi, $fermi1, $weights) = @_;
439  my(@eigk,@k, @freq_list, @qpt_list, @nk_list, @wt_list);
440  my($wt, $first, $Hartree,$nk,$ink,$ink_save, $ncount, $cpt, $spin_polarised);
441
442  $Hartree =27.2114;
443  $first = 1;
444  $nk = -1;
445  $ncount = 0;
446  $spin_polarised = 0;
447  $cpt = 0;
448  #
449  # Pass 1 - read the data
450  while (<>) {
451    #
452    # <>.bands
453    #
454    chop;
455    if (/^K-point +(\d+) +($fnumber) *($fnumber) *($fnumber) *($fnumber)/) {
456      if( ! $first ) {
457        push @{$freq_list[$cpt]}, [@eigk];
458        push @qpt_list,  [@k];
459        push @nk_list,   $nk;
460        push @wt_list,   $wt;
461        print STDERR "Read K-point $nk, ($qpt_list[$ncount-1][0],$qpt_list[$ncount-1][1],$qpt_list[$ncount-1][2])\n" if $verbose;
462      }
463      $first = 0;
464      if( $1 == 1 and $ncount > 1 ) {
465        if( $ncount == $nk+1) {
466          $nk = $ncount;
467        } else {
468          print STDERR "Don't know what to do with duplicate k-point $ncount\n";
469          $nk = -1;
470        }
471      } else {
472        $nk = $1-1;
473      }
474      @k=($2,$3,$4);
475      $wt = $5;
476      @eigk=();
477      $ncount++;
478      $cpt = 0;
479    } elsif (/^Spin component 1/) {
480    } elsif (/^Spin component 2/) {
481      push @{$freq_list[0]}, [@eigk];
482      print STDERR "Read K-point $nk, ($qpt_list[$ncount-1][0],$qpt_list[$ncount-1][1],$qpt_list[$ncount-1][2]), spin 2\n" if $verbose;
483      @eigk=();
484
485      $cpt = 1;
486      $spin_polarised = 1;
487    } elsif (/^ +($fnumber)\s*$/) {
488      push @eigk, $1*$Hartree;
489    } elsif ( /Fermi energy \(in atomic units\) +($fnumber)/ ) {
490      $$fermi = $1 * $Hartree;
491    } elsif ( /Fermi energies \(in atomic units\) +($fnumber) +($fnumber)/ ) {
492      $$fermi = $1 * $Hartree;
493      $$fermi1 = $1 * $Hartree;
494    }
495  }
496  push @{$freq_list[$cpt]}, [@eigk];
497  if( $cpt - 0 > 0 ) { # Force conversion to F.P.
498    #
499    # Already have k-point and weight if this is second spin component
500    #
501    push @qpt_list,  [@k];
502    push @nk_list,   $nk;
503    push @wt_list,   $wt;
504  }
505
506  print STDERR "Read K-point $nk, ($qpt_list[$ncount-1][0],$qpt_list[$ncount-1][1],$qpt_list[$ncount-1][2])\n" if $verbose;
507  #
508  # Reorder list and store into output arrays
509  #
510  $ink_save = -1;
511  for $ink (0..$#nk_list) {
512    $nk = $nk_list[$ink];
513    print STDERR "Reordering K-point $nk, ($qpt_list[$ink][0],$qpt_list[$ink][1],$qpt_list[$ink][2])\n" if $verbose;
514    if( $nk < 0 ) {
515      $ink_save = $ink;
516    } else {
517      $$freqs_u[$nk] = $freq_list[0][$ink];
518      $$freqs_d[$nk] = $freq_list[1][$ink] if $spin_polarised;
519     
520      $$qpts[$nk] =  $qpt_list[$ink];
521      $$weights[$nk] = $wt_list[$ink];
522    }
523  }
524  print STDERR "Spin_polarised=",$spin_polarised,"\n";
525  #
526  # Due to a CASTEP bug some state may be double-labelled. In that case we have a "spare" set
527  # of frequencies in $ink_save.  See if we can guess where to put it.
528  for $ink (1..$#$qpts) {
529    if( $#{$$freqs_u[$ink]} < 1 ) {  # Aha!
530    print STDERR "Placing missing K-point $ink, ($qpt_list[$ink_save][0],$qpt_list[$ink_save][1],$qpt_list[$ink_save][2])\n" if $verbose;
531     
532      $$freqs_u[$ink] = $freq_list[0][$ink_save];
533      $$freqs_d[$ink] = $freq_list[1][$ink_save] if $spin_polarised;
534      $$qpts[$ink] =  $qpt_list[$ink_save];
535    }
536  }
537
538}
539
540sub ftrim {
541  my($val,$prec) = @_;
542  my($str,$w);
543 
544  $w=$prec+2;
545  $str = sprintf "%${w}.${prec}f",$val;
546  $str =~ s/0+$//;
547  $str;
548}
549
550sub graceout {
551    my($dos, $base, $bwidth, $plotfile, $psflag, $xlabel, $ylabel, $title, $fermi, $fexptl ) = @_;
552    my($i,  $setnum);
553    my($xmax, $ymax, $n, $ncpts, $sign);
554
555    $| = 1;
556    $xmax = $base+$bwidth*$#{$dos->[0]};
557
558    print   "\@timestamp def \"",scalar(localtime),"\"\n";
559    print   "\@with g0\n";
560    print   "\@title \"",$title,"\"\n";
561    print   "\@world xmin $base\n";
562    print   "\@world xmax $xmax\n";
563    print   "\@view  ymin 0.35\n";
564    print   "\@view  xmax 0.75\n";
565    print   "\@legend  0.625, 0.825\n";
566    print   "\@autoscale onread xyaxes\n";
567   
568    print   "\@xaxis label \"$xlabel\"\n";
569    print   "\@yaxis label \"$ylabel\"\n";
570
571    $setnum=0;
572    if( $#{$dos->[1]} > 0 ) {
573      print "\@ G0.S",$setnum," line color 2\n";
574      print "\@ G0.S",$setnum," legend \"\\f{Symbol}a\\f{}\"\n";
575    }
576    printf STDOUT "\@target G0.S%d\n",$setnum;
577    print "\@type xy\n";
578
579    $ymax = 0.0;
580    foreach $n (0..$#{$dos->[0]}) {
581      printf STDOUT "%15.6f %15.6f\n", $n*$bwidth+$base, ${$dos->[0]}[$n];
582      $ymax =  ${$dos->[0]}[$n] if  ${$dos->[0]}[$n] > $ymax;
583    }
584    print "&\n";
585    #
586    # 2nd component?
587    #
588    if( $#{$dos->[1]} > 0 ) {
589      $setnum++;
590
591      $sign = 1;
592      if( $opt_mirror ){
593        $sign = -1;
594
595        print "\@altxaxis  on\n";
596        print "\@altxaxis  type zero true\n";
597        print "\@altxaxis  tick off\n";
598        print "\@altxaxis  ticklabel off\n";
599
600      }
601
602      print "\@ G0.S",$setnum," line color 4\n";
603      print "\@ G0.S",$setnum," legend \"\\f{Symbol}b\\f{}\"\n";
604      printf STDOUT "\@target G0.S%d\n",$setnum;
605      print "\@type xy\n";
606
607      foreach $n (0..$#{$dos->[1]}) {
608        printf STDOUT "%15.6f %15.6f\n", $n*$bwidth+$base, $sign*${$dos->[1]}[$n];
609        $ymax =  ${$dos->[1]}[$n] if  ${$dos->[1]}[$n] > $ymax;
610      }
611      print "&\n";
612
613    }
614    $setnum++;
615
616    if( $fermi ne "" ) {
617      print "\@ G0.S",$setnum," line linestyle 3\n";
618      print "\@ G0.S",$setnum," line color 1\n";
619      print "\@ G0.S",$setnum," legend \"\\f{Symbol}e\\f{}\\sF\\N\"\n";
620      printf STDOUT "\@target G0.S%d\n",$setnum;
621      print "\@type xy\n";
622      if( $opt_mirror ) {
623        printf STDOUT "%15.6f %15.6f\n",$fermi,-$ymax;
624      } else {
625        printf STDOUT "%15.6f %15.6f\n",$fermi,0;
626      }
627      printf STDOUT "%15.6f %15.6f\n",$fermi,$ymax;
628      print "&\n";
629    }
630    if( $psflag == 1){
631      print   "\@hardcopy device \"PS\"\n";
632      printf  "\@print to \"%s.ps\"\n\@print\n", $plotfile;
633    }
634    if( $psflag == 2){
635      print   "\@hardcopy device \"EPS\"\n\@device \"EPS\" op \"bbox:tight\"\n";
636      printf  "\@print to \"%s.eps\"\n\@print\n", $plotfile;
637    }
638}
639
640sub gnuplotout {
641    my($dos, $base, $bwidth, $plotfile, $psflag, $xlabel, $ylabel, $title, $fermi, $fexptl ) = @_;
642    my($i,   $setnum);
643    my($xmax, $ymax, $n, $sign);
644
645    $| = 1;
646    $xmax = $base+$bwidth*$#{$dos->[0]};
647    $xmax = $fermi + 1 if ($fermi ne "" && $fermi + 1 > $xmax) ;
648
649    if( $psflag == 1){
650      print   "set terminal postscript landscape color solid\n";
651      printf   "set output \"%s.ps\"\n", $plotfile;
652    }
653    if( $psflag == 2){
654      print   "set terminal postscript eps color solid\n";
655      printf   "set output \"%s.eps\"\n", $plotfile;
656    }
657    print   "set data style lines\n";
658    print   "set termoption enhanced\n";
659    print   "set title \"",$title,"\"\n";
660    print   "set xlabel \"$xlabel\"\n";
661    print   "set ylabel \"$ylabel\"\n";
662    print   "set xrange  [$base:$xmax]\n";
663
664    if( $#{$dos->[1]} > 0 ) {
665      print   "plot '-' title '{/Symbol a}', '-' title '{/Symbol b}' lt 3" ;
666
667      if( $opt_mirror ){
668        print ", '-' title '' lt 1 lc 'black'";
669      }
670    } else {
671      print   "plot '-' title '$ylabel' lt 1 lc 'black'" ;
672    }
673    if( $fermi ne "" ) {
674      print  ", '-' title '{/Symbol e}_{/*0.5 F}' lt 0 lc 'black'" ;
675    }
676    print "\n";
677    #
678    # Write datasets
679    #
680    $ymax = 0;
681    foreach $n (0..$#{$dos->[0]}) {
682      printf "%15.6f %15.6f\n", $n*$bwidth+$base, ${$dos->[0]}[$n];
683      $ymax =  ${$dos->[0]}[$n] if  ${$dos->[0]}[$n] > $ymax;
684    }
685    print "end\n";
686    if( $#{$dos->[1]} > 0 ) {
687      $sign = 1;
688      $sign = -1.0 if( $opt_mirror );
689
690      foreach $n (0..$#{$dos->[1]}) {
691        printf "%15.6f %15.6f\n", $n*$bwidth+$base, $sign*${$dos->[1]}[$n];
692        $ymax =  ${$dos->[1]}[$n] if  ${$dos->[1]}[$n] > $ymax;
693      }
694      print "end\n";
695    }
696    #
697    # Data for fermi level indicator
698    #
699    if( $fermi ne "" ) {
700      if( $opt_mirror ) {
701        printf  "%15.6f %15.6f\n",$base, 0;
702        printf  "%15.6f %15.6f\n",$xmax, 0;
703        print "end\n";
704
705        printf  "%15.6f %15.6f\n",$fermi,-$ymax;
706      } else {
707        printf  "%15.6f %15.6f\n",$fermi,0;
708      }
709      printf  "%15.6f %15.6f\n",$fermi,$ymax;
710    }
711    print "end\n";
712
713    print   "\npause -1\n";
714
715}
716
717
718sub write_data {
719  my($plotfd, $freqs, $qpts, $qflg) = @_;
720  my($n,$q,$freqk);
721  $n = 0;
722  for $q (0 .. $#$freqs) {
723    if( $qflg ) {
724      printf $plotfd "%8d    %8.4f %8.4f %8.4f      ",$q,$$qpts[$q][0],$$qpts[$q][1],$$qpts[$q][2];
725    } else {
726      printf $plotfd "%5d",$n++;
727    }
728
729    for $freqk (@{$$freqs[$q]}) {
730      printf $plotfd "  %12.6f", $freqk;
731    }
732    printf $plotfd "\n";
733  }
734}
735
736#sub dot_product {
737#  my ($a, $b) = @_;
738# 
739#  ${$a}[0]*${$b}[0]+${$a}[1]*${$b}[1]+${$a}[2]*${$b}[2];
740#}
741
742sub dot_product {
743  my ($a, $b) = @_;
744  my ($i,$n,$sum);
745
746  $n = $#{$a};
747  for $i (0..$n) {
748    $sum +=  ${$a}[$i]*${$b}[$i];
749  }
750  $sum;
751}
752
753
754sub read_exptl {
755  my ($file) = @_;
756  my (@expdata, @thisline, @qpoint, @dataset);
757  my ($f, $nsets);
758
759  open EXPT, "<$file" or die "Failed to open experimental data file $file";
760
761  $nsets = 0;
762  $expdata[0] = [0,0,0,0];
763  while (<EXPT>) {
764    @thisline = split;
765    if( $#thisline < 0 ) {
766      $nsets++;
767      @dataset = ();
768    } elsif( $#thisline < 3 ) {
769      die "Need at least 4 items on line (qx,qy,qz,f)\n@thisline";
770    }
771    @qpoint = splice @thisline,0, 3;
772    while ( $f = shift @thisline) {
773      unshift @dataset, [@qpoint, $f] if ($f =~ /$number/);
774    }
775    $expdata[$nsets] = [@dataset];
776  }
777  close EXPT;
778  [@expdata];
779}
780
781sub compute_zpe {
782  my ($freqs, $weights) = @_;
783  my ($nq, $mode);
784  my $zpe  = 0.0;
785  my $cm_to_eV = 1.239842e-4;
786
787  for $nq (0..$#{$freqs}) {
788#    printf STDERR "Q-pt %d weight %f\n",$nq+1,$$weights[$nq];
789    for $mode (0..$#{$$freqs[$nq]}) {
790      $freq = ${$freqs}[$nq]->[$mode];
791#      printf STDERR "Mode %d freq %f\n",$mode+1,$freq;
792      $zpe += $freq*$$weights[$nq];
793    }
794  }
795
796  $zpe *= 0.5*$cm_to_eV;
797  $zpe;
798}
799
800sub compute_dos {
801  my ($freqs, $weights, $dos, $bwidth, $gwidth, $baseptr, $intensities) = @_;
802
803  my ($nq, $q, $mode, $freq, $bin, $base, $g,$h, $ngauss, $nlorentz, $intensity,$gammaby2,$sigma);
804  my (@hist);
805
806  if ( $$baseptr > 1.0e20 ) {
807    $base = $$baseptr;
808
809    for $nq (0..$#{$freqs}) {
810      for $mode (0..$#{$$freqs[$nq]}) {
811        $freq = ${$freqs}[$nq][$mode];
812        $base = $freq if $freq < $base;
813      }
814    }
815    $$baseptr = $base;
816  } else {
817    $base = $$baseptr;
818  }
819
820  for $nq (0..$#{$freqs}) {
821    for $mode (0..$#{$$freqs[$nq]}) {
822      $freq = ${$freqs}[$nq]->[$mode];
823      $bin = ($freq-$base)/$bwidth;
824      $intensity = 1.0;
825      $intensity = 0.0 if abs($freq) < $fzerotol;
826      $intensity =  ${$intensities}[$nq]->[$mode] if ref $intensities ne "";
827      $hist[$bin] += $$weights[$nq]*$intensity;
828    }
829  }
830
831  $ngauss = 3.0*$gwidth/$bwidth;  # 3 sd should do
832  $nlorentz = 25.0*$gwidth/$bwidth;  # 5 widths?
833
834  $gammaby2 = $gwidth/2;
835  $sigma = $gwidth/2.354;
836
837   for $h (0..$#hist) {
838     if( $lorentzian ) { 
839       for $g (-$nlorentz..$nlorentz) {
840         $$dos[$h+$g] += $hist[$h]*$gammaby2/($g**2+$gammaby2**2)/$pi if ($h+$g >= 0);
841       }
842     } else {
843       for $g (-$ngauss..$ngauss) {
844         $$dos[$h+$g] += $hist[$h]*exp( - ($g*$bwidth)**2/(2*$sigma**2))/ (sqrt(2*$pi)*$sigma) if ($h+$g >= 0);
845       }
846     }
847  }
848
849  if ( $verbose ) {
850    for $b (0..$#$dos) {
851      printf "%d %4d %12.4f\n", $b, $hist[$b], $$dos[$b];
852    }
853  }
854}
855
856sub compute_raman {
857  my ($freqs, $weights, $dos, $bwidth, $gwidth, $baseptr, $intensities) = @_;
858
859  my $c = 299792458;
860  my $laser_wavelength = 514.5e-9; # Ar at 514.5 nm
861  my $planck = 6.6260755e-34;
862  my $cm1k = 1.438769; # cm(-1) => K conversion
863  my ($n,$freq,$mode, $factor, $bose_occ, $cross_section, @cross_secs);
864
865  $factor = pow(2*$pi/$laser_wavelength,4)*$planck/(8*$pi**2*45)*1e12;
866
867#  print STDERR "MULT=$factor\n";
868
869  @cross_secs = ();
870  $n=$#{$$freqs[0]};
871  for $mode (0..$n) {
872    $freq = ${$freqs}[0][$mode];
873    $bose_occ = 1.0/(exp($cm1k*$freq/$temperature)-1);
874    if( $freq > $fzerotol ){
875      $cross_section = $factor/$freq * (1 + $bose_occ)*$$intensities[0][$mode];
876    } else {
877      $cross_section = 0;
878    }
879    push @cross_secs, $cross_section;
880#    print STDERR "f=$freq;  occ=$bose_occ; act=$$intensities[0][$mode]; ds/dO= $cross_section\n";
881  }
882  compute_dos($freqs, $weights, $dos, $bwidth, $gwidth, $baseptr, [[@cross_secs]]);
883 
884}