CONFLEX Tutorials

動的反応座標（Dynamic Reaction Coordinate,DRC）法による動力学計算

配座同士の変換や、分子間相互作用で配位している分子の移動は、ポテンシャルエネルギー曲面上において低周波数の固有振動モードに沿って構造が変化することで起こると考えられます（下図）。
ここでは、基準振動解析により求めた振動数と振動モードを利用して動力学計算を行う、動的反応座標（DRC）法の概要と計算例を示します。

Potential Energy
ポテンシャルエネルギー曲面と分子振動の模式図

【DRC法の概要】

熱平衡状態において、エネルギーが各基準振動モードに等分配されるとしますと、i番目の基準振動モードの振幅αi（振動による平衡構造からのずれ）は Amplitude で表されます。ここでRは気体定数、Tは絶対温度、cは光速、ωiはi番目の基準振動数です。
時間tでの熱揺らぎによる平衡構造からのずれΔx(t)は、以下の式で求めることができます。 Delta X ここでqiは基準振動モード、δiは振動モードqiの位相差で、デフォルトでは-90°にしています。

この式を使って、まず平衡構造からの変位ベクトルと初速を求めるのですが、この時、振幅αiは、下記の式を用いて振動数の総和と運動エネルギーを元に求めた値Factorでスケーリングします。 Factor ここでmjは原子jの質量、vijはi番目の振動モードでの原子jの速度を表します。なお、Factor（スケーリングファクター）は、振動数の和の代わりに別の値を指定（DRC_KICK=）することも、また値そのものを指定（DRC_SCALE=）することもできます。

平衡構造から動いた後は、その構造でのグラジエント値から変位と速度を求めることで、トラジェクトリーと各点でのエネルギー値を求めます。

【DRC計算例：α-D-Glucose + H2O】

下図に示した、α-D-Glucoseに水分子が一つ配位した系に、DRC計算を適用します。図の赤い矢印は最低振動モード（振動数=30.09 cm-1）を表していて、このモードを使って計算を行います。

α-D-Glucose + H2Oの極小構造（E=74.7739 kcal/mol）と最低振動モード（振動数=30.09 cm-1）
DRC Fig.2

α-D-Glucose + H2Oの座標データ（aDglucose_H2O.mol）

aDglucose_H2O.mol


 27 26  0  0  0               999 V2000
   -0.5024    2.4148    0.6016 O   0  0  0  0  0  0  0  0  0  0  0  0
   -0.8520   -1.0932   -0.5731 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.4519   -1.5744   -0.2774 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.2020    1.2411   -0.1553 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1949    0.1125    0.1452 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.1431    1.8092   -0.2200 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.5330   -1.9764    1.0847 O   0  0  0  0  0  0  0  0  0  0  0  0
    2.8319   -0.9820   -0.1321 O   0  0  0  0  0  0  0  0  0  0  0  0
   -2.6186    0.5024   -0.2701 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.5391   -0.5233   -0.5592 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.2184    0.7799    0.1723 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.5094   -0.5881   -0.0183 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.2794    2.9953    0.4978 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6445   -2.4621   -0.8896 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.2535    1.5201   -1.2154 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.2141   -0.1177    1.2184 H   0  0  0  0  0  0  0  0  0  0  0  0
    3.0267    1.3889   -0.1951 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0183   -2.7768    1.1469 H   0  0  0  0  0  0  0  0  0  0  0  0
    2.6839   -1.4193    0.7318 H   0  0  0  0  0  0  0  0  0  0  0  0
   -2.6702    0.7196   -1.3421 H   0  0  0  0  0  0  0  0  0  0  0  0
   -2.9791    1.3731    0.2851 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.6027   -0.3262   -1.6358 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.3472    0.6668    1.2558 H   0  0  0  0  0  0  0  0  0  0  0  0
   -3.0727   -1.3725   -0.4004 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.3106    0.8234   -2.7522 O   0  0  0  0  0  0  0  0  0  0  0  0
   -0.3290    0.0085   -2.2451 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.9697    0.7822   -3.4490 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  4  1  0  0  0  0
  1 13  1  0  0  0  0
  2  3  1  0  0  0  0
  2  5  1  0  0  0  0
  3  7  1  0  0  0  0
  3 10  1  0  0  0  0
  3 14  1  0  0  0  0
  4  5  1  0  0  0  0
  4 11  1  0  0  0  0
  4 15  1  0  0  0  0
  5  9  1  0  0  0  0
  5 16  1  0  0  0  0
  6 11  1  0  0  0  0
  6 17  1  0  0  0  0
  7 18  1  0  0  0  0
  8 10  1  0  0  0  0
  8 19  1  0  0  0  0
  9 12  1  0  0  0  0
  9 20  1  0  0  0  0
  9 21  1  0  0  0  0
 10 11  1  0  0  0  0
 10 22  1  0  0  0  0
 11 23  1  0  0  0  0
 12 24  1  0  0  0  0
 25 26  1  0  0  0  0
 25 27  1  0  0  0  0
M  END

[Interfaceから実行する場合]

aDglucose_H2O.molファイルをCONFLEX Interfaceを用いて開きます。 Interface Glucose H2O

Calculationメニューから「CONFLEX」を選択し、開いた計算設定ダイアログのDetail Settingsをクリックします。詳細設定ダイアログが開きます。 Basic Settings

まず、詳細設定ダイアログの「General Settings」ダイアログにある「Calculation Type:」のプルダウンメニューから「Dynamic Reaction Coordinate」を選択します。 General Settings

