A Maxwell-demon-like Mechanism Verified by Molecular Dynamics Simulations


 Here we report a series of molecular dynamics simulations to confirm the feasibility of a molecule-rectifying mechanism which can continuously induce a net particle flow but has to work with an external electric field. Here we also propose an optimized model that can work without the external electric field. The results of a series of simulations show that our new model can also continuously induce net flows without any external forces. The new model can generate a considerable vapor pressure difference of up to 7.073kPa at a temperature of 370K. The new model will be easier to be verified by physical experiments and can be used to develop useful nanodevices. It is generally believed that it is impossible to exploit the kinetic energy of molecules in thermal motion at equilibrium state, but our simulation results may change this view.


Introduction
In the past two decades, many nano-scale directional-flow-inducing mechanisms have been proposed [1][2][3][4][5][6][7][8] . They share a common feature--the net particle flow is generated by breaking the dynamic symmetry of the system with external force. Recently, we also proposed a new net-flow-inducing mechanism in ref. 9 which can induce a continuous net flow without continuous energy input(for brief introduction about this mechanism, see Sec.1.1 in SI). If it really works, that will be of great theoretical and practical significance because it is considered impossible to exploit the kinetic energy of molecules in thermal motion. Due to the importance of salt aqueous solutions in biology and chemistry, their dynamics have been in-depth studied through molecular dynamics simulation(MD), from the solvation structures [10][11][12][13][14] , dynamics near vapor-liquid interface [15][16][17] to surface potential and surface tension, etc 18,19 , providing a lot of microscopic details that are not easy to obtain from physical experiments. We will verify our net-flow-inducing mechanism by MD simulation and provide microscopic details.

Verification of the Model Proposed in ref.9 2.1 Methods
From the P-shaped model shown in supplementary Fig.1, we can modify it and get model A as shown in Fig.1 for MD simulation. The size of the simulation box is lx*ly*lz (2.6*2.6*15nm), and periodic boundary conditions are applied to all the three directions. A water slab with a thickness of about 4nm is set in the middle of the box, approximately located from 4.2 to 8.2nm on the Z-axis. The water slab is composed of 950 water molecules, and its two vapor-liquid interfaces are orthogonal to the Z-axis. Right on the water slab is a 1nm thick hydrophobic porous membrane (HPM) that is located from 8.3 to 9.3 on the Z-axis. The hydrophobic membrane is composed of molecules that have week interaction with water molecules (here we use 147 nitrogen molecules with their coordinates frozen to make the hydrophobic membrane). There is a circular pore with a diameter of 1.6nm in the middle of the membrane. There is only one free ion in the water slab. In order to make the system electrically neutral, a counter ion is fixed on the boundary of the Z-axis of the box. A constant external electric field is applied in the Z-axis direction of the system, so the electrostatic force on the free ions in the water slab is constant and upward. The water model is TIP4P. The time step is 1 fs. The trajectory output is a frame per 1ps. A cutoff radius of 1.2nm is applied to short-range electrostatic and van der Waals interactions. The simulations are conducted at a temperature of 370K in the NVT (constant number of particles, volume, and temperature) ensemble. Nosé-Hoover thermostat with a coupling time of 0.2 ps is used to maintain the temperature of the system. The simulations are conducted with Opls_aa of Gromacs. We will respectively try different types of ions as the free ions in the water slab of model A and different external electric field strengths to see whether the model can induce a directional net flow. The net flow induced by model A is determined by counting the net number of water molecules across the reference plane (NNWM). Here we define the XY-plane at 2.5 on Z-axis as the reference plane(see the dotted line in Fig.1). If a water molecule moves downwards and crosses the reference plane, the value of NNWM will be increased by 1; if moving upwards to cross it, decreased by 1. In other words, if the value of NNWM at the end of a simulation is positive, there will be a downward-moving net vapor flow in the system.

External Electric Field of Appropriate Strength
Model A requires an external electric field with an appropriate strength. We found that the absolute value of the appropriate electric field strength is around 0.3 V/nm (i.e. if it is over 0.3 V/nm, the free ion in the water slab will escape from the liquid surface easily). Unfortunately, 0.3 V/nm is still too weak to concentrate the free ion in the upper half of the water slab. In order to shorten our simulation time, we have to reducing the ion's diffusion coefficient by modifying its mass(for details, see Sec.2.1 in SI). Therefore, we increase the mass of the free ion fi u to 1000 times the original mass of Clion in all following simulation cases of model A(i.e. fi u =35453).

