Fix rigid: Bad principal moments

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I want to simulate bulk CO2 with fix_gcmc command using EPM2 force field.

I encountered some problems when calculating the density profile:

I have tried Both rigid and shake,

  1. when using fix_shake, Lammps crashes with an Error: shake determinant =0.0. So I searched the solutions on Lammps Mailing List and found that “the SHAKE/RATTLE implementation in LAMMPS cannot process linear molecules with more than 2 atoms.”

2) therefore, I use fix_rigid instead. But I failed with the ERROR: Fix rigid: Bad principal moments (…/fix_rigid_small.cpp:2229).

I have looked through the LAMMPS Mailing List but I haven’t found the answers to the above problems.

The version of LAMMPS I am using is lammps-11Aug17.

Any help would be helpful.

Thank you very much

Best regards,

Juan Zhou

#####################INPUT SCRIPT############################

variable mu index -0.99968

variable disp index 0.5

variable press equal 100

variable pressatm equal ${press}*0.9869233

variable fugacoeff equal 0.49048

variable temp index 300.0

variable lbox index 50.0

variable spacing index 5.0

units real

atom_style full

boundary p p p

pair_style lj/cut/coul/long 14.0

pair_modify mix arithmetic tail yes

kspace_style ewald 0.0001

bond_style harmonic

angle_style harmonic

lattice sc ${spacing}

region box block 0 {lbox} 0 {lbox} 0 ${lbox} units box

create_box 2 box &

bond/types 1 &

angle/types 1 &

extra/bond/per/atom 2 &

extra/angle/per/atom 1 &

extra/special/per/atom 2

molecule CO2 CO2.txt

create_atoms 0 box mol CO2 464563 units box

pair_coeff 1 1 0.0558535 2.757

pair_coeff 2 2 0.1598562 3.033

bond_coeff 1 2057.09623 1.149

angle_coeff 1 111.950929 180

mass 1 12.0107

mass 2 15.9994

minimize 1.0e-4 1.0e-6 100 1000

min_style cg

group gas type 1 2

neighbor 2.0 bin

neigh_modify every 1 delay 10 check yes

velocity all create ${temp} 4928459 mom yes rot yes dist gaussian

timestep 1.0

thermo_style custom step temp press density atoms pe ke

thermo 100

fix 1 gas rigid/nvt/small molecule temp {temp} {temp} 100 mol CO2

fix_modify 1 dynamic/dof yes

variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans)

fix 2 gas gcmc 10 1000 0 0 54341 {temp} {mu} {disp} mol CO2 pressure {pressatm} fugacity_coeff ${fugacoeff} &

tfac_insert ${tfac} group gas rigid 1

compute_modify thermo_temp dynamic/dof yes

run 100000

##################### CO2.txt ############################

3 atoms

2 bonds

1 angles

Coords

1 2.0000 1.0000 0.0000

2 3.0221 1.5679 0.0000

3 0.9779 0.4321 0.0000

Types

1 1

2 2

3 2

Charges

1 0.6512

2 -0.3256

3 -0.3256

Bonds

1 1 1 2

2 1 1 3

Angles

1 1 2 1 3

Special Bond Counts

1 2 0 0

2 1 1 0

3 1 1 0

Special Bonds

1 2 3

2 1 3

3 1 2

#Shake Flags

#1 1

#2 2

#3 2

#Shake Atoms

#1 1 2 3

#2 1 2

#3 1 3

#Shake Bond Types

#1 1 1 1

#2 1

#3 1

2

Hi,

I want to simulate bulk CO2 with fix_gcmc command using EPM2 force field.

I encountered some problems when calculating the density profile:

I have tried Both rigid and shake,

1) when using fix_shake, Lammps crashes with an Error: shake
determinant =0.0. So I searched the solutions on Lammps Mailing List and
found that “the SHAKE/RATTLE implementation in LAMMPS cannot process linear
molecules with more than 2 atoms.”

yes. the solution to the linear system required by the SHAKE algorithm
is degenerate for linear molecules. so to do keep CO2 linear with with
SHAKE, you need to implement a special case, where you translate the
linear molecule to an "effective 2 atom" system, solve the constraint
for that and project the solution back to the linear molecule. that is
very complex to do in parallel, and - obviously - nobody wanted this
feature enough to contribute an implementation to LAMMPS.

2) therefore, I use fix_rigid instead. But I failed with the ERROR:
Fix rigid: Bad principal moments (../fix_rigid_small.cpp:2229).

I have looked through the LAMMPS Mailing List but I haven’t found the
answers to the above problems.

you have already quoted the answer to the first problem, and you'll
find the answer to what error messages mean here:
http://lammps.sandia.gov/doc/Section_errors.html

axel.

3

