fit npt problem, volume not changing correctly

Dear All,

I was trying to calculate a surfactant/water system on a Al2O3 surface.

The simulation box was set up with packmol code, where I placed the Al2O3 at the bottom of a box and then inserted a certain number of surfactant (SDS) and H2O molecules. Since the Al2O3 will be treated rigid, I thought the X and Y lengths of the simulation box are also fixed. In order to make the packing process easier, the Z length was intentionally a larger number (162.5 angstrom instead of the actual 112.5 angstrom, which is estimated from the number of SDS/H2O and their density at such T, P condition).

Once the initial configuration file was prepared, I though the ‘fix npt’ will bring the volume to the correct value. So I am expecting that the X and Y lengths of the box be fixed (due to a rigid Al2O3 surface model) and the Z length will be shrinking from 162.5 angstrom to roughly 112.5 angstrom. Of course the final Z value might be a bit different, because the SDS micellization.

However, this was not I observed. With a few steps, the Z length was expanded to be around 172 angstrom and stayed there for the following 130000 steps (time step 1.0 fs). There is a vacuum like space with few SDS or H2O molecules, see attached to snapshots at step 130000.

I did monitor some other values: the X, Y lengths of the simulation box are not changing during the test run; the center of mass of the Al2O3 is not shifting; the temperature is quite nicely equilibrated around the targeted value; the pressure is fluctuating +/- 200.0 around the setup value 1.0atm, implying the pressure control was not very successful yet.

But what brothers me is, why this was happening? Why the ‘fix npt’ was not able to reduce the volume (the Z length) to a reasonable number? Shall I keep it running for much longer and hope it works? Or I am missing any key control? To me, this looks like both the Al2O3 surface and its mirror images are fixed, preventing the volume from shrinking. Is this what was happening?

The in file is attached. The core part is explained below:

---- core lines from the in file --------

This is to define the groups

group sds type 2 3 4 5 6 7

group na type 8
group al2o3 type 1 9
group water type 10 11
group mobile union sds na water

velocity mobile create 300.0 293288 dist gaussian mom yes rot yes

neighbor 0.3 bin

neigh_modify delay 0 every 1 check yes
neigh_modify exclude molecule al2o3

CFG, Oxygen of H2O is colored as Nitrogen purposely;

dump 3 all cfg 1000 cfg.*.cfg mass type xs ys zs vx vy vz fx fy fz
dump_modify 3 element Al C C O S O C Na O N H
dump_modify 3 sort id

compute tempFree mobile temp
compute pressFree all pressure tempFree

thermo 200

thermo_style custom step vol c_tempFree c_pressFree temp press
thermo_modify line one

fix 1 water shake 0.0001 20 0 b 8 a 9

fix 2 mobile npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
fix_modify 3 temp tempFree press pressFree

run 5000000

----- end of the in file --------

I tried both ‘iso’ and ‘aniso’ in the ‘fix npt’ part, no difference.

Also tried to add this command:

fix 1 al2o3 rigid single force * off off off torque * off off off

However, this brings forth a new problem: the simulation box was expanding nonstop along X, Y and Z directions.

I greatly appreciate any comment and thank you in advance for reading this message!

Best wishes,
Paul

—Paul Huang
Chemical Engineering
University of Oklahoma

Box_SDSonly.png

Box_all.png

lmp3.in (4.6 KB)

Dear All,

I was trying to calculate a surfactant/water system on a Al2O3 surface.

The simulation box was set up with packmol code, where I placed the Al2O3
at the bottom of a box and then inserted a certain number of surfactant
(SDS) and H2O molecules. Since the Al2O3 will be treated rigid, I thought
the X and Y lengths of the simulation box are also fixed. In order to make
the packing process easier, the Z length was intentionally a larger number
(162.5 angstrom instead of the actual 112.5 angstrom, which is estimated
from the number of SDS/H2O and their density at such T, P condition).

Once the initial configuration file was prepared, I though the 'fix npt'
will bring the volume to the correct value. So I am expecting that the X
and Y lengths of the box be fixed (due to a rigid Al2O3 surface model) and
the Z length will be shrinking from 162.5 angstrom to roughly 112.5
angstrom. Of course the final Z value might be a bit different, because the
SDS micellization.

Sorry, but this is an incorrect expectation. If there is vacuum there is
no rule that the box dimension along the vacuum should shrink under the NPT
ensemble. Fix npt changes the box dimension based on Pstart,Pstop values
so it can expand or shrink. If your vapor phase molecules are pushing
against the box along Z, then the box will of course expand.

