mirror of
https://github.com/pfloos/quack
synced 2025-05-06 15:14:55 +02:00
wip
This commit is contained in:
Loris Burth
parent
c80aabeab6
commit
5b99b4e44b
116
PyDuck.py
116
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):
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user