merge trunk into efield branch

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

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

src/AndersonMixer.C

@@ -15,14 +15,41 @@
 // AndersonMixer.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: AndersonMixer.C,v 1.12 2009-09-08 14:26:01 fgygi Exp $
 
 #include "AndersonMixer.h"
 #include "blas.h"
 #include <iostream>
+#if USE_MPI
+#include <mpi.h>
+#else
+typedef int MPI_Comm;
+#endif
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
+AndersonMixer::AndersonMixer(const int m, const int nmax,
+  const MPI_Comm* const pcomm) : m_(m), nmax_(nmax), pcomm_(pcomm)
+{
+#if USE_MPI
+  if ( pcomm_ != 0 )
+  {
+    MPI_Comm_rank(*pcomm_,&mype_);
+    MPI_Comm_size(*pcomm_,&npes_);
+  }
+#endif
+
+  assert( nmax >= 0 );
+  x_.resize(nmax_+1);
+  f_.resize(nmax_+1);
+  for ( int n = 0; n < nmax_+1; n++ )
+  {
+    x_[n].resize(m_);
+    f_[n].resize(m_);
+  }
+  restart();
+}
+
+////////////////////////////////////////////////////////////////////////////////
 void AndersonMixer::restart(void)
 {
   n_ = -1;
@@ -52,15 +79,17 @@ void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
     f_[k_][i] = f[i];
   }
 
-  valarray<double> a;
-  valarray<double> b;
+  valarray<double> a,atmp;
+  valarray<double> b,btmp;
   valarray<double> theta;
   if ( n_ > 0 )
   {
     // compute matrix A = F^T F and rhs b = F^T f
     // compute the lower part of A only (i>=j)
     a.resize(n_*n_);
+    atmp.resize(n_*n_);
     b.resize(n_);
+    btmp.resize(n_);
     theta.resize(n_);
     for ( int i = 0; i < n_; i++ )
     {
@@ -85,42 +114,19 @@ void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
       b[i] = bsum;
     }
 
-    if ( pctxt_ != 0 )
+#if USE_MPI
+    if ( pcomm_ != 0 )
     {
-      pctxt_->dsum(n_*n_,1,&a[0],n_*n_);
-      pctxt_->dsum(n_,1,&b[0],n_);
-    }
-
-#if 0
-    // print matrix a and rhs b
-    if ( pctxt_ != 0 )
-    {
-      for ( int ip =0; ip < pctxt_->size(); ip++ )
-      {
-        pctxt_->barrier();
-        if ( pctxt_->mype() == ip )
-        {
-          // print matrix a and rhs b
-          double anrm = 0.0;
-          for ( int i = 0; i < n_; i++ )
-            for ( int j = 0; j <=i; j++ )
-            {
-              cout << pctxt_->mype() << ": "
-                   << "a("<<i<<","<<j<<")=" << a[i+j*n_] << endl;
-              anrm += a[i+j*n_]*a[i+j*n_];
-            }
-            for ( int i = 0; i < n_; i++ )
-              cout << pctxt_->mype() << ": "
-                   <<  "b("<<i<<")=" << b[i] << endl;
-          //cout << " AndersonMixer: n=" << n_ << " anorm = " << anrm << endl;
-        }
-      }
+      MPI_Allreduce(&a[0],&atmp[0],n_*n_,MPI_DOUBLE,MPI_SUM,*pcomm_);
+      a = atmp;
+      MPI_Allreduce(&b[0],&btmp[0],n_,MPI_DOUBLE,MPI_SUM,*pcomm_);
+      b = btmp;
     }
 #endif
 
     // solve the linear system a * theta = b
     // solve on task 0 and bcast result
-    if ( pctxt_ == 0 || pctxt_->onpe0() )
+    if ( pcomm_ == 0 || mype_ == 0 )
     {
       const bool diag = false;
       if ( diag )
@@ -241,27 +247,11 @@ void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
 
     }
 
+#if USE_MPI
     // broadcast theta from task 0
-    if ( pctxt_ != 0 )
-    {
-      if ( pctxt_->onpe0() )
-        pctxt_->dbcast_send(n_,1,&theta[0],n_);
-      else
-        pctxt_->dbcast_recv(n_,1,&theta[0],n_,0,0);
-    }
-
-#if 0
-    for ( int ip = 0; ip < pctxt_->size(); ip++ )
+    if ( pcomm_ != 0 )
     {
-      pctxt_->barrier();
-      if ( pctxt_->mype() == ip )
-      {
-        cout << pctxt_->mype() << ": ";
-        cout << " AndersonMixer: theta = ";
-        for ( int i = 0; i < theta.size(); i++ )
-          cout << theta[i] << " ";
-        cout << endl;
-      }
+      MPI_Bcast(&theta[0],n_,MPI_DOUBLE,0,*pcomm_);
     }
 #endif
 

src/AndersonMixer.h

@@ -15,7 +15,6 @@
 // AndersonMixer.h
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: AndersonMixer.h,v 1.8 2009-03-08 01:10:30 fgygi Exp $
 
 #ifndef ANDERSONMIXER_H
 #define ANDERSONMIXER_H
@@ -23,7 +22,11 @@
 #include <vector>
 #include <valarray>
 #include <cassert>
-#include "Context.h"
+#ifdef USE_MPI
+#include <mpi.h>
+#else
+typedef int MPI_Comm;
+#endif
 
 class AndersonMixer
 {
@@ -34,26 +37,15 @@ class AndersonMixer
   int     nmax_;                 // maximum number of vectors (without current)
   int     n_;                    // number of vectors
   int     k_;                    // index of current vector
-  const   Context* const pctxt_; // pointer to relevant Context, null if local
+  const   MPI_Comm* const pcomm_;// pointer to relevant Context, null if local
+  int     mype_;
+  int     npes_;
 
   std::vector<std::valarray<double> > x_,f_;
 
   public:
 
-  AndersonMixer(const int m, const int nmax, const Context* const pctxt) :
-    m_(m), nmax_(nmax), pctxt_(pctxt)
-  {
-    assert( nmax >= 0 );
-    x_.resize(nmax_+1);
-    f_.resize(nmax_+1);
-    for ( int n = 0; n < nmax_+1; n++ )
-    {
-      x_[n].resize(m_);
-      f_[n].resize(m_);
-    }
-    restart();
-  }
-
+  AndersonMixer(const int m, const int nmax, const MPI_Comm* const pcomm);
   void update(double* x, double* f, double* xbar, double* fbar);
   void restart(void);
 };

src/AngleCmd.h

@@ -32,29 +32,53 @@ class AngleCmd : public Cmd
 
   AngleCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "angle"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "angle"; }
