Issues computing stress in a sdpd particles + fene springs system

LAMMPS LAMMPS Mailing List Mirror
1

Dear all,

I using lammps to simulate polymer melts and solutions based on sdpd particles joined by fene springs. I subject these systems to Reverse Poiseuille Flow using

Apply forcing

variable bodyfx atom mass*{gx}*((y>{Ly}/2.0)-(y<${Ly}/2.0))
fix reverse_periodic all addforce v_bodyfx 0.0 0.0

My simulations seem to be running ok since the velocity profile that I obtain for the solvent particles agree with V_y = rhogx/(2eta)(yLy-y^2). For the melts and solutions I obtain reasonable results (i.e., increased viscosity lowers vmax and a bit of flattening of the profiles due to power law-like thinning). Nevertheless when I am computing the stress of the system I find some strange results. For the solvent I find linear profiles that are proportional to tau_xy= rhogxLy/4*(1-4*y/Ly) but not quite correct. The scaling factor seem to depend on the system size.

I calculate the stress by chunking the domain:

compute chunk_vi all chunk/atom bin/1d y lower 0.01 units reduced

and the stress with:

compute stress_v all stress/atom NULL virial
variable sv_xx atom c_stress_v[1]
variable sv_yy atom c_stress_v[2]
variable sv_zz atom c_stress_v[3]
variable sv_xy atom c_stress_v[4]
fix av_stress_v all ave/chunk 1 {Nfreq} {Nfreq} chunk_vi v_sv_xx v_sv_yy v_sv_zz v_sv_xy file sv.av ave window 5

I divide the output stress by the per-atom-volume = vol/count(all).

Attached is an image with the stress output for a system of 31k particles, including:
Solvent only
Melts (5,10,25,50) particles per chain

Solution 25 particles per chain at 50%

All the stress profiles should agree with each other and with the imposed one tau_xy= rhogxLy/4*(1-4*y/Ly). Since the simulation seems to be working properly (i.e. velocity profiles are correct), I am assuming there is something wrong with the stress calculation.

#Additional settings

bond_style fene # or lj 1 1 1
bond_coeff 1 {H} {r0} 0 0
special_bonds fene
mass 1 {sdpd_mass} set group all meso/rho {sdpd_rho}
pair_style hybrid/overlay sph/rhosum 1 sdpd/taitwater/isothermal {sdpd_temp} {sdpd_eta} 28619
pair_coeff * * sph/rhosum {h} pair_coeff * * sdpd/taitwater/isothermal {sdpd_rho} {sdpd_c} {h}
fix integrate_fix_full all meso

Any advice on what I could be missing is welcomed!

Thanks,

David

ReverseFlow.jpg

2

The USER-SDPD package is a contributed package for which you may best get help from the author. I just did a quick check and the last time that author responded to anything on the lammps-users mailing list was over 5 years ago. however, the USER-SDPD package was added less than 2 years ago. so it is probably best to contact the author directly. his e-mail address is in the README file in the USER-SDPD folder.

axel.

3

Dear Axel,

Thanks for your quick response. I am quite certain this was working properly on some previous version of lammps. I looking at the behavior of the different pair_styles to try to discern what is going wrong with the compute stress/atom. Has there been any modifications to the implementation of this compute recently that might affect how stress is computed with respect to a given pair style/ potential ? I know it is quite the general question but maybe there is some specific change on how the interaction potential is taken into account in the virial implementation of the stress.

Thanks,

David

4

Dear Axel,

Thanks for your quick response. I am quite certain this was working properly on some previous version of lammps. I looking at the behavior of the different pair_styles to try to discern what is going wrong with the compute stress/atom. Has there been any modifications to the implementation of this compute recently that might affect how stress is computed with respect to a given pair style/ potential ?

no. there have been only cosmetic changes to either the implementation of compute stress/atom, the pair style in USER-SDPD and the Pair::ev_tally() function. If there had been any functional change, it would have been noticed, since this infrastructure is used by a large number of simulations and we would have noticed differences in our regression tests of those.
Please note that compute stress/atom is not really doing a computation, but rather registers some info, that triggers that the individual pair styles collect suitable per-atom data into designated, on-demand allocated storage, and then just makes data computed by the individual pair styles available. The pair styles call the Pair::ev_tally() or similar functions which then accumulate the per atom stress tensor.

