From 275bad89e07529c395e846321a7dc5946bf8c76e Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 20 Oct 2023 09:46:50 +0200 Subject: [PATCH 01/28] Flags can be set from input --- test.py | 38 ++++++++++++++++++++++++++++++++------ 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/test.py b/test.py index 954b41c..37a6ed5 100644 --- a/test.py +++ b/test.py @@ -21,6 +21,16 @@ import importlib importlib.reload(orb) +def read_step_file(filename): + steps = {} + with open(filename, 'r') as file: + for line in file: + terms = line.strip().split() + option = terms[0] + value = terms[1] + steps[option] = value + return steps + if __name__ == '__main__': parser = argparse.ArgumentParser(description='Collecting all data tostart the program.') parser.add_argument('-d', '--derivative', dest='deriv', type=str, default='ABGV', @@ -62,30 +72,45 @@ derivative = args.deriv print('call MRA DONE') + steplist = 'step_list.txt' + atomlist = 'atom_list.txt' + steps = read_step_file(steplist) + coordinates, number = nucpot.read_file_with_named_lists(atomlist) + ################## Jobs ########################## - computeNuclearPotential = True - readOrbitals = True + computeNuclearPotential = False + readOrbitals = False + readPotential = False runD_1e = False runD2_1e = False runCoulombGen = False - runCoulomb2e = False #Use it + runCoulomb2e = False runKutzelnigg = False - runKutzSimple = True #I will use it + runKutzSimple = False saveOrbitals = False + savePotential = False runGaunt = False runGaugeA = False runGaugeB = False runGaugeC = False runGaugeD = False runGaugeDelta = False - print('Jobs chosen') + + for key in steps: + flag = (steps[key] == 'True') + locals()[key] = flag + print(key, steps[key], locals()[key]) + + if(computeNuclearPotential): print("hello") ################### Reading Atoms ######################### - atomlist = 'atom_list.txt' # Replace with the actual file name + steplist = 'step_list.txt' + atomlist = 'atom_list.txt' coordinates, number = nucpot.read_file_with_named_lists(atomlist) print("Number of Atoms = ", number) print(coordinates) + exit(-1) ################### Define V potential ###################### V_tree = vp.FunctionTree(mra) @@ -176,3 +201,4 @@ if saveOrbitals: spinorb1.save("spinorb1") #spinorb2.save("spinorb2") + From c3119fd4926d4e09dc7efe396ce984265bcf632b Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 27 Oct 2023 18:16:23 +0200 Subject: [PATCH 02/28] Update to gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 35acbfe..19a55ce 100644 --- a/.gitignore +++ b/.gitignore @@ -136,3 +136,5 @@ dmypy.json *.out *.prof *.txt +*.err +*.tree From 3bc4818135634611340fb5920eb7f7619b554adc Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 27 Oct 2023 18:16:57 +0200 Subject: [PATCH 03/28] 1. Some output improvements 2. Modulation of precision for the D2 method --- one_electron.py | 93 +++++++++++++++++++++++++++++++------------------ 1 file changed, 59 insertions(+), 34 deletions(-) diff --git a/one_electron.py b/one_electron.py index 5987ddc..ef1b7d4 100644 --- a/one_electron.py +++ b/one_electron.py @@ -17,94 +17,119 @@ def analytic_1s(light_speed, n, k, Z): tmp3 = 1 + tmp2**2 return light_speed**2 / np.sqrt(tmp3) -def gs_D_1e(spinorb1, potential, mra, prec, derivative): - print('Hartree-Fock 1e') +def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative): + print('One-electron calculations') error_norm = 1 #compute_last_energy = False light_speed = spinorb1.light_speed - - while error_norm > prec: + old_energy = 0 + + while error_norm > thr: hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi energy = spinorb1.dot(add_psi).real - print('Energy',energy - light_speed**2) mu = orb.calc_dirac_mu(energy, light_speed) tmp = orb.apply_helmholtz(v_psi, mu, prec) - tmp.crop(prec/10) +# tmp = orb.apply_dirac_hamiltonian(v_psi, prec, energy, der = derivative) + tmp.cropLargeSmall(prec) new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = derivative) - new_orbital.crop(prec/10) +# new_orbital = orb.apply_helmholtz(tmp, mu, prec) +# new_orbital.cropLargeSmall(prec) + new_orbital.cropLargeSmall(prec) new_orbital.normalize() delta_psi = new_orbital - spinorb1 #orbital_error = delta_psi.dot(delta_psi).real deltasq = delta_psi.squaredNorm() error_norm = np.sqrt(deltasq) print('Error', error_norm) + print('Delta E', energy - old_energy) + print('Energy',energy - light_speed**2) + old_energy = energy spinorb1 = new_orbital hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi energy = spinorb1.dot(add_psi).real - print('Final Energy',energy - light_speed**2) + energy_1s = analytic_1s(light_speed, 1, -1, 1) + print("Exact energy: ", energy_1s - light_speed**2) + print('Final Energy:',energy - light_speed**2) + print('Delta Energy:',energy - old_energy) + print('Error Energy:',energy - energy_1s) return spinorb1 -def gs_D2_1e(spinorb1, potential, mra, prec, derivative): +def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative): print('Hartree-Fock 1e D2') error_norm = 1 light_speed = spinorb1.light_speed c2 = light_speed * light_speed - - while error_norm > prec: + old_energy = 0 + while error_norm > thr: v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) - vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec) + vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() apV_psi = v_psi.alpha_p(prec, derivative) - ap_psi = spinorb1.alpha_p(prec, derivative) - Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec) + ap_psi = spinorb1.alpha_p(prec*light_speed, derivative) + Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec*light_speed) anticom = apV_psi + Vap_psi - anticom.cropLargeSmall(prec) - #beta_v_psi.cropLargeSmall(prec) - vv_psi.cropLargeSmall(prec) +# anticom.cropLargeSmall(prec) +# beta_v_psi.cropLargeSmall(prec) +# vv_psi.cropLargeSmall(prec) RHS = beta_v_psi + vv_psi + anticom * (0.5/light_speed) - #RHS.cropLargeSmall(prec) + RHS.cropLargeSmall(prec) cke = spinorb1.classicT() cpe = (spinorb1.dot(RHS)).real print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", cke + cpe) + classic_energy = cke + cpe + energy = c2*(np.sqrt(1+2*classic_energy/c2)-1) mu = orb.calc_non_rel_mu(cke+cpe) print("mu =", mu) new_orbital = orb.apply_helmholtz(RHS, mu, prec) + new_orbital.cropLargeSmall(prec) new_orbital.normalize() delta_psi = new_orbital - spinorb1 deltasq = delta_psi.squaredNorm() error_norm = np.sqrt(deltasq) print("Error =", error_norm) + print('Delta E', energy - old_energy) + print('Energy',energy, old_energy) + old_energy = energy spinorb1 = new_orbital hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi - energy = spinorb1.dot(add_psi).real + energy_dirac = spinorb1.dot(add_psi).real - cke = spinorb1.classicT() +# v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) + vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() - beta_pot = (beta_v_psi.dot(spinorb1)).real - pot_sq = (v_psi.dot(v_psi)).real - ap_psi = spinorb1.alpha_p(prec, derivative) - anticom = (ap_psi.dot(v_psi)).real - energy_kutzelnigg = cke + beta_pot + pot_sq/(2*c2) + anticom/light_speed - - print('Kutzelnigg =',cke, beta_pot, pot_sq/(2*c2), anticom/light_speed, energy_kutzelnigg) - print('Quadratic approx =',energy_kutzelnigg - energy_kutzelnigg**2/(2*c2)) - print('Correct from Kutzelnigg =', c2*(np.sqrt(1+2*energy_kutzelnigg/c2)-1)) - print('Final Energy =',energy - light_speed**2) - + apV_psi = v_psi.alpha_p(prec*light_speed, derivative) + ap_psi = spinorb1.alpha_p(prec*light_speed, derivative) + Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec) + anticom = apV_psi + Vap_psi +# anticom.cropLargeSmall(prec) +# beta_v_psi.cropLargeSmall(prec) +# vv_psi.cropLargeSmall(prec) + RHS = beta_v_psi + vv_psi + anticom * (0.5/light_speed) + RHS.cropLargeSmall(prec) + cke = spinorb1.classicT() + cpe = (spinorb1.dot(RHS)).real + classic_energy = cke + cpe + energy = c2*(np.sqrt(1+2*classic_energy/c2)-1) energy_1s = analytic_1s(light_speed, 1, -1, 1) - print('Exact Energy =',energy_1s - light_speed**2) - print('Difference 1 =',energy_1s - energy) - print('Difference 2 =',energy_1s - energy_kutzelnigg - light_speed**2) +# print('Kutzelnigg =',cke, cpe, energy_kutzelnigg) +# print('Quadratic approx =',energy_kutzelnigg - energy_kutzelnigg**2/(2*c2)) + print('Exact Energy = ', energy_1s - light_speed**2) + print('Dirac Energy = ', energy_dirac - light_speed**2) + print('Kutze Energy = ', energy) + print('Error Kutze = ', energy - energy_1s + light_speed**2) + print('Error Dirac = ', energy_dirac - energy_1s) + print('Delta Energy = ', energy - old_energy) + print('Dirac - Kutzelnigg = ', energy_dirac - energy - light_speed**2) return spinorb1 From 393ee7e4016e5a1521776b4348df6d152acf1ebf Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 27 Oct 2023 18:20:13 +0200 Subject: [PATCH 04/28] precSmall = precLarge/100 --- orbital4c/orbital.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orbital4c/orbital.py b/orbital4c/orbital.py index f48ae98..0cfdbe5 100644 --- a/orbital4c/orbital.py +++ b/orbital4c/orbital.py @@ -101,7 +101,7 @@ def cropLargeSmall(self, prec): largeNorm = np.sqrt(self.squaredLargeNorm()) smallNorm = np.sqrt(self.squaredSmallNorm()) precLarge = prec * largeNorm - precSmall = prec * smallNorm + precSmall = prec * largeNorm / 100 # print('precisions', precLarge, precSmall) self['La'].crop(precLarge, True) self['Lb'].crop(precLarge, True) From d57cd7bea0256612e47a7c22a17bb79cffbf2932 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 27 Oct 2023 18:21:02 +0200 Subject: [PATCH 05/28] add threshold --- test.py | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/test.py b/test.py index 37a6ed5..faf9d4a 100644 --- a/test.py +++ b/test.py @@ -35,17 +35,19 @@ def read_step_file(filename): parser = argparse.ArgumentParser(description='Collecting all data tostart the program.') parser.add_argument('-d', '--derivative', dest='deriv', type=str, default='ABGV', help='put the type of derivative') - parser.add_argument('-b', '--box', dest='box', type=int, default=60, + parser.add_argument('-b', '--box', dest='box', type=int, default=100, help='put the box dimension') parser.add_argument('-l', '--light_speed', dest='lux_speed', type=float, default=137.03599913900001, help='light of speed') - parser.add_argument('-o', '--order', dest='order', type=int, default=6, + parser.add_argument('-o', '--order', dest='order', type=int, default=10, help='put the order of Polinomial') parser.add_argument('-p', '--prec', dest='prec', type=float, default=1e-4, help='put the precision') + parser.add_argument('-t', '--threshold', dest='thr', type=float, default=1e-4, + help='put the orbital threshold') parser.add_argument('-e', '--coulgau', dest='coulgau', type=str, default='coulomb', help='put the coulomb or gaunt') - parser.add_argument('-v', '--potential', dest='potential', type=str, default='point_charge', + parser.add_argument('-v', '--potential', dest='potential', type=str, default='coulomb_HFYGB', help='tell me wich model for V you want to use point_charge, coulomb_HFYGB, homogeneus_charge_sphere, gaussian') args = parser.parse_args() @@ -66,6 +68,7 @@ def read_step_file(filename): ################# Call MRA ####################### mra = vp.MultiResolutionAnalysis(box=[-args.box,args.box], order=args.order, max_depth = 30) prec = args.prec + thr = args.thr orb.orbital4c.mra = mra orb.orbital4c.light_speed = light_speed cf.complex_fcn.mra = mra @@ -108,9 +111,6 @@ def read_step_file(filename): atomlist = 'atom_list.txt' coordinates, number = nucpot.read_file_with_named_lists(atomlist) - print("Number of Atoms = ", number) - print(coordinates) - exit(-1) ################### Define V potential ###################### V_tree = vp.FunctionTree(mra) @@ -121,14 +121,21 @@ def read_step_file(filename): V_tree = Peps(f) print("Define V", args.potential, "DONE") com_coordinates = nucpot.calculate_center_of_mass(coordinates) - print("Center of Mass (x, y, z):", com_coordinates) + if(savePotential): + V_tree.saveTree(f"potential") + + if(readPotential): + V_tree.loadTree(f"potential") + + print("Number of Atoms = ", number) + print(coordinates) #############################START WITH CALCULATION################################### spinorb1 = orb.orbital4c() spinorb2 = orb.orbital4c() if readOrbitals: spinorb1.read("spinorb1") - #spinorb2.read("spinorb2") + spinorb2 = spinorb1.ktrs(prec) else: gauss_tree_tot = vp.FunctionTree(mra) gauss_tree_tot.setZero() @@ -153,20 +160,16 @@ def read_step_file(filename): spinorb1.init_small_components(prec/10) spinorb1.normalize() spinorb1.cropLargeSmall(prec) - # print('Spin1', spinorb1) spinorb2 = spinorb1.ktrs(prec) - #print("spinorb1") - #print(spinorb1) - length = 2 * args.box print("Using derivative ", derivative) if runD_1e: - spinorb1 = one_electron.gs_D_1e(spinorb1, V_tree, mra, prec, derivative) + spinorb1 = one_electron.gs_D_1e(spinorb1, V_tree, mra, prec, thr, derivative) if runD2_1e: - spinorb1 = one_electron.gs_D2_1e(spinorb1, V_tree, mra, prec, derivative) + spinorb1 = one_electron.gs_D2_1e(spinorb1, V_tree, mra, prec, thr, derivative) if runCoulombGen: spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec, derivative) From 5e0bb384aa7f6ef6ff5c670e4c944e18dd00aa46 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Wed, 24 Jan 2024 11:00:47 +0100 Subject: [PATCH 06/28] The iteration is now interrupted either when the change in orbital norm is below threshold OR when the energy difference is below the threshold/1000 (three orders of magnitude more than the orbital threshold) --- one_electron.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/one_electron.py b/one_electron.py index ef1b7d4..5b5b422 100644 --- a/one_electron.py +++ b/one_electron.py @@ -25,8 +25,8 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative): light_speed = spinorb1.light_speed old_energy = 0 - - while error_norm > thr: + delta_e = 1 + while (error_norm > thr and delta_e > thr/1000): hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi @@ -45,7 +45,8 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative): deltasq = delta_psi.squaredNorm() error_norm = np.sqrt(deltasq) print('Error', error_norm) - print('Delta E', energy - old_energy) + delta_e = np.abs(energy - old_energy) + print('Delta E', delta_e) print('Energy',energy - light_speed**2) old_energy = energy spinorb1 = new_orbital @@ -65,10 +66,11 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative): def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative): print('Hartree-Fock 1e D2') error_norm = 1 + delta_e = 1 light_speed = spinorb1.light_speed c2 = light_speed * light_speed old_energy = 0 - while error_norm > thr: + while (error_norm > thr and delta_e > thr/1000): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() @@ -95,7 +97,8 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative): deltasq = delta_psi.squaredNorm() error_norm = np.sqrt(deltasq) print("Error =", error_norm) - print('Delta E', energy - old_energy) + delta_e = np.abs(energy - old_energy) + print('Delta E', delta_e) print('Energy',energy, old_energy) old_energy = energy spinorb1 = new_orbital From 1c062b3568277a8dc63a0e5ac0fc48b0b1c164a3 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Wed, 24 Jan 2024 11:50:25 +0100 Subject: [PATCH 07/28] Added precision modulation to the 2e case as well --- one_electron.py | 4 ++-- orbital4c/operators.py | 24 ++++++++++++------------ two_electron.py | 12 ++++++------ 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/one_electron.py b/one_electron.py index 5b5b422..27c4d0d 100644 --- a/one_electron.py +++ b/one_electron.py @@ -74,7 +74,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() - apV_psi = v_psi.alpha_p(prec, derivative) + apV_psi = v_psi.alpha_p(prec*light_speed, derivative) ap_psi = spinorb1.alpha_p(prec*light_speed, derivative) Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec*light_speed) anticom = apV_psi + Vap_psi @@ -113,7 +113,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative): beta_v_psi = v_psi.beta2() apV_psi = v_psi.alpha_p(prec*light_speed, derivative) ap_psi = spinorb1.alpha_p(prec*light_speed, derivative) - Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec) + Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec*light_speed) anticom = apV_psi + Vap_psi # anticom.cropLargeSmall(prec) # beta_v_psi.cropLargeSmall(prec) diff --git a/orbital4c/operators.py b/orbital4c/operators.py index d7a6ad3..c09415e 100644 --- a/orbital4c/operators.py +++ b/orbital4c/operators.py @@ -94,11 +94,11 @@ def setup(self): vp.advanced.add(self.prec, rho, rholist) self.potential = (4.0*np.pi) * self.poisson(rho).crop(self.prec) - def __call__(self, phi): + def __call__(self, phi, prec_mod = 1.0): complex_pot = cf.complex_fcn() complex_pot.real = self.potential complex_pot.imag.setZero() - out = orb.apply_complex_potential(1.0, complex_pot, phi, self.prec) + out = orb.apply_complex_potential(1.0, complex_pot, phi, self.prec * prec_mod) out.cropLargeSmall(self.prec) return out @@ -108,20 +108,20 @@ def __init__(self, mra, prec, Psi): self.poisson = vp.PoissonOperator(mra=self.mra, prec=self.prec) self.potential = None - def __call__(self, phi): + def __call__(self, phi, prec_mod = 1): Kij_array = [] coeff_array = [] for i in range(0, len(self.Psi)): V_ij = cf.complex_fcn() - overlap_density = self.Psi[i].overlap_density(phi, self.prec) - V_ij.real = self.poisson(overlap_density.real).crop(self.prec) - V_ij.imag = self.poisson(overlap_density.imag).crop(self.prec) - tmp = orb.apply_complex_potential(1.0, V_ij, self.Psi[i], self.prec) + overlap_density = self.Psi[i].overlap_density(phi, self.prec * prec_mod) + V_ij.real = self.poisson(overlap_density.real).crop(self.prec * prec_mod) + V_ij.imag = self.poisson(overlap_density.imag).crop(self.prec * prec_mod) + tmp = orb.apply_complex_potential(1.0, V_ij, self.Psi[i], self.prec * prec_mod) Kij_array.append(tmp) coeff_array.append(1.0) - output = orb.add_vector(Kij_array, coeff_array, self.prec) + output = orb.add_vector(Kij_array, coeff_array, self.prec * prec_mod) output *= 4.0 * np.pi - output.cropLargeSmall(self.prec) + output.cropLargeSmall(self.prec * prec_mod) return output class PotentialOperator(Operator): @@ -130,11 +130,11 @@ def __init__(self, mra, prec, potential, real = True): self.potential = potential self.real = real - def __call__(self, phi): + def __call__(self, phi, prec_mod): if(self.real): - result = orb.apply_potential(1.0, self.potential, phi, self.prec) + result = orb.apply_potential(1.0, self.potential, phi, self.prec * prec_mod) else: - result = orb.apply_complex_potential(1.0, self.potential, phi, self.prec) + result = orb.apply_complex_potential(1.0, self.potential, phi, self.prec * prec_mod) result.cropLargeSmall(self.prec) return result diff --git a/two_electron.py b/two_electron.py index e95b5e8..7dca6d3 100644 --- a/two_electron.py +++ b/two_electron.py @@ -646,16 +646,16 @@ def build_RHS_D2(Jop, Vop, spinor, prec, light_speed): beta_VT_psi = VT_psi.beta2() beta_VT_psi.cropLargeSmall(prec) - ap_VT_psi = VT_psi.alpha_p(prec) - ap_psi = spinor.alpha_p(prec) - VT_ap_psi = 0.5 * Jop(ap_psi) - Vop(ap_psi) + ap_VT_psi = VT_psi.alpha_p(prec*light_speed) + ap_psi = spinor.alpha_p(prec*light_speed) + VT_ap_psi = 0.5 * Jop(ap_psi,light_speed) - Vop(ap_psi,light_speed) anticom = VT_ap_psi + ap_VT_psi anticom *= 1.0 / (2.0 * light_speed) - anticom.cropLargeSmall(prec) + anticom.cropLargeSmall(prec*light_speed) - VT_VT_psi = 0.5 * Jop(VT_psi) - Vop(VT_psi) + VT_VT_psi = 0.5 * Jop(VT_psi*c2) - Vop(VT_psi*c2) VT_VT_psi *= 1.0 / (2.0 * c2) - VT_VT_psi.cropLargeSmall(prec) + VT_VT_psi.cropLargeSmall(prec*c2) RHS = beta_VT_psi + anticom + VT_VT_psi RHS.cropLargeSmall(prec) From 49ebf67359b63935055d9392eb83ece349cd7031 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Wed, 24 Jan 2024 16:14:20 +0100 Subject: [PATCH 08/28] moved two function definitions out of a third function --- orbital4c/nuclear_potential.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/orbital4c/nuclear_potential.py b/orbital4c/nuclear_potential.py index 02ac967..787e51e 100644 --- a/orbital4c/nuclear_potential.py +++ b/orbital4c/nuclear_potential.py @@ -116,17 +116,19 @@ def point_charge(position, center , charge): distance = np.sqrt(d2) return charge / distance +def smoothing_HFYGB(charge, prec): + factor = 0.00435 * prec / charge**5 + return factor**(1./3.) + +def uHFYGB(r): + u = erf(r)/r + (1/(3*np.sqrt(np.pi)))*(np.exp(-(r**2)) + 16*np.exp(-4*r**2)) + return u + def coulomb_HFYGB(position, center, charge, precision): d2 = ((position[0] - center[0])**2 + (position[1] - center[1])**2 + (position[2] - center[2])**2) distance = np.sqrt(d2) - def smoothing_HFYGB(charge, prec): - factor = 0.00435 * prec / charge**5 - return factor**(1./3.) - def uHFYGB(r): - u = erf(r)/r + (1/(3*np.sqrt(np.pi)))*(np.exp(-(r**2)) + 16*np.exp(-4*r**2)) - return u factor = smoothing_HFYGB(charge, precision) value = uHFYGB(distance/factor) return charge * value / factor From 5520c64ec73bd5da91abf06571c790a42148cfb8 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Wed, 24 Jan 2024 16:19:11 +0100 Subject: [PATCH 09/28] Adds precision modulation for the D2 method for 2 electrons --- orbital4c/operators.py | 2 +- test.py | 10 ++++++---- two_electron.py | 37 ++++++++++++++++++++++++------------- 3 files changed, 31 insertions(+), 18 deletions(-) diff --git a/orbital4c/operators.py b/orbital4c/operators.py index c09415e..c2690e1 100644 --- a/orbital4c/operators.py +++ b/orbital4c/operators.py @@ -130,7 +130,7 @@ def __init__(self, mra, prec, potential, real = True): self.potential = potential self.real = real - def __call__(self, phi, prec_mod): + def __call__(self, phi, prec_mod = 1): if(self.real): result = orb.apply_potential(1.0, self.potential, phi, self.prec * prec_mod) else: diff --git a/test.py b/test.py index faf9d4a..4a17e0f 100644 --- a/test.py +++ b/test.py @@ -172,16 +172,18 @@ def read_step_file(filename): spinorb1 = one_electron.gs_D2_1e(spinorb1, V_tree, mra, prec, thr, derivative) if runCoulombGen: - spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec, derivative) + print("NOT ACTIVE") +# spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec, derivative) if runCoulomb2e: - spinorb1, spinorb2 = two_electron.coulomb_gs_2e(spinorb1, V_tree, mra, prec, derivative) + spinorb1, spinorb2 = two_electron.coulomb_gs_2e(spinorb1, V_tree, mra, prec, thr, derivative) if runKutzelnigg: - spinorb1, spinorb2 = two_electron.coulomb_2e_D2([spinorb1, spinorb2], V_tree, mra, prec, derivative) + print("NOT ACTIVE") +# spinorb1, spinorb2 = two_electron.coulomb_2e_D2([spinorb1, spinorb2], V_tree, mra, prec, derivative) if runKutzSimple: - spinorb1, spinorb2 = two_electron.coulomb_2e_D2_J([spinorb1, spinorb2], V_tree, mra, prec, derivative) + spinorb1, spinorb2 = two_electron.coulomb_2e_D2_J([spinorb1, spinorb2], V_tree, mra, prec, thr, derivative) if runGaunt: two_electron.calcGauntPert(spinorb1, spinorb2, mra, prec) diff --git a/two_electron.py b/two_electron.py index 7dca6d3..6905d11 100644 --- a/two_electron.py +++ b/two_electron.py @@ -188,58 +188,63 @@ def coulomb_2e_D2(spinors, potential, mra, prec, derivative): print("total energy: ", total_energy) return spinors[0], spinors[1] -def coulomb_2e_D2_J(spinors, potential, mra, prec, derivative): +def coulomb_2e_D2_J(spinors, potential, mra, prec, thr, derivative): print('Hartree-Fock (Coulomb interaction) 2e D2 J only') error_norm = 1.0 + delta_e = 1 + old_energy = 0 compute_last_energy = False P = vp.PoissonOperator(mra, prec) light_speed = spinors[0].light_speed c2 = light_speed**2 Vop = oper.PotentialOperator(mra, prec, potential) - while(error_norm > prec): + while(error_norm > thr and delta_e > thr/1000): Jop = oper.CoulombDirectOperator(mra, prec, spinors) RHS = build_RHS_D2(Jop, Vop, spinors[0], prec, light_speed) cke = spinors[0].classicT() cpe = (spinors[0].dot(RHS)).real print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", cke + cpe) print("Orbital energy: ", c2 * ( -1.0 + np.sqrt(1 + 2 * (cpe + cke) / c2))) + class_e = cke + cpe + delta_e = np.abs(class_e - old_energy) mu = orb.calc_non_rel_mu(cke+cpe) + print("mu = ", mu) new_spinor = orb.apply_helmholtz(RHS, mu, prec) - #print("============= Spinor before Helmholtz =============") - #print(spinors[0]) - #print("============= RHS before Helmh#oltz =============") - #print(RHS) - #print("============= New spinor before crop =============") - #print(new_spinor) new_spinor.cropLargeSmall(prec) new_spinor.normalize() delta_psi = new_spinor - spinors[0] deltasq = delta_psi.squaredNorm() error_norm = np.sqrt(deltasq) print('Orbital_Error norm', error_norm) + print('delta_e', delta_e) spinors[0] = new_spinor spinors[1] = new_spinor.ktrs(prec) + old_energy = class_e Jop = oper.CoulombDirectOperator(mra, prec, spinors) RHS = build_RHS_D2(Jop, Vop, spinors[0], prec, light_speed) cke = spinors[0].classicT() cpe = (spinors[0].dot(RHS)).real + class_e = cke + cpe + delta_e = np.abs(class_e - old_energy) final_orbital_energy = c2 * ( -1.0 + np.sqrt(1 + 2 * (cpe + cke) / c2)) Jorb = Jop(spinors[0]) Jenergy = (spinors[0].dot(Jorb)).real final_total_energy = 2.0 * final_orbital_energy - 0.5 * Jenergy - print("Final classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", cke + cpe) + print("Final classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", cke + cpe, "delta_e", delta_e) print("Final orbital energy: ", final_orbital_energy) print("Final Total energy: ", final_total_energy) return spinors[0], spinors[1] -def coulomb_gs_2e(spinorb1, potential, mra, prec, derivative): +def coulomb_gs_2e(spinorb1, potential, mra, prec, thr, derivative): print('Hartree-Fock (Coulomb interaction)') error_norm = 1 + delta_e = 1 compute_last_energy = False + old_energy = 0 P = vp.PoissonOperator(mra, prec) light_speed = spinorb1.light_speed # for i in range(10): - while (error_norm > prec or compute_last_energy): + while ((error_norm > thr and delta_e > thr/1000) or compute_last_energy): n_22 = spinorb1.overlap_density(spinorb1, prec) # Definition of two electron operators @@ -266,8 +271,13 @@ def coulomb_gs_2e(spinorb1, potential, mra, prec, derivative): eps = hd_V_11.real + JmK.real E_tot_JK = 2*eps - JmK.real + orbital_e = eps - light_speed**2 + total_e = E_tot_JK - (2.0 *light_speed**2) + delta_e = np.abs(total_e - old_energy) print('orbital energy', eps - light_speed**2) print('total energy', E_tot_JK - (2.0 *light_speed**2)) + print('delta_e', delta_e) + if(compute_last_energy): break @@ -293,8 +303,9 @@ def coulomb_gs_2e(spinorb1, potential, mra, prec, derivative): print('Orbital_Error norm', error_norm) spinorb1 = new_orbital spinorb2 = spinorb1.ktrs(prec) - if(error_norm < prec): + if(error_norm < thr and delta_e < thr/1000): compute_last_energy = True + old_energy = total_e return spinorb1, spinorb2 #def coulomb_gs(potential, spinors, mra, prec, der = 'ABGV'): @@ -653,7 +664,7 @@ def build_RHS_D2(Jop, Vop, spinor, prec, light_speed): anticom *= 1.0 / (2.0 * light_speed) anticom.cropLargeSmall(prec*light_speed) - VT_VT_psi = 0.5 * Jop(VT_psi*c2) - Vop(VT_psi*c2) + VT_VT_psi = 0.5 * Jop(VT_psi,c2) - Vop(VT_psi,c2) VT_VT_psi *= 1.0 / (2.0 * c2) VT_VT_psi.cropLargeSmall(prec*c2) From 2b44fa66b9a80ddd0d25d01ac95af5db19287ecf Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Thu, 25 Jan 2024 10:45:03 +0100 Subject: [PATCH 10/28] Improves energy priting D2 two electrons --- two_electron.py | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/two_electron.py b/two_electron.py index 6905d11..8b60963 100644 --- a/two_electron.py +++ b/two_electron.py @@ -203,10 +203,11 @@ def coulomb_2e_D2_J(spinors, potential, mra, prec, thr, derivative): RHS = build_RHS_D2(Jop, Vop, spinors[0], prec, light_speed) cke = spinors[0].classicT() cpe = (spinors[0].dot(RHS)).real - print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", cke + cpe) - print("Orbital energy: ", c2 * ( -1.0 + np.sqrt(1 + 2 * (cpe + cke) / c2))) - class_e = cke + cpe - delta_e = np.abs(class_e - old_energy) + classical_energy = cke + cpe + orbital_energy = c2 * ( -1.0 + np.sqrt(1 + 2 * (classical_energy) / c2)) + print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", classical_energy) + print("Orbital energy: ", orbital_energy) + delta_e = np.abs(orbital_energy - old_energy) mu = orb.calc_non_rel_mu(cke+cpe) print("mu = ", mu) new_spinor = orb.apply_helmholtz(RHS, mu, prec) @@ -219,20 +220,28 @@ def coulomb_2e_D2_J(spinors, potential, mra, prec, thr, derivative): print('delta_e', delta_e) spinors[0] = new_spinor spinors[1] = new_spinor.ktrs(prec) - old_energy = class_e + old_energy = orbital_energy Jop = oper.CoulombDirectOperator(mra, prec, spinors) RHS = build_RHS_D2(Jop, Vop, spinors[0], prec, light_speed) cke = spinors[0].classicT() cpe = (spinors[0].dot(RHS)).real - class_e = cke + cpe - delta_e = np.abs(class_e - old_energy) - final_orbital_energy = c2 * ( -1.0 + np.sqrt(1 + 2 * (cpe + cke) / c2)) + classical_energy = cke + cpe + orbital_energy = c2 * ( -1.0 + np.sqrt(1 + 2 * (classical_energy) / c2)) + delta_e = np.abs(orbital_energy - old_energy) Jorb = Jop(spinors[0]) + hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) + Vorb = Vop(spinors[0]) Jenergy = (spinors[0].dot(Jorb)).real - final_total_energy = 2.0 * final_orbital_energy - 0.5 * Jenergy + hdenergy = (spinors[0].dot(hd_psi)).real + Venergy = -(spinors[0].dot(Vorb)).real + orbital_energy_dirac = hdenergy + Jenergy + Venergy + total_energy = 2.0 * orbital_energy - 0.5 * Jenergy print("Final classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", cke + cpe, "delta_e", delta_e) - print("Final orbital energy: ", final_orbital_energy) - print("Final Total energy: ", final_total_energy) + print("Final orbital energy: ", orbital_energy) + print("Final orbital energy Dirac: ", orbital_energy_dirac) + print("Final Total energy: ", total_energy) + print("Final orbital energy difference: ", orbital_energy - old_energy) + print("Dirac - Kutzelnigg", orbital_energy - orbital_energy_dirac) return spinors[0], spinors[1] def coulomb_gs_2e(spinorb1, potential, mra, prec, thr, derivative): From 7a8821203be6b5ae94cbc12d1896e566a023588c Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Thu, 25 Jan 2024 10:57:48 +0100 Subject: [PATCH 11/28] Bugfix --- two_electron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/two_electron.py b/two_electron.py index 8b60963..f46f637 100644 --- a/two_electron.py +++ b/two_electron.py @@ -229,7 +229,7 @@ def coulomb_2e_D2_J(spinors, potential, mra, prec, thr, derivative): orbital_energy = c2 * ( -1.0 + np.sqrt(1 + 2 * (classical_energy) / c2)) delta_e = np.abs(orbital_energy - old_energy) Jorb = Jop(spinors[0]) - hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) + hd_psi = orb.apply_dirac_hamiltonian(spinors[0], prec, der = derivative) Vorb = Vop(spinors[0]) Jenergy = (spinors[0].dot(Jorb)).real hdenergy = (spinors[0].dot(hd_psi)).real From abe2b9e7ec3870eecf0c15570f66c5aef1ef0de6 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Thu, 25 Jan 2024 15:18:41 +0100 Subject: [PATCH 12/28] prints energy components in D and D2 --- two_electron.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/two_electron.py b/two_electron.py index f46f637..26b8a59 100644 --- a/two_electron.py +++ b/two_electron.py @@ -234,7 +234,8 @@ def coulomb_2e_D2_J(spinors, potential, mra, prec, thr, derivative): Jenergy = (spinors[0].dot(Jorb)).real hdenergy = (spinors[0].dot(hd_psi)).real Venergy = -(spinors[0].dot(Vorb)).real - orbital_energy_dirac = hdenergy + Jenergy + Venergy + print("hd J V", hdenergy - c2, Jenergy, Venergy) + orbital_energy_dirac = hdenergy + 0.5*Jenergy + Venergy - c2 total_energy = 2.0 * orbital_energy - 0.5 * Jenergy print("Final classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", cke + cpe, "delta_e", delta_e) print("Final orbital energy: ", orbital_energy) @@ -315,6 +316,7 @@ def coulomb_gs_2e(spinorb1, potential, mra, prec, thr, derivative): if(error_norm < thr and delta_e < thr/1000): compute_last_energy = True old_energy = total_e + print("hd J V", hd_11 - light_speed**2, JmK, V1) return spinorb1, spinorb2 #def coulomb_gs(potential, spinors, mra, prec, der = 'ABGV'): From 6b3f66c9051f82dda9f2b834d42e247aa2fcf131 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Tue, 19 Mar 2024 18:50:15 +0100 Subject: [PATCH 13/28] alternative exp val --- two_electron.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/two_electron.py b/two_electron.py index 26b8a59..336c8ed 100644 --- a/two_electron.py +++ b/two_electron.py @@ -203,6 +203,8 @@ def coulomb_2e_D2_J(spinors, potential, mra, prec, thr, derivative): RHS = build_RHS_D2(Jop, Vop, spinors[0], prec, light_speed) cke = spinors[0].classicT() cpe = (spinors[0].dot(RHS)).real + cpe_alt = pot_exp_val_D2(Jop, Vop, spinor, prec, light_speed) + print("CPE and alt", cpe, cpe_alt) classical_energy = cke + cpe orbital_energy = c2 * ( -1.0 + np.sqrt(1 + 2 * (classical_energy) / c2)) print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", classical_energy) @@ -682,3 +684,22 @@ def build_RHS_D2(Jop, Vop, spinor, prec, light_speed): RHS = beta_VT_psi + anticom + VT_VT_psi RHS.cropLargeSmall(prec) return RHS + +def pot_exp_val_D2(Jop, Vop, spinor, prec, light_speed): + c2 = light_speed**2 + Jpsi = Jop(spinor) + Vpsi = Vop(spinor) + VT_psi = 0.5 * Jpsi - Vpsi + + beta_VT_psi = VT_psi.beta2() + beta_VT_psi.cropLargeSmall(prec) + + ap_psi = spinor.alpha_p(prec*light_speed) + ap_psi.cropLargeSmall(prec) + + exp1 = spinor.dot(beta_VT_psi).real + exp2 = ap_psi.dot(VT_psi).real + exp3 = VT_psi.dot(VT_psi).real + + return exp1 + exp2/light_speed + exp3/(2*c2) + From 569d6adc9dc952d30af27dead83f059296934b03 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Mon, 6 May 2024 17:12:25 +0200 Subject: [PATCH 14/28] pythonized input --- orbital4c/nuclear_potential.py | 69 +------- test.py | 303 +++++++++++++-------------------- two_electron.py | 4 +- 3 files changed, 124 insertions(+), 252 deletions(-) diff --git a/orbital4c/nuclear_potential.py b/orbital4c/nuclear_potential.py index 787e51e..cbc0d05 100644 --- a/orbital4c/nuclear_potential.py +++ b/orbital4c/nuclear_potential.py @@ -5,30 +5,13 @@ from vampyr import vampyr3d as vp from vampyr import vampyr1d as vp1 -import argparse import numpy as np -import numpy.linalg as LA -import sys, getopt - -def read_file_with_named_lists(atomlist): - charge_list = {"H" : 1, "He": 2, "Ne": 10, "Ar": 18, "Kr": 36, "Xe": 54, "Rn": 86, "Th": 90, "U":92, "Pu": 94} - atom_list = {} - index = 0 - with open(atomlist, 'r') as file: - for line in file: - terms = line.strip().split() - charge = charge_list[terms[0]] - atom_list[index] = [terms[0], charge, float(terms[1]), float(terms[2]), float(terms[3])] - index += 1 - number = len(atom_list) - return atom_list, number - def calculate_center_of_mass(atoms_list): total_mass = 0.0 center_of_mass = [0.0, 0.0, 0.0] - for atom in atoms_list.values(): + for atom in atoms_list: # Assuming each atom has mass 1.0 (modify if necessary) mass = 1.0 total_mass += mass @@ -43,57 +26,9 @@ def calculate_center_of_mass(atoms_list): return center_of_mass -#def pot(coordinates, typenuc, mra, prec, der): -# atomic_potentials = [] -# V_tree = vp.FunctionTree(mra) -# V_tree.setZero() -# for atom, origin in coordinates.items(): -# atom = get_original_list_name(atom) -# print("Atom:", atom) -# fileObj = open("Z.txt", "r") -# charge = "" -# for line in fileObj: -# if not line.startswith("#"): -# line = line.strip().split() -# if len(line) == 2: -# if line[0] == atom: -# charge = float(line[1]) -# print("Charge:", charge) -# fileObj.close() -# print("Origin:", origin) -# print() # Print an empty line for separation -# -# if typenuc == 'point_charge': -# Peps = vp.ScalingProjector(mra,prec/10) -# f = lambda x: point_charge(x, origin, charge) -# V = Peps(f) -# elif typenuc == 'coulomb_HFYGB': -# Peps = vp.ScalingProjector(mra,prec/10) -# f = lambda x: coulomb_HFYGB(x, origin, charge, prec) -# V = Peps(f) -# elif typenuc == 'homogeneus_charge_sphere': -# Peps = vp.ScalingProjector(mra,prec/10) -# f = lambda x: homogeneus_charge_sphere(x, origin, charge, atom) -# V = Peps(f) -# elif typenuc == 'gaussian': -# Peps = vp.ScalingProjector(mra,prec/10) -# f = lambda x: gaussian(x, origin, charge, atom) -# V = Peps(f) -# print("Potential for atom ", atom) -# print(V) -# atomic_potentials.append(V) -## vp.advanced.add(prec, V_tree, atomic_potentials) -# V_tree = atomic_potentials[0] + atomic_potentials[1] -# print('Define V Potential', typenuc, 'DONE') -# return V_tree -# - - - - def nuclear_potential(position, atoms_list, typenuc, mra, prec, der): potential = 0 - for atom in atoms_list.values(): + for atom in atoms_list: charge = atom[1] atomname = atom[0] atom_coordinates = [atom[2], atom[3], atom[4]] diff --git a/test.py b/test.py index 4a17e0f..1d48076 100644 --- a/test.py +++ b/test.py @@ -21,189 +21,124 @@ import importlib importlib.reload(orb) -def read_step_file(filename): - steps = {} - with open(filename, 'r') as file: - for line in file: - terms = line.strip().split() - option = terms[0] - value = terms[1] - steps[option] = value - return steps - -if __name__ == '__main__': - parser = argparse.ArgumentParser(description='Collecting all data tostart the program.') - parser.add_argument('-d', '--derivative', dest='deriv', type=str, default='ABGV', - help='put the type of derivative') - parser.add_argument('-b', '--box', dest='box', type=int, default=100, - help='put the box dimension') - parser.add_argument('-l', '--light_speed', dest='lux_speed', type=float, default=137.03599913900001, - help='light of speed') - parser.add_argument('-o', '--order', dest='order', type=int, default=10, - help='put the order of Polinomial') - parser.add_argument('-p', '--prec', dest='prec', type=float, default=1e-4, - help='put the precision') - parser.add_argument('-t', '--threshold', dest='thr', type=float, default=1e-4, - help='put the orbital threshold') - parser.add_argument('-e', '--coulgau', dest='coulgau', type=str, default='coulomb', - help='put the coulomb or gaunt') - parser.add_argument('-v', '--potential', dest='potential', type=str, default='coulomb_HFYGB', - help='tell me wich model for V you want to use point_charge, coulomb_HFYGB, homogeneus_charge_sphere, gaussian') - args = parser.parse_args() - - assert args.coulgau in ['coulomb', 'gaunt', 'gaunt-test'], 'Please, specify coulgau in a rigth way: coulomb or gaunt' - - assert args.potential in ['point_charge', 'coulomb_HFYGB', 'homogeneus_charge_sphere', 'gaussian'], 'Please, specify V' - - assert args.deriv in ['PH', 'BS', 'ABGV'], 'Please, specify the type of derivative' - - ################# Define Paramters ########################### - light_speed = args.lux_speed - alpha = 1/light_speed - k = -1 - l = 0 - n = 1 - m = 0.5 - - ################# Call MRA ####################### - mra = vp.MultiResolutionAnalysis(box=[-args.box,args.box], order=args.order, max_depth = 30) - prec = args.prec - thr = args.thr - orb.orbital4c.mra = mra - orb.orbital4c.light_speed = light_speed - cf.complex_fcn.mra = mra - derivative = args.deriv - print('call MRA DONE') - - steplist = 'step_list.txt' - atomlist = 'atom_list.txt' - steps = read_step_file(steplist) - coordinates, number = nucpot.read_file_with_named_lists(atomlist) - - ################## Jobs ########################## - computeNuclearPotential = False - readOrbitals = False - readPotential = False - runD_1e = False - runD2_1e = False - runCoulombGen = False - runCoulomb2e = False - runKutzelnigg = False - runKutzSimple = False - saveOrbitals = False - savePotential = False - runGaunt = False - runGaugeA = False - runGaugeB = False - runGaugeC = False - runGaugeD = False - runGaugeDelta = False - - for key in steps: - flag = (steps[key] == 'True') - locals()[key] = flag - print(key, steps[key], locals()[key]) - - if(computeNuclearPotential): print("hello") - - ################### Reading Atoms ######################### - steplist = 'step_list.txt' - atomlist = 'atom_list.txt' - coordinates, number = nucpot.read_file_with_named_lists(atomlist) - - - ################### Define V potential ###################### - V_tree = vp.FunctionTree(mra) - if(computeNuclearPotential): - Peps = vp.ScalingProjector(mra, prec/10) - typenuc = args.potential - f = lambda x: nucpot.nuclear_potential(x, coordinates, typenuc, mra, prec, derivative) - V_tree = Peps(f) - print("Define V", args.potential, "DONE") - com_coordinates = nucpot.calculate_center_of_mass(coordinates) - if(savePotential): - V_tree.saveTree(f"potential") +import fileinput + +def make_starting_guess(): + gauss_tree_tot = vp.FunctionTree(mra) + gauss_tree_tot.setZero() + a_coeff = 3.0 + b_coeff = np.sqrt(a_coeff/np.pi)**3 + AO_list = [] + for atom in coordinates: + gauss = vp.GaussFunc(b_coeff, a_coeff, [atom[2], atom[3], atom[4]]) + gauss_tree = vp.FunctionTree(mra) + vp.advanced.build_grid(out=gauss_tree, inp=gauss) + vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss) + AO_list.append(gauss_tree) + gauss_tree_sum = vp.sum(AO_list) + + La_comp = cf.complex_fcn() + La_comp.copy_fcns(real = gauss_tree_sum) + spinorb1.copy_components(La = La_comp) + spinorb1.init_small_components(prec/10) + spinorb1.normalize() + spinorb1.cropLargeSmall(prec) + spinorb2 = spinorb1.ktrs(prec) + return spinorb1, spinorb2 + +input_blob = "" +for line in fileinput.input(): + input_blob += line +exec(input_blob) + +################# Call MRA ####################### +mra = vp.MultiResolutionAnalysis(box=[-box, box], order=order, max_depth = 30) +orb.orbital4c.mra = mra +orb.orbital4c.light_speed = light_speed +cf.complex_fcn.mra = mra +coordinates = molecule + +################### Define V potential ###################### +V_tree = vp.FunctionTree(mra) + +if(computePotential): + Peps = vp.ScalingProjector(mra, prec/10) + typenuc = potential + f = lambda x: nucpot.nuclear_potential(x, coordinates, typenuc, mra, prec, derivative) + V_tree = Peps(f) + print("Define V", potential, "DONE") + com_coordinates = nucpot.calculate_center_of_mass(coordinates) + +if(readPotential): + V_tree.loadTree(f"potential") + +if(savePotential): + V_tree.saveTree(f"potential") - if(readPotential): - V_tree.loadTree(f"potential") - - print("Number of Atoms = ", number) - print(coordinates) - - #############################START WITH CALCULATION################################### - spinorb1 = orb.orbital4c() - spinorb2 = orb.orbital4c() - if readOrbitals: - spinorb1.read("spinorb1") - spinorb2 = spinorb1.ktrs(prec) - else: - gauss_tree_tot = vp.FunctionTree(mra) - gauss_tree_tot.setZero() - a_coeff = 3.0 - b_coeff = np.sqrt(a_coeff/np.pi)**3 - AO_list = [] - for atom in coordinates.values(): - gauss = vp.GaussFunc(b_coeff, a_coeff, [atom[2], atom[3], atom[4]]) - gauss_tree = vp.FunctionTree(mra) - vp.advanced.build_grid(out=gauss_tree, inp=gauss) - vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss) - AO_list.append(gauss_tree) - if number == 1: - gauss_tree_tot = AO_list[0] - elif number == 2: - gauss_tree_tot = AO_list[0] + AO_list[1] - gauss_tree_tot.normalize() - - La_comp = cf.complex_fcn() - La_comp.copy_fcns(real = gauss_tree_tot) - spinorb1.copy_components(La = La_comp) - spinorb1.init_small_components(prec/10) - spinorb1.normalize() - spinorb1.cropLargeSmall(prec) - spinorb2 = spinorb1.ktrs(prec) - - length = 2 * args.box - print("Using derivative ", derivative) - - if runD_1e: - spinorb1 = one_electron.gs_D_1e(spinorb1, V_tree, mra, prec, thr, derivative) - - if runD2_1e: - spinorb1 = one_electron.gs_D2_1e(spinorb1, V_tree, mra, prec, thr, derivative) - - if runCoulombGen: - print("NOT ACTIVE") -# spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec, derivative) - - if runCoulomb2e: - spinorb1, spinorb2 = two_electron.coulomb_gs_2e(spinorb1, V_tree, mra, prec, thr, derivative) - - if runKutzelnigg: - print("NOT ACTIVE") -# spinorb1, spinorb2 = two_electron.coulomb_2e_D2([spinorb1, spinorb2], V_tree, mra, prec, derivative) - - if runKutzSimple: - spinorb1, spinorb2 = two_electron.coulomb_2e_D2_J([spinorb1, spinorb2], V_tree, mra, prec, thr, derivative) - - if runGaunt: - two_electron.calcGauntPert(spinorb1, spinorb2, mra, prec) - - if runGaugeA: - two_electron.calcGaugePertA(spinorb1, spinorb2, mra, prec) - - if runGaugeB: - two_electron.calcGaugePertB(spinorb1, spinorb2, mra, prec) - - if runGaugeC: - two_electron.calcGaugePertC(spinorb1, spinorb2, mra, prec) - - if runGaugeD: - two_electron.calcGaugePertD(spinorb1, spinorb2, mra, prec) - - if runGaugeDelta: - two_electron.calcGaugeDelta(spinorb1, spinorb2, mra, prec) - - if saveOrbitals: - spinorb1.save("spinorb1") - #spinorb2.save("spinorb2") +print("Number of Atoms = ", len(molecule)) +print(coordinates) + +#############################START WITH CALCULATION################################### +spinorb1 = orb.orbital4c() +spinorb2 = orb.orbital4c() + +if readOrbitals: + spinorb1.read("spinorb1") + spinorb2 = spinorb1.ktrs(prec) +else: + spinorb1, spinorb2 = make_starting_guess() + +run_D_1e = scf and not D2 and not two_electrons +run_D2_1e = scf and D2 and not two_electrons +run_D_2e = scf and not D2 and two_electrons and not ktrs +run_D2_2e = scf and D2 and two_electrons and not ktrs +run_D_2e_ktrs = scf and not D2 and two_electrons and ktrs +run_D2_2e_ktrs = scf and D2 and two_electrons and ktrs + + +length = 2 * box +print("Using derivative ", derivative) + +if run_D_1e: + spinorb1 = one_electron.gs_D_1e(spinorb1, V_tree, mra, prec, thr, derivative) + +if run_D2_1e: + spinorb1 = one_electron.gs_D2_1e(spinorb1, V_tree, mra, prec, thr, derivative) + +if run_D_2e: + print("NOT PROPERLY TESTED") + exit(-1) + spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec, derivative) + +if run_D2_2e: + print("NOT PROPERLY TESTED") + exit(-1) + spinorb1, spinorb2 = two_electron.coulomb_2e_D2([spinorb1, spinorb2], V_tree, mra, prec, derivative) + +if run_D_2e_ktrs: + spinorb1, spinorb2 = two_electron.coulomb_gs_2e(spinorb1, V_tree, mra, prec, thr, derivative) + +if run_D2_2e_ktrs: + spinorb1, spinorb2 = two_electron.coulomb_2e_D2_J([spinorb1, spinorb2], V_tree, mra, prec, thr, derivative) + +if runGaunt: + two_electron.calcGauntPert(spinorb1, spinorb2, mra, prec) + +if runGaugeA: + two_electron.calcGaugePertA(spinorb1, spinorb2, mra, prec) + +if runGaugeB: + two_electron.calcGaugePertB(spinorb1, spinorb2, mra, prec) + +if runGaugeC: + two_electron.calcGaugePertC(spinorb1, spinorb2, mra, prec) + +if runGaugeD: + two_electron.calcGaugePertD(spinorb1, spinorb2, mra, prec) + +if runGaugeDelta: + two_electron.calcGaugeDelta(spinorb1, spinorb2, mra, prec) + +if saveOrbitals: + spinorb1.save("spinorb1") diff --git a/two_electron.py b/two_electron.py index 336c8ed..b76e78f 100644 --- a/two_electron.py +++ b/two_electron.py @@ -203,7 +203,7 @@ def coulomb_2e_D2_J(spinors, potential, mra, prec, thr, derivative): RHS = build_RHS_D2(Jop, Vop, spinors[0], prec, light_speed) cke = spinors[0].classicT() cpe = (spinors[0].dot(RHS)).real - cpe_alt = pot_exp_val_D2(Jop, Vop, spinor, prec, light_speed) + cpe_alt = pot_exp_val_D2(Jop, Vop, spinors[0], prec, light_speed) print("CPE and alt", cpe, cpe_alt) classical_energy = cke + cpe orbital_energy = c2 * ( -1.0 + np.sqrt(1 + 2 * (classical_energy) / c2)) @@ -295,7 +295,9 @@ def coulomb_gs_2e(spinorb1, potential, mra, prec, thr, derivative): V_J_K_spinorb1 = v_psi_1 + JmK_phi1 + print(eps,light_speed) mu = orb.calc_dirac_mu(eps, light_speed) + print(mu) tmp = orb.apply_helmholtz(V_J_K_spinorb1, mu, prec) new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, eps, der = derivative) new_orbital *= 0.5/light_speed**2 From 3cd753eda24648d9823118f6017942b8d1500acf Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Wed, 9 Oct 2024 15:50:31 +0200 Subject: [PATCH 15/28] Adds NR starting guess --- one_electron.py | 8 +-- orbital4c/nuclear_potential.py | 13 +---- orbital4c/orbital.py | 12 +++-- starting_guess.py | 71 +++++++++++++++++++++++++++ test.py | 90 ++++++++++++++++------------------ 5 files changed, 125 insertions(+), 69 deletions(-) create mode 100644 starting_guess.py diff --git a/one_electron.py b/one_electron.py index 27c4d0d..9eb3418 100644 --- a/one_electron.py +++ b/one_electron.py @@ -17,7 +17,7 @@ def analytic_1s(light_speed, n, k, Z): tmp3 = 1 + tmp2**2 return light_speed**2 / np.sqrt(tmp3) -def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative): +def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('One-electron calculations') error_norm = 1 @@ -55,7 +55,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi energy = spinorb1.dot(add_psi).real - energy_1s = analytic_1s(light_speed, 1, -1, 1) + energy_1s = analytic_1s(light_speed, 1, -1, charge) print("Exact energy: ", energy_1s - light_speed**2) print('Final Energy:',energy - light_speed**2) print('Delta Energy:',energy - old_energy) @@ -63,7 +63,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative): return spinorb1 -def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative): +def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('Hartree-Fock 1e D2') error_norm = 1 delta_e = 1 @@ -124,7 +124,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative): cpe = (spinorb1.dot(RHS)).real classic_energy = cke + cpe energy = c2*(np.sqrt(1+2*classic_energy/c2)-1) - energy_1s = analytic_1s(light_speed, 1, -1, 1) + energy_1s = analytic_1s(light_speed, 1, -1, charge) # print('Kutzelnigg =',cke, cpe, energy_kutzelnigg) # print('Quadratic approx =',energy_kutzelnigg - energy_kutzelnigg**2/(2*c2)) diff --git a/orbital4c/nuclear_potential.py b/orbital4c/nuclear_potential.py index cbc0d05..080d9b3 100644 --- a/orbital4c/nuclear_potential.py +++ b/orbital4c/nuclear_potential.py @@ -98,18 +98,7 @@ def homogeneus_charge_sphere(position, center, charge, RMS): return prec * factor -def gaussian(position, center, charge, atomname): - fileObj = open("./orbital4c/param_V.txt", "r") - for line in fileObj: - if not line.startswith("#"): - line = line.strip().split() - if len(line) == 3: - if line[0] == atomname: - epsilon = line[2] - else: - print("Data file not correclty formatted! Please check it!") - fileObj.close() - epsilon = float(epsilon) +def gaussian_potential(position, center, charge, epsilon): d2 = ((position[0] - center[0]) ** 2 + (position[1] - center[1]) ** 2 + (position[2] - center[2]) ** 2) diff --git a/orbital4c/orbital.py b/orbital4c/orbital.py index 0cfdbe5..7e60721 100644 --- a/orbital4c/orbital.py +++ b/orbital4c/orbital.py @@ -440,7 +440,9 @@ def alpha_gradient(orbital, prec): return out def calc_dirac_mu(energy, light_speed): - return np.sqrt((light_speed**4-energy**2)/light_speed**2) + val = (light_speed**4-energy**2)/light_speed**2 + print("calc_dirac_mu", val, energy, light_speed) + return np.sqrt(val) def calc_kutzelnigg_mu(energy_sq, light_speed): c2 = light_speed**2 @@ -448,9 +450,11 @@ def calc_kutzelnigg_mu(energy_sq, light_speed): return np.sqrt(-val) def calc_non_rel_mu(energy): - if energy > 0: - return np.sqrt(2.0 * energy) - elif energy < 0: + if energy < 0: return np.sqrt(-2.0 * energy) + else: + print("Positive energy") + exit(-1) + diff --git a/starting_guess.py b/starting_guess.py new file mode 100644 index 0000000..41a76af --- /dev/null +++ b/starting_guess.py @@ -0,0 +1,71 @@ +from vampyr import vampyr3d as vp +import numpy as np +from scipy.special import eval_genlaguerre +from orbital4c import complex_fcn as cf +from orbital4c import orbital as orb + +def make_starting_guess(mra, prec): + gauss_tree_tot = vp.FunctionTree(mra) + gauss_tree_tot.setZero() + a_coeff = 3.0 + b_coeff = np.sqrt(a_coeff/np.pi)**3 + AO_list = [] + for atom in coordinates: + gauss = vp.GaussFunc(b_coeff, a_coeff, [atom[2], atom[3], atom[4]]) + gauss_tree = vp.FunctionTree(mra) + vp.advanced.build_grid(out=gauss_tree, inp=gauss) + vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss) + AO_list.append(gauss_tree) + gauss_tree_sum = vp.sum(AO_list) + + La_comp = cf.complex_fcn() + La_comp.copy_fcns(real = gauss_tree_sum) + spinorb1.copy_components(La = La_comp) + spinorb1.init_small_components(prec/10) + spinorb1.normalize() + spinorb1.cropLargeSmall(prec) + return spinorb1 + +def make_NR_starting_guess(position, charge, mra, prec): + nr_wf_tree = vp.FunctionTree(mra) + nr_wf_tree.setZero() + n = 1 + l = 0 + Peps = vp.ScalingProjector(mra, prec) + guess = lambda x : wf_hydrogenionic_atom(n,l,[x[0]-position[0], x[1]-position[1], x[2]-position[2]],charge) + nr_wf_tree = Peps(guess) + + La_comp = cf.complex_fcn() + La_comp.copy_fcns(real = nr_wf_tree) + + spinorb1 = orb.orbital4c() + spinorb2 = orb.orbital4c() + spinorb1.copy_components(La = La_comp) + spinorb1.init_small_components(prec/10) + spinorb1.normalize() + spinorb1.cropLargeSmall(prec) + return spinorb1 + +#returns the value of the radial WF in the point r +# 1. the nucleus is assumed infintely heavy (mass of electron and Bohr radius used) +# 2. the nucleus is placed in the origin +# 3. atomic units are assumed a0 = 1 hbar = 1 me = 1 4pie0 = 1 +def radial_wf_hydrogenionic_atom(n,l,r,Z): + rho = 2 * Z * r + slater = np.exp(-rho/2) + polynomial = eval_genlaguerre(n-l-1, 2*l+1, rho) + f1 = np.math.factorial(n-l-1) + f2 = np.math.factorial(n+l) + norm = np.sqrt((2*Z/n)**3 * f1 / (2 * n * f2)) + value = norm * rho**l * polynomial * slater + return value + +def wf_hydrogenionic_atom(n,l,position,Z): + if(l != 0): + print("only s orbitals for now") + exit(-1) + distance = np.sqrt(position[0]**2 + position[1]**2 + position[2]**2) + value = radial_wf_hydrogenionic_atom(n, l, distance, Z) + return value + + diff --git a/test.py b/test.py index 1d48076..b73ea8b 100644 --- a/test.py +++ b/test.py @@ -15,78 +15,70 @@ import numpy.linalg as LA import sys, getopt -import one_electron -import two_electron +import one_electron as oneel +import two_electron as twoel +import starting_guess as sg import importlib importlib.reload(orb) import fileinput -def make_starting_guess(): - gauss_tree_tot = vp.FunctionTree(mra) - gauss_tree_tot.setZero() - a_coeff = 3.0 - b_coeff = np.sqrt(a_coeff/np.pi)**3 - AO_list = [] - for atom in coordinates: - gauss = vp.GaussFunc(b_coeff, a_coeff, [atom[2], atom[3], atom[4]]) - gauss_tree = vp.FunctionTree(mra) - vp.advanced.build_grid(out=gauss_tree, inp=gauss) - vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss) - AO_list.append(gauss_tree) - gauss_tree_sum = vp.sum(AO_list) - - La_comp = cf.complex_fcn() - La_comp.copy_fcns(real = gauss_tree_sum) - spinorb1.copy_components(La = La_comp) - spinorb1.init_small_components(prec/10) - spinorb1.normalize() - spinorb1.cropLargeSmall(prec) - spinorb2 = spinorb1.ktrs(prec) - return spinorb1, spinorb2 - input_blob = "" for line in fileinput.input(): input_blob += line + +print(input_blob) exec(input_blob) +# +# 1. This code works now only for atoms and up to two electrons with KTRS +# 2. Orbital guess obtained by using NR hydrogenionic 1s orbital +# 3. The input file is a Python code mostly containing variable allocation +# 4. Input parsing is executing that python input after reading it +# 5. Nuclear potential selected manually (Gaussian now to reproduce Harrison's results) +# + ################# Call MRA ####################### mra = vp.MultiResolutionAnalysis(box=[-box, box], order=order, max_depth = 30) orb.orbital4c.mra = mra orb.orbital4c.light_speed = light_speed cf.complex_fcn.mra = mra -coordinates = molecule +charge = molecule[0][1] +position = [molecule[0][2],molecule[0][3],molecule[0][4]] +radius = molecule[0][5] +epsilon = molecule[0][6] + +print(charge,position) ################### Define V potential ###################### +Peps = vp.ScalingProjector(mra, prec/10) V_tree = vp.FunctionTree(mra) if(computePotential): - Peps = vp.ScalingProjector(mra, prec/10) typenuc = potential - f = lambda x: nucpot.nuclear_potential(x, coordinates, typenuc, mra, prec, derivative) + f = lambda x: nucpot.gaussian_potential(x, position, charge, epsilon) V_tree = Peps(f) - print("Define V", potential, "DONE") - com_coordinates = nucpot.calculate_center_of_mass(coordinates) - -if(readPotential): +elif(readPotential): V_tree.loadTree(f"potential") if(savePotential): V_tree.saveTree(f"potential") print("Number of Atoms = ", len(molecule)) -print(coordinates) +print(molecule) #############################START WITH CALCULATION################################### spinorb1 = orb.orbital4c() spinorb2 = orb.orbital4c() - if readOrbitals: - spinorb1.read("spinorb1") - spinorb2 = spinorb1.ktrs(prec) + spinorb1.read(orbitalName) else: - spinorb1, spinorb2 = make_starting_guess() + spinorb1 = sg.make_NR_starting_guess(position, charge, mra, prec) +spinorb2 = spinorb1.ktrs(prec) + +if saveGuess: + spinorb1.save("guess1") run_D_1e = scf and not D2 and not two_electrons run_D2_1e = scf and D2 and not two_electrons @@ -100,44 +92,44 @@ def make_starting_guess(): print("Using derivative ", derivative) if run_D_1e: - spinorb1 = one_electron.gs_D_1e(spinorb1, V_tree, mra, prec, thr, derivative) + spinorb1 = oneel.gs_D_1e(spinorb1, V_tree, mra, prec, thr, derivative, charge) if run_D2_1e: - spinorb1 = one_electron.gs_D2_1e(spinorb1, V_tree, mra, prec, thr, derivative) + spinorb1 = oneel.gs_D2_1e(spinorb1, V_tree, mra, prec, thr, derivative, charge) if run_D_2e: print("NOT PROPERLY TESTED") exit(-1) - spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec, derivative) + spinorb1, spinorb2 = twoel.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec, derivative) if run_D2_2e: print("NOT PROPERLY TESTED") exit(-1) - spinorb1, spinorb2 = two_electron.coulomb_2e_D2([spinorb1, spinorb2], V_tree, mra, prec, derivative) + spinorb1, spinorb2 = twoel.coulomb_2e_D2([spinorb1, spinorb2], V_tree, mra, prec, derivative) if run_D_2e_ktrs: - spinorb1, spinorb2 = two_electron.coulomb_gs_2e(spinorb1, V_tree, mra, prec, thr, derivative) + spinorb1, spinorb2 = twoel.coulomb_gs_2e(spinorb1, V_tree, mra, prec, thr, derivative) if run_D2_2e_ktrs: - spinorb1, spinorb2 = two_electron.coulomb_2e_D2_J([spinorb1, spinorb2], V_tree, mra, prec, thr, derivative) + spinorb1, spinorb2 = twoel.coulomb_2e_D2_J([spinorb1, spinorb2], V_tree, mra, prec, thr, derivative) if runGaunt: - two_electron.calcGauntPert(spinorb1, spinorb2, mra, prec) + twoel.calcGauntPert(spinorb1, spinorb2, mra, prec) if runGaugeA: - two_electron.calcGaugePertA(spinorb1, spinorb2, mra, prec) + twoel.calcGaugePertA(spinorb1, spinorb2, mra, prec) if runGaugeB: - two_electron.calcGaugePertB(spinorb1, spinorb2, mra, prec) + twoel.calcGaugePertB(spinorb1, spinorb2, mra, prec) if runGaugeC: - two_electron.calcGaugePertC(spinorb1, spinorb2, mra, prec) + twoel.calcGaugePertC(spinorb1, spinorb2, mra, prec) if runGaugeD: - two_electron.calcGaugePertD(spinorb1, spinorb2, mra, prec) + twoel.calcGaugePertD(spinorb1, spinorb2, mra, prec) if runGaugeDelta: - two_electron.calcGaugeDelta(spinorb1, spinorb2, mra, prec) + twoel.calcGaugeDelta(spinorb1, spinorb2, mra, prec) if saveOrbitals: spinorb1.save("spinorb1") From b090df671482b31ccf95d04aacd590aed8b58a14 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Tue, 19 Nov 2024 11:18:31 +0100 Subject: [PATCH 16/28] Tighter convergence check for 1e case --- one_electron.py | 11 ++++++++--- orbital4c/nuclear_potential.py | 6 +++--- test.py | 10 +++++++++- 3 files changed, 20 insertions(+), 7 deletions(-) diff --git a/one_electron.py b/one_electron.py index 9eb3418..9a1cddd 100644 --- a/one_electron.py +++ b/one_electron.py @@ -18,7 +18,7 @@ def analytic_1s(light_speed, n, k, Z): return light_speed**2 / np.sqrt(tmp3) def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): - print('One-electron calculations') + print('One-electron calculations', prec) error_norm = 1 #compute_last_energy = False @@ -26,7 +26,8 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): light_speed = spinorb1.light_speed old_energy = 0 delta_e = 1 - while (error_norm > thr and delta_e > thr/1000): + idx = 0 + while ((error_norm > thr or delta_e > thr/1000) and idx < 100): hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi @@ -50,6 +51,8 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('Energy',energy - light_speed**2) old_energy = energy spinorb1 = new_orbital + print('Converged? ', error_norm, ' > ', thr, ' ---- ', delta_e, ' > ',thr/1000, ' ---- iteration', idx) + idx += 1 hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) @@ -70,7 +73,8 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): light_speed = spinorb1.light_speed c2 = light_speed * light_speed old_energy = 0 - while (error_norm > thr and delta_e > thr/1000): + idx = 0 + while ((error_norm > thr or delta_e > thr/1000) and idx < 100): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() @@ -102,6 +106,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('Energy',energy, old_energy) old_energy = energy spinorb1 = new_orbital + idx += 1 hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) diff --git a/orbital4c/nuclear_potential.py b/orbital4c/nuclear_potential.py index 080d9b3..69f6ff3 100644 --- a/orbital4c/nuclear_potential.py +++ b/orbital4c/nuclear_potential.py @@ -103,6 +103,6 @@ def gaussian_potential(position, center, charge, epsilon): (position[1] - center[1]) ** 2 + (position[2] - center[2]) ** 2) distance = np.sqrt(d2) - prec = charge / distance - u = erf(np.sqrt(epsilon) * distance) - return prec * u + point_charge_potential = charge / distance + gaussian_screening = erf(np.sqrt(epsilon) * distance) + return point_charge_potential * gaussian_screening diff --git a/test.py b/test.py index b73ea8b..97d1070 100644 --- a/test.py +++ b/test.py @@ -57,7 +57,15 @@ if(computePotential): typenuc = potential - f = lambda x: nucpot.gaussian_potential(x, position, charge, epsilon) + f = 0 + if(potential == "gaussian"): + f = lambda x: nucpot.gaussian_potential(x, position, charge, epsilon) + elif(potential == "coulomb_HFYGB"): + f = lambda x: nucpot.coulomb_HFYGB(x, position, charge, prec/10) + elif(potential == "point_charge"): + f = lambda x: nucpot.point_charge(x, position, charge) + else: + exit(-1) V_tree = Peps(f) elif(readPotential): V_tree.loadTree(f"potential") From 4e5a714646b93cb346d1c2bb79e8106f46266d49 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Tue, 19 Mar 2024 18:52:06 +0100 Subject: [PATCH 17/28] Poor mans KAIN --- two_electron.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/two_electron.py b/two_electron.py index b76e78f..2b6fa2b 100644 --- a/two_electron.py +++ b/two_electron.py @@ -220,8 +220,9 @@ def coulomb_2e_D2_J(spinors, potential, mra, prec, thr, derivative): error_norm = np.sqrt(deltasq) print('Orbital_Error norm', error_norm) print('delta_e', delta_e) - spinors[0] = new_spinor - spinors[1] = new_spinor.ktrs(prec) + spinors[0] = new_spinor + spinors[0] + spinors[0].normalize() + spinors[1] = spinors[0].ktrs(prec) old_energy = orbital_energy Jop = oper.CoulombDirectOperator(mra, prec, spinors) RHS = build_RHS_D2(Jop, Vop, spinors[0], prec, light_speed) From 8dc0a72d1a0975b7b168533c30416d73acc94737 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Tue, 19 Nov 2024 11:20:52 +0100 Subject: [PATCH 18/28] Some changes --- one_electron.py | 1 + 1 file changed, 1 insertion(+) diff --git a/one_electron.py b/one_electron.py index 9a1cddd..9f4f124 100644 --- a/one_electron.py +++ b/one_electron.py @@ -21,6 +21,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('One-electron calculations', prec) error_norm = 1 + Z = coordinates[0][1] #compute_last_energy = False light_speed = spinorb1.light_speed From 9801c9b4a72631a951bf242ecf7bf0019d069ac8 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Tue, 19 Nov 2024 11:57:19 +0100 Subject: [PATCH 19/28] Bug fixes --- one_electron.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/one_electron.py b/one_electron.py index 9f4f124..822263f 100644 --- a/one_electron.py +++ b/one_electron.py @@ -21,7 +21,6 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('One-electron calculations', prec) error_norm = 1 - Z = coordinates[0][1] #compute_last_energy = False light_speed = spinorb1.light_speed @@ -97,6 +96,8 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print("mu =", mu) new_orbital = orb.apply_helmholtz(RHS, mu, prec) new_orbital.cropLargeSmall(prec) + if(idx > 10): + new_orbital = new_orbital + spinorb1 new_orbital.normalize() delta_psi = new_orbital - spinorb1 deltasq = delta_psi.squaredNorm() @@ -106,7 +107,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('Delta E', delta_e) print('Energy',energy, old_energy) old_energy = energy - spinorb1 = new_orbital + spinorb1 = new_orbital idx += 1 hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) From cb63df0a287652c725dac2631b49233c2b41ac8c Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 13 Dec 2024 10:00:46 +0100 Subject: [PATCH 20/28] Uniform printing --- one_electron.py | 60 ++++++++++++++++++++++++++------------------ orbital4c/orbital.py | 3 ++- test.py | 8 +++++- 3 files changed, 45 insertions(+), 26 deletions(-) diff --git a/one_electron.py b/one_electron.py index 822263f..a21019c 100644 --- a/one_electron.py +++ b/one_electron.py @@ -21,13 +21,13 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('One-electron calculations', prec) error_norm = 1 - #compute_last_energy = False light_speed = spinorb1.light_speed + c2 = light_speed**2 old_energy = 0 delta_e = 1 idx = 0 - while ((error_norm > thr or delta_e > thr/1000) and idx < 100): + while ((error_norm > thr or delta_e > prec/10) and idx < 100): hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi @@ -38,7 +38,8 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): tmp.cropLargeSmall(prec) new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = derivative) # new_orbital = orb.apply_helmholtz(tmp, mu, prec) -# new_orbital.cropLargeSmall(prec) + if(idx > 10): + new_orbital = new_orbital + spinorb1 new_orbital.cropLargeSmall(prec) new_orbital.normalize() delta_psi = new_orbital - spinorb1 @@ -51,18 +52,34 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('Energy',energy - light_speed**2) old_energy = energy spinorb1 = new_orbital - print('Converged? ', error_norm, ' > ', thr, ' ---- ', delta_e, ' > ',thr/1000, ' ---- iteration', idx) + print('Converged? ', error_norm, ' > ', thr, ' ---- ', delta_e, ' > ',prec/10, ' ---- iteration', idx) idx += 1 + print(new_orbital) hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi energy = spinorb1.dot(add_psi).real energy_1s = analytic_1s(light_speed, 1, -1, charge) - print("Exact energy: ", energy_1s - light_speed**2) - print('Final Energy:',energy - light_speed**2) - print('Delta Energy:',energy - old_energy) - print('Error Energy:',energy - energy_1s) + + beta_v_psi = v_psi.beta2() + ap_psi = spinorb1.alpha_p(prec, derivative) + cke = spinorb1.classicT() + psi_beta_v_vpsi = spinorb1.dot(beta_v_psi).real + psi_ap_V_psi = ap_psi.dot(v_psi).real + psi_V2_psi = v_psi.dot(v_psi).real + cpe = psi_beta_v_vpsi + psi_ap_V_psi/light_speed + psi_V2_psi/c2 + classic_energy = cke + cpe + print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", classic_energy) + energy_kutzelnigg = c2*(np.sqrt(1+2*classic_energy/c2)-1) + + print('Exact Energy = ', energy_1s - c2) + print('Dirac Energy = ', energy - c2) + print('Kutze Energy = ', energy_kutzelnigg) + print('Error Kutze = ', energy_kutzelnigg - energy_1s + light_speed**2) + print('Error Dirac = ', energy - energy_1s) + print('Delta Energy = ', energy - old_energy) + print('Dirac - Kutzelnigg = ', energy - energy_kutzelnigg - light_speed**2) return spinorb1 @@ -74,7 +91,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): c2 = light_speed * light_speed old_energy = 0 idx = 0 - while ((error_norm > thr or delta_e > thr/1000) and idx < 100): + while ((error_norm > thr or delta_e > prec/10) and idx < 100): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() @@ -95,9 +112,9 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): mu = orb.calc_non_rel_mu(cke+cpe) print("mu =", mu) new_orbital = orb.apply_helmholtz(RHS, mu, prec) - new_orbital.cropLargeSmall(prec) if(idx > 10): new_orbital = new_orbital + spinorb1 + new_orbital.cropLargeSmall(prec) new_orbital.normalize() delta_psi = new_orbital - spinorb1 deltasq = delta_psi.squaredNorm() @@ -108,33 +125,28 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): print('Energy',energy, old_energy) old_energy = energy spinorb1 = new_orbital + print('Converged? ', error_norm, ' > ', thr, ' ---- ', delta_e, ' > ',prec/10, ' ---- iteration', idx) idx += 1 + print(new_orbital) hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi energy_dirac = spinorb1.dot(add_psi).real -# v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) - vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() - apV_psi = v_psi.alpha_p(prec*light_speed, derivative) ap_psi = spinorb1.alpha_p(prec*light_speed, derivative) - Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec*light_speed) - anticom = apV_psi + Vap_psi -# anticom.cropLargeSmall(prec) -# beta_v_psi.cropLargeSmall(prec) -# vv_psi.cropLargeSmall(prec) - RHS = beta_v_psi + vv_psi + anticom * (0.5/light_speed) - RHS.cropLargeSmall(prec) + cke = spinorb1.classicT() - cpe = (spinorb1.dot(RHS)).real + psi_beta_v_vpsi = spinorb1.dot(beta_v_psi).real + psi_ap_V_psi = ap_psi.dot(v_psi).real + psi_V2_psi = v_psi.dot(v_psi).real + cpe = psi_beta_v_vpsi + psi_ap_V_psi/light_speed + psi_V2_psi/c2 classic_energy = cke + cpe energy = c2*(np.sqrt(1+2*classic_energy/c2)-1) energy_1s = analytic_1s(light_speed, 1, -1, charge) - -# print('Kutzelnigg =',cke, cpe, energy_kutzelnigg) -# print('Quadratic approx =',energy_kutzelnigg - energy_kutzelnigg**2/(2*c2)) + + print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", classic_energy) print('Exact Energy = ', energy_1s - light_speed**2) print('Dirac Energy = ', energy_dirac - light_speed**2) print('Kutze Energy = ', energy) diff --git a/orbital4c/orbital.py b/orbital4c/orbital.py index 7e60721..ed9179d 100644 --- a/orbital4c/orbital.py +++ b/orbital4c/orbital.py @@ -442,7 +442,8 @@ def alpha_gradient(orbital, prec): def calc_dirac_mu(energy, light_speed): val = (light_speed**4-energy**2)/light_speed**2 print("calc_dirac_mu", val, energy, light_speed) - return np.sqrt(val) + mu = np.sqrt(val) + return mu def calc_kutzelnigg_mu(energy_sq, light_speed): c2 = light_speed**2 diff --git a/test.py b/test.py index 97d1070..12095ed 100644 --- a/test.py +++ b/test.py @@ -40,7 +40,7 @@ # ################# Call MRA ####################### -mra = vp.MultiResolutionAnalysis(box=[-box, box], order=order, max_depth = 30) +mra = vp.MultiResolutionAnalysis(box=[-box, box], order=order, max_depth = 20) orb.orbital4c.mra = mra orb.orbital4c.light_speed = light_speed cf.complex_fcn.mra = mra @@ -59,10 +59,13 @@ typenuc = potential f = 0 if(potential == "gaussian"): + print("Gaussian potential") f = lambda x: nucpot.gaussian_potential(x, position, charge, epsilon) elif(potential == "coulomb_HFYGB"): + print("Harrison potential") f = lambda x: nucpot.coulomb_HFYGB(x, position, charge, prec/10) elif(potential == "point_charge"): + print("point charge potential") f = lambda x: nucpot.point_charge(x, position, charge) else: exit(-1) @@ -72,6 +75,9 @@ if(savePotential): V_tree.saveTree(f"potential") + +print("Potential") +print(V_tree) print("Number of Atoms = ", len(molecule)) print(molecule) From f24496f89996d1776fc7be4f474bf8e7e42ea832 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 28 Mar 2025 16:01:08 +0100 Subject: [PATCH 21/28] Fix nr of iterations --- one_electron.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/one_electron.py b/one_electron.py index a21019c..9b6d036 100644 --- a/one_electron.py +++ b/one_electron.py @@ -27,7 +27,8 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): old_energy = 0 delta_e = 1 idx = 0 - while ((error_norm > thr or delta_e > prec/10) and idx < 100): +# while ((error_norm > thr or delta_e > prec/10) and idx < 100): + while (idx < 30): hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi @@ -91,7 +92,8 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): c2 = light_speed * light_speed old_energy = 0 idx = 0 - while ((error_norm > thr or delta_e > prec/10) and idx < 100): +# while ((error_norm > thr or delta_e > prec/10) and idx < 100): + while(idx < 30) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() From c0552ab691a6ebd2928c825551b574754dc062f8 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 28 Mar 2025 16:18:34 +0100 Subject: [PATCH 22/28] Bugfix --- one_electron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/one_electron.py b/one_electron.py index 9b6d036..da1179f 100644 --- a/one_electron.py +++ b/one_electron.py @@ -93,7 +93,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): old_energy = 0 idx = 0 # while ((error_norm > thr or delta_e > prec/10) and idx < 100): - while(idx < 30) + while(idx < 30): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() From 831dbca89e9fea2b3313ae3cabd451b45f689db6 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Fri, 25 Apr 2025 13:28:22 +0200 Subject: [PATCH 23/28] Added max_iter --- one_electron.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/one_electron.py b/one_electron.py index da1179f..8b69902 100644 --- a/one_electron.py +++ b/one_electron.py @@ -17,7 +17,7 @@ def analytic_1s(light_speed, n, k, Z): tmp3 = 1 + tmp2**2 return light_speed**2 / np.sqrt(tmp3) -def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): +def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): print('One-electron calculations', prec) error_norm = 1 @@ -28,7 +28,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): delta_e = 1 idx = 0 # while ((error_norm > thr or delta_e > prec/10) and idx < 100): - while (idx < 30): + while (idx < niter): hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi @@ -84,7 +84,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge): return spinorb1 -def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): +def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): print('Hartree-Fock 1e D2') error_norm = 1 delta_e = 1 @@ -93,7 +93,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge): old_energy = 0 idx = 0 # while ((error_norm > thr or delta_e > prec/10) and idx < 100): - while(idx < 30): + while(idx < niter): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) beta_v_psi = v_psi.beta2() From 2ee0b040904e7e793ef8d9e5a2f25df529463534 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Tue, 6 May 2025 10:38:33 +0200 Subject: [PATCH 24/28] missing 0.5 added --- one_electron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/one_electron.py b/one_electron.py index 8b69902..47223f3 100644 --- a/one_electron.py +++ b/one_electron.py @@ -143,7 +143,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): psi_beta_v_vpsi = spinorb1.dot(beta_v_psi).real psi_ap_V_psi = ap_psi.dot(v_psi).real psi_V2_psi = v_psi.dot(v_psi).real - cpe = psi_beta_v_vpsi + psi_ap_V_psi/light_speed + psi_V2_psi/c2 + cpe = psi_beta_v_vpsi + psi_ap_V_psi/light_speed + 0.5 * psi_V2_psi / c2 classic_energy = cke + cpe energy = c2*(np.sqrt(1+2*classic_energy/c2)-1) energy_1s = analytic_1s(light_speed, 1, -1, charge) From 19966a1d3feb154d67631cd2fd2b4e94011ec842 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Tue, 6 May 2025 10:52:53 +0200 Subject: [PATCH 25/28] missing factor 2 --- one_electron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/one_electron.py b/one_electron.py index 47223f3..c219bbb 100644 --- a/one_electron.py +++ b/one_electron.py @@ -69,7 +69,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): psi_beta_v_vpsi = spinorb1.dot(beta_v_psi).real psi_ap_V_psi = ap_psi.dot(v_psi).real psi_V2_psi = v_psi.dot(v_psi).real - cpe = psi_beta_v_vpsi + psi_ap_V_psi/light_speed + psi_V2_psi/c2 + cpe = psi_beta_v_vpsi + psi_ap_V_psi/light_speed + 0.5 * psi_V2_psi / c2 classic_energy = cke + cpe print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", classic_energy) energy_kutzelnigg = c2*(np.sqrt(1+2*classic_energy/c2)-1) From 537cc7481f8273b6d6acd0a778a5a63c09c92482 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Wed, 7 May 2025 10:07:46 +0200 Subject: [PATCH 26/28] Uniform precision requirements for all contributions --- one_electron.py | 11 +++++------ orbital4c/orbital.py | 12 ++++++++++++ 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/one_electron.py b/one_electron.py index c219bbb..bd8674e 100644 --- a/one_electron.py +++ b/one_electron.py @@ -44,7 +44,6 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): new_orbital.cropLargeSmall(prec) new_orbital.normalize() delta_psi = new_orbital - spinorb1 - #orbital_error = delta_psi.dot(delta_psi).real deltasq = delta_psi.squaredNorm() error_norm = np.sqrt(deltasq) print('Error', error_norm) @@ -95,11 +94,11 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): # while ((error_norm > thr or delta_e > prec/10) and idx < 100): while(idx < niter): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) - vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec*c2) + vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec) beta_v_psi = v_psi.beta2() - apV_psi = v_psi.alpha_p(prec*light_speed, derivative) - ap_psi = spinorb1.alpha_p(prec*light_speed, derivative) - Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec*light_speed) + apV_psi = v_psi.alpha_p(prec, derivative) + ap_psi = spinorb1.alpha_p(prec, derivative) + Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec) anticom = apV_psi + Vap_psi # anticom.cropLargeSmall(prec) # beta_v_psi.cropLargeSmall(prec) @@ -137,7 +136,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): energy_dirac = spinorb1.dot(add_psi).real beta_v_psi = v_psi.beta2() - ap_psi = spinorb1.alpha_p(prec*light_speed, derivative) + ap_psi = spinorb1.alpha_p(prec, derivative) cke = spinorb1.classicT() psi_beta_v_vpsi = spinorb1.dot(beta_v_psi).real diff --git a/orbital4c/orbital.py b/orbital4c/orbital.py index ed9179d..76ef70a 100644 --- a/orbital4c/orbital.py +++ b/orbital4c/orbital.py @@ -328,6 +328,18 @@ def dot(self, other): result += component return result +def print_expectation_value_VV(psi, V, prec): + v_psi = apply_potential(-1.0, V, psi, prec) + vv_psi = apply_potential(-1.0, V, v_psi, prec) + V2 = V * V + v2_psi = apply_potential(1.0, V2, psi, prec) + exp1 = v_psi.dot(v_psi).real + exp2 = psi.dot(vv_psi).real + exp3 = psi.dot(v2_psi).real + print(" = ", exp1) + print("< psi | vv_psi> = ", exp2) + print("< psi | v2_psi> = ", exp3) + def apply_dirac_hamiltonian(orbital, prec, shift = 0.0, der = 'ABGV'): beta_phi = orbital.beta(shift) grad_phi = orbital.gradient(der) From 70864eb6490b6062631bb08d07989b330679c297 Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Thu, 8 May 2025 12:17:11 +0200 Subject: [PATCH 27/28] Sample output for copilot --- sample_output.txt | 2544 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 2544 insertions(+) create mode 100644 sample_output.txt diff --git a/sample_output.txt b/sample_output.txt new file mode 100644 index 0000000..012b9b3 --- /dev/null +++ b/sample_output.txt @@ -0,0 +1,2544 @@ +BOX_SIZE: 5 +POT_TYPE: coulomb_HFYGB +TWOELECTRON: False +READORB: False +SAVEORB: False +ORBNAME: spinorb1 +P_PREC: 4 +P_THR: 3 +ORDER: 7 +ATOM: Ne +CHARGE: 10 +XPOS: 0.0 +YPOS: 0.0 +ZPOS: 0.0 +RADIUS: 0.0001 +EXPONENT: 357222173.0 +D1_OR_D2: False +READPOT: False +SAVEPOT: True +COMPPOT: True +DERIVATIVE: BS +INDEX: 0 +output_name: output-0.txt +10 [0.0, 0.0, 0.0] +Harrison potential +Potential + +*FunctionTree: nn + genNodes: 0 + square norm: 7674.11 + root scale: -2 + order: 7 + nodes: 1096 + endNodes: 960 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + scale=6 nodes=64 + scale=7 nodes=64 + scale=8 nodes=64 + scale=9 nodes=64 + scale=10 nodes=64 + scale=11 nodes=64 + scale=12 nodes=64 + scale=13 nodes=64 + scale=14 nodes=64 + scale=15 nodes=64 + +Using derivative BS +One-electron calculations 0.0001 +calc_dirac_mu 100.15201967928093 18728.722104819233 137.035999139 +Error 0.0011010164643532898 +Delta E 18728.722104819233 +Energy -50.142955204777536 +Converged? 0.0011010164643532898 > 0.001 ---- 18728.722104819233 > 1e-05 ---- iteration 0 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998669 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 4.01188e-18 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 4.01188e-18 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 4.01188e-18 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000443528 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000443528 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000443528 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00479649122863 18728.795913350295 137.035999139 +Error 0.0004278224120718559 +Delta E 0.07380853106224095 +Energy -50.069146673715295 +Converged? 0.0004278224120718559 > 0.001 ---- 0.07380853106224095 > 1e-05 ---- iteration 1 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998668 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 2.35429e-16 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 2.35429e-16 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 2.35429e-16 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000443932 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000443932 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000443932 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00811522075041 18728.79424954935 137.035999139 +Error 0.0001693486219423944 +Delta E 0.0016638009437883738 +Energy -50.07081047465908 +Converged? 0.0001693486219423944 > 0.001 ---- 0.0016638009437883738 > 1e-05 ---- iteration 2 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998668 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61946e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61946e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61946e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.00044415 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.00044415 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.00044415 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00818185731136 18728.794216141996 137.035999139 +Error 9.210295817906534e-05 +Delta E 3.340735565871e-05 +Energy -50.07084388201474 +Converged? 9.210295817906534e-05 > 0.001 ---- 3.340735565871e-05 > 1e-05 ---- iteration 3 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.58714e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.58714e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.58714e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444261 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444261 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444261 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00820480913998 18728.7942046354 137.035999139 +Error 4.884694620156261e-05 +Delta E 1.1506595910759643e-05 +Energy -50.07085538861065 +Converged? 4.884694620156261e-05 > 0.001 ---- 1.1506595910759643e-05 > 1e-05 ---- iteration 4 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61313e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61313e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61313e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444318 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444318 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444318 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00820986801382 18728.7942020992 137.035999139 +Error 2.5447606658260972e-05 +Delta E 2.5361987354699522e-06 +Energy -50.07085792480939 +Converged? 2.5447606658260972e-05 > 0.001 ---- 2.5361987354699522e-06 > 1e-05 ---- iteration 5 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61498e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61498e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61498e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444346 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444346 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444346 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821169364812 18728.794201183944 137.035999139 +Error 1.3096052785529658e-05 +Delta E 9.152572602033615e-07 +Energy -50.07085884006665 +Converged? 1.3096052785529658e-05 > 0.001 ---- 9.152572602033615e-07 > 1e-05 ---- iteration 6 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61558e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61558e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61558e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.00044436 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.00044436 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.00044436 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821245796678 18728.794200800763 137.035999139 +Error 6.683247346916536e-06 +Delta E 3.8318103179335594e-07 +Energy -50.07085922324768 +Converged? 6.683247346916536e-06 > 0.001 ---- 3.8318103179335594e-07 > 1e-05 ---- iteration 7 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61585e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61585e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61585e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444367 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444367 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444367 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821280479282 18728.794200626886 137.035999139 +Error 3.39116698156732e-06 +Delta E 1.7387719708494842e-07 +Energy -50.07085939712488 +Converged? 3.39116698156732e-06 > 0.001 ---- 1.7387719708494842e-07 > 1e-05 ---- iteration 8 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61599e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61599e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61599e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444371 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444371 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444371 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821296901069 18728.79420054456 137.035999139 +Error 1.7140566252757329e-06 +Delta E 8.232746040448546e-08 +Energy -50.07085947945234 +Converged? 1.7140566252757329e-06 > 0.001 ---- 8.232746040448546e-08 > 1e-05 ---- iteration 9 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61605e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61605e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61605e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444372 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444372 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444372 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821304864387 18728.794200504635 137.035999139 +Error 8.640968762614571e-07 +Delta E 3.992317942902446e-08 +Energy -50.07085951937552 +Converged? 8.640968762614571e-07 > 0.001 ---- 3.992317942902446e-08 > 1e-05 ---- iteration 10 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61609e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61609e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61609e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444373 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444373 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444373 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821308772568 18728.79420048504 137.035999139 +Error 4.34831546760069e-07 +Delta E 1.9594153854995966e-08 +Energy -50.07085953896967 +Converged? 4.34831546760069e-07 > 0.001 ---- 1.9594153854995966e-08 > 1e-05 ---- iteration 11 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.6161e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.6161e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.6161e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821310726182 18728.794200475248 137.035999139 +Error 2.1856870883680098e-07 +Delta E 9.793438948690891e-09 +Energy -50.07085954876311 +Converged? 2.1856870883680098e-07 > 0.001 ---- 9.793438948690891e-09 > 1e-05 ---- iteration 12 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61611e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61611e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61611e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.0082131169553 18728.794200470387 137.035999139 +Error 1.0977637237617293e-07 +Delta E 4.8603396862745285e-09 +Energy -50.07085955362345 +Converged? 1.0977637237617293e-07 > 0.001 ---- 4.8603396862745285e-09 > 1e-05 ---- iteration 13 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821312180204 18728.794200467957 137.035999139 +Error 5.5105909222448724e-08 +Delta E 2.4301698431372643e-09 +Energy -50.07085955605362 +Converged? 5.5105909222448724e-08 > 0.001 ---- 2.4301698431372643e-09 > 1e-05 ---- iteration 14 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821312427144 18728.79420046672 137.035999139 +Error 2.7652466464816207e-08 +Delta E 1.2369127944111824e-09 +Energy -50.07085955729053 +Converged? 2.7652466464816207e-08 > 0.001 ---- 1.2369127944111824e-09 > 1e-05 ---- iteration 15 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821312550296 18728.7942004661 137.035999139 +Error 1.3872841766913942e-08 +Delta E 6.184563972055912e-10 +Energy -50.07085955790899 +Converged? 1.3872841766913942e-08 > 0.001 ---- 6.184563972055912e-10 > 1e-05 ---- iteration 16 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821312598224 18728.79420046586 137.035999139 +Error 6.958392143350944e-09 +Delta E 2.4010660126805305e-10 +Energy -50.070859558149095 +Converged? 6.958392143350944e-09 > 0.001 ---- 2.4010660126805305e-10 > 1e-05 ---- iteration 17 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821312633774 18728.794200465683 137.035999139 +Error 3.4900941204868376e-09 +Delta E 1.7826096154749393e-10 +Energy -50.070859558327356 +Converged? 3.4900941204868376e-09 > 0.001 ---- 1.7826096154749393e-10 > 1e-05 ---- iteration 18 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +calc_dirac_mu 100.00821312640439 18728.79420046565 137.035999139 +Error 1.750145570944906e-09 +Delta E 3.2741809263825417e-11 +Energy -50.0708595583601 +Converged? 1.750145570944906e-09 > 0.001 ---- 3.2741809263825417e-11 > 1e-05 ---- iteration 19 +Large components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.998667 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 3.61612e-21 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + Small components + alpha +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0 + root scale: -2 + order: 7 + nodes: 8 + endNodes: 8 + nodes per scale: + scale=-2 nodes=8 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + beta +Real part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + + Imag part +*FunctionTree: nn + genNodes: 0 + square norm: 0.000444374 + root scale: -2 + order: 7 + nodes: 456 + endNodes: 400 + nodes per scale: + scale=-2 nodes=8 + scale=-1 nodes=64 + scale=0 nodes=64 + scale=1 nodes=64 + scale=2 nodes=64 + scale=3 nodes=64 + scale=4 nodes=64 + scale=5 nodes=64 + +Classic-like energies: cke = 50.52947961731555 cpe = -100.51588985100338 cke + cpe = -49.98641023368783 +Exact Energy = -50.06674201684291 +Dirac Energy = -50.070859558385564 +Kutze Energy = -50.053115926655934 +Error Kutze = 0.01362609018542571 +Error Dirac = -0.00411754154265509 +Delta Energy = -2.546585164964199e-11 +Dirac - Kutzelnigg = -0.0177436317280808 From 54c1bad94af2330c4ada97f76f05947f266492dc Mon Sep 17 00:00:00 2001 From: Luca Frediani Date: Thu, 8 May 2025 17:17:42 +0200 Subject: [PATCH 28/28] Reintroduces prec-based SCF stop --- one_electron.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/one_electron.py b/one_electron.py index bd8674e..b1b3163 100644 --- a/one_electron.py +++ b/one_electron.py @@ -27,8 +27,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): old_energy = 0 delta_e = 1 idx = 0 -# while ((error_norm > thr or delta_e > prec/10) and idx < 100): - while (idx < niter): + while (idx < niter and (delta_e > prec/100 or error_norm > thr)): hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = derivative) v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) add_psi = hd_psi + v_psi @@ -91,8 +90,7 @@ def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15): c2 = light_speed * light_speed old_energy = 0 idx = 0 -# while ((error_norm > thr or delta_e > prec/10) and idx < 100): - while(idx < niter): + while (idx < niter and (delta_e > prec/100 or error_norm > thr)): v_psi = orb.apply_potential(-1.0, potential, spinorb1, prec) vv_psi = orb.apply_potential(-0.5/c2, potential, v_psi, prec) beta_v_psi = v_psi.beta2()