problem with calculating viscousity

Alireza_Shadloo · May 20, 2018, 4:13pm

dear all
I am new in lammps and i try to calculate the viscousity of water molecular (spce) confined between two harmonic walls. in lammps manual, a Sample LAMMPS input script has been written for viscosity of liquid Ar. so i tried to modify this code for my simulation. so my lammps input script is:

variable T equal 300
variable V equal vol
variable dt equal 1
variable p equal 400
variable s equal 5
variable d equal $p*$s

variable kB equal 1.3806504e-23 # [J/K] Boltzmann

variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable atm2Pa equal 101325.0
variable convert2 equal {atm2Pa}*{atm2Pa}{fs2s}*{A2m}{A2m}*{A2m}
dimension 3
include system.in.init

read_data system.data

include system.in.settings
velocity all create $T 34387 rot yes dist gaussian

fix 55 up-wall nvt temp $T $T 10
fix 5 down-wall nvt temp $T $T 10
fix 4 spce nve

fix 7 spce temp/rescale 100 $T $T 0.02 1.0
fix 77 up-wall temp/rescale 100 $T $T 0.02 1.0
fix 777 down-wall temp/rescale 100 $T T 0.02 1.0 timestep {dt}

thermo $d
thermo_style custom step temp pe ke etotal enthalpy press vol ebond ecoul eangle
run 20000

variable pxy equal pxy

variable pxz equal pxz

variable pyz equal pyz

fix SS spce ave/correlate $s $p $d v_pxy v_pxz v_pyz type auto file S0St.dat ave running

variable scale2 equal {convert2}/({kB}$T)$V*s*{dt}

variable v11 equal trap(f_SS[3])*${scale2}

variable v22 equal trap(f_SS[4])*${scale2}

variable v33 equal trap(f_SS[5])*${scale2}

thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
run 10000

variable v equal (v_v11+v_v22+v_v33)/3.0
variable ndens equal count(spce)/vol

variable k equal (v_k11+v_k22+v_k33)/3.0

variable ndens equal count(spce)/vol

