.. tip:: All input files can be downloaded: :download:`Files `. md ====== .. contents:: :local: This keyword controls the molecular dynamics (MD) simulation details. Options ------------ .. option:: type .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - ``NVE`` Microcanonical ensemble, constant number of particles, volume, and energy * - - ``NVT`` Canonical ensemble, constant number of particles, volume, and temperature * - - ``fep-NVE`` Free energy perturbation calculation for microcanonical ensemble * - - ``fep-NVT`` Free energy perturbation calculation for canonical ensemble * - Default - ``NVE`` Define the MD simulation type. .. option:: dt .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``0.0005`` The time step in ps. For systems not containing hydrogen atoms, ``dt`` can be set to ``0.001`` ps. .. option:: num_steps .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - An integer * - Default - ``10`` The number of MD steps to be run. The total simulation time is ``dt``*\ ``num_steps`` ps. .. option:: output_freq .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - An integer * - Default - ``1`` The number of steps that MD information is output. The trajectory, velocity, and gradient will be output to disk. .. option:: temp .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``273.15`` Define the initial temperature for NVE simulation, and the targeted temperature for NVT simulation. .. option:: temp_nhc_length .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - An integer * - Default - ``5`` The length of Nose-Hoover chain to control the temperature. Usually an integer larger than 2 is enough. .. option:: temp_nhc_tau .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``0.5`` The time-constant of Nose-Hoover chain to control the temperature in ps. .. option:: vel_fn .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A file name * - Default - None If this is given, the initial velocities will be set from this XYZ file (in Angstrom/ps). Otherwise, velocities will be generated from Boltzmann distribution. For example: .. code-block:: bash :linenos: md vel_fn mol-vel.xyz end Initial velocities will be read from ``mol-vel.xyz``. .. option:: max_dr .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``5`` The maximum allowed atomic displacement per step in Angstrom. If an atom mover more than ``max_dr`` in a step, Qbics will stop the MD simulation since the system may crash due to the too fast motion. .. option:: plumed_fn .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A file name * - Default - None If this is given, Qbics will apply enhanced sampling algorithms of PLUMED defined in this file. For PLUMED input file grammar, please refer to its `PLUMED manual `_. Theoretical Background ----------------------- Molecular Dynamics Simulations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Molecular dynamics (MD) simulations are a computational method to study the time evolution of a system of interacting particles. In MD simulations, the forces between particles are calculated using a potential energy function, and the equations of motion are solved to obtain the trajectory of the particles. In Qbics, the equations of motions are solved using the velocity-Verlet algorithm, which is a symplectic integrator, meaning that it conserves the total energy of the system over long time. The forces can be evaluated on several levels of theory, including DFT, MM, NDDO, xTB, DFT/MM, NDDO/MM, and xTB/MM. During MD simulations, Qbics will output several files: - ``x-md-traj.xyz``: The trajectory of the system in XYZ format. - ``x-md-vel.xyz``: The velocity of the system in XYZ format. - ``x-md-grad.xyz``: The gradient of the system in XYZ format. Note that force is the negative of gradient. - ``x-md-log.txt``: The log file of the MD simulation, containing the step, time, energy components, temperature, and other information. - ``x-md-pressure.txt``: The pressures of the MD simulation, containing the step, time, cell volume, and pressure tensor. If the system is not periodic, this file will not be generated. If ``plumed_fn`` is given, Qbics will also output some data defined in this file. Before doing MD, it is strongly recommended to perform a **geometry optimization** to obtain a stable initial configuration. Each ensemble has its own conserved quantity. - For NVE, the total energy is conserved: .. math:: E = K + V Here :math:`E` is the total energy, :math:`K` is the kinetic energy, and :math:`V` is the potential energy. For example, if you do an NVE MD simulation, you can check the total energy in file ``x-md-log.txt``: .. code-block:: bash :linenos: # Time Total Kinetic Potential Temp Cost ps Hartree Hartree Hartree Kelvin second 0 0.00000 -34.44404963 0.02280104 -34.46685067 300.00000000 0.257 20 0.01000 -34.44402464 0.01224859 -34.45627323 161.15837759 0.260 40 0.02000 -34.44400898 0.01020550 -34.45421449 134.27682685 0.222 60 0.03000 -34.44401439 0.01628312 -34.46029751 214.24190576 0.220 80 0.04000 -34.44402295 0.02160569 -34.46562864 284.27246073 0.279 100 0.05000 -34.44403491 0.03125266 -34.47528758 411.20057981 0.248 120 0.06000 -34.44402767 0.02061378 -34.46464145 271.22164167 0.228 140 0.07000 -34.44403668 0.02599444 -34.47003112 342.01656058 0.226 160 0.08000 -34.44402275 0.01897475 -34.46299751 249.65643218 0.224 - For NVT, the total "virtual" energy is conserved: .. math:: Q = K + V + \sum_{j=1}^{n_\text{NHC}} \frac{p_{\eta_j}^2}{2Q_j} + 3NkT\eta_1+kT \sum_{j=2}^{n_\text{NHC}}\eta_j Here :math:`n_\text{NHC}` is the Nose-Hoover chain length ``nhc_temp_length``, :math:`k` is the Boltzmann constant, :math:`T` is the temperature, :math:`\eta_j` and :math:`p_{\eta_j}` are the thermostat variables ("position" and "momentum"), and :math:`Q_j` is the thermostat mass (determined automatically by Qbics). For example, if you do an NVT MD simulation, you can check the ``ConstQ`` in file ``x-md-log.txt``: .. code-block:: bash :linenos: # Time Total Kinetic Potential Temp ConstQ Cost ps Hartree Hartree Hartree Kelvin Hartree second 0 0.00000 -34.44404963 0.02280104 -34.46685067 300.00000000 -34.44404963 0.000 20 0.01000 -34.44404123 0.01734433 -34.46138557 228.20455596 -34.44404256 0.293 40 0.02000 -34.44401539 0.01535660 -34.45937199 202.05137626 -34.44402161 0.230 60 0.03000 -34.44399789 0.01914301 -34.46314090 251.87027197 -34.44401529 0.213 80 0.04000 -34.44401533 0.03373382 -34.47774915 443.84591789 -34.44404524 0.303 100 0.05000 -34.44399074 0.01808475 -34.46207549 237.94640025 -34.44402711 0.222 120 0.06000 -34.44400165 0.03088255 -34.47488420 406.33083745 -34.44404511 0.218 140 0.07000 -34.44397432 0.01350937 -34.45748370 177.74683071 -34.44402030 0.327 160 0.08000 -34.44397272 0.01886731 -34.46284002 248.24278233 -34.44402353 0.222 180 0.09000 -34.44398418 0.03438731 -34.47837149 452.44411562 -34.44404296 0.225 In all cases, if the conserved quantity is not conserved, it means that the simulation is unstable. You may need: - Reduce the time step; - Do a geometry optimization to obtain a stable initial configuration; - Check if the potential energy function (quantum chemistry or force field parameters) is reasonable for your system. However, if ``plumed_fn`` is given, the conserved quantity may **NOT** be conserved, since the enhanced sampling algorithm may change the total energy. Nose-Hoover Chain algorithm ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The Nose-Hoover chain (NHC) algorithm is a method to control the temperature in MD simulations. It is a generalization of the Nose-Hoover thermostat. The NHC is a series of thermostats that are connected to each other. The first thermostat controls the temperature of the system, and the subsequent thermostats control the temperature of the previous thermostat. The length of NHC is controlled by ``temp_nhc_length``. It should be larger than 2. If ``temp_nhc_length`` is set to ``1``, the system may not have a correct canonical distribution of velocities. Input Examples ------------------------- Example: Dynamic Effects for Carbocation Reaction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Consider the following reaction: We want to know which group is more likely to migrate. Thus, we will do a 200 fs NVT MD simulation from the optimized structure of this carbocation to see what will happen. We will repeat the MD for 20 times using a shell script ``run.sh``. Here we use xTB for its speed. You can also use DFT if you have enough time. The input files are shown below: .. tabs:: .. tab:: md-1.inp .. code-block:: bash :linenos: :caption: md-1.inp mol C -0.90849 0.15356 -0.48118 C -0.18840 -0.79035 0.24223 C -0.79661 -1.92884 0.73482 C -2.15434 -2.10431 0.51590 C -2.88201 -1.15101 -0.17642 C -2.26735 -0.01407 -0.67971 H -0.22972 -2.66722 1.27887 H -2.64692 -2.98493 0.89733 H -3.93966 -1.29437 -0.33076 H -2.84788 0.72259 -1.21038 C 1.26411 -0.39651 0.38023 C -0.02903 1.30229 -0.93562 C 1.67317 -0.22550 1.85240 H 1.83624 -1.21169 2.27913 H 2.60525 0.32605 1.95669 H 0.88426 0.23363 2.45602 C -0.71228 2.64784 -0.65267 H -0.02047 3.48874 -0.70508 H -1.46470 2.82054 -1.41711 H -1.26222 2.63862 0.29289 C 1.31204 0.95650 -0.38570 H 2.18545 1.61818 -0.54700 C 2.13967 -1.43173 -0.30477 H 1.61435 -2.09539 -0.97595 C 3.43556 -1.57292 -0.10516 H 3.99105 -2.34108 -0.61298 H 4.00084 -0.94936 0.56664 C 0.38001 1.16971 -2.43593 C 1.13470 2.39986 -2.94354 C -0.80675 0.90408 -3.35872 H 1.04783 0.30446 -2.50076 H 1.96481 2.68234 -2.30028 H 1.54266 2.18219 -3.92659 H 0.47031 3.25225 -3.05019 H -1.34333 0.00481 -3.07250 H -1.49736 1.74331 -3.36799 H -0.44263 0.76667 -4.37283 end basis def2-svp end pseudopotential def2-ecp end scf charge +1 spin2p1 1 end xtb chrg +1 end md type nvt dt 0.001 # 0.001 ps = 1 fs num_steps 200 # 0.001*200 = 0.2 ps. In real simultion, you may need > 1 ps. temp 300 temp_nhc_length 5 temp_nhc_tau 0.5 output_freq 10 # 0.001*10 = 10 fs end task md xtb # md m06 # If you want to use DFT. end .. tab:: run.sh .. code-block:: bash :linenos: :caption: run.sh for i in {1..20}; do fn=md-1-$i cp md-1.inp ${fn}.inp qbics ${fn}.inp -n 8 > ${fn}.out done By ``bash run.sh``, you can run 20 MD simulations. The output files will be ``md-1-1-md-traj.xyz``, etc. We observe that both the isopropyl and vinyl group can migrate. For example: .. tabs:: .. tab:: Trajectory: Isopropyl Group Migration .. figure:: figs/md-1-iso.gif .. tab:: Trajectory: Vinyl Group Migration .. figure:: figs/md-1-vinyl.gif For this test, we can calculate the migration ratio of the isopropyl group to the vinyl group. For 20 trajectories, the result is: .. list-table:: :stub-columns: 1 :widths: 5 5 * - Isopropyl group migration - 7 * - Vinyl group migration - 7 * - Non-reactive - 2 * - Other reactions - 4 Of course, if the simulation time is longer, the non-reactive trajectory may become reactive. Anyway, the NVT MD result suggests that the migration ratio is about 1:1. Why do we do dynamics? Due to the **nonstatistical effect** of this reaction, we cannot predict the migration ratio by a static calculation. The dynamics result is more reliable. Please refer to the paper `J. Am. Chem. Soc. 2021, 143, 1088 `_ for more details.