Hello,
I have a system that involves 2 crystal slabs with atoms in fixed positions that sandwich a mix of oil and water. I would like to run a barostat that changes the box size so the water-oil region is held at 300K and 1atm, while the slab atoms’ positions are remapped to the new box dimensions.
I’ve tried using “fix 1 mobile npt temp 300 300 (100*dt) iso 1 1 (1000*dt)” (where group mobile is the water-oil atoms) but my box rapidly expands due to the high pressure the slabs introduced. Is it possible to remove the pressure contributions from the slabs when performing fix npt, or are there other fixes that would achieve what I’m trying to do?
the problem is that you have an anisotropic system, but are enforcing an isotropic box change.
that will not keep your immobilized atoms fully immobile, but also makes no sense for the system you describe.
so you need to allow only the box dimensions alongside the slabs to expand and then you can decide to either couple those two dimensions to expand those isotropically or leave them to be resized independently.
please note that you will still not have completely immobile slab atoms. either their positions will be scaled with the change of the box (dilate all), or they will remain fixed in space (dilate mobile) but due to the change in box size there will be gaps forming.
when I run your input, there is extremely large pressure, so you may need to correct your geometry first or check your force field parameters. from visualization it looks like the force field parameters for the slab atoms may be a problem here. please note that for immobile atoms you don’t need to compute interactions with themselves.
so you can use neigh_modify exclude group to remove them from the neighbor list (this does not remove the long-range coulomb contribution, though).
with this situation a significant expansion of the simulation cell has to be expected. you cannot have a 1atm pressure otherwise.
also, a time step of 1fs is not giving you a stable and energy conserving time integration for a system with hydrogen atoms. you will have to reduce this to at least 0.5fs but you probably need something smaller for good energy conservation. with that in mind, you also need to adjust the kspace convergence when switching to NPT. the choice of 1.0e-4 is giving an acceptable error on forces, but not on pressure. so you would need 1.0e-6 or so for that.
before running with fix npt, i suggest to do some initial equilibration with fix nvt and monitor pressures in x y z direction independently.
overall, for your kind of system, It may be better to just run with a fix volume and sufficiently large system, since it is quite complex to have consistent and correct forces while keeping atoms immobile and proper energy conservation. with a fixed volume this will be much easier (you may want to increase the system size to have less bias on fluctuations in x and y direction).
for your reference you may want to check out the modified input below and observe how the different modifications I have mentioned affect the simulation
axel.
#-----Initialization-------------------------------------------
units real
dimension 3
boundary p p p
atom_style full
timestep 0.25
neigh_modify delay 0 every 1 check yes
#-----Force Field (OPLS)---------------------------------------
bond_style harmonic
angle_style harmonic
dihedral_style opls
#improper_style opls
kspace_style pppm 1e-6
pair_style lj/cut/coul/long 14.0
read_data water-oil-slab.data
replicate 1 3 1
#-----Freeze CaCO3 atom positions-------------------------------
group caco3 type 1 2 3
group mobile subtract all caco3
neigh_modify exclude group caco3 caco3
#-----Output---------------------------------------------------
variable t_ps equal step*dt/1000
compute mytemp mobile temp
dump 1 all custom 100 dump.md2 id type q x y z mass vx vy vz
dump_modify 1 element C Ca O C C H H O
thermo 10
thermo_style custom step pxx pyy pzz press temp density
thermo_modify temp mytemp
#-----Begin Simulation-----------------------------------------
#fix 1 mobile npt temp 300 300 (100*dt) iso 1 1 (1000dt)
fix 1 mobile nvt temp 300 300 $(100dt)
run 1000
unfix 1
fix 1 mobile npt temp 300 300 (100*dt) x 1 1 (1000dt) y 1 1 $(1000dt)
run 10000