Merge pull request #54 from Christian48596/1e2e

1e2e

H-square.py

@@ -138,7 +138,7 @@ def analytic_1s(light_speed, n, k, Z):
     new_orbital.normalize()
     delta_psi = new_orbital - spinor_H
     orbital_error = (delta_psi.dot(delta_psi)).real
-    print('Error',orbital_error, imag, flush = True)
+    print('Error',orbital_error)
     spinor_H = new_orbital
 
 hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = 'BS')

H2+.py

@@ -15,13 +15,14 @@
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description='Collecting all data tostart the program.')
-    parser.add_argument('-a', '--atype', dest='atype', type=str, default='H',
-                        help='put the atom type')
+    parser.add_argument('-d', '--dtype', dest='dtype', type=str, default='dirac',
+                        help='Dirac or Dirac-square operators')
     parser.add_argument('-v', '--potential', dest='potential', type=str, default='point_charge',
                         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.potential in ['point_charge', 'smoothing_HFYGB', 'coulomb_HFYGB', 'homogeneus_charge_sphere', 'gaussian'], 'Please, specify V'
+    assert args.dtype in ['dirac', 'dirac2'], 'Please, specify Dirac-type operator'
 
 def analytic_1s(light_speed, n, k, Z):
     alpha = 1/light_speed
@@ -39,10 +40,10 @@ def analytic_1s(light_speed, n, k, Z):
 n = 1
 m = 0.5
 Z = 1
-atom = args.atype
+atom = 'H'
 
 energy_1s = analytic_1s(light_speed, n, k, Z)
-print('Exact Energy',energy_1s - light_speed**2, flush = True)
+print('Exact Energy',energy_1s - light_speed**2)
 
 mra = vp.MultiResolutionAnalysis(box=[-60,60], order=6)
 prec = 1.0e-4
@@ -77,6 +78,7 @@ def VH2(x, origin1, origin2, Z1, Z2):
 
 Peps = vp.ScalingProjector(mra,prec/10)
 V_tree = Peps(f)
+print('V_tree', V_tree)
 print('Define V Potential', args.potential, 'DONE')
 
 orb.orbital4c.light_speed = light_speed
@@ -105,75 +107,79 @@ def VH2(x, origin1, origin2, Z1, Z2):
 spinor_H.normalize()
 
 
-default_der = 'BS'
-
-orbital_error = 1
-mc2 = light_speed * light_speed
-
-while orbital_error > prec:
-    v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec) 
-    vv_psi = orb.apply_potential(-0.5/mc2, V_tree, v_psi, prec)
-    beta_v_psi = v_psi.beta2()
-    apV_psi = v_psi.alpha_p(prec, 'BS')
-    ap_psi = spinor_H.alpha_p(prec, 'BS')
-    Vap_psi = orb.apply_potential(-1.0, V_tree, ap_psi, prec)
-    anticom = apV_psi + Vap_psi
-    RHS = beta_v_psi + vv_psi + anticom * (0.5/light_speed)
+defult_der = 'BS'
+
+error_norm = 1
+c2 = light_speed * light_speed
+
+if args.dtype == 'dirac':
+    while error_norm > prec:
+        hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = defult_der)
+        v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
+        add_psi = hd_psi + v_psi
+        energy = spinor_H.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)
+        new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = defult_der)
+        new_orbital.crop(prec/10)
+        new_orbital.normalize()
+        delta_psi = new_orbital - spinor_H
+        deltasq = delta_psi.squaredNorm()
+        error_norm = np.sqrt(deltasq)
+        print('Error =', error_norm)
+        spinor_H = new_orbital
+
+    hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = defult_der)
+    v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
+    add_psi = hd_psi + v_psi
+    energy = spinor_H.dot(add_psi).real
+    print('Final Energy =',energy - light_speed**2)
+
+elif args.dtype == 'dirac2':
+    while error_norm > prec:
+        v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec) 
+        vv_psi = orb.apply_potential(-0.5/c2, V_tree, v_psi, prec)
+        beta_v_psi = v_psi.beta2()
+        apV_psi = v_psi.alpha_p(prec, defult_der)
+        ap_psi = spinor_H.alpha_p(prec, defult_der)
+        Vap_psi = orb.apply_potential(-1.0, V_tree, ap_psi, prec)
+        anticom = apV_psi + Vap_psi
+        RHS = beta_v_psi + vv_psi + anticom * (0.5/light_speed)
+        cke = spinor_H.classicT()
+        cpe = (spinor_H.dot(RHS)).real
+        print("Classic-like energies:", "cke =", cke,"cpe =", cpe,"cke + cpe =", cke + cpe)
+        mu = orb.calc_non_rel_mu(cke+cpe)
+        print("mu =", mu)
+        new_orbital = orb.apply_helmholtz(RHS, mu, prec)
+        new_orbital.normalize()
+        delta_psi = new_orbital - spinor_H
+        deltasq = delta_psi.squaredNorm()
+        error_norm = np.sqrt(deltasq)
+        print("Error =", error_norm)
+        spinor_H = new_orbital
+
+    hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der)
+    v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
+    add_psi = hd_psi + v_psi
+    energy = spinor_H.dot(add_psi).real
+
     cke = spinor_H.classicT()
-    cpe,imag = spinor_H.dot(RHS)
-    print('classic', cke,cpe,cpe+cke)
-    mu = orb.calc_non_rel_mu(cke+cpe)
-    print("mu", mu)
-    new_orbital = orb.apply_helmholtz(RHS, mu, prec)
-    new_orbital.normalize()
-    delta_psi = new_orbital - spinor_H
-    orbital_error, imag = delta_psi.dot(delta_psi)
-    print('Error',orbital_error, imag, flush = True)
-    spinor_H = new_orbital
-
-
-    #hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)
-    #v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
-    #add_psi = hd_psi + v_psi
-    #energy, imag = spinor_H.dot(add_psi)
-    #print('Energy',energy - light_speed**2,imag)
-    #mu = orb.calc_dirac_mu(energy, light_speed)
-    #tmp = orb.apply_helmholtz(v_psi, mu, prec)
-    #tmp.crop(prec/10)
-    #new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = default_der)
-    #new_orbital.crop(prec/10)
-    #new_orbital.normalize()
-    #delta_psi = new_orbital - spinor_H
-    #orbital_error, imag = delta_psi.dot(delta_psi)
-    #print('Error',orbital_error, imag)
-    #spinor_H = new_orbital
+    beta_v_psi = v_psi.beta2()
+    beta_pot = (beta_v_psi.dot(spinor_H)).real
+    pot_sq  = (v_psi.dot(v_psi)).real
+    ap_psi = spinor_H.alpha_p(prec, der)
+    anticom = (ap_psi.dot(v_psi)).real
+    energy_kutzelnigg = cke + beta_pot + pot_sq/(2*mc2) + anticom/light_speed
 
-#hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)
-#v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
-#add_psi = hd_psi + v_psi
-#energy, imag = spinor_H.dot(add_psi)
-#print('Final Energy',energy - light_speed**2)
-
-hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)
-v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
-add_psi = hd_psi + v_psi
-energy, imag = spinor_H.dot(add_psi)
-
-cke = spinor_H.classicT()
-beta_v_psi = v_psi.beta2()
-beta_pot,imag = beta_v_psi.dot(spinor_H)
-pot_sq, imag = v_psi.dot(v_psi)
-ap_psi = spinor_H.alpha_p(prec, 'BS')
-anticom, imag = ap_psi.dot(v_psi)
-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)
+    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, 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)

one_electron.py

@@ -9,121 +9,90 @@
 import numpy.linalg as LA
 import sys, getopt
 
+def analytic_1s(light_speed, n, k, Z):
+    alpha = 1/light_speed
+    gamma = orb.compute_gamma(k,Z,alpha)
+    tmp1 = n - np.abs(k) + gamma
+    tmp2 = Z * alpha / tmp1
+    tmp3 = 1 + tmp2**2
+    return light_speed**2 / np.sqrt(tmp3)
+
 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
