Merge branch 'master' into hse-dev-local

Conflicts:
	src/ExchangeOperator.C
	src/ExchangeOperator.h
	src/Makefile
	src/XCOperator.C
	src/XCPotential.C
	src/Xc.h
	src/release.C

git-svn-id: http://qboxcode.org/svn/qb/branches/hse-dev@1426 cba15fb0-1239-40c8-b417-11db7ca47a34

src/ExchangeOperator.C

@@ -2130,7 +2130,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
     }
     // dc now contains the forces
 
-    // "divergence" corrections (only truly divergent for Coulomb potential)
+    // "divergence" corrections
 
     // correct the energy of state i
     for ( int i = 0; i < sd.nstloc(); i++ )

src/release.C

@@ -19,5 +19,5 @@
 #include "release.h"
 std::string release(void)
 {
-  return std::string("1.57.15");
+  return std::string("wrk");
 }

util/cell.plt

+#!/bin/bash
+# cell.plt: plot cell parameters during an MD simulation
+gnuplot -persist <<EOF
+plot "<grep -h -A3 '<unit_cell' $* | grep a=" u 2 w l, \
+     "<grep -h -A3 '<unit_cell' $* | grep a=" u 3 w l, \
+     "<grep -h -A3 '<unit_cell' $* | grep a=" u 4 w l, \
+     "<grep -h -A3 '<unit_cell' $* | grep b=" u 2 w l, \
+     "<grep -h -A3 '<unit_cell' $* | grep b=" u 3 w l, \
+     "<grep -h -A3 '<unit_cell' $* | grep b=" u 4 w l, \
+     "<grep -h -A3 '<unit_cell' $* | grep c=" u 2 w l, \
+     "<grep -h -A3 '<unit_cell' $* | grep c=" u 3 w l, \
+     "<grep -h -A3 '<unit_cell' $* | grep c=" u 4 w l, 0
+fit ax "<grep -h -A3 '<unit_cell' $* | grep a=" u 0:2 via ax
+fit ay "<grep -h -A3 '<unit_cell' $* | grep a=" u 0:3 via ay
+fit az "<grep -h -A3 '<unit_cell' $* | grep a=" u 0:4 via az
+fit bx "<grep -h -A3 '<unit_cell' $* | grep b=" u 0:2 via bx
+fit by "<grep -h -A3 '<unit_cell' $* | grep b=" u 0:3 via by
+fit bz "<grep -h -A3 '<unit_cell' $* | grep b=" u 0:4 via bz
+fit cx "<grep -h -A3 '<unit_cell' $* | grep c=" u 0:2 via cx
+fit cy "<grep -h -A3 '<unit_cell' $* | grep c=" u 0:3 via cy
+fit cz "<grep -h -A3 '<unit_cell' $* | grep c=" u 0:4 via cz
+print "avg a: ", ax, ay, az
+print "avg b: ", bx, by, bz
+print "avg c: ", cx, cy, cz
+EOF

util/econst.plt

+#!/bin/bash
+# econst.plt: plot <econst> in one or more MD simulations
+# use: econst.plt mdrun1.r [mdrun2.r ...]
+gnuplot -persist <<EOF
+fit a*x+b "<grep -h econst $*" u 0:2 via a,b
+fit c "<grep -h econst $*" u 0:2 via c
+p "<grep -h econst $*" u 2 w l, a*x+b, c
+print "Econst_avg=",c,"   dE/dt=",a
+EOF

util/econst_cmp.plt

+#!/bin/bash
+# econst_cmp.plt: compare <econst> in two MD simulations
+# use: econst_cmp.plt mdrun1.r mdrun2.r
+gnuplot -persist <<EOF
+plot "<grep '<econst>' $1" u 2 w l, "<grep '<econst>' $2" u 2 w l
+EOF

util/econste.plt

