initial revision

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

util/dynmat/D3vector.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2008 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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// D3vector.h
+//
+// double 3-vectors
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: D3vector.h,v 1.8 2008/09/08 15:56:18 fgygi Exp $
+
+#ifndef D3VECTOR_H
+#define D3VECTOR_H
+#include <iostream>
+#include <cmath>
+#include <cassert>
+
+class D3vector
+{
+  public:
+
+  double x, y, z;
+
+  // explicit constructor to avoid implicit conversion from double to D3vector
+  explicit D3vector(const double& xv, const double& yv, const double& zv) :
+    x(xv), y(yv), z(zv) {}
+  explicit D3vector(void) : x(0.0), y(0.0), z(0.0) {}
+
+  explicit D3vector(const double* r) : x(r[0]), y(r[1]), z(r[2]) {}
+
+  double& operator[](const int &i)
+  {
+    assert(i>=0 && i <3);
+    if ( i == 0 ) return x;
+    else if ( i == 1 ) return y;
+    else return z;
+  }
+
+  bool operator==(const D3vector &rhs) const
+  {
+    return x == rhs.x && y == rhs.y && z == rhs.z;
+  }
+
+  bool operator!=(const D3vector &rhs) const
+  {
+    return x != rhs.x || y != rhs.y || z != rhs.z;
+  }
+
+  D3vector& operator += ( const D3vector& rhs )
+  {
+    x += rhs.x; y += rhs.y; z += rhs.z;
+    return *this;
+  }
+
+  D3vector& operator -= ( const D3vector& rhs )
+  {
+    x -= rhs.x; y -= rhs.y; z -= rhs.z;
+    return *this;
+  }
+
+  D3vector& operator *= ( const double& rhs )
+  {
+    x *= rhs; y *= rhs; z *= rhs;
+    return *this;
+  }
+
+  D3vector& operator /= ( const double& rhs )
+  {
+    x /= rhs; y /= rhs; z /= rhs;
+    return *this;
+  }
+
+  friend const D3vector operator + (const D3vector& lhs, const D3vector& rhs )
+  {
+    return D3vector(lhs) += rhs;
+  }
+
+  friend const D3vector operator - ( const D3vector& a, const D3vector& b )
+  {
+    return D3vector(a) -= b;
+  }
+
+  friend D3vector operator - ( const D3vector& a ) // unary minus
+  {
+    return D3vector( -a.x, -a.y, -a.z );
+  }
+
+  friend D3vector operator * ( const double& a, const D3vector& b )
+  {
+    return D3vector(b) *= a;
+  }
+
+  friend D3vector operator * ( const D3vector& a, const double& b )
+  {
+    return D3vector(a) *= b;
+  }
+
+  friend D3vector operator / ( const D3vector& a, const double& b )
+  {
+    return D3vector(a) /= b;
+  }
+
+  // scalar product
+  friend double operator * ( const D3vector& a, const D3vector& b )
+  {
+    return a.x * b.x + a.y * b.y + a.z * b.z ;
+  }
+
+  friend D3vector operator ^ ( const D3vector& a, const D3vector& b )
+  {
+    return D3vector( a.y * b.z - a.z * b.y ,
+                     a.z * b.x - a.x * b.z ,
+                     a.x * b.y - a.y * b.x );
+  }
+
+  friend D3vector rotate ( const D3vector& x, const D3vector& w )
+  {
+    if ( length(x) == 0.0 ) return x; // x has zero length
+    double theta = length( w );       // rotate by zero
+    if ( theta == 0.0 ) return x;
+    D3vector ew = normalized ( w );
+    D3vector v = w ^ x;
+    if ( length( v ) == 0.0 ) return x; // x is parallel to the rotation axis
+    v = normalized( v );
+    D3vector u = v ^ ew;
+    double p = x * u;
+    return  (x*ew)*ew + p*cos(theta)*u + p*sin(theta)*v ;
+  }
+
+  friend double length( const D3vector& a )
+  {
+    return sqrt( a.x * a.x + a.y * a.y + a.z * a.z );
+  }
+
+  friend double norm( const D3vector& a )
+  {
+    return a.x * a.x + a.y * a.y + a.z * a.z;
+  }
+
+  friend D3vector normalized( const D3vector a )
+  {
+    return a / length( a );
+  }
+
+  friend std::ostream& operator << ( std::ostream& os, const D3vector& v )
+  {
+    os << v.x << " " << v.y << " " << v.z;
+    return os;
+  }
+
+  friend std::istream& operator >> ( std::istream& is, D3vector& v )
+  {
+    is >> v.x >> v.y >> v.z ;
+    return is;
+  }
+
+};
+#endif

