Задание №1. Зависимость энергии молекулы от длины связи.

In [1]:
import numpy
import scipy.special
import scipy.misc
import psi4
In [2]:
inp = """
C
C 1 cc
H 1 ch 2 hcc
H 1 ch 2 hcc 3 120.0
H 1 ch 2 hcc 3 -120.0
H 2 ch 1 hcc 3 180.0
H 2 ch 1 hcc 5 120.0
H 2 ch 1 hcc 5 -120.0

cc = 1.52986
ch = 1.08439
hcc = 111.200
"""
In [3]:
def run_psi4(inp):
    m = psi4.geometry(inp)
    psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True})
    ener=psi4.energy('scf/cc-pvtz', molecule = m )
    return ener
In [4]:
psi4.core.set_output_file('output_CC.dat')
In [10]:
ener_list=[]
cc_list=[]
zmat_template = """C
C 1 cc
H 1 ch 2 hcc
H 1 ch 2 hcc 3 120.0
H 1 ch 2 hcc 3 -120.0
H 2 ch 1 hcc 3 180.0
H 2 ch 1 hcc 5 120.0
H 2 ch 1 hcc 5 -120.0

ch = 1.08439
hcc = 111.200
cc = {cc:.5f}
"""

# Начальное значение cc и шаг здесь мы взяли стартом 1.32, это значение отличается от тех задания, чтобы показать параболическую зависимость энергии от длины связи.
cc_start = 1.32986
step = 0.02

# Генерируем 20 вариантов (i = 0..19)
for i in range(21):
    cc_value = cc_start + i * step
    cc_list.append(cc_value)
    # Формируем текст Z-матрицы с новым значением cc
    zmat_text = zmat_template.format(cc=cc_value)
    ener=run_psi4(zmat_text)
    ener_list.append(ener)
In [11]:
cc_list
Out[11]:
[1.32986,
 1.34986,
 1.36986,
 1.38986,
 1.4098600000000001,
 1.4298600000000001,
 1.4498600000000001,
 1.4698600000000002,
 1.48986,
 1.50986,
 1.52986,
 1.54986,
 1.56986,
 1.58986,
 1.60986,
 1.62986,
 1.64986,
 1.6698600000000001,
 1.68986,
 1.70986,
 1.72986]
In [12]:
ener_list
Out[12]:
[-79.12147106144542,
 -79.12862588526878,
 -79.13470987830091,
 -79.13981745732596,
 -79.1440349857867,
 -79.14744150815146,
 -79.15010940862017,
 -79.15210500168824,
 -79.15348906310734,
 -79.15431730770699,
 -79.15464082062108,
 -79.1545064471334,
 -79.15395714637484,
 -79.15303231272178,
 -79.1517680688072,
 -79.1501975331354,
 -79.14835106522409,
 -79.14625649026685,
 -79.14393930553118,
 -79.14142287022013,
 -79.13872858014733]
In [13]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [14]:
x_o=cc_list
y_o=ener_list

#function is  f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,1, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Optimized params:", p1)

#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlim(1,2)
plt.show()

Optimized params: [  0.61075072   1.5541475  -79.15517139]
No description has been provided for this image

Аппроксимация функции оставляет желать лучшего, отсюда и небольшая ошибка, на 0.03 в большую сторону по сравнению с референсом из статьи по GAFF2.

Задание №2 зависимость энергии молекулы от величины валентного угла HCC.

In [10]:
psi4.core.set_output_file('output_HCC.dat')
In [27]:
ener_list=[]
hcc_list=[]
zmat_template = """C
C 1 cc
H 1 ch 2 hcc
H 1 ch 2 hcc 3 120.0
H 1 ch 2 hcc 3 -120.0
H 2 ch 1 hcc 3 180.0
H 2 ch 1 hcc 5 120.0
H 2 ch 1 hcc 5 -120.0

ch = 1.08439
hcc = {hcc:.5f}
cc = 1.52986
"""

# Начальное значение cc и шаг
hcc_start = 109.2
step = 0.2

# Генерируем 20 вариантов (i = 0..19)
for i in range(21):
    hcc_value = hcc_start + i * step
    hcc_list.append(hcc_value)
    # Формируем текст Z-матрицы с новым значением cc
    zmat_text = zmat_template.format(hcc=hcc_value)
    ener=run_psi4(zmat_text)
    ener_list.append(ener)
In [28]:
hcc_list
Out[28]:
[109.2,
 109.4,
 109.60000000000001,
 109.8,
 110.0,
 110.2,
 110.4,
 110.60000000000001,
 110.8,
 111.0,
 111.2,
 111.4,
 111.60000000000001,
 111.8,
 112.0,
 112.2,
 112.4,
 112.60000000000001,
 112.8,
 113.0,
 113.2]