+  const char *help_msg(void) const
   {
     return
     "\n angle\n\n"
-    " syntax: angle name1 name2 name3\n\n"
-    "   The angle command prints the angle defined by three atoms.\n\n";
+    " syntax: angle [-pbc] name1 name2 name3\n\n"
+    "   The angle command prints the angle defined by three atoms.\n"
+    "   If the -pbc option is used, the angle is computed using the\n"
+    "   nearest atoms taking into account periodic boundary conditions.\n\n";
   }
 
   int action(int argc, char **argv)
   {
-    if ( argc != 4 )
+    if ( ! ( argc == 4 || argc == 5 ) )
     {
       if ( ui->onpe0() )
       {
-        cout << " use: angle name1 name2 name3" << endl;
+        cout << " use: angle [-pbc] name1 name2 name3" << endl;
       }
       return 1;
     }
 
-    string name1(argv[1]);
-    string name2(argv[2]);
-    string name3(argv[3]);
+    string name1, name2, name3;
+    bool use_pbc = false;
+
+    if ( argc == 4 )
+    {
+      name1 = argv[1];
+      name2 = argv[2];
+      name3 = argv[3];
+    }
+    if ( argc == 5 )
+    {
+      if ( strcmp(argv[1],"-pbc") )
+      {
+        if ( ui->onpe0() )
+        {
+          cout << " use: angle [-pbc] name1 name2 name3" << endl;
+        }
+        return 1;
+      }
+      use_pbc = true;
+      name1 = argv[2];
+      name2 = argv[3];
+      name3 = argv[4];
+    }
+
     Atom* a1 = s->atoms.findAtom(name1);
     Atom* a2 = s->atoms.findAtom(name2);
     Atom* a3 = s->atoms.findAtom(name3);
@@ -82,22 +106,25 @@ class AngleCmd : public Cmd
       return 1;
     }
 
-    D3vector r12(a1->position()-a2->position());
-    D3vector r32(a3->position()-a2->position());
-    if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 )
+    if ( ui->onpe0() )
     {
-      if ( ui->onpe0() )
+      D3vector r12(a1->position()-a2->position());
+      D3vector r32(a3->position()-a2->position());
+      if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 )
       {
         cout << " AngleCmd: atoms are too close" << endl;
+        return 1;
       }