print “average viscosity: $v [Pa.s/ @ T K, {ndens2} /A^3”

my problem is pxy,pxz and pyz. because i think these values are for all atoms (both wall and water molecules) and i just want to calculate the viscousity of water molecule (no walls).
how can i calculate the pxy pxz and pyz just for water (no water+wall)??

can i use compute stress/atom command instead?? if yes how should i modify the fix ave/correlate???
can i use the flowing code??
:
:

compute mystress spce stress/atom

fix SS spce ave/correlate $s $p $d c_mysterss[4] c_mysterss[5] c_mysterss[6] type auto file S0St.dat ave running
:
:

akohlmey · May 21, 2018, 2:01pm

dear all
I am new in lammps and i try to calculate the viscousity of water
molecular (spce) confined between two harmonic walls. in lammps manual, a
Sample LAMMPS input script has been written for viscosity of liquid Ar.
so i tried to modify this code for my simulation. so my lammps input script
is:

you are trying to do two significant modifications at the same time while
being inexperienced. that is a very bad idea. try one at a time, i.e.
either introduce walls to an Ar system or do a bulk SPC/E water system (or
better both) to learn about the issues related to that.

axel.

Alireza_Shadloo · May 24, 2018, 5:53pm

Dear Axel

Thank you very much for your email. As you suggested in the first step i tried to simulate the water box with no wall. my simulation was favored. for the next step i am trying to calculate the viscosity of Argon with wall. Now i have some obstructions to do it. As i understood from lammps manual the pxy is for entire system of atoms and it can't use for mobile (Argon) group. so what is the way to calculating the pxy of mobile group in confined system??

dear all
I am new in lammps and i try to calculate the viscousity of water molecular (spce) confined between two harmonic walls. in lammps manual, a Sample LAMMPS input script has been written for viscosity of liquid Ar. so i tried to modify this code for my simulation. so my lammps input script is:

you are trying to do two significant modifications at the same time while being inexperienced. that is a very bad idea. try one at a time, i.e. either introduce walls to an Ar system or do a bulk SPC/E water system (or better both) to learn about the issues related to that.

axel.

akohlmey · May 24, 2018, 6:14pm

Dear Axel

Thank you very much for your email. As you suggested in the first step i tried to simulate the water box with no wall. my simulation was favored. for the next step i am trying to calculate the viscosity of Argon with wall. Now i have some obstructions to do it. As i understood from lammps manual the pxy is for entire system of atoms and it can't use for mobile (Argon) group. so what is the way to calculating the pxy of mobile group in confined system??

pressure in LAMMPS is computed from virial stress by computing the total
(sum of the per-atom) stress contributions and dividing by the volume.
please see the documentation of compute stress/atom for an example showing
how to do it for the entire system. now, the volume of a fully periodic
system is easy to determine, as it is an input parameter. the volume of a
subgroup of atoms, is much more difficult to determine (it is not a very
well defined entity). there are different ways to approximate the volume of
your argon subsystem, most of them have been discussed on the mailing list
in the past, so you should dig into the mailing list archive to find those
discussions and determine which method is most suitable to your case.

axel.

Alireza_Shadloo · May 25, 2018, 3:13pm

Dear Axel and lammps users
As you proposed me, i tried to calculate the pressure of system by calculating pres/atom. in the first step i try to calculate it for water and Argon box with no walls again. so as the lammps manual has reported (if the diagonal components of the per-atom stress tensor are summed for all atoms in the system and the sum is divided by dV, where d = dimension and V is the
volume of the system, the result should be -P, where P is the total pressure of the system.) the last 2 columns of thermo output in the below input script should be the same:

compute peratom all stress/atom NULL
compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
variable press equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)

thermo_style custom step temp etotal press v_press

but when i add the above input script to my codes the press and v_press in thermo-style are different. see below:

###########water box-spce#######

variable T equal 300
variable V equal vol
variable dt equal 1
variable p equal 400
variable s equal 5
dimension 3
neighbor 2 bin
neigh_modify delay 5
lattice fcc 0.8
include system.in.init
read_data system.data
include system.in.settings

set type 1 charge -0.8476

set type 2 charge 0.4238

group oxygen type 1

minimize 1.0e-4 1.0e-6 100000 400000

velocity all create $T 34387 rot yes dist gaussian

fix 4 all nve

fix 7 all temp/rescale 100 $T $T 0.02 1.0

timestep ${dt}

thermo $d

akohlmey · May 25, 2018, 10:30pm

please note, that a) SPC/E is a rigid water model, so you need to apply
either fix shake, fix rattle or fix rigid to keep those molecules rigid
and b) fix temp/rescale is a very, *very*, **very** bad choice for a
thermostat, especially when you want to validate settings for a production
calculation.

beyond that, i cannot reproduce your findings with the 11May 2018 LAMMPS
version. here is a correspondingly modified input using the data file from
examples/USER/tally.

axel.

units real
atom_style full
read_data data.spce

pair_style lj/cut/coul/long 12.0 12.0
kspace_style pppm 1.0e-4
pair_coeff 1 1 0.15535 3.166
pair_coeff * 2 0.0000 0.0000

bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
bond_coeff 1 1000.00 1.000
angle_coeff 1 100.0 109.47

special_bonds lj/coul 0.0 0.0 1.0

neighbor 2.0 bin

fix 1 all shake 0.0001 20 0 b 1 a 1
fix 2 all nvt temp 300.0 300.0 100.0

# there are bad choices from the example sent to the mailing list, but
there is no discrepancy for the pressure either
#fix 1 all nve
#fix 2 all temp/rescale 10 300 300 0.02 1.0

# make certain that shake constraints are satisfied when starting the real
simulation
run 0 post no