util/dynmat/Makefile

+#-------------------------------------------------------------------------------
+#
+# Copyright (c) 2008 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/>.
+#
+#-------------------------------------------------------------------------------
+# $Id: Makefile,v 1.65 2009/09/08 05:40:50 fgygi Exp $
+#------------------------------------------------------------------------------
+#
+include $(TARGET).mk
+#------------------------------------------------------------------------------
+SRCDIR = $(HOME)/qb/export/rel1_52_1/qb/src/
+EXEC = dynmat
+# use source files in SRCDIR but not objects, which may use MPI
+vpath %.h $(SRCDIR)
+vpath %.C $(SRCDIR)
+OBJECTS = dynmat.o Matrix.o Context.o
+CXXFLAGS += -DTARGET='"$(TARGET)"' -I$(SRCDIR)
+ $(EXEC):    $(OBJECTS)
+	$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
+#------------------------------------------------------------------------------
+# generate dependencies in makefile: use -Y to avoid library header files
+# that are likely to be different on other platforms.
+ depend :
+	makedepend -Y -D$(PLT) *.[cCh]
+#makedepend -D$(PLT) -I/usr/include/CC *.[cC]
+#------------------------------------------------------------------------------
+#  Cleanup object files
+ clean :
+	rm -f *.o
+	rm -f $(EXEC) $(TESTEXECS)
+#------------------------------------------------------------------------------
+ ctags :
+	etags -o tags *.[Ch]
+#------------------------------------------------------------------------------
+# DO NOT DELETE
+
+Context.o: Context.h
+dynmat.o: Context.h Matrix.h
+Matrix.o: Context.h Matrix.h

util/dynmat/README

+
+   Calculation of vibrational frequencies
+   --------------------------------------
+
+1) Generate a Qbox input script "fdmoves.i" using the script "fdmoves.sh"
+
+   Example:
+   $ fdmoves.sh O1 O2 H1 H2 H3 H4 > fdmoves.i
+
+   Note: the atom names must appear in the same order as they are printed in
+   the <atomset> element in Qbox output. This may not correspond to the order
+   in which they were defined with the atom command.
+
+2) Run Qbox using the following input file:
+
+   A sample file describing a relaxed structure must be available (in the
+   following example: "relaxed_sample.xml")
+
+   # finite difference calculation of the dynamical matrix
+   load relaxed_sample.xml
+   set xc PBE
+   set wf_dyn PSDA
+   set ecutprec 4
+   fdmoves.i
+
+   Save the result in "fdmoves.r"
+
+3) Extract all the <force> elements from the file "fdmoves.r"
+
+   $ grep '<force>' fdmoves.r > force.dat
+
+4) Run the dynmat program using "force.dat"
+
+   use: dynmat force.dat N mass1 .. massN
+   where N is the number of atoms, mass1 .. massN are the masses (C=12)
+   Note: the mass must be repeated for each atom.
+
+   $ ./dynmat force.dat 6 16 16 1 1 1 1 
+
+   The frequencies are printed twice in order of increasing magnitude.
+   The first six frequencies should be zero. Deviations from zero give
+   an estimate of the lower bound on the error in frequencies.
+   The dynamical matrix computed from finite-difference force calculations
+   is not exactly symmetric due to errors in the forces and imperfect 
+   relaxation of the original structure. The dynmat program prints two sets
+   of eigenvalues. The first is computed using the upper triangle of the 
+   dynamical matrix and the second using the lower triangle. The difference
+   between the two sets gives an indication of the error due to the lack
+   of symmetry of the dynamical matrix.
+
+   Note: For calculations involving different isotopes, only step 4 must be 
+   repeated with modified mass arguments. The file force.dat can be reused.

