merge efield branch into trunk

git-svn-id: http://qboxcode.org/svn/qb/trunk@1712 cba15fb0-1239-40c8-b417-11db7ca47a34

src/AngleCmd.h

@@ -110,7 +110,7 @@ class AngleCmd : public Cmd
     {
       D3vector r12(a1->position()-a2->position());
       D3vector r32(a3->position()-a2->position());
-      if ( norm(r12) == 0.0 || norm(r32) == 0.0 )
+      if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 )
       {
         cout << " AngleCmd: atoms are too close" << endl;
         return 1;

src/AngleConstraint.C

@@ -323,8 +323,8 @@ void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2,
 {
   D3vector r12(r1-r2);
   D3vector r32(r3-r2);
-  assert(norm(r12) > 0.0);
-  assert(norm(r32) > 0.0);
+  assert(norm2(r12) > 0.0);
+  assert(norm2(r32) > 0.0);
   D3vector e12(normalized(r12));
   D3vector e32(normalized(r32));
   const double ss = length(e12^e32);
@@ -344,7 +344,7 @@ void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2,
     {
       // choose direction e as e12+e32
       D3vector e(e12+e32);
-      assert(norm(e)>0.0);
+      assert(norm2(e)>0.0);
       e = normalized(e);
       const double r12_inv = 1.0/length(r12);
       const double r32_inv = 1.0/length(r32);

src/AtomSet.C

@@ -15,7 +15,6 @@
 // AtomSet.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: AtomSet.C,v 1.29 2010-04-16 21:40:50 fgygi Exp $
 
 #include "AtomSet.h"
 #include "Species.h"
@@ -576,6 +575,24 @@ D3vector AtomSet::dipole(void) const
 }
 
 ////////////////////////////////////////////////////////////////////////////////
+D3tensor AtomSet::quadrupole(void) const
+{
+  D3tensor sum;
+  for ( int is = 0; is < atom_list.size(); is++ )
+  {
+    double charge = species_list[is]->zval();
+    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
+    {
+      D3vector p = atom_list[is][ia]->position();
+      for ( int idir = 0; idir < 3; idir++ )
+        for ( int jdir = 0; jdir < 3; jdir++ )
+          sum[idir*3+jdir] += charge * p[idir] * p[jdir];
+    }
+  }
+  return sum;
+}
+
+////////////////////////////////////////////////////////////////////////////////
 void AtomSet::set_cell(const D3vector& a, const D3vector& b, const D3vector& c)
 {
   cell_.set(a,b,c);

src/AtomSet.h

@@ -15,7 +15,6 @@
 // AtomSet.h
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: AtomSet.h,v 1.26 2010-05-10 20:52:54 fgygi Exp $
 
 #ifndef ATOMSET_H
 #define ATOMSET_H
@@ -23,6 +22,7 @@
 #include "Context.h"
 #include "Atom.h"
 #include "UnitCell.h"
+#include "D3tensor.h"
 #include <vector>
 #include <string>
 #include <list>
@@ -87,6 +87,7 @@ class AtomSet
   void randomize_positions(double amplitude);
   D3vector vcm(void) const;
   D3vector dipole(void) const;
+  D3tensor quadrupole(void) const;
   void reset_vcm(void);
   void fold_in_ws(void);
   int size(void) const;

src/BOSampleStepper.C

@@ -35,6 +35,7 @@
 #include "CGCellStepper.h"
 #include "AndersonMixer.h"
 #include "MLWFTransform.h"
+#include "D3tensor.h"
 
 #ifdef USE_APC
 #include "apc.h"
@@ -181,8 +182,6 @@ void BOSampleStepper::step(int niter)
   Wavefunction& wf = s_.wf;
   const int nspin = wf.nspin();
 
-  const UnitCell& cell = wf.cell();
-
   const double dt = s_.ctrl.dt;
 
   const string wf_dyn = s_.ctrl.wf_dyn;
@@ -207,7 +206,6 @@ void BOSampleStepper::step(int niter)
                             cell_dyn != "LOCKED" );
   // GS-only calculation:
   const bool gs_only = !atoms_move && !cell_moves;
-  const bool use_confinement = ( s_.ctrl.ecuts > 0.0 );
 
   Timer tm_iter;
 
@@ -376,35 +374,13 @@ void BOSampleStepper::step(int niter)
       const bool compute_forces = true;
       double energy =
         ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
-      double enthalpy = energy;
+      double enthalpy = ef_.enthalpy();
 
       if ( onpe0 )
       {
-        cout.setf(ios::fixed,ios::floatfield);
-        cout.setf(ios::right,ios::adjustfield);
-        cout << "  <ekin>   " << setprecision(8)
-             << setw(15) << ef_.ekin() << " </ekin>\n";
-        if ( use_confinement )
-          cout << "  <econf>  " << setw(15) << ef_.econf() << " </econf>\n";
-        cout << "  <eps>    " << setw(15) << ef_.eps() << " </eps>\n"
-             << "  <enl>    " << setw(15) << ef_.enl() << " </enl>\n"
-             << "  <ecoul>  " << setw(15) << ef_.ecoul() << " </ecoul>\n"
-             << "  <exc>    " << setw(15) << ef_.exc() << " </exc>\n"
-             << "  <esr>    " << setw(15) << ef_.esr() << " </esr>\n"
-             << "  <eself>  " << setw(15) << ef_.eself() << " </eself>\n"
-             << "  <ets>    " << setw(15) << ef_.ets() << " </ets>\n";
-        if ( s_.extforces.size() > 0 )
-          cout << "  <eexf>     " << setw(15) << ef_.eexf() << " </eexf>\n";
-        cout << "  <etotal> " << setw(15) << ef_.etotal() << " </etotal>\n";
-        if ( compute_stress )
-        {
-          const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
-          enthalpy = ef_.etotal() + pext * cell.volume();
-          cout << "  <pv>     " << setw(15) << pext * cell.volume()
-               << " </pv>" << endl;
-          cout << "  <enthalpy> " << setw(15) << enthalpy << " </enthalpy>\n"
-             << flush;
-        }
+        cout << ef_;
+        if ( ef_.el_enth() )
+          cout << *ef_.el_enth();
       }
 
       if ( iter > 0 && ionic_stepper )
@@ -468,7 +444,8 @@ void BOSampleStepper::step(int niter)
         if ( ionic_stepper != 0 )
           ekin_stepper = ionic_stepper->ekin_stepper();
         cout << setprecision(8);
-        cout << "  <econst> " << energy+ekin_ion+ekin_stepper << " </econst>\n";
+        cout << "  <econst> " << enthalpy+ekin_ion+ekin_stepper
+             << " </econst>\n";
         cout << "  <ekin_ion> " << ekin_ion << " </ekin_ion>\n";
         cout << "  <temp_ion> " << temp_ion << " </temp_ion>\n";
       }
@@ -865,7 +842,7 @@ void BOSampleStepper::step(int niter)
         while ( !nonscf_converged && ite < max(nite_,1) )
         {
           energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
-          double enthalpy = energy;
+          double enthalpy = ef_.enthalpy();
 
           if ( ite > 0 )
             etotal_int_m = etotal_int;
@@ -884,7 +861,7 @@ void BOSampleStepper::step(int niter)
 
           eigenvalue_sum = real(s_.wf.dot(dwf));
           if ( onpe0 )
-            cout << "  <eigenvalue_sum> "
+            cout << "  <eigenvalue_sum>  "
                  << eigenvalue_sum << " </eigenvalue_sum>" << endl;
 
           wf_stepper->update(dwf);
@@ -893,16 +870,13 @@ void BOSampleStepper::step(int niter)
           {
             cout.setf(ios::fixed,ios::floatfield);
             cout.setf(ios::right,ios::adjustfield);
-            cout << "  <etotal_int> " << setprecision(8) << setw(15)
+            cout << "  <etotal_int>  " << setprecision(8) << setw(15)
                  << energy << " </etotal_int>\n";
-          }
-          if ( compute_stress )
-          {
-            const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
-            enthalpy = energy + pext * cell.volume();
-            if ( onpe0 )
-              cout << "  <enthalpy_int> " << setw(15)
+            if ( compute_stress || ef_.el_enth() )
+            {
+              cout << "  <enthalpy_int>" << setw(15)
                    << enthalpy << " </enthalpy_int>\n" << flush;
+            }
           }
 
           // compare delta_etotal_int only after first iteration
@@ -1007,6 +981,7 @@ void BOSampleStepper::step(int niter)
       {
         for ( int ispin = 0; ispin < nspin; ispin++ )
         {
+          mlwft[ispin]->update();
           mlwft[ispin]->compute_transform();
         }
 
@@ -1028,20 +1003,36 @@ void BOSampleStepper::step(int niter)
             SlaterDet& sd = *(wf.sd(ispin,0));
             cout << " <mlwfset spin=\"" << ispin
                  << "\" size=\"" << sd.nst() << "\">" << endl;
+            double total_spread[6];
+            for ( int j = 0; j < 6; j++ )
+               total_spread[j] = 0.0;
             for ( int i = 0; i < sd.nst(); i++ )
             {
               D3vector ctr = mlwft[ispin]->center(i);
-              double sp = mlwft[ispin]->spread(i);
+              double spi[6];
+              for (int j=0; j<3; j++)
+              {
+                spi[j] = mlwft[ispin]->spread2(i,j);
+                total_spread[j] += spi[j];
+              }
+
               cout.setf(ios::fixed, ios::floatfield);
               cout.setf(ios::right, ios::adjustfield);
               cout << "   <mlwf center=\"" << setprecision(6)
                    << setw(12) << ctr.x
                    << setw(12) << ctr.y
                    << setw(12) << ctr.z
-                   << " \" spread=\" " << sp << " \"/>"
+                   << " \" spread=\" "
+                   << setw(12) << spi[0]
+                   << setw(12) << spi[1]
+                   << setw(12) << spi[2] << " \"/>"
                    << endl;
             }
 
+            cout << " <total_spread> ";
+            for ( int j = 0; j < 3; j++ )
+              cout << setw(10) << total_spread[j];
+            cout << " </total_spread>" << endl;
             D3vector edipole = mlwft[ispin]->dipole();
             cout << " <electronic_dipole spin=\"" << ispin << "\"> " << edipole
                  << " </electronic_dipole>" << endl;
@@ -1073,6 +1064,8 @@ void BOSampleStepper::step(int niter)
         if ( onpe0 )
         {
           cout << ef_;
+          if ( ef_.el_enth() )
+            cout << *ef_.el_enth();
           cout << "<atomset>" << endl;
           cout << atoms.cell();
           for ( int is = 0; is < atoms.atom_list.size(); is++ )
@@ -1132,6 +1125,9 @@ void BOSampleStepper::step(int niter)
       if ( onpe0 )
       {
         cout << ef_;
+        if ( ef_.el_enth() )
+          cout << *ef_.el_enth();
+
       }
     }
 

src/Basis.C

@@ -307,7 +307,7 @@ bool Basis::resize(const UnitCell& cell, const UnitCell& refcell,
   // defcell: cell used to define which vectors are contained in the Basis
   // if refcell is defined, defcell = refcell
   // otherwise, defcell = cell
-  if ( norm(refcell.b(0)) + norm(refcell.b(1)) + norm(refcell.b(2)) == 0.0 )
+  if ( norm2(refcell.b(0)) + norm2(refcell.b(1)) + norm2(refcell.b(2)) == 0.0 )
   {
     defcell = cell;
   }
@@ -320,7 +320,7 @@ bool Basis::resize(const UnitCell& cell, const UnitCell& refcell,
   const D3vector b1 = defcell.b(1);
   const D3vector b2 = defcell.b(2);
 
-  const double normb2 = norm(b2);
+  const double normb2 = norm2(b2);
   const double b2inv2 = 1.0 / normb2;
 
   const D3vector kp = kpx*b0 + kpy*b1 + kpz*b2;
@@ -380,7 +380,7 @@ bool Basis::resize(const UnitCell& cell, const UnitCell& refcell,
       bool found = false;
       for ( int l = lmin; l <= lmax; l++ )
       {
-        const double two_e = norm(k*b1+l*b2);
+        const double two_e = norm2(k*b1+l*b2);
         if ( two_e < two_ecut )
         {
           lstart = min(l,lstart);
@@ -410,7 +410,7 @@ bool Basis::resize(const UnitCell& cell, const UnitCell& refcell,
         bool found = false;
         for ( int l = lmin; l <= lmax; l++ )
         {
-          const double two_e = norm(h*b0+k*b1+l*b2);
+          const double two_e = norm2(h*b0+k*b1+l*b2);
           if ( two_e < two_ecut )
           {
             lstart = min(l,lstart);
@@ -450,7 +450,7 @@ bool Basis::resize(const UnitCell& cell, const UnitCell& refcell,
         bool found = false;
         for ( int l = lmin; l <= lmax; l++ )
         {
-          const double two_e = norm((kpx+h)*b0 + (kpy+k)*b1 + (kpz+l)*b2);
+          const double two_e = norm2((kpx+h)*b0 + (kpy+k)*b1 + (kpz+l)*b2);
           if ( two_e < two_ecut )
           {
             lstart = min(l,lstart);
@@ -658,10 +658,10 @@ void Basis::update_g(void)
     kpgx_[locsize+i] = kpgt.y;
     kpgx_[locsize+locsize+i] = kpgt.z;
 
-    g2_[i] = norm(gt);
+    g2_[i] = norm2(gt);
     g_[i] = sqrt( g2_[i] );
 
-    kpg2_[i] = norm(kpgt);
+    kpg2_[i] = norm2(kpgt);
     kpg_[i] = sqrt( kpg2_[i] );
 
     gi_[i] = g_[i] > 0.0 ? 1.0 / g_[i] : 0.0;

src/CPSampleStepper.C

@@ -109,7 +109,6 @@ void CPSampleStepper::step(int niter)
   const bool compute_hpsi = true;
   const bool compute_forces = ( atoms_dyn != "LOCKED" );
   const bool compute_stress = ( s_.ctrl.stress == "ON" );
-  const bool use_confinement = ( s_.ctrl.ecuts > 0.0 );
 
   CellStepper* cell_stepper = 0;
   if ( cell_dyn == "SD" )
@@ -135,6 +134,7 @@ void CPSampleStepper::step(int niter)
   ef_.update_vhxc(compute_stress);
   double energy =
     ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
+  double enthalpy = ef_.enthalpy();
 
   mdwf_stepper->compute_wfm(dwf);
 
@@ -151,34 +151,9 @@ void CPSampleStepper::step(int niter)
 
     if ( onpe0 )
     {
-      cout.setf(ios::fixed,ios::floatfield);
-      cout.setf(ios::right,ios::adjustfield);
-      cout << "  <ekin>     " << setprecision(8)
-           << setw(15) << ef_.ekin() << " </ekin>\n";
-      if ( use_confinement )
-      {
-        cout << "  <econf>    " << setw(15) << ef_.econf()
-             << " </econf>\n";
-      }
-      cout << "  <eps>      " << setw(15) << ef_.eps() << " </eps>\n"
-           << "  <enl>      " << setw(15) << ef_.enl() << " </enl>\n"
-           << "  <ecoul>    " << setw(15) << ef_.ecoul() << " </ecoul>\n"
-           << "  <exc>      " << setw(15) << ef_.exc() << " </exc>\n"
-           << "  <esr>      " << setw(15) << ef_.esr() << " </esr>\n"
-           << "  <eself>    " << setw(15) << ef_.eself() << " </eself>\n";
-      if ( s_.extforces.size() > 0 )
-        cout << "  <eexf>     " << setw(15) << ef_.eexf() << " </eexf>\n";
-      cout << "  <etotal>   " << setw(15) << ef_.etotal() << " </etotal>\n"
-           << flush;
-      if ( compute_stress )
-      {
-        const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
-        const double enthalpy = ef_.etotal() + pext * s_.wf.cell().volume();
-        cout << "  <pv>     " << setw(15) << pext * s_.wf.cell().volume()
-             << " </pv>" << endl;
-        cout << "  <enthalpy> " << setw(15) << enthalpy << " </enthalpy>\n"
-           << flush;
-      }
+      cout << ef_;
+      if ( ef_.el_enth() )
+        cout << *ef_.el_enth();
     }
 
     if ( compute_forces )
@@ -279,6 +254,7 @@ void CPSampleStepper::step(int niter)
     ef_.update_vhxc(compute_stress);
     energy =
       ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
+    enthalpy = ef_.enthalpy();
 
     if ( s_.ctxt_.mype() == 0 )
       cout << "</iteration>" << endl;

src/ComputeMLWFCmd.C

@@ -15,7 +15,6 @@
 // ComputeMLWFCmd.C:
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: ComputeMLWFCmd.C,v 1.8 2008-11-14 22:11:30 fgygi Exp $
 
 #include "ComputeMLWFCmd.h"
 #include<iostream>
@@ -48,6 +47,7 @@ int ComputeMLWFCmd::action(int argc, char **argv)
     SlaterDet& sd = *(wf.sd(ispin,0));
     MLWFTransform* mlwft = new MLWFTransform(sd);
 
+    mlwft->update();
     mlwft->compute_transform();
     mlwft->apply_transform(sd);
 

src/Control.h

@@ -15,13 +15,13 @@
 // Control.h:
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: Control.h,v 1.15 2009-03-08 01:11:31 fgygi Exp $
 
 #ifndef CONTROL_H
 #define CONTROL_H
 
 #include <string>
 #include <vector>
+#include "D3vector.h"
 
 struct Control
 {
@@ -69,5 +69,8 @@ struct Control
   double btHF;
 
   double scf_tol;
+
+  D3vector e_field;
+  std::string polarization;
 };
 #endif

src/D3tensor.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2014 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// D3tensor.h
+//
+// double 3x3 tensor
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef D3TENSOR_H
+#define D3TENSOR_H
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <cassert>
+#include "blas.h"
+#include "D3vector.h"
+
+using namespace std;
+
+class D3tensor
+{
+  public:
+
+  double a_[9];
+
+  double* a(void) { return &a_[0]; }
+  const double* a(void) const { return &a_[0]; }
+
+  explicit D3tensor(void) { clear(); }
+
+  explicit D3tensor(double xx, double yy, double zz)
+  { a_[0]=xx; a_[4]=yy; a_[8]=zz; }
+
+  explicit D3tensor(double xx, double yy, double zz,
+                    double xy, double yz, double xz, char& uplo)
+  {
+    a_[0] = xx;
+    a_[4] = yy;
+    a_[8] = zz;
+
+    if ( uplo == 'l' )
+    {
+      a_[1] = xy;
+      a_[2] = xz;
+      a_[5] = yz;
+    }
+    else if ( uplo == 'u' )
+    {
+      a_[3] = xy;
+      a_[6] = xz;
+      a_[7] = yz;
+    }
+    else if ( uplo == 's' )
+    {
+      a_[1] = xy;
+      a_[2] = xz;
+      a_[5] = yz;
+      a_[3] = xy;
+      a_[6] = xz;
+      a_[7] = yz;
+    }
+    else
+      assert(false);
+  }
+
+  explicit D3tensor(const D3vector& diag, const D3vector& offdiag)
+  {
+    a_[0] = diag[0];
+    a_[4] = diag[1];
+    a_[8] = diag[2];
+    a_[1] = offdiag[0];
+    a_[5] = offdiag[1];
+    a_[2] = offdiag[2];
+    a_[3] = offdiag[0];
+    a_[7] = offdiag[1];
+    a_[6] = offdiag[2];
+  }
+
+  explicit D3tensor(const double* a)
+  {
+    for ( int i = 0; i < 9; i++ ) a_[i] = a[i];
+  }
+
+  double& operator[](int i)
+  {
+    assert(i>=0 && i < 9);
+    return a_[i];
+  }
+
+  double operator[](int i) const
+  {
+    assert(i>=0 && i < 9);
+    return a_[i];
+  }
+
+  void setdiag(int i, double b)
+  {
+    assert(i>=0 && i<3);
+    a_[i*4] = b;
+  }
+
+  void setdiag(const D3vector& b)
+  {
+    for ( int i = 0; i < 3; i++ )
+      a_[i*4] = b[i];
+  }
+
+  void setoffdiag(int i, double b)
+  {
+    assert(i>=0 && i<3);
+    if ( i == 0 )
+    {
+      a_[1] = b;
+      a_[3] = b;
+    }
+    else if ( i == 2 )
+    {
+      a_[2] = b;
+      a_[6] = b;
+    }
+    else
+    {
+      a_[5] = b;
+      a_[7] = b;
+    }
+  }
+
+  void setoffdiag(const D3vector& b)
+  {
+    a_[1] = b[0];
+    a_[3] = b[0];
+    a_[5] = b[1];
+    a_[7] = b[1];
+    a_[2] = b[2];
+    a_[6] = b[2];
+  }
+
+  bool operator==(const D3tensor &rhs) const
+  {
+    bool eq = true;
+    for ( int i = 0; i < 9; i++ )
+    {
+      if ( rhs[i] != a_[i] )
+      {
+        eq &= false;
+      }
+    }
+    return eq;
+  }
+
+  bool operator!=(const D3tensor &rhs) const
+  {
+    bool neq = false;
+    for ( int i = 0; i < 9; i++ )
+    {
+      if ( rhs[i] != a_[i] )
+      {
+        neq |= true;
+      }
+    }
+    return neq;
+  }
+
+  D3tensor& operator+=(const D3tensor& rhs)
+  {
+    for ( int i = 0; i < 9; i++ )
+      a_[i] += rhs[i];
+    return *this;
+  }
+
+  D3tensor& operator-=(const D3tensor& rhs)
+  {
+    for ( int i = 0; i < 9; i++ )
+      a_[i] -= rhs[i];
+    return *this;
+  }
+
+  D3tensor& operator*=(const double& rhs)
+  {
+    for ( int i = 0; i < 9; i++ )
+      a_[i] *= rhs;
+    return *this;
+  }
+
+  D3tensor& operator/=(const double& rhs)
+  {
+    for ( int i = 0; i < 9; i++ )
+      a_[i] /= rhs;
+    return *this;
+  }
+
+  friend const D3tensor operator+(const D3tensor& lhs, const D3tensor& rhs)
+  {
+    return D3tensor(lhs) += rhs;
+  }
+
+  friend const D3tensor operator-(const D3tensor& a, const D3tensor& b)
+  {
+    return D3tensor(a) -= b;
+  }
+
+  friend D3tensor operator*(double a, const D3tensor& b)
+  {
+    return D3tensor(b) *= a;
+  }
+
+  friend D3tensor operator*(const D3tensor& a, double b)
+  {
+    return D3tensor(a) *= b;
+  }
+
+  friend D3tensor operator/(const D3tensor& a, double b)
+  {
+    return D3tensor(a) /= b;
+  }
+
+  friend D3tensor operator*(D3tensor& a, D3tensor& b)
+  {
+    D3tensor c;
+    int ithree = 3;
+    double one = 1.0, zero = 0.0;
+    char t = 'n';
+    dgemm ( &t, &t, &ithree, &ithree, &ithree, &one, &a[0], &ithree,
+            &b[0], &ithree, &zero, &c[0], &ithree );
+    return c;
+  }
+
+  friend D3tensor operator-(const D3tensor& a) // unary minus
+  {
+    return D3tensor()-a;
+  }
+
+  double norm2(const D3tensor& a) const
+  {
+    return a_[0]*a_[0] + a_[1]*a_[1] + a_[2]*a_[2] +
+           a_[3]*a_[3] + a_[4]*a_[4] + a_[5]*a_[5] +
+           a_[6]*a_[6] + a_[7]*a_[7] + a_[8]*a_[8];
+  }
+
+  double norm(const D3tensor& a) const
+  {
+    return sqrt(norm2(a));
+  }
+
+  double trace(void) const
+  {
+    return a_[0]+a_[4]+a_[8];
+  }
+
+  void traceless(void)
+  {
+    double b = trace() / 3;
+    a_[0] -= b;
+    a_[4] -= b;
+    a_[8] -= b;
+  }
+
+  void clear(void)
+  {
+    for ( int i = 0; i < 9; i++ )
+      a_[i] = 0.0;
+  }
+
+  void identity(void)
+  {
+    clear();
+    a_[0] = 1.0;
+    a_[4] = 1.0;
+    a_[8] = 1.0;
+  }
+
+  friend std::ostream& operator<<(std::ostream& os, const D3tensor& rhs)
+  {
+    const double * const v  = rhs.a();
+    os.setf(ios::fixed,ios::floatfield);
+    os.setf(ios::right,ios::adjustfield);
+    os.precision(8);
+    os << setw(14) << v[0] << " " << setw(14) << v[3] << " " << setw(14) << v[6]
+       << "\n"
+       << setw(14) << v[1] << " " << setw(14) << v[4] << " "