Dear All,
I am trying to calculate the fluid shear viscosity by RNEMD method using
fix viscosity command. According to Muller-Plathe’s original paper,
swapping the momentum of atoms with the same mass conserve the linear
momentum and the kinetic energy of the system. Thus, if NVE without the
thermostat is used, the total energy of the system and also the
temperature is expect to be a constant with natural fluctuations. In my
simulation, I simulate 6x106 time steps with time step as 0.005\tau.
The swap options are defaults, as swap=1 and vtarget=INF. The results
show a continuous decreasing of the total energy and temperature. Does
anyone have any idea why the energy is not conserved?
Best
Frank
There is a bit of energy leak with each swap. This is my sense of what is happening. Each soon-to-be-swapped particle is moving in response to some local potential. When you change the velocity, the particle is likely going to be moving up a very steep potential gradient - the same gradient that was moving it rapidly in the opposite direction. The integrator can't keep up. You can use a small target velocity and swap more frequently/ swap multiple particles, use a smaller timestep, or use a *weakly* coupled thermostat (I've used Langevin) in the direction normal to the swap direction and the velocity gradient (the vorticity direction) - or some combination of these.
Matt
Also, assuming the energy *is* conserved, the total energy gets "partitioned" in an unintuitive way. When you start to shear the system the temperature will drop. This is because some of the kinetic energy goes into the imposed translation. Since the MP method has two mirrored velocity profile, the temperature is correct - the net translation cancels and the reported temperature is just the fluctuating bit. You either need to start at a higher temperature to get the one you want, or use a thermostat to put in some kinetic energy.
Matt
Dear Matt,
Thanks a lot for your reply. I am trying several suggestions you just gave to see whether I can minimize the total energy drift. Your second reply pointed out a very important issue, which is the nonuniform local temperature profile. The local temperature can be measured by subtracting out the streaming velocity, same as command compute temp/profile. As you said and explained in the MP paper, the momentum swap reduces the entropy which cool the local temperature in two swapping slabs. So there should be a parabolic local temperature profile, more accurate two mirror parabolic profile. And the two swapping slabs behave like two heat sinks, which disposes the viscous heat generated from the shear. Thus, the apparent temperature of the system, which not considering subtracting the streaming velocity, should be a constant.
So, I am wondering when you mentioned the temperature will drop. Do you mean the spatial averaged local temperature, computed by subtracting the streaming velocity or the spatial averaged apparent temperature?
Best
Frank
(assuming constant energy) Both temperatures - the apparent temperature and the local temperature that excludes the streaming velocity. If you start from an equilibrium simulation at some temperature, the kinetic energy that gets "reorganized" into the streaming velocity no longer contributes to the temperature. When you calculate the local temperature in some fast streaming bit, the streaming velocity is first taken out. That velocity that is taken out counted toward the temperature in the simulation before you started the swapping. The larger the imposed shear, the lower the temperature drops. That was my point.
Matt
I agree with the local temperature will definitely decrease. But apparent temperature should not drop, since the apparent temperature is just the computed by total kinetic energy. Swapping process preserves the total kinetic energy, and total energy the same . So, if the apparent temperature drops, that mean part of the kinetic energy converts into the potential energy, which is associated with the structure distortion induced by the shear. I think this is the only possible explanation if we assume the integrator does conserves energy.
Frank
OK, that's right. It's been a while since I've done this.
Matt