but don’t take my work for it. you can download or clone the LAMMPS git repository from github, and you will have access to the entire change history of LAMMPS since 2005, when it was taken under source code management (using subversion at the time), and you can easily check and confirm changes to the source code any which way you like.

If there is a significant difference between running an older version of LAMMPS and the current version, you could help us be a) providing a suitable small and easy to run input deck for testing, and b) narrow down to the release where it broke (you can check out the different releases by checking out the patch_XXX or stable_XXX tags from the git clone/repo).

I know it is quite the general question but maybe there is some specific change on how the interaction potential is taken into account in the virial implementation of the stress.

no, and to make us take a closer look, you first have to convincingly prove that there is a change in LAMMPS causing a reproducible issue, and not - for example - some significant difference in your input, or that you have a miscompiled executable due to a buggy compiler or something else.
without this, your only option is to contact the author of the package, as mentioned before.

Axel.

5

Dear all,

I did not get any feedback yet from the makers of USER/sdpd, so I have been experimenting with SPH package pair styles. In particular, I am having similar issues with sph/taitwater/morris pair_style. I simplified my system to a box of solvent particles. After applying reverse poiseuille flow forcing, the system velocity achieves steady state quite quickly and the values agree with theoretical limits (i.e. vmax=Frho/8/eta(Ly/2)^2):

image.png

Nevertheless, the stress values are quite low compared to the theoretical imposed value. I checked if there was some issue with my chunking and averaging so dumped the whole stress and post processed the chunking and averaging:

image.png

Note that I needed to multiply the stress values by a huge factor in order to get them to resemble the imposed theoretical value Txy= FrhoLy/4*(1-4y/Ly) for y/Ly<0.5 and Txy= FrhoLy/4(4*y/Ly-3) for y/Ly>0.5.

Question: I read in the documentation that compute stress/atom gives the per atom stress values in pressure*volume so I have been dividing my averaged values by the per atom volume, which I estimate as vol/#particles. Is this ok?

I am running longer simulations in the hope that this eventually converges to the right value…

Below you can find the code I used:

dimension 3

units si

atom_style meso

variable n equal 50

variable l equal 1e-3

variable Ly equal 2*${l}

variable F equal 1e-4

variable mass equal 1

variable rho equal 1e3

variable eta equal 1e-3

variable c equal 1.25e-4

variable mu equal {eta}/{rho}

variable dx equal {l}/{n}

variable Lx equal 6*${dx}

variable Lz equal 6*${dx}

variable cutoff equal 3*${dx}

lattice sc ${dx}

variable mass equal {rho}*{dx}^3

variable dt equal 1e-5

region box block 0.0 {Lx} 0.0 {Ly} 0.0 ${Lz} units box

create_box 1 box

create_atoms 1 region box

mass 1 ${mass}

set group all meso/rho ${rho}

pair_style sph/taitwater/morris

pair_coeff * * {rho} {c} {eta} {cutoff}

compute rho_peratom all meso/rho/atom

fix integrate all meso

compute stress all stress/atom NULL

variable xy atom c_stress[4]

dump dump_id all custom 10000 *.dat id type x y z vx vy vz v_xy

dump_modify dump_id first yes

dump_modify dump_id sort id

dump_modify dump_id pad 8

thermo 1000

variable skin equal ${dx}/100

neighbor ${skin} bin

neigh_modify delay 0 check yes

timestep ${dt}

variable f atom ${mass}*F*((y>{Ly}/2.0)-(y<${Ly}/2.0))

fix reverce_periodic all addforce v_f 0 0

run 400000

Thanks!

David

6

Dear Lammps Users,

I discovered that the stress in the system that I was discussing in my previous email becomes unstable at long simulation times. Surprisingly the velocity profiles is only slightly affected by this.

I decided to take a step back and do some tests using a LJ fluid instead. Only to find that I have the same issues with calculated stress. Here is the code for that one:

2-d LJ flow simulation

dimension 3
boundary p p p

atom_style atomic
neighbor 0.3 bin
neigh_modify delay 5

create geometry

variable Ly equal 60
variable Lx equal ${Ly}/4

lattice sc 0.7
region box block 0 {Lx} 0 {Ly} 0 ${Lx}
create_box 3 box
create_atoms 1 box

mass 1 1.0
mass 2 1.0
mass 3 1.0

LJ potentials

pair_style lj/cut 1.12246
pair_coeff * * 1.0 1.0 1.12246

define groups

initial velocities

