Input

.gitignore

@@ -136,3 +136,5 @@ dmypy.json
 *.out
 *.prof
 *.txt
+*.err
+*.tree

one_electron.py

@@ -17,94 +17,140 @@ 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, charge, niter=15):
+    print('One-electron calculations', prec)
 
     error_norm = 1
-    #compute_last_energy = False
 
     light_speed = spinorb1.light_speed
-
-    while error_norm > prec:
+    c2 = light_speed**2
+    old_energy = 0
+    delta_e = 1
+    idx = 0
+    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
         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)
+        if(idx > 10):
+            new_orbital = new_orbital + spinorb1
+        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)
+        delta_e = np.abs(energy - old_energy)
+        print('Delta E', delta_e)
+        print('Energy',energy - light_speed**2)
+        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 = spinorb1.dot(add_psi).real
-    print('Final Energy',energy - light_speed**2)
+    energy_1s = analytic_1s(light_speed, 1, -1, charge)
+
+    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 + 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)
+
+    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
 
 
-def gs_D2_1e(spinorb1, potential, mra, prec, derivative):
+def gs_D2_1e(spinorb1, potential, mra, prec, thr, derivative, charge, niter=15):
     print('Hartree-Fock 1e D2')
     error_norm = 1
+    delta_e = 1
     light_speed = spinorb1.light_speed
     c2 = light_speed * light_speed
-
-    while error_norm > prec:
+    old_energy = 0
+    idx = 0
+    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()
         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)
-        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)
+        if(idx > 10):
+            new_orbital = new_orbital + spinorb1
+        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)
-        spinorb1 = new_orbital
+        delta_e = np.abs(energy - old_energy)
+        print('Delta E', delta_e)
+        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 = spinorb1.dot(add_psi).real
+    energy_dirac = spinorb1.dot(add_psi).real
 
-    cke = spinorb1.classicT()
     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)
-
-    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)
+    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 + 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)
+
+    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)
+    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

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]]
@@ -116,17 +51,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
@@ -161,22 +98,11 @@ 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)
     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

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 = 1):
         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