+#!/bin/bash
+# econste.plt: plot <econst> and <etotal> in one or more MD simulations
+# use: econste.plt mdrun1.r [mdrun2.r ...]
+gnuplot -persist <<EOF
+p "<grep -h econst $*" u 2 w l, "<grep -h '<etotal>' $*" u 2 w l
+EOF

util/econste_cmp.plt

+#!/bin/bash
+# econste_cmp.plt: compare <econst> and <etotal> in two MD simulations
+# use: econste_cmp.plt mdrun1.r mdrun2.r
+gnuplot -persist <<EOF
+p "<grep -h econst $1" u 2 w l, "<grep -h '<etotal>' $1" u 2 w l, \
+  "<grep -h econst $2" u 2 w l, "<grep -h '<etotal>' $2" u 2 w l
+
+EOF

util/egap.sh

+#!/bin/bash
+# compute Egap = E(n+1) - E(n)
+# use: egap.sh n file.r
+declare -i nocc=$1
+shift
+declare -i nl=nocc/5+1
+#echo "nl=" $nl
+declare -i nfrac=nocc-5*\(nocc/5\)
+#echo "nfrac=" $nfrac
+
+grep -h -A$nl '<eigenvalues ' ${*}| \
+ awk -v nl=$nl -v nfrac=$nfrac \
+  ' NR%(nl+2)==nl {e[0] = $5;} \
+    NR%(nl+2)==(nl+1) \
+    { e[1]=$1; e[2]=$2; e[3]=$3; e[4]=$4; e[5]=$5; \
+      e_homo = e[nfrac]; e_lumo = e[nfrac+1]; \
+      print "E_HOMO=",e_homo, "E_LUMO=",e_lumo, "Eg=",e_lumo-e_homo}' -

util/etotal.plt

+#!/bin/bash
+# etotal.plt: plot <etotal> in one or more simulations
+# use: etotal.plt run1.r [run2.r ...]
+gnuplot -persist <<EOF
+plot "<grep -h '<etotal>' $*" u 2 w l
+EOF

util/etotal_cmp.plt

+#!/bin/bash
+# etotal_cmp.plt: compare <etotal> in two MD simulations
+# use: etotal_cmp.plt mdrun1.r mdrun2.r
+gnuplot -persist <<EOF
+set grid
+plot "<grep '<etotal>' $1" u 2 w l, "<grep '<etotal>' $2" u 2 w l
+EOF

util/etotal_int.plt

+#!/bin/bash
+# etotal_int.plt: plot <etotal_int> in one or more simulations
+# use: etotal_int.plt run1.r  [run2.r ...]
+gnuplot -persist <<EOF
+plot "<grep -h etotal $*" u 2 w l
+EOF

util/etotal_int_cmp.plt

+#!/bin/bash
+# etotal_int_cmp.plt: compare <etotal_int> in two simulations
+# use: etotal_int_cmp.plt run1.r run2.r
+gnuplot -persist <<EOF
+plot "<grep '<etotal_int>' $1" u 2 w l, "<grep '<etotal_int>' $2" u 2 w l
+EOF

util/force.plt

+#!/bin/bash
+# force.plt: plot all components of ionic forces in one or more simulations
+# use: force.plt mdrun1.r [mdrun2.r ...]
+gnuplot -persist <<EOF
+p "<grep -h force $*" u 2, "" u 3, "" u 4, 0
+EOF

util/get_atomset.sh

+#!/bin/bash
+# get_atomset: extract atomset from a Qbox sample
+#
+# use: get_atomset sample.xml
+#
+nlines=$(grep /atomset -m 1 -n $1 | cut -f1 -d: - )
+head -$nlines $1 
+echo "</fpmd:sample>"

util/qbox_dos.py

