Gay-Berne molecule interacting hydrocarbones(CHn)

Dear users

In simulation of Liquid-crystal Polymer system, mesogen (rigid body, ellipsoidal molecule) must interact with other polymer molecules via coarse-graining.
(for instance, -CH2-(mesogen)-CH2- system can be modeled by hybrid pair-style of Gay-Berne and LJ)
However, mesogens freely rotate, in respect to other hydrocarbons, because there is no constraint on the mesogen orientation in polymer chains and only center of masses are bonded.

As far as I know, LAMMPS does not provides angle defined between orientation of ellipsoidal-atom_style and atomic-atom_style, while it is a key to solve the problem.

A few researches are found in regards to this matter, and they suggest constraint pseudo-angle which is defined between mesogen’s long axis vs. backbone of hydrocarbons.
(Mark R. Wilson’s paper:J. Chem. Phys.107(20), 22 November 1997 provides the detail)

I think solution of this kind, may expand LAMMPS’s applicability greatly.

If anyone ever had similar issues, or have any suggestions in this matter, any comments or suggestions are deeply appreciated.

Thanks.

Dear users

In simulation of Liquid-crystal Polymer system, mesogen (rigid body,
ellipsoidal molecule) must interact with other polymer molecules via
coarse-graining.
(for instance, -CH2-(mesogen)-CH2- system can be modeled by hybrid
pair-style of Gay-Berne and LJ)
However, mesogens freely rotate, in respect to other hydrocarbons, because
there is no constraint on the mesogen orientation in polymer chains and only
center of masses are bonded.

As far as I know, LAMMPS does not provides angle defined between orientation
of ellipsoidal-atom_style and atomic-atom_style, while it is a key to solve
the problem.
A few researches are found in regards to this matter, and they suggest
constraint pseudo-angle which is defined between mesogen's long axis vs.
backbone of hydrocarbons.
(Mark R. Wilson's paper:J. Chem. Phys.107(20), 22 November 1997 provides the
detail)

I think solution of this kind, may expand LAMMPS's applicability greatly.

well, there you have just defined a project for yourself to get fame
and fortune. just grab a text editor and implement it. LAMMPS has a
very liberal approach to include contributed code from everybody, for
as long as it doesn't need any significant or problematic changes to
its core functionality.

cheers,
    axel.

More speficically you could add an angle style
that induces forces and torques based on the
orientation of one (or more) of the ellipsoids
in the 3-body interaction.

Steve

Thanks Steve and Axel for all the comments.

Now, I’m currently working on this matter right now,
but I’m facing some issues that you may help.

First of all, is angle_style can still be utilized when only two atoms (Asphere and sphere) interact?
basically, it seems plausible since there are two vectors (one is between center of masses, and another is based on quaternion of asphere atom).

Secondly, if the answer is yes, how can ev_tally() in angle_style can be used?

I’m bit afraid if the answer on the first question is no, and there are too much point to edit.

Thanks again, and any types of advices are deeply appreciated.

Hayoung.

2013/6/19 Steve Plimpton <[email protected]>

Thanks Steve and Axel for all the comments.

Now, I'm currently working on this matter right now,
but I'm facing some issues that you may help.

First of all, is angle_style can still be utilized when only two atoms
(Asphere and sphere) interact?

an angle style will add forces. LAMMPS doesn't care how. so if you
only add forces to two of the three particles that make up an angle,
that is just fine.

basically, it seems plausible since there are two vectors (one is between
center of masses, and another is based on quaternion of asphere atom).

i don't see the reasoning of that. an angle is always defined between
three particles; particularly the center of masses of the three. if
you have only two particles, this would be a bond. of course, you can
do a bond style interaction that goes beyond just updating the
particle forces (between the center of masses) just as well as an
angle style interaction.

Secondly, if the answer is yes, how can ev_tally() in angle_style can be
used?

you likely need to write your own variable of ev_tally(). have a look
at pair.cpp. it has a lot of them to cater to the many different ways
how LAMMPS handles pairwise nonbonded interactions (and how it treats
implicit bonds/angles/dihedrals in manybody potentials).

I'm bit afraid if the answer on the first question is no, and there are too
much point to edit.

there ain't no escape from the blues. :wink:

axel.

To Axel, Steve, and users of LAMMPS,

Now currently working on this humble project, and here are some difficulties I’m facing now.I’m going to give brief explanation, and if anyone find some holes in it please give me as commentary.

% -------------
my problem here is to fix the orientation of the asphere molecule, with respect to neighboring (actually, bonded) sphere atoms.
(because in current LAMMPS results provide unrealistic spinning of the asphere molecule because there is no way to constrain this rotation.)
In other words, the angle (generated between orientation of asphere molecule and the vector connecting two center of masses) is fixed.
(I’d like to call it Pseudoangle, following Mark R. Wilson’s paper)

Since there is only two molecules in this style, I decided to use bond_style.
one vector between two atoms, and quaternion fixed on the asphere is called to compute the pseudoangle in between them,
and fictitious (as you may call it) energy, force, and torque for these two atoms are easily assessed in compute() function.

but when it comes to single() and equilibrium_distance() function, which is requisite for this type of bond_style,
there is no way to implement this pseudoangle style, because the instances for this functions are quite restricted.

So here are the questions;

  1. there is no way to define equilibrium length of bonding, because length of the bond does not determine the energies in this style. So is it possible to left it returns 0, as long as fix_SHAKE unused?
  2. What is the single() function is for? energy and forces are computed via compute() function, so I cannot figure out when this function is called.
    Also, can you give me a reasoning on why the format (only ids of atoms, and square of bond) is so restricted unlike that of compute()?
    % --------------

if the script I’m working on may help you readers understand, I’m willing to provide it in detail.

Thanks for reading lengthy questions. Any comments would be appreciated.

Hayoung.

2013/7/1 Axel Kohlmeyer <[email protected]…1125…>

To Axel, Steve, and users of LAMMPS,

Now currently working on this humble project, and here are some difficulties
I'm facing now.
I'm going to give brief explanation, and if anyone find some holes in it
please give me as commentary.

% -------------
my problem here is to fix the orientation of the asphere molecule, with
respect to neighboring (actually, bonded) sphere atoms.
(because in current LAMMPS results provide unrealistic spinning of the
asphere molecule because there is no way to constrain this rotation.)
In other words, the angle (generated between orientation of asphere molecule
and the vector connecting two center of masses) is fixed.
(I'd like to call it Pseudoangle, following Mark R. Wilson's paper)

Since there is only two molecules in this style, I decided to use
bond_style.
one vector between two atoms, and quaternion fixed on the asphere is called
to compute the pseudoangle in between them,
and fictitious (as you may call it) energy, force, and torque for these two
atoms are easily assessed in compute() function.

but when it comes to single() and equilibrium_distance() function, which is
requisite for this type of bond_style,
there is no way to implement this pseudoangle style, because the instances
for this functions are quite restricted.

So here are the questions;
1. there is no way to define equilibrium length of bonding, because length
of the bond does not determine the energies in this style. So is it possible
to left it returns 0, as long as fix_SHAKE unused?

yes. the corresponding function is currently used in exactly two places:
1) fix shake and
2) for tip4p pair and pppm styles.

2. What is the single() function is for? energy and forces are computed via
compute() function, so I cannot figure out when this function is called.
Also, can you give me a reasoning on why the format (only ids of atoms, and
square of bond) is so restricted unlike that of compute()?

the single functions (they also exist for pair and others) exist to
compute the forces and energies for just a single bond (or pair). this
can be used in certain fixes or computes. the format is based on the
needs (which are rather specific).

for the time being, it may be sufficient to have it return zeroes just as well.

axel.

You can see where single() is called by grepping for
“pair->single”.

If you don’t want your pair style to define it (some don’t),
then don’t have a single() method in your pair class,
and set single_enable = 0 in the constructor.

Steve

You can see where single() is called by grepping for
"pair->single".

the question was about Bond::single() and you cannot disable that via
single_enable = 0. it is a pure function.

  virtual double single(int, double, int, int, double &) = 0;

axel.