LAMMPS Pressure Units Relation to Size for Polydisperse system

[almehnol]

Hello all,

I am running molecular dynamics simulations in lammps and trying to calculate how polydispersity affects the pressure. To start, I ran MD simulations with the repulsive part of the Morse potential with a small damping factor, with the goal of replicating Brownian Dynamics.

If I run the simulations (in LJ units) with a size of 1 and a volume fraction of 30%, I get close to the Carnahan Starling value as expected. However, when I say double the size, the pressure gets reduced by 1/8. This makes sense given the LJ units are P*sigma^3/epsilon. This observation holds across multiple pressure calculation commands:

compute myPress all pressure NULL virial
thermo_style custom step temp c_myPress[1] c_myPress[2] c_myPress[3] c_myPress[4] c_myPress[5] c_myPress[6] pe ke
thermo_style custom step temp c_myPress[1] c_myPress[2] c_myPress[3] c_myPress[4] c_myPress[5] c_myPress[6] press pe ke c_msd[4]

My questions are:

1. Is my interpretation that if we want to get a dimensionless pressure that depends only on volume fraction, you would scale with the particle size cubed?
2. How would this translate to a bidisperse case? It is not immediately obvious to me how to extend this because I have calculated the stress on each atom and then taken the average for each type and output that, but the relation no longer holds (in fact sometimes the trend reverses and the larger particles have a higher stress output). I’ve tried also normalizing by the Voronoi volume but couldn’t make sense of that output either.

I realize I might be overlooking something very simple, and thank you for your help and patience. Below is my code in case its helpful:

create_atoms 1 random n_part_1 seed1 simulation_box
create_atoms 2 random n_part_2 seed2 simulation_box

group group1 type 1
group group2 type 2

variable size1 equal .51
variable size2 equal 21

variable mass1 equal size1^3
variable mass2 equal size2^3
mass 1 mass1
mass 2 mass2

variable pot_depth equal 10.0
variable pot_steep equal 30.0

variable cutoff_11 equal size1
variable cutoff_12 equal (size1+size2)/2.0
variable cutoff_22 equal size2

variable pot_steep_11 equal pot_steep/size1
variable pot_steep_12 equal 2*pot_steep/(size1+size2)
variable pot_steep_22 equal pot_steep/size2

pair_coeff 1 1 pot_depth pot_steep_11 cutoff_11 cutoff_11
pair_coeff 1 2 pot_depth pot_steep_12 cutoff_12 cutoff_12
pair_coeff 2 2 pot_depth pot_steep_22 cutoff_22 cutoff_22

variable dt equal 0.0005
timestep dt

variable damp equal 10dt
variable damp1 equal dampsize1^2
variable damp2 equal damp*size2^2

minimize 1.0e-10 1.0e-12 100000 1000000

compute myPress all pressure NULL virial
compute perAtomStress all stress/atom NULL virial
compute myStress1 group1 reduce sum c_perAtomStress[1] c_perAtomStress[2] c_perAtomStress[3]
compute myStress2 group2 reduce sum c_perAtomStress[1] c_perAtomStress[2] c_perAtomStress[3]

variable press1 equal c_myStress1[1]/n_part_1
variable press2 equal c_myStress2[1]/n_part_2

compute msd1 group1 msd com yes
compute msd2 group2 msd com yes

thermo 100
thermo_style custom step temp press v_press1 v_press2 c_myPress[1] c_myPress[2] c_myPress[3] c_msd1[4] c_msd2[4]

fix nve all nve
fix mylangevin1 group1 langevin 1.0 1.0 damp1 seed
fix mylangevin2 group2 langevin 1.0 1.0 damp2 seed2

variable prod_steps equal 50000
run prod_steps

[Germain]

Hi @almehnol,

There is a lot going on here but very few related to LAMMPS. Here are some comments that might help you a bit with the physics.

@almehnol

I am running molecular dynamics simulations in lammps and trying to calculate how polydispersity affects the pressure.

Polydispersity of what? Hard spheres? Or a single system with Morse potential and several parameters? I suppose there is already some literature on both. This is the place you are more likely to find answers to your questions.

@almehnol

However, when I say double the size, the pressure gets reduced by 1/8. This makes sense given the LJ units are P*sigma^3/epsilon.

No, this might make sense because the pressure is computed as \frac{s_{ii}}{V} with s the pressure tensor computed by LAMMPS. But since you do not tell how you “double the size” this would only make sense if you did not double the quantity of matter. Pressure is an intensive quantity. Systems in the same thermodynamic states with different sizes have the same pressure value.

@almehnol

This observation holds across multiple pressure calculation commands:

compute myPress all pressure NULL virial
thermo_style custom step temp c_myPress[1] c_myPress[2] c_myPress[3] c_myPress[4] c_myPress[5] c_myPress[6] pe ke
thermo_style custom step temp c_myPress[1] c_myPress[2] c_myPress[3] c_myPress[4] c_myPress[5] c_myPress[6] press pe ke c_msd[4]

There is no differences in the pressure computation between these commands. You just add more values at the end of the printing. The press value should give you the same value as the compute MyPress all pressure NULL virial but with the kinetic term in addition.