However, this was not I observed. With a few steps, the Z length was
expanded to be around 172 angstrom and stayed there for the following
130000 steps (time step 1.0 fs). There is a vacuum like space with few SDS
or H2O molecules, see attached to snapshots at step 130000.

I did monitor some other values: the X, Y lengths of the simulation box
are not changing during the test run; the center of mass of the Al2O3 is
not shifting; the temperature is quite nicely equilibrated around the
targeted value; the pressure is fluctuating +/- 200.0 around the setup
value 1.0atm, implying the pressure control was not very successful yet.

But what brothers me is, why this was happening? Why the 'fix npt' was not
able to reduce the volume (the Z length) to a reasonable number? Shall I
keep it running for much longer and hope it works? Or I am missing any key
control? To me, this looks like both the Al2O3 surface and its mirror
images are fixed, preventing the volume from shrinking. Is this what was
happening?

I doubt it will change if you run it longer. If you want the box to shrink,
use a larger set of Pstart,Pstop values along x, or even change the
boundary condition along x to be shrink-wrapped (you won't be able to use
npt of course after changing the boundary condition).

Ray

Ray,

You are right that the box will not shrink even there is only one molecule in the ‘vacuum’ phase. This was exactly what I observed.

As suggested by Steve, the Pdamp I used might be too small. I will correct this in the next test.

I am also afraid that no matter how long I keep the job running, I will not see the shrink, the extent I was hoping for.

Try to set up a new configuration now, instead of using a larger Z length, this time I change the vdw radius setting in the Topo tool so that I can generate a correct size of box with desired number of species, and the bond/angle/dihedeal information is correct.

Thanks again!

Paul

Ray,

You are right that the box will not shrink even there is only one molecule
in the 'vacuum' phase. This was exactly what I observed.

As suggested by Steve, the Pdamp I used might be too small. I will correct
this in the next test.

please note, that the documentation usually makes suggestions for the
*production* part of a calculation. during equilibration, especially
the initial parts, all bets are off, and it can happen quite
frequently that does a "wrong" thing is "the right thing(tm)". since
equilibrium is a state function, it doesn't matter how you get there.
:wink:

e.g. a berendsen barostat and a dissipative thermostat (langevin,
temp/csld) are much better suited to push a system quickly to a
desired volume/temperature and to encourage equipartitioning of the
kinetic energy. you can always switch back to nose-hoover, once the
system is close to equilibrium.

how much effort and how many steps are needed, depends very much on
the individual system.

axel.

Hi Axel,

Do you have any suggestion on minimization setting? I regenerated the data file, using a reasonable Z length and packed all SDS/H2O/Al2O3. It is of course crowded and I need to minimize the box before the MD runs.

However, the minimization always stops after satisfying the criterion, “quadratic factors are zero”. Part of the log file is pasted below. The MD fails at the first iteration, due to the ‘nan’ error.

I guess I should further minimize the box before any MD iteration. But after trying different min_style, cg or sd, and changing ‘min_modify line’, from ‘forcezero’ to ‘backtrack’, or ‘quadratic’, none is working.

Right now, the minimize setting in the in file is:

minimize 1.0e-16 1.0e-18 10000 100000
min_style cg
min_modify line force zero

Thanks and have a great Independence Day!

Paul

---- part of minimization log file -----

Step Temp E_pair E_mol TotEng Press
0 0 1.1550966e+23 84502095 1.1550966e+23 1.7947248e+22
55 0 -1.0402947e+16 85801454 -1.0402947e+16 -1.1662403e+14
Loop time of 155.612 on 48 procs for 55 steps with 169200 atoms

Minimization stats:
Stopping criterion = quadratic factors are zero
Energy initial, next-to-last, final =
1.15509664449e+23 -1.04029465185e+16 -1.04029465185e+16
Force two-norm initial, final = 5.0781e+25 1.03526e+30
Force max component initial, final = 2.6018e+25 7.32042e+29
Final line search alpha, max atom move = 6.02198e-49 4.40835e-19
Iterations, force evaluations = 55 734

----end, part of minimization log file -----

---- MD failure —

Step Volume tempFree pressFre Temp Press
55 1765243.2 380.32666 -nan 314.43706 -nan
ERROR on proc 37: Non-numeric box dimensions - simulation unstable (…/pppm.cpp:1882)

---- end, MD failure —