次に、詳細設定ダイアログの「Dynamic Reaction Coordinate」ダイアログでDRC計算の設定を行います。
最低振動モードを利用するので、「Mode:」は「1」であることを確認してください。「Snapshot:」は「10」とします。これは、何ステップ毎に構造を出力させるかを指定しています。「Kick:」は「3.225」とします。これは、Factorの計算式における分子（振動数の総和）を、指定した値（3.225 kcal/mol）に変更します。 DRC Settings 設定が終わりましたら、詳細設定ダイアログのSubmitをクリックします。DRC計算が行われます。

[コマンドラインから実行する場合]

計算設定は、aDglucose_H2O.iniファイルにキーワードを記述することで行います。

DMB.iniファイル

DRC
SNAPSHOT=10
DRC_SMODE=1
DRC_NSTEP=0.0002
DRC_KICK=3.225

「DRC」キーワードは、DRC計算を行うことを意味します。
「SNAPSHOT=」は何ステップ毎に構造を出力させるかの指定、「DRC_SMODE=」はDRC計算に用いる振動モードが一つの場合の番号の指定、「DRC_NSTEP=」は時間刻み幅（単位はps）、「DRC_KICK=」は振幅のスケーリングファクターを計算する際の振動エネルギーの指定（単位はkcal/mol）をそれぞれ表しています。

aDglucose_H2O.molとaDglucose_H2O.iniをフォルダに格納し、下記コマンドを実行してください。計算が始まります。 

C:\CONFLEX\bin\flex9a_win_x64.exe   -par   C:\CONFLEX\par   aDglucose_H2Oenter

上記は、Windowsの場合です。他のOSにおける実行コマンドについては、本書の「実行方法」を参照してください。

計算結果

計算を実行しますと、時間経過に伴うエネルギー値の変化がaDglucose_H2O.bsoに、また10 step毎の構造の変化がaDglucose_H2O-D.sdfにそれぞれ出力されます。
構造とエネルギー値の変化を下図に示します。2 psで、初期構造より安定な配位構造が得られていることがわかります。 α-D-Glucose + H2OのDRC計算によるトラジェクトリーとエネルギー変化

aDglucose_H2O.bso :

 !=====================================================================================================================!
 !                                                                                                                     !
 ! DYNAMIC REACTION COORDINATES(Simple Molecular Dynamics): DRC                                                        !
 !                                                                                                                     !
 !---------------------------------------------------------------------------------------------------------------------!
 !                                                                                                                     !
 ! SNAPSHOT PRINTING FOR EACH N-TH STEP:        10                                                                     !
 !    SCALING FACTOR FOR ALL VIBRATIONS:         2.3331                                                                !
 !  TEMPERATURE OF VIBRATIONAL DYNAMICS:        25.00  [DEGREE CELSIUS]                                                !
 !  STARTING TIME OF DYNAMIC SIMULATION:         0.0000 [PS]                                                           !
 !  TERMINAL TIME OF DYNAMIC SIMULATION:         4.4347 [PS]                                                           !
 !  INTERVAL TIME OF DYNAMIC SIMULATION:         0.0002 [PS]                                                           !
 !     SINGLE VIBRATIONAL MODE SELECTED:                                                                               !
 !                          NO.     WAVE NUMBER     PHASE        PERIOD          VIB ENERGY      KINECTIC              !
 !                                  [CM**-1]        [DEGREE]     [PS]            [KCAL/MOL]      [KCAL/MOL]            !
 !                           1         30.09      -90.0000        1.1087         8.60222E-02     0.59249               !
 !                                                                                                                     !
 ! *WARNING: DRC_PHASE - INITIALIZATION FAILED, THEN ALL PHASES ARE SET TO ZERO.                                       !
 !                                                                                                                     !
 !=====================================================================================================================!


 ENERGY PROFILE:

     NSTEP     TIME        ENERGY              DELTA E          GRMS             XRMS             VMAX            ATOM     KINETIC           TOTAL E  
               [PS]       [KCAL/MOL]          [KCAL/MOL]       [KCAL/MOL/ANGS]  [ANGS]           [ANGS]                   [KCAL/MOL]        [KCAL/MOL]
         0     0.0000   74.773940376484         0.000000        2.5863909E-07
         1     0.0002   74.773944532084        4.1555993E-06    4.0413004E-04    4.3581107E-04    2.0981112E-03             3.225000          3.225004    
         2     0.0004   74.773957054786        1.6678301E-05    1.0907419E-03    2.4898240E-09    2.1673713E-08     26      3.224996          3.225013    
         3     0.0006   74.773978042534        3.7666050E-05    2.1430837E-03    9.7631357E-09    6.5026510E-08     26      3.224983          3.225021    
         4     0.0008   74.774007609620        6.7233136E-05    3.5643141E-03    2.1742388E-08    1.0712481E-07     26      3.224962          3.225030    
         5     0.0010   74.774045882484        1.0550600E-04    5.3337169E-03    3.8194172E-08    1.4712886E-07     26      3.224933          3.225038    
         6     0.0012   74.774092994323        1.5261784E-04    7.4197651E-03    5.8797683E-08    1.8425723E-07     26      3.224894          3.225047    
         7     0.0014   74.774149078954        2.0870247E-04    9.7836141E-03    8.3151274E-08    2.1778532E-07     26      3.224847          3.225056    
         8     0.0016   74.774214264425        2.7388794E-04    1.2380669E-02    1.1078037E-07    2.4705933E-07     26      3.224791          3.225065    
         9     0.0018   74.774288666916        3.4829043E-04    1.5161733E-02    1.4114677E-07    2.7150898E-07     26      3.224726          3.225074    
        10     0.0020   74.774372385439        4.3200895E-04    1.8074090E-02    1.7365925E-07    2.9065857E-07     26      3.224652          3.225084    
  
・・・・・