How to control the uniform velocity of dpd particles placed in parallel plates which are moved with constant velocity

gsourav577 · September 9, 2024, 6:27am

I am trying to simulate dilute polymer solution with uniform velocity U approaching a stationary sphere. For this purpose, I am initially using dpd particles placed in parallel plates and then moving the parallel plates with certain velocity so that dpd particles achieve velocity of plate. After that I will put beads with spring in that system then sphere made of dpd particles.
But, I am stuck in the initial phase for creating uniform flow using parallel plate. Though the flow generated by moving plates is uniform but it does not matches the plate velocity, the uniform velocity is more than plate velocity. I am attaching the input script used in my simulation
##start of script

Initialization

units lj
dimension 3
boundary p p s
atom_style dpd

Ensure ghost atoms store velocities for DPD pair style

comm_modify vel yes

Create simulation box

region simbox block -30 30 -4 4 -12 12
create_box 2 simbox

Set masses for atom types

mass 1 1.0 # Mass for plate particles (Group A)
mass 2 1.0 # Mass for fluid particles (Group B)

Define regions for plates and fluid

region lower_plate block -30 30 -4 4 -12 -10
region upper_plate block -30 30 -4 4 10 12
region FLUID block -30 30 -4 4 -10 10
region FLUID2 block -30 30 -4 4 -9.9999999 9.9999999

Create atoms for plates with high density

lattice fcc 6
create_atoms 1 region lower_plate
create_atoms 1 region upper_plate

Create atoms for fluid with lower density

lattice fcc 4
create_atoms 2 region FLUID

Group definitions

group lowerplate region lower_plate
group upperplate region upper_plate
group fluid region FLUID
group walls union lowerplate upperplate
group fluid_v region FLUID2

Define DPD settings

pair_style dpd 1.0 2.0 12345 #
pair_coeff 1 1 0.0 4.5 2.0 # No interaction between wall particles
pair_coeff 1 2 18.75 4.5 2.0 # Interaction between wall and fluid particles, initially 9.68
pair_coeff 2 2 18.75 4.5 2.0 # Interaction between fluid particles

Freeze wall particles

fix freeze_walls walls setforce 0.0 0.0 0.0

Reflective boundary conditions for fluid particles at plate boundaries

fix wall_lower fluid wall/reflect zlo EDGE
fix wall_upper fluid wall/reflect zhi EDGE

Move the plates with a certain velocity along the x-axis

variable plate_velocity equal 0.35 # Set the plate velocity
velocity lowerplate set v_plate_velocity 0.0 0.0
velocity upperplate set v_plate_velocity 0.0 0.0

Neighbor list settings

neighbor 2 bin # Increase the skin distance
neigh_modify delay 0 every 1 check yes

Time integration using NVT for thermostatting

fix 2 fluid nvt temp 1.0 1.0 0.5 # Nose-Hoover thermostat applied to fluid

Initial velocities for fluid particles

velocity fluid_v create 1.0 12345

Chunks along z-axis (cuboidal binning)

compute cuboidChunks fluid chunk/atom bin/3d x lower 1 y lower 1 z lower 1 units box

Compute temperature for fluid region

compute fluid_temp fluid temp/region FLUID

Fixes for averaging density and velocity in cuboidal bins

fix aveCuboid fluid ave/chunk 100 10 1000 cuboidChunks density/number vx vy vz norm sample file cuboid_density_velocity.txt

Output settings

thermo 100
thermo_style custom step temp c_fluid_temp

Dump temperature of the fluid region only

dump 1 fluid custom 100 fluid_temp.lammpstrj id type x y z vx vy vz fx fy fz
dump 2 all custom 100 dumpvel.lammpstrj id type x y z vx vy vz fx fy fz

restart setting
restart 5000 restart.lammpstrj

Run settings

timestep 0.0025 # Reduce the time step for stability
run 100000
##end of script
I have created two regions FLUID and FLUID2, FLUID2 region is subset of FLUID region with y range being slightly small. I have used this region to create initial random velocity in the FLUID. As If I create random velocity in FLUID region, dpd particles diffuses in wall. The wall particles are going with a velocity of about 3.8 (This is ave/chunk calculated value).
What is causing this issue, I have checked with low velocity also (plate velocity 0.05) it is also showing same phenomena. Is there any alternate way to get uniform controlled velocity?
P.S. I am new to lammps and DPD simulations in lammps

Germain · September 9, 2024, 7:52am

Hi @gsourav577,