Results and Discussion
Each of the simulation cases of model A has a specific number, such as A1, A2, A3, etc. In case A1, the free ion we adopt is Cl -, and the external electric field strength is -0.3 V/nm. After minimizing the system's energy, the simulation ran for 400 ns. The result shows that the Clion did not escape during the simulation, and the external electric field generated an obvious density gradient of Clion along the Z-axis in the water slab. The probability of the free Clion being present above 7.6 on the Z-axis is up to 90.01%. We also found that a few water molecules condensed around the fixed counter ion located at the boundary. In fact, the system has reached a general equilibrium in the first 1ns. We do not delete any part of the data in order to show the trend of the whole process. At 400 ns, the value of NNWM is 750 with an average flow rate of 1.875/ns, indicating that there is a net vapor flow as the arrows shown in Fig.1. There is a net evaporation on the lower liquid surface of the water slab, and the vapor molecules move across the lower boundary of the box and flow into the pore of the hydrophobic membrane. This is in good agreement with the prediction of ref. 9. The line chart a1

Comparison with pure water system
We conducted another simulation of a pure water system (case PW). On the base of model A, we remove the free ion, the counter ion, the hydrophobic porous membrane and the external electric field. Only the water slab is left. The box size, number of water molecules, temperature and other simulation settings are the same as described in Sec.2.1. The simulation ran for 400 ns, and the line chart pw Fig.2 shows the variation of NNWM value pw n as the function of time.  Table 1 lists the statistical results of each simulation case. We can see that the statistical mean of NNWM of case PW is closest to the horizontal axis, while the statistical means of NNWM of other simulation cases of model A are all positive and farther away from the horizontal axis, and their standard deviation values are also significantly larger than that of case PW. Where, N-Case denotes the number of each case, Type-i denotes the type of the free ion in the water slab; Str-EEF denotes the strength of the external electric field(V/nm); Time-sim denotes the effective simulation time(ns), and during Time-sim The free ion did not escape from the liquid surface; V-NNWM denotes the NNWM value at the last moment of Time-sim and F-rate denotes the average flow rate(water molecules/ns); Max and Min denote the maximum and minimum of the NNWM samples, respectively; V-mean And V-std denote the statistical mean and standard deviation of all NNWM samples, respectively. The sampling frequency is a sample per ps.

More Cases to Support Model A
We also tried ions of Na + K + and Fas the free ion in model A respectively. They all did not escape during simulation. We can see the results of Cases of A2, A3 and A4 in The absolute value of the external electric field strength in these cases are all greater than 0.3 V/nm, and the free ions all escaped from the liquid surface during simulation. In case A5, the Fion escaped at 458 ns, we only adopt the data from 0~455 ns; in case A6, the Clion escaped at 377 ns, we adopt data from 0~370 ns; in case A7, the Na + ion escaped at 105ns, we adopt data from 0~100 ns(see Fig.6~Fig.8 in SI). Although we regard these cases as failures, the results also show that there is a net vapor flow before the free ion escapes, and the flow rate seems to be positively correlated with the strength of the external electric field.

