Merge branch 'develop'

.gitignore

+*.[oa]

build/centos6.mk

@@ -60,11 +60,11 @@ ifeq ($(FFT),FFTW3)
  PLTFLAGS += -DFFTWMEASURE
 #PLTFLAGS += -DFFTW_TRANSPOSE
  PLTFLAGS += -DFFTW3_2D
- FFTWDIR=$(HOME)/software/fftw/fftw-3.3.4
- FFTWINCLUDEDIR=$(FFTWDIR)/api
- FFTWLIBDIR=$(FFTWDIR)/.libs
- INCLUDE += -I$(FFTWINCLUDEDIR)
- LIBPATH += -L$(FFTWLIBDIR)
+#FFTWDIR=$(HOME)/software/fftw/fftw-3.3.4
+#FFTWINCLUDEDIR=$(FFTWDIR)/api
+#FFTWLIBDIR=$(FFTWDIR)/.libs
+#INCLUDE += -I$(FFTWINCLUDEDIR)
+#LIBPATH += -L$(FFTWLIBDIR)
  LIBS += -lfftw3
 endif
 

src/AlphaRSH.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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// AlphaRSH.h
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALPHARSH_H
+#define ALPHARSH_H
+
+#include<iostream>
+#include<iomanip>
+#include<sstream>
+#include<stdlib.h>
+
+#include "Sample.h"
+
+class AlphaRSH : public Var
+{
+  Sample *s;
+
+  public:
+
+  const char *name ( void ) const { return "alpha_RSH"; };
+
+  int set ( int argc, char **argv )
+  {
+    if ( argc != 2 )
+    {
+      if ( ui->onpe0() )
+      cout << " alpha_RSH takes only one value" << endl;
+      return 1;
+    }
+
+    double v = atof(argv[1]);
+    if ( v < 0.0 )
+    {
+      if ( ui->onpe0() )
+        cout << " alpha_RSH must be non-negative" << endl;
+      return 1;
+    }
+
+    s->ctrl.alpha_RSH = v;
+
+    return 0;
+  }
+
+  string print (void) const
+  {
+     ostringstream st;
+     st.setf(ios::left,ios::adjustfield);
+     st << setw(10) << name() << " = ";
+     st.setf(ios::right,ios::adjustfield);
+     st << s->ctrl.alpha_RSH;
+     return st.str();
+  }
+
+  AlphaRSH(Sample *sample) : s(sample)
+  {
+    s->ctrl.alpha_RSH = 0;
+  }
+};
+#endif

src/AngleConstraint.C

@@ -80,7 +80,7 @@ vector<vector<double> > &rp) const
   D3vector r3(pr3);
   D3vector g1,g2,g3;
   grad_sigma(r1,r2,r3,g1,g2,g3);
-  const double a = bond_angle(r1,r2,r3);
+  // const double a = bond_angle(r1,r2,r3);
   double ng = g1*g1 + g2*g2 + g3*g3;
   assert(ng>=0.0);
 
@@ -358,9 +358,9 @@ void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2,
     // angle is large enough. Use finite differences
     //cout << " ========= grad_sigma using finite differences" << endl;
 
-    const double r12_inv = 1.0/length(r12);
-    const double r32_inv = 1.0/length(r32);
-    const double a = bond_angle(r1,r2,r3);
+    // const double r12_inv = 1.0/length(r12);
+    // const double r32_inv = 1.0/length(r32);
+    // const double a = bond_angle(r1,r2,r3);
 
     const double l12 = length(r1-r2);
     const double l32 = length(r3-r2);

src/AtomSet.C

