第一原理分子動力学プログラム STATE Senri Wiki
開始行:
*自由エネルギーの計算 [#a3b50620]
ここでは最適化された構造を用いthermochemistryモジュールの...
- Pythonスクリプト (ase-vib+thermo.py)
from ase.io import read
from ase.io import write
from ase.calculators.espresso import Espresso
from ase.optimize import BFGS
from ase.vibrations import Vibrations
from ase.thermochemistry import IdealGasThermo
input_data = {'control': {'outdir' : './tmp/',
'prefix' : 'h2'
},
'system': {'ecutwfc': 80,
'ecutrho': 800},
'electrons': {'diagonalization'...
'conv_thr' : 1.0e-12,
'mixing_beta' : 0.3}
}
pseudopotentials = {'H':'H.pbe-kjpaw_psl.1.0.0.UPF'}
calc = Espresso(pseudopotentials=pseudopotentials,input_...
tprnfor=True)
h2=read('scf.in',format='espresso-in')
h2.calc = calc
e_h2 = h2.get_potential_energy()
vib = Vibrations(h2)
vib.run()
vib.summary(log='H2_vib_summary.txt')
vib.write_mode(-1)
vib_energies = vib.get_energies()
thermo = IdealGasThermo(vib_energies=vib_energies,
potentialene...
atoms=h2,
geometry='li...
symmetrynumb...
G = thermo.get_gibbs_energy(temperature=298.15, pressure...
ここで計算に使用する構造は以下のようにQEフォーマットのscf...
h2=read('scf.in',format='espresso-in')
先ず全エネルギーのSCF計算を行っています。
e_h2 = h2.get_potential_energy()
その後、自由エネルギーの計算に必要な(調和)振動数の計算...
vib = Vibrations(h2)
vib.run()
vib.summary(log='H2_vib_summary.txt')
vib.write_mode(-1)
ここで得られた振動数と零点エネルギーは以下のようにH2_vib_...
---------------------
# meV cm^-1
---------------------
0 0.0i 0.0i
1 0.0 0.0
2 0.0 0.0
3 4.5 35.9
4 4.5 35.9
5 535.2 4316.3
---------------------
Zero-point energy: 0.272 eV
そして自由エネルギーに必要な計算を実行しています。
thermo = IdealGasThermo(vib_energies=vib_energies,
potentialene...
atoms=h2,
geometry='li...
symmetrynumb...
symmetrynumberは[['''Essentials of Computational Chemistr...
最終的に298.15 K、101325 Pa (101.325 kPa (~0.1 MPa))にお...
G = thermo.get_gibbs_energy(temperature=298.15, pressure...
Enthalpy components at T = 298.15 K:
===============================
E_pot -31.747 eV
E_ZPE 0.268 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.026 eV
Cv_vib (0->T) 0.000 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -31.389 eV
===============================
Entropy components at T = 298.15 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0012188 eV/K 0.363 eV
S_rot 0.0001341 eV/K 0.040 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0000000 eV/K 0.000 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0013518 eV/K 0.403 eV
=================================================
Free energy components at T = 298.15 K and P = 101325.0 ...
=======================
H -31.389 eV
-T*S -0.403 eV
-----------------------
G -31.792 eV
=======================
** 参考 [#ddba03f3]
- [[Thermochemistry>https://wiki.fysik.dtu.dk/ase/ase/the...
終了行:
*自由エネルギーの計算 [#a3b50620]
ここでは最適化された構造を用いthermochemistryモジュールの...
- Pythonスクリプト (ase-vib+thermo.py)
from ase.io import read
from ase.io import write
from ase.calculators.espresso import Espresso
from ase.optimize import BFGS
from ase.vibrations import Vibrations
from ase.thermochemistry import IdealGasThermo
input_data = {'control': {'outdir' : './tmp/',
'prefix' : 'h2'
},
'system': {'ecutwfc': 80,
'ecutrho': 800},
'electrons': {'diagonalization'...
'conv_thr' : 1.0e-12,
'mixing_beta' : 0.3}
}
pseudopotentials = {'H':'H.pbe-kjpaw_psl.1.0.0.UPF'}
calc = Espresso(pseudopotentials=pseudopotentials,input_...
tprnfor=True)
h2=read('scf.in',format='espresso-in')
h2.calc = calc
e_h2 = h2.get_potential_energy()
vib = Vibrations(h2)
vib.run()
vib.summary(log='H2_vib_summary.txt')
vib.write_mode(-1)
vib_energies = vib.get_energies()
thermo = IdealGasThermo(vib_energies=vib_energies,
potentialene...
atoms=h2,
geometry='li...
symmetrynumb...
G = thermo.get_gibbs_energy(temperature=298.15, pressure...
ここで計算に使用する構造は以下のようにQEフォーマットのscf...
h2=read('scf.in',format='espresso-in')
先ず全エネルギーのSCF計算を行っています。
e_h2 = h2.get_potential_energy()
その後、自由エネルギーの計算に必要な(調和)振動数の計算...
vib = Vibrations(h2)
vib.run()
vib.summary(log='H2_vib_summary.txt')
vib.write_mode(-1)
ここで得られた振動数と零点エネルギーは以下のようにH2_vib_...
---------------------
# meV cm^-1
---------------------
0 0.0i 0.0i
1 0.0 0.0
2 0.0 0.0
3 4.5 35.9
4 4.5 35.9
5 535.2 4316.3
---------------------
Zero-point energy: 0.272 eV
そして自由エネルギーに必要な計算を実行しています。
thermo = IdealGasThermo(vib_energies=vib_energies,
potentialene...
atoms=h2,
geometry='li...
symmetrynumb...
symmetrynumberは[['''Essentials of Computational Chemistr...
最終的に298.15 K、101325 Pa (101.325 kPa (~0.1 MPa))にお...
G = thermo.get_gibbs_energy(temperature=298.15, pressure...
Enthalpy components at T = 298.15 K:
===============================
E_pot -31.747 eV
E_ZPE 0.268 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.026 eV
Cv_vib (0->T) 0.000 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -31.389 eV
===============================
Entropy components at T = 298.15 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0012188 eV/K 0.363 eV
S_rot 0.0001341 eV/K 0.040 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0000000 eV/K 0.000 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0013518 eV/K 0.403 eV
=================================================
Free energy components at T = 298.15 K and P = 101325.0 ...
=======================
H -31.389 eV
-T*S -0.403 eV
-----------------------
G -31.792 eV
=======================
** 参考 [#ddba03f3]
- [[Thermochemistry>https://wiki.fysik.dtu.dk/ase/ase/the...
ページ名: