Usage issue of fix shake

Please reply to the list, not to me, thanks.

Using fix shake while minimising is an error, so make sure it is disabled. I see now that shake generates an error in my version of LAMMPS if used with minimize, so I guess it could be something else going wrong. Again, what does your system look like after minimisation? Is it still reasonable? How large is the final energy? Please send an input script and data files so people can take a look.

Thank you. I think it is the problem of minimization. This is what I did to keep both minimization and fix shake.
Firstly, let it run for a few steps to minimize and dump the xyz information. Then update the datafile with the new xyz coordinates and add fix shake in script.

Best wishes
Xiaoyu

Thank you. I think it is the problem of minimization. This is what I did
to keep both minimization and fix shake.
Firstly, let it run for a few steps to minimize and dump the xyz
information. Then update the datafile with the new xyz coordinates and add
fix shake in script.

​why so complicated? LAMMPS has a "write_data" command.​ that does this in
one go and has the additional benefit of producing a properly formatted
data file.

in any case, unless you provide a simple and small(!) test example that
reproduces the segfault and you have verified that the segfault still
happens with the latest LAMMPS patchlevel (26 Oct 2015 today) nobody will
take a serious look. chances are high that you are doing something wrong,
but since nobody here can do mind reading or has a crystal ball, we need
hard and verifiable data, not vague descriptions and guesses.

axel.

I do have some issue about my model. Actually I do not want the bond to be rigid. But problems keep coming with all efforts I have done. I can show you some of my command lines. I also attach the datafile.