compute mobile all temp
velocity all create 1.0 482748 temp mobile
fix 1 all nve
#fix 2 flow temp/rescale 200 1.0 1.0 0.02 1.0
#fix_modify 2 temp mobile

Couette flow

#velocity lower set 0.0 0.0 0.0
#velocity upper set 3.0 0.0 0.0
#fix 3 boundary setforce 0.0 0.0 0.0
#fix 4 all enforce2d

Poiseuille flow

#velocity boundary set 0.0 0.0 0.0
#fix 3 lower setforce 0.0 0.0 0.0
#fix 4 upper setforce 0.0 NULL 0.0
#fix 5 upper aveforce 0.0 -1.0 0.0
#fix 6 flow addforce 0.5 0.0 0.0
#fix 7 all enforce2d
variable L equal ly
variable f atom 0.1*((y>{L}/2.0)-(y<{L}/2.0))
fix reverse_periodic all addforce v_f 0 0

compute pTemp all temp/profile 1 1 0 y 5

compute stress all stress/atom NULL
compute pstress all stress/atom pTemp
variable xy atom c_stress[4]
variable pxy atom c_stress[4]

Run

timestep 0.003
thermo 500
thermo_modify temp mobile

#dump 1 all atom 500 dump.flow
dump dump_id all custom 1000 *.lammpstrj id type x y z vx vy vz v_xy v_pxy v_f fx fy fz
dump_modify dump_id first yes
dump_modify dump_id sort id
dump_modify dump_id pad 8

dump 2 all image 1000 image.*.ppm type type &
zoom 1.6 adiam 1.0
dump_modify 2 pad 5

#dump 3 all movie 100 movie.mpg type type &

zoom 1.6 adiam 1.5

#dump_modify 3 pad 5

#run 100000

run 100000 every 1000 “velocity all zero linear” &
“velocity all zero angular”

Can someone explain what am I missing? Stress should converge to the theoretical value fn(L/2-y) [Fedosov J. of Chem. Phys- 2010]

Thanks in advance,

David

image.png

image.png

7

Dear all,

Can anyone confirm or explain the units of viscosity in the sdpd/taitwater/isothermal documentation site? They are given as dynamic viscosity = mass*distance/time. I think this is wrong. Take SI units for example, eta = [Pa s] = N/m^2 s = kg/m/s so mass/distance/time. Is this just a typo or am I missing something?

Best,

David

image.png

image.png

8

For your LJ simulations you need to explain more clearly
what you are asking. It’s unlikely there is a problem computing
stress or pressure incorrectly, either on a per-atom
basis or globally.

What do you mean that stress is unstable at long times?

Steve

image.png

image.png

9

For this Q, I think it can only be answered by the SDPD contributors.

Steve

image.png

image.png

10

Dear Steve,

Thanks for your email. I am still investigating this issue. Nevertheless, regarding my last question I think it is clear that there is a typo in the documentation regarding the units of dynamic viscosity.

" viscosity = dynamic viscosity of the fluid (mass*distance/time units) "

I will post as soon as I find a solution to the stress issue.

Best,

David

image.png

image.png

11

Dear Steve,

I would like to run a question through you and the other lammps users. This is directly related to my problem but it might be useful for other users as well to better understand how a new pair style should be properly created for lammps.

In my case, several force contributions are considered for the energy and momentum of particles. The total force which for the pair_lg_cut.cpp looks like:
f[i][0] += delx * fpair

f[i][1] += …
f[i][2] += …

In my case it is given by:

f[i][0] += delx * fpair + (velx + delx * delVdotDelR / rsq) * fvisc + f_random[0];
f[i][1] += …
f[i][2] += …

However when ev_tally is called the call for both cases remains:

ev_tally (i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz);

I checked some other pair_styles and all contributions to the force between particles are included in fpair but that is not the case for the sdpd-package. In my case, if I turn viscosity off (fvisc=0) when I compute the stress/atom I get the right stress values for my simulation.

Am I correct to assume that this is the problem that I am encountering with the stress calculation?

Best,

David

image.png

image.png

12

Hi - not clear on what you refer to as your pair style vs what is done in the sdpd package.
If you look at vanilla DPD in src/pair_dpd.cpp you’ll see that ev_tally() has flair
as an arg. But fpair includes all 3 terms in the DPD eq, including the viscous and
random terms. There is nothing magic about the name fpair.
For true pairwise forces (there are other ev_tally variants for many body forces)
you’re just supposed to pass it the total equal/opposite force on the 2 particles.

Steve

image.png

image.png