[lammps-users] NPT with rigid bodies in fluid

I am new at MD and I have a lot of conceptual and technical problems with my simulation.

My simulation has a upper and lower surfaces that don’t cover the length of the simulation cell in X, meaning that to the right and left there are two “reservoir regions”, like in Uzi Landmann method of GCMD (Gran canonical MD). The rest of the simulation is filled with a fluid (I artificially create a air bubble in the center, that disappears but that’s another issue)

My problems are as follows:

I am using a NPT ensemble and letting the X dimension extend and contract to maintain a P=1atm, and not controlling the other 2 dimensions.

If I do the NPT over all the system (fix 1 all npt 0.826 0.826 0.025 yz .0024 0.0024 NULL NULL NULL NULL 0.5) (I want a T=100k, P=1atm,Tdamp=50fs,Pdamp=1000fs)
then the solid is stretched to the point where the fluid atoms flow in between the slabs, which doesn’t make sense.

If I do the NPT over the fluid (fix 1 fluid npt 0.826 0.826 0.025 yz .0024 0.0024 NULL NULL NULL NULL 0.5 dilate partial)
then my simulation box stretches a lot to the point where in my screen I only see a few atoms of the 20,000 and the two solids, most of the atoms are near the vertical walls where I can’t see them. It seems to contract a little and then it stretchs a lot.

My concern is that if I do NPT only on the liquid, the liquid molecules won’t interact with the solid. But if I NPT over all, then the solid gets deformed and the solid atoms interact between them, but I think I can correct that by using neigh_modify over them. To stretch in the X direction, is LAMMPS calculating the pressure only on the walls? which in my case shouldn’t include the virial term between the frozen atoms, right? What about the pressure on the other 4 faces, where part of the atoms are frozen, how is that pressure calculated? I don’t know if my problem is that the parameters are unrealistic for a platinum solid and an argon flui, or if I’m converting them wrong.

I really appreciate all the help I can get with this. I am confused.
Lina Merchan

This is my code:

3-d LJ diamond crystal simulation

dimension 3
boundary p p p
#units lj m=40g/mol sigma=0.34nm
atom_style atomic
neighbor 0.3 bin
neigh_modify delay 5

create geometry

lattice fcc 2.65
region box block 0 35 0 15.9 -12.9 12.9
create_box 2 box

atom regions

region upper block 10.3 24.7 12.3 INF INF INF units box
region out-upper block 8.3 26.7 10.3 INF INF INF units box side out
region lower block 10.3 24.7 INF 3.5 INF INF units box
region out-lower block 8.3 26.7 INF 5.5 INF INF units box side out
region boundary union 2 upper lower
region bubble sphere 16.7 5.5 0 2.5 side out
region ifluid intersect 3 out-upper out-lower bubble
region fluid intersect 2 out-upper out-lower

#create surface with bubble

create_atoms 2 region boundary

lattice fcc 0.62

create_atoms 1 region ifluid

group fluid region fluid
group boundary region boundary
group upper region upper
group lower region lower

mass 1 1.0 #argon 39.948/40
mass 2 5.0 #platinum 195.084/40

LJ potentials

pair_style lj/cut 3.5
pair_coeff 1 1 1 1 3.5 #argon
pair_coeff 2 2 50 0.73 3.5 #platinum
pair_coeff 1 2 1.2 0.86 3.5 #platinum

#define groups

set fluid atom 1
set boundary atom 2

initial velocities

velocity fluid create .826 482748 loop geom

fix 1 fluid npt .826 .826 0.025 yz .0024 0.0024 NULL NULL NULL NULL 0.5 dilate partial
fix 2 boundary setforce 0.0 0.0 0.0


timestep 0.001
thermo 500
neigh_modify exclude type 2 2

dump 1 all atom 50 dump.linaArgonPlatHemiNPT
run 5000


The main issue is that your Argon is vaporizing. If you change this line:

fix 1 fluid npt .826 .826 0.025 yz .0024 0.0024 NULL NULL NULL NULL 0.5 dilate partial


fix 1 fluid npt .826 .826 0.025 yz .0024 0.0024 NULL NULL NULL NULL 5.0 dilate partial

it still expands a lot, but much slower and more controlled. I think that if run long enough, it might eventually contract back down to where it should be so that you can run your equilibrium MD simulation.

The issue is that your initial Ar configuration is at relatively high T and density and initially in an unstable high energy state. It needs to gradually relax and still remain in the liquid state. So you need a larger Pdamp value to slow the expansion and hopefully eventually contract.

I think that doing NPT only on the liquid is better than trying to do it on the entire system.

As far as how the virial is computed, LAMMPS includes all of the pair interaction forces that are computed, so in your case, that means all of the 1-1 and 1-2 interactions, but none of the 2-2 interactions.

Best wishes,