To calculate the pressure of a pure CO2 (TraPPE) system at NVT ensemble:
When I start with an initial configuration at liquid phase, after 5 ns simulation, sometimes, I get pressure values with 10-50 % error compared to the experimental values. However, when I start with an initial configuration at gaseous state (with known density) after 1 ns simulation, I can calculate the pressure very close to the experimental value but it takes a long CPU time to complete the run.
I even tried the saturation points and so the same P& T for both phases. According to the literature (MC simulation), TraPPE can predict the saturation vapour pressure for CO2.
Can anyone describe what changes do I need to make in my MD code, so I can predict precise values of pressure at liquid states?
Please remember that there are no mind readers here, so how should we know what changes you need to make since we only have some very vague descriptions of what you did?
To even begin to help we would need more specific information, such as:
what values of density you tried
what values of pressure you observed
what convergence you saw (such as a graph of pressure vs simulation time)
what settings you used (timestep, thermostat parameters, rigidity algorithm)
And most importantly, why you expect to be able to replicate MC results with an MD algorithm. Yes, in principle they sample the same distribution, but (1) they may not actually sample the same distribution if settings for one or the other technique are wrong (2) even if they sample the same distribution, the efficiency by which independent realisations are generated will differ between the techniques and give you problems with one or the other.
I started with an initial configuration with density of 18.0 mol/L and equilibrated the system with fix rigid/npt for 2ns at 300 K and 94 atm. After this step, I got an average density of 18.24 mol/L. Then I changed the size of the box to match the average volume that I calculated during the npt run and switched to fix rigid/nvt for 10 ns simulation. However, in the end, the pressure that I get is 78.72 atm! which is too understimated.
I was expecting values around 94 atm, the same that I fixed during the npt.
I even saw more deviation for other cases. This occurs for any configuration at the liquid state but for configurations at gaseous state I can get values very close to those that I had fixed during the npt (but in longer simulation time).
I tested other forcefields but I still see the same issue.
I followed two approaches to analyze the pressure data:
First method (as shown in the graphs sent before), I used the instantaneous pressure data that were printed every 100 fs, then averaged 10000 of these samples so that each 1 million simulation is equivalent to 1 block, so in total I have 10 blocks for 10 ns simulation time. The final value for pressure was then the average of the 10 blocks that is equal to 78.7194 atm with standard deviation of 0.95.
In the second method, I used the fix ave/time command to average the pressure values over every 20 fs and 50 samples per block (so the size of block is 1000 fs), and all are stored in a dat file until the end of 10 ns simulation.
In this case the average of all 10000 pressure data points (attachment) is 78.8063 with standard deviation of all data equal to 28.92.
@Hasan2023, for the NPT portion of your simulation, I see in the input script that you did an equilibration of 50000 timesteps before running the 2ns from which you seem to have gathered the data for calculating the equilibrium volume. Did you make a plot of volume versus time for these 50000 timesteps? I was thinking that maybe it was not enough for your system to equilibrate.
If it indeed your system has not equilibrated before you started gathering data for calculating average volume, it would make sense that this given average volume you are calculating is the incorrect one for the thermodynamic state of (300K, 94 atm). And then it would make sense why the NVT simulation at that volume and T = 300K is not giving you a value closer to 94 atm.