Please have a look at the formatting and posting guidelines on this topic to help people read your input file and get some general advice on how to formulate your question and investigate your issue.

If you are new to lammps and DPD, the best place to start is the relevant literature, the manual and more importantly, get feedback to your questions from an experimented advisor/colleague face to face.

After a quick look at your input, I can already see that you are using NVT with DPD pair_style which is unnecessary. The random and dissipative forces of DPD take care of thermostatting your system. This is why you provide an input temperature. Since those forces depend on relative velocities of interacting particles, using a fix like setforce 0 0 0 seems ill advised to me. The way I’ve seen such surfaces thermalised using DPD was through tethering forces like in this article. The NVE integrator should be enough.

Overall I think you need to carefully think about your setup and look for suitable parameters in the literature.

gsourav577 · February 9, 2025, 5:49am

Hi @Germain, thank you for your response. You are right DPD is itself thermostat and NVT is not needed. Actually, I am trying to reproduce Poiseuille flow from the Paper, but using NVE integrator alone does not fix the temperature of fluid. In brief, dpd particles (solvent) is sandwiched in parallel plate with frozen dpd particles and a force field is used to drive the motion. The parameters chosen are from the paper itself.
If I use fix nvt then temperature remains in control.

My input script is

# Initialization
units           lj
dimension       3
boundary        p p s
atom_style      bond

# Ensure ghost atoms store velocities for DPD pair style
comm_modify     vel yes
#comm_modify     cutoff 2.5

# Define variables for upper and lower z-values of the fluid region
variable        fluid_z_up equal 15.25          # upper z value of fluid
variable        fluid_z_down equal -15.25       # lower z value of fluid
variable        xl equal -30
variable        xh equal 30
variable        yl equal -1.5
variable        yh equal 1.5
# Define new variables for adjusted z-values for simulation box
variable        simbox_z_up equal ${fluid_z_up}+2
variable        simbox_z_down equal ${fluid_z_down}-2

# Define variables for adjusted z-values for plates
variable        lower_plate_z_up equal ${fluid_z_down}
variable        lower_plate_z_down equal ${fluid_z_down}-2
variable        upper_plate_z_down equal ${fluid_z_up}
variable        upper_plate_z_up equal ${fluid_z_up}+2

# Create simulation box with z boundaries based on fluid_z_up and fluid_z_down
region          simbox block v_xl v_xh v_yl v_yh v_simbox_z_down v_simbox_z_up
create_box      2 simbox

# Set masses for atom types
mass            1 1.0  # Mass for plate particles (Group A)
mass            2 1.0  # Mass for fluid particles (Group B)

# Define regions for plates, fluid, and spherical wall
region          lower_plate block v_xl v_xh v_yl v_yh v_lower_plate_z_down v_lower_plate_z_up
region          upper_plate block v_xl v_xh v_yl v_yh v_upper_plate_z_down v_upper_plate_z_up
region          fluid_actual block v_xl v_xh v_yl v_yh v_fluid_z_down v_fluid_z_up


# Create atoms for plates with high density
lattice         fcc 6
create_atoms    1 region lower_plate
create_atoms    1 region upper_plate

# Create atoms for fluid with lower density
lattice         fcc 4
create_atoms    2 region fluid_actual


#group by type
group           walls type 1
group           fluid type 2

# Define DPD settings
pair_style      dpd 1.0 1.0 12345  #  standard

#Set dpd coefficient
pair_coeff      1 1 0.0 4.5 1.0  # No interaction between wall particles
pair_coeff      1 2 9.68 4.5 1.0 # Interaction between wall and fluid particles,
pair_coeff      2 2 18.75 4.5 1.0 # Interaction between fluid particles

# Freeze wall particles
fix             freeze_walls walls setforce 0.0 0.0 0.0

# Reflective boundary conditions for fluid particles at plate boundaries
fix             wall_lower fluid wall/reflect zlo v_fluid_z_down
fix             wall_upper fluid wall/reflect zhi v_fluid_z_up

# Apply a constant force to fluid particles to simulate gravity
variable        grav_x equal 0.02  # Adjust this value to set the strength of the gravitational field
fix             gravity fluid addforce v_grav_x 0.0 0.0

# Neighbor list settings
neighbor        0.3 bin  # Increase the skin distance
neigh_modify    delay 5 check yes

# Initial velocities for fluid particles
velocity        fluid create 1.0 12345


# Time integration using NVT for thermostatting
fix             2 all nve  # Plain velocity verlet integration