+#!/usr/bin/python
+# qbox_dos.py: extract electronic DOS from Qbox output
+# generate DOS plot in gnuplot format
+# use: qbox_dos.py [-last] emin emax width file.r
+# emin, emax: bounds of plot in [eV]
+# width: gaussian broadening in [eV]
+# the DOS is accumulated separately for each spin
+# With the -last option, only the last <eigenset> is used to compute the DOS
+
+import xml.sax
+import sys
+import math
+
+if (len(sys.argv) != 5) and (len(sys.argv) != 6) :
+  print "use: ",sys.argv[0]," [-last] emin emax width file.r"
+  sys.exit()
+
+iarg = 1
+lastonly = False
+if (sys.argv[iarg] == "-last") :
+  lastonly = True
+  iarg += 1
+
+emin = float(sys.argv[iarg])
+iarg += 1
+emax = float(sys.argv[iarg])
+iarg += 1
+width = float(sys.argv[iarg])
+iarg += 1
+infile = sys.argv[iarg]
+
+ndos = 501
+de = (emax - emin)/(ndos-1)
+
+# normalized gaussian distribution in one dimension
+# f(x) = 1/(sqrt(pi)*width) * exp(-(x/width)^2 )
+def gauss(x, width):
+  return (1.0/(math.sqrt(math.pi)*width)) * math.exp(-(x/width)**2)
+  
+# Qbox output handler to extract and process data
+class QboxOutputHandler(xml.sax.handler.ContentHandler):
+  def __init__(self):
+    self.nspin = 1
+    self.readData = 0
+    self.dos_up = [0] * ndos
+    self.dos_dn = [0] * ndos
+
+  def startElement(self, name, attributes):
+    if (name == "eigenset") and (lastonly):
+      self.dos_up = [0] * ndos
+      self.dos_dn = [0] * ndos
+    if name == "eigenvalues":
+      self.n = attributes["n"]
+      self.spin = int(attributes["spin"])
+      self.kpoint = attributes["kpoint"]
+      self.weight = float(attributes["weight"])
+      self.readData = 1
+      self.buffer = ""
+      if self.spin == 1:
+        self.nspin = 2
+
+  def characters(self, data):
+    if self.readData:
+      self.buffer += data
+
+  def endElement(self, name):
+    if name == "eigenvalues":
+      self.readData = 0
+      self.accumulate_dos()
+
+  def accumulate_dos(self):
+    self.e = self.buffer.split()
+    if self.spin == 0:
+      for i in range(len(self.e)):
+        for j in range(ndos):
+          ej = emin + j * de
+          self.dos_up[j] += gauss(float(self.e[i])-ej, width ) * self.weight
+    if self.spin == 1:
+      for i in range(len(self.e)):
+        for j in range(ndos):
+          ej = emin + j * de
+          self.dos_dn[j] += gauss(float(self.e[i])-ej, width ) * self.weight
+
+  def print_dos(self):
+    print "# ",infile," spin=0 width=",width
+    for j in range(ndos):
+      ej = emin + j * de
+      print ej, self.dos_up[j]
+    if self.nspin == 2:
+      print
+      print
+      print "# ",infile," spin=1 width=",width
+      for j in range(ndos):
+        ej = emin + j * de
+        print ej, self.dos_dn[j]
+
+parser = xml.sax.make_parser()
+handler = QboxOutputHandler()
+parser.setContentHandler(handler)
+parser.parse(infile)
+handler.print_dos()

util/qbox_maxforce.sh

+#!/bin/bash
+# get the largest force component in a file of "<force> fx fy fz </force>"
+# use: qbox_maxforce nat file.r
+grep '<force>' $2 | tail -$1 | \
+awk '{
+  if ( mx*mx < $2*$2 ) mx = $2;
+  if ( my*my < $3*$3 ) my = $3;
+  if ( mz*mz < $4*$4 ) mz = $4;
+} END {printf("%9.2e %9.2e %9.2e\n", mx, my, mz)}' - 

util/qbox_msd.py