-      return 1;
-    }
 
-    const double sp = normalized(r12) * normalized(r32);
-    const double c = max(-1.0,min(1.0,sp));
-    const double a = (180.0/M_PI)*acos(c);
-    if ( ui->onpe0() )
-    {
+      if ( use_pbc )
+      {
+        const UnitCell& cell = s->wf.cell();
+        cell.fold_in_ws(r12);
+        cell.fold_in_ws(r32);
+      }
+      const double sp = normalized(r12) * normalized(r32);
+      const double c = max(-1.0,min(1.0,sp));
+      const double a = (180.0/M_PI)*acos(c);
       cout.setf(ios::fixed,ios::floatfield);
       cout << " angle " << name1 << "-" << name2  << "-" << name3
            << ": "

src/AtomCmd.h

@@ -34,8 +34,8 @@ class AtomCmd : public Cmd
 
   AtomCmd(Sample *sample) : s(sample) {};
 
-  char *name(void) const { return "atom"; }
-  char *help_msg(void) const
+  const char *name(void) const { return "atom"; }
+  const char *help_msg(void) const
   {
     return
     "\n atom\n\n"

src/AtomSet.C

@@ -496,6 +496,7 @@ void AtomSet::randomize_velocities(double temp)
       v[is][3*ia+2] = width * xi2;
     }
   }
+  reset_vcm();
   set_velocities(v);
 }
 

src/AtomSetHandler.C

@@ -15,9 +15,6 @@
 // AtomSetHandler.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: AtomSetHandler.C,v 1.13 2008-09-08 15:56:18 fgygi Exp $
-
-#if USE_XERCES
 
 #include "AtomSetHandler.h"
 #include "AtomSet.h"