# Compute temperature for fluid region
compute         fluid_temp fluid temp/region fluid_actual

# Output settings
thermo          1000
thermo_style    custom step temp c_fluid_temp pe ke etotal press
dump            2    all custom 1000 dumpvelforce.lammpstrj id type x y z vx vy vz fx fy fz 
restart         100000 restart.lammpstrj 

# Run settings
timestep        0.01  # Reduce the time step for stability
run             10000000

Attached is the file for input
PouissulllenewNVE.txt (3.5 KB)
I have also attached temperature vs time step for ‘fluid’ group
PouissulllenewNVE.txt.c_fluid_tempdt0.01 This is for dt=0.01

And
PouissulllenewNVE.txt.c_fluid_tempdt0.001
For dt = 0.001

Thanks for patience and reading this post
LAMMPS (17 Apr 2024) version used

Germain · February 9, 2025, 2:21pm

Some quick comments regarding your inputs compared to the procedure of the paper:

# Create atoms for plates with high density
lattice         fcc 6
create_atoms    1 region lower_plate
create_atoms    1 region upper_plate

There is nothing mentioning a different density between walls and fluid. Are you sure about this parameter?

#Set dpd coefficient
pair_coeff      1 1 0.0 4.5 1.0  # No interaction between wall particles
pair_coeff      1 2 9.68 4.5 1.0 # Interaction between wall and fluid particles,
pair_coeff      2 2 18.75 4.5 1.0 # Interaction between fluid particles

# Freeze wall particles
fix             freeze_walls walls setforce 0.0 0.0 0.0

Since DPD uses forces and velocities in order to compute dissipative and random factor, I would recommend using exactly the same parameters as the paper (a_{wall}=5). The less parameters you modify, the less error sources you introduce.
Since dpd potential is applied to all particles including walls, it is ill advised to me to use setforce 0 to prevent particles from moving. Simply integrate the fluid group.

# Reflective boundary conditions for fluid particles at plate boundaries
fix             wall_lower fluid wall/reflect zlo v_fluid_z_down
fix             wall_upper fluid wall/reflect zhi v_fluid_z_up

If you use wall in the z direction, this makes no sense and is unnecessary.

# Apply a constant force to fluid particles to simulate gravity
variable        grav_x equal 0.02  # Adjust this value to set the strength of the gravitational field
fix             gravity fluid addforce v_grav_x 0.0 0.0

In the paper, the gravitational field is only applied after reaching thermodynamic equilibrium.
Note that the paper also modify the velocities of the particles near walls. The fix_move command might help you do that but I would recommend digging into this only after solving the temperature increase issue.

I can’t help you much further but there are already things to check for in your script. It might be good to look for direct insights from colleagues closer to your research topic and/or experienced in DPD simulation to get direct feedback on your simulation results.

srtee · February 10, 2025, 6:17am

Good job having a simple script that’s relatively easy to debug (and with a .txt version!). That allowed me to quickly observe your system’s behaviour with LAMMPS-GUI and make a few notes:

The neigh_modify delay 5 setting is very aggressive – on your timestep 0.01 system every single neighbourlist build is marked dangerous. With a non-equilibrium system your particles are changing neighbours very often and the thermostatting occurs exclusively via fluid-wall interactions, so it’s no wonder the simulation isn’t stable.

On the other hand it appears to be just fine for your timestep 0.001 system (although that’s a very conservative timestep as you know). Still – you should stick to the defaults unless you know that (1) you can neighbourlist more aggressively (less often) and (2) everything else in your simulation is okay.

Your temperature compute compute temp/region will not give you an accurate temperature reading. There are at least two factors to consider:

At steady flow, the entire group of fluid particles has a non-zero center-of-mass motion relative to the box. This “counts” as extra kinetic energy for “normal” temperature computes like compute temp/region. The command compute temp/com accounts for CoM motion and removes that component of kinetic energy. It registers much less temperature rise for timestep 0.001.
At steady Poiseuille flow, there is a non-uniform velocity profile across the box. Thus, even in the CoM frame, particles nearer each wall have an average negative velocity while particles in the middle have an average positive velocity. All of the non-zero streaming velocities make further spurious contributions to the temperature (although relatively less than the CoM motion). To be ultra-accurate you can use compute temp/profile to calculate the temperature as peculiar motion relative to the average velocity in each slice of the box – but this introduces a free parameter in terms of how many slices you use (too many and the profile becomes too noisy; too few and it becomes over-homogeneous).