util/dynmat/dynmat.C

+//
+// dynmat.C: compute and diagonalize a dynamical matrix
+// Forces are read from Qbox output
+// The Qbox output should correspond to a sequence of calculations
+// using symmetric finite displacements for all atoms in the x,y,z directions
+//
+// use: diag input_file
+// input_file: forces from Qbox XML output file (collected with grep)
+
+#include <cassert>
+#include <cstdlib>
+#include <cmath>
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <fstream>
+#include <valarray>
+#include <vector>
+using namespace std;
+
+#include "Context.h"
+#include "Matrix.h"
+
+int main(int argc, char **argv)
+{
+  if ( argc < 5 )
+  {
+    cout << "use: dynmat force.dat h nat1 mass1 [nat2 mass2] ..." << endl;
+    return 1;
+  }
+  char* infilename = argv[1];
+  ifstream infile(infilename);
+
+  const double h = atof(argv[2]);
+  // cout << "h=" << h << endl;
+  const double Ha2cm1 = 219474.65;
+
+  vector<double> mass;
+  for ( int iarg = 3; iarg < argc; iarg++, iarg++ )
+  {
+    // read number of atoms and mass for each species
+    const int n = atoi(argv[iarg]);
+    const double m = atof(argv[iarg+1]);
+    // cout << n << " atoms of mass " << m << endl;
+    for ( int i = 0; i < n; i++ )
+      mass.push_back(1822.89*m);
+  }
+  const int nat = mass.size();
+
+  // read forces
+  double f[2][nat][3][nat][3];
+  double fx,fy,fz;
+  string dum;
+  for ( int ia1 = 0; ia1 < nat; ia1++ )
+  {
+    for ( int idir1 = 0; idir1 < 3; idir1++ )
+    {
+      for ( int ishift = 0; ishift < 2; ishift++ )
+      {
+        for ( int ia2 = 0; ia2 < nat; ia2++ )
+        {
+          infile >> dum >> fx >> fy >> fz;
+          while (infile.get() != '\n');
+          f[ishift][ia1][idir1][ia2][0] = fx; 
+          f[ishift][ia1][idir1][ia2][1] = fy; 
+          f[ishift][ia1][idir1][ia2][2] = fz; 
+        }
+      }
+    }
+  }
+
+  // compute matrix elements: centered finite difference
+  // a[3*ia1+idir1][3*ia2+idir2] = 
+  //  ( f[1][ia1][idir1][ia2][idir2] - f[0][ia1][idir1][ia2][idir2] ) / ( 2 h )
+
+  Context ctxt;
+  const int n = 3*nat;
+  DoubleMatrix a(ctxt,n,n);
+
+  for ( int ia1 = 0; ia1 < nat; ia1++ )
+  {
+    for ( int idir1 = 0; idir1 < 3; idir1++ )
+    {
+      for ( int ia2 = 0; ia2 < nat; ia2++ )
+      {
+        for ( int idir2 = 0; idir2 < 3; idir2++ )
+        {
+          int i = 3*ia1+idir1;
+          int j = 3*ia2+idir2;
+          double aij = ( f[1][ia1][idir1][ia2][idir2] - 
+                         f[0][ia1][idir1][ia2][idir2] ) / ( 2.0 * h );
+          aij = aij / sqrt(mass[ia1]*mass[ia2]);
+          a[j*a.n()+i] = aij;
+        }
+      }
+    }
+  }
+
+  // cout << a;
+  valarray<double> w(a.m());
+  DoubleMatrix asave(a);
+  a.syevd('u',w);
+
+  cout << " frequencies:" << endl;
+  for ( int i = 0; i < a.n(); i++ )
+  {
+    if ( w[i] < 0.0 )
+      cout << setw(8) << (int) (Ha2cm1 * sqrt(-w[i])) << " I" << endl;
+    else
+      cout << setw(8) << (int) (Ha2cm1 * sqrt(w[i])) << endl;
+  }
+
+  a = asave;
+  a.syevd('l',w);
+  cout << " frequencies:" << endl;
+  for ( int i = 0; i < a.n(); i++ )
+  {
+    if ( w[i] < 0.0 )
+      cout << setw(8) << (int) (Ha2cm1 * sqrt(-w[i])) << " I" << endl;
+    else
+      cout << setw(8) << (int) (Ha2cm1 * sqrt(w[i])) << endl;
+  }
+}

