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: Axel Kohlmeyer <akohlmey@...24...>
Date: Tue, 17 Apr 2018 14:18:22 -0400

On Tue, Apr 17, 2018 at 11:41 AM, Pengyu Huang <hughh.py@...24...> wrote:
> 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.

>
> 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@lists.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@...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@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.