Dear LAMMPS users,
I would like to apparently having two group of water molecules with the same interaction because I will have to treat the two groups of water molecules slightly differently at some point,
, i.e, both type 1 and type 3 are O atoms and both type 2 and type 4 are H atoms.
Is it possible with the pair_style hybrid/overlay command in LAMMPS? For example, with something like:
pair_style hybrid/overlay lj/cut/tip4p/long 1 2 1 1 0.15 8.5 8.5 lj/cut/tip4p/long 3 4 1 1 0.15 8.5 8.5
pair_coeff 1 1 lj/cut/tip4p/long 0.155 3.1536
pair_coeff * 2 lj/cut/tip4p/long 0.0000 0.0000
pair_coeff 3 3 lj/cut/tip4p/long 0.155 3.1536
pair_coeff 1 3 lj/cut/tip4p/long 0.155 3.1536
pair_coeff * 4 lj/cut/tip4p/long 0.0000 0.0000
kspace_style pppm/tip4p 1.0e-5
The script above doesn’t work because the log.lammps output is:
ERROR: Expected integer parameter in input script or data file (…/pair_hybrid_overlay.cpp:58)
I believe it is not due to any problem in the data file because if I assign other pair potential as the second sub-style there is no running errors.
Best Regards
Lunna
The pair_coeff syntax is different when you use the same pair style twice. How should LAMMPS otherwise be able to match the coefficients with the styles?
Axel
Hi,
there might be another problem. Even if you correct the syntax you’re
simulations will provide incorrect results, because the kspace solver will
recognize only either type 1 or type 3 as TIP4P oxygen, but not both.
Maybe your should think about whether you can define the groups in some other
way. Will defining groups by atom IDs or the initial configuration work for you?
Best, Rolf
"Maybe your should think about whether you can define the groups in some other
way. Will defining groups by atom IDs or the initial configuration work for you?"
Dear Rolf/LAMMPS Users,
I was thinking about something similar, however, I have not yet figured out a way to do that.
I wish to switch the tip4p interaction for one group of water molecules gradually from off to on while keeping the other group fully interacting tip4p so I will need to have two independent but the same interactions.
I might be wrong but I feel that even if I group via atom type/id I can not have two group of water molecules, because the format of the tip4 pair_style in LAMMPS required otype and htype but not group ids.
lj/cut/tip4p/long args = otype htype btype atype qdist cutoff (cutoff2)
Best
Lunna
Hi,
if you only need to switch of vdW interactions, a combination of the pair styles lj/cut
and tip4p/long might be solution to your problem with
pair_style hybrid/overlay lj/cut 8.5 tip4p/long …
Otherwise I don’t know how what you are trying to do is possible with LAMMPS
without changing the source code.
Best, Rolf
Hi,
if you only need to switch of vdW interactions, a combination of the pair
styles lj/cut
and tip4p/long might be solution to your problem with
pair_style hybrid/overlay lj/cut 8.5 tip4p/long ..
Otherwise I don't know how what you are trying to do is possible with LAMMPS
without changing the source code.
one can also do TIP4P water as an explicit 4-site model using fix
rigid and lj/cut/coul/long and regular PPPM. the downside is that
modeling rigid water with fix rigid requires a much shorter time step
than when using fix shake, but - as the saying goes - there ain't no
such thing as a free lunch... 
axel.
Dear Axel/LAMMPS Users,
Thank you so much for the help!
"one can also do TIP4P water as an explicit 4-site model using fix
rigid and lj/cut/coul/long and regular PPPM. the downside is that
modeling rigid water with fix rigid requires a much shorter time step
than when using fix shake, but - as the saying goes - there ain’t no
such thing as a free lunch… "
I am trying to use “fix rigid + explicit 4-site tip4p models”.
After separating the O mass to O site and O charge to the massless M site, it seems that LAMMPS does not allow me to use massless particles because it complains about "Invalid mass value (…/atom.cpp).
Hence I have assigned a tiny mass (e.g. 0.00001) to the massless M site, which runs fine and in short run produces reasonable potential energy comparable to the original tip4p model in LAMMPS if I do not use “fix recenter”, while if I use “fix recenter” following the fix rigid/npt command, the dynamics is wrong, (e.g. potential energy/radial distribution function is wrong and the temperature is also wrong and not oscillating at the desired T=298K)
Does it mean that I cannot constrain the centre of mass of the whole system if I use fix rigid? Is there a way to use both fix rigid and constrain the centre of mass of the whole system?
Best Wishes
Lunna
Dear Axel/LAMMPS Users,
Thank you so much for the help!
"one can also do TIP4P water as an explicit 4-site model using fix
rigid and lj/cut/coul/long and regular PPPM. the downside is that
modeling rigid water with fix rigid requires a much shorter time step
than when using fix shake, but - as the saying goes - there ain't no
such thing as a free lunch... "
I am trying to use "fix rigid + explicit 4-site tip4p models".
After separating the O mass to O site and O charge to the massless M site,
it seems that LAMMPS does not allow me to use massless particles because it
complains about "Invalid mass value (../atom.cpp).
Hence I have assigned a tiny mass (e.g. 0.00001) to the massless M site,
yes, that is the way to do it.
which runs fine and in short run produces reasonable potential energy
comparable to the original tip4p model in LAMMPS if I do not use "fix
recenter", while if I use "fix recenter" following the fix rigid/npt
command, the dynamics is wrong, (e.g. potential energy/radial distribution
function is wrong and the temperature is also wrong and not oscillating at
the desired T=298K)
why do you need to use fix recenter in the first place? this is almost
always an indication of something going wrong. it is dangerous, too,
because it can mask if your system is becoming subject to the
so-called "flying icecube syndrome".
the cleanest way to keep a system from building up a COM drift is to
use a weak restraint via fix spring in tether mode. in some cases also
using fix momentum might do the trick and would be preferable to fix
recenter. the latter is mostly meaningful if you would have to
perfectly recenter the trajectory for analysis in postprocessing, but
this is very rarely needed.
axel.
Dear Axel/LAMMPS Users
"the cleanest way to keep a system from building up a COM drift is to
use a weak restraint via fix spring in tether mode. in some cases also
using fix momentum might do the trick and would be preferable to fix
recenter. the latter is mostly meaningful if you would have to
perfectly recenter the trajectory for analysis in postprocessing, but
this is very rarely needed"
Thanks so much for the great help on this thread!
Therefore I have been using “fix momentum” together with my explicit 4-site tip4p water model with lj/cut/coul/long and regular pppm with “fix rigid molecule”
“fix momentum” helps a lot.
However COM drift still happen within ~4ns for my explicit 4-site tip4p water models, while no COM drift occurs for the default tip4p water models in LAMMPS within ~18ns. Before the drift, pe and pressure of the two models agree so I assume the dynamics of the explicit tip4p model I built is around correct before the COM drift.
I know that “fix momentum” does not guarantee prevention of COM drift and fix recenter is too strong to apply.
Then I cannot use “fix spring tether” because LAMMPS has noted that “One exception is for rigid bodies, which should not be used with the fix spring command, if the rigid body will cross a periodic boundary.” in its fix spring page. Since I am using pbc, certainly my water molecules move across boundary every now and then. And initial testing with fix spring tether mode did go wrong because thermodynamics properties e.g. pe are incorrect.
Hence I wonder is there any other way I can try to prevent COM drift all together with “fix rigid” ?
Apart from “fix rigid”, is there any inherent reason that “fix shake” + “explicit 4-site tip4p model” won’t work because in this case I might be able to use “fix spring tether” ?
In this case, I have been thinking about fix shake the virtual bond length between O site and massless M site, and the virtual angle HOM. However my initial test simulation broke down in early stage and I wonder if it is due to any intrinsic reason on LAMMPS setting that make this combination difficult (though I cannot find any from the fix shake LAMMPS page), or it might due to my poor parameter setting? For the latter I will do more tuning, while for the first reason I can give up before spending more tests on it.
Best Wishes
Lunna
Dear Axel/LAMMPS Users
"the cleanest way to keep a system from building up a COM drift is to
use a weak restraint via fix spring in tether mode. in some cases also
using fix momentum might do the trick and would be preferable to fix
recenter. the latter is mostly meaningful if you would have to
perfectly recenter the trajectory for analysis in postprocessing, but
this is very rarely needed"
Thanks so much for the great help on this thread!
Therefore I have been using "fix momentum" together with my explicit 4-site
tip4p water model with lj/cut/coul/long and regular pppm with "fix rigid
molecule"
"fix momentum" helps a lot.
However COM drift still happen within ~4ns for my explicit 4-site tip4p
water models, while no COM drift occurs for the default tip4p water models
in LAMMPS within ~18ns. Before the drift, pe and pressure of the two models
agree so I assume the dynamics of the explicit tip4p model I built is around
correct before the COM drift.
I know that "fix momentum" does not guarantee prevention of COM drift and
fix recenter is too strong to apply.
Then I cannot use "fix spring tether" because LAMMPS has noted that "One
exception is for rigid bodies, which should not be used with the fix spring
command, if the rigid body will cross a periodic boundary." in its fix
spring page. Since I am using pbc, certainly my water molecules move across
boundary every now and then. And initial testing with fix spring tether mode
did go wrong because thermodynamics properties e.g. pe are incorrect.
Hence I wonder is there any other way I can try to prevent COM drift all
together with "fix rigid" ?
you are asking the wrong question. you should find out *why* you get a
COM drift in the first place. there has to be something in your setup
that encourages it and that you have to eliminate or reduce. tinkering
with the symptoms rarely solves the real problem.
axel.
Dear Axel/LAMMPS Users
"you are asking the wrong question. you should find out why you get a
COM drift in the first place. there has to be something in your setup
that encourages it and that you have to eliminate or reduce. tinkering
with the symptoms rarely solves the real problem."
Regarding the reason of my COM drift, I have noticed that the COM has changed dramatically (might be too large to call it a drift) when the volume expanded to very large/infinity during the “fix rigid/npt” (298K and 1 bar).
In such case the pressure is still oscillating around 1 bar with a different oscillation pattern from before the volume change, and the pressure becomes almost constant at 1 (which I assume is reasonable for very large system volume it stuck itself with)
In real units, I have tried a lot of combinations of Tdamp (10-100) and Pdamp (100-1000) for “fix rigid/npt” using time step 0.5 fs/1.0 fs. The same problem of volume expanding happened at different times but none of the combinations was completely free from it.
The commands from my input file that could have affected the dynamics are listed below:
bond_style none
angle_style none
velocity all create 298 28362 dist gaussian mom yes rot yes
fix 1 all rigid/npt molecule temp 298 298 100 iso 1 1 1000
fix 2 all momentum 1 linear 1 1 1 angular
I have started from an unphysical initial configuration (with the bond/angles at the desired values but molecules are not in the right orientations) but since the system has not complained and has produced correct energies and distributions functions for the first ~2ns , so I assume that the “volume expansion problem” might have nothing to do with my initial configuration.
Many thanks and I am really sorry for my never-ending question list…
Best Wishe
Lunna
Dear Axel/LAMMPS Users
"you are asking the wrong question. you should find out *why* you get a
COM drift in the first place. there has to be something in your setup
that encourages it and that you have to eliminate or reduce. tinkering
with the symptoms rarely solves the real problem."
Regarding the reason of my COM drift, I have noticed that the COM has
changed dramatically (might be too large to call it a drift) when the volume
expanded to very large/infinity during the "fix rigid/npt" (298K and 1 bar).
In such case the pressure is still oscillating around 1 bar with a different
oscillation pattern from before the volume change, and the pressure becomes
almost constant at 1 (which I assume is reasonable for very large system
volume it stuck itself with)
In real units, I have tried a lot of combinations of Tdamp (10-100) and
Pdamp (100-1000) for "fix rigid/npt" using time step 0.5 fs/1.0 fs. The same
problem of volume expanding happened at different times but none of the
combinations was completely free from it.
this sounds like an issue with your simulation protocol. you should
not start with fix rigid/npt right away. usually know the desired
target density or at least some approximation initially and then
should start equilibration with a fixed volume integration. this will
prevent excessive times required to recompress the system to get back
to the desired volume.
axel.