Continuity of the Net Vapor Flows
We specially extended the simulation Cases of A1 and A2 to 2000 ns (2 μs). Neither the Clion nor the Na + ion escaped from the liquid surface during the simulations. The results show that the NNWM value in each of the two cases continuously increase with time. Case A1-1 is the extended one of case A1. In Table 1, it can be seen that the value of NNWM is up to 2767 at 2000 ns, and the average flow rate is 1.3835/ns; case A2-1 is the extended one of case A2, the value of NNWM is up to 2458 at 2000 ns with an average flow rate of 1.229/ns. The difference between the average flow rates in the two cases is 12.6%. We are not sure whether this difference is related to the difference in adsorption concentrations between Cland Na + ions at the vapor-liquid interface 17 . Also, the statistical mean values and standard deviation values of the two cases are much greater than other cases'; their standard deviation values are even more than 16 times greater than that of case PW. The trend in each case has obviously converged to a stable flow rate, showing that the net vapor flow is continuous (For details, See Fig.9 and Fig.10 in SI). Also, in each case, although more than 2000 vapor molecules condensed on the liquid surface at the pore, the water slab remained almost the same as its initial state because all the condensed water entered the bulk of the water slab. Moreover, because the surface tension prevented the free ions from keeping moving upwards, the external electric field, in fact, did not keep performing work on the ions; therefore, the continuous net vapor flow was not driven by the work from the external electric field, which is also in line with the argument in ref. 9.
From the qualitative perspective, the results of all the above cases can strongly support the argument that model A can induce continuous net vapor flows. But it is not easy to verify model A by physical experiments. First, the masses of common ions are usually much smaller than that we adopted in these cases, the concentration gradient generated by the external electric field will be much smaller, it will require much longer time for experiments. In addition, the model is not expandable. When we tried putting two Clions in the water slab, they repelled each other and could not be concentrated near the liquid surface at the pore. Therefore, we need new models that can be easily verified by physical experiments. Ref.20 pointed out that asymmetrically charged materials could be made by surface modification. We believe that the external electric field can be replaced with some fixed charges, which can also concentrate the free ions. With this idea, we propose the following new model.

Verification of Our New Model 3.1 Methods
As shown in Fig.3, its box size, number of water molecules, location of the water slab and hydrophobic porous membrane, the diameter of the pore, water model, temperature and other simulation parameters are the same as model A (see Sec.2.1). In model B, it does not require the external electric field. The four ions placed on the internal wall of the pore act as fixed charges, and they are geometrically symmetrical on the same XY-plane at 8.5 on Z-axis. There are four free counter ions in the water slab, so the net charge of the system is zero. We will respectively try different types of ions as free ions in the water slab and different temperatures to see whether model B can induce net flows. The energy of the system will be minimized before each of the following simulation cases runs. Note that, unlike model A, we do not need to modify the mass of any ion in model B; we just adopt the default values of ion masses in Gromacs.

Results and Discussion
Each of the following simulation cases for model B has a specific number, such as B1, B2, and so on. In case B1, we adopt four Clions as the fixed charges, and adopt four Na + ions as the free ions in the water slab. The simulation ran for 400 ns. In the first few ps of the simulation, we could see that the Na + ions moved towards the fixed charges, pulling dozens of water molecules into the pore. The results show that, although the average velocity of the Na + ions is up to 580.19 m/s, the electric field of the fixed charge can still concentrate the free Na + near the liquid surface at the pore. The probability of the free Na + ions being present above 8.3 on Z-axis is 66.26%(i.e. the Na + ions are present in the pore most of the time). At 400 ns, the NNWM value is 399(for details, see Fig.11 in SI). Table 2 provides more statistical results for this case. The result shows that there is a net vapor flow moving downward in the system as the arrows shown in Fig.3. Where, N-Case indicates the number of each simulation case. Type-i indicate the type of the free ions in the water slab; Time-sim denotes the simulation time(ns), and V-NNWM denotes the NNWM value at the last moment of Time-sim. F-rate denotes the average flow rate(water molecules/ns). Max and Min denote the maximum and minimum of the NNWM samples, respectively; V-mean And V-std denote the statistical mean and standard deviation of all NNWM samples, respectively. The sampling frequency is a sample per ps.
We especially extended the simulation of Case B1 to 4000 ns (4 μs). In Table 2, Case B1-1 is the extended one of Case B1. We can see that the NNWM value is up to 4352 at 4000 ns with an average flow rate of 1.088/ns. The statistical mean and the standard deviation value is very large, the standard deviation is even more than 30 times larger than that of the case PW in Table 1. Fig.4 clearly shows the trend that the value of NNWM increases with time. It has converged to a stable flow rate and shows no signs of weakening over time. Also, the water slab remained almost the same as its initial state because all the condensed water molecules entered the bulk. Moreover, model B works without the external electric field, so the free ions will never escape from the liquid surface. It strongly supports that Model B can continuously induce a net flow without any external forces.

