Bug fix in Poisson PCG Mehrstellen r.h.s.

src/PCGSolver.cc

@@ -176,8 +176,15 @@ bool PCGSolver<OperatorType, ScalarDataType, PrecondDataType>::solve(
     /* compute Ax */
     oper_.apply(gf_phi, lhs);
     // set r = b
-    pb::GridFunc<ScalarDataType> res(gf_rhs);
-    oper_.transform(res);
+    pb::GridFunc<ScalarDataType> rhs(gf_rhs);
+
+    // transform r.h.s. to account for dielectric model
+    oper_.transform(rhs);
+
+    // apply Mehrstelllen r.h.s. if appropriate
+    pb::GridFunc<ScalarDataType> res(finegrid, bc_[0], bc_[1], bc_[2]);
+    oper_.rhs(rhs, res);
+
     // compute r = r - Ax
     res -= lhs;
 

src/pb/Mgm.h

@@ -35,12 +35,10 @@ bool Mgm(T1& A, T2& vh, const GridFunc<T3>& rho, const short cogr,
 
     // Compute r.h.s. from rho
     GridFunc<T3> res(rho);
-    A.transform(res);
+    A.transform(res); // to account for dielectric model
     GridFunc<T3> rhs(finegrid, bcx, bcy, bcz);
-    A.rhs(res, rhs);
+    A.rhs(res, rhs); // to account for possible Mehrstellen operator
 
-    // Hartree units
-    // work GridFunc<T3>
     GridFunc<T3> lhs(finegrid, bcx, bcy, bcz);
 
     short bcwork[3] = { bcx, bcy, bcz };
@@ -54,7 +52,6 @@ bool Mgm(T1& A, T2& vh, const GridFunc<T3>& rho, const short cogr,
     nb_sweeps                 = 0;
     for (short i = 0; i < max_sweeps; i++)
     {
-
         A.apply(vh, lhs);
         // res=rhs-lhs;
         res.diff(rhs, lhs);

tests/Chebyshev/test.py

@@ -70,21 +70,24 @@
   print("Chebyshev-MVP test FAILED for taking too many iterations")
   sys.exit(1)
 
+energy_ref = -108.264614049136
+tol = 1.e-6
 print("Check energy...")
 last_energy = eval(energies[-1])
 print("Energy = {}".format(last_energy))
-if last_energy>-108.0868:
+if abs(last_energy-energy_ref)>tol:
   print("Last energy = {}".format(last_energy))
   sys.exit(1)
 
 tol = 1.e-5
+entropy_ref = 0.17653976
 print("Check entropy...")
 for line in lines:
   if line.count(b'-TS'):
     words=line.split()
     entropy = eval(words[3])
     print("Entropy = {}".format(entropy))
-    entropy_diff = entropy+0.17819
+    entropy_diff = entropy+entropy_ref
     if abs(entropy_diff)>tol:
         print("Check entropy test FAILED. Entropy difference = {}".format(abs(entropy_diff)))
         sys.exit(1)

tests/Cl2_ONCVPSP_LDA/test.py

@@ -36,7 +36,7 @@
 lines=output.split(b'\n')
 
 tol = 4.e-6
-Fz  = 1.2e-3
+Fz  = -7.33e-04
 for line in lines:
   num_matches = line.count(b'%%')
   if num_matches:
@@ -51,7 +51,7 @@
       for i in range(5,7):
         force = eval(words[i])
         if abs(force)>tol:
-          print("force = {}".format(force))
+          print("Force larger than tol, force = {}".format(force))
           sys.exit(1)
       #check value of force in z direction
       if abs(eval(words[7])-Fz)>2.e-5:

tests/Fatom/test.py

@@ -65,7 +65,7 @@
       print("ERROR Eigenvalue 0 = {}".format(eval(eigenvalues[0])))
       sys.exit(1)
     for ii in range(3):
-      if abs(eval(eigenvalues[1+ii])+0.409)>tole:
+      if abs(eval(eigenvalues[1+ii])+0.410)>tole:
         print("ERROR Eigenvalue {} = {}".format(1+ii,eval(eigenvalues[1+ii])))
         sys.exit(1)
     sys.exit(0)

tests/MVP/test.py

@@ -87,12 +87,13 @@
 
 print(eigenvalues)
 tol = 1.e-4
-eigenvalue0 = -0.208
+eigenvalue0 = -0.210
 if abs(eigenvalues[0]-eigenvalue0)>tol:
   print("Expected eigenvalue 0 to be {}".format(eigenvalue0))
   sys.exit(1)
-eigenvalue50 = 0.208
+eigenvalue50 = 0.205
 if abs(eigenvalues[50]-eigenvalue50)>tol:
+  print("Eeigenvalue 50 = {}".format(eigenvalues[50]))
   print("Expected eigenvalue 50 to be {}".format(eigenvalue50))
   sys.exit(1)
 

tests/ShortSighted/test.py

@@ -67,8 +67,8 @@
 print("Check energies...")
 tol = 1.e-3
 count = 0
-energy1_ref = -83.904
-energy4_ref = -83.871
+energy1_ref = -83.929
+energy4_ref = -83.896
 
 for line in lines:
   num_matches1 = line.count(b'IONIC')
@@ -93,7 +93,7 @@
 tol = 1.e-1
 count = 0
 temperature1_ref = 948.253
-temperature4_ref = 916.029
+temperature4_ref = 916.336
 
 for line in lines:
   num_matches1 = line.count(b'Kinetic')

tests/SpinO2/test.py

@@ -43,15 +43,15 @@
     words=line.split()
     energy = eval(words[5][:-1])
 
-ref_energy = -31.805
+ref_energy = -31.808
 print("energy = {}".format(energy))
 if abs(ref_energy-energy) > 1.e-3:
-  print("Incorrect energy!")
+  print("Expected energy = {}".format(ref_energy))
   sys.exit(1)
 
 #make sure forces are below tolerance
 tol = 6.e-4
-Fz  = -1.06e-2
+Fz  = -0.96e-2
 for line in lines:
   #find output lines with forces
   if line.count(b'##'):

tests/SpinO2LDA/test.py

@@ -43,15 +43,15 @@
     words=line.split()
     energy = eval(words[5][:-1])
 
-ref_energy = -31.6105
+ref_energy = -31.6130
 print("energy = {}".format(energy))
 if abs(ref_energy-energy) > 1.e-3:
-  print("Incorrect energy!")
+  print("Incorrect energy, expected {}".format(ref_energy))
   sys.exit(1)
 
 #make sure forces are below tolerance
 tol = 4.e-4
-Fz  = 1.e-2
+Fz  = 1.06e-2
 for line in lines:
   #find output lines with forces
   if line.count(b'##'):

tests/SpreadPenalty/test.py

@@ -79,7 +79,7 @@
 
 #we tolerate an energy difference since the initial wave functions
 #are very delocalized and the spread penalty remains active all along
-energy_ref = -17.16448
+energy_ref = -17.1660
 tol = 5.e-4
 if abs(energy-energy_ref) > tol:
   print("Test failed: last energy value incorrect!")