1) Don't use the Berendsen thermostat.
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2699284/
https://sourceforge.net/p/lammps/mailman/message/34356144/
https://lammps.sandia.gov/threads/msg23729.html
2) There was a discussion recently about how to mix rigid molecules
together with non-rigid molecules in LAMMPS under NPT conditions.
That discussion sounds relevant to what you are trying to do:
https://lammps.sandia.gov/threads/msg77508.html
That discussion led to a LAMMPS example which mixes TIP5P with
ethanol. TIP5P is a popular model for water molecules which is
typically held rigid. (In the LAMMPS example below, the atoms in
these molecules were assigned to a group named "gRigid".) The
remaining molecules (in this case ethanol) are flexible. (The atoms
in these molecules were assigned to a group named "gFlexible".) I
tried several different protocols, applying different integrators to
each group. Most of the protocols I tried resulted in equilibrated
structures that were too big or too small, or too hot or too cold, or
failed to equilibrate.)
The protocol which worked the best is explained in this file
(scroll down to the "fxRigid", "fxFlexible", and "fxLanvegin" lines:
https://github.com/jewettaij/moltemplate/blob/master/examples/all_atom/force_field_OPLSAA/waterTIP5P%2Bethanol/run.in.npt
If you want the complete example (which includes all of the LAMMPS
input files and data files), go here:
https://github.com/jewettaij/moltemplate/tree/master/examples/all_atom/force_field_OPLSAA/waterTIP5P%2Bethanol
3) Just pointing something out:
Rather than holding molecules rigid, users typically build their
system using a force-field containing angle, dihedral, and improper
forces which allows the molecules to bend in a realistic way
(hopefully). That might be easier to get working than using an exotic
combination of rigid and flexible integrators.
There are several molecule builders on the LAMMPS prepost page
(https://lammps.sandia.gov/prepost.html) (Disclaimer: I wrote one of
them.) For what it's worth, there are a couple examples with benzene
here:
https://github.com/jewettaij/moltemplate/tree/master/examples/all_atom/force_field_OPLSAA/ethylene%2Bbenzene
https://github.com/jewettaij/moltemplate/tree/master/examples/all_atom/force_field_AMBER/ethylene%2Bbenzene
Cheers
Andrew
Dear all,
A word of caution about the rigid TIP5P - flexible ethanol system mentioned below. While relatively short runs went OK, longer ones inevitable ended up in crashes even when I reduced the time step to 1 fs from the original 2 fs. I guess there is something in mixing two thermostats when water and ethanol molecules interact so strongly they form hydrogen bonds.
Regards
Tatiana Kuznetsova
Dear all,
A word of caution about the rigid TIP5P - flexible ethanol system
mentioned below. While relatively short runs went OK, longer ones
inevitable ended up in crashes even when I reduced the time step to 1 fs
what kind of crash?
from the original 2 fs. I guess there is something in mixing two
unlike with fix shake (where 2fs is usually ok with hydrogen
containing systems), using fix rigid on water often requires a much
shorter time step to integrate the equations of motion in a stable
fashion, which is due to the rotational degrees of freedom. even 1fs
seems large to me. i remember having to use 0.5fs or 0.25fs depending
on thermodynamic conditions and choice of potentials.
thermostats when water and ethanol molecules interact so strongly they
form hydrogen bonds.
i don't think it is the thermostats causing problems, but rather the
fact that the integration of the equations of motion for rotation and
translation require different resolution. also, there are sometimes
issues with TIP5P with certain force fields, as (partially) charged
particles can come very close to the "satellite charge sites" or the
hydrogens and thus cause a "coulomb catastrophe" with large time
steps. CHARMM for example has this issue and thus uses a modified
TIP3P potential with LJ interactions added on the water hydrogen
(which are not present in the original TIP3P).
axel.
The crash was preceded by a spike in temperature, followed by the out of range atoms
590500 304.40907 -1.4897406 0.46266251 -1.0270781 -73.866713 20943.356 -2.0939991 0.42788286 0.16836921 0.0080063807
590600 303.87182 -1.4879322 0.46184595 -1.0260862 -465.75036 20954.904 -2.1006582 0.42873019 0.17507641 0.008919429
590700 1448.5642 -1.2839186 2.2016306 0.91771202 17051.72 20963.335 -1.9184155 0.44819917 0.17724456 0.0090531777
590800 946.76699 -0.74159264 1.4389636 0.697371 1935.4828 21093.792 -1.964391 0.99188896 0.21886863 0.012040794
590900 673.22368 -0.1213171 1.0232131 0.90189601 -4184.8012 21503.296 -1.8613126 1.5016576 0.22496153 0.013376439
ERROR on proc 7: Out of range atoms - cannot compute PPPM (../pppm_gpu.cpp:223)
Last command: run 2000000
ERROR on proc 5: Out of range atoms - cannot compute PPPM (../pppm_gpu.cpp:223)
Last command: run 2000000
The crash was preceded by a spike in temperature, followed by the out of
range atoms
which is exactly what i would expect with the two scenarios (or a
combination of the two) that i described.
axel.
Thanks for reporting the problem Tatiana, and for your feedback Axel.
Indeed I also thought it was weird that it seemed to run well at 2fs.
It looks like I didn't run the tests long enough. (These 10000
timestep tests were laughably too short to detect these kinds of
problems.)
I confess I don't have time to debug this example carefully this week,
but I will reduce the timestep used in the example from 2fs to 0.25fs
and add a warning that even this timestep may be too large. (While
testing the example today, I also ran into a moltemplate bug. Sigh.
Fixing that too.)
If anyone else has tried simulating TIP5P mixed with other molecules
and has feedback, please let me know.
Cheers
Andrew
P.S. A coulomb catastrophe is definitely possible with TIP5P. The
two negatively charged ("lone-pair") particles in the TIP5P model
provide no steric Lennard-Jones (1/r^12) repulsion with other atoms.
If there is a more recent version of TIP5P which places steric
particles at the lone-pair positions, feel free to send me an email
where to find it and I will include this model with moltemplate.