rRESPA with shake – LAMMPS General Discussion

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

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

#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

#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

#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 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 200

#undump trjfile

#thermo 20000

Run the simulation

#run 1000000

Read more here: Source link