diff --git a/src/PCGSolver.cc b/src/PCGSolver.cc index 7a5fd157..254b1a23 100644 --- a/src/PCGSolver.cc +++ b/src/PCGSolver.cc @@ -176,8 +176,15 @@ bool PCGSolver::solve( /* compute Ax */ oper_.apply(gf_phi, lhs); // set r = b - pb::GridFunc res(gf_rhs); - oper_.transform(res); + pb::GridFunc rhs(gf_rhs); + + // transform r.h.s. to account for dielectric model + oper_.transform(rhs); + + // apply Mehrstelllen r.h.s. if appropriate + pb::GridFunc res(finegrid, bc_[0], bc_[1], bc_[2]); + oper_.rhs(rhs, res); + // compute r = r - Ax res -= lhs; diff --git a/src/pb/Mgm.h b/src/pb/Mgm.h index 4b5ad68d..a821261a 100644 --- a/src/pb/Mgm.h +++ b/src/pb/Mgm.h @@ -35,12 +35,10 @@ bool Mgm(T1& A, T2& vh, const GridFunc& rho, const short cogr, // Compute r.h.s. from rho GridFunc res(rho); - A.transform(res); + A.transform(res); // to account for dielectric model GridFunc rhs(finegrid, bcx, bcy, bcz); - A.rhs(res, rhs); + A.rhs(res, rhs); // to account for possible Mehrstellen operator - // Hartree units - // work GridFunc GridFunc lhs(finegrid, bcx, bcy, bcz); short bcwork[3] = { bcx, bcy, bcz }; @@ -54,7 +52,6 @@ bool Mgm(T1& A, T2& vh, const GridFunc& 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); diff --git a/tests/Chebyshev/test.py b/tests/Chebyshev/test.py index 5cb9a284..28718ee3 100644 --- a/tests/Chebyshev/test.py +++ b/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) diff --git a/tests/Cl2_ONCVPSP_LDA/test.py b/tests/Cl2_ONCVPSP_LDA/test.py index e7ff345a..dcb9bff3 100755 --- a/tests/Cl2_ONCVPSP_LDA/test.py +++ b/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: diff --git a/tests/Fatom/test.py b/tests/Fatom/test.py index db960c3a..def0b3a9 100755 --- a/tests/Fatom/test.py +++ b/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) diff --git a/tests/MVP/test.py b/tests/MVP/test.py index 57d9522e..8a55cc5b 100644 --- a/tests/MVP/test.py +++ b/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) diff --git a/tests/ShortSighted/test.py b/tests/ShortSighted/test.py index 68a617aa..25f7290b 100755 --- a/tests/ShortSighted/test.py +++ b/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') diff --git a/tests/SpinO2/test.py b/tests/SpinO2/test.py index 5d5e823c..8d08bfd9 100755 --- a/tests/SpinO2/test.py +++ b/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'##'): diff --git a/tests/SpinO2LDA/test.py b/tests/SpinO2LDA/test.py index 77557ff0..26010c0a 100755 --- a/tests/SpinO2LDA/test.py +++ 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'##'): diff --git a/tests/SpreadPenalty/test.py b/tests/SpreadPenalty/test.py index ea70f074..059f3320 100755 --- a/tests/SpreadPenalty/test.py +++ 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!")