Using Bond, Angle and Dihedral Information in Ellipsoids

Hi,

I am performing coarse-grained MD simulations on polymer chains, with Gay-Berne potential for non-bonded interactions in LAMMPS (30 Jul 2016).

While running nvt/asphere or npt/asphere, I am setting a temperature of 300K, but the temperature of my system, is hovering around double that value, 600K. Could this be a problem of considering the mass of the atoms twice, once as a particle and once as an ellipsoid?

Here are my implementation details:-

To provide bond, angle and dihedral information in the ellipsoidal beads, I am using atom_style hybrid ellipsoid molecular.

The pair_style is the following: pair_style gayberne 1 1 2 25

In the data file, I am defining the following:
1000 atoms
1000 ellipsoids

990 bonds
980 angles
970 dihedrals

1 atom types
1 bond types
1 angle types
1 dihedral types

Masses
1 288.30569
PairIJ Coeffs
1 1 0.6 1.74 9.5 9.5 0.3 9.5 9.5 0.3

In the atoms section, I am setting the Ellipsoid flag of all atoms to 1 as there is only one type of beads in the system, and all beads are anisotropic. The density of the ellipsoid is set as the value obtained by dividing the mass by volume of the ellipsoid.
Sample line of my atoms section:
Atoms
Atom_type Atom_id x y z Ellipsoidflag Density MoleculeID
1 1 28.11866053 54.654925277 22.118563583 1 14.62 1

Sample line of my Ellipsoids section:
Ellipsoids

1 12.650 1.725 1.725 0.357795920 0 0.006946311 0.933773971

Please suggest if I am making a fundamental mistake in implementation of ellipsoidal beads in polymer chains for using the Gay-Berne potential.

Thanks,
Sandipan

Hi,

I am performing coarse-grained MD simulations on polymer chains, with Gay-Berne potential for non-bonded interactions in LAMMPS (30 Jul 2016).

While running nvt/asphere or npt/asphere, I am setting a temperature of 300K, but the temperature of my system, is hovering around double that value, 600K. Could this be a problem of considering the mass of the atoms twice, once as a particle and once as an ellipsoid?

no.

how do you compute the temperature? the default temperature compute will consider all particles as point particles.

axel.

I am getting the temperature from the thermo output:
thermo_style custom step press temp density etotal pe ke epair

I had ran the following fix:
fix 3 all npt/asphere temp 300 300 0.5 iso 1 1 5

The timestep is 0.005 fs.

Note: I have performed all-atom and coarse-grained simulations on the same polymer with standard LJ potential (using atom_style full in that case), and got results which match with experimental properties. This problem of getting double the desired temperature is occurring only while running simulations with ellipsoid beads.

that is what i was referring to. you are not computing the correct temperature since your system has more degrees of freedom than a system with point particles. the thermo output cannot know about this and will have to be set up to use a different compute with thermo_modify or you need to ignore its temp output and look at the output of compute temp/asphere instead.

LAMMPS is very flexible and supports many kinds of systems and types of interactions, but if you deviate from what is common and what is assumed for the default settings, you need to instruct LAMMPS accordingly or else you will get bogus results. in may cases the choices are ambiguous and thus LAMMPS makes no choice at all and it is expected of the user to make the suitable choices and modifications to the input instead. that is the price you have to pay for the flexibility.

axel.

Thanks Axel.

I will try thermo_modify or compute temp/asphere, whichever will suit my purpose better and post the results here.

It is perfectly understandable that the user needs to give specific instructions to lammps about the system in order to get proper results.

Regards,
Sandipan

The problem of getting double the desired temperature of the system has been solved with the following:

fix 3 all npt/asphere temp 300 300 0.5 iso 1 1 5
compute 1 all temp/asphere
compute_modify 1 extra 1 (this step is to reduce the number of dof from 6 to 5 because the lengths of the 2 minor axes are the same i.e. the ellipsoid is symmetric about the major axis)
thermo_modify temp 1

Thanks Axel.