Fix rigid: Bad Principal moments

Hello,
I am simulating CO2 molecules with fix nvt/rigid/small molecule. For the initial case without any static framework, I am able to simulate it correctly with the defined fix.

The problem occurs when I define a fixed framework. After defining a static framework I get the error “Fix rigid: Bad Principal moments”. With the small number of random molecules (less than 10), this command is working fine with the framework. The problem occurs when the molecules exceed 10. The dump file after minimization suggests that the atoms are minimized appropriately. The framework is frozen with setforces to zero.

There are no overlapping atoms as well.

I tried to attach the file with the thread but couldn’t do so. Following is the datafile script:

dimension         3
units             real
boundary          p p p
atom_style        full

pair_style      hybrid/overlay lj/cut/coul/long 12.0 lennard/mdf 11.0 12.0

kspace_style    pppm 1.0e-5 #ewald
bond_style      harmonic
angle_style     harmonic
dihedral_style  none
improper_style  none

# ...

read_data      Framework.data

region         box block 0 24 31.0 50.0 0 105
molecule       CO2 CO2_EPM2_upd.mol
create_atoms   0 random 300 885631 box mol CO2 154798

set type          1                 charge  0.999
set type          2                 charge  1.668
set type          3                 charge -0.889
set type          4                 charge  0.6512
set type          5                 charge -0.3256

group   framework  type     1   2   3
group   co2    type     4   5

# ======
# Force-fields
# ========

pair_coeff	1	1	lj/cut/coul/long	2.6504000	2.6185000		
pair_coeff	1	2	lj/cut/coul/long	0.2053280	3.0969000		
pair_coeff	1	3	lj/cut/coul/long	6.2702000	1.1966000		
pair_coeff	1	4	lj/cut/coul/long	0.0702167	3.2901500		
pair_coeff	1	5	lj/cut/coul/long	0.1187568	3.4281500		
pair_coeff	2	2	lj/cut/coul/long	0.4780000	2.3705000		
pair_coeff	2	3	lennard/mdf	        0.0000009	0.0000000		
pair_coeff	2	4	lj/cut/coul/long	0.1634631	2.5637500		
pair_coeff	2	5	lj/cut/coul/long	0.2764637	2.7017500		
pair_coeff	3	3	lj/cut/coul/long	0.0002730	4.7441500		
pair_coeff	3	4	lj/cut/coul/long	0.0881482	2.9239500		
pair_coeff	3	5	lj/cut/coul/long	0.1490842	3.0619500		
pair_coeff	4	4	lj/cut/coul/long	0.0559000	2.7570000		
pair_coeff	4	5	lj/cut/coul/long	0.0945432	2.8950000		
pair_coeff	5	5	lj/cut/coul/long	0.1599000	3.0330000		


pair_modify 	tail yes mix arithmetic

# but meant to be rigid?
bond_coeff      1 11500 1.149
angle_coeff     1 11500 180

# ===================

# Simulation Protocols

dump         molecule_structure_init		all custom 10 traj-nvt-init.lammpstrj id mol type x y z ix iy iz 

velocity          framework    set     0.0 0.0 0.0
fix		  freeze    framework    setforce  0.0 0.0 0.0


minimize          1.0e-6    1.0e-7   1000   10000

velocity          all       create   ${Tinit}  5871232    rot no     dist gaussian
fix               integrate_nvt  co2  rigid/nvt/small molecule temp ${T} ${T} 50
fix               rmMom      co2  momentum   1  linear  1  1  1

undump molecule_structure_init
# ==================================================================================================================

timestep     	  ${dt}
reset_timestep    0

dump         molecule_structure		all custom 1000 traj-nvt.lammpstrj id mol type x y z ix iy iz 

compute	          rmsd  co2  msd  com  yes
fix               rmsd  co2  ave/time  ${sampint}   ${corrlen}   ${nevery}  c_rmsd[1]   c_rmsd[2]   c_rmsd[3]   c_rmsd[4]    file  co2_${count}.msd

compute           pp_co2 co2  stress/atom NULL
compute           p_co2 co2 reduce sum c_pp_co2[1] c_pp_co2[2] c_pp_co2[3]
variable          press_co2 equal -(c_p_co2[1]+c_p_co2[2]+c_p_co2[3])/(3*vol)
#fix               press_co2   co2  ave/time     ${sampint}   ${corrlen}   ${nevery}  v_press_co2  file  co2_${count}.press   mode  vector

compute 	myTemp co2 temp

# ==================================================================================================================

thermo_style      custom     step  pe  density  temp  c_myTemp  press  v_press_co2
thermo       	  ${ndelay}		

neighbor          2.0 bin
neigh_modify      exclude group framework framework
neigh_modify      every 1 delay 0 check yes

run          	  ${t_nvt}


write_restart     eq.co2

Following is the log error:

ERROR: Fix rigid: Bad principal moments (src/RIGID/fix_rigid_small.cpp:2322)
Last command: run          	  ${t_nvt}

I would be happy if you help me out with this.

Please let me know if you need additional information.

The problem is that you are using a minimization. That distorts or bends your supposed to be linear CO2 molecules and thus produces moments of inertia that show them as near-linear. Handling those as non-linear (i.e. prolate or oblate or quasi-spherical) would lead to numerical instablities or require an impractically short timestep, that LAMMPS chooses to error out instead.

Minimization with rigid objects is a “known to be hard” problem. The best way to approach this from within LAMMPS would be to do some kind of simulated annealing scheme, where you assign some temperature with the velocity command and then use fix viscous to drain energy from the system for a while and then switch to regular MD.
Or you just start with the MD right away but use the built in langevin thermostat of fix rigid with a rather short relaxation constant, i.e. large damping and randomization.

2 Likes