Turbomole を用いた量子化学計算の基本的な手順：エチレンの DFT 計算を例として
- Turbomole は量子化学計算で用いられるプログラムパッケージのひとつです
- このノートでは、エチレン分子を例にして、Turbomole を使って構造最適化、振動数解析、励起状態計算に取り組むときの手順を紹介します
基本的な手順
- Turbomole の入力ファイル（
controlなど）は、代表的な量子化学計算パッケージの Gaussian などと比べると、かなり複雑です - したがって基本的には、Turbomole の
defineツールを使って、計算の条件などを対話的に設定します defineツールで基本的な設定をした後に、内容ごとに異なるプログラムを用いて計算を実行します- たとえば、一点計算は
dscf、構造最適化はjobex、振動数解析はaoforce、励起状態計算はescfなど
- たとえば、一点計算は
分子研スパコンを使っている場合
- ログイン後、
moduleコマンドで、Turbomoleの環境設定をおこないます
module load turbomole初期構造の準備
- エチレン分子の初期構造を
ethylene.xyzというファイルに保存します
ethylene.xyz
6
C2H4 molecule
C 0.00 0.00 0.00
C -1.31 0.00 0.00
H -1.85 -0.94 0.00
H -1.85 0.94 0.00
H 0.54 0.94 0.00
H 0.54 -0.94 0.00- Turmobole に付属する
x2tツールを使って、xyz形式から、Turbomole 用の座標形式に変換します
x2t ethylene.xyz | tee coordcoord
$coord
0.00691780724847 -0.00000000001683 -0.00000000001018 c
-2.48245903052942 0.00000000001673 0.00000000001010 c
-3.55380171099544 -1.72798218518559 -0.00000000000268 h
-3.55380171101130 1.72798218518850 -0.00000000000269 h
1.07826048771524 1.72798218518535 0.00000000000274 h
1.07826048773116 -1.72798218518815 0.00000000000271 h
$user-defined bonds
$end構造最適化
- 構造最適化では、分子に掛かる力の方向に少しずつ動かして、安定な（＝ポテンシャルエネルギーが小さい）構造を探索します
define の起動
- ターミナルで
defineを実行して、対話型セッションを開始します
define ***********************************************************
* *
* D E F I N E *
* *
* TURBOMOLE'S INTERACTIVE INPUT PROGRAM *
* *
* Quantum Chemistry Group University of Karlsruhe *
* *
***********************************************************
DATA WILL BE WRITTEN TO THE NEW FILE control
IF YOU WANT TO READ DEFAULT-DATA FROM ANOTHER control-TYPE FILE,
THEN ENTER ITS LOCATION/NAME OR OTHERWISE HIT >return<.- INPUT: Enter
INPUT TITLE OR
ENTER & TO REPEAT DEFINITION OF DEFAULT INPUT FILE- INPUT:
ethylene- このインプットのタイトルを ethylene にします
座標の設定
- まずは、計算に用いる分子の立体構造を指定します
symmetry group of the molecule : c1
the group has the following generators :
c1(z)
1 symmetry operations found
SPECIFICATION OF MOLECULAR GEOMETRY ( #ATOMS=0 SYMMETRY=c1 )
YOU MAY USE ONE OF THE FOLLOWING COMMANDS :
sy <group> <eps> : DEFINE MOLECULAR SYMMETRY (default for eps=3d-1)
desy <eps> : DETERMINE MOLECULAR SYMMETRY AND ADJUST
COORDINATES (default for eps=1d-6)
syndi <eps> : LIKE DESY, BUT FIND ONLY GROUPS WITH NON-
DEGENERATE IRREPS (D2h AND SUBGROUPS)
susy : ADJUST COORDINATES FOR SUBGROUPS
ai : ADD ATOMIC COORDINATES INTERACTIVELY
a <file> : ADD ATOMIC COORDINATES FROM FILE <file>
aa <file> : ADD ATOMIC COORDINATES IN ANGSTROEM UNITS FROM FILE <file>
sub : SUBSTITUTE AN ATOM BY A GROUP OF ATOMS
i : INTERNAL COORDINATE MENU
ired : REDUNDANT INTERNAL COORDINATES
pbc_ired : PERIODIC REDUNDANT INTERNAL COORDINATES
red_info : DISPLAY REDUNDANT INTERNAL COORDINATES
ff : UFF-FORCEFIELD CALCULATION
m : MANIPULATE GEOMETRY
frag : Define Fragments for BSSE calculation
w <file> : WRITE MOLECULAR COORDINATES TO FILE <file>
r <file> : RELOAD ATOMIC AND INTERNAL COORDINATES FROM FILE <file>
name : CHANGE ATOMIC IDENTIFIERS
del : DELETE ATOMS
fix : FIX ATOMS
dis : DISPLAY MOLECULAR GEOMETRY
banal : CARRY OUT BOND ANALYSIS
* : TERMINATE MOLECULAR GEOMETRY SPECIFICATION
AND WRITE GEOMETRY DATA TO CONTROL FILE- INPUT:
a coord- 座標ファイルを読み込みます
CARTESIAN COORDINATES FOR 6 ATOMS HAVE SUCCESSFULLY
BEEN ADDED.
DEFINITIONS OF INTERNAL COORDINATES HAVE N O T BEEN READ.- INPUT:
desy- 対称性を自動で設定します
molecule is planar
2 symmetry operations found :
given set of 6 coordinate triples separated
into 2 bin(s) for a given accuracy of 0.5000D-06
no atom resides at center of mass
try inversion ...
4 symmetry operations found :
8 symmetry operations found :
thrsym = 1.000000000000000E-006
c2 axis orthogonal to cn axis : 5
c2 axis orthogonal to sn axis : 4
mirror plane orthogonal to cn/sn axis : 2
mirror plane parallel to cn/sn axis : 2
c3 axis different from cn/sn axis : 0
RESULTING SCHOENFLIES SYMBOL = d2h
symmetry group of the molecule : d2h
the group has the following generators :
c2(z)
c2(x)
mirror plane sigma(xy)
8 symmetry operations found
THE MOLECULAR GEOMETRY WILL BE SYMMETRIZED !
SPECIFICATION OF MOLECULAR GEOMETRY ( #ATOMS=6 SYMMETRY=d2h )- INPUT:
ired- 内部座標を自動で設定します
START LOOP OVER SUBUNITS
time in mkbmet+rdiag cpu: 0.01 sec wall: 0.01 sec ratio: 1.0
Lowest Eigenvalue of BmBt is: 0.3939705452
Lowest Eigenvalue of BmBt-Mat 0.3939705452
time in Decoupling cpu: 0.00 sec wall: 0.00 sec ratio: 1.0
time in Preparation cpu: 0.00 sec wall: 0.00 sec ratio: 1.0
time in Loop Subunit cpu: 0.01 sec wall: 0.01 sec ratio: 1.0
time in Unpack cpu: 0.00 sec wall: 0.00 sec ratio: 1.0
time in Total Time cpu: 0.01 sec wall: 0.01 sec ratio: 1.0
SPECIFICATION OF MOLECULAR GEOMETRY ( #ATOMS=6 SYMMETRY=d2h )- INPUT:
*- 座標の指定を終わります
GEOMETRY DATA WILL BE WRITTEN TO FILE coord基底関数の設定
- 次に、計算に用いる基底関数を設定します
SUPPLYING BASIS SETS TO 6 ATOMS
ATOMIC ATTRIBUTE DEFINITION MENU ( #atoms=6 #bas=6 #ecp=0 )
b : ASSIGN ATOMIC BASIS SETS
bb : b RESTRICTED TO BASIS SET LIBRARY
bl : LIST ATOMIC BASIS SETS ASSIGNED
bm : MODIFY DEFINITION OF ATOMIC BASIS SET
bp : SWITCH BETWEEN 5d/7f AND 6d/10f
lib : SELECT BASIS SET LIBRARY
ecp : ASSIGN EFFECTIVE CORE POTENTIALS
ecpb : ecp RESTRICTED TO BASIS SET LIBRARY
ecpi : GENERAL INFORMATION ABOUT EFFECTIVE CORE POTENTIALS
ecpl : LIST EFFECTIVE CORE POTENTIALS ASSIGNED
ecprm: REMOVE EFFECTIVE CORE POTENTIAL(S)
c : ASSIGN NUCLEAR CHARGES (IF DIFFERENT FROM DEFAULTS)
cem : ASSIGN NUCLEAR CHARGES FOR EMBEDDING
m : ASSIGN ATOMIC MASSES (IF DIFFERENT FROM DEFAULTS)
iso : ASSIGN ISOTOPE FOR NUCLEAR COUPLING CALCULATION
dis : DISPLAY MOLECULAR GEOMETRY
dat : DISPLAY ATOMIC ATTRIBUTES YET ESTABLISHED
h : EXPLANATION OF ATTRIBUTE DEFINITION SYNTAX
* : TERMINATE THIS SECTION AND WRITE DATA OR DATA REFERENCES TO control
GOBACK=& (TO GEOMETRY MENU !)- INPUT:
b all 6-31G*- 基底関数として、全ての原子に 6-31G(d) を用います
SUPPLYING BASIS SETS TO 6 ATOMS
#
# BASIS SET LIBRARY FOR CARBON
# ECPs, HONDO-BASIS SETS FROM basen AND
# FULLY OPTIMIZED BASIS SETS FROM newbas MERGED 02/6/93
#
# abbreviation hondo refers to the version 7.0 of HONDO
#
########################################################################
# HF limit : E(3P) = -37.688619 a.u. (C. Froese Fischer, 1977)
########################################################################
# Roothaan parameters for C(3P) in symmetry I:
# a = 3/4 b = 3/2
########################################################################
#
c 6-31G*
c 6-31G**
# ROHF(equiv) energy is -37.67686564828 a.u. (virial theorem = 2.000362
# UHF(noneq) energy is -37.68057667956 a.u. (virial theorem = 2.000594
# P.C. Hariharan and J.A. Pople, Theor. Chim. Acta 28, 213 (1973)
# obtained from EMSL Basis Set Exchange Library 11/4/08 6:45 AM
# antique valence double zeta
#
# BASIS SET LIBRARY FOR HYDROGEN
# ECPs, HONDO-BASIS SETS FROM basen AND
# FULLY OPTIMIZED BASIS SETS FROM newbas MERGED 02/6/93
#
# abbreviation hondo refers to the version 7.0 of HONDO
#
########################################################################
# HF limit : E(2S) = -0.5 a.u.
########################################################################
# Roothaan parameters for H(2S):
# a = 0 b = 0
########################################################################
#
h 6-31G*
# HF(equiv) energy is -0.49823291073 a.u. (virial theorem = 2.023350085
# W.J. Hehre, R. Ditchfield, and J.A. Pople, J. Chem. Phys. 56, 2257 (19
# obtained from EMSL Basis Set Exchange Library 11/4/08 6:44 AM
# antique valence double zeta- INPUT:
bl- 指定した基底関数を確認します
------------------------------------------------------------------------
LIST OF BASIS SETS DEFINED YET
------------------------------------------------------------------------
INDEX | BASIS SET NICKNAME
------------------------------------------------------------------------
1 | c 6-31G*
2 | h 6-31G*
------------------------------------------------------------------------
NOTE THAT YOU MAY USE bm #<i> IF YOU WANT TO MODIFY THE
BASIS SET WITH THE INDEX <i> BY MEANS OF THE bm COMMAND- INPUT:
*- 基底関数の設定を終わります
BASIS SETS WILL BE WRITTEN TO FILE basis BY DEFAULT
+--------------------------------------------------+
| basis set information |
+--------------------------------------------------+
we will work with the 1s 3p 5d 7f 9g ... basis set
...i.e. with spherical basis functions...
type atoms prim cont basis
---------------------------------------------------------------------------
c 2 27 14 6-31G* [3s2p1d|10s4p1d]
h 4 4 2 6-31G* [2s|4s]
---------------------------------------------------------------------------
total: 6 70 36
---------------------------------------------------------------------------
total number of primitive shells : 19
total number of contracted shells : 20
total number of cartesian basis functions : 38
total number of SCF-basis functions : 36
ATOMIC COORDINATES ATOM SHELLS CHARGE PSEUDO MASS
1.24468842 0.00000000 0.00000000 c 6 6. 0 12.011
-1.24468842 0.00000000 0.00000000 c 6 6. 0 12.011
-2.31603110 -1.72798219 0.00000000 h 2 1. 0 1.008
-2.31603110 1.72798219 0.00000000 h 2 1. 0 1.008
2.31603110 1.72798219 0.00000000 h 2 1. 0 1.008
2.31603110 -1.72798219 0.00000000 h 2 1. 0 1.008
we will work with the 1s 3p 5d
SYMMETRY HAS BEEN CHANGED
there are 8 real representations : ag b1g b2g b3g au b1u b2u b3u電子配置の設定
- 次に、分子軌道を初期推定する方法や電子配置を設定します
OCCUPATION NUMBER & MOLECULAR ORBITAL DEFINITION MENU
CHOOSE COMMAND
infsao : OUTPUT SAO INFORMATION
atb : Switch for writing MOs in ASCII or binary format
eht : PROVIDE MOS && OCCUPATION NUMBERS FROM EXTENDED HUECKEL GUESS
use <file> : SUPPLY MO INFORMATION USING DATA FROM <file>
man : MANUAL SPECIFICATION OF OCCUPATION NUMBERS
hcore : HAMILTON CORE GUESS FOR MOS
flip : FLIP SPIN OF A SELECTED ATOM
& : MOVE BACK TO THE ATOMIC ATTRIBUTES MENU
THE COMMANDS use OR eht OR * OR q(uit) TERMINATE THIS MENU !!!
FOR EXPLANATIONS APPEND A QUESTION MARK (?) TO ANY COMMAND- INPUT:
eht- 拡張 Huckel 法を用いて、初期 MO の推定をおこないます
PROVIDING EHT AOS FOR THE FOLLOWING SET OF ATOMS :
1 c 2 c
for the 6 electrons of the actual atom you have
to provide at least basis functions for the AO's : 2s 1p 0d 0f
PROVIDING EHT AOS FOR THE FOLLOWING SET OF ATOMS :
3 h 4 h 5 h 6 h
for the 1 electrons of the actual atom you have
to provide at least basis functions for the AO's : 1s 0p 0d 0f
DO YOU WANT THE DEFAULT PARAMETERS FOR THE EXTENDED HUECKEL CALCULATION ?
DEFAULT=y HELP=?- INPUT: Enter or
y
ENTER THE MOLECULAR CHARGE (DEFAULT=0)- INPUT: Enter or
0- 電荷はゼロ（中性状態）
NUMBER OF ELECTRONS IN YOUR MOLECULE IS 16
AUTOMATIC OCCUPATION NUMBER ASSIGNMENT ESTABLISHED !
FOUND CLOSED SHELL SYSTEM !
HOMO/LUMO-SEPARATION : 0.232438
ORBITAL SYMMETRY ENERGY DEFAULT
(SHELL) TYPE OCCUPATION
5 1b2u -0.60494 2
6 3ag -0.51365 2
7 1b1u -0.50851 2
8 1b1g -0.49943 2
9 1b2g -0.26699 0
10 2b2u 0.42819 0
11 3b3u 0.62853 0
DO YOU ACCEPT THIS OCCUPATION ? DEFAULT=y- INPUT: Enter or
y
PROVIDING 'derivative' DEFAULT PARAMETERS ...
PROVIDING FORCE RELAXATION DEFAULT PARAMETERS ...
mo occupation :
irrep mo's occupied
ag 9 3
b1g 5 1
b2g 3 0
b3g 1 0
au 1 0
b1u 3 1
b2u 5 1
b3u 9 2
number of basis functions : 36
number of occupied orbitals : 8計算手法の選択
- 次に、計算手法を設定します
GENERAL MENU : SELECT YOUR TOPIC
scf : SELECT NON-DEFAULT SCF PARAMETER
mp2 : OPTIONS AND DATA GROUPS FOR rimp2 and mpgrad
cc : OPTIONS AND DATA GROUPS FOR ricc2
pnocc : OPTIONS AND DATA GROUPS FOR pnoccsd
ex : EXCITED STATE AND RESPONSE OPTIONS
prop : SELECT TOOLS FOR SCF-ORBITAL ANALYSIS
drv : SELECT NON-DEFAULT INPUT PARAMETER FOR EVALUATION
OF ANALYTICAL ENERGY DERIVATIVES
(GRADIENTS, FORCE CONSTANTS)
rex : SELECT OPTIONS FOR GEOMETRY UPDATES USING RELAX
stp : SELECT NON-DEFAULT STRUCTURE OPTIMIZATION PARAMETER
e : DEFINE EXTERNAL ELECTROSTATIC FIELD
dft : DFT Parameters
ri : RI Parameters
rijk : RI-JK-HF Parameters
rirpa : RIRPA Parameters
gw : OPTIONS AND DATA GROUPS FOR GW (escf)
senex : seminumeric exchange parameters
hybno : hybrid Noga/Diag parameters
dsp : DFT dispersion correction
nmr : NMR shift parameters
ncoup : NMR coupling parameters
epr : EPR parameters
trunc : USE TRUNCATED AUXBASIS DURING ITERATIONS
marij : MULTIPOLE ACCELERATED RI-J
fde : Frozen Density Embedding
dis : DISPLAY MOLECULAR GEOMETRY
list : LIST OF CONTROL FILE
& : GO BACK TO OCCUPATION/ORBITAL ASSIGNMENT MENU
* or q : END OF DEFINE SESSION- INPUT:
dft- 密度汎関数法を使います
STATUS OF DFT_OPTIONS:
DFT is NOT used
functional b-p
gridsize m3
ENTER DFT-OPTION TO BE MODIFIED
func : TO CHANGE TYPE OF FUNCTIONAL
grid : TO CHANGE GRIDSIZE
on: TO SWITCH ON DFT
Just <ENTER>, q or '*' terminate this menu.- INPUT:
func- 汎関数を設定します
SURVEY OF AVAILABLE EXCHANGE-CORRELATION ENERGY FUNCTIONALS
FUNCTIONAL | TYPE | EXCHANGE | CORRELATION | REFERENCES
---------------------------------------------------------------------
s-vwn | LDA | S | VWN(V) | 1-3
s-vwn_Gaussian | LDA | S | VWN(III) | 1-3
（省略）
hse06 | RSH |short-range PBE | |
cam-b3lyp | RSH |b(B88)+aHF |0.81LYP+0.19VWN5| 27
wb97x | RSH | from LibXC | |
wb97x-d | RSH | from LibXC | |
wb97x-v | RSH | from LibXC | |
wb97x-3c | RSH | from LibXC | |
wb97m-v | RSH | from LibXC | |
lc-wpbe | RSH | from LibXC | |
m11 | RSH | from LibXC | |
（省略）- INPUT:
wb97x-d- ωB97X-D 汎関数を選択します
STATUS OF DFT_OPTIONS:
DFT is NOT used
functional wb97x-d
gridsize m3- INPUT:
*- 基底関数の設定を終わります
設定の保存と終了
- INPUT:
*defineのセッションを終了し、設定をファイルに保存します
***********************************************************
* *
* e n d o f *
* D E F I N E *
* *
***********************************************************basis control coord ethylene ethylene.xyz mosbasis
$basis
*
c 6-31G*
# c (10s4p1d) / [3s2p1d] {631/31/1}
*
6 s
3047.5249000 0.18347000000E-02
457.36951000 0.14037300000E-01
103.94869000 0.68842600000E-01
29.210155000 0.23218440000
9.2866630000 0.46794130000
3.1639270000 0.36231200000
3 s
7.8682724000 -0.11933240000
1.8812885000 -0.16085420000
0.54424930000 1.1434564000
1 s
0.16871440000 1.0000000000
3 p
7.8682724000 0.68999100000E-01
1.8812885000 0.31642400000
0.54424930000 0.74430830000
1 p
0.16871440000 1.0000000000
1 d
0.80000000000 1.0000000000
*
h 6-31G*
# h (4s) / [2s] {31}
*
3 s
18.731137000 0.33494600000E-01
2.8253937000 0.23472695000
0.64012170000 0.81375733000
1 s
0.16127780000 1.0000000000
*
$endcontrol
$title
ethylene
$symmetry d2h
$redundant file=coord
$user-defined bonds file=coord
$coord file=coord
$optimize
internal on
redundant on
cartesian off
global off
basis off
$atoms
c 1-2 \
basis =c 6-31G*
h 3-6 \
basis =h 6-31G*
$basis file=basis
$scfmo file=mos
$closed shells
ag 1-3 ( 2 )
b1g 1 ( 2 )
b1u 1 ( 2 )
b2u 1 ( 2 )
b3u 1-2 ( 2 )
$scfiterlimit 30
$scfconv 7
$thize 0.10000000E-04
$thime 5
$scfdamp start=0.300 step=0.050 min=0.100
$scfdump
$scfintunit
unit=30 size=0 file=twoint
$scfdiis
$maxcor 500 MiB per_core
$scforbitalshift automatic=.1
$drvopt
cartesian on
basis off
global off
hessian on
dipole on
nuclear polarizability
$interconversion off
qconv=1.d-7
maxiter=25
$coordinateupdate
dqmax=0.3
interpolate on
statistics 5
$forceupdate
ahlrichs numgeo=0 mingeo=3 maxgeo=4 modus=<g|dq> dynamic fail=0.3
threig=0.005 reseig=0.005 thrbig=3.0 scale=1.00 damping=0.0
$forceinit on
diag=default
$energy file=energy
$grad file=gradient
$forceapprox file=forceapprox
$rundimensions
natoms=6
nbf(CAO)=38
nbf(AO)=36
$last step define
$endcoord
$coord natoms= 6
1.24468841888868 0.00000000000000 0.00000000000000 c
-1.24468841888868 0.00000000000000 0.00000000000000 c
-2.31603109935414 -1.72798218518704 0.00000000000000 h
-2.31603109935414 1.72798218518704 0.00000000000000 h
2.31603109935414 1.72798218518704 0.00000000000000 h
2.31603109935414 -1.72798218518704 0.00000000000000 h
$user-defined bonds
$redundant
number_of_atoms 6
degrees_of_freedom 3
internal_coordinates 3
frozen_coordinates 0
# definitions of redundant internals
1 k 1.0000000000000 stre 1 2 val= 2.48938
2 k -0.8164965809277 bend 6 5 1 val= 2.69814
0.4082482904639 bend 2 5 1
0.4082482904639 bend 2 6 1
3 k 1.0000000000000 stre 1 5 val= 2.03315
3 non zero eigenvalues of BmBt
1 4.018807994 1 0
1
2 0.916751311 2 0
2
3 0.393970545 3 0
3
$endmos
$scfmo expanded format(4d20.14)
# molecular orbitals of project :
# ---> ethylene <---
1 ag eigenvalue=-.11331129956760D+02 nsaos=9
-.99547155775189D+00-.31012736005475D-01-.99119903000270D-02-.12297534249664D-01
-.10414474590208D-010.84976680810906D-04-.21945580448892D-030.11447629849686D-01
0.21676513838550D-01
2 ag eigenvalue=-.85737982670084D+00 nsaos=9
0.21163521349115D+00-.32572717358856D+00-.42468307399236D+00-.18042015327744D-01
-.71921833271127D-02-.12090612486699D-030.88386319815307D-03-.68524462552372D-01
-.12164506647763D+00
3 ag eigenvalue=-.51365079688816D+00 nsaos=9
0.65552781497607D-01-.73044929086838D-01-.14582474645070D+000.47901078613018D+00
0.34040045867238D+000.34562250547799D-03-.37160976879735D-020.17756795596056D+00
0.32139850371599D+00
1 b1g eigenvalue=-.49943058403148D+00 nsaos=5
0.42056355721658D+000.22857163258478D+000.10229844069096D-020.30437651055092D+00
0.50965883830392D+00
1 b1u eigenvalue=-.50850575435471D+00 nsaos=3
0.54519790029590D+000.42373241876959D+000.24117758841640D-02
1 b2u eigenvalue=-.60493754850913D+00 nsaos=5
0.35185047969971D+000.23593419373656D+000.45909390681048D-02-.20198499940200D+00
-.37259482480880D+00
1 b3u eigenvalue=-.11334517389930D+02 nsaos=9
0.99189143112463D+000.43175565559743D-010.34200003097057D-01-.19796018232524D-01
-.16165856786627D-010.52236243843878D-03-.64441039437352D-030.65729510179553D-02
0.15251512581716D-01
2 b3u eigenvalue=-.70399743338664D+00 nsaos=9
-.15943982078143D+000.26462979445791D+000.28808943158560D+000.15407788265258D+00
0.12366023098259D+00-.20177475908427D-02-.97846891429849D-03-.18752997497078D+00
-.32281838229220D+00
$end計算の実行
defineでの設定後、jobexコマンドで構造最適化を実行します
jobex | tee jobex.log JOBEX
GENERAL SHELL-SCRIPT FOR GEOMETRY OPTIMIZATION IN
T U R B O M O L E
Usage: jobex <ARGUMENTS>
where the most frequently used arguments are:
----------------------------------------------------------------------
PARAMETER | ARGUMENT | DEFAULT
----------------------------------------------------------------------
view this help | -h |
| |
switch to RI versions | -ri |
(RIDFT, RDGRAD) | |
| !
switch to RIPER | -riper |
(includes periodic DFT) | |
| |
load directory | -l <directory path> |
containing executables | | $TURBODIR/bin/`sysname`
| |
select comp level | -level <arg> | -level scf
| <arg> = scf |
| -or- mp2 |
| -or- cc2 |
| -or- NumGrad |
| -or- rirpa |
| -or- uff |
| |
switch to RI-JK versions | -rijk |
of HF and CPHF | (works only with |
(RIDFT, RICC2) | '-level cc2') |
| |
select number of CPU's | -np <number> | depends on the
you'd like to use for | | machine, usually
parallel calculation | | -np 2
----------------------------------------------------------------------
and the more specialised functions are:
----------------------------------------------------------------------
PARAMETER | ARGUMENT | DEFAULT
----------------------------------------------------------------------
max. number of cycles | -c <integer> | -c 100
| |
cycle to start with | -start <integer> | 0
| |
convergency criteria (1) | |
total energy | -energy <integer> | -energy 6
cart. gradient norm | -gcart <integer> | -gcart 3
expnt. gradient norm | -gexp <integer> | -gexp 3
| |
select first step | -dscf dscf step | -dscf
| -grad gradient step | (2)
| |
program version | |
special versions | ecp -or- ncs -or- para | ' '
non-standard extension | -modus <extension> (3) |
| |
molecular dynamics | -md | (4)
| |
MD master file | -mdfile <filename> | -mdfile mdmaster
| |
MD shell script | -mdscript <filename> |
| |
excited state geometry | -ex |
optimization | |
| |
transition state | -trans | (5)
geometry optimization | |
| |
structure optimization | -relax | (6)
with the program relax | |
| |
All output directly | -v |
to terminal | |
| |
All output directed | -outfile <filename> |
to specified file | |
| |
All output files are | -keep |
preserved | |
| |
Avoid search for new | -noired |
internal redundant | |
coords if the old ones | |
fail | |
| |
Write xyz format movie | -movie |
file as the standard | jobex -movie > movie.xyz |
output | |
----------------------------------------------------------------------
(1) <integer> corresponds to a final threshold of 10**(-<integer>).
(2) If 'nextstep' file exists, the contents of this file are the first step.
(3) Non-standard filenames for DSCF and GRAD may be specified
by argument -modus <extension>, where the programs are named
'dscf_<extension>' and 'grad_<extension>': e.g to run files 'dscf_oh'
and 'grad_oh', use "-modus oh". It is assumed that there is only one
version of the relaxation programs STATPT or RELAX.
(4) Option -md uses the MD program FROG instead of STATPT to update
coordinates with gradients from 'gradient' file.
Option -mdfile looks for MD commands in <filename> instead of in mdmaster.
Option -mdscript calls a shell script <filename> before the FROG step.
(5) Transition state geometry optimization. An initial Hessian matrix must
be provided (or results of the lowest eigenvalue search). As default, the
aoforce program is run in the first cycle.
(6) Structure optimization using the relax program.
----------------------------------------------------------------------
THIS SHELL SCRIPT CAN BE CHECKED BY RUNNING jobex -check
----------------------------------------------------------------------分子研スパコンの場合
- ジョブファイル（
jobex.job）を準備して、jsubでキューイングシステムに投入します
jobex.job
#!/bin/bash
#PBS -l select=1:ncpus=8:mpiprocs=1:ompthreads=8:jobtype=core
#PBS -l walltime=24:00:00
module load turbomole
cd ${PBS_O_WORKDIR}
jobex | tee jobex.logjsub jobex.job結果の確認
- 計算の詳細は
job.lastに出力されます
job.last
OPTIMIZATION CYCLE 2
Sat May 4 20:01:54 JST 2024
STARTING grad VIA YOUR QUEUING SYSTEM!
（省略）
next step = grad
energy change : actual value = -0.4000E-09 threshold = 0.1000E-05
geom. gradient : actual value = 0.1242E-05 threshold = 0.1000E-02
CONVERGENCY CRITERIA FULFILLED IN CYCLE 2- エネルギーの変化は
energy、最適化された構造はcoordで確認できます
energy
$energy SCF SCFKIN SCFPOT
1 -78.03136095675 78.09016578910 -156.12152674585
2 -78.03136095717 78.09016339148 -156.12152434865
3 -78.03136095985 78.09016467965 -156.12152563950
$endcoord
$coord
1.24468674539171 0.00000000000000 0.00000000000000 c
-1.24468674539171 0.00000000000000 0.00000000000000 c
-2.31603595919819 -1.72798164334127 0.00000000000000 h
-2.31603595919819 1.72798164334127 0.00000000000000 h
2.31603595919819 1.72798164334127 0.00000000000000 h
2.31603595919819 -1.72798164334127 0.00000000000000 h
$user-defined bonds
$redundant
number_of_atoms 6
degrees_of_freedom 3
internal_coordinates 3
frozen_coordinates 0
# definitions of redundant internals
1 k 1.0000000000000 stre 1 2 val= 2.48937
2 k -1.1547005383792 bend 6 5 1 val= 2.69839
0.5773502691897 bend 2 5 1
0.5773502691897 bend 2 6 1
3 k 2.0000000000000 stre 1 5 val= 2.03315
3 non zero eigenvalues of BmBt
1 4.018807994 1 0
1
2 0.916751311 2 0
2
3 0.393970545 3 0
3
$endcoord形式の座標データは、t2xコマンドで、通常のxyz形式に変換できます
t2x coord | tee ethylene_opt.xyzethylene_opt.xyz
6
Energy =
C 0.6586599 0.0000000 0.0000000
C -0.6586599 0.0000000 0.0000000
H -1.2255934 -0.9144085 0.0000000
H -1.2255934 0.9144085 0.0000000
H 1.2255934 0.9144085 0.0000000
H 1.2255934 -0.9144085 0.0000000振動数解析
- 最適化された立体構造を用いて、分子の振動状態を解析します
- 振動数解析をおこなうことで、赤外 or ラマンスペクトルを予測することができます
計算の実行
- 構造最適化を実行したときと同じディレクトリで
aoforceコマンドを使用して振動数解析を実行します
aoforce | tee force.out分子研スパコンの場合
- ジョブファイル（
aoforce.job）を準備して、jsubでキューイングシステムに投入します
jobex.job
#!/bin/bash
#PBS -l select=1:ncpus=8:mpiprocs=1:ompthreads=8:jobtype=core
#PBS -l walltime=24:00:00
module load turbomole
cd ${PBS_O_WORKDIR}
aoforce | tee force.outjsub aoforce.job結果の確認
- 計算の詳細は
force.outに出力されます
force.out
（省略）
mode 7 8 9 10 11 12
frequency 897.15 1094.43 1098.32 1154.55 1352.75 1496.52
symmetry b2u b1u b2g au b1g ag
IR YES YES NO NO NO NO
|dDIP/dQ| (a.u.) 0.0003 0.0075 0.0000 0.0000 0.0000 0.0000
intensity (km/mol) 0.22 100.73 0.00 0.00 0.00 0.00
intensity ( % ) 0.21 100.00 0.00 0.00 0.00 0.00
（省略）- 振動スペクトルは
vibspectrumファイルで確認できます
vibspectrum
$vibrational spectrum
# mode symmetry wave number IR intensity selection rules
# cm**(-1) km/mol IR RAMAN
1 -0.00 0.00000 - -
2 -0.00 0.00000 - -
3 0.00 0.00000 - -
4 0.00 0.00000 - -
5 0.00 0.00000 - -
6 0.00 0.00000 - -
7 b2u 897.15 0.21646 YES NO
8 b1u 1094.43 100.72648 YES NO
9 b2g 1098.32 0.00000 NO YES
10 au 1154.55 0.00000 NO NO
11 b1g 1352.75 0.00000 NO YES
12 ag 1496.52 0.00000 NO YES
13 b3u 1611.16 5.94160 YES NO
14 ag 1856.79 0.00000 NO YES
15 b3u 3318.48 25.16488 YES NO
16 ag 3342.06 0.00000 NO YES
17 b1g 3392.88 0.00000 NO YES
18 b2u 3419.06 40.17169 YES NO
$end励起状態計算
- 励起状態計算では、電子励起状態とそれに関連する特性を調べることができます
define セッションの開始
defineを実行して設定セッションを開始します
define***********************************************************
* *
* D E F I N E *
* *
* TURBOMOLE'S INTERACTIVE INPUT PROGRAM *
* *
* Quantum Chemistry Group University of Karlsruhe *
* *
***********************************************************
FILE control ALREADY EXISTS
I WILL PLUG IN THE NEW DATA.
DATA WILL BE TAKEN FROM control BY DEFAULT
DEFAULT TITLE OF THIS PROJECT IS :
ethylene
HIT >return< TO ACCEPT DEFAULT TITLE OR
INPUT TITLE- INPUT: Enter
SYMMETRY d2h AND CARTESIAN COORDINATES FOR 6 ATOMS
HAVE BEEN READ FROM THE DEFAULT INPUT FILE control .
DEFINITIONS OF INTERNAL COORDINATES HAVE N O T BEEN READ.
SPECIFICATION OF BOND TOPOLOGY HAS BEEN READ.
DO YOU WANT TO CHANGE THE GEOMETRY DATA ? DEFAULT=n GOBACK=&- INPUT: Enter
ATOMIC ATTRIBUTE DATA (BASES,CHARGES,MASSES,ECPS) HAVE BEEN
TAKEN FROM THE DEFAULT INPUT FILE control.
DO YOU WANT TO CHANGE THESE DATA ? DEFAULT=n GOBACK=&- INPUT: Enter
+--------------------------------------------------+
| basis set information |
+--------------------------------------------------+
we will work with the 1s 3p 5d 7f 9g ... basis set
...i.e. with spherical basis functions...
type atoms prim cont basis
---------------------------------------------------------------------------
c 2 27 14 6-31G* [3s2p1d|10s4p1d]
h 4 4 2 6-31G* [2s|4s]
---------------------------------------------------------------------------
total: 6 70 36
---------------------------------------------------------------------------
total number of primitive shells : 19
total number of contracted shells : 20
total number of cartesian basis functions : 38
total number of SCF-basis functions : 36
ATOMIC COORDINATES ATOM SHELLS CHARGE PSEUDO MASS
1.24468675 0.00000000 0.00000000 c 6 6. 0 12.011
-1.24468675 0.00000000 0.00000000 c 6 6. 0 12.011
-2.31603596 -1.72798164 0.00000000 h 2 1. 0 1.008
-2.31603596 1.72798164 0.00000000 h 2 1. 0 1.008
2.31603596 1.72798164 0.00000000 h 2 1. 0 1.008
2.31603596 -1.72798164 0.00000000 h 2 1. 0 1.008
we will work with the 1s 3p 5d
there are 8 real representations : ag b1g b2g b3g au b1u b2u b3u
mo occupation :
irrep mo's occupied
ag 9 3
b1g 5 1
b2g 3 0
b3g 1 0
au 1 0
b1u 3 1
b2u 5 1
b3u 9 2
number of basis functions : 36
number of occupied orbitals : 8
reading orbital data $scfmo from file mos
orbital characterization : scfconv=7
MOLECULAR ORBITAL DATA (OCCUPATION NUMBERS,MOS) HAVE BEEN
TAKEN FROM THE DEFAULT INPUT FILE control .
DO YOU WANT TO CHANGE THESE DATA ? DEFAULT=n GOBACK=&- INPUT: Enter
orbitals $scfmo will be written to file mos
ADJUSTING DATA GROUP $drvopt FOR CURRENT REQUIREMENTS
FORCE CONSTANT INITIALIZATION $forceinit WILL BE ENABLED
DO YOU WANT TO DELETE DATA GROUPS LIKE
$energy
$grad
$hessian
$hessian (projected)
$dipole
$dipgrad
$vibrational normal modes
$vibrational spectrum
LEFT OVER FROM PREVIOUS CALCULATIONS ? DEFAULT(n)- INPUT: Enter
mo occupation :
irrep mo's occupied
ag 9 3
b1g 5 1
b2g 3 0
b3g 1 0
au 1 0
b1u 3 1
b2u 5 1
b3u 9 2
number of basis functions : 36
number of occupied orbitals : 8励起状態計算の設定
- 次に、励起状態の計算条件を設定します
GENERAL MENU : SELECT YOUR TOPIC
scf : SELECT NON-DEFAULT SCF PARAMETER
mp2 : OPTIONS AND DATA GROUPS FOR rimp2 and mpgrad
cc : OPTIONS AND DATA GROUPS FOR ricc2
pnocc : OPTIONS AND DATA GROUPS FOR pnoccsd
ex : EXCITED STATE AND RESPONSE OPTIONS
prop : SELECT TOOLS FOR SCF-ORBITAL ANALYSIS
drv : SELECT NON-DEFAULT INPUT PARAMETER FOR EVALUATION
OF ANALYTICAL ENERGY DERIVATIVES
(GRADIENTS, FORCE CONSTANTS)
rex : SELECT OPTIONS FOR GEOMETRY UPDATES USING RELAX
stp : SELECT NON-DEFAULT STRUCTURE OPTIMIZATION PARAMETER
e : DEFINE EXTERNAL ELECTROSTATIC FIELD
dft : DFT Parameters
ri : RI Parameters
rijk : RI-JK-HF Parameters
rirpa : RIRPA Parameters
gw : OPTIONS AND DATA GROUPS FOR GW (escf)
senex : seminumeric exchange parameters
hybno : hybrid Noga/Diag parameters
dsp : DFT dispersion correction
nmr : NMR shift parameters
ncoup : NMR coupling parameters
epr : EPR parameters
trunc : USE TRUNCATED AUXBASIS DURING ITERATIONS
marij : MULTIPOLE ACCELERATED RI-J
fde : Frozen Density Embedding
dis : DISPLAY MOLECULAR GEOMETRY
list : LIST OF CONTROL FILE
& : GO BACK TO OCCUPATION/ORBITAL ASSIGNMENT MENU
* or q : END OF DEFINE SESSION- INPUT:
ex
MAIN MENU FOR RESPONSE CALCULATIONS
OPTION | STATUS | DESCRIPTION
-------------------------------------------------------------------
rpas | off | RPA SINGLET EXCITATIONS (TDHF OR TDDFT)
ciss | off | TDA SINGLET EXCITATIONS (CI SINGLES)
rpat | off | RPA TRIPLET EXCITATIONS (TDHF OR TDDFT)
cist | off | TDA TRIPLET EXCITATIONS (CI SINGLES)
polly | off | STATIC POLARIZABILITY
dynpol | off | DYNAMIC POLARIZABILITY
single | off | SINGLET STABILITY ANALYSIS
triple | off | TRIPLET STABILITY ANALYSIS
nonrel | off | NON-REAL STABILITY ANALYSIS
bse | off | BETHE-SALPETER EX.
cbse | off | corr-aug. BETHE-SALPETER EX.
ENTER <OPTION> TO SWITCH ON/OFF OPTION, * OR q TO QUIT- INPUT:
rpas- TD-DFT を選択します
MAIN MENU FOR RESPONSE CALCULATIONS
OPTION | STATUS | DESCRIPTION
-------------------------------------------------------------------
rpas | on | RPA SINGLET EXCITATIONS (TDHF OR TDDFT)
ciss | off | TDA SINGLET EXCITATIONS (CI SINGLES)
rpat | off | RPA TRIPLET EXCITATIONS (TDHF OR TDDFT)
cist | off | TDA TRIPLET EXCITATIONS (CI SINGLES)
polly | off | STATIC POLARIZABILITY
dynpol | off | DYNAMIC POLARIZABILITY
single | off | SINGLET STABILITY ANALYSIS
triple | off | TRIPLET STABILITY ANALYSIS
nonrel | off | NON-REAL STABILITY ANALYSIS
bse | off | BETHE-SALPETER EX.
cbse | off | corr-aug. BETHE-SALPETER EX.
ENTER <OPTION> TO SWITCH ON/OFF OPTION, * OR q TO QUIT- INPUT:
*- 次の設定に進みます
STATE SELECTION MENU
IRREP | #STATES | #SELECTED |
----------------------------------------------------------------
ag | 42 | 0 |
b1g | 34 | 0 | Lz
b2g | 22 | 0 | Ly
b3g | 14 | 0 | Lx
au | 14 | 0 |
b1u | 20 | 0 | z
b2u | 34 | 0 | y
b3u | 44 | 0 | x
SELECT IRREP AND NUMBER OF STATES
ENTER ? FOR HELP, * OR Q TO QUIT, & TO GO BACK- INPUT:
all 3- 各対称性について最初の３つの励起状態を計算してみます
STATE SELECTION MENU
IRREP | #STATES | #SELECTED |
----------------------------------------------------------------
ag | 42 | 3 |
b1g | 34 | 3 | Lz
b2g | 22 | 3 | Ly
b3g | 14 | 3 | Lx
au | 14 | 3 |
b1u | 20 | 3 | z
b2u | 34 | 3 | y
b3u | 44 | 3 | x
SELECT IRREP AND NUMBER OF STATES
ENTER ? FOR HELP, * OR Q TO QUIT, & TO GO BACK- INPUT:
*- 次の設定に進みます
GENERAL OPTIONS FOR RESPONSE CALCULATIONS
OPTION | VALUE | DESCRIPTION
-------------------------------------------------------
rpacor | 200 | MEMORY (MB) (RECOMMENDED: 200)
rpaconv | 5 | CONVERGENCE THRESHOLD (DEFAULT: 5)
-------------------------------------------------------
RECOMMENDED MEMORY CAN BE DIVIDED BY 2 FOR non-hybrid DFT!
ENTER <OPTION> <VALUE> TO CHANGE OPTION, * OR q TO QUIT- INPUT:
*- 次の設定に進みます
SET SCF DENSITY CONVERGENCE THRESHOLD $denconv TO 1d-7
(RECOMMENDED FOR RESPONSE CALCULATIONS) (y/n, DEFAULT:y)?- INPUT: Enter
- デフォルトの設定のままにします
ジョブの保存と終了
- INPUT:
*- 設定を終了します
***********************************************************
* *
* e n d o f *
* D E F I N E *
* *
***********************************************************control
$title
ethylene
$symmetry d2h
$redundant file=coord
$user-defined bonds file=coord
$coord file=coord
$optimize
internal on
redundant on
cartesian off
global off
basis off
$atoms
c 1-2 \
basis =c 6-31G*
h 3-6 \
basis =h 6-31G*
$basis file=basis
$scfmo file=mos
$scfiterlimit 30
$scfconv 7
$thize 0.10000000E-04
$thime 5
$scfdamp start=0.300 step=0.050 min=0.100
$scfdump
$scfintunit
unit=30 size=0 file=twoint
$scfdiis
$scforbitalshift automatic=.1
$drvopt
cartesian on
basis off
global off
hessian on
dipole on
nuclear polarizability
$interconversion off
qconv=1.d-7
maxiter=25
$coordinateupdate
dqmax=0.3
interpolate on
statistics 5
$forceupdate
ahlrichs numgeo=0 mingeo=3 maxgeo=4 modus=<g|dq> dynamic fail=0.3
threig=0.005 reseig=0.005 thrbig=3.0 scale=1.00 damping=0.0
$forceinit on
diag=default
$energy file=energy
$grad file=gradient
$forceapprox file=forceapprox
$rundimensions
natoms=6
nbf(CAO)=38
nbf(AO)=36
$last SCF energy change = -.26793003E-08
$charge from dscf
-0.000 (not to be modified here)
$dipole from force
x -0.00000000000000 y 0.00000000000000 z 0.00000000000000 a.u.
| dipole | = 0.0000000000 debye
$optinfo file=optinfo
$hessapprox file=hessapprox
$orbital_max_rnorm 0.10769779602030E-06
$dipgrad file=dipgrad
$hessian (projected) file=hessian
$vibrational normal modes file=vib_normal_modes
$vibrational reduced masses
1.3258967772 1.8945238567 3.6287625946 2.0471565005 4.2506770105
3.5143354754 1.0426731556 1.1607969330 1.5196113978 1.0079700000
1.5224873469 1.2529577899 1.1117148104 2.9267262232 1.0481533486
1.0779906478 1.1164337846 1.1179469536
$nvibro 18
$vibrational spectrum file=vibspectrum
$closed shells
ag 1-3 ( 2 )
b1g 1 ( 2 )
b1u 1 ( 2 )
b2u 1 ( 2 )
b3u 1-2 ( 2 )
$scfinstab rpas
$soes all 3
#$maxcor 500 MiB per_core
$denconv 1d-7
$last step define
$end計算の実行
dscfコマンドで一点計算をおこなったあと、escfコマンドで励起状態を計算します
dscf | tee dscf.logescf | tee escf.log分子研スパコンの場合
- ジョブファイル（
escf.job）を準備して、jsubでキューイングシステムに投入します
jobex.job
#!/bin/bash
#PBS -l select=1:ncpus=8:mpiprocs=1:ompthreads=8:jobtype=core
#PBS -l walltime=24:00:00
module load turbomole
cd ${PBS_O_WORKDIR}
dscf | tee dscf.log
escf | tee escf.logjsub escf.job計算の確認
- 計算の詳細は
escf.logに出力されます - 励起スペクトルは
exspectrumで確認することができます
exspectrum
# Excitation spectrum ethylene; written at 2024-05-04 21:13:05
# Exc. energy (Eh) energy (eV) energy (cm-1) energy (nm) Osc.(vel) Osc.(len)
1 ag 0.589441 16.03950 0.12936731D+06 77.299 0.00000000 0.00000000
2 ag 0.643592 17.51302 0.14125205D+06 70.795 0.00000000 0.00000000
3 ag 0.706412 19.22245 0.15503951D+06 64.500 0.00000000 0.00000000
1 b1g 0.494028 13.44318 0.10842659D+06 92.228 0.00000000 0.00000000
2 b1g 0.681824 18.55339 0.14964315D+06 66.826 0.00000000 0.00000000
3 b1g 0.708695 19.28457 0.15554056D+06 64.292 0.00000000 0.00000000
1 b2g 0.388259 10.56506 0.85212985D+05 117.353 0.00000000 0.00000000
2 b2g 0.433395 11.79328 0.95119227D+05 105.131 0.00000000 0.00000000
3 b2g 0.570503 15.52418 0.12521096D+06 79.865 0.00000000 0.00000000
1 b3g 0.356027 9.68798 0.78138813D+05 127.977 0.00000000 0.00000000
2 b3g 0.378958 10.31198 0.83171732D+05 120.233 0.00000000 0.00000000
3 b3g 0.809799 22.03577 0.17773044D+06 56.265 0.00000000 0.00000000
1 au 0.479360 13.04405 0.10520737D+06 95.050 0.00000000 0.00000000
2 au 0.507747 13.81651 0.11143768D+06 89.736 0.00000000 0.00000000
3 au 0.929827 25.30187 0.20407338D+06 49.002 0.00000000 0.00000000
1 b1u 0.377417 10.27005 0.82833559D+05 120.724 0.00326119 0.00074012
2 b1u 0.664826 18.09085 0.14591251D+06 68.534 0.13719402 0.13363109
3 b1u 0.712697 19.39348 0.15641898D+06 63.931 0.16650119 0.19466603
1 b2u 0.544350 14.81251 0.11947098D+06 83.702 1.09170693 1.06128111
2 b2u 0.584721 15.91108 0.12833151D+06 77.923 0.11403575 0.10592582
3 b2u 0.671742 18.27903 0.14743033D+06 67.829 0.43679806 0.47637986
1 b3u 0.304073 8.27424 0.66736232D+05 149.844 0.44155754 0.45235991
2 b3u 0.534956 14.55690 0.11740934D+06 85.172 0.73906903 0.73111966
3 b3u 0.622605 16.94194 0.13664596D+06 73.182 0.34006115 0.32566921さいごに
defineで対話的にインプットファイルを準備しながら計算に取り組んでいくという、他の量子化学計算パッケージとは異なるスタイルには慣れが必要かも知れません