Merge pull request #49 from ilfreddy/JKcoulomb

Reactivating operators + D2 for 2 electrons

ReVAMPyR.py

@@ -108,8 +108,7 @@
 if args.readOrbitals == 'Yes':
     spinorb1 = orb.orbital4c()
     spinorb1.read(args.atype)
-    spinorb2 = spinorb1.ktrs()
-    spinorb2.cropLargeSmall(prec)
+    spinorb2 = spinorb1.ktrs(prec)
     spinorb2.normalize()
     print('Read alpha spinorbital and defined beta spinorbital using KTRS DONE')
 elif args.readOrbitals == 'No':
@@ -128,8 +127,7 @@
     spinorb1.copy_components(La=complexfc)
     spinorb1.init_small_components(prec/10, der)
     spinorb1.normalize()
-    spinorb2 = spinorb1.ktrs()
-    spinorb2.cropLargeSmall(prec)
+    spinorb2 = spinorb1.ktrs(prec)
     spinorb2.normalize()
     print('Define alpha and beta spinorbitals DONE')
 

orbital4c/complex_fcn.py

@@ -65,6 +65,18 @@ def __sub__(self, other):
         output.imag = self.imag - other.imag
         return output
 
+    def real_mul(self, coeff):
+        output = complex_fcn()
+        output.real = self.real * coeff
+        output.imag = self.imag * coeff
+        return output
+
+    def imag_mul(self, coeff):
+        output = complex_fcn()
+        output.real = self.imag * (-1.0 * coeff)
+        output.imag = self.real * coeff
+        return output
+
     def __rmul__(self, other):
         output = complex_fcn()
         output.real = self.real * np.real(other) - self.imag * np.imag(other)
@@ -247,9 +259,9 @@ def apply_potential(factor, potential, func, prec):
 def apply_helmholtz(func, mu, light_speed, prec):
     out_func = complex_fcn()
     H = vp.HelmholtzOperator(func.mra, mu, prec)
-    if(func.real.squaredNorm() > 0):
+    if(func.real.squaredNorm() > 1e-12):
         vp.advanced.apply(prec, out_func.real, H, func.real)
-    if(func.imag.squaredNorm() > 0):
+    if(func.imag.squaredNorm() > 1e-12):
         vp.advanced.apply(prec, out_func.imag, H, func.imag)
     return out_func
 
@@ -320,3 +332,17 @@ def vector_gradient(vector, der = "BS"):
         tensor.append(vector[i].gradient(der))
     return tensor
 
+# Note: some thresholding of the contributions should be considered here.
+def add_vector(func_array, coeff_array, prec):
+    output = complex_fcn()
+    real_array = []
+    imag_array = []
+    for i in range(len(func_array)):
+        real_array.append(( coeff_array[i].real, func_array[i].real))
+        real_array.append((-coeff_array[i].imag, func_array[i].imag))
+        imag_array.append(( coeff_array[i].real, func_array[i].imag))
+        imag_array.append(( coeff_array[i].imag, func_array[i].real))
+    vp.advanced.add(prec, output.real, real_array)
+    vp.advanced.add(prec, output.imag, imag_array)
+    return output
+

orbital4c/operators.py

@@ -16,7 +16,7 @@ def __init__(self, mra, guessorb, c, origin, prec):
         self.complexfc = None
 
         if self.guessorb == 'slater':
-            print('cazzo')
+            print("Slater-type not yet implemented")
         elif guessorb == 'gaussian':
 ################################   DEFINE GAUSSIAN FUNCTION AS GUESS  ################################
             a_coeff = 3.0
@@ -46,42 +46,109 @@ def __call__(self, component):
         phi.normalize()
         return phi
 
-class CoulombDirectOperator():
+class Operator():
     def __init__(self, mra, prec, Psi):
         self.mra = mra
         self.prec = prec
         self.Psi = Psi