@almehnol

My questions are:

1. Is my interpretation that if we want to get a dimensionless pressure that depends only on volume fraction, you would scale with the particle size cubed?

This makes no sense. Pressure has a dimension. Using reduced units doesn’t change that, you just express it as a ratio compare to another value.

@almehnol

How would this translate to a bidisperse case? It is not immediately obvious to me how to extend this because I have calculated the stress on each atom and then taken the average for each type and output that, but the relation no longer holds (in fact sometimes the trend reverses and the larger particles have a higher stress output). I’ve tried also normalizing by the Voronoi volume but couldn’t make sense of that output either.

This is a scientific question that should be discussed with people in the concerned field of research. There is a lot of discussion on local stress in the scientific literature and I suppose there is also in the field of polydisperse systems (I only about them from their use in forming glassy systems). But overall this is not related to LAMMPS, which is the principal aim of this forum. You will get better insights from colleagues, advisor or expert from the field whom shall orient you toward better resources on the topic than ;ost people on this forum.

[almehnol]

Hi @Germain,

Thank you for taking the time to respond. Sorry if some of my post was unclear or irrelevant. I’ll respond to your points to clarify some things.

@Germain

Polydispersity of what? Hard spheres? Or a single system with Morse potential and several parameters? I suppose there is already some literature on both. This is the place you are more likely to find answers to your questions.

Size polydispersity of hard (which I try to achieve with a strongly repulsive Morse potential that I cutoff at the particle size) spheres. I’m aware of an abundance of literature on the topic. I am having a hard time reproducing some of it in LAMMPS (going beyond the monodisperse case where I set the size to 1), which is why I posed the question.

@Germain

No, this might make sense because the pressure is computed as \frac{s_{ii}}{V}siiV\frac{s_{ii}}{V} with sss the pressure tensor computed by LAMMPS. But since you do not tell how you “double the size” this would only make sense if you did not double the quantity of matter. Pressure is an intensive quantity. Systems in the same thermodynamic states with different sizes have the same pressure value.

If I double both the particle diameter and box size (meaning same volume fraction and number of particles), I get a 1/8 times smaller pressure output. I adjust the mass by diameter cubed and the damping factor by diameter squared. For what it’s worth, I get the same behavior doing both fix brownian and fix langevin. It is of course likely I’m overlooking something obvious that is leading to this output. I agree that the pressure should be defined by a thermodynamic state variable (for hard spheres just volume fraction). The pressure I want is P/nKT.

@Germain

There is no differences in the pressure computation between these commands. You just add more values at the end of the printing. The press value should give you the same value as the compute MyPress all pressure NULL virial but with the kinetic term in addition.

I figured this would be the case and it indeed is what I saw.

@Germain

This makes no sense. Pressure has a dimension. Using reduced units doesn’t change that, you just express it as a ratio compare to another value.

I apologize for the confusion. I am trying to calculate P/nKT, which is dimensionless. I should have perhaps posed the question much differently. Maybe I should have just asked how to convert LAMMPS pressure output to P/nkT depending on the size of the particle(s) in LAMMPS LJ units.

@Germain

This is a scientific question that should be discussed with people in the concerned field of research. There is a lot of discussion on local stress in the scientific literature and I suppose there is also in the field of polydisperse systems (I only about them from their use in forming glassy systems). But overall this is not related to LAMMPS, which is the principal aim of this forum. You will get better insights from colleagues, advisor or expert from the field whom shall orient you toward better resources on the topic than ;ost people on this forum.

I agree perhaps this part of my question extends beyond the scope of this forum.

[Germain]

@almehnol

The pressure I want is P/nKT. […] I am trying to calculate P/nKT, which is dimensionless. I should have perhaps posed the question much differently. Maybe I should have just asked how to convert LAMMPS pressure output to P/nkT depending on the size of the particle(s) in LAMMPS LJ units.

Sorry but \frac{P}{nK_BT} obviously has a dimension which is [L^3]. The dimensionless compressibility factor is \frac{PV}{nK_BT}. Let me stress this: the value you say you aim at computing cannot be a pressure value.

The rest is something you should sort out depending on the technical details of your method, the model you use and the equation you expect to follow.

[almehnol]

@Germain

Sorry but \frac{P}{nK_BT}PnKBT\frac{P}{nK_BT} obviously has a dimension which is [L^3][L3][L^3]. The dimensionless compressibility factor is \frac{PV}{nK_BT}PVnKBT\frac{PV}{nK_BT}.

I wrote the symbol n to denote number density, N/V. So dimensionless. I’ve seen this called the suspension pressure in colloids papers, but it is more commonly called compressibility in the MD community. That being said, I think a lot of my confusion may stem from me not having a grasp on subtle definitions and differences in convention when it comes to normalization and non-dimensionalization in the MD and colloids literature. I will take this into consideration, and this is almost certainly beyond the scope of the forum. Thanks for your help.

[Germain]

I’ve never seen such a use of the term but this would explain that.

We had these very dense discussions about LJ units some time ago if that can help, here and there.