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
Nose-Hoover
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).
Kevin