The Pressure Difference between the upper and lower Vapor Regions
If model B really can induce a net vapor flow, the saturated vapor pressure of the lower liquid surface of the water slab should always be greater than that of the upper one, resulting in a pressure difference between them. The following simulation case can help us measure it. Based on the configuration and settings of Case B1, we add a virtual wall on each of the Z-axis boundaries (i.e. periodic boundary conditions are applied to X and Y directions), which can prevent vapor molecules from crossing the Z-axis boundaries, making the lower and upper vapor regions independent from each other, as shown in Fig.5. The simulation ran for 1000 ns(1 μs). Our statistics show that the the vapor pressure of the upper vapor region and that of the lower one are respectively 83.872kPa and 90.945kPa; the lower one is 108.43% of the upper one with a pressure difference being 0 p =7.073kPa. Our analysis shows that the pressure difference between the upper and lower vvvvvapor regions is closely related to the ion concentration difference between the upper and lower surfaces of the water slab; also, the four Clions fixed on the inner wall of the pore(the fixed charges) also contribute to the ion concentration difference; i.e. the fixed ions and the free ions also take part in inducing the net flow(for details, see Sec.3.1 in SI). All above results clearly show that the asymmetric dynamic mechanism of Model B can generate a pressure difference between the upper and lower vapor regions without any external forces. The vapor pressure difference of 7.073 kPa is a considerable value. When a pipe is bridged between the two vapor regions, the vapor will flow from the region of high pressure to the region of low pressure. The vapor pressure difference is caused by the inherent property of dynamic asymmetry of the mechanism, so the vapor flow will never stop. Does the mechanism of Model B sound like a Maxwell demon? Table (2) also lists the results of other simulation cases of model B. In cases of B2~B5, we respectively tried ions of K + , Li + , Cland Fas free ions in model B. The results show that model B works with any of these ion as free ion(for details, see Sec.3.2 in SI). Halide and metal ions have large hydration enthalpies 22 . The water molecules in their solvation shell partially lose their degrees of freedom, so their chemical potential is lower than those far from the ions. When an ions is close to the vapor-liquid interface, the liquid surface layer will not be able to completely screen the electric field of the ion 18 , causing an additive potential to the local liquid surface, where water molecules need greater kinetic energy to escape, so the evaporation rate will decrease accordingly. We believe that model B is widely applicable to salt solutions. The results of the above simulation cases strongly support our view.

More cases to Support Model B
We modified relevant parameters based on case B1 to conduct the following simulations, so we can compare the results of the following cases to that of case B1. We tried water models of TIP3P and TIP5P respectively in cases of B6 and B7. The results show that model B also works with these two water models(for details, see Sec.3.3 in Si). We tried thermostats of Berendsen and V-rescale respectively In cases of B8 and B9. The results show that model B also works with these two thermostats(for details, see Sec.3.4 in SI). We also tried different temperatures. We set the system at 400k in Case B10 and 430K in case B11 respectively. The results show that higher temperatures result in larger net flows(for details, see Sec.3.5 in SI). We also tried different force fields. We adopted gromos53a6 in case B12 and charmm27 in case B13 respectively. The results show that model B can also induce net flows(for details, see Sec.3.6 in SI).
In all simulation cases of model B, we took the default values in Gromacs for masses of free ions, so we believe that the results are of qualitative and quantitative significance. Compared with model A, model B is obviously better because it gets rid of the external electric field and the free ions in the water slab will never escape. Also, model B is expandable. A number of units of model B can be put together to form a larger-scale system, which make it more practical and easier to be verified by physical experiments. If we replace the fixed charges with electrodes in model B, then we will be able to regulate the net flow; it can be a very useful nanodevice in the future. Finally, we have a question: if the net particle flow is continuous, is the time reversal symmetry of the system broken? That is to say, if we suddenly reverse the signs of the velocities of all particles, will the net flow continue and move in an opposite direction?

Conclusion
The kinetic energy of the molecules in thermal motion at thermal equilibrium state is regarded as useless energy that is impossible to exploited, but our job in this paper probably changes this view. The results of the simulations of the two models all show that the Maxwell-demon-like molecule-rectifying mechanism proposed in ref.9 is feasible. This is the first mechanism which captures the intermolecular interaction energy released during condensation of vapor molecules to generate a driving force. Its verification will bring us a deeper understanding of not only thermodynamics, statistical mechanics, but also even the entire natural sciences. Theory, simulation, and experiment are regarded as a trilogy of new discoveries. We have completed the first two parts. We believe that our simulation results are instructive for future physical experiments and we hope they also arouse interest of those physicists with experimental conditions.