!pip install -q condacolab
import condacolab
condacolab.install()
condacolab.check()
✨🍰✨ Everything looks OK! ✨🍰✨ Everything looks OK!
%%capture
! mamba install -c anaconda intel-openmp
! mamba install -c psi4 psi4
! mamba install -c conda-forge psiresp
!pip install py3Dmol
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import psi4
import numpy as np
psi4.core.set_output_file('output1.dat')
#real-like ethan
import py3Dmol
m = psi4.geometry("""
C
C 1 1.52986
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 5 60
H 2 1.08439 1 111.200 5 180
H 2 1.08439 1 111.200 5 -60
""")
view = py3Dmol.view()
view.addModel(m.save_string_xyz_file(),'xyz')
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
<py3Dmol.view at 0x7f4328d742e0>
#bound length
psi4.set_memory('500 MB')
psi4.core.set_output_file('output1.dat')
f = open('psi4-1.out','w')
f.write('length\tenergy\n')
bond_E = []
for length in np.round(np.arange(1.25,1.75,0.02),2):
inp = psi4.geometry(f"""
C
C 1 {length}
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 5 60
H 2 1.08439 1 111.200 5 180
H 2 1.08439 1 111.200 5 -60
""")
energy = psi4.energy('scf/cc-pvdz')
f.write(str(length) + '\t' + str(energy) + '\n')
bond_E.append(energy)
f.close()
x_o=np.arange(1.2,1.7,0.02)
y_o=bond_E
#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, -70] # 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.1,1.8)
plt.show()
Optimized params: [ 0.76062806 1.4992685 -79.23695634]
<ipython-input-6-941ef5b07d0f>:37: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5) <ipython-input-6-941ef5b07d0f>:37: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
Вышло не очень. Впрочем, мы же помним, что фунция там чуть другая - глянем дальше.)
#angle
psi4.set_memory('500 MB')
psi4.core.set_output_file('output2.dat')
f = open('psi4-2.out','w')
f.write('angle\tenergy\n')
angle_E = []
for angle in np.round(np.arange(109.2,113.4,0.1),2):
inp = psi4.geometry(f"""
C
C 1 1.52986
H 1 1.08439 2 {angle}
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 5 60
H 2 1.08439 1 111.200 5 180
H 2 1.08439 1 111.200 5 -60
""")
energy = psi4.energy('scf/cc-pvdz')
f.write(str(angle) + '\t' + str(energy) + '\n')
angle_E.append(energy)
f.close()
x_o=np.arange(109.2,113.4,0.1)
y_o=angle_E
#function is f(x)=ax^2 + bx + c
fitfunc = lambda p, x: p[0]*x**2 + p[1]*x + p[2] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1, 1, -70] # 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(109,113.6)
plt.show()
Optimized params: [ 3.81547897e-05 -8.48091748e-03 -7.87633605e+01]
<ipython-input-7-34296aab55be>:37: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5) <ipython-input-7-34296aab55be>:37: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
Неплохо вышло. Все же квадратный многочлен... Оптимальный угол - 112° против 109.5° для sp3 углерода из статьи о GAFF. Вероятно, это следствие меньшего размера атома водорода - отталкивание меньше, угол больше.
#torsion
psi4.set_memory('500 MB')
psi4.core.set_output_file('output3.dat')
f = open('psi4-3.out','w')
f.write('torsion\tenergy\n')
torsion_E = []
for torsion in np.round(np.arange(-180,181,12),1):
inp = psi4.geometry(f"""
C
C 1 1.52986
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 5 {torsion}
H 2 1.08439 1 111.200 5 {torsion+120}
H 2 1.08439 1 111.200 5 {torsion-120}
""")
energy = psi4.energy('scf/cc-pvdz')
f.write(str(torsion) + '\t' + str(energy) + '\n')
torsion_E.append(energy)
f.close()
x_o=np.arange(-180,181,12)
y_o=torsion_E
#function is f(x)=k * cos(ax + b) + c
fitfunc = lambda p, x: p[0]*np.cos(np.radians(p[1] * np.pi * x - p[2]* np.pi)) + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1, 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, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
plt.xlim(-185,185)
plt.show()
Optimized params: [ 2.66937519e-03 9.55015104e-01 3.31413461e-04 -7.92319747e+01]
<ipython-input-8-119318b04d1d>:36: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5) <ipython-input-8-119318b04d1d>:36: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
Получилось три максимума - как и ожидалось. Фиттирование прошло довольно успешно.
#bound length high scale
psi4.set_memory('500 MB')
psi4.core.set_output_file('output4.dat')
f = open('psi4-4.out','w')
f.write('length\tenergy\n')
range = np.arange(0.5,5,0.1)
bond_E = []
for length in np.round(range,2):
inp = psi4.geometry(f"""
C
C 1 {length}
H 1 1.08439 2 111.200
H 1 1.08439 2 111.200 3 120
H 1 1.08439 2 111.200 3 -120
H 2 1.08439 1 111.200 5 60
H 2 1.08439 1 111.200 5 180
H 2 1.08439 1 111.200 5 -60
""")
energy = psi4.energy('scf/cc-pvdz')
f.write(str(length) + '\t' + str(energy) + '\n')
bond_E.append(energy)
f.close()
x_o=range
y_o=bond_E
#function is f(x)=k(1-e^(-a(r-re))^2 + b - morze potential
fitfunc = lambda p, x: p[0] * (1 - np.exp(- p[1] * (x - p[2]))) ** 2 + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1, 1, 1, -70] # 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(0,5.5)
plt.show()
/usr/local/lib/python3.9/site-packages/scipy/optimize/minpack.py:475: RuntimeWarning: Number of calls to function has reached maxfev = 1000. warnings.warn(errors[info][0], RuntimeWarning) <ipython-input-9-931348a190ef>:39: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5) <ipython-input-9-931348a190ef>:39: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
Optimized params: [-1.12279437e+02 -1.66081076e-03 -2.66717473e+02 -4.29604023e+01]
Плоховато...
x_o=range
y_o=bond_E
#function is f(x)=k(1-e^(-a(r-re))^2 + b - morze potential
fitfunc = lambda p, x: p[0] * (1 - np.exp(- p[1] * (x - p[2]))) ** 2 + p[3] # Target function
errfunc = lambda p, x, y: fitfunc(p, x) - y # Error function
p0 = [1, 2, 1, -70] # 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(0,5.5)
plt.show()
Optimized params: [ 0.30617538 1.86664453 1.41970257 -79.2589307 ]
<ipython-input-10-3bfd762f4c2e>:13: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "ro" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5) <ipython-input-10-3bfd762f4c2e>:13: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "r-" (-> color='r'). The keyword argument will take precedence. plt.plot(x_o, y_o, "ro", x_o,fitfunc(p1,x_o),"r-",c='blue',alpha=0.5)
А подбор начальных значений помог, забавно. Локальный минимум)
Идеальный радиус связи близок к литературному - 1.4342 против 1.526 из статьи о GAFF. Разница, вероятно, вызвана некоторыми приближениями, возможно, а еще и влиянием атомов водорода.
%%capture
!mamba install -c conda-forge psiresp
import psiresp
import qcfractal.interface as ptl
# set up molecule
ethan = psiresp.Molecule.from_smiles("CC")
from psiresp.testing import FractalSnowflake
server = FractalSnowflake()
client = ptl.FractalClient(server, verify=False)
/usr/local/lib/python3.9/site-packages/psiresp/charge.py:282: FutureWarning: `symmetric_atoms_are_equivalent` will be set to False by default for now, as it is a new feature. It will be set to True by default in the future warnings.warn(
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-12-fa77f3ca7ff3> in <module> 6 7 from psiresp.testing import FractalSnowflake ----> 8 server = FractalSnowflake() 9 client = ptl.FractalClient(server, verify=False) 10 /usr/local/lib/python3.9/site-packages/psiresp/testing.py in __init__(self, max_workers, storage_project_name, max_active_services, logging, start_server, reset_database) 109 start_server: bool = True, 110 reset_database: bool = False,): --> 111 storage = TemporaryPostgres(database_name=storage_project_name) 112 storage_uri = storage.database_uri(safe=False, database="") 113 super().__init__(max_workers=max_workers, /usr/local/lib/python3.9/site-packages/psiresp/testing.py in __init__(self, database_name, tmpdir, quiet, logger) 94 self.config = FractalConfig(database=config_data) 95 self.psql = PostgresHarness(self.config) ---> 96 self.psql.initialize_postgres() 97 self.psql.init_database() 98 # self.psql.start() /usr/local/lib/python3.9/site-packages/qcfractal/postgres_harness.py in initialize_postgres(self) 332 init_status = self._run([shutil.which("initdb"), "-D", self.config.database_path]) 333 if "Success." not in init_status["stdout"]: --> 334 raise ValueError(f"Could not initialize the PostgreSQL server. Error below:\n\n{init_status['stderr']}") 335 336 # Change any configurations ValueError: Could not initialize the PostgreSQL server. Error below: initdb: error: cannot be run as root initdb: hint: Please log in (using, e.g., "su") as the (unprivileged) user that will own the server process.
Ну ладно, нет так нет...
from IPython.core.display import display, HTML
html_custom = '<h4>Полученные выше файлы лежат в этой <a href="https://kodomo.fbb.msu.ru/~maxim2002/term8/pr6_files">папке</a>.</h4>'
display(HTML(html_custom))