Quantum Chemistry Refinement
This notebook demonstrates using Auto3D structures as starting points for DFT calculations:
Input generation - Gaussian, ORCA, Psi4 formats
Basis set selection - appropriate for your accuracy needs
Property calculations - energies, frequencies, NMR
Result analysis - parsing output files
Why Use Auto3D for QM?
Better starting geometries → fewer optimization steps
Pre-filtered conformers → avoid redundant calculations
NNP energies → initial ranking before expensive QM
[ ]:
import os
import tempfile
from pathlib import Path
import Auto3D
from Auto3D import Auto3DOptions, main
from rdkit import Chem
from rdkit.Chem import AllChem
print(f"Auto3D version: {Auto3D.__version__}")
1. Generate Optimized Structures
[ ]:
# Molecules for QM calculation
molecules = {
"caffeine": "CN1C=NC2=C1C(=O)N(C(=O)N2C)C",
"aspirin": "CC(=O)OC1=CC=CC=C1C(=O)O",
}
with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
for name, smi in molecules.items():
f.write(f"{smi} {name}\n")
input_file = f.name
# Generate tight geometry
if __name__ == "__main__":
config = Auto3DOptions(
path=input_file,
k=1,
optimizing_engine="AIMNET",
use_gpu=True,
convergence_threshold=0.003, # Tight for QM starting point
)
sdf_output = main(config)
print(f"Output: {sdf_output}")
2. Gaussian Input Generation
[ ]:
def mol_to_gaussian(mol, job_type="opt freq", method="B3LYP", basis="6-31G(d)",
charge=0, multiplicity=1, nproc=8, mem="16GB"):
"""
Generate Gaussian input file from RDKit molecule.
Args:
mol: RDKit molecule with 3D coordinates
job_type: Gaussian job type (opt, freq, nmr, etc.)
method: DFT functional or method
basis: Basis set
charge: Molecular charge
multiplicity: Spin multiplicity
nproc: Number of processors
mem: Memory allocation
"""
name = mol.GetProp("_Name") if mol.HasProp("_Name") else "molecule"
conf = mol.GetConformer()
# Header
lines = [
f"%NProcShared={nproc}",
f"%Mem={mem}",
f"%Chk={name}.chk",
f"# {method}/{basis} {job_type}",
"",
f"{name} - Auto3D optimized geometry",
"",
f"{charge} {multiplicity}",
]
# Coordinates
for atom in mol.GetAtoms():
pos = conf.GetAtomPosition(atom.GetIdx())
symbol = atom.GetSymbol()
lines.append(f"{symbol:2s} {pos.x:12.6f} {pos.y:12.6f} {pos.z:12.6f}")
lines.append("") # Blank line at end
return "\n".join(lines)
# Generate Gaussian inputs
if 'sdf_output' in dir():
mols = list(Chem.SDMolSupplier(sdf_output, removeHs=False))
for mol in mols:
if mol is None:
continue
name = mol.GetProp("_Name").split("_")[0]
charge = Chem.GetFormalCharge(mol)
# Geometry optimization + frequency
gjf_opt = mol_to_gaussian(mol, "opt freq", "B3LYP", "6-311+G(2d,p)", charge)
print(f"\n=== {name}.gjf ===")
print(gjf_opt[:500] + "...")
3. ORCA Input Generation
[ ]:
def mol_to_orca(mol, job_type="Opt Freq", method="B3LYP", basis="def2-TZVP",
charge=0, multiplicity=1, nproc=8):
"""
Generate ORCA input file from RDKit molecule.
"""
name = mol.GetProp("_Name") if mol.HasProp("_Name") else "molecule"
conf = mol.GetConformer()
# Header
lines = [
f"# ORCA input - {name}",
f"! {method} {basis} {job_type} D3BJ", # Include dispersion
f"%pal nprocs {nproc} end",
"",
f"* xyz {charge} {multiplicity}",
]
# Coordinates
for atom in mol.GetAtoms():
pos = conf.GetAtomPosition(atom.GetIdx())
symbol = atom.GetSymbol()
lines.append(f" {symbol:2s} {pos.x:12.6f} {pos.y:12.6f} {pos.z:12.6f}")
lines.append("*")
return "\n".join(lines)
# Generate ORCA input
if 'sdf_output' in dir():
mol = next(Chem.SDMolSupplier(sdf_output, removeHs=False))
if mol:
orca_input = mol_to_orca(mol)
print("=== ORCA Input ===")
print(orca_input)
4. Psi4 Input (Python Script)
[ ]:
def mol_to_psi4(mol, method="B3LYP", basis="cc-pVDZ", charge=0, multiplicity=1):
"""
Generate Psi4 Python script from RDKit molecule.
"""
name = mol.GetProp("_Name") if mol.HasProp("_Name") else "molecule"
conf = mol.GetConformer()
# Build geometry string
geom_lines = []
for atom in mol.GetAtoms():
pos = conf.GetAtomPosition(atom.GetIdx())
symbol = atom.GetSymbol()
geom_lines.append(f" {symbol:2s} {pos.x:12.6f} {pos.y:12.6f} {pos.z:12.6f}")
geometry = "\n".join(geom_lines)
script = f'''
"""Psi4 calculation for {name}"""
import psi4
psi4.set_memory("16 GB")
psi4.set_num_threads(8)
mol = psi4.geometry("""
{charge} {multiplicity}
{geometry}
""")
# Geometry optimization
psi4.set_options({{
"basis": "{basis}",
"reference": "rhf" if {multiplicity} == 1 else "uhf",
"scf_type": "df",
"d_convergence": 1e-8,
}})
E_opt = psi4.optimize("{method}", molecule=mol)
print(f"Optimized energy: {{E_opt:.6f}} Hartree")
# Frequency calculation
E_freq, wfn = psi4.frequency("{method}", molecule=mol, return_wfn=True)
freqs = wfn.frequencies().to_array()
print(f"Frequencies (cm-1): {{freqs}}")
'''
return script.strip()
if 'sdf_output' in dir():
mol = next(Chem.SDMolSupplier(sdf_output, removeHs=False))
if mol:
psi4_script = mol_to_psi4(mol)
print("=== Psi4 Script ===")
print(psi4_script[:800] + "...")
5. Method/Basis Set Recommendations
[ ]:
print("""
QM Method Recommendations for Drug-Like Molecules
=================================================
GEOMETRY OPTIMIZATION:
Fast: B3LYP/6-31G(d) (~1 min/mol)
Standard: B3LYP-D3/6-311G(d,p) (~5 min/mol)
Accurate: ωB97X-D/def2-TZVP (~30 min/mol)
SINGLE-POINT ENERGY:
Standard: B3LYP-D3/def2-TZVP
Accurate: DLPNO-CCSD(T)/cc-pVTZ (ORCA)
Gold std: CCSD(T)/CBS extrapolation
THERMOCHEMISTRY:
Method: B3LYP-D3/6-311+G(2d,2p) with frequency
Scale: 0.967 for B3LYP frequencies
NMR CALCULATIONS:
1H/13C: mPW1PW91/6-311+G(2d,p) with GIAO
or: B3LYP/6-311+G(2d,p) with GIAO
DISPERSION:
Always include for drug-like molecules:
- D3(BJ) for Gaussian/ORCA
- DFT-D3 for Psi4
- Built-in for ωB97X-D, ωB97M-V
""")
6. Batch Input Generation
[ ]:
def generate_batch_inputs(sdf_path, output_dir, software="gaussian"):
"""
Generate QM inputs for all molecules in SDF.
"""
output_dir = Path(output_dir)
output_dir.mkdir(exist_ok=True)
mols = list(Chem.SDMolSupplier(sdf_path, removeHs=False))
generated = []
for mol in mols:
if mol is None:
continue
name = mol.GetProp("_Name").replace(" ", "_")
charge = Chem.GetFormalCharge(mol)
if software.lower() == "gaussian":
content = mol_to_gaussian(mol, charge=charge)
ext = ".gjf"
elif software.lower() == "orca":
content = mol_to_orca(mol, charge=charge)
ext = ".inp"
elif software.lower() == "psi4":
content = mol_to_psi4(mol, charge=charge)
ext = ".py"
else:
raise ValueError(f"Unknown software: {software}")
output_file = output_dir / f"{name}{ext}"
with open(output_file, 'w') as f:
f.write(content)
generated.append(str(output_file))
return generated
print("Batch generation function available.")
print("Usage: generate_batch_inputs('structures.sdf', 'qm_inputs/', 'gaussian')")
7. Output Parsing (Example)
[ ]:
def parse_gaussian_log(log_path):
"""
Parse Gaussian log file for key results.
"""
results = {
"converged": False,
"energy": None,
"frequencies": [],
"imaginary_freqs": 0,
}
with open(log_path) as f:
for line in f:
if "Normal termination" in line:
results["converged"] = True
elif "SCF Done" in line:
results["energy"] = float(line.split()[4])
elif "Frequencies --" in line:
freqs = [float(x) for x in line.split()[2:]]
results["frequencies"].extend(freqs)
results["imaginary_freqs"] += sum(1 for f in freqs if f < 0)
return results
print("Output parsing functions available.")
print("\nFor comprehensive parsing, consider:")
print(" - cclib: pip install cclib")
print(" - ASE: ase.io.read()")
print(" - GaussView or Avogadro for visualization")
Summary
This tutorial covered:
Auto3D → QM workflow - NNP-optimized starting geometries
Gaussian input file generation
ORCA input file generation
Psi4 Python script generation
Method selection - appropriate basis sets and functionals
Batch processing - generating inputs for multiple molecules
Key benefits:
Fewer SCF cycles with better starting geometry
Pre-ranking with NNP energies
Conformer filtering before expensive QM
[ ]:
# Cleanup
if 'input_file' in dir() and os.path.exists(input_file):
os.unlink(input_file)