I am trying to simulate an evaporation process of a liquid droplet surrounded by nitrogen particles. In my cases, the liquid argon droplet is expected to absorb heat from the nitrogen for evaporation. The Langevin thermostat is used to maintain the temperature of the nitrogen particles.
After the equilibrium process, the whole system (including nitrogen and argon) is at 100 K. During the latter evaporation process, the temperature of the nitrogen is abruptly increased to 120 K. With this 20 K difference between nitrogen and argon, the argon droplet is expected to evaporate.
Here are my codes:
read_restart 0401_1bar_den_ok.restart #Read after-equilibrium files
region dele sphere 0 0 0 200 side out #Deletion zone definition
fix 301 all nve/omp
fix 304 n2 langevin 120 120 0.5 199409 tally yes # tally yes to output
fix rigid all shake 0.0001 10 0 b 1 #constrain nitrogen bonds
fix deletion ar evaporate 100 100000 dele 199409 #perform deletion
fix 4 ar recenter 0 0 0 units box #recenter to prevent drifting
timestep 0.005
compute artemp ar temp/com
compute_modify artemp dynamic/dof yes
compute n2temp n2 temp/com
fix heat_flux all ave/time 10000 1 10000 f_304 file langevin_heat_flux.out
thermo 10000
thermo_style custom step etotal temp c_artemp c_n2temp f_deletion f_heat_flux
run 1000000
It can be seen that the nitrogen temperature (c_n2temp) is around 120 K, the temperature of the argon droplet (c_artemp) is increasing as it absorbs heat from the nitrogen surrounding it. After checking the trajectory files, I am pretty sure that the evaporation process is going as I expected.
However, the energy exchange between the nitrogen and the idealised heat reservoir (f_heat_flux) remains negative. This suggests that heat is being continuously extracted from the nitrogen group. It sounds strange to me that since the liquid droplet will be absorbing heat from the nitrogen group for evaporation, the output (f_heat_flux) should be positive, indicating that heat is being added. so, from a thermodynamic equilibrium perspective, the temperature of the nitrogen can be relatively constant.
Remember that Langevin dynamics combines a dissipative and a random term. Correct me if I’m wrong but I guess you use real units given your system. So when I see:
You are setting a damping factor of 0.5\ \mathrm{fs} which is extremely aggressive as it makes the nitrogen dynamics very viscous. This might explain why there is so much energy taken out of the system in between outputs. In the general case, I wouldn’t trust results with such a setting except with a valid reason behind it.
This makes sense. Thanks for the precision, this jumped to my eyes first.
Now I see you applied rigid restriction to your molecules, have you considered using the langevin option from the rigid integrator? With the restriction of the degrees of freedom, the computation of the correct temperature for the Langevin fix can also be an issue in the thermostatting of your system. You might want to do some test for comparison.
My guess is that this recentering is causing the argon atoms to clash with nitrogen particles, since the argon is being unphysically displaced. The thermostat then tries to continually remove this heat from the system resulting in the negative heat flux you’re seeing.
Note that fix recenter will take a shift keyword together with a group name (in this case shift all would work) to shift the entire “embedding” group. Do try that and see if it helps.
a) with fix 4 ar recenter 0 0 0 shift all units box
b) Turn off recenter fix
c) Turn off Langevin fix
The starting point to perform this three conditons are:
a) see shift all effects.
b) as your suggest, recenter fix may influence the evaporation dynamics. so i just turn it off to see the effects.
c) to check without energy added into the system, how system will evolve
Results are shown below:
a)with fix 4 ar recenter 0 0 0 shift all units box
Following is my hypothesis,
it can be seen in c), in a NVE ensemble, with the evaporation goes on, the temperature of total system is dropping with a relatively constant total energy. It is kind of meet the expectations. So, i believe relavent settings may be correct without Langevin thermostat.
However, in a) b), i would like to conclude the recenter fix have little effects on the dynamics, sorry by the way. But i noticed the sum of TotEng and f_heat_flux is relatively constant during the simulation. The total energy increased can be regard as input of Langevin thermostat. So, can i regard negative value of Langevin tally yes output as a, may be, code error?
Do you have any suggestions? Any help will be fine since it is really a fun question for me.
No need to apologise – you are doing a good job translating our advice into concrete experiments and reporting back, which is helping you to understand the system better.
I would suggest a run with the following modifications:
no fix ... evaporate;
fix langevin with an increasing temperature over the run (for example 120 to 180; you can set the final temperature to be different from the initial!).
Then we will know for sure from monitoring the system’s total energy and temperature that energy is being pumped into the system by the Langevin fix. Comparing this to the Langevin fix’s tally will then let us confirm if the sign (positive or negative) is accurate.