On Fri, Mar 30, 2018 at 11:17 AM, Juan Zhou

  1. therefore, I use fix_rigid instead. But I failed with the ERROR:
    Fix rigid: Bad principal moments (…/fix_rigid_small.cpp:2229).

you’ll
find the answer to what error messages mean here:
http://lammps.sandia.gov/doc/Section_errors.html

…as for what to do about it, (unless somebody has any better suggestions), I would add another ghost atom to your system which does not lie on the same line as the other 3 atoms. Put it somewhere reasonably near the center of your molecule, but not on the line.
Create a new atom type for this atom (and a new “group” for all the atoms which belong to this type). You can then use “neigh_modify exclude” to disable the computation of pairwise (nonbonded) forces between these atoms and all of the atoms in your system. (You can also do the same thing by setting the corresponding epsilon parameters of the pair_coeffs for this atom tyoe to 0, which you should also do.)

Using the “neigh_modify exclude” command reduces computation time because forces involving these atoms are never calculated. So adding the extra atom should hopefully not slow down the simulation very much.

If anybody else has a better idea, please post.

Cheers
Andrew

4

Hi,

Many thanks for your suggestion!

I have tried to add another type of atom as you suggested. But I encountered the same error ERROR: Fix rigid: Bad principal moments (…/fix_rigid_small.cpp:2236). (When I set the bond_coeff and angle_coeff of the ghost atom as a non-zero value, the density in log file remains zero.) Input file and CO2.txt attached.

It has taken me a lot of time to simulate bulk CO2. I would really appreciate it if somebody could help me solve this problem.

in.gcmc_CO2 (2.08 KB)

CO2.txt (535 Bytes)

5

On Fri, Mar 30, 2018 at 11:17 AM, Juan Zhou
> 2) therefore, I use fix_rigid instead. But I failed with the
> ERROR:
> Fix rigid: Bad principal moments (../fix_rigid_small.cpp:2229).

you'll
find the answer to what error messages mean here:
http://lammps.sandia.gov/doc/Section_errors.html

...as for what to do about it, (unless somebody has any better suggestions),
I would add another ghost atom to your system which does not lie on the same
line as the other 3 atoms. Put it somewhere reasonably near the center of
your molecule, but not on the line.

sorry, but this makes no sense. it is contrary to the root of the
problem. the error stems from the fact, that the molecule defined is
not *exactly* linear (only very, very nearly linear), thus fix
rigid/small cannot detect it as such and therefore will not treat it
as linear, which will result in what is effectively the same issue as
with fix shake. adding another atom does not solve either problem. fix
shake won't allow "shake clusters" where not all bonds are connected
to the same atom and for fix rigid, placing an atom away from the line
of connected atoms, will modify the principal moments to be even
further away from the exactly linear geometry. it will have an effect
only, when it is modifying the principal moments so much, that it will
result in completely *bogus* principal moments, and then the simulated
molecule would not move like a CO2 molecule would move.

[...]

If anybody else has a better idea, please post.

it is *trivial* to enter a CO2 geometry that is *exactly* linear. even
more so, LAMMPS provides an example for it and for how to use it with
fix gcmc and using a fix rigid variant.

axel.

6

Hi,

Many thanks for your suggestion!

I have tried to add another type of atom as you suggested. But I encountered
the same error ERROR: Fix rigid: Bad principal moments
(../fix_rigid_small.cpp:2236). (When I set the bond_coeff and angle_coeff of
the ghost atom as a non-zero value, the density in log file remains zero.)
Input file and CO2.txt attached.

It has taken me a lot of time to simulate bulk CO2. I would really
appreciate it if somebody could help me solve this problem.

you wasted a lot of time by not looking over the materials that LAMMPS provides.
there is an example included in the LAMMPS distribution for doing GCMC with CO2.
if you want to do bulk CO2, you should simply modify it in such a way,
that you first have a "full" bulk system of molecules using the create_atoms
command at the desired density, properly equilibrate that to the
desired temperature,
and then add fix gcmc to make it a GCMC simulation.

this should be rather obvious and straightforward to do.
the fact, that feel you have already spent too much time on it,
does not entitle you for any special treatment, or somebody doing your
work for you.
in fact, this is the normal condition of doing research.
for as long as you don't understand properly what you are doing,
things will always take a long time.
so if you cannot set up a complex simulation, you have to figure out
ways to understand.
e.g. by simplifying what you do and by more carefully studying
documentation and examples.

finally, for using fix gcmc, you are strongly recommended to upgrade
your LAMMPS version, as there have been multiple bugfixes applied to
fix gcmc since the version you are using.
that doesn't affect the specific, and rather trivial problem that you
are currently struggling with, but more subtle and less
straightforward to detect issues.

axel.

7

Thank you so much for your suggestions.

I will try it and study the documentation and examples more carefully.

Best wishes,
Juan