Merge branch 'develop'

src/LineMinimizer.cpp

@@ -179,8 +179,7 @@ double LineMinimizer::next_alpha(double alpha, double f, double fp)
     if ( debug_print )
       cout << "LineMinimizer: U2, increase alpha" << endl;
 
-    const double delta = 1.1;
-    double new_alpha = std::min(alpha+delta*(alpha-alpha_low), alpha_max_);
+    double new_alpha = std::min(alpha+delta_*(alpha-alpha_low), alpha_max_);
     if ( new_alpha == alpha_max_ )
       done_ = true;
     return new_alpha;

src/LineMinimizer.h

@@ -27,7 +27,7 @@ class LineMinimizer
          alpha_m,alpha_low,alpha_high;
   bool first_use, done_, fail_, bracketing, use_psi;
   bool debug_print;
-  double alpha_start_, alpha_max_, sigma1_, sigma2_;
+  double alpha_start_, alpha_max_, sigma1_, sigma2_, delta_;
   int nstep_, nstep_max_;
 
   double psi(double alpha, double f) { return f - f0 - alpha * fp0 * sigma1_; }
@@ -37,7 +37,7 @@ class LineMinimizer
   public:
 
   LineMinimizer(void) : sigma1_(0.1), sigma2_(0.5), alpha_start_(0.1),
-   alpha_max_(1.0), first_use(true), done_(false), fail_(false),
+   alpha_max_(1.0), delta_(1.1), first_use(true), done_(false), fail_(false),
    bracketing(false), use_psi(true), nstep_(0), nstep_max_(5),
    debug_print(false) {}
   void reset(void) { first_use = true; done_ = false; fail_ = false;
@@ -53,6 +53,7 @@ class LineMinimizer
   void set_sigma2(double s) { sigma2_ = s; }
   void set_alpha_start(double a) { alpha_start_ = a; }
   void set_alpha_max(double a) { alpha_max_ = a; }
+  void set_delta(double d) { delta_ = d; }
   void set_nstep_max(int n) { nstep_max_ = n; }
   void set_debug_print(void) { debug_print = true; }
 

src/MDWavefunctionStepper.cpp

@@ -280,7 +280,9 @@ double MDWavefunctionStepper::ekin_eh(void)
       }
     }
   }
-  wf_.sd_context().dsum(1,1,&ekin_e,1);
+  double tsum;
+  MPI_Allreduce(&ekin_e,&tsum,1,MPI_DOUBLE,MPI_SUM,MPIdata::comm());
+  ekin_e = tsum;
   tmap_["ekin_e"].stop();
   return ekin_e;
 }

src/RSHFunctional.cpp

@@ -655,10 +655,8 @@ void RSHFunctional::RSH_exchange(const double rho, const double grad,
   RSH_enhance(s,kF,w,&fxhse,&dfx_ds,&dfx_dkf);
 
   // calculate exchange energy
-  // ex = (1 - a) ex,SR + ex,LR
-  //    = (1 - a) ex,SR + ex,PBE - ex,SR
-  //    = ex,PBE - a ex,SR
-  *ex = exLDA * ( fxpbe - a_ex * fxhse );
+  *ex = exLDA * ( ( 1.0 - alpha_RSH_ ) * fxpbe +
+                  ( alpha_RSH_ - beta_RSH_ ) * fxhse );
 
   // calculate potential
   *vx1 = third4 * exLDA * ( fxpbe - s2 * fs - a_ex * ( fxhse - s * dfx_ds

src/release.cpp

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

util/qbox_vavg.py

+#!/usr/bin/env python3
+#
+# qbox_vavg.py:
+# use: qbox_vavg.py file.cube
+# generate plot data of averaged local potential in x,y,z directions
+# read data from cube file generated by Qbox with the plot -vlocal command
+# write plot data on standard output in gnuplot format
+# current version works in orthorhombic cells only
+
+import numpy as np
+import sys
+
+def usage():
+  print("use: ",sys.argv[0]," cube_file")
+  sys.exit()
+
+argc=len(sys.argv)
+if ( argc != 2 ):
+  usage()
+
+filename = sys.argv[1]
+
+f = open(filename, 'r')
+lines = f.readlines()
+
+# number of atoms and origin on third line
+buf = lines[2].split()
+nat = int(buf[0])
+origin = [ float(buf[1]), float(buf[2]), float(buf[3]) ]
+#print("nat = ", nat)
+#print("origin = ", origin)
+
+# read grid parameters
+buf = lines[3].split()
+np0 = int(buf[0])
+da0 = [ float(buf[1]), float(buf[2]), float(buf[3]) ]
+#print(np0,da0)
+buf = lines[4].split()
+np1 = int(buf[0])
+da1 = [ float(buf[1]), float(buf[2]), float(buf[3]) ]
+#print(np1,da1)
+buf = lines[5].split()
+np2 = int(buf[0])
+da2 = [ float(buf[1]), float(buf[2]), float(buf[3]) ]
+#print(np2,da2)
+
+# check that the cell is orthorhombic
+assert da0[1] == 0.0 and da0[2] == 0
+assert da1[0] == 0.0 and da1[2] == 0
+assert da2[0] == 0.0 and da2[1] == 0
+
+# read vlocal data starting at line 6+nat
+v = np.empty(0)
+for line in lines[6+nat:]:
+  vals = ([float(val) for val in line.split()])
+  v = np.append(v,vals)
+
+# fastest increasing index in cube file is z
+v = v.reshape(np0,np1,np2)
+
+vx = np.sum(v,(1,2))/np0
+vy = np.sum(v,(0,2))/np1
+vz = np.sum(v,(0,1))/np2
+
+#print("Vavg(x) min/max = ", '%8f %8f' %(min(vx), max(vx)))
+#print("Vavg(y) min/max = ", '%8f %8f' %(min(vy), max(vy)))
+#print("Vavg(z) min/max = ", '%8f %8f' %(min(vz), max(vz)))
+
+# output position in bohr units
+x = [origin[0]+da0[0]*i for i in range(np0)]
+y = [origin[1]+da1[1]*i for i in range(np1)]
+z = [origin[2]+da2[2]*i for i in range(np2)]
+
+# print average potential on stdout in gnuplot format
+print("# Vavg(x)")
+for i in range(np0):
+  print('%.8f %.8f' %(x[i], vx[i]))
+
+print("\n\n# Vavg(y)")
+for i in range(np1):
+  print('%.8f %.8f' %(y[i], vy[i]))
+
+print("\n\n# Vavg(z)")
+for i in range(np2):
+  print('%.8f %.8f' %(z[i], vz[i]))