+
+    def matrix_true(self, Phi):
+        n_orbitals = len(Phi)
+        mat = np.zeros((n_orbitals, n_orbitals), complex)
+        for i in range(n_orbitals):
+            si = Phi[i]
+            Osi = self(si)
+            for j in range(i+1):
+                sj = Phi[j]
+                val = sj.dot(Osi)
+                mat[j][i] = val
+                if (i != j):
+                    mat[i][j] = np.conjugate(mat[j][i]) 
+        return mat
+
+    def matrix(self, Phi):
+        n_orbitals = len(Phi)
+        mat = np.zeros((n_orbitals, n_orbitals), complex)
+        si = Phi[0]
+        Osi = self(si)
+        sj = Phi[0]
+        val = sj.dot(Osi)
+        mat[0][0] = val
+        mat[1][1] = val
+        return mat
+
+
+class CoulombDirectOperator(Operator):
+    def __init__(self, mra, prec, Psi):
+        super().__init__(mra, prec, Psi)
         self.poisson = vp.PoissonOperator(mra=self.mra, prec=self.prec)
         self.potential = None
         self.setup()
 
     def setup(self):
-        rho = self.Psi[0].density(self.prec)
-        for i in range(1, len(self.Psi)):
-            rho += self.Psi[i].density(self.prec)
-        rho.crop(self.prec)
-        self.potential = (4.0*np.pi)*self.poisson(rho).crop(self.prec)
-
-    def __call__(self, Phi):
-        return self.potential * Phi
+        rho = vp.FunctionTree(self.mra)
+        rholist = []
+        for i in range(0, len(self.Psi)):
+            dens = self.Psi[i].overlap_density(self.Psi[i], self.prec)
+            rholist.append((1.0, dens.real))
+        vp.advanced.add(self.prec, rho, rholist)
+        self.potential = (4.0*np.pi) * self.poisson(rho).crop(self.prec)
 
+    def __call__(self, phi):
+        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.cropLargeSmall(self.prec)
+        return out
 
-class CoulombExchangeOperator():
+class CoulombExchangeOperator(Operator):
     def __init__(self, mra, prec, Psi):
-        self.mra = mra
-        self.prec = prec
-        self.Psi = Psi
-        self.poisson = vp.PoissonOperator(mra=mra, prec=self.prec)
+        super().__init__(mra, prec, Psi)
+        self.poisson = vp.PoissonOperator(mra=self.mra, prec=self.prec)
         self.potential = None
 
+    def __call__(self, phi):
+        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)
+            Kij_array.append(tmp)
+            coeff_array.append(1.0)
+        output = orb.add_vector(Kij_array, coeff_array, self.prec) 
+        output *= 4.0 * np.pi
+        output.cropLargeSmall(self.prec)
+        return output
 
-    def __call__(self, Phi):
-        V_j0 = self.poisson(self.Psi[0].exchange(Phi, self.prec))
-        self.potential = (V_j0 * self.Psi[0])
-        for i in range(1, len(self.Psi)):
-            V_ji = self.poisson(self.Psi[i].exchange(Phi, self.prec))
-            self.potential += (V_ji * self.Psi[i])
-        self.potential *= 4.0*np.pi
-        return self.potential
+class PotentialOperator(Operator):
+    def __init__(self, mra, prec, potential, real = True):
+        super().__init__(mra, prec, [])
+        self.potential = potential
+        self.real = real
+
+    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.cropLargeSmall(self.prec)
+        return result
+
+class FockOperator(Operator):
+    def __init__(self, mra, prec, operators, factors, der = "ABGV"):
+        super().__init__(mra, prec, [])
+        self.der = der
+        self.operators = operators
+        self.factors = factors
 
+    def __call__(self, phi):
+        Fphi = orb.apply_dirac_hamiltonian(phi, self.prec, shift = 0.0, der = self.der)
+        for i in range(len(self.operators)):
+            Fphi += self.factors[i] * self.operators[i](phi)
+        Fphi.cropLargeSmall(self.prec)
+        return Fphi
 

orbital4c/orbital.py

@@ -101,8 +101,8 @@ def cropLargeSmall(self, prec):
         largeNorm = np.sqrt(self.squaredLargeNorm())
         smallNorm = np.sqrt(self.squaredSmallNorm())
         precLarge = prec * largeNorm
-        precSmall = prec * smallNorm
-        print('precisions', precLarge, precSmall)
+        precSmall = prec * largeNorm
+#        print('precisions', precLarge, precSmall)
         self['La'].crop(precLarge, True)
         self['Lb'].crop(precLarge, True)
         self['Sa'].crop(precSmall, True)
@@ -268,7 +268,9 @@ def alpha_p(self, prec, der = "ABGV"):
         apx = orb_grad[0].alpha(0, prec)
         apy = orb_grad[1].alpha(1, prec)
         apz = orb_grad[2].alpha(2, prec)
