I am trying to do NVT simulation of elliptical particles with off-centre point dipoles attached to it. I am attaching the data file and the LAMMPS input file. The temperature I am trying to fix is 2.0 but I am not able to do so.
I used atom_style hybrid. Used 2 types of particles (ellipsoid without dipole and point particle with point dipole. ) and used fixed rigid to treat them as single particles. Need to know if I am doing it correctly or not.
TIA.
This makes sense. I don't think the code for "fix rigid" understands
non-point-like particles (ellipses, dipoles). I don't know of any
other way to accomplish what you are trying to do using the official
LAMMPS code.
The most advanced simulations using ellipsoids in LAMMPS are
implemented using the "MOLC" code developed by Matteo Ricci and
Ottello Roscioni. It does something very similar to what you are
trying to do. The paper they wrote includes their LAMMPS code, but
their code lacks documentation. The paper is here.
They have not abandoned this. If you email Otello, he would probably
be willing to help you use their code because they want to promote
this method.
If you read their paper, it will probably be confusing because it
discusses a complex coarse-graining procedure for turning all-atom
molecules into coarse-grained ellipsoidal molecules. Fortunately
their code is general and works for all ellipsoidal simulations. You
can use their LAMMPS code without following their complicated
coarse-graining procedure. Several pair_styles (and also bond_styles)
have been created in order to make those kinds of simulations
possible. But the code has not been submitted for inclusion with the
LAMMPS distribution yet. (I think Matteo got bogged down writing the
documentation and it never happened.)
The name of the source code files are:
bond_ellipsoid.cpp
bond_ellipsoid.h
bond_molc.cpp
bond_molc.h
pair_coul_cut_offcentre.cpp
pair_coul_cut_offcentre.h
pair_coul_long_offcentre.cpp
pair_coul_long_offcentre.h
pppm_offcentre.cpp
pppm_offcentre.h
... and the name of the
bond_style ellipsoid
bond_style molc
pair_style coul/cut/offcentre
pair_style coul/long/offcentre
kspace_style pppm/offcentre
You will need to download the LAMMPS source code, copy these files
into the "src/" subdirectory, and recompile LAMMPS.
In your simulations, you would probably want to use something like:
pair_style hybrid/overlay gayberne 1.0 1.0 -3.0 14.
coul/long/offcentre 14. 4 &
kspace_style pppm/offcentre 1e-4 4 &
Unfortunately I suspect their code will probably not compile because
there have been some changes to LAMMPS since their paper came out.
(The syntax of the "numeric()" and "inumeric()" functions has
changed.) If you cannot figure out how to get LAMMPS to compile, then
there is an old version of LAMMPS that includes their code and
compiles. It is located here.
You will need to download the second archive file, unpack the LAMMPS
subdirectory and compile that version of LAMMPS.
Alternatively you can try emailing Otello for the latest source code.
There is also an example at the link above which demonstrates how to
run these kinds of ellipsoidal simulations. Unfortunately it is not a
simple example. It requires installing "moltemplate" to create the
LAMMPS data files and input scripts you need to run the example. You
must download the first archive file from the link above and unpack
it. The example is located in the directory named
"supplemental_figure_6_ellipsoidal_biphenyl".
Incidentally, that example is part of a paper
(https://doi.org/10.1016/j.jmb.2021.166841). I can send it to you if
you don't have easy download access.
--- VISUALIZE your results ---
No matter how you implement this simulation (using fix rigid or using
the new "MOLC" code from Matteo and Otello), you must visualize your
trajectory. This is the first thing you should try when things are
not working.
I ran a short simulation using a modified version of your in.ellipse
file (attached) and I visualized the trajectory file with OVITO.
Although the distance between the two particles remains the same, the
orientations of both the ellipsoid and dipole are free to spin. That
is presumably not what you wanted, and it would easily explain why
energy is not conserved and your temperature is increasing. (The
point dipoles were probably getting too close to each other. I am
surprised the simulation did not crash.)
To visualize your trajectory files, I think you need to create a
single trajectory file containing both the ellipsoid and dipole
orientations. (See attached "in.ellipsoid" file.) Then display both
the dipole and ellipsoid degrees of freedom. Then select two of the
particles which are part of the same molecule. Then add a
modification "Invert Selection" Then add modification "delete
selected particles".
that statement is not correct, which can be easily confirmed by looking at the source code.
you will see that for rigid bodies both per particle force and per particle torque are considered during the initial_integrate() and final_integrate() functions
and that ellipsoids and spheriods (but not point dipoles) have extra degrees of freedom that are removed from the system in the dof() function.
This makes sense. I don't think the code for "fix rigid" understands
non-point-like particles (ellipses, dipoles).
that statement is not correct, which can be easily confirmed by looking at the source code.
you will see that for rigid bodies both per particle force and per particle torque are considered during the initial_integrate() and final_integrate() functions
and that ellipsoids and spheriods (but not point dipoles) have extra degrees of freedom that are removed from the system in the dof() function.
axel.
That's good to know. Soumik is using this command to integrate the
equations of motion:
fix fixrigid all rigid/nvt/small molecule temp 2.0 2.0 0.1
I did not check the source code. But, for what it's worth, when I ran
his simulation and attempted to look at the trajectory with OVITO, it
seemed as though the dipole and ellipsoid orientations were not
frozen. They seem to be spinning (more-or-less) freely. But the
distance between the dipole and the ellipsoid was preserved. (But I
could definitely be making a mistake when using OVITO.)
I confess, part of the reason I replied was to promote the MOLC code.
(It's really cool.) But if there's a way to run this simulation with
fix/rigid instead, then that might be better. The MOLC pair_styles
and bond_styles are hard to use without documentation.
Andrew
P.S. If anyone is curious, I am using my customized, updated version
of his "in.ellipsoid" file, where I changed the dump command to
include the dipole information to make visualization easier. I
attached that input script to this message. The remainder of the
"in.ellipsoids" file is the same as what he posted. I'm using these
instructions to display the ellipsoids and dipole directions in OVITO:
(You have to alter the default size of the particles and reduce the
dipole lengths to see anything. It also helps if you only display one
of the molecules instead of the entire system.)
Hi - I’m coming to this thread late, but I think the current rix rigid should work properly with bodies consisting of a mix of point, sphere, or ellipsoid particles,
any of which could have a point dipole assigned to them. And the point dipole should rotate with the body. If that isn’t working is is likely a bug.
Andrew, the simplest test would be a two particle rigid dimer (point particles) with a dipole assigned to one or both of the point particles.
Can you run that test with OVITO or dump out the body and dipole orientations to see if the bodies remain rigid?
Just for a simple 2-body system.
Hi - I’m coming to this thread late, but I think the current rix rigid should work properly with bodies consisting of a mix of point, sphere, or ellipsoid particles,
any of which could have a point dipole assigned to them. And the point dipole should rotate with the body. If that isn’t working is is likely a bug.
Andrew, the simplest test would be a two particle rigid dimer (point particles) with a dipole assigned to one or both of the point particles.
Can you run that test with OVITO or dump out the body and dipole orientations to see if the bodies remain rigid?
Just for a simple 2-body system.
Hi Steve, Axel, and Soumik
I think I figured out the source of Soumik’s problem with temperature, and it’s not what I thought. The molecules are rigid. The ellipsoid and dipole do indeed maintain the same orientation and position as you would expect for a rigid molecule. (This is true when I reduce the size of the system to a single rigid molecule, or when I simulate the original full-size system.)
My apologies for the false alarm. (When I attempted to view the trajectory in OVITO, it looks like I guessed the meanings of the quaternion columns incorrectly. It’s easy to do. Incidentally, the corrected OVITO settings are included in a PNG image which is part of an attachment included with this post.)
So, why is the temperature weird? to the documentation, for aspherical particles, one needs to use “fix nvt/asphere” for correct thermostating of ellipsoidal particles. But I can’t think of a way to combine that integrator with fix rigid. Perhaps this is the source of the problem.
I tried simulating Soumik’s system using several different thermostats including fix langevin, as well as Soumik’s original choice (“fxrigid ED rigid/nvt/small molecule temp 2.0 2.0 0.1”). When I used his original thermostat (“fxrigid ED rigid/nvt/small molecule temp 2.0 2.0 0.1”), it takes a while for the system to equilibrate. The temperature is reported to be significantly higher than the target temperature (approximately 2.7 instead of 2.0). But I might not have waited long enough for equilbration (only 2000 timesteps).
When in doubt, use the simplest possible thermostat: fix langevin. When I use fix langevin (together with “fxrigid all rigid/nve/small molecule”), the reported temperature was also too high. (Again, approximately 2.8 instead of 2.0) (I get similar behavior when I use “fix rigid/small molecule langevin 2.0 2.0 0.5 123456”)
However, unless I am mistaken, using “fix langevin” should eventually produce a system with the correct kinetic energy once thermal equilibrium is reached (even if the reported temperature is higher than it should be, and the rotational motion is underdamped). So that’s probably what I would do if I were in Soumik’s shoes. (This is called “Langevin Method #1” in the attached “run.in” file.) I’m curious if someone has a better idea.
Andrew
P.S. If anyone cares, the correct OVITO settings are included at the end of this email in the attached .gz file. The simplified “run.in” LAMMPS input script contains the integrator and thermostat settings. You can try different thermostats and see what happens. (You can change the number of molecules in the system by editing the “system.lt” file and following the instructions in the comments and README files.)
So, why is the temperature weird? to the documentation, for aspherical particles, one needs to use “fix nvt/asphere” for correct thermostating of ellipsoidal particles. But I can’t think of a way to combine that integrator with fix rigid. Perhaps this is the source of the problem.
For individual particles fix nvt/asphere is good. But for rigid bodies fix rigid/nvt/small should do the equivalent thing. And there is a langevin option to use that as a thermostat instead of NoseHoover. This will thermostat both the translational and rotational DOF of the collection of bodies.
You also have to be careful about how you measure the temp of the system. You should use compute temp/body.
the temperature reported in the thermo output is by default from “compute temp”, but that compute doesn’t know about extended particles.
so you have to do:
compute ext all compute temp/asphere
thermo_modify temp ext
to get the correct temperature reported in thermo output.
When using any rigid fixes, they will automatically compute their own internal temperatures for translation and rotation and also count and report the removed degrees of freedom due to extended particles. When using “external” thermostats like fix langevin they will only consider point particles and use an instance of compute temp internally to get the current temperature of the group of particles they are operating on.
the temperature reported in the thermo output is by default from “compute temp”, but that compute doesn’t know about extended particles.
so you have to do:
compute ext all compute temp/asphere
thermo_modify temp ext
to get the correct temperature reported in thermo output.
Thanks Axel. For what it’s worth, I tried playing with using “compute temp/asphere” and “compute_modify extra/dof”, but I still end up with temperatures that are reported to be (30%-50%) too high, which is roughly similar to what Soumik reported.
(Steve’s suggestion to use “compute temp/body” doesn’t work because “Compute temp/body requires atom style body”. Regarding your suggestion, I can’t use group “all” with “compute temp/asphere” because only half of the particles are aspherical. Instead I tried using the “spheroid” group for this which has only aspherical particles. But this means I’m currently calculating the kinetic energy of only the spheroids without considering the kinetic energy from the point dipoles, which could be contributing to the problem.)
I’m giving up. I’ll let Soumik follow up if he needs more help.
In my opinion, this is an edge-case issue that does not affect many users. But I strongly suspect that the system is behaving correctly (equilibrating to the correct temperature) even if the temperature is reported incorrectly.
An updated version of the “run.in” input script is attached which includes various settings I tried. (It is part of a .gz archive which is attached. The relevant lines from “run.in” are 34-38,68) (If I had more time to play with this, I would compute the temperature from the kinetic energy to verify this is true, and also read the docs and source code to dig deeper into what is going on.)
Andrew
When using any rigid fixes, they will automatically compute their own internal temperatures for translation and rotation and also count and report the removed degrees of freedom due to extended particles. When using “external” thermostats like fix langevin they will only consider point particles and use an instance of compute temp internally to get the current temperature of the group of particles they are operating on.
P.S.
Thanks. That makes sense.
I realize this is besides the point you were making here, but or the record, it doesn’t seem to matter if I use:
“fix rigid/nvt/small”
“fix rigid/small molecule langevin”
“fix langevin” + “fix rigid/small/nve”
…In all cases, the reported temperature is too high (between 30% and 50%).
(See lines 52, 63-64, 68 of “run.in” in attached .gz file.) If the temperature is just being reported incorrectly, then it makes sense that it does not matter how I integrate the equations of motion. Incidentally, the attached .gz file contains a picture of what each molecule looks like.
