Dear LAMMPS developers,
I am running the simulation for a system with flexible spc and flexible co2, and silica surfaces. I found the different pressures in the bulk water region between the Verlet integrator and r-RESPA integrator. I used a timestep of 0.25 fs for Verlet and the following set up for r-RESPA:
run_style respa 3 4 2 bond 1 angle 1 dihedral 1 improper 1 pair 2 kspace 3,
so that the bond and angle forces are calculated every 0.25 fs, pair forces at 1 fs and kspace at 2 fs. For kspace calculation (pair_style lj/cut/coul/long), pppm with accuracy 1e-6 was used.
The pressure in bulk water region (far away from the interfaces) with Verlet is around 12 MPa, while one with r-RESPA gives 115 MPa, which is 10x larger. The pressures in bulk co2 region from the two systems are statistically similar. I am aware when using large timestep for flexible water with Verlet integrator, the calculation of pressure would be incorrect, as I have tested something similar in http://lammps.sandia.gov/threads/msg66029.html (Although I got correct pressures for rigid spc at 1 fs…). So, would this means that the way of the pressure calculation in r-RESPA (at a larger timestep?) cause the problem?
Thank you for your time.
Regards,
Pengyu
Just forgot to mention in my previous email, the outermost timestep is 2 fs in r-RESPA, so that the bond and angle forces are calculated every 0.25 fs, pair forces at 1 fs and kspace at 2 fs.
Regards,
Pengyu
you also forgot to mention the LAMMPS version you are using and
whether you get the same issue, when you run r-RESPA with 0.25fs time
step and only one level or a multiplier of 1 (so that all levels run
with the same size timestep).
finally, you mention a multi-component system so, you should make
tests with a homogeneous system first.
axel.
Thank you, Axel, for your reply.
Sorry that I should have provided more details and further tests. The LAMMPS version I am using is the previous stable version 11Aug2017. I also did a quick check, it seems like the r-RESPA has not been changed.
I did a quick test with different integrators for the flexible SPC running in the NVT ensemble with “fix nvt” at T = 333.15K, restarting from a pre-equilibrated water box (which has 8287 SPC molecules) in the NPT ensemble at a targeted pressure of 20 MPa and T = 333.15 K. I have compared three integration methods: 1. r-RESPA with timestep = 0.25 fs and only one level for all forces calculations; 2 Verlet with timestep = 0.25 fs; 3 r-RESPA with timestep = 2 fs, with bond and angles calculated at 0.25 fs, pair at 1 fs and kspace at 2 fs (run_style respa 3 4 2 bond 1 angle 1 dihedral 1 improper 1 pair 2 kspace 3). In all cases, the Tdamp was 100 times of the timestep (or outermost timestep).
I have attached the plots for the total energy and the pressure against with time. The simulation was run 25 ps for the Cases 1 and 2, and 50 ps for the Case 3. It may be a bit short, but I think it already captures the problem I described. The Case 1 and 2 had similar results, while the Case 3 had higher total energy and pressure than the Cases 1 and 2. The pressure for Cases 1 and 2 were around 16 -18 MPa, respectively, while the Case 3 had an average pressure of 120 MPa. All cases had similar degrees of fluctuations, given the standard deviations were around 35-40 MPa. If there is no problem with the r-RESPA, would this means that I cannot use it for flexible water to get the correct pressure?
Regards,
Pengyu
Thank you, Axel, for your reply.
Sorry that I should have provided more details and further tests. The LAMMPS
version I am using is the previous stable version 11Aug2017. I also did a
quick check, it seems like the r-RESPA has not been changed.
that is not so easy to say. the respa.cpp and respa.h files only
contain the scheduling of the calculations. there are plenty of RESPA
specific functions scattered throughout other source files, which can
have an impact on the output.
I did a quick test with different integrators for the flexible SPC running
in the NVT ensemble with "fix nvt" at T = 333.15K, restarting from a
pre-equilibrated water box (which has 8287 SPC molecules) in the NPT
ensemble at a targeted pressure of 20 MPa and T = 333.15 K. I have compared
three integration methods: 1. r-RESPA with timestep = 0.25 fs and only one
level for all forces calculations; 2 Verlet with timestep = 0.25 fs; 3
r-RESPA with timestep = 2 fs, with bond and angles calculated at 0.25 fs,
pair at 1 fs and kspace at 2 fs (run_style respa 3 4 2 bond 1 angle 1
dihedral 1 improper 1 pair 2 kspace 3). In all cases, the Tdamp was 100
times of the timestep (or outermost timestep).
I have attached the plots for the total energy and the pressure against with
time. The simulation was run 25 ps for the Cases 1 and 2, and 50 ps for the
Case 3. It may be a bit short, but I think it already captures the problem I
described. The Case 1 and 2 had similar results, while the Case 3 had higher
total energy and pressure than the Cases 1 and 2. The pressure for Cases 1
and 2 were around 16 -18 MPa, respectively, while the Case 3 had an average
pressure of 120 MPa. All cases had similar degrees of fluctuations, given
the standard deviations were around 35-40 MPa. If there is no problem with
the r-RESPA, would this means that I cannot use it for flexible water to get
the correct pressure?
so far you have confirmed, that the r-RESPA formulation in LAMMPS
falls back on the verlet code even when using the various levels to
compute forces and stress.
the next step would be to determine systematically, where the pressure
shift is coming from and whether it is dependent on the timestep.
i suggest you start with using an outer time step of 0.5 fs and
run_style respa 2 2 bond 1 angle 1 dihedral 1 improper 1 pair 1 kspace 2
as well as with
run_style respa 2 2 bond 1 angle 1 dihedral 1 improper 1 pair 2 kspace 2
then time step of 1fs with respa 2 4, then time step 2fs with respa 2 8
this should allow you to identify whether the pressure shift is
dependent on either pair or kspace or both and whether it happens
suddenly or gradually.
axel.
Thanks for your suggestion, Axel.
Please find attached the results of total energy and pressure against time from more tests done with r-RESPA.
All bond and angles were calculated at every 0.25 fs and only outermost timestep was varied up to 1 fs for the pair and kspace calculation, as you suggested (I did not try 2 fs though). The results showed that the pressure shift was from the pair forces calculation with larger timesteps. The moving averages are also plotted so that it is easier to compare the changes in pressure. What do you think about the results? Would it be any problem with r-RESPA formulation in LAMMPS? Or simply we cannot simulate the flexible water with a larger timestep for correct pressure calculation?
Regards,
Pengyu
Thanks for your suggestion, Axel.
Please find attached the results of total energy and pressure against time
from more tests done with r-RESPA.
All bond and angles were calculated at every 0.25 fs and only outermost
timestep was varied up to 1 fs for the pair and kspace calculation, as you
suggested (I did not try 2 fs though). The results showed that the pressure
shift was from the pair forces calculation with larger timesteps. The moving
averages are also plotted so that it is easier to compare the changes in
pressure. What do you think about the results? Would it be any problem with
r-RESPA formulation in LAMMPS? Or simply we cannot simulate the flexible
water with a larger timestep for correct pressure calculation?
it is not my job to make your conclusions for you.
since you can single out the pairwise interactions there may be other
factors related to those, that you have not checked out yet, and that
may have an impact here. for example, the neighborlist update settings
or the cutoff(s).
axel.
Thanks Axel.
I did not see any dangerous builds from those tests so I think it may be unlikely due to the neighbor list update. Anyway, I will have a look at the code when I have time and see if I can figure out myself. Thanks for your help again.
Regards,
Pengyu
Hi Axel,
Sorry to bother you again. I have found a paper about pressure calculation for r-REPSA method: https://onlinelibrary.wiley.com/doi/epdf/10.1002/jcc.24731. It suggests that the pressure calculation should include the contribution from the forces calculated at the intermediate timestep as well (between innermost and outermost timesteps). It seems like the compute pressure command only calculated with the forces at every outermost timestep (correct me if I am wrong). I am not sure whether the method of pressure calculation they proposed would resolve the problem. If I would like to implement that method in LAMMPS, would it be lots of work? since I have little experience with c++.
Thank you for your time.
Regards,
Pengyu