@@ -503,6 +503,25 @@ void AtomSet::randomize_velocities(double temp)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
+D3vector AtomSet::rcm(void) const
+{
+  D3vector mrsum;
+  double msum = 0.0;
+  for ( int is = 0; is < atom_list.size(); is++ )
+  {
+    double mass = species_list[is]->mass();
+    for ( int ia = 0; ia < atom_list[is].size(); ia++ )
+    {
+      D3vector r = atom_list[is][ia]->position();
+      mrsum += mass * r;
+      msum += mass;
+    }
+  }
+  if ( msum == 0.0 ) return D3vector(0,0,0);
+  return mrsum / msum;
+}
+
+////////////////////////////////////////////////////////////////////////////////
 D3vector AtomSet::vcm(void) const
 {
   D3vector mvsum;
@@ -542,6 +561,141 @@ void AtomSet::reset_vcm(void)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
+void AtomSet::reset_rotation(void)
+{
+  // check for special case of zero or one atom
+  if ( size() < 2 ) return;
+
+  D3vector rc = rcm();
+  D3vector vc = vcm();
+  vector<vector<double> > rt;
+  get_positions(rt);
+  vector<vector<double> > vt;
+  get_velocities(vt);
+
+  // compute angular momentum w.r.t. the center of mass
+  D3vector L;
+  for ( int is = 0; is < vt.size(); is++ )
+  {
+    double mass = species_list[is]->mass();
+    for ( int ia = 0; ia < na(is); ia++ )
+    {
+      D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
+      r -= rc;
+      D3vector v(vt[is][3*ia+0],vt[is][3*ia+1],vt[is][3*ia+2]);
+      v -= vc;
+      L += mass * ( r ^ v );
+    }
+  }
+
+  // compute inertia tensor a and vector omega
+
+  // check for special case of all atoms aligned, for which the
+  // inertia tensor has a zero eigenvalue
+
+  // check if all atoms are aligned
+  // collect all positions in a single vector of D3vectors
+  vector<D3vector> rv(size());
+  int iat = 0;
+  for ( int is = 0; is < vt.size(); is++ )
+    for ( int ia = 0; ia < na(is); ia++ )
+      rv[iat++] = D3vector(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
+
+  // normalized direction e = (rv[1]-rv[0])
+  D3vector e = normalized(rv[1]-rv[0]);
+  bool aligned = true;
+  for ( int i = 2; (i < size()) && aligned; i++ )
+  {
+    D3vector u = normalized(rv[i]-rv[0]);
+    aligned &= length(u^e) < 1.e-6;
+  }
+
+  D3vector omega;
+  // compute the inertia tensor
+  // treat separately the case of all atoms aligned
+  if ( aligned )
+  {
+    // inertia tensor reduces to a scalar
+    double a = 0.0;
+    for ( int is = 0; is < vt.size(); is++ )
+    {
+      double mass = species_list[is]->mass();
+      for ( int ia = 0; ia < na(is); ia++ )
+      {
+        D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
+        r -= rc;
+        a += mass * norm2(r);
+      }
+    }
+    omega = L / a;
+  }
+  else
+  {
+    double a00,a01,a02,a10,a11,a12,a20,a21,a22;
+    a00=a01=a02=a10=a11=a12=a20=a21=a22=0.0;
+
+    for ( int is = 0; is < vt.size(); is++ )
+    {
+      double mass = species_list[is]->mass();
+      for ( int ia = 0; ia < na(is); ia++ )
+      {
+        D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
+        r -= rc;
+        a00 += mass * ( r.y * r.y + r.z * r.z );
+        a11 += mass * ( r.x * r.x + r.z * r.z );
+        a22 += mass * ( r.x * r.x + r.y * r.y );
+        a01 -= mass * r.x * r.y;
+        a02 -= mass * r.x * r.z;
+        a12 -= mass * r.y * r.z;
+      }
+    }
+    a10 = a01;
+    a20 = a02;
+    a21 = a12;
+
+    // inverse b of the inertia tensor a
+    // determinant
+    double det = a00 * ( a11 * a22 - a21 * a12 ) -
+                 a01 * ( a10 * a22 - a20 * a12 ) +
+                 a02 * ( a10 * a21 - a20 * a11 );
+    // the determinant must be positive
+    assert(det>1.e-8);
+    double b00,b01,b02,b10,b11,b12,b20,b21,b22;
+    b00 = ( -a12*a21 + a11*a22 ) / det;
+    b10 = (  a12*a20 - a10*a22 ) / det;
+    b20 = ( -a11*a20 + a10*a21 ) / det;
+
+    b01 = (  a02*a21 - a01*a22 ) / det;
+    b11 = ( -a02*a20 + a00*a22 ) / det;
+    b21 = (  a01*a20 - a00*a21 ) / det;
+
+    b02 = ( -a02*a11 + a01*a12 ) / det;
+    b12 = (  a02*a10 - a00*a12 ) / det;
+    b22 = ( -a01*a10 + a00*a11 ) / det;
+
+    // omega = inverse(a) * L
+    omega.x = b00 * L.x + b01 * L.y + b02 * L.z;
+    omega.y = b10 * L.x + b11 * L.y + b12 * L.z;
+    omega.z = b20 * L.x + b21 * L.y + b22 * L.z;
+  }
+
+  // correct velocities: v = v - omega ^ r
+  for ( int is = 0; is < vt.size(); is++ )
+  {
+    for ( int ia = 0; ia < na(is); ia++ )
+    {
+      D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
+      r -= rc;
+      D3vector dv = omega ^ r;
+      vt[is][3*ia+0] -= dv.x;
+      vt[is][3*ia+1] -= dv.y;
+      vt[is][3*ia+2] -= dv.z;
+    }
+  }
+  set_velocities(vt);
+}
+
+////////////////////////////////////////////////////////////////////////////////
 void AtomSet::fold_in_ws(void)
 {
   vector<vector<double> > p;

src/AtomSet.h

@@ -85,10 +85,12 @@ class AtomSet
   void rescale_velocities(double fac);
   void randomize_velocities(double temp);
   void randomize_positions(double amplitude);
+  D3vector rcm(void) const;
   D3vector vcm(void) const;
   D3vector dipole(void) const;
   D3tensor quadrupole(void) const;
   void reset_vcm(void);
+  void reset_rotation(void);
   void fold_in_ws(void);
   int size(void) const;
  };

src/BOSampleStepper.C

@@ -862,7 +862,7 @@ void BOSampleStepper::step(int niter)
         // if ( onpe0 && nite_ > 0 )
         //   cout << " delta_ehart = " << delta_ehart << endl;
         int ite = 0;
-        double energy, etotal_int;
+        double etotal_int;
 
         double eigenvalue_sum, eigenvalue_sum_m = 0.0;
         // if nite == 0: do 1 iteration, no screening in charge mixing
@@ -871,7 +871,7 @@ void BOSampleStepper::step(int niter)
         while ( !nonscf_converged && ite < max(nite_,1) )
         {
           tmap["energy"].start();
-          energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
+          ef_.energy(true,dwf,false,fion,false,sigma_eks);
           tmap["energy"].stop();
 
           // compute the sum of eigenvalues (with fixed weight)

src/Base64Transcoder.C

@@ -120,7 +120,7 @@ int Base64Transcoder::decode(const int nchars, const char* const from,
   // nchars: number of chars in array "from"
   // the number of bytes successfully translated is returned
 
-  byte a0,a1,a2,a3,b0,b1,b2,b3;
+  byte a2,a3,b0,b1,b2,b3;
   int c;
   const char* fptr = from;
   const char* const fptr_end = from+nchars+1;
@@ -142,7 +142,7 @@ int Base64Transcoder::decode(const int nchars, const char* const from,
 #endif
       break;
     }
-    a0 = (byte) c;
+    // a0 = (byte) c;
     b0 = (byte) dtable[c];
 
     do
@@ -158,7 +158,7 @@ int Base64Transcoder::decode(const int nchars, const char* const from,
 #endif
       break;
     }
-    a1 = (byte) c;
+    // a1 = (byte) c;
     b1 = (byte) dtable[c];
 
     do

src/BetaRSH.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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// BetaRSH.h
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef BETARSH_H
+#define BETARSH_H
+
+#include<iostream>
+#include<iomanip>
+#include<sstream>
+#include<stdlib.h>
+
+#include "Sample.h"
+
+class BetaRSH : public Var
+{
+  Sample *s;
+
+  public:
+
+  const char *name ( void ) const { return "beta_RSH"; };
+
+  int set ( int argc, char **argv )
+  {
+    if ( argc != 2 )
+    {
+      if ( ui->onpe0() )
+      cout << " beta_RSH takes only one value" << endl;
+      return 1;
+    }
+
+    double v = atof(argv[1]);
+    if ( v < 0.0 )
+    {
+      if ( ui->onpe0() )
+        cout << " beta_RSH must be non-negative" << endl;
+      return 1;
+    }
+
+    s->ctrl.beta_RSH = v;
+
+    return 0;
+  }
+
+  string print (void) const
+  {
+     ostringstream st;
+     st.setf(ios::left,ios::adjustfield);
+     st << setw(10) << name() << " = ";
+     st.setf(ios::right,ios::adjustfield);
+     st << s->ctrl.beta_RSH;
+     return st.str();
+  }
+
+  BetaRSH(Sample *sample) : s(sample)
+  {
+    s->ctrl.beta_RSH = 0.25;
+  }
+};
+#endif

src/Bisection.C

@@ -401,7 +401,8 @@ void Bisection::compute_transform(const SlaterDet& sd, int maxsweep, double tol)
   // adiag_ is resized by jade
 
   // diagonalize projectors
-  int nsweep = jade(maxsweep,tol,amat_,*u_,adiag_);
+  // int nsweep = jade(maxsweep,tol,amat_,*u_,adiag_);
+  jade(maxsweep,tol,amat_,*u_,adiag_);
 #ifdef TIMING
   if ( ctxt_.onpe0() )
     cout << "Bisection::compute_transform: nsweep=" << nsweep

src/BisectionCmd.h

@@ -65,6 +65,26 @@ class BisectionCmd : public Cmd
     nLevels[1]=atoi(argv[2]);
     nLevels[2]=atoi(argv[3]);
 
+    if ( epsilon < 0.0 )
+    {
+      if ( ui->onpe0() )
+      {
+        cout << " BisectionCmd: threshold must be non-negative" << endl;
+      }
+      return 1;
+    }
+
+    if ( nLevels[0] < 0 || nLevels[0] > 5 ||
+         nLevels[1] < 0 || nLevels[1] > 5 ||
+         nLevels[2] < 0 || nLevels[2] > 5 )
+    {
+      if ( ui->onpe0() )
+      {
+        cout << " BisectionCmd: levels must be in [0,5]" << endl;
+      }
+      return 1;
+    }
+
     tm.reset();
     for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
     {

src/BlHF.h

@@ -46,10 +46,10 @@ class BlHF : public Var
     int v0 = atoi(argv[1]);
     int v1 = atoi(argv[2]);
     int v2 = atoi(argv[3]);
-    if ( v0 < 0 || v1 < 0 || v2 < 0 )
+    if ( v0 < 0 || v1 < 0 || v2 < 0 || v0 > 5 || v1 > 5 || v2 > 5 )
     {
       if ( ui->onpe0() )
-        cout << " blHF values must be > 0" << endl;
+        cout << " blHF values must be in [0,5]" << endl;
       return 1;
     }
 

src/CPSampleStepper.C

@@ -134,7 +134,6 @@ 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);
 
@@ -254,7 +253,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();
+    ef_.enthalpy();
 
     if ( s_.ctxt_.mype() == 0 )
       cout << "</iteration>" << endl;

src/ChargeDensity.C

@@ -21,6 +21,7 @@
 #include "Wavefunction.h"
 #include "FourierTransform.h"
 #include "SlaterDet.h"
+#include "blas.h" // dasum
 
 #include <iomanip>
 #include <algorithm> // fill
@@ -207,6 +208,18 @@ void ChargeDensity::update_rhor(void)
       for ( int i = 0; i < rhor_size; i++ )
         prhor[i] = rhotmp[i].real() * omega_inv;
     }
+
+    // integral of the charge density
+    tmap["charge_integral"].start();
+    int ione=1;
+    int n = rhor_size;
+    double sum = dasum(&n,prhor,&ione);
+    sum *= omega / vft_->np012();
+
+    // sum on all indices except spin: sum along columns of spincontext
+    wf_.spincontext()->dsum('c',1,1,&sum,1);
+    tmap["charge_integral"].stop();
+    total_charge_[ispin] = sum;
   }
 }
 

src/ConstraintSet.C

@@ -27,6 +27,7 @@
 #include <iostream>
 #include <iomanip>
 #include <cstdlib> // atof
+#include <cstring> // strcmp
 using namespace std;
 
 const int constraints_maxiter = 50;

src/Control.h

@@ -54,6 +54,9 @@ struct Control
 
   std::string xc;
   double alpha_PBE0;
+  double alpha_RSH;
+  double beta_RSH;
+  double mu_RSH;
   std::string spin;
   int delta_spin;
 

src/EnergyFunctional.C

@@ -409,7 +409,6 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
       const double weight = wf.weight(ikp);
       const SlaterDet& sd = *(wf.sd(ispin,ikp));
       const Basis& wfbasis = sd.basis();
-      const D3vector kp = wfbasis.kpoint();
       // factor fac in next lines: 2.0 for G and -G (if basis is real) and
       // 0.5 from 1/(2m)
       const double fac = wfbasis.real() ? 1.0 : 0.5;

src/ExchangeOperator.h

@@ -33,8 +33,6 @@ class ExchangeOperator
   double eex_;
   // constant of support function for exchange integration
   double rcut_;
-  // mixing coefficient for exchange energy and dwf accumulation
-  double HFCoeff_;
 
   // HF stress tensor
   std::valarray<double> sigma_exhf_;
@@ -74,10 +72,8 @@ class ExchangeOperator
   // G vectors
   valarray<double> qpG21_;
   valarray<double> qpG22_;
-  valarray<double> qpG2i1_;
-  valarray<double> qpG2i2_;
-  valarray<double> G2_;
-  valarray<double> G2i_;
+  valarray<double> int_pot1_;
+  valarray<double> int_pot2_;
 
   // numbers of states
   int nLocalStates_;
@@ -173,18 +169,23 @@ class ExchangeOperator
   vector<DoubleMatrix*> uc_;
   vector<long int> localization_;
 
+  // screened interaction potential paramters
+  double alpha_sx_, beta_sx_, mu_sx_;
+  // interaction potential. g2 is the wave vector squared.
+  double vint(double g2);
+  // derivative of the interaction potential w.r.t. g2
+  double dvint(double g2);
+  double vint_div_scal(double rc);
+
   public:
 
   // constructor
-  ExchangeOperator(Sample& s_, double HFCoeff);
+  // screened interaction potential: alpha*erf(mu*r)/r + beta*erfc(mu*r)/r
+  ExchangeOperator(Sample& s_, double alpha_sx, double beta_sx, double mu_sx);
 
   // destructor
   ~ExchangeOperator();
 
-  // parameters
-  void setmixCoeff(double value) { HFCoeff_ = value; };
-  double HFCoeff() { return HFCoeff_; };
-
   // exchange energy and forces computation
   double eex() { return eex_; };
   double update_energy(bool compute_stress);

src/ExponentialIntegral.C

+////////////////////////////////////////////////////////////////////////////////
+//
+// 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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// ExponentialIntegral.C
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Calculate exponential integral using the algorithm of
+// Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include "ExponentialIntegral.h"
+#include <cmath>
+#include <algorithm>
+#include <vector>
+
+using namespace std;
+
+namespace util
+{
+
+// Calculate the exponential integral E_1(x):
+//
+//        inf
+//          /     -t
+//         |     e
+// E (x) = | dt -----
+//  1      |      t
+//        /
+//         x
+//
+// Input:  x - position at which exponential integral is evaluated (x > 0)
+// Return: E1(x)
+double E1(const double x)
+{
+
+  if (x < series_cutoff)
+  {
+    // use series expansion
+    return E1_series(x);
+  }
+  else
+  {
+    // use gauss_laguerre expression
+    return exp(-x) * gauss_laguerre(x);
+  }
+}
+
+// Series expansion of the exponential integral
+//
+//                          n_cut
+//                          -----     n  n
+//                           \    (-1)  x
+// E (x) = -gamma - ln(x) -   )   --------
+//  1                        /     n * n!
+//                          -----
+//                          n = 1
+//
+// where gamma is the Euler constant.
+// n_cut is set to 25
+// Input:  x - position at which exponential integral is evaluated (x > 0)
+// Return: approximation by series expansion for E_1(x)
+double E1_series(const double x)
+{
+  // Euler constant
+  const double EULER_GAMMA = 0.57721566490153286060651209008241;
+
+  // Cutoff for series expansion
+  const int itmax = 25;
+
+  // initialize summation result
+  double res = 0.0;
+
+  // perform the summation
+  for (int it = itmax; it > 1; it--)
+  {
+    // calculate 1/n
+    const double fact = 1.0 / it;
+    // add next term of summation
+    res = x * fact * (fact - res);
+  }
+
+  // add everything up
+  return -EULER_GAMMA - log(x) + x * (1.0 - res);
+}
+
+// The Gauss Laguerre expansion of the exponential integral can be written as
+//
+//             N
+// E (x0)    -----     a
+//  1         \         n
+// ------ =    )   ---------
+//   -x0      /     x  + x0
+//  e        -----   n
+//            n=1
+//
+// where the a_n and x_n are determined by least quadrature (see reference)
+// Input: x0 - point at which Gaussian Laguerre quadrature is calculated
+// Return: E_1(x0) / exp(-x0) in this approximation
+double gauss_laguerre(const double x0)
+{
+  // initialize constants a_n and x_n
+  const double size = 15;
+  const double a[] = { 0.2182348859400869e+00, 0.3422101779228833e+00,
+      0.2630275779416801e+00, 0.1264258181059305e+00, 0.4020686492100091e-01,
+      0.8563877803611838e-02, 0.1212436147214252e-02, 0.1116743923442519e-03,