@@ -104,20 +101,6 @@ void AtomSetHandler::endElement(const XMLCh* const uri,
   istringstream stst(content);
   if ( locname == "unit_cell")
   {
-    event_type event = unit_cell;
-    as_.context().ibcast_send(1,1,(int*)&event,1);
-    // notify listening nodes
-    double buf[9];
-    buf[0] = as_.cell().a(0).x;
-    buf[1] = as_.cell().a(0).y;
-    buf[2] = as_.cell().a(0).z;
-    buf[3] = as_.cell().a(1).x;
-    buf[4] = as_.cell().a(1).y;
-    buf[5] = as_.cell().a(1).z;
-    buf[6] = as_.cell().a(2).x;
-    buf[7] = as_.cell().a(2).y;
-    buf[8] = as_.cell().a(2).z;
-    as_.context().dbcast_send(9,1,buf,9);
   }
   else if ( locname == "atom")
   {
@@ -143,21 +126,6 @@ void AtomSetHandler::endElement(const XMLCh* const uri,
       throw;
     }
 
-    // notify listening nodes and broadcast atom info
-    event_type event = atom;
-    as_.context().ibcast_send(1,1,(int*)&event,1);
-    as_.context().string_bcast(current_atom_name,0);
-    as_.context().string_bcast(current_atom_species,0);
-    double buf[3];
-    buf[0] = current_atom_position.x;
-    buf[1] = current_atom_position.y;
-    buf[2] = current_atom_position.z;
-    as_.context().dbcast_send(3,1,buf,3);
-    buf[0] = current_atom_velocity.x;
-    buf[1] = current_atom_velocity.y;
-    buf[2] = current_atom_velocity.z;
-    as_.context().dbcast_send(3,1,buf,3);
-
   }
   else if ( locname == "position" )
   {
@@ -193,7 +161,7 @@ StructureHandler* AtomSetHandler::startSubHandler(const XMLCh* const uri,
     }
 
     // delegate to SpeciesHandler
-    current_species = new Species(as_.context(),current_species_name);
+    current_species = new Species(current_species_name);
     return new SpeciesHandler(*current_species);
   }
   else
@@ -211,7 +179,7 @@ void AtomSetHandler::endSubHandler(const XMLCh* const uri,
   // cout << " AtomSetHandler::endSubHandler " << locname << endl;
   if ( locname == "species" )
   {
-    SpeciesReader sp_reader(current_species->context());
+    SpeciesReader sp_reader;
 
     // check if only the uri was provided
     if ( current_species->uri() != "" )
@@ -219,22 +187,9 @@ void AtomSetHandler::endSubHandler(const XMLCh* const uri,
       // href was found in species definition
       // attempt to read the species from that uri
 
-      try
-      {
-        sp_reader.readSpecies(*current_species,current_species->uri());
-      }
-      catch ( const SpeciesReaderException& e )
-      {
-        cout << " SpeciesReaderException caught in AtomSetHandler" << endl;
-      }
+      sp_reader.uri_to_species(current_species->uri(),*current_species);
     }
 
-    // notify listening nodes and broadcast species info
-    event_type event = species;
-    as_.context().ibcast_send(1,1,(int*)&event,1);
-    as_.context().string_bcast(current_species_name,0);
-    sp_reader.bcastSpecies(*current_species);
-
     // cout << "AtomSetHandler::endSubHandler: adding Species:"
     //      << current_species_name << endl;
     try
@@ -250,5 +205,3 @@ void AtomSetHandler::endSubHandler(const XMLCh* const uri,
   }
   delete last;
 }
-
-#endif

src/AtomsDyn.h

@@ -33,7 +33,7 @@ class AtomsDyn : public Var
 
   public:
 
-  char *name ( void ) const { return "atoms_dyn"; };
+  const char *name ( void ) const { return "atoms_dyn"; };
 
   int set ( int argc, char **argv )
   {

src/BOSampleStepper.h

@@ -23,6 +23,7 @@
 #include "SampleStepper.h"
 #include "EnergyFunctional.h"
 #include "Sample.h"
+#include "ChargeDensity.h"
 #include "Wavefunction.h"
 
 class WavefunctionStepper;

src/Basis.h

@@ -21,22 +21,27 @@
 
 #include "D3vector.h"
 #include "UnitCell.h"
-#include "Context.h"
 #include <vector>
 
+#ifdef USE_MPI
+#include <mpi.h>
+#else
+  typedef int MPI_Comm;
+#endif
+
 class Basis
 {
   private:
 
-  Context ctxt_;
-  int nprow_, myrow_;
+  MPI_Comm comm_;
+  int npes_, mype_;
 
   UnitCell cell_;         // cell dimensions
   UnitCell refcell_;      // reference cell dimensions
   D3vector kpoint_;       // k-point in units of b0,b1,b2
   double ecut_;           // energy cutoff of wavefunctions in Rydberg
-  int idxmin_[3];          // minimum index in each direction
-  int idxmax_[3];          // maximum index in each direction
+  int idxmin_[3];         // minimum index in each direction
+  int idxmax_[3];         // maximum index in each direction
   int size_;              // basis size
   int nrods_;             // total number of rods
   std::vector<int> localsize_; // localsize_[ipe]
@@ -57,7 +62,7 @@ class Basis
   std::vector<double> kpg2_;  // 2-norm of g vectors g2[localsize]
   std::vector<double> g2i_;   // inverse square norm of g vec g2i[localsize]
   std::vector<double> kpg2i_; // inverse square norm of k+g vec kpg2i[localsize]
-  int np_[3];            // cache for the function np
+  int np_[3];                 // cache for the function np
   std::vector<double> gx_;    // g vec components gx[j*localsize+i], j=0,1,2
   std::vector<double> kpgx_;  // k+g vec components kpgx[j*localsize+i], j=0,1,2
   std::vector<int> isort_loc; // index array to access locally sorted vectors
@@ -67,9 +72,11 @@ class Basis
 
   public:
 
-  const Context& context(void) const; // context on which Basis is defined
+  MPI_Comm comm(void) const;     // MPI_Comm on which Basis is defined
+  int mype(void) const { return mype_; }
+  int npes(void) const { return npes_; }
 
-  const UnitCell& cell() const;   // cell dimensions
+  const UnitCell& cell() const;  // cell dimensions
   const UnitCell& refcell() const;// reference cell dimensions
   const D3vector kpoint() const; // k-point in units of b0,b1,b2
   int np(int i) const;           // good size of FFT grid in direction i
@@ -135,7 +142,7 @@ class Basis
   double memsize(void) const;
   double localmemsize(void) const;
 
-  Basis(const Context &ctxt, D3vector kpoint);
+  Basis(MPI_Comm comm, D3vector kpoint);
   //Basis(const Basis &b);
   ~Basis(void);
 

src/BasisMapping.C

@@ -15,7 +15,6 @@
 //  BasisMapping.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: BasisMapping.C,v 1.4 2008-09-08 15:56:18 fgygi Exp $
 
 #include "Basis.h"
 #include "Context.h"
@@ -26,13 +25,10 @@
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
-BasisMapping::BasisMapping (const Basis &basis) : ctxt_(basis.context()),
- basis_(basis)
-
+BasisMapping::BasisMapping (const Basis &basis) : basis_(basis)
 {
-  assert(ctxt_.npcol() == 1);
-  nprocs_ = ctxt_.size();
-  myproc_ = ctxt_.myproc();
+  nprocs_ = basis_.npes();
+  myproc_ = basis_.mype();
 
   np0_ = basis.np(0);
   np1_ = basis.np(1);
@@ -428,11 +424,11 @@ void BasisMapping::transpose_fwd(const complex<double> *zvec,
 #if USE_MPI
   int status = MPI_Alltoallv((double*)&sbuf[0],&scounts[0],&sdispl[0],
       MPI_DOUBLE,(double*)&rbuf[0],&rcounts[0],&rdispl[0],MPI_DOUBLE,
-      ctxt_.comm());
+      basis_.comm());
   if ( status != 0 )
   {
     cout << " BasisMapping: status = " << status << endl;
-    ctxt_.abort(2);
+    MPI_Abort(basis_.comm(),2);
   }
 #else
   assert(sbuf.size()==rbuf.size());
@@ -513,7 +509,7 @@ void BasisMapping::transpose_bwd(const complex<double> *ct,
 #if USE_MPI
   int status = MPI_Alltoallv((double*)&rbuf[0],&rcounts[0],&rdispl[0],
       MPI_DOUBLE,(double*)&sbuf[0],&scounts[0],&sdispl[0],MPI_DOUBLE,
-      ctxt_.comm());
+      basis_.comm());
   assert ( status == 0 );
 #else
   assert(sbuf.size()==rbuf.size());

src/BasisMapping.h

@@ -15,7 +15,6 @@
 //  BasisMapping.h
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: BasisMapping.h,v 1.4 2008-09-08 15:56:18 fgygi Exp $
 
 #ifndef BASISMAPPING_H
 #define BASISMAPPING_H
@@ -24,13 +23,11 @@
 #include <vector>
 
 class Basis;
-class Context;
 
 class BasisMapping
 {
   private:
 
-  const Context& ctxt_;
   const Basis& basis_;
   int nprocs_, myproc_;
 
@@ -58,8 +55,6 @@ class BasisMapping
   int nvec(void) const { return nvec_; }
   int zvec_size(void) const { return nvec_ * np2_; }
 
-  const Context& context(void) const { return ctxt_; }
-
   // map a function c(G) to zvec_
   void vector_to_zvec(const std::complex<double> *c,
                       std::complex<double> *zvec);

src/Bisection.C

@@ -110,6 +110,9 @@ Bisection::Bisection(const SlaterDet& sd, int nlevels[3])
       if ( np_[idim] % base != 0 ) np_[idim] += base/2;
     }
   }
+  while (!sd.basis().factorizable(np_[0])) np_[0] += (1<<nlevels[0]);
+  while (!sd.basis().factorizable(np_[1])) np_[1] += (1<<nlevels[1]);
+  while (!sd.basis().factorizable(np_[2])) np_[2] += (1<<nlevels[2]);
 
   // number of grid points of augmented grid for normalization
   ft_ = new FourierTransform(sd.basis(),np_[0],np_[1],np_[2]);
@@ -428,9 +431,9 @@ void Bisection::compute_transform(const SlaterDet& sd, int maxsweep, double tol)
     double time = (*i).second.real();
     double tmin = time;
     double tmax = time;
-    wf.context().dmin(1,1,&tmin,1);
-    wf.context().dmax(1,1,&tmax,1);
-    if ( wf.context().myproc()==0 )
+    gcontext_.dmin(1,1,&tmin,1);
+    gcontext_.dmax(1,1,&tmax,1);
+    if ( gcontext_.myproc()==0 )
     {
       cout << "<timing name=\""
            << setw(15) << (*i).first << "\""
@@ -672,6 +675,7 @@ void Bisection::trim_amat(const vector<double>& occ)
   // set to zero the matrix elements of the matrices amat_[k] if they couple
   // states with differing occupation numbers
 
+  const double trim_tol = 1.e-6;