4
1
mirror of https://github.com/pfloos/quack synced 2025-01-03 01:56:09 +01:00
This commit is contained in:
Abdallah Ammar 2024-08-29 20:48:16 +02:00
parent 8d7fb2a292
commit 9aebbe74f0
13 changed files with 410 additions and 307 deletions

3
.gitignore vendored
View File

@ -1,2 +1,5 @@
*.o
*.
__pycache__
.ninja_deps

6
tests/.gitignore vendored Normal file
View File

@ -0,0 +1,6 @@
FeatherBench.db
FeatherBench.json
*.xyz

View File

@ -4,21 +4,21 @@ import sqlite3
from molecule import save_molecules_to_json, load_molecules_from_json
from molecule import create_database, add_molecule_to_db
from swift_set import swiftset
from feather_bench import FeatherBench
# Save molecules to JSON
save_molecules_to_json(swiftset, 'swiftset.json')
save_molecules_to_json(FeatherBench, 'FeatherBench.json')
# Load molecules from JSON
loaded_molecules = load_molecules_from_json('swiftset.json')
loaded_molecules = load_molecules_from_json('FeatherBench.json')
print(loaded_molecules)
# Create a database and add molecules
db_name = 'swiftset.db'
db_name = 'FeatherBench.db'
create_database(db_name)
for molecule in swiftset:
for molecule in FeatherBench:
add_molecule_to_db(db_name, molecule)

93
tests/feather_bench.py Normal file
View File

@ -0,0 +1,93 @@
from molecule import Molecule
He = Molecule(
name="He",
multiplicity=1,
geometry=[
{"element": "He", "x": 0.0, "y": 0.0, "z": 0.0}
],
properties={
"properties_rhf":{
"6-31g": {
"RHF energy": -2.855160426884076,
"RHF HOMO energy": -0.914126628614305,
"RHF LUMO energy": 1.399859335225087,
"RHF dipole moment": 0.000000000000000,
"RMP2 correlation energy": -0.011200122910187,
"CCD correlation energy": -0.014985063408247,
"DCD correlation energy": -0.014985062907429,
"CCSD correlation energy": -0.015001711549550,
"drCCD correlation energy": -0.018845374502248,
"rCCD correlation energy": -0.016836324636164,
"crCCD correlation energy": 0.008524677369855,
"lCCD correlation energy": -0.008082420815100,
"pCCD correlation energy": -0.014985062519068,
"RCIS singlet excitation energy": 1.911193619935257,
"RCIS triplet excitation energy": 1.455852629402236,
"phRRPA correlation energy": -0.018845374129105,
"phRRPAx correlation energy": -0.015760565121283,
"crRRPA correlation energy": -0.008868581132405,
"ppRRPA correlation energy": -0.008082420815100,
"RG0F2 correlation energy": -0.011438430540374,
"RG0F2 HOMO energy": -0.882696116247871,
"RG0F2 LUMO energy": 1.383080391811630,
"evRGF2 correlation energy": -0.011448483158486,
"evRGF2 HOMO energy": -0.881327878713477,
"evRGF2 LUMO energy": 1.382458968133448,
"RG0W0 correlation energy": -0.019314094399756,
"RG0W0 HOMO energy": -0.870533880190454,
"RG0W0 LUMO energy": 1.377171287010956,
"evRGW correlation energy": -0.019335511771724,
"evRGW HOMO energy": -0.868460640957913,
"evRGW LUMO energy": 1.376287581471769,
"RG0T0pp correlation energy": -0.008082420815100,
"RG0T0pp HOMO energy": -0.914126628614305,
"RG0T0pp LUMO energy": 1.399859335225087,
"evRGTpp correlation energy": -0.008082420815100,
"evRGTpp HOMO energy": -0.914126628614305,
"evRGTpp LUMO energy": 1.399859335225087
}
},
"properties_uhf":{
"6-31g": {
}
},
"properties_ghf":{
"6-31g": {
}
},
"properties_rohf":{
"6-31g": {
}
}
}
)
# ---
#H2O = Molecule(
# name="H2O",
# multiplicity=1,
# geometry=[
# {"element": "O", "x": 0.0000, "y": 0.0000, "z": 0.0000},
# {"element": "H", "x": 0.7571, "y": 0.0000, "z": 0.5861},
# {"element": "H", "x": -0.7571, "y": 0.0000, "z": 0.5861}
# ],
# properties={
# "cc-pvdz": {
# }
#)
# ---
FeatherBench = [
He,
#H2O
]

