書記が物理やるだけ#421 密度汎関数法の計算-2
問題
解答
DFTにおいて，全エネルギーの最も重要な近似項は交換-相関エネルギーであり，LDAとGGAはこの項を電子密度の情報を用いて表現するための基本的な近似手法である。
LDA（Local Density Approximation: 局所密度近似）では，系の各点 r における電子密度が均一な電子ガス (Homogeneous Electron Gas, HEG) のそれに局所的に等しいと仮定する。これにより計算コストを抑えることができる。
GGA（Generalized Gradient Approximation: 一般化勾配近似）では，系のエネルギー密度が，局所的な電子密度だけでなく，その空間的な変化の度合いすなわち電子密度の勾配にも依存すると仮定することで精度を上げる。
B3LYPは，GGAに分類される汎関数よりもさらに高度な近似であるハイブリッド汎関数の代表例である。B3LYPの交換-相関エネルギーは，主に以下の4つの異なる項の線形結合として定義される；ハートリー・フォック交換項，LDA交換項，GGA交換項(Becke 88)，GGA相関項 (Lee-Yang-Parr)。ここでHF交換項を20%導入することで，LDAや純粋なGGAでは困難だった分子系の熱化学的性質や構造予測において精度を上げることができる。
以下のPythonコードを実行することで，水分子のB3LYP/6-31Gレベルでの平衡構造最適化とエネルギー計算が実行できる。
import pyscf.gto
import pyscf.scf
import pyscf.geomopt
# 1. 分子と基底関数の定義
# 初期構造（水分子）
mol = pyscf.gto.M(
atom = '''
O 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 -0.7579180000 0.5866170000
H 0.0000000000 0.7579180000 0.5866170000
''',
basis = '6-31g', # 基底関数系
charge = 0,
spin = 0,
)
# 2. DFT計算の設定 (RKS: Restricted Kohn-Sham)
mf = pyscf.scf.RKS(mol)
mf.xc = 'b3lyp'
# 3. 構造最適化の実行
print('--- Geometry Optimization Start (B3LYP/6-31G) ---')
# 'mf'オブジェクトを渡すことで、指定されたDFTレベルで最適化が行われる
# 注意: geometricライブラリが別途インストールされている必要があります (pip install geometric)
optimized_mol = pyscf.geomopt.optimize(mf)
print('--- Geometry Optimization Finished ---')
# 4. 最適化された構造での最終エネルギー計算
final_mf = pyscf.scf.RKS(optimized_mol)
final_mf.xc = 'b3lyp'
final_energy = final_mf.kernel() # SCF計算の実行
print('\n--- Final Optimized Structure & Energy ---')
print(f'Final SCF Energy (Hartree): {final_energy:.10f}')
# 最適化された原子座標の表示
print('\nOptimized Atomic Coordinates (Bohr/Angstrom):')
# PySCFのMoleオブジェクトからデカルト座標を取得
print(optimized_mol.atom_coords())結果：
geometric-optimize called with the following command line:
/opt/anaconda3/lib/python3.9/site-packages/ipykernel_launcher.py -f /Users/ryosukeokubo/Library/Jupyter/runtime/kernel-5ac82df7-4714-4ead-bae5-d7cabe32e79a.json
())))))))))))))))/
())))))))))))))))))))))))),
*)))))))))))))))))))))))))))))))))
#, ()))))))))/ .)))))))))),
#%%%%, ()))))) .))))))))*
*%%%%%%, )) .. ,))))))).
*%%%%%%, ***************/. .)))))))
#%%/ (%%%%%%, /*********************. )))))))
.%%%%%%# *%%%%%%, *******/, **********, .))))))
.%%%%%%/ *%%%%%%, ** ******** .))))))
## .%%%%%%/ (%%%%%%, ,****** /)))))
%%%%%% .%%%%%%# *%%%%%%, ,/////. ****** ))))))
#% %% .%%%%%%/ *%%%%%%, ////////, *****/ ,)))))
#%% %%% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
#%%%%. %%%%%# /%%%%%%* #%%%%%% /////) ****** ))))),
#%%%%##% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
## %%% .%%%%%%/ *%%%%%%, ////////. *****/ ,)))))
#%%%%# /%%%%%%/ (%%%%%% /)/)// ****** ))))))
## .%%%%%%/ (%%%%%%, ******* ))))))
.%%%%%%/ *%%%%%%, **. /******* .))))))
*%%%%%%/ (%%%%%% ********/*..,*/********* *))))))
#%%/ (%%%%%%, *********************/ )))))))
*%%%%%%, ,**************/ ,))))))/
(%%%%%% () ))))))))
#%%%%, ()))))) ,)))))))),
#, ()))))))))) ,)))))))))).
()))))))))))))))))))))))))))))))/
())))))))))))))))))))))))).
())))))))))))))),
-=# geomeTRIC started. Version: 1.1 #=-
Current date and time: 2025-11-26 11:22:26
#========================================================#
#| Arguments passed to driver run_optimizer(): |#
#========================================================#
customengine <pyscf.geomopt.geometric_solver.PySCFEngine object at 0x7fd72c037dc0>
input /var/folders/5x/frnl8pn17pl4mfm2qfl31l140000gn/T/tmpwdbjcndt/0cadc1e7-7c3f-4311-8854-7d1a20b1c47b
logIni /opt/anaconda3/lib/python3.9/site-packages/pyscf/geomopt/log.ini
maxiter 100
----------------------------------------------------------
Custom engine selected.
Bonds will be generated from interatomic distances less than 1.20 times sum of covalent radii
--- Geometry Optimization Start (B3LYP/6-31G) ---
9 internal coordinates being used (instead of 9 Cartesians)
Internal coordinate system (atoms numbered from 1):
Distance 1-2
Distance 1-3
Angle 2-1-3
Translation-X 1-3
Translation-Y 1-3
Translation-Z 1-3
Rotation-A 1-3
Rotation-B 1-3
Rotation-C 1-3
<class 'geometric.internal.Distance'> : 2
<class 'geometric.internal.Angle'> : 1
<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: 100
> 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
O 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
H 0.000000 -0.757918 0.586617 0.000000 0.000000 0.000000
H 0.000000 0.757918 0.586617 0.000000 0.000000 0.000000
WARN: Mole.unit (angstrom) is changed to Bohr
converged SCF energy = -76.3850100867294
--------------- RKS_Scanner gradients ---------------
x y z
0 O -0.0000000000 -0.0000000000 0.0117583333
1 H -0.0000000000 0.0205089252 -0.0058829593
2 H 0.0000000000 -0.0205089252 -0.0058829593
----------------------------------------------
cycle 1: E = -76.3850100867 dE = -76.385 norm(grad) = 0.0323838
Step 0 : Gradient = 1.870e-02/2.134e-02 (rms/max) Energy = -76.3850100867
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.60000e-01 5.58157e-01 5.58157e-01
Geometry optimization cycle 2
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
O 0.000000 0.000000 0.015848 0.000000 0.000000 0.015848
H 0.000000 -0.798747 0.578814 0.000000 -0.040829 -0.007803
H 0.000000 0.798747 0.578814 0.000000 0.040829 -0.007803
converged SCF energy = -76.3860634307406
--------------- RKS_Scanner gradients ---------------
x y z
0 O -0.0000000000 0.0000000000 0.0010165518
1 H 0.0000000000 -0.0030742784 -0.0005113316
2 H 0.0000000000 0.0030742784 -0.0005113316
----------------------------------------------
cycle 2: E = -76.3860634307 dE = -0.00105334 norm(grad) = 0.00452313
Step 1 : Displace = 3.515e-02/4.158e-02 (rms/max) Trust = 1.000e-01 (=) Grad = 2.611e-03/3.117e-03 (rms/max) E (change) = -76.3860634307 (-1.053e-03) Quality = 0.784
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 1.86595e-01 5.58157e-01 5.71017e-01
Geometry optimization cycle 3
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
O 0.000000 0.000000 0.011060 0.000000 -0.000000 -0.004788
H 0.000000 -0.791522 0.581290 0.000000 0.007225 0.002476
H 0.000000 0.791522 0.581290 0.000000 -0.007225 0.002476
converged SCF energy = -76.3861160263633
--------------- RKS_Scanner gradients ---------------
x y z
0 O 0.0000000000 0.0000000000 0.0005515418
1 H -0.0000000000 0.0000150349 -0.0002793299
2 H 0.0000000000 -0.0000150349 -0.0002793299
----------------------------------------------
cycle 3: E = -76.3861160264 dE = -5.25956e-05 norm(grad) = 0.000678749
Step 2 : Displace = 6.821e-03/7.620e-03 (rms/max) Trust = 1.414e-01 (+) Grad = 3.919e-04/5.515e-04 (rms/max) E (change) = -76.3861160264 (-5.260e-05) Quality = 1.077
Hessian Eigenvalues: 4.99990e-02 5.00000e-02 5.00000e-02 ... 1.65291e-01 5.58157e-01 5.97449e-01
Geometry optimization cycle 4
Cartesian coordinates (Angstrom)
Atom New coordinates dX dY dZ
O 0.000000 -0.000000 0.010403 -0.000000 -0.000000 -0.000656
H 0.000000 -0.790984 0.581740 0.000000 0.000539 0.000450
H 0.000000 0.790984 0.581740 -0.000000 -0.000539 0.000450
converged SCF energy = -76.3861166637285
--------------- RKS_Scanner gradients ---------------
x y z
0 O -0.0000000000 -0.0000000000 0.0001247536
1 H 0.0000000000 0.0000669948 -0.0000659979
2 H -0.0000000000 -0.0000669948 -0.0000659979
----------------------------------------------
cycle 4: E = -76.3861166637 dE = -6.37365e-07 norm(grad) = 0.00018235
Step 3 : Displace = 6.822e-04/7.375e-04 (rms/max) Trust = 2.000e-01 (+) Grad = 1.053e-04/1.248e-04 (rms/max) E (change) = -76.3861166637 (-6.374e-07) Quality = 1.129
Hessian Eigenvalues: 4.99990e-02 5.00000e-02 5.00000e-02 ... 1.65291e-01 5.58157e-01 5.97449e-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: 2.637 seconds
--- Geometry Optimization Finished ---
converged SCF energy = -76.3861166637351
--- Final Optimized Structure & Energy ---
Final SCF Energy (Hartree): -76.3861166637
Optimized Atomic Coordinates (Bohr/Angstrom):
[[ 2.89382989e-14 -1.83696495e-14 1.96593584e-02]
[ 2.53322948e-14 -1.49474228e+00 1.09932835e+00]
[ 4.91135093e-15 1.49474228e+00 1.09932835e+00]]最適化プロセスはサイクル 4 で収束し，平衡構造における基底状態エネルギー，最終的な原子核の位置が得られた。
本記事のもくじはこちら：
