Вычисление параметров для молекулярной механики¶
подготовка
In [1]:
from pathlib import Path
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import psi4
Matplotlib is building the font cache; this may take a moment.
In [2]:
# Главная папка проекта
project_dir = Path.home() / "Documents" / "jupyter" / "mm_ethane"
# Подпапки
results_dir = project_dir / "results"
logs_dir = results_dir / "psi4_logs"
xyz_dir = results_dir / "xyz"
plots_dir = results_dir / "plots"
orca_dir = results_dir / "orca_inputs" # опционально, если захочешь ещё и ORCA-входы
scratch_dir = project_dir / "scratch"
for folder in [project_dir, results_dir, logs_dir, xyz_dir, plots_dir, orca_dir, scratch_dir]:
folder.mkdir(parents=True, exist_ok=True)
# Psi4 scratch
os.environ["PSI_SCRATCH"] = str(scratch_dir)
# Настройки Psi4
psi4.set_memory("4 GB")
psi4.set_num_threads(4)
print("Проект:", project_dir)
print("Результаты:", results_dir)
print("Scratch:", scratch_dir)
Проект: Memory set to 3.725 GiB by Python driver. /Users/romashka/Documents/jupyter/mm_ethane Результаты: /Users/romashka/Documents/jupyter/mm_ethane/results Scratch: /Users/romashka/Documents/jupyter/mm_ethane/scratch Threads set to 4 by Python driver.
Начало основного блока заданий¶
In [3]:
R0_CC = 1.52986 # исходная длина C-C, Å
R_CH = 1.08439 # длина C-H, Å
A_HCC = 111.200 # угол H-C-C, градусы
In [4]:
def make_ethane_zmat(rcc, rch=R_CH, ahcc=A_HCC):
return f"""
0 1
C
C 1 rCC
H 1 rCH 2 aHCC
H 1 rCH 2 aHCC 3 120.0
H 1 rCH 2 aHCC 3 -120.0
H 2 rCH 1 aHCC 3 180.0
H 2 rCH 1 aHCC 5 120.0
H 2 rCH 1 aHCC 5 -120.0
rCC = {rcc:.5f}
rCH = {rch:.5f}
aHCC = {ahcc:.3f}
units angstrom
symmetry c1
"""
In [5]:
print(make_ethane_zmat(R0_CC))
0 1 C C 1 rCC H 1 rCH 2 aHCC H 1 rCH 2 aHCC 3 120.0 H 1 rCH 2 aHCC 3 -120.0 H 2 rCH 1 aHCC 3 180.0 H 2 rCH 1 aHCC 5 120.0 H 2 rCH 1 aHCC 5 -120.0 rCC = 1.52986 rCH = 1.08439 aHCC = 111.200 units angstrom symmetry c1
Здесь мы подготовили основу всего задания: шаблон молекулы, где можно менять только одну длину связи.
In [7]:
import py3Dmol
m0 = psi4.geometry(make_ethane_zmat(R0_CC))
view = py3Dmol.view(width=500, height=350)
view.addModel(m0.save_string_xyz_file(), "xyz")
view.setStyle({"stick": {}})
view.zoomTo()
view.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
Убедились, что геометрия собрана правильно.
In [8]:
METHOD_BASIS = "scf/cc-pvtz"
# Для быстрого теста можно временно поставить:
# METHOD_BASIS = "hf/6-31g"
HARTREE_TO_KCAL = 627.509474
In [9]:
def safe_label_from_r(r):
return f"ethane_cc_r_{r:.5f}".replace(".", "p")
In [10]:
def run_psi4(zmat_text, label, method_basis=METHOD_BASIS):
out_file = logs_dir / f"{label}.out"
xyz_file = xyz_dir / f"{label}.xyz"
zmat_file = xyz_dir / f"{label}.zmat.txt"
result = {
"label": label,
"method_basis": method_basis,
"out_file": out_file.relative_to(project_dir).as_posix(),
"xyz_file": xyz_file.relative_to(project_dir).as_posix(),
"zmat_file": zmat_file.relative_to(project_dir).as_posix(),
"status": "ok"
}
try:
psi4.core.clean()
psi4.core.set_output_file(str(out_file), False)
psi4.set_options({
"reference": "rhf",
"maxiter": 200,
"fail_on_maxiter": True
})
mol = psi4.geometry(zmat_text)
energy = psi4.energy(method_basis, molecule=mol)
xyz_text = mol.save_string_xyz_file()
xyz_file.write_text(xyz_text, encoding="utf-8")
zmat_file.write_text(zmat_text, encoding="utf-8")
result["energy_hartree"] = float(energy)
except Exception as e:
result["energy_hartree"] = np.nan
result["status"] = f"failed: {e}"
finally:
psi4.core.clean()
return result
In [11]:
# Проверка на одной точке
test_label = safe_label_from_r(R0_CC)
test_result = run_psi4(make_ethane_zmat(R0_CC), test_label)
test_result
Out[11]:
{'label': 'ethane_cc_r_1p52986',
'method_basis': 'scf/cc-pvtz',
'out_file': 'results/psi4_logs/ethane_cc_r_1p52986.out',
'xyz_file': 'results/xyz/ethane_cc_r_1p52986.xyz',
'zmat_file': 'results/xyz/ethane_cc_r_1p52986.zmat.txt',
'status': 'ok',
'energy_hartree': -79.15464083093346}
Теперь вместо ручного запуска каждого случая мы сможем просто прогнать цикл по разным длинам связи и собрать все энергии в одну таблицу.
In [12]:
step = 0.02
r_values = R0_CC + step * np.r_[np.arange(-10, 0), np.arange(1, 11)]
r_values = np.round(r_values, 5)
r_values
Out[12]:
array([1.32986, 1.34986, 1.36986, 1.38986, 1.40986, 1.42986, 1.44986,
1.46986, 1.48986, 1.50986, 1.54986, 1.56986, 1.58986, 1.60986,
1.62986, 1.64986, 1.66986, 1.68986, 1.70986, 1.72986])
In [13]:
records = []
for rcc in r_values:
label = safe_label_from_r(rcc)
zmat = make_ethane_zmat(rcc)
res = run_psi4(zmat, label)
res["r_cc_angstrom"] = rcc
records.append(res)
df = pd.DataFrame(records).sort_values("r_cc_angstrom").reset_index(drop=True)
df["dE_hartree"] = df["energy_hartree"] - df["energy_hartree"].min()
df["dE_kcal_mol"] = df["dE_hartree"] * HARTREE_TO_KCAL
summary_csv = results_dir / "scan_summary.csv"
df.to_csv(summary_csv, index=False)
df
Out[13]:
| label | method_basis | out_file | xyz_file | zmat_file | status | energy_hartree | r_cc_angstrom | dE_hartree | dE_kcal_mol | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | ethane_cc_r_1p32986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p32986.out | results/xyz/ethane_cc_r_1p32986.xyz | results/xyz/ethane_cc_r_1p32986.zmat.txt | ok | -79.121471 | 1.32986 | 0.033035 | 20.730026 |
| 1 | ethane_cc_r_1p34986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p34986.out | results/xyz/ethane_cc_r_1p34986.xyz | results/xyz/ethane_cc_r_1p34986.zmat.txt | ok | -79.128626 | 1.34986 | 0.025881 | 16.240306 |
| 2 | ethane_cc_r_1p36986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p36986.out | results/xyz/ethane_cc_r_1p36986.xyz | results/xyz/ethane_cc_r_1p36986.zmat.txt | ok | -79.134710 | 1.36986 | 0.019797 | 12.422542 |
| 3 | ethane_cc_r_1p38986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p38986.out | results/xyz/ethane_cc_r_1p38986.xyz | results/xyz/ethane_cc_r_1p38986.zmat.txt | ok | -79.139817 | 1.38986 | 0.014689 | 9.217487 |
| 4 | ethane_cc_r_1p40986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p40986.out | results/xyz/ethane_cc_r_1p40986.xyz | results/xyz/ethane_cc_r_1p40986.zmat.txt | ok | -79.144035 | 1.40986 | 0.010471 | 6.570947 |
| 5 | ethane_cc_r_1p42986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p42986.out | results/xyz/ethane_cc_r_1p42986.xyz | results/xyz/ethane_cc_r_1p42986.zmat.txt | ok | -79.147442 | 1.42986 | 0.007065 | 4.433321 |
| 6 | ethane_cc_r_1p44986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p44986.out | results/xyz/ethane_cc_r_1p44986.xyz | results/xyz/ethane_cc_r_1p44986.zmat.txt | ok | -79.150109 | 1.44986 | 0.004397 | 2.759187 |
| 7 | ethane_cc_r_1p46986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p46986.out | results/xyz/ethane_cc_r_1p46986.xyz | results/xyz/ethane_cc_r_1p46986.zmat.txt | ok | -79.152105 | 1.46986 | 0.002401 | 1.506933 |
| 8 | ethane_cc_r_1p48986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p48986.out | results/xyz/ethane_cc_r_1p48986.xyz | results/xyz/ethane_cc_r_1p48986.zmat.txt | ok | -79.153489 | 1.48986 | 0.001017 | 0.638420 |
| 9 | ethane_cc_r_1p50986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p50986.out | results/xyz/ethane_cc_r_1p50986.xyz | results/xyz/ethane_cc_r_1p50986.zmat.txt | ok | -79.154317 | 1.50986 | 0.000189 | 0.118688 |
| 10 | ethane_cc_r_1p54986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p54986.out | results/xyz/ethane_cc_r_1p54986.xyz | results/xyz/ethane_cc_r_1p54986.zmat.txt | ok | -79.154506 | 1.54986 | 0.000000 | 0.000000 |
| 11 | ethane_cc_r_1p56986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p56986.out | results/xyz/ethane_cc_r_1p56986.xyz | results/xyz/ethane_cc_r_1p56986.zmat.txt | ok | -79.153957 | 1.56986 | 0.000549 | 0.344691 |
| 12 | ethane_cc_r_1p58986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p58986.out | results/xyz/ethane_cc_r_1p58986.xyz | results/xyz/ethane_cc_r_1p58986.zmat.txt | ok | -79.153032 | 1.58986 | 0.001474 | 0.925032 |
| 13 | ethane_cc_r_1p60986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p60986.out | results/xyz/ethane_cc_r_1p60986.xyz | results/xyz/ethane_cc_r_1p60986.zmat.txt | ok | -79.151768 | 1.60986 | 0.002738 | 1.718356 |
| 14 | ethane_cc_r_1p62986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p62986.out | results/xyz/ethane_cc_r_1p62986.xyz | results/xyz/ethane_cc_r_1p62986.zmat.txt | ok | -79.150198 | 1.62986 | 0.004309 | 2.703881 |
| 15 | ethane_cc_r_1p64986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p64986.out | results/xyz/ethane_cc_r_1p64986.xyz | results/xyz/ethane_cc_r_1p64986.zmat.txt | ok | -79.148351 | 1.64986 | 0.006155 | 3.862557 |
| 16 | ethane_cc_r_1p66986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p66986.out | results/xyz/ethane_cc_r_1p66986.xyz | results/xyz/ethane_cc_r_1p66986.zmat.txt | ok | -79.146257 | 1.66986 | 0.008250 | 5.176922 |
| 17 | ethane_cc_r_1p68986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p68986.out | results/xyz/ethane_cc_r_1p68986.xyz | results/xyz/ethane_cc_r_1p68986.zmat.txt | ok | -79.143939 | 1.68986 | 0.010567 | 6.630977 |
| 18 | ethane_cc_r_1p70986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p70986.out | results/xyz/ethane_cc_r_1p70986.xyz | results/xyz/ethane_cc_r_1p70986.zmat.txt | ok | -79.141423 | 1.70986 | 0.013084 | 8.210063 |
| 19 | ethane_cc_r_1p72986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p72986.out | results/xyz/ethane_cc_r_1p72986.xyz | results/xyz/ethane_cc_r_1p72986.zmat.txt | ok | -79.138729 | 1.72986 | 0.015778 | 9.900755 |
Автоматически пробегаем по 20 длинам связи, для каждой запускаем Psi4, сохраняем файлы и собираем все энергии в таблицу df. Дополнительно считаем относительную энергию, чтобы график было легче читать.
In [14]:
# Собираем все xyz в один общий файл
multi_xyz_file = results_dir / "ethane_cc_scan_20.xyz"
with multi_xyz_file.open("w", encoding="utf-8") as fout:
for xyz_rel_path in df["xyz_file"]:
xyz_path = project_dir / xyz_rel_path
fout.write(xyz_path.read_text(encoding="utf-8").strip())
fout.write("\n")
print("Общий xyz-файл:", multi_xyz_file)
Общий xyz-файл: /Users/romashka/Documents/jupyter/mm_ethane/results/ethane_cc_scan_20.xyz
Определяем константу связи через квадратичный fit около минимума¶
In [15]:
r_min_scan = df.loc[df["energy_hartree"].idxmin(), "r_cc_angstrom"]
fit_mask = np.abs(df["r_cc_angstrom"] - r_min_scan) <= 0.08
fit_df = df.loc[fit_mask].copy()
fit_df
Out[15]:
| label | method_basis | out_file | xyz_file | zmat_file | status | energy_hartree | r_cc_angstrom | dE_hartree | dE_kcal_mol | |
|---|---|---|---|---|---|---|---|---|---|---|
| 8 | ethane_cc_r_1p48986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p48986.out | results/xyz/ethane_cc_r_1p48986.xyz | results/xyz/ethane_cc_r_1p48986.zmat.txt | ok | -79.153489 | 1.48986 | 0.001017 | 0.638420 |
| 9 | ethane_cc_r_1p50986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p50986.out | results/xyz/ethane_cc_r_1p50986.xyz | results/xyz/ethane_cc_r_1p50986.zmat.txt | ok | -79.154317 | 1.50986 | 0.000189 | 0.118688 |
| 10 | ethane_cc_r_1p54986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p54986.out | results/xyz/ethane_cc_r_1p54986.xyz | results/xyz/ethane_cc_r_1p54986.zmat.txt | ok | -79.154506 | 1.54986 | 0.000000 | 0.000000 |
| 11 | ethane_cc_r_1p56986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p56986.out | results/xyz/ethane_cc_r_1p56986.xyz | results/xyz/ethane_cc_r_1p56986.zmat.txt | ok | -79.153957 | 1.56986 | 0.000549 | 0.344691 |
| 12 | ethane_cc_r_1p58986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p58986.out | results/xyz/ethane_cc_r_1p58986.xyz | results/xyz/ethane_cc_r_1p58986.zmat.txt | ok | -79.153032 | 1.58986 | 0.001474 | 0.925032 |
| 13 | ethane_cc_r_1p60986 | scf/cc-pvtz | results/psi4_logs/ethane_cc_r_1p60986.out | results/xyz/ethane_cc_r_1p60986.xyz | results/xyz/ethane_cc_r_1p60986.zmat.txt | ok | -79.151768 | 1.60986 | 0.002738 | 1.718356 |
In [16]:
a, b, c = np.polyfit(fit_df["r_cc_angstrom"], fit_df["energy_hartree"], 2)
r_eq_fit = -b / (2 * a)
e0_fit = a * r_eq_fit**2 + b * r_eq_fit + c
# 1) Если писать потенциал как E = 1/2 * k * (r - r0)^2
k_half_hartree = 2 * a
k_half_kcal = k_half_hartree * HARTREE_TO_KCAL
# 2) Если писать потенциал как E = k * (r - r0)^2 (частая форма в MM)
k_mm_hartree = a
k_mm_kcal = k_mm_hartree * HARTREE_TO_KCAL
report_text = f"""
Метод: {METHOD_BASIS}
Квадратичный fit: E(r) = a*r^2 + b*r + c
a = {a:.8f} Hartree/Å^2
b = {b:.8f} Hartree/Å
c = {c:.8f} Hartree
r_eq (из fit) = {r_eq_fit:.5f} Å
E0 (из fit) = {e0_fit:.8f} Hartree
Если использовать запись:
E = 1/2 * k * (r - r0)^2
то:
k = {k_half_hartree:.8f} Hartree/Å^2
k = {k_half_kcal:.4f} kcal/(mol·Å^2)
Если использовать запись:
E = k * (r - r0)^2
то:
k = {k_mm_hartree:.8f} Hartree/Å^2
k = {k_mm_kcal:.4f} kcal/(mol·Å^2)
"""
force_report = results_dir / "force_constant_report.txt"
force_report.write_text(report_text, encoding="utf-8")
print(report_text)
Метод: scf/cc-pvtz Квадратичный fit: E(r) = a*r^2 + b*r + c a = 0.51672280 Hartree/Å^2 b = -1.58676880 Hartree/Å c = -77.93642438 Hartree r_eq (из fit) = 1.53542 Å E0 (из fit) = -79.15459940 Hartree Если использовать запись: E = 1/2 * k * (r - r0)^2 то: k = 1.03344560 Hartree/Å^2 k = 648.4969 kcal/(mol·Å^2) Если использовать запись: E = k * (r - r0)^2 то: k = 0.51672280 Hartree/Å^2 k = 324.2485 kcal/(mol·Å^2)
In [17]:
x_dense = np.linspace(fit_df["r_cc_angstrom"].min() - 0.01,
fit_df["r_cc_angstrom"].max() + 0.01, 300)
y_dense = np.polyval([a, b, c], x_dense)
y_dense_rel = (y_dense - df["energy_hartree"].min()) * HARTREE_TO_KCAL
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.scatter(df["r_cc_angstrom"], df["dE_kcal_mol"], label="Все точки скана")
ax.scatter(fit_df["r_cc_angstrom"], fit_df["dE_kcal_mol"], label="Точки для fit")
ax.plot(x_dense, y_dense_rel, label="Квадратичный fit")
ax.axvline(r_eq_fit, linestyle="--", label=f"r0 = {r_eq_fit:.4f} Å")
ax.set_xlabel("Длина связи C–C, Å")
ax.set_ylabel("Относительная энергия, kcal/mol")
ax.set_title("Квадратичный fit около минимума")
ax.grid(alpha=0.3)
ax.legend()
fig.tight_layout()
fit_plot = plots_dir / "cc_scan_fit.png"
fig.savefig(fit_plot, dpi=250, bbox_inches="tight")
plt.show()
Мы приближаем участок около минимума параболой, потому что вблизи равновесия потенциальная поверхность хорошо описывается гармоническим потенциалом. Именно из коэффициента при квадрате и получают силовую константу связи.
Ссылки на все файлы результатов¶
In [51]:
from IPython.display import HTML, display
In [52]:
links_html = """
<h3>Файлы результатов</h3>
<ul>
<li><a href="results/scan_summary.csv">Сводная таблица scan_summary.csv</a></li>
<li><a href="results/force_constant_report.txt">Текстовый отчёт по fit</a></li>
<li><a href="results/plots/cc_scan_relative_energy.png">График зависимости энергии от длины связи</a></li>
<li><a href="results/plots/cc_scan_fit.png">График с квадратичным fit</a></li>
<li><a href="results/ethane_cc_scan_20.xyz">Общий файл со всеми 20 геометриями</a></li>
</ul>
"""
snippet_file = project_dir / "results_links_snippet.html"
snippet_file.write_text(links_html, encoding="utf-8")
display(HTML(links_html))
Вычисления для валентного угла HCC¶
In [23]:
# 1.1. Z-matrix для скана по углу H-C-C
def make_ethane_angle_zmat(ahcc, rcc=R0_CC, rch=R_CH):
return f"""
0 1
C
C 1 rCC
H 1 rCH 2 aHCC
H 1 rCH 2 aHCC 3 120.0
H 1 rCH 2 aHCC 3 -120.0
H 2 rCH 1 aHCC 3 180.0
H 2 rCH 1 aHCC 5 120.0
H 2 rCH 1 aHCC 5 -120.0
rCC = {rcc:.5f}
rCH = {rch:.5f}
aHCC = {ahcc:.3f}
units angstrom
symmetry c1
"""
In [25]:
# Удобная подпись файлов
def safe_angle_label(angle):
return f"ethane_hcc_{angle:.3f}".replace(".", "p")
In [26]:
# Значения угла
angle_values = np.arange(109.2, 113.2 + 0.0001, 0.2)
angle_values
Out[26]:
array([109.2, 109.4, 109.6, 109.8, 110. , 110.2, 110.4, 110.6, 110.8,
111. , 111.2, 111.4, 111.6, 111.8, 112. , 112.2, 112.4, 112.6,
112.8, 113. , 113.2])
In [27]:
# Запуск расчётов
angle_records = []
for angle in angle_values:
label = safe_angle_label(angle)
zmat = make_ethane_angle_zmat(angle)
res = run_psi4(zmat, label)
res["hcc_angle_deg"] = angle
angle_records.append(res)
df_angle = pd.DataFrame(angle_records).sort_values("hcc_angle_deg").reset_index(drop=True)
df_angle["dE_hartree"] = df_angle["energy_hartree"] - df_angle["energy_hartree"].min()
df_angle["dE_kcal_mol"] = df_angle["dE_hartree"] * HARTREE_TO_KCAL
angle_csv = results_dir / "angle_hcc_scan_summary.csv"
df_angle.to_csv(angle_csv, index=False)
df_angle
Out[27]:
| label | method_basis | out_file | xyz_file | zmat_file | status | energy_hartree | hcc_angle_deg | dE_hartree | dE_kcal_mol | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | ethane_hcc_109p200 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_109p200.out | results/xyz/ethane_hcc_109p200.xyz | results/xyz/ethane_hcc_109p200.zmat.txt | ok | -79.151884 | 109.2 | 0.003238 | 2.031595 |
| 1 | ethane_hcc_109p400 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_109p400.out | results/xyz/ethane_hcc_109p400.xyz | results/xyz/ethane_hcc_109p400.zmat.txt | ok | -79.152269 | 109.4 | 0.002853 | 1.790090 |
| 2 | ethane_hcc_109p600 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_109p600.out | results/xyz/ethane_hcc_109p600.xyz | results/xyz/ethane_hcc_109p600.zmat.txt | ok | -79.152629 | 109.6 | 0.002492 | 1.563745 |
| 3 | ethane_hcc_109p800 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_109p800.out | results/xyz/ethane_hcc_109p800.xyz | results/xyz/ethane_hcc_109p800.zmat.txt | ok | -79.152966 | 109.8 | 0.002155 | 1.352582 |
| 4 | ethane_hcc_110p000 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_110p000.out | results/xyz/ethane_hcc_110p000.xyz | results/xyz/ethane_hcc_110p000.zmat.txt | ok | -79.153278 | 110.0 | 0.001843 | 1.156625 |
| 5 | ethane_hcc_110p200 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_110p200.out | results/xyz/ethane_hcc_110p200.xyz | results/xyz/ethane_hcc_110p200.zmat.txt | ok | -79.153566 | 110.2 | 0.001555 | 0.975895 |
| 6 | ethane_hcc_110p400 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_110p400.out | results/xyz/ethane_hcc_110p400.xyz | results/xyz/ethane_hcc_110p400.zmat.txt | ok | -79.153830 | 110.4 | 0.001291 | 0.810416 |
| 7 | ethane_hcc_110p600 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_110p600.out | results/xyz/ethane_hcc_110p600.xyz | results/xyz/ethane_hcc_110p600.zmat.txt | ok | -79.154069 | 110.6 | 0.001052 | 0.660212 |
| 8 | ethane_hcc_110p800 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_110p800.out | results/xyz/ethane_hcc_110p800.xyz | results/xyz/ethane_hcc_110p800.zmat.txt | ok | -79.154284 | 110.8 | 0.000837 | 0.525308 |
| 9 | ethane_hcc_111p000 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_111p000.out | results/xyz/ethane_hcc_111p000.xyz | results/xyz/ethane_hcc_111p000.zmat.txt | ok | -79.154475 | 111.0 | 0.000647 | 0.405730 |
| 10 | ethane_hcc_111p200 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_111p200.out | results/xyz/ethane_hcc_111p200.xyz | results/xyz/ethane_hcc_111p200.zmat.txt | ok | -79.154641 | 111.2 | 0.000480 | 0.301501 |
| 11 | ethane_hcc_111p400 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_111p400.out | results/xyz/ethane_hcc_111p400.xyz | results/xyz/ethane_hcc_111p400.zmat.txt | ok | -79.154782 | 111.4 | 0.000339 | 0.212650 |
| 12 | ethane_hcc_111p600 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_111p600.out | results/xyz/ethane_hcc_111p600.xyz | results/xyz/ethane_hcc_111p600.zmat.txt | ok | -79.154899 | 111.6 | 0.000222 | 0.139202 |
| 13 | ethane_hcc_111p800 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_111p800.out | results/xyz/ethane_hcc_111p800.xyz | results/xyz/ethane_hcc_111p800.zmat.txt | ok | -79.154992 | 111.8 | 0.000129 | 0.081185 |
| 14 | ethane_hcc_112p000 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_112p000.out | results/xyz/ethane_hcc_112p000.xyz | results/xyz/ethane_hcc_112p000.zmat.txt | ok | -79.155060 | 112.0 | 0.000062 | 0.038626 |
| 15 | ethane_hcc_112p200 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_112p200.out | results/xyz/ethane_hcc_112p200.xyz | results/xyz/ethane_hcc_112p200.zmat.txt | ok | -79.155103 | 112.2 | 0.000018 | 0.011555 |
| 16 | ethane_hcc_112p400 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_112p400.out | results/xyz/ethane_hcc_112p400.xyz | results/xyz/ethane_hcc_112p400.zmat.txt | ok | -79.155121 | 112.4 | 0.000000 | 0.000000 |
| 17 | ethane_hcc_112p600 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_112p600.out | results/xyz/ethane_hcc_112p600.xyz | results/xyz/ethane_hcc_112p600.zmat.txt | ok | -79.155115 | 112.6 | 0.000006 | 0.003990 |
| 18 | ethane_hcc_112p800 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_112p800.out | results/xyz/ethane_hcc_112p800.xyz | results/xyz/ethane_hcc_112p800.zmat.txt | ok | -79.155084 | 112.8 | 0.000038 | 0.023556 |
| 19 | ethane_hcc_113p000 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_113p000.out | results/xyz/ethane_hcc_113p000.xyz | results/xyz/ethane_hcc_113p000.zmat.txt | ok | -79.155028 | 113.0 | 0.000094 | 0.058729 |
| 20 | ethane_hcc_113p200 | scf/cc-pvtz | results/psi4_logs/ethane_hcc_113p200.out | results/xyz/ethane_hcc_113p200.xyz | results/xyz/ethane_hcc_113p200.zmat.txt | ok | -79.154947 | 113.2 | 0.000175 | 0.109538 |
In [28]:
# Построение графика энергии от угла
plt.figure(figsize=(7, 4.5))
plt.plot(df_angle["hcc_angle_deg"], df_angle["dE_kcal_mol"], "o-")
plt.xlabel("Угол H–C–C, градусы")
plt.ylabel("Относительная энергия, kcal/mol")
plt.title("Зависимость энергии этана от валентного угла H–C–C")
plt.grid(alpha=0.3)
plt.tight_layout()
angle_plot = plots_dir / "angle_hcc_scan.png"
plt.savefig(angle_plot, dpi=250, bbox_inches="tight")
plt.show()
Получаем зависимость энергии этана от валентного угла HCC.
In [29]:
# Сохранение полученных коэффицентов
fitfunc = lambda p, x: p[0] * (p[1] - x)**2 + p[2]
errfunc = lambda p, x, y: fitfunc(p, x) - y
In [31]:
# Fit для угла
from scipy import optimize
x_angle = df_angle["hcc_angle_deg"].to_numpy()
y_angle = df_angle["energy_hartree"].to_numpy()
p0 = [1e-3, x_angle[np.argmin(y_angle)], y_angle.min()]
p_angle, success_angle = optimize.leastsq(errfunc, p0[:], args=(x_angle, y_angle))
k_angle_hartree_deg2 = p_angle[0]
theta0_deg = p_angle[1]
E0_angle = p_angle[2]
print("Optimized params for angle scan:")
print("k (Hartree/deg^2) =", k_angle_hartree_deg2)
print("theta0 (deg) =", theta0_deg)
print("E0 (Hartree) =", E0_angle)
Optimized params for angle scan: k (Hartree/deg^2) = 0.00030639022711007074 theta0 (deg) = 112.45159201706947 E0 (Hartree) = -79.15512082302072
In [32]:
# График с fitted-кривой
x_dense = np.linspace(x_angle.min(), x_angle.max(), 400)
plt.figure(figsize=(7, 4.5))
plt.plot(x_angle, (y_angle - y_angle.min()) * HARTREE_TO_KCAL, "o", label="Расчётные точки")
plt.plot(x_dense, (fitfunc(p_angle, x_dense) - y_angle.min()) * HARTREE_TO_KCAL, "-", label="Fit")
plt.xlabel("Угол H–C–C, градусы")
plt.ylabel("Относительная энергия, kcal/mol")
plt.title("Fit зависимости энергии от угла H–C–C")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
angle_fit_plot = plots_dir / "angle_hcc_fit.png"
plt.savefig(angle_fit_plot, dpi=250, bbox_inches="tight")
plt.show()
In [33]:
# Сохранение коэффициентов в текстовый файл
angle_report = results_dir / "angle_fit_report.txt"
angle_report.write_text(
f"""Fit для угла H-C-C
k = {k_angle_hartree_deg2:.8f} Hartree/deg^2
theta0 = {theta0_deg:.5f} deg
E0 = {E0_angle:.8f} Hartree
""",
encoding="utf-8"
)
print("Коэффициенты сохранены в:", angle_report)
Коэффициенты сохранены в: /Users/romashka/Documents/jupyter/mm_ethane/results/angle_fit_report.txt
In [ ]:
## Блок сравнения
Для сравнения с GAFF были использованы параметры из файла gaff.dat для типов атомов, соответствующих этану: c3-c3 для связи C–C и c3-c3-hc для угла H–C–C. Из файла gaff.dat были взяты значения K_r = 303.1 kcal/(mol·Å²) и r_eq = 1.535 Å для связи c3-c3, а также K_θ = 46.4 kcal/(mol·rad²) и θ_eq = 110.05° для угла c3-c3-hc.
In [34]:
DEG2RAD = np.pi / 180.0
RAD2DEG = 180.0 / np.pi
In [35]:
# Перевод угловой константы
k_angle_kcal_deg2 = k_angle_hartree_deg2 * HARTREE_TO_KCAL
k_angle_kcal_rad2 = k_angle_kcal_deg2 * (RAD2DEG ** 2)
print("k_angle =", k_angle_kcal_deg2, "kcal/(mol·deg^2)")
print("k_angle =", k_angle_kcal_rad2, "kcal/(mol·rad^2)")
k_angle = 0.192262770252581 kcal/(mol·deg^2) k_angle = 631.161443056022 kcal/(mol·rad^2)
In [37]:
k_bond_hartree = k_mm_hartree
k_bond_kcal = k_mm_hartree * HARTREE_TO_KCAL
r0_bond = r_eq_fit
print("k_bond_hartree =", k_bond_hartree)
print("k_bond_kcal =", k_bond_kcal)
print("r0_bond =", r0_bond)
k_bond_hartree = 0.5167228014968882 k_bond_kcal = 324.2484533711187 r0_bond = 1.535415893216224
Аппроксимация зависимости энергии этана от длины связи C–C параболой позволила определить равновесную длину связи 𝑟0 = 1.5354 Å и силовую константу 𝑘 =324.25 kcal/(mol·Å²), характеризующую жёсткость ковалентной связи.
In [38]:
comparison_text = f"""
Сравнение с GAFF
1) Связь c3-c3
GAFF:
K_r = 300.9 kcal/(mol·Å^2)
r_eq = 1.5375 Å
Мой результат:
K_r = {k_bond_kcal:.4f} kcal/(mol·Å^2)
r_eq = {r0_bond:.5f} Å
2) Угол c3-c3-hc
GAFF:
K_theta = 46.4 kcal/(mol·rad^2)
theta_eq = 110.05 deg
Мой результат:
K_theta = {k_angle_kcal_rad2:.4f} kcal/(mol·rad^2)
theta_eq = {theta0_deg:.5f} deg
"""
comparison_file = results_dir / "gaff_comparison.txt"
comparison_file.write_text(comparison_text, encoding="utf-8")
print(comparison_text)
Сравнение с GAFF 1) Связь c3-c3 GAFF: K_r = 300.9 kcal/(mol·Å^2) r_eq = 1.5375 Å Мой результат: K_r = 324.2485 kcal/(mol·Å^2) r_eq = 1.53542 Å 2) Угол c3-c3-hc GAFF: K_theta = 46.4 kcal/(mol·rad^2) theta_eq = 110.05 deg Мой результат: K_theta = 631.1614 kcal/(mol·rad^2) theta_eq = 112.45159 deg
Для связи C–C результаты хорошо согласуются с GAFF. Равновесная длина связи практически совпадает, а константа жёсткости отличается умеренно. Для угла H–C–C согласование хуже. Равновесный угол отличается не слишком сильно, но константа получилась намного больше, что, вероятно, связано с жёстким одномерным сканом без переоптимизации остальных координат и с тем, что параметры GAFF являются обобщёнными для широкого класса молекул.
Операции для торсионного угла C–C¶
In [40]:
# Вспомогательная функция, чтобы угол всегда лежал в диапазоне [-180, 180]
def wrap_deg(angle):
return ((angle + 180) % 360) - 180
In [41]:
# Z-matrix для скана по торсионному углу
# Здесь мы поворачиваем всю одну метильную группу относительно другой вокруг связи C-C
def make_ethane_torsion_zmat(phi_deg, rcc=R0_CC, rch=R_CH, ahcc=111.200):
d1 = wrap_deg(180.0 + phi_deg)
d2 = wrap_deg(120.0 + phi_deg)
d3 = wrap_deg(-120.0 + phi_deg)
return f"""
0 1
C
C 1 rCC
H 1 rCH 2 aHCC
H 1 rCH 2 aHCC 3 120.0
H 1 rCH 2 aHCC 3 -120.0
H 2 rCH 1 aHCC 3 d1
H 2 rCH 1 aHCC 5 d2
H 2 rCH 1 aHCC 5 d3
rCC = {rcc:.5f}
rCH = {rch:.5f}
aHCC = {ahcc:.3f}
d1 = {d1:.3f}
d2 = {d2:.3f}
d3 = {d3:.3f}
units angstrom
symmetry c1
"""
In [42]:
# Просто имена файлов
def safe_torsion_label(phi):
sign = "m" if phi < 0 else "p"
return f"ethane_torsion_{sign}{abs(int(phi)):03d}"
In [43]:
# Значения торсионного угла
torsion_values = np.arange(-180, 180 + 12, 12)
torsion_values
Out[43]:
array([-180, -168, -156, -144, -132, -120, -108, -96, -84, -72, -60,
-48, -36, -24, -12, 0, 12, 24, 36, 48, 60, 72,
84, 96, 108, 120, 132, 144, 156, 168, 180])
In [44]:
# Запуск расчётов
torsion_records = []
for phi in torsion_values:
label = safe_torsion_label(phi)
zmat = make_ethane_torsion_zmat(phi)
res = run_psi4(zmat, label)
res["torsion_deg"] = phi
torsion_records.append(res)
df_tors = pd.DataFrame(torsion_records).sort_values("torsion_deg").reset_index(drop=True)
df_tors["dE_hartree"] = df_tors["energy_hartree"] - df_tors["energy_hartree"].min()
df_tors["dE_kcal_mol"] = df_tors["dE_hartree"] * HARTREE_TO_KCAL
torsion_csv = results_dir / "torsion_scan_summary.csv"
df_tors.to_csv(torsion_csv, index=False)
df_tors.head()
Out[44]:
| label | method_basis | out_file | xyz_file | zmat_file | status | energy_hartree | torsion_deg | dE_hartree | dE_kcal_mol | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | ethane_torsion_m180 | scf/cc-pvtz | results/psi4_logs/ethane_torsion_m180.out | results/xyz/ethane_torsion_m180.xyz | results/xyz/ethane_torsion_m180.zmat.txt | ok | -79.156791 | -180 | 0.000758 | 0.475817 |
| 1 | ethane_torsion_m168 | scf/cc-pvtz | results/psi4_logs/ethane_torsion_m168.out | results/xyz/ethane_torsion_m168.xyz | results/xyz/ethane_torsion_m168.zmat.txt | ok | -79.155663 | -168 | 0.001886 | 1.183264 |
| 2 | ethane_torsion_m156 | scf/cc-pvtz | results/psi4_logs/ethane_torsion_m156.out | results/xyz/ethane_torsion_m156.xyz | results/xyz/ethane_torsion_m156.zmat.txt | ok | -79.154552 | -156 | 0.002997 | 1.880561 |
| 3 | ethane_torsion_m144 | scf/cc-pvtz | results/psi4_logs/ethane_torsion_m144.out | results/xyz/ethane_torsion_m144.xyz | results/xyz/ethane_torsion_m144.zmat.txt | ok | -79.153882 | -144 | 0.003667 | 2.301203 |
| 4 | ethane_torsion_m132 | scf/cc-pvtz | results/psi4_logs/ethane_torsion_m132.out | results/xyz/ethane_torsion_m132.xyz | results/xyz/ethane_torsion_m132.zmat.txt | ok | -79.153914 | -132 | 0.003634 | 2.280676 |
Здесь мы рассчитали энергию этана при разных торсионных углах вокруг связи C–C. В отличие от длины связи и валентного угла, здесь получается периодическая зависимость энергии.
In [45]:
# Построение графика торсионного профиля
plt.figure(figsize=(8, 4.8))
plt.plot(df_tors["torsion_deg"], df_tors["dE_kcal_mol"], "o-")
plt.xlabel("Торсионный угол, градусы")
plt.ylabel("Относительная энергия, kcal/mol")
plt.title("Зависимость энергии этана от торсионного угла")
plt.grid(alpha=0.3)
plt.tight_layout()
torsion_plot = plots_dir / "torsion_scan_ethane.png"
plt.savefig(torsion_plot, dpi=250, bbox_inches="tight")
plt.show()
Построили график зависимости энергии от торсионного угла
Нахождение минимумов функций¶
In [46]:
# Убираем дублирование точки 180°, потому что она эквивалентна -180°
df_tors_unique = df_tors[df_tors["torsion_deg"] < 180].copy().reset_index(drop=True)
x = df_tors_unique["torsion_deg"].to_numpy()
y = df_tors_unique["dE_kcal_mol"].to_numpy()
min_indices = []
for i in range(len(y)):
prev_i = (i - 1) % len(y)
next_i = (i + 1) % len(y)
if y[i] <= y[prev_i] and y[i] <= y[next_i]:
min_indices.append(i)
torsion_minima = df_tors_unique.iloc[min_indices][["torsion_deg", "dE_kcal_mol"]]
n_minima = len(torsion_minima)
print("Число уникальных минимумов:", n_minima)
torsion_minima
Число уникальных минимумов: 3
Out[46]:
| torsion_deg | dE_kcal_mol | |
|---|---|---|
| 8 | -84 | 4.155530e-09 |
| 18 | 36 | 1.515966e-10 |
| 28 | 156 | 0.000000e+00 |
In [47]:
plt.figure(figsize=(8, 4.8))
plt.plot(df_tors_unique["torsion_deg"], df_tors_unique["dE_kcal_mol"], "o-", label="Профиль")
plt.scatter(torsion_minima["torsion_deg"], torsion_minima["dE_kcal_mol"], s=80, label="Минимумы")
plt.xlabel("Торсионный угол, градусы")
plt.ylabel("Относительная энергия, kcal/mol")
plt.title("Минимумы торсионного профиля этана")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
torsion_min_plot = plots_dir / "torsion_scan_ethane_minima.png"
plt.savefig(torsion_min_plot, dpi=250, bbox_inches="tight")
plt.show()
In [50]:
torsion_report = results_dir / "torsion_report.txt"
torsion_report.write_text(
f"""Торсионный скан этана
Число уникальных минимумов: {n_minima}
Координаты минимумов:
{torsion_minima.to_string(index=False)}
""",
encoding="utf-8"
)
print("Отчёт сохранён:", torsion_report)
Отчёт сохранён: /Users/romashka/Documents/jupyter/mm_ethane/results/torsion_report.txt
Мы автоматически нашли минимумы функции и сохранили их в отдельный файл. У этана получилось три уникальных минимума за 360°.
Для торсионного угла вокруг связи C–C был выполнен скан в диапазоне от -180° до 180° с шагом 12°. По полученным энергиям построен торсионный профиль этана. Функция имеет периодический характер и содержит три уникальных минимума за полный оборот, что соответствует трём энергетически выгодным staggered-конформациям.
In [ ]: