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
Run
timestep 0.001
thermo 500
neigh_modify exclude type 2 2
dump 1 all atom 50 dump.linaArgonPlatHemiNPT
run 5000