.. _examples-retis-2d: RETIS in a 2D potential ======================= In this example, you will explore a rare event with the Replica Exchange Transition Interface Sampling (RETIS) algorithm. In this example, we will implement a new potential for |pyretis| and we will create three new order parameters which we will use to investigate a transition between two stable states in the potential. Using the RETIS algorithm, we will compute the rate constant for the transition between the two stable states. This example is structured as follows: .. contents:: :local: Creating the potential function ------------------------------- We will consider a relatively simple 2D potential where a single particle is moving. The potential energy (:math:`V_\text{pot}`) is given by .. math:: V_\text{pot}(x, y) = (x^2 + y^2)^2 - 10 \exp(-30 (x - 0.2)^2 -3 (y - 0.4)^2) -10 \exp(-30 (x + 0.2)^2 -3 (y + 0.4)^2) where :math:`x` and :math:`y` gives the positions. The potential energy is shown in :numref:`fig2dpot_ex_retis` below. .. _fig2dpot_ex_retis: .. figure:: /_static/examples/retis2d-hysteresis/potential-2d-hyst.png :alt: The 2D potential example :align: center The potential energy as a function of the position. We have two stable states, at (x, y) = (-0.2, -0.4) and (x, y) = (0.2, 0.4) (indicated with the circles), separated by a saddle point at (x, y) = (0, 0). We begin by creating a new potential function to use with |pyretis|, see the :ref:`user guide ` for a short introduction on how custom potential functions can be created. In short, we will now have to define a new class which calculates the potential energy and the corresponding force. We begin by simply creating a new class for the new potential function. Here, we define the parameters our new potential should accept. We implement the potential slightly more general than the potential given above so that we can change the parameters more easily, if we wish to do so, later. The potential we will be implementing is given by: .. math:: V_\text{pot}(x, y) = \gamma_1 (x^2 + y^2)^2 + \gamma_2 \exp(\alpha_1 (x - x_0)^2 + \alpha_2 (y - y_0)^2) + \gamma_3 \exp(\beta_1 (x + x_0)^2 + \beta_2(y + y_0)^2) where :math:`\gamma_1`, :math:`\gamma_2`, :math:`\gamma_3`, :math:`\alpha_1`, :math:`\alpha_2`, :math:`\beta_1`, :math:`\beta_2`, :math:`x_0` and :math:`y_0` are potential parameters. Next, we define a new method as part of our new potential class which will be responsible for calculating the potential energy and a method responsible for calculating the force. Finally, we add a method that will calculate both the potential and the force at the same time. The full Python code for the new potential class can be found below. .. pyretis-collapse-block:: :heading: Show/hide the new potential class .. literalinclude:: /_static/examples/retis2d-hysteresis/potential.py :language: python This new potential function can be included in a |pyretis| simulation by referencing it in the input file as follows (assuming that it has been stored in a file named ``potential.py``): .. literalinclude:: /_static/examples/retis2d-hysteresis/potential.txt :language: rst Creating the order parameters ----------------------------- In this case, it is not so obvious how to define the order parameter. Three simple possibilities are (see :numref:`fig2dpot_order` for an illustration): 1. The x coordinate of the particle, where we include the potential energy in order to define the two stable states. 2. The y coordinate of the particle, where we include the potential energy in order to define the two stable states. 3. The projection of the distance vector from the position of the particle (x, y) to (-0.2, -0.4) onto the straight line between (-0.2, -0.4) and (0.2, 0.4). Here, we also include the potential energy in order to define the two stable states. .. _fig2dpot_order: .. figure:: /_static/examples/retis2d-hysteresis/orderparameters.png :alt: Order parameters for the 2D example. :align: center Illustration of the three different order parameters with location of interfaces (dotted lines). From left to right: using the x position as the order parameter, using the y position as the order parameter and, finally, in the rightmost figure, using a combination of x and y as the order parameter. The white line in the rightmost figure show the line which we project the order parameter onto. The blue paths are valid paths for the final path ensembles and in the bottom row, the corresponding order parameters are shown for these paths. The light dotted interface lines represent the locations where we include the potential energy in the definition of the order parameters. The two marked points (cross and triangle) in the paths in the bottom row show the second last and last points in the paths, respectively. The two marked points (cross and circle) in the paths in the top row show the first and last points in the paths, respectively. Since we include the potential energy in the order parameter definition, we will have to create customized order parameters. There is a short :ref:`introduction ` to how the order parameters are implemented in |pyretis| in the user guide and there is also a more lengthy recipe for :ref:`creating custom order parameters ` available. Here, we will first consider the x (or y) order parameter and then the more involved combination. Creating the x (or y) positional order parameter ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The order parameter we are now creating will just be the position of the particle with two extra conditions: 1. If the position is below a certain value (``inter_a``) close to the first interface and if the potential energy is above a certain threshold value (``energy_a``), then the order parameter will not be allowed to come closer to the first interface. 2. If the position is closer to the last interface than a certain value (``inter_b``) and if the potential energy is above a certain threshold value (``energy_b``), then the order parameter will not be allowed to come closer to the last interface. Thus, in addition to the position of the particle, we need to consider the 4 additional parameters ``inter_a``, ``inter_b``, ``energy_a``, ``energy_b``. The Python code for this new order parameter is given below: .. pyretis-collapse-block:: :heading: Show/hide the new order parameter .. literalinclude:: /_static/examples/retis2d-hysteresis/order-x.py :language: python This new order parameter can be included in a |pyretis| simulation by adding the following order parameter section (assuming that the order parameter is stored in a file ``order.py``): .. literalinclude:: /_static/examples/retis2d-hysteresis/order-x.txt :language: rst Creating the projection-of-positions order parameter ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For this order parameter, we will consider the distance from the particle to the stable A state and project the distance vector onto the line joining states A and B. In addition, as in the previously defined order parameters, we will include the energy in the definition of the states. The Python code for this new order parameter is given below: .. pyretis-collapse-block:: :heading: Show/hide the new order parameter .. literalinclude:: /_static/examples/retis2d-hysteresis/order-xy.py :language: python Note that this order parameter is less general than the previous ones, for instance, we make explicit use of the location of the two minimums. If you are feeling adventurous, please try to add these locations as parameters for the order parameter. This new order parameter can be included in a |pyretis| simulation by adding the following order parameter section (assuming that the order parameter is stored in a file ``order.py``): .. literalinclude:: /_static/examples/retis2d-hysteresis/order-xy.txt :language: rst Running the RETIS simulation with |pyretis| ------------------------------------------- We have now defined the potential and the order parameters we are going to use and we are ready to run the RETIS simulations. Running the RETIS simulation using x as the order parameter ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In order to run the RETIS simulation, the following steps are suggested: 1. Create a new folder, for instance ``2D-hysteresis``. 2. Place the order parameter script in this folder and name it ``order.py``. 3. Download the initial configuration :download:`initial.xyz ` and place it in the same directory. 4. Create a new sub-directory, ``retis-x``. 5. Add the following input script :download:`retis-x.rst ` to the ``retis-x`` folder. 6. Run the simulation using (in the ``retis-x`` folder): .. code-block:: pyretis pyretisrun -i retis-x.rst -p Some examples from the analysis can be found below: .. _fig2dpot_retisx: .. figure:: /_static/examples/retis2d-hysteresis/retis-x-many.png :alt: Example paths harvested. :align: center Example of paths harvested for the RETIS simulation. All paths shown are accepted paths for the last path ensemble. .. _fig2dpot_prob_x: .. figure:: /_static/examples/retis2d-hysteresis/prob-x.png :alt: Crossing probability. :align: center The crossing probability after 1 000 000 RETIS cycles. The numerical value is :math:`3.77 \times 10^{-7} \pm 6.6 \%`. The initial flux is :math:`1.97 \pm 0.7 \%`. Running the RETIS simulation using y as the order parameter ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In order to run the RETIS simulation, the following steps are suggested: 1. Create a new folder, for instance ``2D-hysteresis``. 2. Place the order parameter script in this folder and name it ``order.py``. 3. Download the initial configuration :download:`initial.xyz ` and place it in the same directory. 4. Create a new sub-directory, ``retis-y``. 5. Add the following input script :download:`retis-y.rst ` to the ``retis-y`` folder. 6. Run the simulation using (in the ``retis-y`` folder): .. code-block:: pyretis pyretisrun -i retis-y.rst -p Some examples from the analysis can be found below: .. _fig2dpot_retisy: .. figure:: /_static/examples/retis2d-hysteresis/retis-y-many.png :alt: Example paths harvested. :align: center Example of paths harvested for the RETIS simulation. All paths shown are accepted paths for the last path ensemble. .. _fig2dpot_prob_y: .. figure:: /_static/examples/retis2d-hysteresis/prob-y.png :alt: Crossing probability. :align: center The crossing probability after 1 000 000 RETIS cycles. The numerical value is :math:`4.94 \times 10^{-7} \pm 5.6 \%`. The initial flux is :math:`1.66 \pm 0.4 \%`. Running the RETIS simulation using the projection-of-position as the order parameter ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In order to run the RETIS simulation, the following steps are suggested: 1. Create a new folder, for instance ``2D-hysteresis``. 2. Create a new sub-directory, ``retis-xy`` and move into this directory. 3. Place the order parameter script in this folder and name it ``order.py``. 4. Download the initial configuration :download:`initial2.xyz ` and place it in the same directory and rename it to ``initial.xyz``. 5. Add the following input script :download:`retis-xy.rst ` to the ``retis-xy`` folder. 6. Run the simulation using (still in the ``retis-xy`` folder): .. code-block:: pyretis pyretisrun -i retis-xy.rst -p Some examples from the analysis can be found below: .. _fig2dpot_retisxy: .. figure:: /_static/examples/retis2d-hysteresis/retis-xy-many.png :alt: Example paths harvested. :align: center Example of paths harvested for the RETIS simulation. All paths shown are accepted paths for the last path ensemble. .. _fig2dpot_prob_xy: .. figure:: /_static/examples/retis2d-hysteresis/prob-xy.png :alt: Crossing probability. :align: center The crossing probability after 1 000 000 RETIS cycles. The numerical value is :math:`5.03 \times 10^{-7} \pm 9.0 \%`. The initial flux is :math:`1.54 \pm 0.6 \%`. Comparing results from the analysis ----------------------------------- By using the three different order parameters, we obtain different estimates for the rate constant as summarized in the table below. .. _table_2d_hysteresis: .. table:: Comparison of the rate constant for the 2D potential :class: table-striped table-hover +-------------------+-------------+-------------+ | Order parameter | Rate constant / 1e-7 | + +-------------+-------------+ | | Lower bound | Upper bound | +===================+=============+=============+ | x | 7.3 | 8.3 | +-------------------+-------------+-------------+ | y | 7.7 | 8.6 | +-------------------+-------------+-------------+ | xy | 7.3 | 8.3 | +-------------------+-------------+-------------+