Software
Contents
17. Software#
Los software de estructura electrónica se dedican a calcular diversas propiedades de las moléculas utilizando la teoría que se ve en química cuántica. En general, esto nos sirve para predecir la reactividad química, desde poder predecir si una reacción procederá o no, hasta cosas más avanzadas como cinéticas de reacción y pKas o potenciales redox. En este notebook estaremos usando el software PySCF, sin embargo, existe una gran cantidad de software, desde los libres y gratuitos hasta los privados y de pago.
Note
Algunos programas de estructura electrónica son:
Gaussian (https://gaussian.com)
PySCF (https://pyscf.org)
Psi4 (http://www.psicode.org)
NWChem (http://www.nwchem-sw.org)
QChem (http://www.q-chem.com)
TeraChem (http://www.petachem.com)
deMon2k (http://www.demon-software.com)
Orca (https://orcaforum.cec.mpg.de)
Molcas (http://www.molcas.org)
ADF (https://www.scm.com)
Quantum Espresso (https://www.quantum-espresso.org)
17.1. Aprendiendo a usar PySCF#
Warning
Si está utilizando Google Colab o la ejecución en línea, debe de ejecutar al inicio el siguiente código
!pip install pyscf
!pip install geomeTRIC
Para usar PySCF, puede importarlo como si de una librería se tratase, es decir
import pyscf
Importe PySCF en la siguiente celda
# importe pyscf
# Descomentar estas líneas si está en modo online
#!pip install pyscf
#!pip install geomeTRIC
import pyscf
El siguiente paso es declarar las coordenadas de los átomos que forman la molécula. Para ello se pueden usar visualizadores como Avogadro o IQmol. También es posible obtener valores experimentales o calculados de https://cccbdb.nist.gov/ . En este caso utilizaremos los resultados experimentales de benceno.
Use las siguientes líneas para declarar la geometría y seleccionar una base
benzene = pyscf.gto.Mole(atom = """
C 0.0000 1.3970 0.0000
C 1.2098 0.6985 0.0000
C 1.2098 -0.6985 0.0000
C 0.0000 -1.3970 0.0000
C -1.2098 -0.6985 0.0000
C -1.2098 0.6985 0.0000
H 0.0000 2.4810 0.0000
H 2.1486 1.2405 0.0000
H 2.1486 -1.2405 0.0000
H 0.0000 -2.4810 0.0000
H -2.1486 -1.2405 0.0000
H -2.1486 1.2405 0.0000
""",basis = "6-31G")
benzene = benzene.build()
#Geometría
benzene = pyscf.gto.Mole(atom = """
C 0.0000 1.3970 0.0000
C 1.2098 0.6985 0.0000
C 1.2098 -0.6985 0.0000
C 0.0000 -1.3970 0.0000
C -1.2098 -0.6985 0.0000
C -1.2098 0.6985 0.0000
H 0.0000 2.4810 0.0000
H 2.1486 1.2405 0.0000
H 2.1486 -1.2405 0.0000
H 0.0000 -2.4810 0.0000
H -2.1486 -1.2405 0.0000
H -2.1486 1.2405 0.0000
""",basis = "6-31G")
benzene = benzene.build()
En la mayoría de los software es común (pero no obligatorio) que antes de mandar el cálculo de una molécula se asigne una cantidad de memoria RAM, por ejemplo 2000 MB.
En PySCF esto se hace mediante la instrucción
mol.max_memory = 2000
Cada vez que modifique algún parámetro de la molécula, debe volver a construirla antes de usarla.
Asigne memoria a su cálculo
# Establezca memoria (y reconstruya la molécula)
# Establezca memoria (en MB) (y reconstruya la molécula)
benzene.max_memory = 2000
benzene = benzene.build()
Note
Los software usan la RAM asignada para guardar vectores y matrices como lo ha hecho en las prácticas anteriores. Si la memoria es suficiente, el programa guardará todo y el cálculo será más rápido, si no la hay, el cálculo será más lento. A más memoria los cálculos tienden a ser igual o más rápidos
Warning
La cantidad de memoria que puede asignar al cálculo depende de la cantidad de RAM que tenga su computadora. Recomendamos asignar menos memoria del total disponible ya que la memoria se reparte con los demás programas de su computadora.
Para realizar un cálculo de energía de una molécula con la geometría y la base especificada arriba, es necesario especificar un método. Por ejemplo, para hacer Hartree-Fock, escribiremos
from pyscf improt scf
rhf = scf.RHF(mol)
E = rhf.kernel()
hemos de cambiar mol por el nombre que le asignamos a nuestra molécula.
Realice un cálculo con el método HF y la base 6-31G
# Benceno HF
# Benceno HF
from pyscf import scf
rhf = scf.RHF(benzene)
E = rhf.kernel()
converged SCF energy = -230.623508002015
Calcule la energía de benceno con MP2 y 6-31G
from pyscf import mp
mp2 = mp.MP2(mol)
Ecorr = mp2.kernel()
# Benceno MP2/6-31G
# Benceno MP2/6-31G
from pyscf import mp
mp2 = mp.MP2(benzene)
Ecorr = mp2.kernel()
converged SCF energy = -230.623508002015
E(MP2) = -231.147363985909 E_corr = -0.52385598389374
Calcule la energía de benceno con el funcional M062X y la base 6-31G
Use el código
from pyscf import dft
rks = dft.RKS(mol)
rks.xc = 'M062X'
E = rks.kernel()
# Benceno M062X/6-31G
# Benceno M062X/6-31G
from pyscf import dft
rks = dft.RKS(benzene)
rks.xc = 'M062X'
E = rks.kernel()
converged SCF energy = -232.091626867516
Usualmente la geometría especificada no es necesariamente la geometría real. Es posible pedir al software que mueva los átomos hasta encontrar las coordenadas que representen un mínimo de energía con el método y base usados. Esto se llama optimización de geometrías, para lo cual hay que importar el optimizador y pasar la instancia del método. Por ejemplo, para optimizar con M062X, utilizaremos el rks que definimos antes
from pyscf.geomopt.geometric_solver import optimize
mol = geometric_opt(rks)
Optimice la geometría de benceno con el método M062X y base 6-31G e imprima su energía.
Warning
Este cálculo puede tardar entre segundos y unos minutos dependiendo del procesador de cada computadora.
# Optimizacion de Geometría
from pyscf.geomopt.geometric_solver import optimize
mol = optimize(rhf)
geometric-optimize called with the following command line:
/home/jfhlewyee/anaconda3/envs/jb/lib/python3.10/site-packages/ipykernel_launcher.py -f /tmp/tmpys1usha8.json --HistoryManager.hist_file=:memory:
())))))))))))))))/
())))))))))))))))))))))))),
*)))))))))))))))))))))))))))))))))
#, ()))))))))/ .)))))))))),
#%%%%, ()))))) .))))))))*
*%%%%%%, )) .. ,))))))).
*%%%%%%, ***************/. .)))))))
#%%/ (%%%%%%, /*********************. )))))))
.%%%%%%# *%%%%%%, *******/, **********, .))))))
.%%%%%%/ *%%%%%%, ** ******** .))))))
## .%%%%%%/ (%%%%%%, ,****** /)))))
%%%%%% .%%%%%%# *%%%%%%, ,/////. ****** ))))))
#% %% .%%%%%%/ *%%%%%%, ////////, *****/ ,)))))
#%% %%% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
#%%%%. %%%%%# /%%%%%%* #%%%%%% /////) ****** ))))),
#%%%%##% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
## %%% .%%%%%%/ *%%%%%%, ////////. *****/ ,)))))
#%%%%# /%%%%%%/ (%%%%%% /)/)// ****** ))))))
## .%%%%%%/ (%%%%%%, ******* ))))))
.%%%%%%/ *%%%%%%, **. /******* .))))))
*%%%%%%/ (%%%%%% ********/*..,*/********* *))))))
#%%/ (%%%%%%, *********************/ )))))))
*%%%%%%, ,**************/ ,))))))/
(%%%%%% () ))))))))
#%%%%, ()))))) ,)))))))),
#, ()))))))))) ,)))))))))).
()))))))))))))))))))))))))))))))/
())))))))))))))))))))))))).
())))))))))))))),
-=# geomeTRIC started. Version: 1.0 #=-
Current date and time: 2022-12-13 13:54:32
Custom engine selected.
Bonds will be generated from interatomic distances less than 1.20 times sum of covalent radii
36 internal coordinates being used (instead of 36 Cartesians)
Internal coordinate system (atoms numbered from 1):
Distance 1-2
Distance 1-6
Distance 1-7
Distance 2-3
Distance 2-8
Distance 3-4
Distance 3-9
Distance 4-5
Distance 4-10
Distance 5-6
Distance 5-11
Distance 6-12
Angle 2-1-7
Angle 6-1-7
Angle 1-2-8
Angle 3-2-8
Angle 2-3-9
Angle 4-3-9
Angle 3-4-10
Angle 5-4-10
Angle 4-5-11
Angle 6-5-11
Angle 1-6-12
Angle 5-6-12
Out-of-Plane 1-2-6-7
Out-of-Plane 2-1-3-8
Out-of-Plane 3-2-4-9
Out-of-Plane 4-3-5-10
Out-of-Plane 5-4-6-11
Out-of-Plane 6-1-5-12
Dihedral 6-1-2-3
Dihedral 6-1-2-8
Dihedral 7-1-2-3
Dihedral 7-1-2-8
Dihedral 2-1-6-5
Dihedral 2-1-6-12
Dihedral 7-1-6-5
Dihedral 7-1-6-12
Dihedral 1-2-3-4
Dihedral 1-2-3-9
Dihedral 8-2-3-4
Dihedral 8-2-3-9
Dihedral 2-3-4-5
Dihedral 2-3-4-10
Dihedral 9-3-4-5
Dihedral 9-3-4-10
Dihedral 3-4-5-6
Dihedral 3-4-5-11
Dihedral 10-4-5-6
Dihedral 10-4-5-11
Dihedral 4-5-6-1
Dihedral 4-5-6-12
Dihedral 11-5-6-1
Dihedral 11-5-6-12
Translation-X 1-12
Translation-Y 1-12
Translation-Z 1-12
Rotation-A 1-12
Rotation-B 1-12
Rotation-C 1-12
<class 'geometric.internal.Distance'> : 12
<class 'geometric.internal.Angle'> : 12
<class 'geometric.internal.OutOfPlane'> : 6
<class 'geometric.internal.Dihedral'> : 24
<class 'geometric.internal.TranslationX'> : 1
<class 'geometric.internal.TranslationY'> : 1
<class 'geometric.internal.TranslationZ'> : 1
<class 'geometric.internal.RotationA'> : 1
<class 'geometric.internal.RotationB'> : 1
<class 'geometric.internal.RotationC'> : 1
> ===== Optimization Info: ====
> Job type: Energy minimization
> Maximum number of optimization cycles: 300
> Initial / maximum trust radius (Angstrom): 0.100 / 0.300
> Convergence Criteria:
> Will converge when all 5 criteria are reached:
> |Delta-E| < 1.00e-06
> RMS-Grad < 3.00e-04
> Max-Grad < 4.50e-04
> RMS-Disp < 1.20e-03
> Max-Disp < 1.80e-03
> === End Optimization Info ===
Geometry optimization cycle 1
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
C 0.000000 1.397000 0.000000 0.000000 0.000000 0.000000
C 1.209800 0.698500 0.000000 0.000000 0.000000 0.000000
C 1.209800 -0.698500 0.000000 0.000000 0.000000 0.000000
C 0.000000 -1.397000 0.000000 0.000000 0.000000 0.000000
C -1.209800 -0.698500 0.000000 0.000000 0.000000 0.000000
C -1.209800 0.698500 0.000000 0.000000 0.000000 0.000000
H 0.000000 2.481000 0.000000 0.000000 0.000000 0.000000
H 2.148600 1.240500 0.000000 0.000000 0.000000 0.000000
H 2.148600 -1.240500 0.000000 0.000000 0.000000 0.000000
H 0.000000 -2.481000 0.000000 0.000000 0.000000 0.000000
H -2.148600 -1.240500 0.000000 0.000000 0.000000 0.000000
H -2.148600 1.240500 0.000000 0.000000 0.000000 0.000000
converged SCF energy = -230.623508002017
--------------- SCF_Scanner gradients ---------------
x y z
0 C 0.0000000000 0.0015348049 0.0000000000
1 C 0.0013014809 0.0007803543 0.0000000000
2 C 0.0013014809 -0.0007803543 -0.0000000000
3 C 0.0000000000 -0.0015348049 0.0000000000
4 C -0.0013014809 -0.0007803543 -0.0000000000
5 C -0.0013014809 0.0007803543 0.0000000000
6 H 0.0000000000 0.0080211277 -0.0000000000
7 H 0.0069607311 0.0040195330 -0.0000000000
8 H 0.0069607311 -0.0040195330 0.0000000000
9 H 0.0000000000 -0.0080211277 0.0000000000
10 H -0.0069607311 -0.0040195330 0.0000000000
11 H -0.0069607311 0.0040195330 -0.0000000000
----------------------------------------------
cycle 1: E = -230.623508002 dE = -230.624 norm(grad) = 0.0200258
Step 0 : Gradient = 5.781e-03/8.038e-03 (rms/max) Energy = -230.6235080020
Hessian Eigenvalues: 2.30000e-02 2.30000e-02 2.30000e-02 ... 4.61250e-01 4.61268e-01 4.61268e-01
Geometry optimization cycle 2
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
C -0.000000 1.385947 -0.000000 -0.000000 -0.011053 -0.000000
C 1.200366 0.693005 0.000000 -0.009434 -0.005495 0.000000
C 1.200366 -0.693005 -0.000000 -0.009434 0.005495 -0.000000
C -0.000000 -1.385947 -0.000000 -0.000000 0.011053 -0.000000
C -1.200366 -0.693005 -0.000000 0.009434 0.005495 -0.000000
C -1.200366 0.693005 -0.000000 0.009434 -0.005495 -0.000000
H -0.000000 2.457996 0.000000 -0.000000 -0.023004 0.000000
H 2.128777 1.229045 0.000000 -0.019823 -0.011455 0.000000
H 2.128777 -1.229045 -0.000000 -0.019823 0.011455 -0.000000
H -0.000000 -2.457996 -0.000000 -0.000000 0.023004 -0.000000
H -2.128777 -1.229045 -0.000000 0.019823 0.011455 -0.000000
H -2.128777 1.229045 -0.000000 0.019823 -0.011455 -0.000000
WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.
converged SCF energy = -230.624434899035
--------------- SCF_Scanner gradients ---------------
x y z
0 C -0.0000000000 -0.0016273417 -0.0000000000
1 C -0.0013627097 -0.0007867159 0.0000000000
2 C -0.0013627097 0.0007867159 0.0000000000
3 C 0.0000000000 0.0016273417 -0.0000000000
4 C 0.0013627097 0.0007867159 0.0000000000
5 C 0.0013627097 -0.0007867159 -0.0000000000
6 H 0.0000000000 -0.0009777321 0.0000000000
7 H -0.0008418080 -0.0004912628 0.0000000000
8 H -0.0008418080 0.0004912628 0.0000000000
9 H 0.0000000000 0.0009777321 0.0000000000
10 H 0.0008418080 0.0004912628 0.0000000000
11 H 0.0008418080 -0.0004912628 0.0000000000
----------------------------------------------
cycle 2: E = -230.624434899 dE = -0.000926897 norm(grad) = 0.00457296
Step 1 : Displace = 1.797e-02/2.300e-02 (rms/max) Trust = 1.000e-01 (=) Grad = 1.320e-03/1.627e-03 (rms/max) E (change) = -230.6244348990 (-9.269e-04) Quality = 0.814
Hessian Eigenvalues: 2.30000e-02 2.30000e-02 2.30000e-02 ... 4.61250e-01 4.61268e-01 5.51445e-01
Geometry optimization cycle 3
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
C -0.000000 1.388631 0.000000 -0.000000 0.002683 0.000000
C 1.202301 0.694220 -0.000000 0.001935 0.001215 -0.000000
C 1.202301 -0.694220 0.000000 0.001935 -0.001215 0.000000
C -0.000000 -1.388631 0.000000 -0.000000 -0.002683 0.000000
C -1.202301 -0.694220 -0.000000 -0.001935 -0.001215 0.000000
C -1.202301 0.694220 -0.000000 -0.001935 0.001215 0.000000
H -0.000000 2.461850 -0.000000 -0.000000 0.003854 -0.000000
H 2.131768 1.230763 -0.000000 0.002991 0.001718 -0.000000
H 2.131768 -1.230763 -0.000000 0.002991 -0.001718 -0.000000
H -0.000000 -2.461850 0.000000 -0.000000 -0.003854 0.000000
H -2.131768 -1.230763 -0.000000 -0.002991 -0.001718 -0.000000
H -2.131768 1.230763 -0.000000 -0.002991 0.001718 -0.000000
converged SCF energy = -230.624474764293
--------------- SCF_Scanner gradients ---------------
x y z
0 C 0.0000000000 0.0002963729 0.0000000000
1 C 0.0001293122 0.0000674271 -0.0000000000
2 C 0.0001293122 -0.0000674271 -0.0000000000
3 C 0.0000000000 -0.0002963729 0.0000000000
4 C -0.0001293122 -0.0000674271 0.0000000000
5 C -0.0001293122 0.0000674271 -0.0000000000
6 H 0.0000000000 -0.0000415785 -0.0000000000
7 H -0.0000544982 -0.0000200601 0.0000000000
8 H -0.0000544982 0.0000200601 -0.0000000000
9 H 0.0000000000 0.0000415785 -0.0000000000
10 H 0.0000544982 0.0000200601 -0.0000000000
11 H 0.0000544982 -0.0000200601 0.0000000000
----------------------------------------------
cycle 3: E = -230.624474764 dE = -3.98653e-05 norm(grad) = 0.000526966
Step 2 : Displace = 3.063e-03/3.854e-03 (rms/max) Trust = 1.414e-01 (+) Grad = 1.521e-04/2.964e-04 (rms/max) E (change) = -230.6244747643 (-3.987e-05) Quality = 0.957
Hessian Eigenvalues: 2.30000e-02 2.30000e-02 2.30000e-02 ... 4.61250e-01 4.61268e-01 5.83445e-01
Geometry optimization cycle 4
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
C -0.000000 1.387632 -0.000000 -0.000000 -0.000998 -0.000000
C 1.202712 0.694143 0.000000 0.000412 -0.000077 0.000000
C 1.202712 -0.694143 0.000000 0.000412 0.000077 0.000000
C -0.000000 -1.387632 0.000000 -0.000000 0.000998 0.000000
C -1.202712 -0.694143 -0.000000 -0.000412 0.000077 -0.000000
C -1.202712 0.694143 -0.000000 -0.000412 -0.000077 -0.000000
H -0.000000 2.460934 -0.000000 -0.000000 -0.000916 -0.000000
H 2.132109 1.231018 0.000000 0.000340 0.000255 0.000000
H 2.132109 -1.231018 0.000000 0.000340 -0.000255 0.000000
H -0.000000 -2.460934 0.000000 -0.000000 0.000916 0.000000
H -2.132109 -1.231018 -0.000000 -0.000340 -0.000255 -0.000000
H -2.132109 1.231018 -0.000000 -0.000340 0.000255 -0.000000
converged SCF energy = -230.624474239419
--------------- SCF_Scanner gradients ---------------
x y z
0 C 0.0000000000 -0.0003543213 0.0000000000
1 C 0.0001292611 0.0001025532 -0.0000000000
2 C 0.0001292611 -0.0001025532 0.0000000000
3 C -0.0000000000 0.0003543213 0.0000000000
4 C -0.0001292611 -0.0001025532 -0.0000000000
5 C -0.0001292611 0.0001025532 0.0000000000
6 H -0.0000000000 -0.0000153288 -0.0000000000
7 H 0.0000502132 -0.0000106425 0.0000000000
8 H 0.0000502132 0.0000106425 -0.0000000000
9 H 0.0000000000 0.0000153288 -0.0000000000
10 H -0.0000502132 0.0000106425 0.0000000000
11 H -0.0000502132 -0.0000106425 0.0000000000
----------------------------------------------
cycle 4: E = -230.624474239 dE = 5.24875e-07 norm(grad) = 0.000609096
Step 3 : Displace = 6.517e-04/9.983e-04 (rms/max) Trust = 2.000e-01 (+) Grad = 1.758e-04/3.543e-04 (rms/max) E (change) = -230.6244742394 (+5.249e-07) Quality = -1.328
Hessian Eigenvalues: 2.30000e-02 2.30000e-02 2.30000e-02 ... 4.61250e-01 4.61268e-01 5.83445e-01
Converged! =D
#==========================================================================#
#| If this code has benefited your research, please support us by citing: |#
#| |#
#| Wang, L.-P.; Song, C.C. (2016) "Geometry optimization made simple with |#
#| translation and rotation coordinates", J. Chem, Phys. 144, 214108. |#
#| http://dx.doi.org/10.1063/1.4952956 |#
#==========================================================================#
Time elapsed since start of run_optimizer: 14.463 seconds
17.2. Referencias#
Q. Sun, T.C. Berkelbach, N.S. Blunt, G.H. Booth, S. Guo, Z. Li, J. Liu, J.D. McClain, E.R. Sayfutyarova, S. Sharma, S. Wouters, y G.K.-L. Chan, PySCF: the Python-based simulations of chemistry framework, Wiley Interdiscip. Rev. Comput. Mol. Sci. 8, e1340 (2018).
Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N.S. Blunt, N.A. Bogdanov, G.H. Booth, J. Chen, Z.-H. Cui, J.J. Eriksen, Y. Gao, S. Guo, J. Hermann, M.R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J.D. McClain, M. Motta, B. Mussard, H.Q. Pham, A. Pulkin, W. Purwanto, P.J. Robinson, E. Ronca, E.R. Sayfutyarova, M. Scheurer, H.F. Schurkus, J.E.T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L.K. Wagner, X. Wang, A. White, J.D. Whitfield, M.J. Williamson, S. Wouters, J. Yang, J.M. Yu, T. Zhu, T.C. Berkelbach, S. Sharma, A.Y. Sokolov, y G.K.-L. Chan, Recent developments in the PySCF program package, J. Chem. Phys. 153, 024109 (2020).
Q. Sun, Libcint: An efficient general integral library for Gaussian basis functions, J. Comput. Chem. 36, 1664 (2015).