velocity all create 300 432567 dist uniform
timestep 2.0

compute spa all stress/atom NULL
compute ppa all reduce sum c_spa[1] c_spa[2] c_spa[3]
variable press equal -(c_ppa[1]+c_ppa[2]+c_ppa[3])/(3.0*vol)

thermo_style custom step temp press v_press
thermo 10

run 50

here is the output i get:

LAMMPS (11 May 2018)
  using 1 OpenMP thread(s) per MPI task
Reading data file ...
  orthogonal box = (0.02645 0.02645 0.02641) to (35.5328 35.5328 35.4736)
  1 by 1 by 1 MPI processor grid
  reading atoms ...
  4500 atoms
  scanning bonds ...
  2 = max bonds/atom
  scanning angles ...
  1 = max angles/atom
  reading bonds ...
  3000 bonds
  reading angles ...
  1500 angles
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj: 0 0 0
  special bond factors coul: 0 0 0
  2 = max # of 1-2 neighbors
  1 = max # of 1-3 neighbors
  1 = max # of 1-4 neighbors
  2 = max # of special neighbors
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj: 0 0 1
  special bond factors coul: 0 0 1
  2 = max # of 1-2 neighbors
  1 = max # of 1-3 neighbors
  2 = max # of special neighbors
Finding SHAKE clusters ...
  0 = # of size 2 clusters
  0 = # of size 3 clusters
  0 = # of size 4 clusters
  1500 = # of frozen angles
PPPM initialization ...
  using 12-bit tables for long-range coulomb (../kspace.cpp:321)
  G vector (1/distance) = 0.218482
  grid = 15 15 15
  stencil order = 5
  estimated absolute RMS force accuracy = 0.0319435
  estimated relative force accuracy = 9.61968e-05
  using double precision FFTs
  3d grid and FFT values/proc = 8000 3375
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 14
  ghost atom cutoff = 14
  binsize = 7, bins = 6 6 6
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d/newton
      bin: standard
Setting up Verlet run ...
  Unit style : real
  Current step : 0
  Time step : 1
Per MPI rank memory allocation (min/avg/max) = 26.54 | 26.54 | 26.54 Mbytes
Step Temp E_pair E_mol TotEng Press
       0 0 -16692.358 0 -16692.358 -1289.8319
Loop time of 9.53674e-07 on 1 procs for 0 steps with 4500 atoms

PPPM initialization ...
  using 12-bit tables for long-range coulomb (../kspace.cpp:321)
  G vector (1/distance) = 0.218482
  grid = 15 15 15
  stencil order = 5
  estimated absolute RMS force accuracy = 0.0319435
  estimated relative force accuracy = 9.61968e-05
  using double precision FFTs
  3d grid and FFT values/proc = 8000 3375
Setting up Verlet run ...
  Unit style : real
  Current step : 0
  Time step : 2
Per MPI rank memory allocation (min/avg/max) = 32.9 | 32.9 | 32.9 Mbytes
Step Temp Press v_press
       0 300 -918.21727 -918.21727
      10 242.20787 -1456.1649 -1456.1649
      20 247.70596 -1304.2404 -1304.2404
      30 251.54535 -828.92395 -828.92395
      40 258.60457 -611.08964 -611.08964
      50 266.04811 -492.00946 -492.00946
Loop time of 2.25743 on 1 procs for 50 steps with 4500 atoms

Performance: 3.827 ns/day, 6.271 hours/ns, 22.149 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads

Alireza_Shadloo · May 28, 2018, 10:25am

Dear Axel
I did your instructions step by step and my problem was solved. But i did not know why the fix temp/rescale is a very, very, very bad choice for a thermostat.

beyond that, How can i perform my simulation in constant temperature when i am using fix vne??? is there any good choice for a thermostat???