+#!/usr/bin/python
+# qbox_msd.py: compute the mean-square displacement of atoms of a given species
+# generate plot of <r^2(t)> in gnuplot format
+# use: qbox_msd.py species file.r
+
+import xml.sax
+import sys
+import math
+
+if len(sys.argv) != 3:
+  print "use: ",sys.argv[0]," species file.r"
+  sys.exit()
+
+species = sys.argv[1]
+infile = sys.argv[2]
+
+# Qbox output handler to extract and process data
+class QboxOutputHandler(xml.sax.handler.ContentHandler):
+  def __init__(self):
+    self.step = 0
+    self.readData = 0
+    self.readAtom = 0
+    self.tau0 = []
+
+  def startElement(self, name, attributes):
+    if name == "atomset":
+      self.tau=[]
+    elif name == "atom":
+      self.readAtom = (attributes["species"] == species)
+    elif (name == "position") & self.readAtom:
+        self.readData = 1
+        self.buffer = ""
+
+  def characters(self, data):
+    if self.readData:
+      self.buffer += data
+
+  def endElement(self, name):
+    if (name == "position") & self.readAtom:
+      self.readData = 0
+      pos = self.buffer.split()
+      x = float(pos[0])
+      y = float(pos[1])
+      z = float(pos[2])
+      self.tau.append([x,y,z])
+    elif name == "atomset":
+      if self.step == 0:
+        # copy initial positions to tau0
+        self.tau0 = self.tau
+      # compute square displacement
+      disp2sum = 0.0
+      for i in range(len(self.tau)):
+        dx = self.tau[i][0]-self.tau0[i][0]
+        dy = self.tau[i][1]-self.tau0[i][1]
+        dz = self.tau[i][2]-self.tau0[i][2]
+        disp2sum += dx*dx + dy*dy + dz*dz
+      print self.step,'%12.6f'%(disp2sum/len(self.tau))
+      self.step += 1
+
+print "# ",species,infile
+parser = xml.sax.make_parser()
+handler = QboxOutputHandler()
+parser.setContentHandler(handler)
+parser.parse(infile)
+print
+print

util/qbox_repair_h2o.py

+#!/usr/bin/python
+# qbox_repair_h2o.py: repair broken h2o molecules in a Qbox sys file
+# move hydrogen atoms across periodic boundaries to repair molecules
+# use: qbox_repair_h2o.py file.sys
+
+import sys
+import math
+
+olist = []
+hlist = []
+
+def distance(a,b,sx,sy,sz):
+  return math.sqrt((a[3]-b[3]-sx)**2+(a[4]-b[4]-sy)**2+(a[5]-b[5]-sz)**2)
+
+def fold_in_ws(atom):
+  x = atom[3]
+  y = atom[4]
+  z = atom[5]
+  while x > 0.5*a_cell + 1.e-5:
+    x -= a_cell
+  while x < -0.5*a_cell - 1.e-5:
+    x += a_cell
+  while y > 0.5*b_cell + 1.e-5:
+    y -= b_cell
+  while y < -0.5*b_cell - 1.e-5:
+    y += b_cell
+  while z > 0.5*c_cell + 1.e-5:
+    z -= c_cell
+  while z < -0.5*c_cell - 1.e-5:
+    z += c_cell
+  atom[3] = x
+  atom[4] = y
+  atom[5] = z
+
+f = open(sys.argv[1])
+for line in f:
+  l = line.split()
+  if ( l[0] == "set" ) & ( l[1] == "cell" ):
+    a_cell = float(l[2])
+    b_cell = float(l[6])
+    c_cell = float(l[10])
+    print line
+  elif ( l[0] == "species" ):
+    print line
+  elif ( l[0] == "atom" ) & ( l[2] == "oxygen" ):
+    olist.append([l[0],l[1],l[2],float(l[3]),float(l[4]),float(l[5])])
+  elif ( l[0] == "atom" ) & ( l[2] == "hydrogen" ):
+    hlist.append([l[0],l[1],l[2],float(l[3]),float(l[4]),float(l[5])])
+
+for o in olist:
+  fold_in_ws(o)
+for h in hlist:
+  fold_in_ws(h)
+
+for h in hlist:
+  # find nearest oxygen atom in olist
+  mindist = 1.e10
+  for o in olist:
+    # compute minimal distance o-h
+    for sx in [-a_cell,0,a_cell]:
+      for sy in [-b_cell,0,b_cell]:
+        for sz in [-c_cell,0,c_cell]:
+          d = distance(o,h,sx,sy,sz)
+          # print "dist(",sx,sy,sz,") = ",d
+          if d < mindist:
+            mindist = d
+            sx_min = sx
+            sy_min = sy
+            sz_min = sz
+            nearest_o = o;
+  # print "min shift is: ",sx_min,sy_min,sz_min
+  if ( sx_min != 0 ) | ( sy_min != 0 ) | ( sz_min != 0 ):
+    print "# current ",h[1]," at ", h[3],h[4],h[5]
+    print "# nearest O is at ", nearest_o[3],nearest_o[4],nearest_o[5]
+    print "# move ",h[1]," by ", sx_min, sy_min, sz_min
+    h[3] += sx_min    
+    h[4] += sy_min    
+    h[5] += sz_min    
+
+for o in olist:
+  print o[0],o[1],o[2],'%10.5f'%o[3],'%10.5f'%o[4],'%10.5f'%o[5]
+for h in hlist:
+  print h[0],h[1],h[2],'%10.5f'%h[3],'%10.5f'%h[4],'%10.5f'%h[5]

