Geometry optimization ¶
In [1]:
Copied!
from pyscf import gto, scf, mp
from pyscf.geomopt.geometric_solver import optimize
import py3Dmol
from pyscf import gto, scf, mp
from pyscf.geomopt.geometric_solver import optimize
import py3Dmol
In [2]:
Copied!
mol = gto.Mole()
mol.atom = """
H 1.090669 0.975611 0.000000
H -0.435544 1.078064 0.889793
H -0.435544 1.078064 -0.889793
C 0.048224 0.665550 0.000000
O 0.044194 -0.754448 0.000000
H -0.874480 -1.043721 0.000000
"""
mol.unit = "Angstrom"
mol.basis = "cc-pvtz"
mol.build()
mol = gto.Mole()
mol.atom = """
H 1.090669 0.975611 0.000000
H -0.435544 1.078064 0.889793
H -0.435544 1.078064 -0.889793
C 0.048224 0.665550 0.000000
O 0.044194 -0.754448 0.000000
H -0.874480 -1.043721 0.000000
"""
mol.unit = "Angstrom"
mol.basis = "cc-pvtz"
mol.build()
Out[2]:
<pyscf.gto.mole.Mole at 0x7aab94c83290>
In [3]:
Copied!
mf = scf.RHF(mol).run()
mp2 = mp.MP2(mf).run()
mf = scf.RHF(mol).run()
mp2 = mp.MP2(mf).run()
converged SCF energy = -115.089057168863 E(MP2) = -115.546749317889 E_corr = -0.457692149025912 E(SCS-MP2) = -115.548492148244 E_corr = -0.459434979380392
In [4]:
Copied!
opt = mp2.Gradients().optimizer(solver="geometric")
opt = mp2.Gradients().optimizer(solver="geometric")
In [5]:
Copied!
mol_eq = opt.kernel()
mol_eq = opt.kernel()
geometric-optimize called with the following command line:
/home/alex/miniconda3/envs/pyscf-dev/lib/python3.12/site-packages/ipykernel_launcher.py --f=/home/alex/.local/share/jupyter/runtime/kernel-v2-13676DEPx0bbU6O5e.json
())))))))))))))))/
())))))))))))))))))))))))),
*)))))))))))))))))))))))))))))))))
#, ()))))))))/ .)))))))))),
#%%%%, ()))))) .))))))))*
*%%%%%%, )) .. ,))))))).
*%%%%%%, ***************/. .)))))))
#%%/ (%%%%%%, /*********************. )))))))
.%%%%%%# *%%%%%%, *******/, **********, .))))))
.%%%%%%/ *%%%%%%, ** ******** .))))))
## .%%%%%%/ (%%%%%%, ,****** /)))))
%%%%%% .%%%%%%# *%%%%%%, ,/////. ****** ))))))
#% %% .%%%%%%/ *%%%%%%, ////////, *****/ ,)))))
#%% %%% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
#%%%%. %%%%%# /%%%%%%* #%%%%%% /////) ****** ))))),
#%%%%##% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
## %%% .%%%%%%/ *%%%%%%, ////////. *****/ ,)))))
#%%%%# /%%%%%%/ (%%%%%% /)/)// ****** ))))))
## .%%%%%%/ (%%%%%%, ******* ))))))
.%%%%%%/ *%%%%%%, **. /******* .))))))
*%%%%%%/ (%%%%%% ********/*..,*/********* *))))))
#%%/ (%%%%%%, *********************/ )))))))
*%%%%%%, ,**************/ ,))))))/
(%%%%%% () ))))))))
#%%%%, ()))))) ,)))))))),
#, ()))))))))) ,)))))))))).
()))))))))))))))))))))))))))))))/
())))))))))))))))))))))))).
())))))))))))))),
-=# geomeTRIC started. Version: 1.0 #=-
Current date and time: 2024-03-23 21:20:54
Custom engine selected.
Bonds will be generated from interatomic distances less than 1.20 times sum of covalent radii
18 internal coordinates being used (instead of 18 Cartesians)
Internal coordinate system (atoms numbered from 1):
Distance 1-4
Distance 2-4
Distance 3-4
Distance 4-5
Distance 5-6
Angle 1-4-2
Angle 1-4-3
Angle 1-4-5
Angle 2-4-3
Angle 2-4-5
Angle 3-4-5
Angle 4-5-6
Dihedral 1-4-5-6
Dihedral 2-4-5-6
Dihedral 3-4-5-6
Translation-X 1-6
Translation-Y 1-6
Translation-Z 1-6
Rotation-A 1-6
Rotation-B 1-6
Rotation-C 1-6
<class 'geometric.internal.Distance'> : 5
<class 'geometric.internal.Angle'> : 7
<class 'geometric.internal.Dihedral'> : 3
<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
H 1.090669 0.975611 0.000000 0.000000 0.000000 0.000000
H -0.435544 1.078064 0.889793 0.000000 0.000000 -0.000000
H -0.435544 1.078064 -0.889793 0.000000 0.000000 0.000000
C 0.048224 0.665550 0.000000 0.000000 0.000000 0.000000
O 0.044194 -0.754448 0.000000 0.000000 0.000000 0.000000
H -0.874480 -1.043721 0.000000 0.000000 0.000000 0.000000
converged SCF energy = -115.089057168865
E(MP2_Scanner) = -115.546749322394 E_corr = -0.457692153528722
E(SCS-MP2_Scanner) = -115.548492145681 E_corr = -0.459434976815418
--------------- MP2_Scanner gradients ---------------
x y z
0 H 0.0037113235 0.0009762788 0.0000000000
1 H -0.0019602440 0.0009230030 0.0032228600
2 H -0.0019602440 0.0009230030 -0.0032228600
3 C 0.0009918653 0.0009070122 -0.0000000000
4 O 0.0041423323 -0.0016285445 0.0000000000
5 H -0.0049250330 -0.0021007524 -0.0000000000
----------------------------------------------
cycle 1: E = -115.546749322 dE = -115.547 norm(grad) = 0.00975584
Step 0 : Gradient = 3.983e-03/5.354e-03 (rms/max) Energy = -115.5467493224 Hessian Eigenvalues: 2.30000e-02 5.00000e-02 5.00000e-02 ... 3.50939e-01 4.24598e-01 5.48032e-01
Geometry optimization cycle 2
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 1.082552 0.971601 0.000000 -0.008117 -0.004010 0.000000
H -0.432970 1.076985 0.884536 0.002574 -0.001079 -0.005257
H -0.432970 1.076985 -0.884536 0.002574 -0.001079 0.005257
C 0.045743 0.662878 0.000000 -0.002481 -0.002672 0.000000
O 0.044989 -0.752475 0.000000 0.000795 0.001973 0.000000
H -0.869825 -1.036854 0.000000 0.004655 0.006867 0.000000
converged SCF energy = -115.089435766108
E(MP2_Scanner) = -115.546859964084 E_corr = -0.457424197976251
E(SCS-MP2_Scanner) = -115.548558979997 E_corr = -0.459123213889003
--------------- MP2_Scanner gradients ---------------
x y z
0 H -0.0001877465 -0.0004649158 0.0000000000
1 H 0.0000457271 -0.0000837216 -0.0001258365
2 H 0.0000457271 -0.0000837216 0.0001258365
3 C -0.0005118004 0.0008954399 -0.0000000000
4 O 0.0007734022 -0.0007348676 -0.0000000000
5 H -0.0001653095 0.0004717867 -0.0000000000
----------------------------------------------
cycle 2: E = -115.546859964 dE = -0.000110642 norm(grad) = 0.00165925
Step 1 : Displace = 6.317e-03/9.054e-03 (rms/max) Trust = 1.000e-01 (=) Grad = 6.774e-04/1.067e-03 (rms/max) E (change) = -115.5468599641 (-1.106e-04) Quality = 0.959 Hessian Eigenvalues: 2.30000e-02 5.00000e-02 5.00000e-02 ... 3.66299e-01 4.14907e-01 5.44078e-01
Geometry optimization cycle 3
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 1.082897 0.973519 0.000000 0.000345 0.001918 -0.000000
H -0.432286 1.077334 0.884344 0.000684 0.000349 -0.000192
H -0.432286 1.077334 -0.884344 0.000684 0.000349 0.000192
C 0.046431 0.662731 0.000000 0.000688 -0.000147 -0.000000
O 0.043206 -0.752154 0.000000 -0.001783 0.000321 -0.000000
H -0.870442 -1.039645 0.000000 -0.000617 -0.002791 -0.000000
converged SCF energy = -115.08949678528
E(MP2_Scanner) = -115.546863613983 E_corr = -0.457366828702596
E(SCS-MP2_Scanner) = -115.548563758329 E_corr = -0.459066973048177
--------------- MP2_Scanner gradients ---------------
x y z
0 H -0.0001328572 0.0000353638 -0.0000000000
1 H 0.0000557015 -0.0001064424 -0.0001281383
2 H 0.0000557015 -0.0001064424 0.0001281383
3 C 0.0001528682 0.0005858007 0.0000000000
4 O -0.0001884534 -0.0003703000 -0.0000000000
5 H 0.0000570395 -0.0000379797 0.0000000000
----------------------------------------------
cycle 3: E = -115.546863614 dE = -3.6499e-06 norm(grad) = 0.000790234
Step 2 : Displace = 1.683e-03/2.861e-03 (rms/max) Trust = 1.414e-01 (+) Grad = 3.226e-04/6.054e-04 (rms/max) E (change) = -115.5468636140 (-3.650e-06) Quality = 0.916
Hessian Eigenvalues: 2.30000e-02 5.00000e-02 5.00000e-02 ... 3.49406e-01 4.33789e-01 5.43616e-01
Geometry optimization cycle 4
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
H 1.082943 0.973054 0.000000 0.000046 -0.000465 0.000000
H -0.432390 1.077525 0.884539 -0.000104 0.000190 0.000195
H -0.432390 1.077525 -0.884539 -0.000104 0.000190 -0.000195
C 0.046194 0.662504 0.000000 -0.000237 -0.000227 0.000000
O 0.043388 -0.751885 0.000000 0.000182 0.000268 -0.000000
H -0.870226 -1.039603 -0.000000 0.000216 0.000042 -0.000000
converged SCF energy = -115.089498861653
E(MP2_Scanner) = -115.546864093532 E_corr = -0.45736523187897
E(SCS-MP2_Scanner) = -115.548562835516 E_corr = -0.459063973862505
--------------- MP2_Scanner gradients ---------------
x y z
0 H 0.0000030778 -0.0000065800 0.0000000000
1 H 0.0000037381 -0.0000278095 0.0000000211
2 H 0.0000037381 -0.0000278095 -0.0000000211
3 C 0.0000082656 0.0001563074 -0.0000000000
4 O -0.0000343972 -0.0000868525 0.0000000000
5 H 0.0000155777 -0.0000072558 -0.0000000000
----------------------------------------------
cycle 4: E = -115.546864094 dE = -4.79549e-07 norm(grad) = 0.000187482
Step 3 : Displace = 3.290e-04/4.659e-04 (rms/max) Trust = 2.000e-01 (+) Grad = 7.654e-05/1.565e-04 (rms/max) E (change) = -115.5468640935 (-4.795e-07) Quality = 1.218 Hessian Eigenvalues: 2.30000e-02 5.00000e-02 5.00000e-02 ... 3.49406e-01 4.33789e-01 5.43616e-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: 24.833 seconds
In [6]:
Copied!
xyz = mol_eq.atom_coords(unit="Angstrom")
atom_symbols = [mol_eq.atom_pure_symbol(i) for i in range(mol_eq.natm)]
xyz_str = f"{len(atom_symbols)}\n\n" + "".join(
[
f"{atom_symbols[i]} {xyz[i][0]:.12f} {xyz[i][1]:.12f} {xyz[i][2]:.12f}\n"
for i in range(len(atom_symbols))
]
)
print(xyz_str)
xyz = mol_eq.atom_coords(unit="Angstrom")
atom_symbols = [mol_eq.atom_pure_symbol(i) for i in range(mol_eq.natm)]
xyz_str = f"{len(atom_symbols)}\n\n" + "".join(
[
f"{atom_symbols[i]} {xyz[i][0]:.12f} {xyz[i][1]:.12f} {xyz[i][2]:.12f}\n"
for i in range(len(atom_symbols))
]
)
print(xyz_str)
6 H 1.082942730243 0.973054295538 0.000000000000 H -0.432389937115 1.077524913433 0.884539060760 H -0.432389937115 1.077524913433 -0.884539060759 C 0.046194206046 0.662504379510 0.000000000000 O 0.043388186215 -0.751885455981 0.000000000000 H -0.870226248270 -1.039603045920 -0.000000000000
In [7]:
Copied!
view = py3Dmol.view(width=400, height=400)
view.addModel(xyz_str, "xyz")
view.setStyle({"stick": {}})
view.zoomTo()
view.show()
view = py3Dmol.view(width=400, height=400)
view.addModel(xyz_str, "xyz")
view.setStyle({"stick": {}})
view.zoomTo()
view.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
The final structure in Z-matrix would be
H
H 1 1.7577
H 1 1.7577 2 60.43
C 1 1.0823 2 36.02 3 323.22
O 4 1.4144 1 106.79 2 238.84
H 5 0.9578 4 107.59 1 180.00