In [29]:
ener_list
Out[29]:
[-79.15188374120805,
 -79.15226860400742,
 -79.15262930725052,
 -79.15296581594981,
 -79.15327809447629,
 -79.15356610613043,
 -79.15382981348127,
 -79.15406917825445,
 -79.15428416118151,
 -79.15447472228726,
 -79.1546408206485,
 -79.15478241444777,
 -79.15489946116351,
 -79.15499191724419,
 -79.15505973829067,
 -79.15510287913138,
 -79.15512129348916,
 -79.15511493448462,
 -79.15508375409306,
 -79.15502770340468,
 -79.15494673277927]
In [31]:
x_o=hcc_list
y_o=ener_list

#function is  f(x)=k(b-x)^2 + a
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function

p0 = [1,1, -79] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Optimized params:", p1)

#Plot it
plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlim(108,114)
plt.show()

Optimized params: [ 3.06390208e-04  1.12451592e+02 -7.91551208e+01]
No description has been provided for this image
In [ ]:
Задание №3 Зависимость энергии молекулы от торсионного угла С-С.
In [17]:
psi4.core.set_output_file('output_torsion_CC.dat')
In [18]:
ener_list = []
tor_list = []

zmat_template = """C
C 1 cc
H 1 ch 2 hcc
H 1 ch 2 hcc 3 120.0
H 1 ch 2 hcc 3 -120.0
H 2 ch 1 hcc 3 {t1}
H 2 ch 1 hcc 3 {t2}
H 2 ch 1 hcc 3 {t3}

cc = 1.52986
ch = 1.08439
hcc = 111.200
"""

tor_start = -180
step = 6

for i in range(61):
    tor_value = tor_start + i * step
    tor_list.append(tor_value)
    t1 = tor_value
    t2 = tor_value + 120.0
    t3 = tor_value - 120.0
    zmat_text = zmat_template.format(t1=t1, t2=t2, t3=t3)
    ener=run_psi4(zmat_text)
    ener_list.append(ener)
In [19]:
tor_list
Out[19]:
[-180,
 -174,
 -168,
 -162,
 -156,
 -150,
 -144,
 -138,
 -132,
 -126,
 -120,
 -114,
 -108,
 -102,
 -96,
 -90,
 -84,
 -78,
 -72,
 -66,
 -60,
 -54,
 -48,
 -42,
 -36,
 -30,
 -24,
 -18,
 -12,
 -6,
 0,
 6,
 12,
 18,
 24,
 30,
 36,
 42,
 48,
 54,
 60,
 66,
 72,
 78,
 84,
 90,
 96,
 102,
 108,
 114,
 120,
 126,
 132,
 138,
 144,
 150,
 156,
 162,
 168,
 174,
 180]
In [51]:
ener_list
Out[51]:
[-79.26001266654873,
 -79.25989321402565,
 -79.2595460653074,
 -79.25900385168322,
 -79.25831799035231,
 -79.25755427824926,
 -79.25678696770814,
 -79.2560916895934,
 -79.25553784349881,
 -79.25518128802292,
 -79.25505826169737,
 -79.25518128801822,
 -79.25553784348088,
 -79.25609168958417,
 -79.25678696769788,
 -79.25755427823658,
 -79.2583179903536,
 -79.2590038516943,
 -79.2595460653112,
 -79.25989321403887,
 -79.26001266652985,
 -79.25989321401934,
 -79.25954606531306,
 -79.2590038516832,
 -79.25831799035998,
 -79.257554278248,
 -79.25678696768703,
 -79.25609168959686,
 -79.25553784346853,
 -79.25518128800182,
 -79.25505826170816,
 -79.25518128802287,
 -79.25553784347883,
 -79.25609168961161,
 -79.2567869677126,
 -79.25755427824416,
 -79.25831799033709,
 -79.25900385167328,
 -79.25954606533094,
 -79.2598932140219,
 -79.2600126665324,
 -79.25989321403509,
 -79.25954606532528,
 -79.25900385170411,
 -79.25831799034388,
 -79.25755427825439,
 -79.2567869676932,
 -79.25609168959164,
 -79.25553784348762,
 -79.25518128802449,
 -79.25505826171215,
 -79.25518128801319,
 -79.25553784346846,
 -79.256091689632,
 -79.25678696771524,
 -79.25755427825327,
 -79.25831799034069,
 -79.25900385169271,
 -79.25954606530128,
 -79.25989321403172,
 -79.26001266654788]
In [25]:
x_o=tor_list
y_o=ener_list

fitfunc = lambda p, x: p[0] * (1 + np.cos(p[1] * np.pi * np.radians(x) - p[2])) + p[3]  # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y  # Error function

p0 = [0, 1, -1, -1] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o, y_o))
print("Optimized params:", p1)

#Plot it
plt.plot(x_o, fitfunc(p1, x_o), c='red')
plt.plot(x_o, y_o)
plt.xlim(-180,180)
plt.show()

Optimized params: [ 2.47680387e-03  9.55027970e-01  7.53376390e-06 -7.92600213e+01]
No description has been provided for this image

У функции зависимости энергии молекулы от торсионного угла наблюдается 3 минимума (минимум при -180 и 180 это один и тот же минимум).