From 5b99b4e44b6fb2aa223de96806d409c5cf2bb8eb Mon Sep 17 00:00:00 2001 From: Loris Burth Date: Wed, 19 Mar 2025 15:42:12 +0100 Subject: [PATCH] wip --- PyDuck.py | 116 +++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 85 insertions(+), 31 deletions(-) diff --git a/PyDuck.py b/PyDuck.py index 487bff5..f1c1f70 100755 --- a/PyDuck.py +++ b/PyDuck.py @@ -24,9 +24,9 @@ parser = argparse.ArgumentParser( # Initialize all the options for the script parser.add_argument('-b', '--basis', type=str, required=True, - help='Name of the file containing the basis set in the $QUACK_ROOT/basis/ directory. If cap is used the basis in psi4 style has to be provided in the directory cap_data/basis directory.') + help='Name of the file containing the basis set in the $QUACK_ROOT/basis/ directory. If cap is used the basis in psi4 style has to be provided in the directory data/basis_psi4 directory and same basis in nwchem format in data/basis_nwchem.') parser.add_argument('--use_local_basis', default=False, action='store_true', - help='If True, basis is loaded from local storage. Needed for CAP.') + help='If True, basis is loaded from local storage. Needed for CAP. From file in data/basis_nwchem/$basis') parser.add_argument('--bohr', default='Angstrom', action='store_const', const='Bohr', help='By default QuAcK assumes that the xyz files are in Angstrom. Add this argument if your xyz file is in Bohr.') parser.add_argument('-c', '--charge', type=int, default=0, @@ -50,14 +50,12 @@ parser.add_argument('--working_dir', type=str, default=QuAcK_dir, parser.add_argument('-x', '--xyz', type=str, required=True, help='Name of the file containing the nuclear coordinates in xyz format in the $QUACK_ROOT/mol/ directory without the .xyz extension') parser.add_argument("--use_cap", action="store_true", default=False, - help="If true cap integrals are calculated by opencap and written to a file. The basis has to be provided in the cap_data/basis dir in psi4 style and the onsets in the cap_data/onsets dir.") + help="If true cap integrals are calculated by opencap and written to a file. The basis has to be provided in the data/basis_psi4 dir in psi4 style and in nwchem format in the data/basis_nwchem format and the onsets in the data/onsets dir.") # Parse the arguments args = parser.parse_args() +working_dir = args.working_dir input_basis = args.basis -# Basisset path -pyscf_path = os.path.dirname(pyscf.__file__) -basis_path = os.path.join(pyscf_path, "gto", "basis", input_basis + ".dat") unit = args.bohr charge = args.charge frozen_core = args.frozen_core @@ -68,7 +66,6 @@ print_2e = args.print_2e formatted_2e = args.formatted_2e mmap_2e = args.mmap_2e aosym_2e = args.aosym_2e -working_dir = args.working_dir # Read molecule f = open(working_dir+'/mol/'+xyz, 'r') lines = f.read().splitlines() @@ -81,16 +78,32 @@ for line in lines: pos = (float(tmp[1]), float(tmp[2]), float(tmp[3])) list_pos_atom.append([atom, pos]) f.close() +# Create PySCF molecule +if args.use_local_basis: + print("blub") + atoms = list(set(atom[0] for atom in list_pos_atom)) + basis_dict = {atom: gto.basis.parse_nwchem.load( + working_dir + "/data/basis_nwchem/" + input_basis, atom) for atom in atoms} + basis = basis_dict + mol = gto.M( + atom=list_pos_atom, + basis=basis, + charge=charge, + spin=multiplicity - 1 + # symmetry = True # Enable symmetry + ) +else: + mol = gto.M( + atom=list_pos_atom, + basis=input_basis, + charge=charge, + spin=multiplicity - 1 + # symmetry = True # Enable symmetry + ) -# Definition of the molecule -mol = gto.M( - atom=list_pos_atom, - basis=input_basis, - charge=charge, - spin=multiplicity - 1 - # symmetry = True # Enable symmetry -) +# Example usage: +print(basis_dict) # Fix the unit for the lengths mol.unit = unit # @@ -135,35 +148,79 @@ subprocess.call(['rm', '-f', working_dir + '/int/nBas.dat']) f = open(working_dir+'/int/nBas.dat', 'w') f.write(" {} ".format(str(norb))) f.close() +print("before cap") + + +def create_psi4_basis(basis_dict): + """ + Converts a dictionary representation of a basis set (pyscf internal format) into a Psi4-formatted string. + + Parameters: + basis_dict (dict): A dictionary where keys are element symbols and values are lists of primitives. + Each primitive is represented as [angular momentum, [exponent, coefficient(s)]]. + + Returns: + str: The Psi4-formatted basis set string. + """ + l_mapping = {0: 'S', 1: 'P', 2: 'D', 3: 'F', 4: 'G', 5: 'H'} + basis_str = "****\n" + for element, shells in basis_dict.items(): + basis_str += f"{element} 0\n" + + for shell in shells: + l_value = shell[0] + l_letter = l_mapping.get(l_value, str(l_value)) + primitives = shell[1:] + num_primitives = len(primitives) + + # Determine number of contractions + max_contractions = max(len(p) - 1 for p in primitives) + + for contraction_idx in range(max_contractions): + basis_str += f"{l_letter} {num_primitives} 1.00\n" + for primitive in primitives: + exponent = primitive[0] + coefficient = primitive[1 + contraction_idx] if len( + primitive) > (1 + contraction_idx) else 0.0 + basis_str += f" {exponent: .6E} {coefficient: .6E}\n" + basis_str += "****\n" + basis_filename_psi4 = working_dir + "data/basis_psi4/" + input_basis + with open(basis_filename_psi4, "w") as file: + file.write(basis_str.strip()) + return basis_filename_psi4 + # CAP definition if args.use_cap: - f = open(working_dir+'/cap_data/onsets/'+args.xyz, 'r') + f = open(working_dir+'/data/onsets/'+args.xyz, 'r') lines = f.read().splitlines() for line in lines: tmp = line.split() onset_x = float(tmp[0]) onset_y = float(tmp[1]) onset_z = float(tmp[2]) + f.close() + print("onsets read") # xyz file with open(working_dir + "/mol/" + xyz, "r") as f: lines = f.readlines() + f.close() + print("xyz read") num_atoms = int(lines[0].strip()) atoms = [line.strip() for line in lines[2:2+num_atoms]] sys_dict = { "molecule": "inline", "geometry": "\n".join(atoms), # XYZ format as a string - "basis_file": working_dir + "/cap_data/basis/" + input_basis, + "basis_file": create_psi4_basis(basis_dict), "bohr_coordinates": unit == 'Bohr' } + print("sys dict constructed") cap_system = pyopencap.System(sys_dict) - print("opencap") - print(cap_system.get_overlap_mat("pyscf")) - print("pyscf") - print(ovlp) - # if not(cap_system.check_overlap_mat(ovlp, "pyscf")): - # raise Exception( - # "Provided cap basis does not match to the pyscf basis.") + print("cap system created") + print("Before overlap check") + if not(cap_system.check_overlap_mat(ovlp, "pyscf")): + raise Exception( + "Provided cap basis does not match to the pyscf basis.") cap_dict = {"cap_type": "box", "cap_x": onset_x, "cap_y": onset_y, @@ -171,13 +228,9 @@ if args.use_cap: "Radial_precision": "16", "angular_points": "590", "thresh": 10} + print("Before cap system") pc = pyopencap.CAP(cap_system, cap_dict, norb) - pc.get_ao_cap(ordering="pyscf") - pc.renormalize_cap(ovlp, "pyscf") cap_ao = pc.get_ao_cap(ordering="pyscf") - print(cap_ao) - print("S after normalizing") - print(cap_system.get_overlap_mat("pyscf")) def write_matrix_to_file(matrix, size, file, cutoff=1e-15): @@ -205,8 +258,9 @@ subprocess.call(['rm', '-f', working_dir + '/int/y.dat']) write_matrix_to_file(y, norb, working_dir+'/int/y.dat') subprocess.call(['rm', '-f', working_dir + '/int/z.dat']) write_matrix_to_file(z, norb, working_dir+'/int/z.dat') -subprocess.call(['rm', '-f', working_dir + '/int/CAP.dat']) -write_matrix_to_file(cap_ao, norb, working_dir+'/int/CAP.dat') +if args.use_cap: + subprocess.call(['rm', '-f', working_dir + '/int/CAP.dat']) + write_matrix_to_file(cap_ao, norb, working_dir+'/int/CAP.dat') def write_tensor_to_file(tensor, size, file_name, cutoff=1e-15):