[lammps-users] Inaccurate Viscosity measurement with RNEMD for a polymer with charged groups

To anybody with advice:

I am attempting to measure the viscosity on a polymer [polyethylene with 5.4 mol% methyl methacrylic acid groups (randomly distributed)] with the RNEMD method developed by Muller-Plathe.

I have performed a thorough trouble-shooting study, and I am still unable to measure a viscosity within the same order of magnitude as is done experimentally.

First, I was able to reproduce the numbers measured by Muller-Plathe in a later publication where their group looked at temperature controllers, charged groups, and molecules with containing atoms with different masses. I chose to focus on their work for SPC/E water. I altered many variables (the controller used, the cutoff, using ewald summation, the timestep, the aspect ratio, switching the force field to TIP3P, number of slabs, and how the velocity profile is printed). The only three factors which had a significant effect on the value measured were the force field, the aspect ratio, and how the velocity profile is printed.

It appears you HAVE to have dimensions of LxLx3L, a cube never reaches steady-state. According to the paper, if the velocity swaps every 100 steps, you need to record the velocity at steps 101, 201, etc… Lammps has no way of doing this, except printing EVERY step. So, I had lammps print every step, and wrote a code to transfer the velocities from step 101, 201, etc…

With respect to force field, I investigated the force field I am using for my polymer (dreiding with non-bonded interactions for charged entities from CHARMM) by looking at the viscosity of ethylene glycol (viscosity of 16.1 cP). The viscosity I measured was approximately 12.6 cP (which is about the same % difference with experiment as was SPC/E water). So it appears the force field I am using for the polymer is acceptable for RNEMD (modulus, Tg, Tm, and density are all accurate for the polymer).

I built the polymer with a new aspect ratio (LxLx3L) and was able to get the correct modulus and density as before. When I ran the RNEMD algorithm, I got a much better velocity profile, BUT, the viscosity still calculates out to be ~1 P (exp: ~1300 P). I thought MAYBE at an elevated temperature (400 K), the temp controller has a major effect. So I redid water at a series of temperatures with NVE and NVT. The NVE dropped to 0 K very quickly. NVE with temp rescaling at the same frequency as the velocity swapping and NVT give about the same viscosity as expected (~75% of the experimental value). So it seems (at least for water), elevated temperatures don’t have an effect on the accuracy of the algorithm.

There MAY be shear thinning (swap frequency is every 100 fs with a timestep of 1.0 fs), but 3 orders of magnitude seems like a serious drop. I am going to rerun the polymer with NVE and temp rescaling and a lower swap frequency (every 1 ps), but I would LOVE to hear any advice from the community.

Some quick details:
Lammps Aug 18, 2010
pair_style lj/charmm/coul/long 10.0 12.0 (dreiding united atom force field, but charmm for acid groups)
kspace_style pppm 1.0E-6
neighbor 2.0 bin
fix 1 all nvt temp 400.0 400.0 100.0
fix 3 all ave/spatial 1 1 1 y lower 6.54010000 vx file vel.profile units box
fix m2 all viscosity 100 x y 20
timestep 1.0

I think that is about as much info as I can provide. Also, I did get the same value as the Muller-Plathe paper for water, so I know I am processing the numbers correctly (conversions and all).


I'll let others comment on RNEMD for polymers. Maybe
Matt Petersen has ideas. Since
you are only swapping velocities of one atom in the polymer
molecule, I would imagine it would take a long time to

According to the paper, if the velocity swaps every 100
steps, you need to record the velocity at steps 101, 201, etc... Lammps has
no way of doing this, except printing EVERY step

You can do this via the stagger() function in a variable
and dump modify every command. From the variable doc page:
The stagger(x,y) function uses the current timestep to generate a new
timestep. X,y > 0 and x > y is required. The generated timesteps
increase in a staggered fashion, as the sequence
x,x+y,2x,2x+y,3x,3x+y,etc. For any current timestep, the next timestep
in the sequence is returned. Thus if stagger(1000,100) is used in a
variable by the dump_modify every command, it will generate the
sequence of output timesteps:


I don't really have much in the way of advice (haven't used RNEMD for polymers), but I have a few comments that might be useful.

Regarding the length of the periodic cell. I don't think you necessarily need a 3x longer z dimension. The minimum size needed will depend on your system. You could need less, or more.

Regarding sampling frequency. The only problem I encountered when testing various sampling schemes was when I collected just after the swap, but before the next integration. Any other scheme worked more or less the same.

Regarding shear thinning. When I swapped very frequently I found that the velocity profiles became curved. The MP method assumes linear response. I don't know if the curved profiles were the result of non-linear behavior, or how you might measure the viscosity from them. But I think you will see a curved profile if there is shear thinning.


Some additional thoughts. Are you polymer chains well mixed
and intertwined. If not and you shear too fast, then you
could get an artificially low viscosity b/c the chains are not
interacting. This would look like shear thinning. You might
also try a conventional NEMD run with a tilting box
to see if you get similar results.