Volume Stops shrinking under NPT

LAMMPS LAMMPS General Discussion
1

Dear All,

I’m modeling a fluid on a substrate and trying to get it to its correct density on the substrate using the NPT ensemble. I’m applying the barostat just to the z-direction to keep the substrate intact. Unfortunately, the volume stops shrinking at some point after 2ps, and I’m left with a wrong density and enormous negative pressures. I would appreciate your input on the problem. The code, output, and image are as follows;

Input:

clear
units real
boundary p p p
atom_style full
read_data fluid.data
group fluid type 1:30 32:69

lattice fcc   ---
region substrate block 0.0 61.625 0.0 61.625 0.0 3.625 units box
create_atoms	31	region substrate
group substrate type 31

# ==================================== Interactions =========================
bond_style   harmonic
angle_style  harmonic
dihedral_style  harmonic
improper_style  harmonic

pair_style hybrid lj/cut/coul/long 12.0 10.0

kspace_style pppm 0.0001
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 1.0 angle yes dihedral yes

include "My_setting.input"
# ===========================================================================
neighbor 2.0 bin
neigh_modify every 1 delay 0 check no

reset_timestep  0
timestep 0.1 
thermo 1000 
thermo_style custom step temp etotal press vol density

velocity fluid create 300.0 4928459 rot yes dist gaussian
fix     1 all npt temp 300 300.0 1.0 z 1.0 1.0 1000.0 drag 1.0 dilate substrate
fix 	freeze substrate setforce 0.0 0.0 0.0


dump   1 all custom 10000 Trajectories/dump.pos.* id mol type q x y z vx vy vz
run   1000000
write_data  Trajectories/data.NPTequilibration*

output:

Step Temp TotEng Press Volume Density
0 276.88822 146025.02 44182.343 1139292.2 0.31560197
10000 299.88122 -28540.013 -533.47346 1098571.1 0.3273005
20000 299.26042 -28799.956 -729.75996 1052442.4 0.34164613
30000 300.40534 -29392.287 -509.83968 1009404.2 0.35621296
40000 301.06122 -29949.315 -364.91893 968431.91 0.37128357
50000 299.82676 -30543.019 -745.2195 927746.8 0.38756572
60000 300.84502 -30828.993 -719.37462 887383.9 0.40519426
70000 299.17647 -31676.042 -694.40367 848425.96 0.42379993
80000 300.05681 -32579.691 -800.94226 810019.63 0.443894
90000 300.84153 -33046.346 -811.81568 773491.34 0.46485699
100000 300.08891 -33858.211 -833.10227 738268.88 0.4870351
110000 300.98323 -34795.409 -791.73706 703303.46 0.51124853
120000 300.01058 -35274.475 -1044.2651 670032.65 0.53663484
130000 300.09914 -36184.233 -1344.0993 637688.12 0.56385378
140000 299.81756 -36971.788 -1298.8514 606641.95 0.59271018
150000 300.28012 -37782.043 -1276.5986 577068.92 0.62308477
160000 300.46901 -38205.26 -1340.7503 550242.96 0.65346199
170000 300.17108 -38355.818 -1544.6005 529375.93 0.67922027
180000 300.00045 -38559.119 -1445.429 517689.18 0.69455355
190000 299.66369 -38358.241 -1671.447 516771.75 0.6957866
200000 299.96648 -38388.54 -1626.7359 516583.59 0.69604004
210000 300.04314 -38441.914 -1023.9868 516317.46 0.69639879
220000 299.3715 -38874.461 -1210.1786 518287.48 0.69375177
230000 300.89808 -38803.557 -1083.6293 516885.71 0.69563319
240000 299.9742 -39186.729 -1421.647 517446.16 0.69487975

image

2

I think the whole simulation setup and model design is fundamentally flawed.
If you have a liquid on a substrate, you cannot have periodic boundaries and thus cannot use fix npt.

Also, looking at the density of the entire system is not likely to give you are meaningful number since it is not going to be homogeneous.

Furthermore, with a free surface, the density of the liquid will equilibrate to whatever density the force field determines. You probably also need a (harmonic or reflective) wall to keep atoms from leaving the box.

3

Based on what you said, I should equilibrate the fluid first, then add a substrate to the system. That step requires unfolding and refolding the fluid back to a bigger box but, I’m guessing that would provide a better result with the fluid having the correct density at the starting point of the simulation. Does that sound reasonable?

Thank you very much for your time.

Best regards,
Amir

4

No. That sounds overly complicated.

Why can’t you just create the liquid at an estimated density right away and then minimize it before assigning an initial velocity and then equilibrate it? directly with the substrate and a harmonic wall to keep atoms from escaping. If set up properly, it will be just as if you prepare an initial bulk system, but you don’t have to worry about breaking bulk system apart and re-equilibrating it when placing next to the substrate. These processes tend to create long-lasting shockwaves that may distort your results.

Or if you have to create the system at a lower density, you can just use a correspondingly larger box with the substrate on one side and the harmonic wall on the other which has its location defined by an equal style variable. That variable would be be using the current timestep to compute the position and could thus be moved toward the substrate and you can move it until you reach the desired density.

Then for the final part of the equilibration, you move the wall back and carefully let the system relax to whatever density is consistent with your force field. That won’t be a constant because of the interaction with the substrate on one side and the free surface on the other.

1 Like
5

Thank you so much. Moving a harmonic wall towards the fluid is an amazing solution. I will do as you suggested. Thank you for your time and help.

Best regards,
Amir