Dynamics simulation by dynamic reaction coordinate (DRC) method
Low-frequency molecular vibrations on a potential energy surface are thought to provide conformational conversions and intermolecular movements (Figure 1). Here, we introduce dynamics simulations by dynamic reaction coordinate (DRC) method that uses frequencies and vibrational modes obtained by a normal mode analysis.
[Overview of DRC method]
When an energy is equally distributed to each normal mode in a thermal equilibrium state, the amplitude of i-th mode (αi) is represented as
where the R is gas constant, the T is absolute temperature, the c is speed of light, and the ωi is frequency of i-th normal mode.
The deviation from the equilibrium structure due to thermal fluctuations at time t (Δx(t)) is expressed as
where the qi is i-th normal mode and the δi is phase difference of i-th mode (-90 degrees as default). This equation is used to calculate displacement vectors and initial velocities of each atom in the equilibrium structure. The amplitude αi is scaled by the Factor value calculated by the sum of frequencies and the sum of kinetic energies.
where the mj is mass of atom j, the vij is velocity of atom j due to the i-th normal mode. In the calculation of Factor (scaling factor) value, you can specify another value by using [DRC_KICK=] keyword instead of the sum of frequencies. Moreover, you can directly specify the Factor value by using [DRC_SCALE=] keyword.
After each atom in the equilibrium structure moves according to the displacement vectors, the displacements and velocities are calculated by the gradients on each atom.
[Example calculation: α-D-Glucose + H2O]
We apply DRC method to a system consists of α-D-Glucose and water molecules. The red arrows in the below figure show the lowest normal mode (frequency: 30.09 cm-1), and this mode is used for the DRC calculation.
Structure data of α-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
[Execution by Interface]
Open the aDglucose_H2O.mol file by CONFLEX Interface.
Select [CONFLEX] in Calculation menu, and click in the calculation setting dialog displayed. A detail setting dialog will be displayed.
First, in [General Settings] dialog on the detail setting dialog, select [Dynamic Reaction Coordinate] in the pull-down menu of [Calculation Type:].
Settings of the DRC calculation are made in [Dynamic Reaction Coordinate] dialog.
Make sure [Mode:] is [1] because we use the lowest normal mode. Set [Snapshot:] to 10, and the coordinate data is outputted every 10 steps. Set [Kick:] to 3.225, and the numerator in the equation of the Factor changes from the sum of frequencies to 3.225 kcal/mol.
When the calculation settings are complete, click . The calculation will start.
[Execution by command line]
The calculation settings are defined by describing keywords in the aDglucose_H2O.ini file.
DMB.ini file
DRC SNAPSHOT=10 DRC_SMODE=1 DRC_NSTEP=0.0002 DRC_KICK=3.225
[DRC] keyword means that to preform a DRC calculation.
By [SNAPSHOT=10], the coordinate data is outputted every 10 steps.
By [DRC_SMODE=1], the lowest normal mode is applied to the DRC calculation.
By [DRC_NSTEP=0.0002], the time step is set to 0.0002 ps.
By [DRC_KICK=3.225], the numerator in the equation of the Factor changes from the sum of frequencies to 3.225 kcal/mol.
Store the two files of aDglucose_H2O.mol and aDglucose_H2O.ini in an one folder, and execute below command. The calculation will start.
C:\CONFLEX\bin\flex9a_win_x64.exe -par C:\CONFLEX\par aDglucose_H2Oenter
The above command is for Windows OS. For the other OS, please refer to [How to execute CONFLEX].
Calculation results
When the calculation finished, the energies on each step are outputted to aDglucose_H2O.bso, and the structures every 10 steps are printed to aDglucose_H2O-D.sdf.
The figure below shows the structure and energy changes on DRC trajectory. You can see that the structure at 2 ps is more stable than the initial one (E=74.7739 kcal/mol).
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
・・・・・