Dynamics simulation using the dynamic reaction coordinate (DRC) method
Low-frequency molecular vibrations on a potential energy surface are thought to facilitate conformational changes and intermolecular motions (Figure 1). Here, we introduce dynamic simulations using the dynamic reaction coordinate (DRC) method, which incorporates frequencies and vibrational modes obtained from normal mode analysis.
[Overview of DRC method]
When energy is equally distributed among all normal modes in a state of thermal equilibrium, the amplitude of the i-th mode (αi) can be expressed as follows:
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 by 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 the Factor (scaling factor) value, you can specify an alternative value using the [DRC_KICK=] keyword instead of using the sum of frequencies. Additionally, the Factor value can be directly specified using the [DRC_SCALE=] keyword.
After each atom in the equilibrium structure is displaced according to the displacement vectors, the subsequent displacements and velocities are calculated based on the gradients acting on each atom.
[Example calculation: α-D-Glucose + H2O]
We apply the DRC method to a system consisting of α-D-glucose and water molecules. The red arrows in the figure below indicate the lowest normal mode (frequency: 30.09 cm-1), which 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 from Interface]
Open the aDglucose_H2O.mol file using CONFLEX Interface.
Select [CONFLEX] from the Calculation menu, and then click in the calculation setting dialog that appears. A detailed settings dialog will be displayed.
First, in [General Settings] dialog of the detailed settings dialog, select [Dynamic Reaction Coordinate] from the [Calculation Type:] pull-down menu.
The settings for the DRC calculation are configured in the [Dynamic Reaction Coordinate] dialog.
Make sure that [Mode:] is set to [1], as we are using the lowest normal mode.
Set [Snapshot:] to 10 so that coordinate data is output every 10 steps.
Set [Kick:] to 3.225, which changes the numerator in the Factor equation from the sum of frequencies to 3.225 kcal/mol.
After completing the calculation settings, click to start the calculation.
[Execution from command line]
The calculation settings are defined by specifying 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 specifies that a DRC calculation should be performed.
[SNAPSHOT=10] outputs coordinate data every 10 steps.
[DRC_SMODE=1] applies the lowest normal mode to the DRC calculation.
[DRC_NSTEP=0.0002] sets the time step to 0.0002 ps.
[DRC_KICK=3.225] changes the numerator in the Factor equation from the sum of frequencies to 3.225 kcal/mol.
Store the two files of aDglucose_H2O.mol and aDglucose_H2O.ini in a single folder, and execute the following command to start the calculation.
C:\CONFLEX\bin\conflex-10a.exe -par C:\CONFLEX\par aDglucose_H2Oenter
The command above is for Windows OS. For other OS, please refer to [How to execute CONFLEX].
Calculation results
When the calculation is finished, the energies at each step are output to aDglucose_H2O.bso, and the structures at every 10 steps are saved in aDglucose_H2O-D.sdf.
The figure below shows the structural and energy changes along the DRC trajectory.
It can be seen 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, and 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
・・・・・