util/qbox_replicate.sh

+#!/bin/bash
+#
+# qbox_replicate: replicate the unit cell in the a0,a1,a2 directions.
+# use: qbox_replicate cell.sys n0 n1 n2 > newcell.sys
+#
+if (( $# != 4 ))
+  then echo "use: qbox_replicate cell.sys n0 n1 n2 > newcell.sys"
+  exit
+fi
+gawk -v n0=$2 -v n1=$3 -v n2=$4 \
+    ' / cell/  {a0x=$3;a0y=$4;a0z=$5; \
+               a1x=$6;a1y=$7;a1z=$8; \
+               a2x=$9;a2y=$10;a2z=$11; \
+      print "set cell ", \
+        n0*$3,n0*$4,n0*$5, n1*$6,n1*$7,n1*$8, n2*$9,n2*$10,n2*$11} \
+      /ref_cell/  { \
+      print "set ref_cell ", \
+        n0*$3,n0*$4,n0*$5, n1*$6,n1*$7,n1*$8, n2*$9,n2*$10,n2*$11} \
+      /species/  {print} \
+      /atom/ {x=$4 - (n0-1)*a0x/2 - (n1-1)*a1x/2 - (n2-1)*a2x/2; \
+              y=$5 - (n0-1)*a0y/2 - (n1-1)*a1y/2 - (n2-1)*a2y/2; \
+              z=$6 - (n0-1)*a0z/2 - (n1-1)*a1z/2 - (n2-1)*a2z/2; \
+              for ( i=0; i<n0; i++ )
+              for ( j=0; j<n1; j++ )
+              for ( k=0; k<n2; k++ )
+                printf("atom %s_%d%d%d %s %12.6f %12.6f %12.6f\n", \
+                       $2,i,j,k,$3, \
+                  x+i*a0x+j*a1x+k*a2x, \
+                  y+i*a0y+j*a1y+k*a2y, \
+                  z+i*a0z+j*a1z+k*a2z) \
+             }' $1

util/qbox_translate.sh

+#!/bin/bash
+# qbox_translate: translate all atoms 
+# use: qbox_replicate cell.sys dx dy dz > newcell.sys
+#
+if (( $# != 4 ))
+  then echo "use: qbox_translate cell.sys dx dy dz > newcell.sys"
+  exit
+fi
+awk -v dx=$2 -v dy=$3 -v dz=$4 \
+    ' /^#/  {print} \
+      / cell/  {print} \
+      /ref_cell/  {print}
+      /species/  {print} \
+      /atom/ {x=$4 + dx; y=$5 + dy; z=$6 + dz; \
+              printf("atom %s %s %12.6f %12.6f %12.6f\n", \
+                       $2,$3,x,y,z) \
+             }' $1

util/qbox_xyz.sh

+#!/bin/bash
+# qbox_xyz.sh: get atomic positions from an MD simulation in xyz format
+# use: qbox_xyz.sh  mdrun.r > file.xyz
+#
+xsltproc - $1 << EOF
+<?xml version="1.0"?>
+<xsl:stylesheet version="1.0" xmlns:xsl="http://www.w3.org/1999/XSL/Transform"
+xmlns:fpmd="http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0">
+<xsl:output method="text" indent="yes"/>
+<xsl:strip-space elements="*"/>
+<xsl:template match="/">
+  <xsl:apply-templates select="//iteration/atomset"/>
+</xsl:template>
+
+<xsl:template match="iteration/atomset">
+  <xsl:value-of select="count(child::atom)"/> <xsl:text>
+</xsl:text>
+<xsl:number count="iteration"/> <xsl:text>
+</xsl:text>
+  <xsl:apply-templates select="atom"/>
+</xsl:template>
+
+<xsl:template match="atom">
+  <xsl:variable name="sym" select="substring(@name,1,2)"/>
+  <xsl:variable name="symbol" select="translate(\$sym,'0123456789_-:.',' ')"/>
+  <xsl:value-of select="\$symbol"/> <xsl:text> </xsl:text>
+  <xsl:variable name="pos" select="normalize-space(position)"/>
+  <xsl:variable name="x" select="substring-before(\$pos,' ')"/>
+  <xsl:variable name="y" select="substring-before(substring-after(\$pos,' '),' ')"/>
+  <xsl:variable name="z" select="substring-after(substring-after(\$pos,' '),' ')"/>
+  <xsl:value-of select="\$x * 0.529177"/> <xsl:text> </xsl:text>
+  <xsl:value-of select="\$y * 0.529177"/> <xsl:text> </xsl:text>
+  <xsl:value-of select="\$z * 0.529177"/> <xsl:text>
+</xsl:text>
+</xsl:template>
+
+<xsl:template match="*"/>
+</xsl:stylesheet>
+EOF

util/sample_to_move.sh

+#!/bin/bash
+# sample_to_move.sh
+# generate a set of Qbox commands to move atoms to the positions
+# in a given sample.
+# use: sample_to_move.sh sample.xml
+#
+xsltproc - $1 << EOF
+<?xml version="1.0"?>
+<xsl:stylesheet version="1.0" xmlns:xsl="http://www.w3.org/1999/XSL/Transform"
+xmlns:fpmd="http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0">
+<xsl:output method="text" indent="yes"/>
+<xsl:strip-space elements="*"/>
+<xsl:template match="/fpmd:sample">
+  <xsl:apply-templates/>
+</xsl:template>
+<xsl:template match="atomset">
+  <xsl:apply-templates/>
+</xsl:template>
+<xsl:template match="atom">
+  <xsl:text>move </xsl:text>
+  <xsl:value-of select="@name"/> 
+  <xsl:text> to </xsl:text>
+  <xsl:value-of select="position"/> <xsl:text>
+</xsl:text>
+</xsl:template>
+<xsl:template match="*"/>
+</xsl:stylesheet>
+EOF

util/sample_to_sys.sh

+#!/bin/bash
+# sample_to_sys.sh: extract atomic positions from a sample and generate
+# a Qbox input file to create atoms at these positions
+# use: sample_to_sys.sh sample.xml > sample.sys
+#
+xsltproc - $1 << EOF
+<?xml version="1.0"?>
+<xsl:stylesheet version="1.0" xmlns:xsl="http://www.w3.org/1999/XSL/Transform"
+xmlns:fpmd="http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0">
+<xsl:output method="text" indent="yes"/>
+<xsl:strip-space elements="*"/>
+<xsl:template match="/fpmd:sample">
+  <xsl:apply-templates/>
+</xsl:template>
+<xsl:template match="atomset">
+  <xsl:apply-templates/>
+</xsl:template>
+<xsl:template match="unit_cell">
+  <xsl:text>set cell </xsl:text>
+  <xsl:value-of select="@a"/> <xsl:text> </xsl:text>
+  <xsl:value-of select="@b"/> <xsl:text> </xsl:text>
+  <xsl:value-of select="@c"/>  <xsl:text>
+</xsl:text>
+</xsl:template>
+<xsl:template match="species">
+  <xsl:text>species </xsl:text>
+  <xsl:value-of select="@name"/> <xsl:text> </xsl:text>
+  <xsl:value-of select="@href"/> <xsl:text>
+</xsl:text>
+</xsl:template>
+<xsl:template match="atom">
+  <xsl:text>atom </xsl:text>
+  <xsl:value-of select="@name"/> <xsl:text> </xsl:text>
+  <xsl:value-of select="@species"/> <xsl:text> </xsl:text>
+  <xsl:value-of select="position"/>  <xsl:text> </xsl:text> 
+  <xsl:value-of select="velocity"/> <xsl:text> 
+</xsl:text>
+</xsl:template>
+<xsl:template match="*"/>
+</xsl:stylesheet>
+EOF

util/sigma.plt

+#!/bin/bash
+# sigma.plt: plot stress tensor components in one or more simulations
+# use: sigma.plt mdrun1.r [mdrun2.r ...]
+gnuplot -persist <<EOF
+plot "<grep -h sigma_xx $*" u 2 w l, "<grep -h sigma_yy $*" u 2 w l, \
+"<grep -h sigma_zz $*" u 2 w l, "<grep -h sigma_xy $*" u 2 w l, \
+"<grep -h sigma_yz $*" u 2 w l, "<grep -h sigma_xz $*" u 2 w l, 0
+fit sxx "<grep -h '<sigma_xx' $*" u 0:2 via sxx
+fit syy "<grep -h '<sigma_yy' $*" u 0:2 via syy
+fit szz "<grep -h '<sigma_zz' $*" u 0:2 via szz
+fit sxy "<grep -h '<sigma_xy' $*" u 0:2 via sxy
+fit syz "<grep -h '<sigma_yz' $*" u 0:2 via syz
+fit sxz "<grep -h '<sigma_xz' $*" u 0:2 via sxz
+print "avg sxx: ", sxx
+print "avg syy: ", syy
+print "avg szz: ", szz
+print "avg sxy: ", sxy
+print "avg syz: ", syz
+print "avg sxz: ", sxz
+EOF

util/temp_ion.plt

+#!/bin/bash
+# temp_ion.plt: plot <temp_ion> in one or more MD simulations
+# use: temp_ion.plt mdrun1.r [mdrun2.r ...]
+gnuplot -persist <<EOF
+fit a*x+b "<grep -h temp_ion $*" u 0:2 via a,b
+fit c "<grep -h temp_ion $*" u 0:2 via c
+p "<grep -h temp_ion $*" u 2 w l, a*x+b, c
+ print "Tavg=",c,"   dT/dt=",a
+EOF

util/volume.plt

+#!/bin/bash
+# volume.plt: plot the unit cell volume during one or more MD simulations
+# use: volume.plt mdrun1.r [mdrun2.r ...]
+gnuplot -persist <<EOF
+set grid
+plot "<grep -h '<unit_cell_volume>' $*" u 2 w l
+EOF