LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] problem with calculating viscousity
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] problem with calculating viscousity


From: Axel Kohlmeyer <akohlmey@...24...>
Date: Fri, 25 May 2018 18:30:25 -0400



On Fri, May 25, 2018 at 11:13 AM, Alireza Shadloo via lammps-users <lammps-users@lists.sourceforge.net> wrote:
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:

​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

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 1.9583     | 1.9583     | 1.9583     |   0.0 | 86.75
Bond    | 0.00083184 | 0.00083184 | 0.00083184 |   0.0 |  0.04
Kspace  | 0.09287    | 0.09287    | 0.09287    |   0.0 |  4.11
Neigh   | 0.18012    | 0.18012    | 0.18012    |   0.0 |  7.98
Comm    | 0.005444   | 0.005444   | 0.005444   |   0.0 |  0.24
Output  | 0.003695   | 0.003695   | 0.003695   |   0.0 |  0.16
Modify  | 0.015152   | 0.015152   | 0.015152   |   0.0 |  0.67
Other   |            | 0.00103    |            |       |  0.05

Nlocal:    4500 ave 4500 max 4500 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:    21131 ave 21131 max 21131 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:    2.60198e+06 ave 2.60198e+06 max 2.60198e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 2601983
Ave neighs/atom = 578.218
Ave special neighs/atom = 2
Neighbor list builds = 4
Dangerous builds = 1
Total wall time: 0:00:02


 

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

######################################################################
compute    peratom all stress/atom NULL virial
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
#########################################################################
    run        100000

and the output is:
    Step       Temp       TotEng         Press       v_press
       7          300    194.41663    15293.182     11560.68
    2000          300   -854.19867    1051.9721   -2680.5301
    4000          300   -856.42596   -589.09325   -4321.5954
    6000          300   -861.18595   -229.62695   -3962.1291
    8000          300   -870.58185    1378.5822   -2353.9199
   10000          300   -854.92008   -326.32136   -4058.8235
   12000          300   -850.47454    1177.1738   -2555.3283
   14000          300   -840.29667    3618.6325   -113.86961
   16000          300   -831.11535    672.74308    -3059.759
   18000          300   -831.23039    -2774.582   -6507.0842
   20000          300   -886.10002   -762.82538   -4495.3275
   22000          300   -858.96214     1507.502   -2225.0001
   24000          300   -858.68609   -3163.5358   -6896.0379
   26000          300   -872.12528    511.16238   -3221.3397
   28000          300   -867.83787    2817.4805   -915.02158
   30000          300   -864.69485   -1677.6946   -5410.1967
   32000          300   -886.64287   -2048.6087   -5781.1108
   34000          300   -886.10274    3377.1385   -355.36359
   36000          300   -868.25313   -2487.6172   -6220.1193
   38000          300   -883.09155   -3482.9154   -7215.4175
   40000          300   -834.51806   -4225.3875   -7957.8896
   42000          300   -868.85648     2481.121   -1251.3811
   44000          300   -822.86436    60.369273   -3672.1328
   46000          300    -832.7916   -309.43043   -4041.9325
   48000          300   -852.65089   -2237.4508   -5969.9529
   50000          300   -855.38901    2862.0323   -870.46982
   52000          300   -850.83745    -1744.503   -5477.0051

why the pressurs are not the same??can you give me some suggestions?
best
Alireza










On Thursday, May 24, 2018, 10:44:09 PM GMT+4:30, Axel Kohlmeyer <akohlmey@...24...> wrote:




On Thu, May 24, 2018 at 1:53 PM, Alireza Shadloo via lammps-users <lammps-users@...396...sourceforge.net> wrote:
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.

 


On Monday, May 21, 2018, 6:31:50 PM GMT+4:30, Axel Kohlmeyer <akohlmey@...24...> wrote:




On Sun, May 20, 2018 at 12:13 PM, <lammps-users@...396... sourceforge.net> wrote:
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.​


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
:
:

------------------------------ ------------------------------ ------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
______________________________ _________________
lammps-users mailing list
lammps-users@...396... sourceforge.net
https://lists.sourceforge.net/ lists/listinfo/lammps-users




--
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.

------------------------------ ------------------------------ ------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
______________________________ _________________
lammps-users mailing list
lammps-users@...396... sourceforge.net
https://lists.sourceforge.net/ lists/listinfo/lammps-users




--
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...6297....sourceforge.net
https://lists.sourceforge.net/lists/listinfo/lammps-users




--
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.