util/dynmat/fdmoves.sh

+#!/bin/bash
+# use: fdmoves.sh atom1 atom2  .. atomN
+# Note: atom names must appear in the same order as in the Qbox <atomset>
+#
+# generate moves for finite difference calculations of the dynamical matrix
+
+# h: atomic displacement in (a.u.)
+h=0.01
+# nitscf: number of scf iterations for each displacement
+nitscf=20
+
+for atom in ${*}
+do
+  # moves along x direction
+  echo move $atom by  $h 0 0
+  echo run 0 $nitscf
+  echo move $atom by -$h 0 0
+
+  echo move $atom by -$h 0 0
+  echo run 0 $nitscf
+  echo move $atom by  $h 0 0
+
+  # moves along y direction
+  echo move $atom by  0  $h  0
+  echo run 0 $nitscf
+  echo move $atom by  0 -$h  0
+
+  echo move $atom by  0 -$h  0
+  echo run 0 $nitscf
+  echo move $atom by  0  $h  0
+
+  # moves along z direction
+  echo move $atom by  0  0  $h
+  echo run 0 $nitscf
+  echo move $atom by  0  0 -$h
+
+  echo move $atom by  0  0 -$h
+  echo run 0 $nitscf
+  echo move $atom by  0  0  $h
+
+done

util/dynmat/x8664_gcc_serial.mk

+#-------------------------------------------------------------------------------
+#
+# Copyright (c) 2009 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/>.
+#
+#-------------------------------------------------------------------------------
+#
+#  x8664_gcc_serial.mk
+#
+#-------------------------------------------------------------------------------
+# $Id: x8664_gcc_serial.mk,v 1.2 2009/12/04 03:16:17 fgygi Exp $
+#
+ PLT=Linux_x8664
+#-------------------------------------------------------------------------------
+ XERCESCDIR=$(HOME)/software/xml/xerces-c-src_2_8_0
+ FFTWDIR=$(HOME)/software/fftw/Linux_x8664/fftw-2.1.3/fftw
+ BLASDIR=$(HOME)/software/atlas/ATLAS/Linux_P4E64SSE3/lib
+ LAPACKDIR=$(HOME)/software/lapack/LAPACK
+
+ PLTOBJECTS = readTSC.o
+
+ CXX=/usr/bin/g++
+ LD=$(CXX)
+
+ PLTFLAGS += -DIA32 -DUSE_FFTW -D_LARGEFILE_SOURCE \
+             -D_FILE_OFFSET_BITS=64 -DADD_ \
+             -DAPP_NO_THREADS -DXML_USE_NO_THREADS -DUSE_XERCES
+
+
+ INCLUDE = -I$(FFTWDIR) -I$(XERCESCDIR)/include
+
+ CXXFLAGS= -g -Wunused -D$(PLT) $(INCLUDE) $(PLTFLAGS) $(DFLAGS)
+
+ LIBPATH = -L$(GCCDIR)/lib -L$(FFTWDIR)/.libs -L/usr/X11R6/lib \
+           -L$(LAPACKDIR) -L$(BLASDIR) \
+           -L$(XERCESCDIR)/lib
+
+ LIBS =  $(PLIBS) -lfftw \
+         -llapack -lf77blas -latlas -lm \
+         -Xlinker -Bstatic \
+          -lc -lgfortran -static-libgcc -lxerces-c -luuid \
+         -Xlinker -Bdynamic
+
+ LDFLAGS = $(LIBPATH) $(LIBS)
+
+ PLAT=Linux_x8664
+
+#-------------------------------------------------------------------------------