I am attempting to apply a pressure of 50 MPa to a 40 Å × 40 Å graphene sheet (containing 680 atoms) using the fix aveforce command. However, I’ve encountered some confusion regarding the implementation.
According to published literature, for a similar setup, the simulation was run by moving the piston 62 Å over a 10 ns period. When I convert the 50 MPa pressure to LAMMPS units, I obtain a value of 3.60E-2 (kcal/mol)/ų, which corresponds to a force of 57.57 (kcal/mol)/Å. To achieve this, I used the command: fix force piston aveforce NULL NULL 0.085
However, this causes the piston to move too quickly, preventing the simulation from reaching the desired 10 ns runtime.
Questions:
Could there be a mistake in my unit conversion for the force in LAMMPS?
Would using fix setforce be more appropriate for this scenario?
As an alternative, I resorted to using the velocity command to move the piston the required 62 Å within 10 ns, as this allowed me to estimate the travel time more precisely.
I would greatly appreciate any insights or suggestions you may have.
Are you sure that the pressure is imposed in your “published literature”? From that description, it seems that the displacement of the piston is imposed instead of the pressure.
Yes, their study included simulations at varying pressures (50–250 MPa), each run for 10 ns.
Thank you for your input. You’re absolutely right. Verifying calculations is an important part of the work. I’ve already checked my approach but wanted to seek additional perspectives in case I missed something.
Please double-check that you are using the exact same setup as in the literature, especially the number of carbon sheets used as a “piston”.
In my own experience such “pistons” do not work very well. As you have discovered, a few carbon sheets does not have enough mass density, and therefore enough inertia, to reliably compress most liquids, like room temperature water or ionic liquids. It is much easier for the pressure fluctuations in the liquid to deform and propel the carbon than vice versa.
I would perform a bulk NPT simulation of your liquid of interest and then fix the distance between walls that replicates said density in a suitable bulk sub-region of the confined liquid.
To assist fellow LAMMPS users and ensure accuracy, I’m sharing my pressure-to-force unit conversion calculations below. I’d greatly appreciate your insights and feedback to help identify any potential errors and improve these calculations.
Converting 50 MPa to (kcal/mol)/ų:
Start with 50 MPa = 5.0 × 10⁷ Pa = 5.0 × 10⁷ N/m² = 5.0 × 10⁷ J/m³.
So, as I understand you preffered to use the NPT for the bulk water and then insert carbon walls. but doesn’t it beak and disturb the relaxed structure?
The protocol you describe – bulk NPT simulation and then insert carbon walls – will certainly require re-equilibration to establish the electrolyte ordering near the walls.
I was not clear in my earlier post: I would use a bulk NPT simulation to establish the benchmark density, and then run an electrolyte + walls simulation from scratch.
I would include a fix ave/chunk measurement of the density profile across that simulation and determine (1) that a “bulk-like” region exists in the middle of the box where density fluctuations are minimal (2) that the wall spacing is set to a value such that the density in said “bulk-like” region is equal (within uncertainty) to the benchmark value.
But yes, my original post didn’t rule out the protocol you suggested which would indeed require re-equilibration.