18
tests/inp/options.RHF Normal file
View File

@ -0,0 +1,18 @@
# HF: maxSCF thresh DIIS guess mix shift stab search
10000 0.0000001 5 1 0.0 0.0 F F
# MP: reg
F
# CC: maxSCF thresh DIIS
64 0.0000001 5
# spin: TDA singlet triplet
F T T
# GF: maxSCF thresh DIIS lin eta renorm reg
256 0.00001 5 F 0.0 0 F
# GW: maxSCF thresh DIIS lin eta TDA_W reg
256 0.00001 5 F 0.0 F F
# GT: maxSCF thresh DIIS lin eta TDA_T reg
256 0.00001 5 F 0.0 F F
# ACFDT: AC Kx XBS
F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA
F F F F T

205
tests/lunch_bench.py Normal file
View File

@ -0,0 +1,205 @@
import sys
import os
import shutil
from pathlib import Path
import subprocess
import platform
from datetime import datetime
import argparse
from molecule import get_molecules_from_db
from molecule import generate_xyz
from utils import print_col
current_date = datetime.now()
quack_root = os.getenv('QUACK_ROOT')
# User Name
user_name = os.getlogin()
# Operating System
os_name = platform.system()
os_release = platform.release()
os_version = platform.version()
# CPU Information
machine = platform.machine()
processor = platform.processor()
# System Architecture
architecture = platform.architecture()[0]
# Python Version
python_version_full = platform.python_version_tuple()
PYTHON_VERSION = "{}.{}".format(python_version_full[0], python_version_full[1])
print(f"The current date and time is {current_date.strftime('%Y-%m-%d %H:%M:%S')}")
print(f"User Name: {user_name}")
print(f"Operating System: {os_name} {os_release} ({os_version})")
print(f"CPU: {processor} ({machine})")
print(f"System Architecture: {architecture}")
print(f"QUACK_ROOT: {quack_root}")
print(f"Python version: {python_version_full}\n\n")
parser = argparse.ArgumentParser(description="Benchmark Data Sets")
parser.add_argument(
'-s', '--set_type',
choices=['light', 'medium', 'heavy'],
default='light',
help="Specify the type of data set: light (default), medium, or heavy."
)
args = parser.parse_args()
if args.set_type == 'light':
bench = 'FeatherBench'
bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
elif args.set_type == 'medium':
bench = 'BalanceBench'
bench_title = "\n\nSelected Medium Benchmark: {}\n\n".format(bench)
elif args.set_type == 'heavy':
bench = 'TitanBench'
bench_title = "\n\nSelected Heavy Benchmark: {}\n\n".format(bench)
else:
bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
print(bench_title.center(150, '-'))
print("\n\n")
# ---
class Quack_Job:
def __init__(self, mol, multip, basis, geom, methd):
self.mol = mol
self.multip = multip
self.basis = basis
self.geom = geom
self.methd = methd
def prep_inp(self):
# geometry
generate_xyz(self.geom, filename="{}.xyz".format(self.mol))
# input files
for inp in ["methods", "options"]:
inp_file = "{}.{}".format(inp, self.methd.upper())
if os.path.exists("inp/{}".format(inp_file)):
shutil.copy("inp/{}".format(inp_file), "../mol/{}".format(inp_file))
else:
print_col("File 'inp/{}' does not exist.".format(inp_file), "red")
sys.exit(1)
def run(file_out, mol, bas, multip):
os.chdir('..')
print(f" :$ cd ..")
for file_in in ["methods", "options"]:
command = ['cp', 'tests/{}.RHF'.format(file_in), 'input/{}'.format(file_in)]
print(f" :$ {' '.join(command)}")
result = subprocess.run(command, capture_output=True, text=True)
if result.returncode != 0:
print("Error moving file: {}".format(result.stderr))
command = [
'python{}'.format(PYTHON_VERSION), 'PyDuck.py',
'-x', '{}'.format(mol),
'-b', '{}'.format(bas),
'-m', '{}'.format(multip)
]
print(f" :$ {' '.join(command)}")
with open(file_out, 'w') as fobj:
result = subprocess.run(command, stdout=fobj, stderr=subprocess.PIPE, text=True)
if result.stderr:
print("Error output:", result.stderr)
os.chdir('tests')
print(f" :$ cd tests")
# ---
def main():
work_path = Path('{}/tests/work'.format(quack_root))
if not work_path.exists():
work_path.mkdir(parents=True, exist_ok=True)
print(f"Directory '{work_path}' created.\n")
for mol in molecules:
mol_name = mol.name
mol_mult = mol.multiplicity
mol_geom = mol.geometry
mol_data = mol.properties
print_col(" Molecule: {} (2S+1 = {})".format(mol_name, mol_mult), "blue")
for mol_prop_name, mol_prop_data in mol_data.items():
print_col(" Testing {}".format(mol_prop_name), "cyan")
methd = mol_prop_name[len('properties_'):]
if(len(mol_prop_data) == 0):
print_col(" {} is empty. Skipping...".format(mol_prop_name), "cyan")
print()
continue
for basis_name, basis_data in mol_prop_data.items():
print_col(" Basis set = {}".format(basis_name), "yellow")
if(len(basis_data) == 0):
print_col(" {} is empty. Skipping...".format(basis_name), "yellow")
print()
continue
work_methd = Path('{}/{}'.format(work_path, methd))
if not work_methd.exists():
work_methd.mkdir(parents=True, exist_ok=True)
#print(f"Directory '{work_methd}' created.\n")
New_Quack_Job = Quack_Job(mol_name, mol_mult, basis_name, mol_geom, methd)
New_Quack_Job.prep_inp()
# for name, val in basis_data.items():
# print(f" name = {name}")
# print(f" val = {val}")
print()
print()
print()
quit()
# # create input files
# class_methd.gen_input()
#
# file_out = "{}/{}/{}_{}_{}.out".format(work_path, prop, mol_name, mol_mult, bas)
#
# print(" testing {} for {}@{} (2S+1 = {})".format(prop, mol_name, bas, mol_mult))
# print(" file_out: {}".format(file_out))
#
# class_methd.run_job(file_out, mol_name, bas, mol_mult)
db_name = '{}.db'.format(bench)
molecules = get_molecules_from_db(db_name)
main()

