I am trying to simulate a box composed of polymer, water and graphene sheets. Polymer and water are sandwiched between two graphene sheets in the middle and also, there are water and polymer again at the two ends of the box (the sequence is as: water/polymer+graphene + water/polymer + graphene + water/polymer). I have assembled the initial configuration in materials studio first with a lower density (so that polymers have freedom to move around and get their optimum configuration). Now, I am using lammps to run my simulation. Because I need a non-periodic boundary condition in the z direction, I used p p f boundary condition and used fix wall/reflect to prevent movement of particles in the z direction. So, now my issue is that with the current setting of my simulation condition, the box size does not change, which is necessary for me to get the proper density for my structure over time. So, I was wondering if anybody could provide me with any suggestion on this matter. I have attached my input file here, for your reference.
Your help is highly appreciated,
echo both
units real
dimension 3
boundary p p f
atom_style full
#restart 400 restart
# Pair_Coefficient ..............................................................................................................................................
pair_style lj/class2/coul/long 10 15
dielectric 1.0
bond_style class2
angle_style class2
dihedral_style class2
improper_style class2
read_data System.data
#read_restart System.data
#info coeffs out log
#...................................................................................................................................................................................
kspace_style pppm 0.0001
kspace_modify slab 3.0
neighbor 2 bin
neigh_modify every 1 delay 0 check yes
#velocity all create 300.0 497 dist gaussian
#minimize 1.0e-6 1.0e-6 5000 5000
#Exports ......................................................................................................................................
timestep 1
thermo 1000
thermo_style custom step cpu temp pe ke epair ecoul etotal density lx ly lz
dump 1 all xyz 3000 System.xyz
#reset_timestep
reset_timestep 0
#nvt run
#fix 1 all nvt temp 300.0 300.0 50.0
#run 1100000 upto
#write_restart restart-1.equil
#unfix 1
#rise temperature
#fix 2 all nvt temp 300.0 600.0 50.0
#run 1000000 upto
#undump 1
#write_restart restart-1.equil
#unfix 2
#..............................................................................................................................................
fix xwalls all wall/reflect zlo EDGE zhi EDGE
#run at high temp for 3000ps
fix 3 all npt temp 600 600 50 iso 1 1 1000
run 100000 upto
undump 1
unfix 3
#..............................................................................................................................................
I have never done something quite like you are trying but I had an idea from reading the description you gave:
Why not use a fix deform (fix deform command — LAMMPS documentation) together with the fix nvt (or any other set of fixes that culminate in equations of motion developed for sampling the NVT ensemble (or simply for carrying out the dynamics at constant volume)) for the atom dynamics? From reading the wall/reflect command, it seems that it would accept that the dimension of the simulation box changes with time (the keyword EDGE you summoened would iteratively be adapted to the edge of the simulation domain in the current timestep). Else, I suppose it is also possible to define a dynamic region using the region command (region command — LAMMPS documentation) + use the wall/region command (fix wall/region command — LAMMPS documentation).
Then you could keep your boundary command as you used it (p p f) and output Pxx, Pyy and Pzz in the log.lammps to evaluate whether or not each of them match your desired, target pressure. If you use the wall/region though it might be a problem as the forces coming from the interaction of the wall with the atoms are included in the computation of pressure (if I remember correctly).
Alternatively to the fix nvt + fix deform, I think it should also be possible to use NPT for the two dimensions only and fix deform in the z direction (the way you are writing in the current script shouldnt work I think because the box is not periodic in the z direction and the iso keyword should take into consideration Pzz as well as deformation in the z direction). Then you would need to keep tabs in the Pzz only.
Dear @ceciliaalvares , thank you very much for your comprehensive response. There are a few points that I’d like to discuss here:
I tried your approach for fix nvt + fix deform. However, the problem with fix deform is that I am not sure what should be the final size of my box. In fact, the box size should naturally shrink to get to a size which everything is wrapped together, properly. Also, when using NPT in x and y dimension it introduces wrinkles on the graphene sheet, which I do not want.
Also, for the fix wall/reflect command, it is mentioned in LAMMPS documentary that this fix is not invoked during energy minimization. So, I was wondering what is the proper sequence for mentioning this fix and energy minimization.
Lastly, when I run energy minimization at the initial stage on my input file with p p f boundary codition, I get an error like what I have attached here, which I presume could be due to the loss of atoms out of my simulation box from the fixed face. So, I am not sure how should I prevent this in my simulation.
I appreciate your consideration,
LAMMPS (20 Nov 2019)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
using 1 OpenMP thread(s) per MPI task
units real
dimension 3
boundary p p f
atom_style full
restart 400 restart
# Pair_Coefficient ..............................................................................................................................................
pair_style lj/class2/coul/long 10 15
dielectric 1.0
bond_style class2
angle_style class2
dihedral_style class2
improper_style class2
read_data System.data
Reading data file ...
orthogonal box = (0.551562 1.95093 0.044454) to (98.9516 87.1678 300.044)
4 by 4 by 10 MPI processor grid
reading atoms ...
54164 atoms
scanning bonds ...
4 = max bonds/atom
scanning angles ...
6 = max angles/atom
scanning dihedrals ...
14 = max dihedrals/atom
reading bonds ...
55442 bonds
reading angles ...
103990 angles
reading dihedrals ...
154971 dihedrals
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
4 = max # of 1-2 neighbors
7 = max # of 1-3 neighbors
18 = max # of 1-4 neighbors
20 = max # of special neighbors
special bonds CPU = 2.14571 secs
read_data CPU = 3.00426 secs
#read_restart System.data
#info coeffs out log
#...................................................................................................................................................................................
kspace_style pppm 0.0001
kspace_modify slab 3.0
neighbor 2 bin
neigh_modify every 1 delay 0 check yes
velocity all create 300.0 497 dist gaussian
minimize 1.0e-6 1.0e-6 5000 5000
PPPM initialization ...
using 12-bit tables for long-range coulomb (../kspace.cpp:323)
G vector (1/distance) = 0.152693
grid = 24 24 144
stencil order = 5
estimated absolute RMS force accuracy = 0.0321897
estimated relative force accuracy = 9.69385e-05
using single precision MKL FFT
3d grid and FFT values/proc = 12463 3636
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 17
ghost atom cutoff = 17
binsize = 8.5, bins = 12 11 36
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/class2/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Setting up cg style minimization ...
Unit style : real
Current step : 0
ERROR on proc 150: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 9: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 119: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 49: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 19: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 39: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 29: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 0: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 30: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 10: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 20: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 59: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 79: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 109: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 40: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 80: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 50: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 100: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 110: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 60: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 70: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 120: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 130: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 140: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 129: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 149: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 159: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
ERROR on proc 139: Bond atom missing in image check (../domain.cpp:765)
Last command: minimize 1.0e-6 1.0e-6 5000 5000
Fix wall/reflect is not a good choice for what you want to do. Personally, I would use fix wall/harmonic to have a soft and repulsive force.
Using kspace is a further complication since internally, the system is changed to be periodic again, but the box is grown to 3x the size.
Since you just want to shrink the system and are not after perfect production quality results, I would recommend to switch to a cutoff coulomb and run a large (say 15 angstrom) cutoff.
You should not use fix deform here, but rather have a “moving” wall, where the wall position changed by using an equal style variable. Based on that variable you can determine the volume up to a small epsilon and thus estimate the density.
If you want to do something more complex, you could use the fix controller command — LAMMPS documentation and then use fix compute stress/atom in combination with compute reduce to determine the total stress and then divide by the volume estimated as mentioned above to get the pressure of the confined system. With fix controller you can determine the wall position (make sure that your computation contains boundaries, so the position stays within a reasonable value and couple not too strong).
You will see on the forum if you search for "pressure non-periodic that this topic is fairly often discussed. Take a look through those posts and try some of those suggestions.
Since you are simulating a system of polymers between graphene, you should adjust the system size by moving the graphene walls. The graphene-confined fluid interaction contributes significantly to the z-z pressure.