-        return -1j * (apx + apy + apz)
+        result = -1j * (apx + apy + apz)
+        result.cropLargeSmall(prec)
+        return result
 
     def classicT(self, der = 'ABGV'):
         orb_grad = self.gradient(der)
@@ -280,15 +282,15 @@ def classicT(self, der = 'ABGV'):
     def alpha_vector(self, prec):
         return [self.alpha(0, prec), self.alpha(1, prec), self.alpha(2, prec)]
 
-    def ktrs(self):   #Kramers´ Time Reversal Symmetry
+    def ktrs(self, prec):   #Kramers´ Time Reversal Symmetry
         out_orb = orbital4c()
         tmp = self.complex_conj()
         ktrs_order = np.array([1, 0, 3, 2])
         ktrs_coeff = np.array([-1,  1,  -1,  1])
         for idx in range(4):
             coeff = ktrs_coeff[idx]
             comp = ktrs_order[idx]
-            out_orb.comp_array[idx] = coeff * tmp.comp_array[comp]
+            out_orb.comp_array[idx] = tmp.comp_array[comp].real_mul(coeff)
         return out_orb
 
 #Beta c**2
@@ -368,6 +370,16 @@ def apply_complex_potential(factor, potential, orbital, prec):
 #            print('Warning: adding two empty trees')
 #    return out_orb
 
+def add_vector(orbital_array, coeff_array, prec):
+    output = orbital4c()
+    for comp in output.comp_dict:
+        func_array = []
+        for orbital in orbital_array:
+            func_array.append(orbital[comp])
+        output[comp] = cf.add_vector(func_array, coeff_array, prec)
+    return output
+
+
 def apply_helmholtz(orbital, mu, prec):
     out_orbital = orbital4c()
     for comp in orbital.comp_dict.keys():
@@ -430,6 +442,11 @@ def alpha_gradient(orbital, prec):
 def calc_dirac_mu(energy, light_speed):
     return np.sqrt((light_speed**4-energy**2)/light_speed**2)
 
+def calc_kutzelnigg_mu(energy_sq, light_speed):
+    c2 = light_speed**2
+    val = energy_sq/c2 - c2
+    return np.sqrt(-val)
+
 def calc_non_rel_mu(energy):
     return np.sqrt(-2.0 * energy)
 

test.py

@@ -69,7 +69,7 @@
     Z = args.charge
     atom = args.atype
     ################# Call MRA #######################
-    mra = vp.MultiResolutionAnalysis(box=[-args.box,args.box], order=args.order)
+    mra = vp.MultiResolutionAnalysis(box=[-args.box,args.box], order=args.order, max_depth = 30)
     prec = args.prec
     origin = [args.cx, args.cy, args.cz]
 
@@ -78,18 +78,20 @@
     cf.complex_fcn.mra = mra
     print('call MRA DONE')
 
-    computeNuclearPotential = False
-    readOrbitals = True
-    runCoulomb = False
-    saveOrbitals = False
-    runGaunt = True
-    runGaugeA = True
-    runGaugeB = True
-    runGaugeC = True
-    runGaugeD = True
-    runGaugeDelta = True
+    computeNuclearPotential = True
+    readOrbitals            = True
+    runCoulomb              = False
+    runCoulombGen           = False
+    runKutzelnigg           = False
+    runKutzSimple           = True
+    saveOrbitals            = False
+    runGaunt                = False 
+    runGaugeA               = False 
+    runGaugeB               = False 
+    runGaugeC               = False 
+    runGaugeD               = False 
+    runGaugeDelta           = False 
     default_der = args.deriv
-
     ################### Define V potential ######################
     if(computeNuclearPotential):
         if args.potential == 'point_charge':
@@ -130,13 +132,23 @@
         spinorb1.copy_components(La=complexfc)
         spinorb1.init_small_components(prec/10)
         spinorb1.normalize()
-        spinorb2 = spinorb1.ktrs() #does this go out of scope?
+        spinorb1.cropLargeSmall(prec)
+        spinorb2 = spinorb1.ktrs(prec) #does this go out of scope?
 
     length = 2 * args.box
 
-    if runCoulomb:
-        spinorb1, spinorb2 = two_electron.coulomb_gs_2e(spinorb, V_tree, mra, prec)
+    if runCoulombGen:
+        spinorb1, spinorb2 = two_electron.coulomb_gs_gen([spinorb1, spinorb2], 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)