rRESPA with shake

LAMMPS LAMMPS General Discussion
1

Dear lammps user,
I am using lammps March 3 , 2020 version. I have an ionic system which consist of uranyl , nitrate and water. Water was made rigid molecule by shake. I want to implement rRESPA to speed up the calcalculation where bond , angle and hybrid-1 will be computed by level 1 and hybrid-2 , kspace will be computed by level 2. Without rRESPA in NVE ensemble energy conservation is fine but with rRESPA system is not stable. I have attached my input file

#LAMMPS Input file

aqueous uranyl nitrate

Intialization

units real
dimension 3
boundary p p p
atom_style full
atom_modify map array sort 500 2.0

Atom Definition

#read_restart restart.500000.equil

read_data data.unw
replicate 1 1 1

bond_style harmonic
bond_coeff 1 1000.95000000000 1.80000000000000 #U-OU
bond_coeff 2 600.57000000000 1.26000000000000 #N-ON
bond_coeff 3 700.000000000000 1.08000000000000 #O-H

angle_style harmonic
angle_coeff 1 300.280000000000 180.000000000000 #Ou-U-Ou
angle_coeff 2 300.280000000000 120.000000000000 #On-N-On
angle_coeff 3 700.000000000000 109.470000000000 #H-O-H

#pair_style lj/long/coul/long long long 10.0
pair_style hybrid/overlay lj/cut 15.0 coul/long 15.0
kspace_style pppm 1.0e-4
#kspace_style ewald/disp 1.0e-4

pair_coeff * * coul/long

pair_coeff 1 1 lj/cut 2 0.0273 3.3500
pair_coeff 1 2 lj/cut 2 0.1094 3.1000
pair_coeff 1 3 lj/cut 2 0.0661 3.2350
pair_coeff 1 4 lj/cut 2 0.0640 3.1450
pair_coeff 1 5 lj/cut 2 0.0651 3.2580
pair_coeff 1 6 lj/cut 2 0.0000 1.6750
pair_coeff 2 2 lj/cut 2 0.4384 2.8500
pair_coeff 2 3 lj/cut 2 0.2650 2.9850
pair_coeff 2 4 lj/cut 2 0.2566 2.8950
pair_coeff 2 5 lj/cut 2 0.2610 3.0080
pair_coeff 2 6 lj/cut 2 0.0000 1.4250
pair_coeff 3 3 lj/cut 1 0.1601 3.1200
pair_coeff 3 4 lj/cut 1 0.1551 3.0300
pair_coeff 3 5 lj/cut 1 0.1578 3.1430
pair_coeff 3 6 lj/cut 1 0.0000 1.5600
pair_coeff 4 4 lj/cut 1 0.1501 2.9400
pair_coeff 4 5 lj/cut 1 0.1527 3.0530
pair_coeff 4 6 lj/cut 1 0.0000 1.4700
pair_coeff 5 5 lj/cut 1 0.1554 3.1660
pair_coeff 5 6 lj/cut 1 0.0000 1.5830
pair_coeff 6 6 lj/cut 1 0.0000 0.0000

group uranyl type 1:2
group nitrate type 3:4
group water type 5:6

Settings

#velocity uranyl create 300.0 4928459 rot yes dist gaussian
#velocity nitrate create 300.0 4928458 rot yes dist gaussian
#velocity water create 300.0 4928457 rot yes dist gaussian

#velocity guest scale 200.0

fix 1 all nve
fix 2 all temp/rescale 10 300.0 301.0 0.1 1.0
fix 3 water shake 0.0001 10 0 b 3 a 3
unfix 3

minimize 1.0e-6 1.0e-6 1000 1000

fix 3 water shake 0.0001 10 0 b 3 a 3

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

timestep 10.0
run_style respa 2 10 hybrid 1 2

Output

#Equilibration 100 ps

#restart 200000 restart.*.equil

thermo_style custom step temp etotal pe ke ecoul elong evdwl ebond eangle vol lx ly lz

thermo 1

Run the simulation

run 100

#Step 1 1 ns vacf and msd dump
reset_timestep 0
unfix 2
restart 500000 restart.*.equil

dump 1 all xyz 10 HISTORY.xyz
dump_modify 1 element U OU N ON OW HW
#dump trjfile all custom 1000 dump.lammpstrj id mol type x y z
#dump dumpXYZ guest xyz 1000 waterpos1ps.xyz
dump dump1XYZ uranyl custom 10 uranylpos10fs.xyz id x y z
dump dump2XYZ nitrate custom 10 nitratepos10fs.xyz id x y z
dump dump3XYZ water custom 10 waterpos10fs.xyz id x y z
dump dumpV1XYZ uranyl custom 10 uranylvel10fs.dat id vx vy vz
dump dumpV2XYZ nitrate custom 10 nitratevel10fs.dat id vx vy vz
dump dumpV3XYZ water custom 10 watervel10fs.dat id vx vy vz
#dump dumpFXYZ guest custom 1000 ethenefrc1ps.dat id fx fy fz

thermo 2

Run the simulation

run 200

#undump trjfile

#thermo 20000

Run the simulation

#run 1000000

2

No surprise here. Your input is quite bogus.

  • for this kind of setup you will never get a 1:10 split between the levels because of kspace which has to cover all atoms. There will already by an impact on energy conservation for a 1:2 split and you may get away with 1:3 or 1:4 if you have a strong thermostat, but already the fact that for the fast atoms short-range coulomb is changing while long-range is not, is giving you inconsistencies.
  • you are misunderstanding how the hybrid division is happening. It is between coul/long and lj/cut, not between two groups of lj/cut. This also leads to bogus results. Because of having lj/cut and coul/long as separate sub-styles you are slowing down your calculation, because you have to loop over the same atoms in the neighbor lists twice compared to using lj/cut/coul/long.
  • your pair_coeff statements are obviously wrong. Their three arguments will be read as epsilon, sigma, cutoff. The first argument is obviously meant for assigning to different substyles, but for that you would have to define lj/cut in the pair style line twice.

Overall, since you are already using fix shake, there is very little to gain from r-RESPA. An outer timestep of 10fs is quite impossible in my opinion.

In general, what you are trying to do is an ill-conceived idea and incorrectly applied.
I have already explained how unlikely any speedup will be for this kind of setup in a previous response.