From 27e09bb5b4b71b9fec6cbb948ee749c43c38ed17 Mon Sep 17 00:00:00 2001 From: Christian Tantardini <85489709+Christian48596@users.noreply.github.com> Date: Tue, 29 Aug 2023 14:54:55 +0200 Subject: [PATCH 1/6] implementation of general nuclear potential --- Z.txt | 3 + atom_list.txt | 3 + orbital4c/nuclear_potential.py | 1 - test.py | 104 ++++++++++++++++++++++----------- 4 files changed, 77 insertions(+), 34 deletions(-) create mode 100644 Z.txt create mode 100644 atom_list.txt diff --git a/Z.txt b/Z.txt new file mode 100644 index 0000000..0d2e6d5 --- /dev/null +++ b/Z.txt @@ -0,0 +1,3 @@ +H 1 +He 2 +Pu 94 \ No newline at end of file diff --git a/atom_list.txt b/atom_list.txt new file mode 100644 index 0000000..482e6eb --- /dev/null +++ b/atom_list.txt @@ -0,0 +1,3 @@ + H 0.1 0.2 0.3 + H 0.2 0.4 0.3 +Pu 0.3 0.6 0.2 \ No newline at end of file diff --git a/orbital4c/nuclear_potential.py b/orbital4c/nuclear_potential.py index c7f7d7d..3a6d288 100644 --- a/orbital4c/nuclear_potential.py +++ b/orbital4c/nuclear_potential.py @@ -2,7 +2,6 @@ import copy as cp from scipy.special import erf - def point_charge(position, center , charge): d2 = ((position[0] - center[0])**2 + (position[1] - center[1])**2 + diff --git a/test.py b/test.py index d1026fa..3c6b105 100644 --- a/test.py +++ b/test.py @@ -22,20 +22,20 @@ if __name__ == '__main__': parser = argparse.ArgumentParser(description='Collecting all data tostart the program.') - parser.add_argument('-a', '--atype', dest='atype', type=str, default='He', - help='put the atom type') + #parser.add_argument('-a', '--atype', dest='atype', type=str, default='He', + # help='put the atom type') parser.add_argument('-d', '--derivative', dest='deriv', type=str, default='ABGV', help='put the type of derivative') - parser.add_argument('-z', '--charge', dest='charge', type=float, default=2.0, - help='put the atom charge') + #parser.add_argument('-z', '--charge', dest='charge', type=float, default=2.0, + # help='put the atom charge') parser.add_argument('-b', '--box', dest='box', type=int, default=60, help='put the box dimension') - parser.add_argument('-cx', '--center_x', dest='cx', type=float, default=0.0, - help='position of nucleus in x') - parser.add_argument('-cy', '--center_y', dest='cy', type=float, default=0.0, - help='position of nucleus in y') - parser.add_argument('-cz', '--center_z', dest='cz', type=float, default=0.0, - help='position of nucleus in z') + #parser.add_argument('-cx', '--center_x', dest='cx', type=float, default=0.0, + # help='position of nucleus in x') + #parser.add_argument('-cy', '--center_y', dest='cy', type=float, default=0.0, + # help='position of nucleus in y') + #parser.add_argument('-cz', '--center_z', dest='cz', type=float, default=0.0, + # help='position of nucleus in z') 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, @@ -48,9 +48,9 @@ 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.atype != 'H', 'Please consider only atoms with more than one electron' + #assert args.atype != 'H', 'Please consider only atoms with more than one electron' - assert args.charge > 1.0, 'Please consider only atoms with more than one electron' + #assert args.charge > 1.0, 'Please consider only atoms with more than one electron' assert args.coulgau in ['coulomb', 'gaunt', 'gaunt-test'], 'Please, specify coulgau in a rigth way: coulomb or gaunt' @@ -66,18 +66,46 @@ l = 0 n = 1 m = 0.5 - Z = args.charge - atom = args.atype + + def read_file_with_named_lists(atomlist): + atom_lists = {} + + with open(atomlist, 'r') as file: + for line in file: + terms = line.strip().split() + atom = terms[0] + origin = terms[1:] + origin = [float(element) for element in origin] + + if atom in atom_lists: + # Append an identifier to make the key unique + identifier = len(atom_lists[atom]) + 1 + unique_key = f"{atom}_{identifier}" + atom_lists[unique_key] = origin + else: + atom_lists[atom] = origin + + return atom_lists + + def get_original_list_name(key): + return key.split('_')[0] + + atomlist = './atom_list.txt' # Replace with the actual file name + coordinates = read_file_with_named_lists(atomlist) + + #origin = [args.cx, args.cy, args.cz] + #Z = args.charge + #atom = args.atype ################# Call MRA ####################### mra = vp.MultiResolutionAnalysis(box=[-args.box,args.box], order=args.order, max_depth = 30) prec = args.prec - origin = [args.cx, args.cy, args.cz] - orb.orbital4c.mra = mra orb.orbital4c.light_speed = light_speed cf.complex_fcn.mra = mra + default_der = args.deriv print('call MRA DONE') - + + ################## Jobs ########################## computeNuclearPotential = True readOrbitals = True runCoulomb = False @@ -90,26 +118,36 @@ runGaugeB = False runGaugeC = False runGaugeD = False - runGaugeDelta = False - default_der = args.deriv + runGaugeDelta = False + print('Jobs chosen') + ################### Define V potential ###################### if(computeNuclearPotential): - if args.potential == 'point_charge': - Peps = vp.ScalingProjector(mra,prec/10) - f = lambda x: nucpot.point_charge(x, origin, Z) - V_tree = Peps(f) - elif args.potential == 'coulomb_HFYGB': + V_tot = vp.FunctionTree(mra) + V_tot.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 + f = lambda x: nucpot.point_charge(x, origin, charge) Peps = vp.ScalingProjector(mra,prec/10) - f = lambda x: nucpot.coulomb_HFYGB(x, origin, Z, prec) - V_tree = Peps(f) - elif args.potential == 'homogeneus_charge_sphere': - Peps = vp.ScalingProjector(mra,prec/10) - f = lambda x: nucpot.homogeneus_charge_sphere(x, origin, Z, atom) - V_tree = Peps(f) - elif args.potential == 'gaussian': - Peps = vp.ScalingProjector(mra,prec/10) - f = lambda x: nucpot.gaussian(x, origin, Z, atom) V_tree = Peps(f) + V_tot += V_tree + + #Peps = vp.ScalingProjector(mra,prec/10) + #f = lambda x: nucpot.point_charge(x, origin, charge) + #V_tree = Peps(f) print('Define V Potential', args.potential, 'DONE') #############################START WITH CALCULATION################################### From b009e2d86d470b25d71be6eb20798da6a6e6431f Mon Sep 17 00:00:00 2001 From: Christian Tantardini <85489709+Christian48596@users.noreply.github.com> Date: Tue, 29 Aug 2023 15:00:09 +0200 Subject: [PATCH 2/6] possibility to have many nucleus with different nuclear model charge --- test.py | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/test.py b/test.py index 3c6b105..76052ad 100644 --- a/test.py +++ b/test.py @@ -140,14 +140,23 @@ def get_original_list_name(key): fileObj.close() print("Origin:", origin) print() # Print an empty line for separation - f = lambda x: nucpot.point_charge(x, origin, charge) - Peps = vp.ScalingProjector(mra,prec/10) - V_tree = Peps(f) + if args.potential == 'point_charge': + Peps = vp.ScalingProjector(mra,prec/10) + f = lambda x: nucpot.point_charge(x, origin, Z) + V_tree = Peps(f) + elif args.potential == 'coulomb_HFYGB': + Peps = vp.ScalingProjector(mra,prec/10) + f = lambda x: nucpot.coulomb_HFYGB(x, origin, Z, prec) + V_tree = Peps(f) + elif args.potential == 'homogeneus_charge_sphere': + Peps = vp.ScalingProjector(mra,prec/10) + f = lambda x: nucpot.homogeneus_charge_sphere(x, origin, Z, + V_tree = Peps(f) + elif args.potential == 'gaussian': + Peps = vp.ScalingProjector(mra,prec/10) + f = lambda x: nucpot.gaussian(x, origin, Z, atom) + V_tree = Peps(f) V_tot += V_tree - - #Peps = vp.ScalingProjector(mra,prec/10) - #f = lambda x: nucpot.point_charge(x, origin, charge) - #V_tree = Peps(f) print('Define V Potential', args.potential, 'DONE') #############################START WITH CALCULATION################################### From 667665bd3f4e8014618c34dcdf9536f44065d1b3 Mon Sep 17 00:00:00 2001 From: Christian Tantardini <85489709+Christian48596@users.noreply.github.com> Date: Wed, 30 Aug 2023 17:51:45 +0200 Subject: [PATCH 3/6] one_electron file to compute D and shadows of D2 for one electron --- one_electron.py | 151 ++++++++++++++++++++++++++++++++++++++++++++++++ test.py | 17 ++++-- 2 files changed, 163 insertions(+), 5 deletions(-) create mode 100644 one_electron.py diff --git a/one_electron.py b/one_electron.py new file mode 100644 index 0000000..69a0d5e --- /dev/null +++ b/one_electron.py @@ -0,0 +1,151 @@ +from orbital4c import complex_fcn as cf +from orbital4c import orbital as orb +from orbital4c import operators as oper +from scipy.constants import hbar +from scipy.linalg import eig, inv +from scipy.special import legendre, laguerre, erf, gamma +from vampyr import vampyr3d as vp +import numpy as np +import numpy.linalg as LA +import sys, getopt + +def gs_D_1e(spinorb1, potential, mra, prec, der): + print('Hartree-Fock 1e') + + error_norm = 1 + compute_last_energy = False + + P = vp.PoissonOperator(mra, prec) + light_speed = spinorb1.light_speed + + while (error_norm > prec or compute_last_energy): + n_22 = spinorb1.overlap_density(spinorb1, prec) + + # Definition of two electron operators + B22 = P(n_22.real) * (4 * np.pi) + + # Definiton of Dirac Hamiltonian for spinorb 1 + hd_psi_1 = orb.apply_dirac_hamiltonian(spinorb1, prec, 0.0, der) + hd_11 = spinorb1.dot(hd_psi_1) + print("hd_11", hd_11) + + # Applying nuclear potential to spinorb 1 + v_psi_1 = orb.apply_potential(-1.0, potential, spinorb1, prec) + V1 = spinorb1.dot(v_psi_1) + + hd_V_11 = hd_11 + V1 + + eps = hd_V_11.real + + print('orbital energy', eps - light_speed**2) + + if(compute_last_energy): + break + + mu = orb.calc_dirac_mu(eps, light_speed) + tmp = orb.apply_helmholtz(V1, mu, prec) + new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, eps, der) + new_orbital *= 0.5/light_speed**2 + print("============= Spinor before Helmholtz =============") + print(spinorb1) + print("============= New spinor before crop =============") + print(new_orbital) + new_orbital.normalize() + new_orbital.cropLargeSmall(prec) + + # Compute orbital error + delta_psi = new_orbital - spinorb1 + deltasq = delta_psi.squaredNorm() + error_norm = np.sqrt(deltasq) + print('Orbital_Error norm', error_norm) + spinorb1 = new_orbital + if(error_norm < prec): + compute_last_energy = True + return spinorb1 + + +def gs_1e_D2(spinorb1, potential, mra, prec, der): + print('Hartree-Fock 1e D2') + + error_norm = 1.0 + compute_last_energy = False + + P = vp.PoissonOperator(mra, prec) + + light_speed = spinorb1.light_speed + c2 = light_speed**2 + + while(error_norm > prec): + print("Applying V") + Vpsi = orb.apply_potential(-1.0, potential, spinorb1, prec) + VVpsi = orb.apply_potential(-0.5/mc2, potential, Vpsi, prec) + beta_Vpsi = Vpsi.beta2() + apV_psi = Vpsi.alpha_p(prec, der) + ap_psi = spinorb1.alpha_p(prec, der) + Vap_psi = orb.apply_potential(-1.0, potential, ap_psi, prec) + anticom = apV_psi + Vap_psi + anticom *= 1.0 / (2.0 * light_speed) + RHS = beta_Vpsi + VVpsi + anticom * (0.5/light_speed) + + anticom.cropLargeSmall(prec) + beta_Vpsi.cropLargeSmall(prec) + VVpsi.cropLargeSmall(prec) + + print("anticom") + print(anticom) + print("beta_Vpsi") + print(beta_Vpsi) + print("VV_psi") + print(VVpsi) + RHS = beta_Vpsi + anticom + VVpsi + RHS.cropLargeSmall(prec) + print("RHS") + print(RHS) + + cke = spinorb1.classicT() + cpe = (spinorb1.dot(RHS)).real + print("Classic-like energies: ", cke, cpe, cke + cpe) + print("Orbital energy: ", c2 * ( -1.0 + np.sqrt(1 + 2 * (cpe + cke) / c2))) + mu = orb.calc_non_rel_mu(cke+cpe) + print("mu: ", mu) + + + new_spinorb1 = orb.apply_helmholtz(RHS, mu, prec) + print("normalization") + new_spinorb1.normalize() + print("crop") + new_spinorb1.cropLargeSmall(prec) + new_spinorb1.append(new_spinorb1) + new_spinorb1.append(new_spinorb1.ktrs(prec)) + + # Compute orbital error + delta_psi = new_spinorb1 - spinorb1 + deltasq = delta_psi.squaredNorm() + error_norm = np.sqrt(deltasq) + print('Orbital_Error norm', error_norm) + spinorb1 = new_spinorb1 + + hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der) + Vpsi = orb.apply_potential(-1.0, potential, spinor, prec) + add_psi = hd_psi + Vpsi + energy = (spinor.dot(add_psi)).real + + cke = spinorb1.classicT() + beta_Vpsi = Vpsi.beta2() + beta_pot = (beta_Vpsi.dot(spinorb1)).real + pot_sq = (Vpsi.dot(Vpsi)).real + ap_psi = spinorb1.alpha_p(prec, der) + anticom = (ap_psi.dot(Vpsi)).real + energy_kutzelnigg = cke + beta_pot + pot_sq/(2*mc2) + anticom/light_speed + + print('Kutzelnigg',cke, beta_pot, pot_sq/(2*mc2), anticom/light_speed, energy_kutzelnigg) + print('Quadratic approx',energy_kutzelnigg - energy_kutzelnigg**2/(2*mc2)) + print('Correct from Kutzelnigg', mc2*(np.sqrt(1+2*energy_kutzelnigg/mc2)-1)) + print('Final Energy',energy - light_speed**2) + + energy_1s = analytic_1s(light_speed, n, k, Z) + + print('Exact Energy',energy_1s - light_speed**2) + print('Difference 1',energy_1s - energy) + print('Difference 2',energy_1s - energy_kutzelnigg - light_speed**2) + return spinorb1 diff --git a/test.py b/test.py index 76052ad..6026797 100644 --- a/test.py +++ b/test.py @@ -108,7 +108,8 @@ def get_original_list_name(key): ################## Jobs ########################## computeNuclearPotential = True readOrbitals = True - runCoulomb = False + runCoulomb1e = False + runCoulomb2e = False runCoulombGen = False runKutzelnigg = False runKutzSimple = True @@ -184,18 +185,24 @@ def get_original_list_name(key): length = 2 * args.box + if runD_1e: + spinorb1 = one_electron.gs_D_1e(spinorb1, V_tree, mra, prec) + + if runD2_1e: + spinorb1 = one_electron.gs_1e_D2(spinorb1, V_tree, mra, prec) + if runCoulombGen: spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec) - + + if runCoulomb2e: + spinorb1, spinorb2 = two_electron.coulomb_gs_2e(spinorb1, V_tree, mra, prec) + if runKutzelnigg: spinorb1, spinorb2 = two_electron.coulomb_2e_D2([spinorb1, spinorb2], V_tree, mra, prec, 'ABGV') if runKutzSimple: spinorb1, spinorb2 = two_electron.coulomb_2e_D2_J([spinorb1, spinorb2], V_tree, mra, prec, der = 'ABGV') - if runCoulomb: - spinorb1, spinorb2 = two_electron.coulomb_gs_2e(spinorb1, V_tree, mra, prec) - if runGaunt: two_electron.calcGauntPert(spinorb1, spinorb2, mra, prec) From fc6800813590512b911497771596c678d0eddd0b Mon Sep 17 00:00:00 2001 From: Christian Tantardini <85489709+Christian48596@users.noreply.github.com> Date: Fri, 1 Sep 2023 15:15:40 +0200 Subject: [PATCH 4/6] reading of atoms creating a center of mass and a spinorbital guess from the combination of gaussians --- atom_list.txt | 5 +- one_electron.py | 11 ++- orbital4c/nuclear_potential.py | 91 ++++++++++++++++++++++- test.py | 130 ++++++++++----------------------- two_electron.py | 2 +- 5 files changed, 136 insertions(+), 103 deletions(-) diff --git a/atom_list.txt b/atom_list.txt index 482e6eb..6b01b2b 100644 --- a/atom_list.txt +++ b/atom_list.txt @@ -1,3 +1,2 @@ - H 0.1 0.2 0.3 - H 0.2 0.4 0.3 -Pu 0.3 0.6 0.2 \ No newline at end of file + H 0.1 0.2 -0.7 + H 0.1 0.2 1.3 diff --git a/one_electron.py b/one_electron.py index 69a0d5e..2ef4e6a 100644 --- a/one_electron.py +++ b/one_electron.py @@ -19,10 +19,10 @@ def gs_D_1e(spinorb1, potential, mra, prec, der): light_speed = spinorb1.light_speed while (error_norm > prec or compute_last_energy): - n_22 = spinorb1.overlap_density(spinorb1, prec) + n_11 = spinorb1.overlap_density(spinorb1, prec) # Definition of two electron operators - B22 = P(n_22.real) * (4 * np.pi) + B11 = P(n_11.real) * (4 * np.pi) # Definiton of Dirac Hamiltonian for spinorb 1 hd_psi_1 = orb.apply_dirac_hamiltonian(spinorb1, prec, 0.0, der) @@ -43,7 +43,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, der): break mu = orb.calc_dirac_mu(eps, light_speed) - tmp = orb.apply_helmholtz(V1, mu, prec) + tmp = orb.apply_helmholtz(v_psi_1, mu, prec) new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, eps, der) new_orbital *= 0.5/light_speed**2 print("============= Spinor before Helmholtz =============") @@ -64,7 +64,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, der): return spinorb1 -def gs_1e_D2(spinorb1, potential, mra, prec, der): +def gs_D2_1e(spinorb1, potential, mra, prec, der): print('Hartree-Fock 1e D2') error_norm = 1.0 @@ -115,8 +115,7 @@ def gs_1e_D2(spinorb1, potential, mra, prec, der): new_spinorb1.normalize() print("crop") new_spinorb1.cropLargeSmall(prec) - new_spinorb1.append(new_spinorb1) - new_spinorb1.append(new_spinorb1.ktrs(prec)) + # Compute orbital error delta_psi = new_spinorb1 - spinorb1 diff --git a/orbital4c/nuclear_potential.py b/orbital4c/nuclear_potential.py index 3a6d288..c8fb9c9 100644 --- a/orbital4c/nuclear_potential.py +++ b/orbital4c/nuclear_potential.py @@ -1,6 +1,93 @@ +from scipy.constants import hbar +from scipy.linalg import eig, inv +from scipy.special import legendre, laguerre, erf, gamma +from scipy.special import gamma +from vampyr import vampyr3d as vp +from vampyr import vampyr1d as vp1 + +import argparse import numpy as np -import copy as cp -from scipy.special import erf +import numpy.linalg as LA +import sys, getopt + +def read_file_with_named_lists(atomlist): + atom_lists = {} + + with open(atomlist, 'r') as file: + for line in file: + terms = line.strip().split() + atom = terms[0] + origin = terms[1:] + origin = [float(element) for element in origin] + + if atom in atom_lists: + # Append an identifier to make the key unique + identifier = len(atom_lists[atom]) + 1 + unique_key = f"{atom}_{identifier}" + atom_lists[unique_key] = origin + else: + atom_lists[atom] = origin + total_atom_lists = len(atom_lists) + return atom_lists, total_atom_lists + +def get_original_list_name(key): + return key.split('_')[0] + + +def calculate_center_of_mass(coordinates): + total_mass = 0.0 + center_of_mass = [0.0, 0.0, 0.0] + + for atom, origin in coordinates.items(): + # Assuming each atom has mass 1.0 (modify if necessary) + mass = 1.0 + total_mass += mass + + # Update the center of mass coordinates + for i in range(3): + center_of_mass[i] += origin[i] * mass + + # Calculate the weighted average to get the center of mass + for i in range(3): + center_of_mass[i] /= total_mass + + return center_of_mass + + +def pot(V_tree, coordinates, typenuc, mra, prec, der): + for atom, origin in coordinates.items(): + atom = get_original_list_name(atom) + print("Atom:", atom) + fileObj = open("/Users/cta018/vampyr-dev/ReMRChem/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) + V_tree += V + print('Define V Potential', typenuc, 'DONE') def point_charge(position, center , charge): d2 = ((position[0] - center[0])**2 + diff --git a/test.py b/test.py index 6026797..dca4009 100644 --- a/test.py +++ b/test.py @@ -15,6 +15,7 @@ import numpy.linalg as LA import sys, getopt +import one_electron import two_electron import importlib @@ -22,20 +23,10 @@ if __name__ == '__main__': parser = argparse.ArgumentParser(description='Collecting all data tostart the program.') - #parser.add_argument('-a', '--atype', dest='atype', type=str, default='He', - # help='put the atom type') parser.add_argument('-d', '--derivative', dest='deriv', type=str, default='ABGV', help='put the type of derivative') - #parser.add_argument('-z', '--charge', dest='charge', type=float, default=2.0, - # help='put the atom charge') parser.add_argument('-b', '--box', dest='box', type=int, default=60, help='put the box dimension') - #parser.add_argument('-cx', '--center_x', dest='cx', type=float, default=0.0, - # help='position of nucleus in x') - #parser.add_argument('-cy', '--center_y', dest='cy', type=float, default=0.0, - # help='position of nucleus in y') - #parser.add_argument('-cz', '--center_z', dest='cz', type=float, default=0.0, - # help='position of nucleus in z') 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, @@ -48,10 +39,6 @@ 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.atype != 'H', 'Please consider only atoms with more than one electron' - - #assert args.charge > 1.0, 'Please consider only atoms with more than one electron' - assert args.coulgau in ['coulomb', 'gaunt', 'gaunt-test'], 'Please, specify coulgau in a rigth way: coulomb or gaunt' assert args.potential in ['point_charge', 'smoothing_HFYGB', 'coulomb_HFYGB', 'homogeneus_charge_sphere', 'gaussian'], 'Please, specify V' @@ -67,35 +54,6 @@ n = 1 m = 0.5 - def read_file_with_named_lists(atomlist): - atom_lists = {} - - with open(atomlist, 'r') as file: - for line in file: - terms = line.strip().split() - atom = terms[0] - origin = terms[1:] - origin = [float(element) for element in origin] - - if atom in atom_lists: - # Append an identifier to make the key unique - identifier = len(atom_lists[atom]) + 1 - unique_key = f"{atom}_{identifier}" - atom_lists[unique_key] = origin - else: - atom_lists[atom] = origin - - return atom_lists - - def get_original_list_name(key): - return key.split('_')[0] - - atomlist = './atom_list.txt' # Replace with the actual file name - coordinates = read_file_with_named_lists(atomlist) - - #origin = [args.cx, args.cy, args.cz] - #Z = args.charge - #atom = args.atype ################# Call MRA ####################### mra = vp.MultiResolutionAnalysis(box=[-args.box,args.box], order=args.order, max_depth = 30) prec = args.prec @@ -107,12 +65,13 @@ def get_original_list_name(key): ################## Jobs ########################## computeNuclearPotential = True - readOrbitals = True - runCoulomb1e = False - runCoulomb2e = False + readOrbitals = False + runD_1e = True + runD2_1e = False runCoulombGen = False + runCoulomb2e = False runKutzelnigg = False - runKutzSimple = True + runKutzSimple = False saveOrbitals = False runGaunt = False runGaugeA = False @@ -121,44 +80,27 @@ def get_original_list_name(key): runGaugeD = False runGaugeDelta = False print('Jobs chosen') - + + + ################### Reading Atoms ######################### + atomlist = '/Users/cta018/vampyr-dev/ReMRChem/atom_list.txt' # Replace with the actual file name + coordinates, total_atom_lists = nucpot.read_file_with_named_lists(atomlist) + ################### Define V potential ###################### if(computeNuclearPotential): - V_tot = vp.FunctionTree(mra) - V_tot.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 args.potential == 'point_charge': - Peps = vp.ScalingProjector(mra,prec/10) - f = lambda x: nucpot.point_charge(x, origin, Z) - V_tree = Peps(f) - elif args.potential == 'coulomb_HFYGB': - Peps = vp.ScalingProjector(mra,prec/10) - f = lambda x: nucpot.coulomb_HFYGB(x, origin, Z, prec) - V_tree = Peps(f) - elif args.potential == 'homogeneus_charge_sphere': - Peps = vp.ScalingProjector(mra,prec/10) - f = lambda x: nucpot.homogeneus_charge_sphere(x, origin, Z, - V_tree = Peps(f) - elif args.potential == 'gaussian': - Peps = vp.ScalingProjector(mra,prec/10) - f = lambda x: nucpot.gaussian(x, origin, Z, atom) - V_tree = Peps(f) - V_tot += V_tree - print('Define V Potential', args.potential, 'DONE') + V_tree = vp.FunctionTree(mra) + V_tree.setZero() + typenuc = args.potential + V_tree = nucpot.pot(V_tree, coordinates, typenuc, mra, prec, der= 'BS') + + ################### Define Center of Mass ################### + if total_atom_lists >= 2: + # Calculate the center of mass + com_coordinates = nucpot.calculate_center_of_mass(coordinates) + print("Center of Mass (x, y, z):", com_coordinates) + else: + print("There are not enough atoms (less than 2) to calculate the center of mass.") + com_coordinates = coordinates #############################START WITH CALCULATION################################### @@ -168,15 +110,21 @@ def get_original_list_name(key): spinorb1.read("spinorb1") spinorb2.read("spinorb2") else: + gauss_tree_tot = vp.FunctionTree(mra) + gauss_tree_tot.setZero() a_coeff = 3.0 b_coeff = np.sqrt(a_coeff/np.pi)**3 - gauss = vp.GaussFunc(b_coeff, a_coeff, origin) - gauss_tree = vp.FunctionTree(mra) - vp.advanced.build_grid(out=gauss_tree, inp=gauss) - vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss) - gauss_tree.normalize() + for atom, origin in coordinates.items(): + gauss = vp.GaussFunc(b_coeff, a_coeff, origin) + gauss_tree = vp.FunctionTree(mra) + vp.advanced.build_grid(out=gauss_tree, inp=gauss) + vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss) + gauss_tree.normalize() + gauss_tree_tot += gauss_tree + + gauss_tree_tot.normalize() complexfc = cf.complex_fcn() - complexfc.copy_fcns(real=gauss_tree) + complexfc.copy_fcns(real=gauss_tree_tot) spinorb1.copy_components(La=complexfc) spinorb1.init_small_components(prec/10) spinorb1.normalize() @@ -186,10 +134,10 @@ def get_original_list_name(key): length = 2 * args.box if runD_1e: - spinorb1 = one_electron.gs_D_1e(spinorb1, V_tree, mra, prec) + spinorb1 = one_electron.gs_D_1e(spinorb1, V_tree, mra, prec, der= 'BS') if runD2_1e: - spinorb1 = one_electron.gs_1e_D2(spinorb1, V_tree, mra, prec) + spinorb1 = one_electron.gs_D2_1e(spinorb1, V_tree, mra, prec, der= 'BS') if runCoulombGen: spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], V_tree, mra, prec) diff --git a/two_electron.py b/two_electron.py index 3609c4f..9614955 100644 --- a/two_electron.py +++ b/two_electron.py @@ -645,7 +645,7 @@ def build_RHS_D2(Jop, Vop, spinor, prec, light_speed): ap_psi = spinor.alpha_p(prec) VT_ap_psi = 0.5 * Jop(ap_psi) - Vop(ap_psi) anticom = VT_ap_psi + ap_VT_psi - anticom *= 1.0 / (2.0 * light_speed + anticom *= 1.0 / (2.0 * light_speed) anticom.cropLargeSmall(prec) VT_VT_psi = 0.5 * Jop(VT_psi) - Vop(VT_psi) From 05cccc3bfe25f81275fafdfc61c319a5f0e599cc Mon Sep 17 00:00:00 2001 From: Christian Tantardini <85489709+Christian48596@users.noreply.github.com> Date: Fri, 1 Sep 2023 16:17:28 +0200 Subject: [PATCH 5/6] apply_potential does not work --- one_electron.py | 7 ++++--- orbital4c/operators.py | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/one_electron.py b/one_electron.py index 2ef4e6a..708a92c 100644 --- a/one_electron.py +++ b/one_electron.py @@ -16,6 +16,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, der): compute_last_energy = False P = vp.PoissonOperator(mra, prec) + #Vop = oper.PotentialOperator(mra, prec, potential) light_speed = spinorb1.light_speed while (error_norm > prec or compute_last_energy): @@ -30,8 +31,8 @@ def gs_D_1e(spinorb1, potential, mra, prec, der): print("hd_11", hd_11) # Applying nuclear potential to spinorb 1 - v_psi_1 = orb.apply_potential(-1.0, potential, spinorb1, prec) - V1 = spinorb1.dot(v_psi_1) + Vpsi1 = orb.apply_potential(-1.0, potential, spinorb1, prec) + V1 = spinorb1.dot(Vpsi1) hd_V_11 = hd_11 + V1 @@ -43,7 +44,7 @@ def gs_D_1e(spinorb1, potential, mra, prec, der): break mu = orb.calc_dirac_mu(eps, light_speed) - tmp = orb.apply_helmholtz(v_psi_1, mu, prec) + tmp = orb.apply_helmholtz(Vpsi1, mu, prec) new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, eps, der) new_orbital *= 0.5/light_speed**2 print("============= Spinor before Helmholtz =============") diff --git a/orbital4c/operators.py b/orbital4c/operators.py index cfaaa0c..d7a6ad3 100644 --- a/orbital4c/operators.py +++ b/orbital4c/operators.py @@ -134,7 +134,7 @@ def __call__(self, phi): if(self.real): result = orb.apply_potential(1.0, self.potential, phi, self.prec) 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) result.cropLargeSmall(self.prec) return result From b6c21a32d6e9fe5301ef2c70163911fff7986f5e Mon Sep 17 00:00:00 2001 From: Christian Tantardini <85489709+Christian48596@users.noreply.github.com> Date: Tue, 5 Sep 2023 16:47:16 +0200 Subject: [PATCH 6/6] removed atom_list.txt from chache --- atom_list.txt | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 atom_list.txt diff --git a/atom_list.txt b/atom_list.txt deleted file mode 100644 index 6b01b2b..0000000 --- a/atom_list.txt +++ /dev/null @@ -1,2 +0,0 @@ - H 0.1 0.2 -0.7 - H 0.1 0.2 1.3