LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] r-RESPA for pressure calculation
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] r-RESPA for pressure calculation


From: Pengyu Huang <hughh.py@...24...>
Date: Wed, 18 Apr 2018 01:41:02 +1000

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

On Tue, Apr 17, 2018 at 11:23 PM, Axel Kohlmeyer <akohlmey@...24...> wrote:
On Mon, Apr 16, 2018 at 10:38 PM, Pengyu Huang <hughh.py@...24...> wrote:
> 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.

>
> Regards,
> Pengyu
>
>
>
>
>
>
>
>
>
> On Tue, Apr 17, 2018 at 2:21 AM, Axel Kohlmeyer <akohlmey@...24...> wrote:
>>
>> 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.
>>
>> On Mon, Apr 16, 2018 at 9:40 AM, Pengyu Huang <hughh.py@...24...> wrote:
>> > 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
>> >
>> > On Mon, Apr 16, 2018 at 11:37 PM, Pengyu Huang <hughh.py@...24...>
>> > wrote:
>> >>
>> >> 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
>> >
>> >
>> >
>> >
>> > ------------------------------------------------------------------------------
>> > Check out the vibrant tech community on one of the world's most
>> > engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>> > _______________________________________________
>> > lammps-users mailing list
>> > lammps-users@...396...sourceforge.net
>> > https://lists.sourceforge.net/lists/listinfo/lammps-users
>> >
>>
>>
>>
>> --
>> Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
>> College of Science & Technology, Temple University, Philadelphia PA, USA
>> International Centre for Theoretical Physics, Trieste. Italy.
>
>



--
Dr. Axel Kohlmeyer  akohlmey@...43...4...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.

Attachment: E_for_r-RESPA_test2.jpg
Description: JPEG image

Attachment: Moving_average_E_for_r-RESPA_test2.jpg
Description: JPEG image

Attachment: P_for_r-RESPA_test2.jpg
Description: JPEG image

Attachment: Moving_average_P_for_r-RESPA_test2.jpg
Description: JPEG image