Hello LAMMPS community
I am trying to load data from standard LAMMPS data.spce file and then replicate to fill computational cell. However, my goal is to have a spherical computational cell I am planning to create with region and create_box commands. How to properly populate a spherical computation cell with SPC/E water molecules?
Many thanks
If you have a well equilibrated(!) periodic box of water, you can use the replicate command to build a sufficiently large system that that your sphere will fit in entirely. Then you can define a sphere region with “side out” setting and use the “delete_atoms” command to delete all water molecules (and their bonds/angles!) in that region. Finish with reset_atom_ids and reset_mol_ids, and perhaps the set command to reset all image flags to 0 and write it out with write_data.
Thank you for your response, Axel!
I am also curious if using a spherical box breaks pppm solver or not since it uses 3D grid for summation.
Best regards
You cannot have a spherical box. Even if you define a spherical region and then create a box from it, the box will be a cube. PPPM won’t work any other way regardless.
In other words, you can either have an isolated system but then cannot use a long-range solver (well, except for MSM, but that has its own set of problems). or you will be simulating an infinite grid of spheres that will “see” each other through PBC (as required by Ewald and PPPM). You can mitigate the error to some degree by increasing the vacuum area, but the computational cost (PPPM’s FFTs will cost you per grid point, whether there are atoms or not).
Thank you so much for clarification!
After I removed molecules from the outside of spherical region as you suggested, I try to integrate further with nvt but get an error message:
Out of range atoms - cannot compute PPPM
My script (LAMMPS version 09 Jan 2020) :
units real
atom_style full
boundary p p p
timestep 1.0
read_data data.spce
replicate 2 2 2
pair_style lj/cut/coul/long 9.0 9.0
pair_modify tail yes
kspace_style pppm 1.0e-4
pair_coeff 1 1 0.15539421659476232 3.16555789
pair_coeff * 2 0.0000 0.0000
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
bond_coeff 1 5000.00 1.000
angle_coeff 1 500.0 109.47
special_bonds lj/coul 0.0 0.0 0.5
neighbor 2.0 bin
neigh_modify every 1 delay 10 check yes
dump d_init all custom 500 spce_init.dump id type mol x y z
fix 1 all shake 0.0001 20 0 b 1 a 1
fix 2 all nvt temp 300.0 300.0 100.0
velocity all create 300 432567 dist uniform
thermo_style one
thermo 50
run 100
unfix 1
unfix 2
undump d_init
region to_del sphere 35.0 35.0 35.0 20 side out
delete_atoms region to_del bond yes mol yes reset_ids
set group all image 0 0 0
dump mydump all custom 500 spce.dump id type mol x y z
fix 3 all nvt temp 300.0 300.0 100.0
run 1000
I am not sure why this might happen. Based on the log file, the grid for PPPM is recalculated successfully after I remove the atoms. I get the error only when I start integrating.
It looks like you forgot to re-issue fix shake. it is needed for SPC/E as a model but also you cannot get a stable time integration at a 1fs timestep without.