Dears,
I am a beginner in molecular dynamics and I need your help to solve this problem.
I am running a simulation using water/surfactant/oil/surfactant/water with 28438 atoms at 323K and 98 atm, 3ns NPT to shrink the box (5x5x12 nm3), obtaining the right density and 2 ns NVT.
At NPT, the pressure is very good, near 100 atm without big fluctuations. But at NVT the pressure goes down, remaining in equilibrium in 50 atm.
My question is: how should I change the script to maintain the pressure close to 98 atm during NVT steps?
I tried do solve this problem reading the LAMMPS Manual and also here the Mailing Lists. But I was not able.
My Script:
pair_modify mix arithmetic tail yes kspace_style pppm 1.0e-4 minimize 1.0e-4 1.0e-6 100 1000 reset_timestep 0 fix SHAKE all shake 0.0001 20 0 b 1 3 9 11 a 1 neighbor 2.0 bin timestep 1.0 variable TK equal 323.15 variable PBAR equal 98.6923 velocity all create ${TK} 123456 rot yes mom yes dist gaussian variable nsteps equal 2000000 # 2000 # variable nptsteps equal 3000000 # 2000 # variable nhalf equal ${nsteps}/2 variable nprint equal ${nsteps}/2000 variable ndump equal ${nsteps}/100 variable xyzdump equal ${nsteps}/20 variable padval equal 7 # 5 variable nrdf equal 100 # 10 # variable Tdamp equal 100 # 10 # variable Pdamp equal 500 # 50 # # fix 1 all nvt temp ${TK} ${TK} ${Tdamp} fix 1 all npt temp ${TK} ${TK} ${Tdamp} iso ${PBAR} ${PBAR} ${Pdamp} variable xPress equal c_thermo_press[1] variable yPress equal c_thermo_press[2] variable zPress equal c_thermo_press[3] variable IFT equal (0.5*lz*(v_zPress-0.5*(v_xPress+v_yPress))*0.0101325) dump 2 all image ${xyzdump} nvt.*.jpg type type zoom 1.0 adiam 1.5 size 720 1280 dump_modify 2 pad ${padval} thermo_style custom step pe ke etotal temp press density vol lx ly lz pxx & pyy pzz pxy pxz pyz v_IFT thermo ${nprint} run ${nptsteps} unfix 1 fix 2 all nvt temp ${TK} ${TK} ${Tdamp} # fix 2 all npt temp ${TK} ${TK} ${Tdamp} iso ${PBAR} ${PBAR} ${Pdamp} run ${nsteps} write_restart restart.*.lmp write_data nvt.*.data
Thank you so much for your attention.
Kind regards,
For a single phase system of fixed particle numbers, it is a thermodynamic consequence that specification of T and P fixes V and vice versa. In our group, for accurate property calculations at specified (NPT), we first do 10 independent (different starting configurations) NPT runs and use the resulting average of the box volumes to do an NVT run. Each NPT run consists of an equilibration stage followed by a production run. If the standard deviation of the volumes is too large, then you must go back and make appropriate changes to your NPT run. This procedure avoids the requirement for a barostat in a direct NPT run, resulting in improved precision.
Of course, you must also use the correct formulation for the property you want to calculate in the NVT run. For example. volumetric properties are generally not a problem, but the appropriate expression must be used for the chemical potential.
I hope this helps.
You can't force an NVT simulation to result in a specified pressure. It will equilibrate to a value based on the temperature and the density of the system.
But I would also expect that first equilibrating with NPT and then running NVT should result in a more or less equal value of pressure (Though the common practice is usually in reverse, first NVT then NPT). I can think of two reasons why you don't get equal pressure:
1) The system has not equilibrated in the NPT run. Have you looked at some measure to ensure reaching equilibrium? Easiest to check is the potential energy. It should plateau with fluctuations around equilibrium values. The fact that the pressure has reached equilibrium doesn't mean the whole system has equilibrated. Because the barostat takes care of that.
2) It might be that something has happened in the second part of the run. Phase segregation for example, which would change the equilibrium properties of the system. You can check that by looking at the trajectories of the two segments. Is the configuration the same in both?
Mohsen has some good questions and advice. I would also add that 100 atm is fairly small and I would expect some fluctuations if you are trying to hold it. The script seems ok from what I can see. A few more thoughts.
1. How bad are the fluctuations for such a small pressure? You could consider changing the pressure damping parameter to try to get the system to react faster or slower depending on how quickly it is equilibrating.
2. As alluded to, the potential may have something to do with it. If you tried the script with a simple copper potential, does this work? Is the pressure more stable?
3. System size can affect this. If you have the luxury of doubling the size of the system, does this discrepancy get better? The other thing would be the timestep. It seems reasonable, but I don't know your potential.
Best of luck.
In NVT and NPT, N is the number of molecules, V is the volume of the system, T is the temperature and P is the pressure. P and V are inversely related. So, NVT is like cooking in a pressure cooker, while NPT is like cooking with the lead open. If you define both for an ensemble, the system will never equilibrate. You cannot basically put the lead on and off at the same time.
For a single phase system of fixed particle numbers, it is a thermodynamic consequence that specification of T and P fixes V and vice versa. In our group, for accurate property calculations at specified (NPT), we first do 10 independent (different starting configurations) NPT runs and use the resulting average of the box volumes to do an NVT run. Each NPT run consists of an equilibration stage followed by a production run. If the standard deviation of the volumes is too large, then you must go back and make appropriate changes to your NPT run. This procedure avoids the requirement for a barostat in a direct NPT run, resulting in improved precision.
Of course, you must also use the correct formulation for the property you want to calculate in the NVT run. For example. volumetric properties are generally not a problem, but the appropriate expression must be used for the chemical potential.
I hope this helps.
Dears, thank you so much. I changed the volume and the kspace and now the pressure in close to 100 atm.
I really appriciated your answers.
- Badges
- Science topic