Using Fix NPT to adjust pressure with a fix boundary along one axis

LAMMPS LAMMPS General Discussion
,
1

Hello:

I’m trying to simulate a bulk liquid in the presence of a rigid surface. Meaning that my system which is made out of liquid and the rigid surface is periodic along x and y axis but is rigid along z axis. I have done NVT and created a restart file. Following that, I want to do a NPT simulation, like having a piston along z axis to adjust the pressure inside. This piston is made out of carbon atoms. I have attached my code here as well. The error I get is that “Cannot use fix nvt/npt/nph on a non-periodic dimension (src/fix_nh.cpp:392)”.

variable temperature equal 298.00 # in K

variable HCtype equal 1
variable CTtype equal 2
variable YCtype equal 3
variable YNtype equal 4
variable HStype equal 5 #Hsurface
variable OBtype equal 6 #silanol O
variable SiBtype equal 7 #O
variable OxStype equal 8 #silanol Si
variable SiStype equal 9 #Si
variable Ctype equal 10 #hydrophobic wall

LJ parameters are from Rossky’s 1994 JCP paper

Note the unit conversion from kJ/mol to kcal/mol

variable HCeps equal 0.01570 # in kcal/mol
variable HCsig equal 2.649532788 # in Angstroms
variable CTeps equal 0.10940
variable CTsig equal 3.399669508
variable YCeps equal 0.13410
variable YCsig equal 3.545776898
variable YNeps equal 0.13310
variable YNsig equal 3.011237667
variable HSeps equal 0.0
variable HSsig equal 0.0
variable OBeps equal 0.15504
variable OBsig equal 3.15400
variable SiBeps equal 0.12753
variable SiBsig equal 3.79500
variable OxSeps equal 0.15504
variable OxSsig equal 3.15400
variable SiSeps equal 0.12753
variable SiSsig equal 3.79500
variable Ceps equal 0.066
variable Csig equal 2.53355

variable timestep equal 1.00 # in fs
variable dumptime equal 10 # output quantities every 5 time steps
variable thermotime equal 10 # output thermo quantities every 5 time steps

processors * * *

units real

boundary p p p # fixed boundaries along z to be simulated using a piston

atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj 0 0 0.5 coul 0 0 0.83

pair_style hybrid lj/cut/coul/long 14.0 lj/cut 14.0
pair_modify pair lj/cut shift yes # for the repulsive wall
kspace_style pppm 1.0e-5
#kspace_modify slab 3.0

read_restart tmp.restart.45000000
neighbor 3.0 bin

pair_coeff {HCtype} {HCtype} lj/cut/coul/long {HCeps} {HCsig}
pair_coeff {CTtype} {CTtype} lj/cut/coul/long {CTeps} {CTsig}
pair_coeff {YCtype} {YCtype} lj/cut/coul/long {YCeps} {YCsig}
pair_coeff {YNtype} {YNtype} lj/cut/coul/long {YNeps} {YNsig}
pair_coeff {HStype} {HStype} lj/cut/coul/long {HSeps} {HSsig}
pair_coeff {OBtype} {OBtype} lj/cut/coul/long {OBeps} {OBsig}
pair_coeff {SiBtype} {SiBtype} lj/cut/coul/long {SiBeps} {SiBsig}
pair_coeff {OxStype} {OxStype} lj/cut/coul/long {OxSeps} {OxSsig}
pair_coeff {SiStype} {SiStype} lj/cut/coul/long {SiSeps} {SiSsig}

pair_coeff * {Ctype} lj/cut {Ceps} ${Csig} 2.8438137 # Rmin = 2**(1/6)*Csig
pair_modify mix arithmetic

group sio2 type 5 6 7 8 9
group fluid type 2 1 3 4
group piston type 10
group pfluid union fluid piston

neigh_modify delay 10 every 1 check yes
neigh_modify exclude type 6 5

timestep ${timestep}

compute temp_pres all temp
compute fluid_temp fluid temp
compute Pres all pressure temp_pres

fix rigidwall piston rigid single force * off off on torque * off off off
fix 3 pfluid npt temp {temperature} {temperature} 10 z 1 1 100 dilate piston
fix_modify 3 temp fluid_temp
fix 2 sio2 setforce 0.0 0.0 0.0

velocity fluid create ${temperature} 502232 dist gaussian mom yes rot yes

dump dmp all custom ${dumptime} dump.lammpstrj id mol type x y z vx vy vz fx fy fz
dump_modify dmp flush yes sort id format line “%d %d d .6f .6f .6f .6f .6f .6f .6f .6f .6f”

thermo ${thermotime}
thermo_style custom step temp epair emol ecoul evdwl elong c_Pres
thermo_modify flush yes

run 1000 # run simulation for 250 ps

2

Your code is not quoted correctly and thus not properly readable.

It also seems to be an example of mindless excessive use of variables, which make it specifically hard to read on top of the improper formatting.

It seems to me you may be confusing “rigid” and “immobile” here.

As far as I can see, your quoted input cannot create this error since it has fully periodic dimensions.

Your choice of “delay” seems rather large for a 1fs timestep.

These commands in combination make very little sense.

To simulate a piston like wall, I would use ‘p p f’ boundaries, fix aveforce with an equal style variable and then control the added force using fix controller.

3

I see. I’ll use fix average force with only z-axis force to be defined as a variable. Then using fix controller. Bu I still need the following line thought??

fix 3 fluid npt temp 300 300 100 z 1 1 100 dilate piston

4

no.

5

The “fix controller” command can simultaneously control the pressure by adjusting the piston position and the temperature by adjusting the velocity?

6

You still need to use time integration on the mobile atoms so that they move, but your don’t want fix npt.
You can either use fix nve plus a thermostat fix or use an integrator that combines time integration with thermostatting, e.g. fix nvt. For details see the list of fixes in the fix command documentation.

7

I still have two questions:

  1. If “set aveforce” is used to enables motion along z axis only then it’s parameters are 0 0 fz, where fz is to be calculated accordingly. Then how do I use a feedback loop to do this?

  2. “fix controller” command is providing this feedback loop?

8

Please read the documentation and study the examples for fix controller.

9

If I’m understanding correctly, I use NVT plus a fix controller to adjust the pressure? Does this pressure corresponds to the “thermodynamic pressure” of the system inside the confinement?

10

You are confusing this forum with an adviser. You have been given more than sufficient information and pointers to further information to implement what you were describing. It is not the job of the forum to confirm every little detail and do your thinking and science for you. Where in the end would then be your achievement?