Merge pull request #52 from Christian48596/1e2e

implementation of general nuclear potential

Z.txt

@@ -0,0 +1,3 @@
+H   1
+He  2
+Pu 94

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)
+    #Vop = oper.PotentialOperator(mra, prec, potential)
+    light_speed = spinorb1.light_speed
+
+    while (error_norm > prec or compute_last_energy):
+        n_11 = spinorb1.overlap_density(spinorb1, prec)
+
+        # Definition of two electron operators
+        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)
+        hd_11 = spinorb1.dot(hd_psi_1)
+        print("hd_11", hd_11)
+
+        # Applying nuclear potential to spinorb 1 
+        Vpsi1 = orb.apply_potential(-1.0, potential, spinorb1, prec)
+        V1 = spinorb1.dot(Vpsi1)
+
+        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(Vpsi1, 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_D2_1e(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)
+
+
+        # 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

orbital4c/nuclear_potential.py

@@ -1,7 +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 +

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