some articles have proposed using fix nve in constant temperaturer is better than fix nvt when your are calculating fluid properties(viscosity, thermal conductivity and …). do you have any suggestion???
I know my second question has no explicit answer but i prefer to know other ideas about it.

akohlmey · May 28, 2018, 1:49pm

Dear Axel
I did your instructions step by step and my problem was solved. But i did
not know why the fix temp/rescale is a very, *very*, **very** bad choice
for a thermostat.

do a little research on that. i've explained this many times on this
mailing list. if you just look at how it works, it should be self-evident.
i would not even want to call fix temp/rescale a thermostat. it is a brutal
hack. its only justification is how simple it is to implement. there are
some uses outside of production calculations, but even then i would argue,
that a dissipative thermostat algorithm has all the advantages and fewer of
the disadvantages.

beyond that, How can i perform my simulation in constant temperature when
i am using fix vne??? is there any good choice for a thermostat???

LAMMPS has a whole zoo of thermostats. most of them have one or more
journal articles explaining how they work and what they are good for.

some articles have proposed using fix nve in constant temperaturer is
better than fix nvt when your are calculating fluid properties(viscosity,
thermal conductivity and ...). do you have any suggestion???

i would read that as using fix nve and *no* thermostat, as any thermostat
does a manipulation of the system and thus has an impact on computed
properties. however, for a well done, well equilibrated, and otherwise only
very weakly coupled system, the impact from the thermostat should be
negligible compared to the statistical error of thermodynamic properties.

axel.

Alireza_Shadloo · May 29, 2018, 1:15pm

Dear Axel
I have returned to my main simulation after doing your Advice:" simulation of a bulk SPC/E water system". now it is time to simulate and calculate the viscousity of water molecular confined between two harmonic walls.
According to your comment and some articles i decided to apply the Nose-hoover thermostst on the walls instead of water molecules becuase as you said “any thermostat does a manipulation of the system and thus has an impact on computed properties”. So it is better to use thermostat on the walls and not on the confined water to keep the system in the constatnt temperature.(entire system= wall atoms with thermostat + water molecules with no thermostst).

But unfortunately the temperature of system does not fix in desired temperature or increases if i exert the force on water to run them. it seems fixe nvt do not work correctly and can not dissipated the energy conveyed to the fluid molecules. what is wrong?
I am using pppm solver so i have to use p p p as a boundary while i do not have any periodic condition in y direction !!!
I have used moltemplate to make input files.

variable T equal 300
variable V equal vol
variable dt equal 1
variable p equal 5
variable s equal 100
variable d equal $p*$s
dimension 3
boundary p p p
####### read input files #####################
include system.in.init
read_data system.data
include system.in.settings
set type 1 charge -0.8476 # o charged
set type 2 charge 0.4238 # H charged

########### harmonic wall ##########################
fix harmonic-up up-wall spring/self 1000.0
fix harmonic-down down-wall spring/self 1000.0

akohlmey · May 29, 2018, 2:24pm

what you are asking below are not really LAMMPS questions but questions about how to do research with MD. there are a lot of things not done properly, but i don’t have the time fill in for your adviser. just to give you a list of obvious issues:

you claim wanting to use the SPC/E water model, but that is a rigid model, yet i see no fix enforcing rigidity
have you checked, whether your choice of time step is suitable for your system. my guess would be that it is not
have you checked, whether your system can reproduce the desired temperature without the added force for a flow?
there is no equilibration
obviously, you need to adjust the added force to reach a dynamic equilibrium. have you tried other, smaller forces
obviously, you have not read the documentation (or literature, or textbooks) well enough to have found the slab option, which requires your system to be oriented differently, but specifically introduces a 1-d poisson solver to decouple coulomb interactions between periodic images.

… and this is just from a quick glance. all of these are issues intrinsic to MD simulations and that you should be tutored in this by your adviser or some competent collaborator. nobody here has the time to give you personal lessons via e-mails (which also is hugely inefficient and time consuming compared to personal tutoring).