reading of atoms creating a center of mass and a spinorbital guess fr…

…om the combination of gaussians

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
+ H   0.1   0.2  -0.7
+ H   0.1   0.2   1.3

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

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 +

test.py

@@ -15,27 +15,18 @@
 import numpy.linalg as LA
 import sys, getopt
 
+import one_electron
 import two_electron
 
 import importlib
 importlib.reload(orb)
 
 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)

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)