My simulation box is like a adsorption of water by using inorganic membrane like zeolite. I have tried several models simulation the dynamics by LAMMPS. For those with not very strong polarity of the membrane ( Al/Si=0 or 0.6 O=-0.21 or -0.4), it works well. While for those with strong polarity (Al/Si=1.5, O=-1.3), it collapses at very first few steps. ( with the error of "out of atoms” or “missing bonds” or “missing atoms”)

My method is tethering those membrane atoms and make the water molecules moving freely. I use lj/cut/coul/long as pair style for whole system and using the SPC/E potential for water.

I tried to shrink the tilmestep, but it did not work. When I turn off the coul interaction, the system could run well. Can anyone tell me is the algorithm i using is wrong or those parameters are wrong. I tried ewald and pppm both, but neither worked. The parameters are proved to work in others’ essays. So i do not think they are wrong. But they are using a combination of buckinghum and lj as the potential model. Could that be the problem? And if I want to do simulation in my own method, is there a way to fix the problem?

pair_style lj/cut/coul/long 12
bond_style harmonic
angle_style harmonic
kspace_style pppm 1.0e-4

create geometry

lattice fcc 1.0
read_data data.water_mem

#LJ Potentials
#1 = Si_zeolite
#2 = Al_zeolite
#3 = O_zeolite
#4 = Na_zeolite
#5 = O_water
#6 = H_water

pair_coeff 1 1 0.00000 0.000
pair_coeff 2 2 0.00000 0.000
pair_coeff 3 3 0.19870 3.400
pair_coeff 4 4 0.03974 3.200
pair_coeff 5 5 0.15540 3.170
pair_coeff 6 6 0.00000 0.000
pair_modify mix arithmetic

define groups and regions

group Si_zeolite type 1
group Al_zeolite type 2
group O_zeolite type 3
group Na type 4
group O_water type 5
group H_water type 6

group zeo union Si_zeolite Al_zeolite O_zeolite
group mem union Si_zeolite Al_zeolite O_zeolite Na
group mem1 id <> 1 352
group mem2 id <> 401 752
group water union O_water H_water

region center block 24.2775 113.9088 INF INF INF INF
region 1 block INF 8.0 INF INF INF INF
region 2 block 130.1863 INF INF INF INF INF
region out union 2 1 2
region 1m block 8.0 24.2775 INF INF INF INF
region 2m block 113.9088 130.1863 INF INF INF INF
region membrane union 2 1m 2m

variable Wcount_c equal count(O_water,center)
variable Wcount_m equal count(O_water,membrane)
variable Wcount_o equal count(O_water,out)

velocity mem create 293 289859
velocity water create 293 289859

fix 1 zeo spring/self 5000

fix 21 all nve
fix 22 mem temp/rescale 1 {TEMP} {TEMP} 0.05 1.0
fix 23 water temp/rescale 1 {TEMP} {TEMP} 0.05 1.0

compute Tmem mem temp
compute Twater water temp

dump 3 all custom {NDUMP} dump.{SYS}.${TEMP}.lammpstrj id type x y z

restart {NRESTART} restart.*.{SYS}.${TEMP}

output density profiles

fix 31 water ave/spatial 1 {NPROFOUT} {NPROFOUT} x lower 1.00555 density/mass density/number ave running file water.{SYS}.{TEMP}
fix 32 zeo ave/spatial 1 {NPROFOUT} {NPROFOUT} x lower 1.00555 density/mass density/number ave running file zeo.{SYS}.{TEMP}

output molecule number profiles

fix 41 O_water ave/time 1 {NPROFOUT} {NPROFOUT} v_Wcount_c v_Wcount_m v_Wcount_o ave one file water.number.${TEMP}

Run

thermo_style custom step dt temp c_Tmem c_Twater pe v_Wcount_c v_Wcount_m v_Wcount_o cpu
timestep {TIMESTEP} thermo {NTHERMOOUT}

minimize 1e-4 1e-6 1000 10000
reset_timestep 0
run ${NRUN}

data.water_mem (332 KB)

I do have some issue about my model. Actually I do not want the bond to be
rigid. But problems keep coming with all efforts I have done. I can show
you some of my command lines. I also attach the datafile.

My simulation box is like a adsorption of water by using inorganic
membrane like zeolite. I have tried several models simulation the dynamics
by LAMMPS. For those with not very strong polarity of the membrane (
Al/Si=0 or 0.6 O=-0.21 or -0.4), it works well. While for those with strong
polarity (Al/Si=1.5, O=-1.3), it collapses at very first few steps. ( with
the error of "out of atoms” or “missing bonds” or “missing atoms”)

My method is tethering those membrane atoms and make the water molecules
moving freely. I use lj/cut/coul/long as pair style for whole system and
using the SPC/E potential for water.

I tried to shrink the tilmestep, but it did not work. When I turn off the
coul interaction, the system could run well. Can anyone tell me is the
algorithm i using is wrong or those parameters are wrong. I tried ewald and
pppm both, but neither worked. The parameters are proved to work in others’
essays. So i do not think they are wrong. But they are using a combination
of buckinghum and lj as the potential model. Could that be the problem? And
if I want to do simulation in my own method, is there a way to fix the
problem?

​you are arbitrarily changing parts of your force field. as i keep
mentioning, MD is GI-GO (= garbage in, garbage out). force fields are
carefully parameterized so that there is a balance between the difference
contributions. if you are changing parts of that arbitrarily or mixing and
matching different pieces, chances are that things are going out of whack
eventually and then you cannot run a stable time integration anymore
because your force field as such doesn't support a stable structure.

you are not making a convincing case that you are actually using parameters
from other publications correctly. at any rate, what you describe seems to
be more of a science problem than a LAMMPS problem. if you can prove(!)
that you have set up a simulation *exactly* as in a corresponding reference
paper with consistent and correct parameters and that LAMMPS cannot run
this input, we'll take a second look. for now we have to assume that you
are just using bad parameters and then there is no surprise that your
simulation cannot work.

pair_style lj/cut/coul/long 12
bond_style harmonic
angle_style harmonic
kspace_style pppm 1.0e-4

​this is incomplete and thus useless.​

​axel.​

That is exactly I am trying to do right now. Those are all old publications decades ago. They were using dated algorithms and they did not described their simulation part very clearly. So it is not easy to reproduce their results. And their works were used in their research field which is not very similar to mine so that is why I wanted to combine their parameters with my model. But so far, it seems that my model is not compatible with their parameters.

Thank you.

That is exactly I am trying to do right now. Those are all old
publications decades ago.

​please note that i *am* one of those "decades old" people that did
simulations decades ago when i was in grad school and published results and
that i was using parameters that were decades ago at that time. not so much
has changed since. ever parameterizations from the 1970s get the
fundamentals right, the differences to newer parameters sets are rather
small and tweaks to improve details.

They were using dated algorithms
​.

dated algorithms for what? most of what you are using in terms of
algorithms *is* 20 years old and older and i'd rather call it "tried and
tested". if there would have been significantly better models as a
comparable computational cost, that is what people would be using, but most
of the classical MD technology in use today dates from the 1980s to 1990s.
old doesn't mean automatically bad.


and they did not described their simulation part very clearly. So it is
not easy to reproduce their results.

​ ...and i'd also contend that many publications nowadays are much more
"cavalier" and much less specific and detailed about how to set up the
presented simulations and how to prepare systems than what they were 20-30
years ago. people had to be much more thorough and careful, since the
amount of computational resources were much less and one had to make very
certain at every step, that things were correct and worth the effort.

And their works were used in their research field which is not very
similar to mine so that is why I wanted to combine their parameters with my
model. But so far, it seems that my model is not compatible with their
parameters.

​that is a *fundamental* property of classical force fields​. *all*
parameters have to be balanced with each other. as i already mentioned,
once cannot simply mix and match parameters that have different underlying
strategies for deriving partial charges and solvation energies. that is why
people that parameterize new components for force fields like CHARMM or
Amber still use outdated quantum chemical methodology, in order to make all
additional parameters consistent with the existing parameters. if you'd
change, e.g. to modern state-of-the-art density functionals, then you'd
basically have to redo *all* existing parameters to make them consistent
with each other. a lot of the empirical balancing of non-bonded parameters
or the force constant for bonded interactions would have to be recomputed
as well.

what you describe about how your system is behaving for some choices of
parameters, is rather drastic and an indication of correspondingly
drastically inadequate parameters.

axel.

Okay, all right. If there is any further problems, I will email the lammps mail list. Thank you all the same for your advice.

Take care
Xiaoyu