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: Alireza Shadloo <alirezashadloo@...16...>
Date: Fri, 25 May 2018 15:13:30 +0000 (UTC)

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

######################################################################
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@lists.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@...12...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.