i generated a lammps system using charmmgui for a ~10000 atoms protein with 40 angstroms of water around it so that my protein has plenty of space to swim around.
im trying to use fix wall/harmonic at the edges of the simulation box to avoid parts of the protein wrapping around periodic boundary conditions. one problem with the input generated by charmmgui is that some water molecules have coordinate(s) outside the simulation box.
i tried delete_atoms to remove water molecules that might either straddle, or are too close to, the edges of the simulation box so i can put a fix wall/harmonic there. however i’m getting:
ERROR on proc 0: Bond atom missing in image check (src/domain.cpp:794)
even though i use both bonds yes and mol yes in my delete_atoms command.
is this a bug ?
since my original system is >200mb, i made a minimal example with a small protein (PDB 6DNY) with 68 atoms and 20A of water around it. the commands i added to default charmmgui generated minimization script are enclosed by # <<< and # >>>. full inputs in zip file below.
6dny-min.in:
echo none
variable dcdfreq index 50
variable outputname index step4.0_minimization
units real
boundary p p p
newton off
pair_style lj/charmmfsw/coul/long 10 12
pair_modify mix arithmetic
kspace_style pppm 1e-6
atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmmfsw
special_bonds charmm
improper_style harmonic
timestep 1
fix cmap all cmap charmmff.cmap
fix_modify cmap energy yes
read_data step3_input.data fix cmap crossterm CMAP
neighbor 1 bin
neigh_modify delay 5 every 1
velocity all create 303.15 42273 dist gaussian
include restraints/constraint_angletype
fix 1 all nvt temp 303.15 303.15 100.0
shell sed -e "s/\$bb/1.0/g" -e "s/\$sc/0.1/g" step3_input.col > restraints/${outputname}.col
fix restraint all colvars restraints/${outputname}.col output ${outputname}
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
variable edge equal 4
region edges block $(xlo+v_edge) $(xhi-v_edge) $(ylo+v_edge) $(yhi-v_edge) $(zlo+v_edge) $(zhi-v_edge) side out
delete_atoms region edges compress no bond yes mol yes
write_data 6dny-min-after-delete_atoms.data
#change_box all x delta -8 8 y delta -8 8 z delta -8 8 units box
#reset_atoms image all
#fix 4 all wall/harmonic &
xlo $(xlo+4) 1000.0 0.0 4.0 xhi $(xhi-4) 1000.0 0.0 4.0 &
ylo $(ylo+4) 1000.0 0.0 4.0 yhi $(yhi-4) 1000.0 0.0 4.0 &
zlo $(zlo+4) 1000.0 0.0 4.0 zhi $(zhi-4) 1000.0 0.0 4.0 &
pbc yes units box
#fix_modify 4 energy yes
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
thermo ${dcdfreq}
thermo_style custom step time xlo xhi ylo yhi zlo zhi etotal pe ke temp press ebond eangle edihed eimp evdwl ecoul elong temp vol
dump 1 all dcd ${dcdfreq} ${outputname}.dcd
dump_modify 1 unwrap yes
min_style cg
minimize 0.0 1.0e-8 10000 20000
write_dump all custom ${outputname}.dump id type x y z vx vy vz ix iy iz
write_data ${outputname}.data
i’m aware of the “WARNING: System is not charge neutral, net charge = -5” but this only happens in the minimal example since some of the ions are close to the edges and get deleted by coincidence. this doesn’t happen in my original much larger system.
It would probably be simpler to just edit the data file and change the lo/hi values of the box to include all atoms in any non-periodic direction. If you want to add walls, it may be a good idea to add a little extra space, so you don’t get too large forces on the atoms close to the box boundaries when adding the wall(s). You can safely do that, since the image flags for all atoms are zero (the default when not specified) and thus the changing of the box dimensions does not change the positions (as it would for non-zero image flags).
i prefer not to manually find the min/max of x,y,z and edit the data file because my systems are typically >1M atoms. since this is a common bug from charmmgui generating atom positions outside the simulation box, i want a reproducible walkaround in lammps code to add to default charmmgui generated minimization script.
first i tried
change_box all x final $(min(x[i])-3) $(max(x[i])+3) &
y final $(min(y[i])-3) $(max(y[i])+3) &
z final $(min(z[i])-3) $(max(z[i])+3) &
units box
but i got:
ERROR: Invalid special function in variable formula (src/variable.cpp:4470)
then i tried
compute xmin all reduce min x
change_box all x final $(c_xmin-3) units box
but i got
ERROR: Variable formula compute cannot be invoked before initialization by a run (src/variable.cpp:1523)
so then
run 0
compute xmin all reduce min x
change_box all x final $(c_xmin-3) units box
and i still got
WARNING: No fixes with time integration, atoms won’t move (src/verlet.cpp:60)
compute xmin all reduce min x
change_box all x final $(c_xmin-3) units box
ERROR: Variable formula compute cannot be invoked before initialization by a run (src/variable.cpp:1523)
im ok with atoms not moving. it’s actually what i want to catch the atoms outside the simulation box before they get non-zero image flags
But also, if you know CHARMMGUI makes a 40A box, why not just manually set the box size to 44A or some similar margin? It’s a sed one-liner:
sed -i 's/.*\([xyz]\)lo.*/-2 42 \1lo \1hi/' myfile.data
(or whatever your start and stop numbers are in each dimension) and you know beforehand exactly what box size to use in LAMMPS too.
(You could also use change_box in LAMMPS to the same effect – the trick is to compute property/atom xu yu zubeforerun 0 and then use set to set the wrapped coordinates in the new, bigger box to their unwrapped previous values.)
For this to work you need to put compute xmin ...beforerun 0. I think.
The way you’ve set it up in the post quote, LAMMPS does a run 0, including calculating all computes* and variables it knows about, and then you tell it to compute xmin ... – by which point you need another run 0.
*The compute has to have a “consumer” that uses it to make an output, so you also need something like
thermo_style custom c_xmin
to tell LAMMPS it really, really needs to run that compute, I think.
that’s what i did based on @mkanski’s tip about bound group function of variable. For future googlers, here’s the solution:
change_box all &
x final $(bound(all,xmin)-5) $(bound(all,xmax)+5) &
y final $(bound(all,ymin)-5) $(bound(all,ymax)+5) &
z final $(bound(all,zmin)-5) $(bound(all,zmax)+5) &
units box
variable xmin equal $(xlo)
variable xmax equal $(xhi)
variable ymin equal $(ylo)
variable ymax equal $(yhi)
variable zmin equal $(zlo)
variable zmax equal $(zhi)
fix 4 all wall/harmonic &
xlo v_xmin 1000 0 5 xhi v_xmax 1000 0 5 &
ylo v_ymin 1000 0 5 yhi v_ymax 1000 0 5 &
zlo v_zmin 1000 0 5 zhi v_zmax 1000 0 5 &
pbc yes units box
However this was a futile attempt to implement walls and equilibrate with NPT because “it is odd to do NPT on a dimension where you also have a wall. The two commands are competing with each other.” as @sjplimp said years ago. I got the classic
when running the charmmgui equilibration script with fix npt.
In the end, i went back to the original charmmgui minimization script with no changes even if some water and ion particles are outside simulation box and get non-zero image flags.
My protein is solvated with a margin of 40A water around it. My goal was to avoid the protein wrapping itself around PBCs
# -- DUMP ALL EXCEPT WATER MOLECULES (PROTEIN AND IONS) --
group solv type 13 41
group protein_ions subtract all solv
dump 1 protein_ions dcd 1 2q3z-extra-equil.dcd
# ----- FIX WALL PROTEIN ONLY -----
group protein_only subtract protein_ions k na ca mg cl
variable xmin equal $(xlo)
variable xmax equal $(xhi)
variable ymin equal $(ylo)
variable ymax equal $(yhi)
variable zmin equal $(zlo)
variable zmax equal $(zhi)
fix 4 protein_only wall/harmonic &
xlo v_xmin 1000 0 5 xhi v_xmax 1000 0 5 &
ylo v_ymin 1000 0 5 yhi v_ymax 1000 0 5 &
zlo v_zmin 1000 0 5 zhi v_zmax 1000 0 5 &
pbc yes units box
Please note that fix wall needs to have variable positions to stay at the edges of the simulation box as fix npt continuously changes the volume to adjust pressure.
open 4pyg-3mM.psf coords 4pyg-3mM.dcd
rainbow residues target c
surface #1 transparency 62 color gray replace true
surface enclose :277,335,358 color white replace false
show :ASP,GLU atoms
color @cal deep pink
