15 Otimização de Geometrias
15.1 Cálculo de forças e gradientes
Na aproximação de Born-Oppenheimer, a energia total de um sistema molecular depende apenas das coordenadas nucleares. A energia total é uma função das coordenadas nucleares, e o gradiente da energia em relação às coordenadas nucleares é chamado de força.
Lembrando que a energia total de um sistema molecular é dada por: \[ E(\{\mathbf{R}_A\}) = \frac{1}{2}\sum_{A=1}^M\sum_{B\neq A}^M \frac{Z_A Z_B}{|\mathbf{R}_A - \mathbf{R}_B|} + E_\text{elec}(\{\mathbf{R}_A\}) \] onde desprezamos a energia cinética nuclear, enquanto a energia eletrônica \(E_\text{elec}(\mathbf{R})\) sendo uma função das coordenadas nucleares \(\mathbf{R} = \{\mathbf{R}_A\}\), e dada por \[ E_\text{elec}(\{\mathbf{R}_A\}) = \min_{\Psi} \langle \Psi | \hat{H}_\text{elec} | \Psi \rangle \] que pode ser calculada usando métodos de química quântica, como Hartree-Fock ou DFT.
O gradiente da energia molecular em relação à coordenada de cada núcleo \(A\) éestá associado à força \(\mathbf{F}_A\) sobre o núcleo \(A\), e é dado por: \[\begin{aligned} \mathbf{F}_A &= -\nabla_{\mathbf{R}_A} E(\{\mathbf{R}_A\}) \\ &= - \sum_{B\neq A}^M \frac{Z_A Z_B (\mathbf{R}_A - \mathbf{R}_B)}{|\mathbf{R}_A - \mathbf{R}_B|^3} -\nabla_{\mathbf{R}_A} E_\text{elec}(\{\mathbf{R}_A\}) \end{aligned} \] onde devemos calcular o gradiente da energia eletrônica em relação às coordenadas nucleares.
O teorema de Hellmann-Feynman nos diz que o gradiente da energia eletrônica em relação às coordenadas nucleares pode ser calculado diretamente a partir da densidade eletrônica \(\rho(\mathbf{r})\) e do potencial nuclear \(V_\text{nuc}(\mathbf{r})\): \[ \nabla_{\mathbf{R}_A} E_\text{elec}(\{\mathbf{R}_A\}) = \min_{\Psi} \langle \Psi | \nabla_{\mathbf{R}_A}\hat{H}_\text{elec} | \Psi \rangle \] e usando a definição de densidade eletrônica, podemos escrever: \[\begin{aligned} \nabla_{\mathbf{R}_A} E_\text{elec}(\{\mathbf{R}_A\}) &= \int \rho(\mathbf{r}) \nabla_{\mathbf{R}_A} V_\text{nuc}(\mathbf{r}) d\mathbf{r}\\ &= \int \rho(\mathbf{r}) \frac{Z_A (\mathbf{r} - \mathbf{R}_A)}{|\mathbf{r} - \mathbf{R}_A|^3} d\mathbf{r} \end{aligned}\]
de modo que a força experimentada pelo núcleo \(A\) pode ser calculada durante a execução de um cálculo de energia eletrônica pelo ciclo SCF, tanto para métodos de Hartree-Fock quanto para DFT.
Assim, para otimizar a geometria de uma molécula, podemos usar o gradiente da energia em relação às coordenadas nucleares para ajustar as posições dos átomos de forma a minimizar a energia total do sistema. No equilíbrio, as forças sobre todos os núcleos devem ser zero, ou seja, o gradiente da energia em relação às coordenadas nucleares deve ser nulo: \[ \nabla_{\mathbf{R}_A} E(\{\mathbf{R}_A\}) = 0 \quad \forall A \] no equilíbrio, a energia total é mínima em relação às coordenadas nucleares.
15.2 Otimização de geometria de moléculas usando PySCF
15.2.1 Modificando a geometria de uma molécula manualmente
15.2.2 Otimizando com o módulo geometric do PySCF
Para executar este notebook no Google Colab, você pode clicar no link: 
No PySCF, você pode realizar uma otimização geométrica usando o módulo geomopt.
Nesse caso devemos instalar o pacote geometric que é uma dependência do PySCF para otimização geométrica. Você pode instalar o pacote geometric usando o seguinte comando:
pip install geometricA seguir, apresentamos um exemplo de como otimizar a geometria da molécula de água (H2O) usando DFT com o funcional B3LYP:
from pyscf import gto, scf
mol_h2o = gto.M(
atom = [['O',(0, 0, 0)], ['H',(0.7586, 0.5043, 0)], ['H',(-0.7586, 0.5043, 0)]],
basis = '6-31++g**')
rks_h2o = scf.RKS(mol_h2o, xc='b3lyp')
e_h2o_dft = rks_h2o.kernel()
from pyscf.geomopt.geometric_solver import optimize
mol_h2o_eq = optimize(rks_h2o)O output da otimização geométrica é bem extenso, mas ao final você verá algo como:
Geometry optimization cycle 5
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
O -0.000000 -0.051815 0.000000 0.000000 0.000093 0.000000
H 0.769480 0.530647 -0.000000 0.000140 0.000097 0.000000
H -0.769480 0.530647 0.000000 -0.000140 0.000097 -0.000000
converged SCF energy = -76.4332283554975
--------------- RKS_Scanner gradients ---------------
x y z
0 O 0.0000000000 -0.0000029586 0.0000000000
1 H -0.0000002211 -0.0000023781 -0.0000000000
2 H 0.0000002211 -0.0000023781 0.0000000000
----------------------------------------------
cycle 5: E = -76.4332283555 dE = -2.89951e-08 norm(grad) = 4.49016e-06
Step 4 : Displace = 1.139e-04/1.395e-04 (rms/max) Trust = 2.828e-01 (+) Grad = 2.592e-06/2.959e-06 (rms/max) E (change) = -76.4332283555 (-2.900e-08) Quality = 0.985
Hessian Eigenvalues: 4.99324e-02 5.00000e-02 5.00000e-02 ... 1.50727e-01 5.37869e-01 6.75229e-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: 6.947 secondsAo final da otimização, podemos imprimir as coordenadas atômicas otimizadas da molécula de água (H2O) em Angstroms e Bohr:
print(mol_h2o_eq.tostring())
print('Atomic coordinates (Ang):')
print(mol_h2o_eq.atom_coords(unit='Ang'))
print('Atomic coordinates (Bohr):')
print(mol_h2o_eq.atom_coords(unit='Bohr'))cujo output será algo como:
O -0.00000000 -0.05181462 0.00000000
H 0.76948042 0.53064706 -0.00000000
H -0.76948042 0.53064706 0.00000000
Atomic coordinates (Ang):
[[-2.60657630e-14 -5.18146164e-02 3.98461609e-14]
[ 7.69480415e-01 5.30647061e-01 -9.36708137e-15]
[-7.69480415e-01 5.30647061e-01 5.80331331e-14]]
Atomic coordinates (Bohr):
[[-4.92571533e-14 -9.79154342e-02 7.52983312e-14]
[ 1.45410724e+00 1.00277761e+00 -1.77012184e-14]
[-1.45410724e+00 1.00277761e+00 1.09666728e-13]]Podemos calcular o comprimento das ligações O-H e o ângulo H-O-H usando as coordenadas atômicas otimizadas:
import numpy as np
coords = mol_h2o_eq.atom_coords(unit='Ang')
lOH1 = np.linalg.norm(coords[0] - coords[1])
lOH2 = np.linalg.norm(coords[0] - coords[2])
print(f'O-H1 bond length (Ang): {lOH1:.3f} Ang')
print(f'O-H2 bond length (Ang): {lOH2:.3f} Ang')
vec1 = coords[1] - coords[0]
vec2 = coords[2] - coords[0]
cos_theta = np.dot(vec1, vec2) / (np.linalg.norm(vec1) * np.linalg.norm(vec2))
theta_rad = np.arccos(cos_theta)
theta_deg = np.degrees(theta_rad)
print(f'H-O-H bond angle (deg): {theta_deg:.3f} deg')cujo output será algo como:
O-H1 bond length (Ang): 0.965 Ang
O-H2 bond length (Ang): 0.965 Ang
H-O-H bond angle (deg): 105.752 deg15.3 Otimização de geometria de moléculas usando VASP e ASE
Para executar este notebook no seu computador, você pode clicar no link: Abrir no GitHub
Usando o ASE vimos que podemos criar a geometria de uma molécula e calcular a energia total usando o VASP.
Criando uma molécula de água (H2O) usando o ASE:
from ase import Atoms, Atom
h2omol = Atoms([Atom('O', [0, 0, 0]),
Atom('H', [0.0, -0.760265, 0.588373]),
Atom('H', [0.0, 0.760265, 0.588373])])
h2omol.center(vacuum=4.0) # caixa com 4 Angstroms de vácuo
h2omol.pbc = True # condição de contorno periódicaEm seguida podemos criar a calculadora do VASP para otimização da geometria usando ibrion=2 (conjugate gradient), isif=0 (relaxa somente os átomos, mantendo a célula fixa), nsw=50 (número máximo de passos de otimização) e ediffg=-1e-3 (tolerância para convergência das forças):
from ase.calculators.vasp import Vasp
calc = Vasp(directory='H2O_relaxed',
xc='PBE', # funcional GGA
encut=350, # safe default for PAW-PBE sets
kpts=[1, 1, 1],gamma=True, # k-points
ibrion=2, # CG ionic relax
isif=0, # relaxa somente os átomos, mantendo a célula fixa
nsw=50, # número máximo de passos de relaxação
lreal='Auto', # projeção de orbitais no espaço real
lwave=True, lcharg=True,lvtot=True, # mantem WAVECAR/CHGCAR/LOCPOT para post-processing
)
h2omol.calc = calcAo calcular a energia total da molécula de água, o VASP irá otimizar a geometria da molécula:
E_h2o = h2omol.get_potential_energy()
print(f'Energia total da molécula de água: {E_h2o:.3f} eV')cujo output será algo como:
Energia total da molécula de água: -14.221 eVO número de passos iônicos executados durante a otimização da geometria pode ser obtido usando:
print("Número de passos iônicos executados:", calc.get_number_of_iterations())cujo output nesse caso foi:
Número de passos iônicos executados: 2E as forças sobre os átomos da molécula de água podem ser obtidas usando:
forces = h2omol.get_forces()
print("Forças sobre os átomos (eV/Ang):")
print(forces)cujo output será algo como:
Forças sobre os átomos (eV/Ang):
[[ 0. 0. 0.00071941]
[-0. 0.00046482 -0.0003597 ]
[-0. -0.00046482 -0.0003597 ]]A nova geometria da molécula de água pode ser obtida usando:
print("Nova geometria da molécula de água (Ang):")
print(h2omol.get_positions())cujo output será algo como:
Nova geometria da molécula de água (Ang):
[[4. 4.760265 3.99497349]
[4. 3.99123747 4.59088625]
[4. 5.52929253 4.59088625]]Os comprimentos das ligações O-H e o ângulo H-O-H podem ser calculados usando as coordenadas atômicas otimizadas:
lOH1 = h2omol.get_distance(0, 1)
lOH2 = h2omol.get_distance(0, 2)
thetaOH1H2 = h2omol.get_angle(1, 0, 2)
print(f"O comprimento da ligação O-H1 é {lOH1:.3f} Å")
print(f"O comprimento da ligação O-H2 é {lOH2:.3f} Å")
print(f"O ângulo H1-O-H2 é {thetaOH1H2:.2f}°")cujo output será algo como:
O comprimento da ligação O-H1 é 0.973 Å
O comprimento da ligação O-H2 é 0.973 Å
O ângulo H1-O-H2 é 104.46°15.4 Cálculo de tensões
Para um sólido, a energia total depende do vetor de rede \(\mathbf{a}, \mathbf{b}, \mathbf{c}\), e o gradiente da energia em relação aos vetores de rede é chamado de tensão. A tensão é uma matriz 3x3, e é dada por: \[ \sigma_{ij} = \frac{1}{V} \frac{\partial E}{\partial \epsilon_{ij}} \] onde \(V\) é o volume da célula unitária, e \(\epsilon_{ij}\) é a deformação aplicada à célula unitária. A tensão pode ser calculada usando métodos de química quântica, como DFT, e é usada para otimizar a geometria de sólidos.
15.5 Otimização de geometria de sólidos usando VASP e ASE
Para executar este notebook no seu computador, você pode clicar no link: Abrir no GitHub
Criando um cristal de NaCl usando o ASE:
from ase.build import bulk
Nacl_crystal = bulk("NaCl", crystalstructure="rocksalt", a=5.64, cubic=True)
print(Nacl_crystal)
print("Cell:", Nacl_crystal.get_cell())
print("Positions:\n", Nacl_crystal.get_positions())cujo output será algo como:
Atoms(symbols='NaClNaClNaClNaCl', pbc=True, cell=[5.64, 5.64, 5.64])
Cell: Cell([5.64, 5.64, 5.64])
Positions:
[[0. 0. 0. ]
[2.82 0. 0. ]
[0. 2.82 2.82]
[2.82 2.82 2.82]
[2.82 0. 2.82]
[0. 0. 2.82]
[2.82 2.82 0. ]
[0. 2.82 0. ]]
Em seguida podemos criar a calculadora do VASP para otimização da geometria usando ibrion=2 (conjugate gradient), isif=7 (relaxa somente a célula), nsw=50 (número máximo de passos de otimização) e ediffg=-1e-3 (tolerância para convergência das forças):
calc = Vasp(directory='NaCl_relaxed',
xc='PBE', # funcional GGA
encut=350, # safe default for PAW-PBE sets
kpts=[1, 1, 1],gamma=True, # k-points
ibrion=2, # CG ionic relax
isif=7, # relaxa somente a célula
nsw=50, # número máximo de passos de relaxação
ediffg=-1e-3, # critério de convergência para relaxação (forças)
lreal='Auto', # projeção de orbitais no espaço real
lwave=True, lcharg=True,lvtot=True, # mantem WAVECAR/CHGCAR/LOCPOT para post-processing
)
Nacl_crystal.calc = calcA otimização da geometria do cristal de NaCl pode ser realizada usando o método get_potential_energy() do objeto Nacl_crystal, que irá calcular a energia total e otimizar a geometria do cristal:
E_nacl = Nacl_crystal.get_potential_energy()
print(f'Energia total do cristal de NaCl: {E_nacl:.3f} eV')cujo output será algo como:
Energia total do cristal de NaCl: -25.451 eVO número de passos de otimização da célula executados durante a otimização da geometria do cristal de NaCl pode ser obtido usando:
print("Número de passos de otimização da célula executados:", calc.get_number_of_iterations())cujo output nesse caso foi:
Número de passos de otimização da célula executados: 4A nova estrutura do cristal de NaCl pode ser obtida usando:
print("Nova estrutura do cristal de NaCl (Ang):")
print(Nacl_crystal.get_cell())cujo output será algo como:
Nova estrutura do cristal de NaCl (Ang):
Cell([5.911942143216201, 5.911942143216201, 5.911942143216201])