View File

@ -3,22 +3,18 @@ import json
import sqlite3
class Molecule:
def __init__(self, name, multiplicity, geometry, energies):
def __init__(self, name, multiplicity, geometry, properties):
self.name = name
self.multiplicity = multiplicity
self.geometry = geometry # List of tuples (atom, x, y, z)
self.energies = energies # Dictionary of dictionaries: {method: {basis: energy}}
def get_energy(self, method, basis):
"""Retrieve energy for a specific method and basis set."""
return self.energies.get(method, {}).get(basis, None)
self.geometry = geometry
self.properties = properties
def to_dict(self):
return {
"name": self.name,
"multiplicity": self.multiplicity,
"geometry": self.geometry,
"energies": self.energies,
"properties": self.properties,
}
@staticmethod
@ -27,7 +23,7 @@ class Molecule:
name=data["name"],
multiplicity=data["multiplicity"],
geometry=data["geometry"],
energies=data["energies"]
properties=data["properties"]
)
def save_molecules_to_json(molecules, filename):
@ -45,7 +41,7 @@ def create_database(db_name):
conn = sqlite3.connect(db_name)
cursor = conn.cursor()
cursor.execute('''CREATE TABLE IF NOT EXISTS molecules
(name TEXT, multiplicity INTEGER, geometry TEXT, energies TEXT)''')
(name TEXT, multiplicity INTEGER, geometry TEXT, properties TEXT)''')
conn.commit()
conn.close()
@ -53,7 +49,7 @@ def add_molecule_to_db(db_name, molecule):
conn = sqlite3.connect(db_name)
cursor = conn.cursor()
geometry_str = json.dumps(molecule.geometry)
energies_str = json.dumps(molecule.energies)
energies_str = json.dumps(molecule.properties)
cursor.execute("INSERT INTO molecules VALUES (?, ?, ?, ?)",
(molecule.name, molecule.multiplicity, geometry_str, energies_str))
conn.commit()
@ -62,25 +58,48 @@ def add_molecule_to_db(db_name, molecule):
def get_molecules_from_db(db_name):
conn = sqlite3.connect(db_name)
cursor = conn.cursor()
cursor.execute("SELECT name, multiplicity, geometry, energies FROM molecules")
cursor.execute("SELECT name, multiplicity, geometry, properties FROM molecules")
rows = cursor.fetchall()
molecules = []
for row in rows:
name, multiplicity, geometry_str, energies_str = row
geometry = json.loads(geometry_str)
energies = json.loads(energies_str) # energies is a dictionary of dictionaries
molecules.append(Molecule(name, multiplicity, geometry, energies))
properties = json.loads(energies_str)
molecules.append(Molecule(name, multiplicity, geometry, properties))
conn.close()
return molecules
def write_geometry_to_xyz(molecule, filename):
def generate_xyz(elements, filename="output.xyz", verbose=False):
"""
Generate an XYZ file from a list of elements.
Parameters:
elements (list): A list of dictionaries, where each dictionary represents
an atom with its element and x, y, z coordinates.
filename (str): The name of the output XYZ file. Default is 'output.xyz'.
"""
# Get the number of atoms
num_atoms = len(elements)
# Open the file in write mode
with open(filename, 'w') as f:
# First line: number of atoms
f.write(f"{len(molecule.geometry)}\n")
# Second line: empty comment line
f.write("\n")
# Remaining lines: atom positions
for atom, x, y, z in molecule.geometry:
f.write(f"{atom} {x:.6f} {y:.6f} {z:.6f}\n")
# Write the number of atoms
f.write(f"{num_atoms}\n")
# Write a comment line (can be left blank or customized)
f.write("XYZ file generated by generate_xyz function\n")
# Write the element and coordinates
for atom in elements:
element = atom['element']
x = atom['x']
y = atom['y']
z = atom['z']
f.write(f"{element} {x:.6f} {y:.6f} {z:.6f}\n")
if(verbose):
print(f"XYZ file '{filename}' generated successfully!")

View File

@ -1,76 +0,0 @@
from molecule import Molecule
He = Molecule(
name="He",
multiplicity=1,
geometry=[
{"element": "He", "x": 0.0, "y": 0.0, "z": 0.0}
],
properties={
"6-31g": {
"RHF energy": -2.855160426884076,
"RHF HOMO energy": -0.914126628614305,
"RHF LUMO energy": 1.399859335225087,
"RHF dipole moment": 0.000000000000000,
"RMP2 correlation energy": -0.011200122910187,
"CCD correlation energy": -0.014985063408247,
"DCD correlation energy": -0.014985062907429,
"CCSD correlation energy": -0.015001711549550,
"drCCD correlation energy": -0.018845374502248,
"rCCD correlation energy": -0.016836324636164,
"crCCD correlation energy": 0.008524677369855,
"lCCD correlation energy": -0.008082420815100,
"pCCD correlation energy": -0.014985062519068,
"RCIS singlet excitation energy": 1.911193619935257,
"RCIS triplet excitation energy": 1.455852629402236,
"phRRPA correlation energy": -0.018845374129105,
"phRRPAx correlation energy": -0.015760565121283,
"crRRPA correlation energy": -0.008868581132405,
"ppRRPA correlation energy": -0.008082420815100,
"RG0F2 correlation energy": -0.011438430540374,
"RG0F2 HOMO energy": -0.882696116247871,
"RG0F2 LUMO energy": 1.383080391811630,
"evRGF2 correlation energy": -0.011448483158486,
"evRGF2 HOMO energy": -0.881327878713477,
"evRGF2 LUMO energy": 1.382458968133448,
"RG0W0 correlation energy": -0.019314094399756,
"RG0W0 HOMO energy": -0.870533880190454,
"RG0W0 LUMO energy": 1.377171287010956,
"evRGW correlation energy": -0.019335511771724,
"evRGW HOMO energy": -0.868460640957913,
"evRGW LUMO energy": 1.376287581471769,
"RG0T0pp correlation energy": -0.008082420815100,
"RG0T0pp HOMO energy": -0.914126628614305,
"RG0T0pp LUMO energy": 1.399859335225087,
"evRGTpp correlation energy": -0.008082420815100,
"evRGTpp HOMO energy": -0.914126628614305,
"evRGTpp LUMO energy": 1.399859335225087
}
}
)
# ---
H2O = Molecule(
name="H2O",
multiplicity=1,
geometry=[
{"element": "O", "x": 0.0000, "y": 0.0000, "z": 0.0000},
{"element": "H", "x": 0.7571, "y": 0.0000, "z": 0.5861},
{"element": "H", "x": -0.7571, "y": 0.0000, "z": 0.5861}
],
properties={
"cc-pvdz": {
}
)
# ---
SwiftSet = [He, H2O]

View File

@ -1,204 +0,0 @@
import os
from pathlib import Path
import subprocess
import platform
from datetime import datetime
from molecule import get_molecules_from_db
current_date = datetime.now()
quack_root = os.getenv('QUACK_ROOT')
# User Name
user_name = os.getlogin()
# Operating System
os_name = platform.system()
os_release = platform.release()
os_version = platform.version()
# CPU Information
machine = platform.machine()
processor = platform.processor()
# System Architecture
architecture = platform.architecture()[0]
# Python Version
python_version_full = platform.python_version_tuple()
PYTHON_VERSION = "{}.{}".format(python_version_full[0], python_version_full[1])
print(f"The current date and time is {current_date.strftime('%Y-%m-%d %H:%M:%S')}")
print(f"User Name: {user_name}")
print(f"Operating System: {os_name} {os_release} ({os_version})")
print(f"CPU: {processor} ({machine})")
print(f"System Architecture: {architecture}")
print(f"QUACK_ROOT: {quack_root}")
print(f"Python version: {python_version_full}\n\n")
# ---
mp2 = "# MP2 MP3\n F F\n"
cc = "# CCD pCCD DCD CCSD CCSD(T)\n F F F F F\n"
rcc = "# drCCD rCCD crCCD lCCD\n F F F F\n"
ci = "# CIS CIS(D) CID CISD FCI\n F F F F F\n"
rpa = "# phRPA phRPAx crRPA ppRPA\n F F F F\n"
gf = "# G0F2 evGF2 qsGF2 ufGF2 G0F3 evGF3\n F F F F F F\n"
gw = "# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW\n F F F F F F\n"
gtpp = "# G0T0pp evGTpp qsGTpp ufG0T0pp\n F F F F\n"
gteh = "# G0T0eh evGTeh qsGTeh\n F F F\n"
tests = "# Rtest Utest Gtest\n F F F\n"
# ---
hf_opt = "# HF: maxSCF thresh DIIS guess mix shift stab search\n 256 0.00001 5 1 0.0 0.0 F F\n"
mp_opt = "# MP: reg\n F\n"
cc_opt = "# CC: maxSCF thresh DIIS\n 64 0.00001 5\n"
tda_opt = "# spin: TDA singlet triplet\n F T T\n"
gf_opt = "# GF: maxSCF thresh DIIS lin eta renorm reg\n 256 0.00001 5 F 0.0 0 F\n"
gw_opt = "# GW: maxSCF thresh DIIS lin eta TDA_W reg\n 256 0.00001 5 F 0.0 F F\n"
gt_opt = "# GT: maxSCF thresh DIIS lin eta TDA_T reg\n 256 0.00001 5 F 0.0 F F\n"
acfdt_opt = "# ACFDT: AC Kx XBS\n F F T\n"
bse_opt = "# BSE: phBSE phBSE2 ppBSE dBSE dTDA\n F F F F T\n"
list_opt = [hf_opt, mp_opt, cc_opt, tda_opt, gf_opt, gw_opt, gt_opt, acfdt_opt, bse_opt]
# ---
class class_RHF:
def gen_input():
f = open("methods", "w")
f.write("# RHF UHF GHF ROHF\n")
f.write(" T F F F\n")
f.write("{}{}{}{}{}{}{}{}{}{}".format(mp2, cc, rcc, ci, rpa, gf, gw, gtpp, gteh, tests))
f.close()
f = open("options", "w")
for opt in list_opt:
f.write("{}".format(opt))
f.close()
def run_job(file_out, mol, bas, multip):
os.chdir('..')
print(f" :$ cd ..")
for file_in in ["methods", "options"]:
command = ['cp', 'tests/{}'.format(file_in), 'input/{}'.format(file_in)]
print(f" :$ {' '.join(command)}")
result = subprocess.run(command, capture_output=True, text=True)
if result.returncode != 0:
print("Error moving file: {}".format(result.stderr))
command = [
'python{}'.format(PYTHON_VERSION), 'PyDuck.py',
'-x', '{}'.format(mol),
'-b', '{}'.format(bas),
'-m', '{}'.format(multip)
]
print(f" :$ {' '.join(command)}")
with open(file_out, 'w') as fobj:
result = subprocess.run(command, stdout=fobj, stderr=subprocess.PIPE, text=True)
if result.stderr:
print("Error output:", result.stderr)
os.chdir('tests')
print(f" :$ cd tests")
# ---
class class_UHF:
def gen_input():
f = open("methods", "w")
f.write("# RHF UHF GHF ROHF\n")
f.write(" F T F F\n")
f.write("{}{}{}{}{}{}{}{}{}{}".format(mp2, cc, rcc, ci, rpa, gf, gw, gtpp, gteh, tests))
f.close()
# ---
class class_GHF:
def gen_input():
f = open("methods", "w")
f.write("# RHF UHF GHF ROHF\n")
f.write(" F F T F\n")
f.write("{}{}{}{}{}{}{}{}{}{}".format(mp2, cc, rcc, ci, rpa, gf, gw, gtpp, gteh, tests))
f.close()
# ---
class class_ROHF:
def gen_input():
f = open("methods", "w")
f.write("# RHF UHF GHF ROHF\n")
f.write(" F F F T\n")
f.write("{}{}{}{}{}{}{}{}{}{}".format(mp2, cc, rcc, ci, rpa, gf, gw, gtpp, gteh, tests))
f.close()
# ---
class_map = {
"RHF": class_RHF,
"UHF": class_UHF,
"GHF": class_GHF,
"ROHF": class_ROHF,
}
def main():
work_path = Path('{}/tests/work'.format(quack_root))
if not work_path.exists():
work_path.mkdir(parents=True, exist_ok=True)
print(f"Directory '{work_path}' created.\n")
for mol in molecules:
mol_name = mol.name
mol_mult = mol.multiplicity
for methd in list_methd:
if methd not in mol.energies:
print(f"Method {methd} does not exist for {mol_name}.")
continue
for bas, _ in mol.energies[methd].items():
work_methd = Path('{}/{}'.format(work_path, methd))
if not work_methd.exists():
work_methd.mkdir(parents=True, exist_ok=True)
print(f"Directory '{work_methd}' created.\n")
class_methd = class_map.get(methd)
# create input files
class_methd.gen_input()
file_out = "{}/{}/{}_{}_{}.out".format(work_path, methd, mol_name, mol_mult, bas)
print(" testing {} for {}@{} (2S+1 = {})".format(methd, mol_name, bas, mol_mult))
print(" file_out: {}".format(file_out))
class_methd.run_job(file_out, mol_name, bas, mol_mult)
print("\n")
print("\n\n")
print(" --- --- --- ---")
print("\n\n\n")
db_name = 'molecules.db'
molecules = get_molecules_from_db(db_name)
list_methd = ["RHF", "UHF", "GHF", "ROHF"]
main()

39
tests/utils.py Normal file
View File

@ -0,0 +1,39 @@
def print_col(text, color):
if(color == "black"):
print("\033[30m{}\033[0m".format(text))
elif(color == "red"):
print("\033[31m{}\033[0m".format(text))
elif(color == "green"):
print("\033[32m{}\033[0m".format(text))
elif(color == "yellow"):
print("\033[33m{}\033[0m".format(text))
elif(color == "blue"):
print("\033[34m{}\033[0m".format(text))
elif(color == "magenta"):
print("\033[35m{}\033[0m".format(text))
elif(color == "cyan"):
print("\033[36m{}\033[0m".format(text))
elif(color == "white"):
print("\033[37m{}\033[0m".format(text))
else:
print("{}".format(text))