-        print("Potential", potential)
-        Vpsi1 = orb.apply_potential(-1.0, potential, spinorb1, prec)
-        V1 = spinorb1.dot(Vpsi1)
-
-        hd_V_11 = hd_11 + V1
+    #compute_last_energy = False
 
-        eps = hd_V_11.real
-
-        print('orbital energy', eps - light_speed**2)
-
-        if(compute_last_energy):
-            break
+    light_speed = spinorb1.light_speed
 
-        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)
+    while error_norm > prec:
+        hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = 'ABGV')
+        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)
+        new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = 'ABGV')
+        new_orbital.crop(prec/10)
         new_orbital.normalize()
-        new_orbital.cropLargeSmall(prec)       
-
-        # Compute orbital error
         delta_psi = new_orbital - spinorb1
+        #orbital_error = delta_psi.dot(delta_psi).real
         deltasq = delta_psi.squaredNorm()
         error_norm = np.sqrt(deltasq)
-        print('Orbital_Error norm', error_norm)
+        print('Error', error_norm)
         spinorb1 = new_orbital
-        if(error_norm < prec):
-            compute_last_energy = True
+
+    hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der = 'ABGV')
+    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)
     return spinorb1
 
 
-def gs_D2_1e(spinorb1, potential, mra, prec, der):
+def gs_D2_1e(spinorb1, potential, mra, prec, der = 'ABGV'):
     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
+    mc2 = light_speed**2
 
-    while(error_norm > prec):
-        print("Applying V")
+    while error_norm > prec:
         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)
+        cke = spinorb1.classicT()
+        cpe = (spinorb1.dot(RHS)).real
+
         VVpsi.cropLargeSmall(prec)
-
-        print("anticom")
-        print(anticom)
-        print("beta_Vpsi")
-        print(beta_Vpsi)
-        print("VV_psi")
-        print(VVpsi)
-        RHS = beta_Vpsi + anticom + VVpsi
+        beta_Vpsi.cropLargeSmall(prec)
+        anticom.cropLargeSmall(prec)
         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)))
+        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)))
+
         mu = orb.calc_non_rel_mu(cke+cpe)
-        print("mu: ", mu)
+        print("mu =", mu)
 
-
         new_spinorb1 = orb.apply_helmholtz(RHS, mu, prec)
-        print("normalization")
+        #print("normalization")
         new_spinorb1.normalize()
-        print("crop")
+        #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)
+        print("Error =", error_norm)
         spinorb1 = new_spinorb1
 
     hd_psi = orb.apply_dirac_hamiltonian(spinorb1, prec, der)

orbital4c/nuclear_potential.py

@@ -52,9 +52,11 @@ def calculate_center_of_mass(coordinates):
         center_of_mass[i] /= total_mass
 
     return center_of_mass
+
 
-
-def pot(V_tree, coordinates, typenuc, mra, prec, der):
+def pot(coordinates, typenuc, mra, prec, der):
+    V_tree = vp.FunctionTree(mra)
+    V_tree.setZero()
     for atom, origin in coordinates.items():
         atom = get_original_list_name(atom)
         print("Atom:", atom)
@@ -88,6 +90,7 @@ def pot(V_tree, coordinates, typenuc, mra, prec, der):
             V = Peps(f)
         V_tree += V
     print('Define V Potential', typenuc, 'DONE')
+    return V_tree
 
 def point_charge(position, center , charge):
     d2 = ((position[0] - center[0])**2 +

orbital4c/orbital.py

@@ -448,5 +448,9 @@ def calc_kutzelnigg_mu(energy_sq, light_speed):
     return np.sqrt(-val)
 
 def calc_non_rel_mu(energy):
-    return np.sqrt(-2.0 * energy)
+    if energy > 0:
+        return np.sqrt(2.0 * energy)
+    elif energy < 0:
+        return np.sqrt(-2.0 * energy)
+