Hello all,
I would like to use tip4p/ice water potential which is a modified version of tip4p. I see that lj/cut/coul/long/tip4p pair style will be included as a patch in the near future. Could you give an approximate time frame of when it will happen?
Also, does anyone know how one can alternatively implement tip4p potential using standard pair styles (i.e. lj/cut/coul/long)? Assigning LJ parameters for oxygen and hydrogens, and placing the positive charges on the hydrogens and the negative on the 4th tip4p site did not work for me. The system blew up immediately, probably due to the 4th site being massless.
TIP4P should be released in the next couple weeks. I don't think
you can do TIP4P via standard pair styles. The 4th fictitious site
has to be treated differently.
I noticed that in the previous version of lammps (Mar06) the usage of
fix_temp/rescale required the input of the number of steps after which the
temperature was rescaled. In the file fix_temp_rescale.h line 36 was
int nevery;
However, this line was taken out in the same file for the latest version.
Even though "nevery" is not defined in the .h file, the .cpp doesn't give an
error of undefined variable when the code is compiled. Where is the variable
nevery stored now? The reason I'm asking is because I created my own fix
following the same structure of the previous fix_temp_rescale.h and .cpp and
I get an error when trying to run it on the new version.
The other aspect of this change you need to be aware of
for a fix you write is: if it is an end-of-step fix, then you
must define Nevery as an argument in your arg list for
the fix. It will be used by modify.cpp to call the fix. That's
why the variable is no longer stored by the fix itself.
To be consistent with all the other end-of-step fixes, you
should have Nevery be the 1st arg of the fix (after the style
of the fix).
Also, I'm setting a temperature increase in a region in space, but I was
wondering how I could set a flag so that when an atom crosses the boundary
of this region it is no longer affected by the heating in that region.
In general, you'll see this condition whenever the code iterates
through the list of atoms owned by the current processor. Most fixes
operate on groups of atoms (defined by the group command - defaulting to
group all) and since the group membership of an atom is defined using a
single 32 bit bitmask that code is a quick way of checking whether the
atom belongs to the group the fix is defined for.
There are two possible ways to implement the answer to your second
question. You can either simply check whether the coordinates of the atom
are outside the boundary (the slow easy way of doing it), or maintaining
per-atom data with your fix (a little more complicated - see fix_msd.cpp
or fix_spring_self.cpp for an example of maintaining per atom data). But
you'll still need to check if any atoms have reentered the region - so
the first option might be the best way.
I ran a simulation on a bead-spring model of polymer.
LJ/cut potential was used for pairs with e=1 and sigma=1. Bond potential
strength was set to 1000.
The result for E_pair, Tot Energy and E_mol came out fine. However, when
I tried to calculate the pair distance from the dump file, I found a few
pairs having really short distance, about 0.1. My question is, what is the
reason that such short pairs can appear? Are they reasonable? Why didn't they
blow up the energy (by estimation in this case, V_lj = 4*(1/0.1)^12 = 4E12 )?
I've attached input file test.bs, data.100 and output file log.lammps,
test.dump. I appreciate any help.
Thanks.
You may want to look at the special_bonds command. By
default it is set to 0,0,0 which means there is no pair
potential between 1-4 atoms (atoms 2 hops away in a bead-
spring chain). Since you only have bond and angle terms
defined, that could mean that 1-4 atoms get arbitrarily close
together.
You are also using angle harmonic with a default theta of
180, which may not be a good choice for that potential. You
might look at the cosine or cosine/squared potentials for
the 180 case.
Thanks, it seems to be my problem. I also appreciate your advice on angle
potential. It is very helpful to know that.
Now I added
special_bond 1.0 1.0 1.0
into my previous attached script file, but LAMMPS encountered SIGSEGV(11)
error after change. The print out messages are:
I would use special_bonds 0,0,1 for this system, not 1,1,1.
The reason being that you have defined bonds and angles,
so presumably you don't want additional pair terms between
those atoms (the 0,0), but you did not define dihedrals so
you do want a pairwise term between 1-4 atoms (the 1).
The latter should prevent 1-4 atoms getting too close,
which was the issue in your first email.