Problems in fix rigid: Bad principal moments

LAMMPS LAMMPS Beginners
1

Hello everyone,

First of all, I would like to apologise if this is a trivial mistake or can be easily solved by some functionality that I’m not aware of.

So, essencially I have some CO2 adsorption GCMC results from RASPA that I am trying to use as a starting point in my MD simulations. From the topics that I have seen in this forum, most likely my issue is with the starting PDB file that is given by RASPA, which doesn’t provide enough accuracy to ensure a linear CO2 molecule for fix rigid.
I have already tried to use fix viscous to drain energy from the system, as well as some other suggestions provided in this forum.
The reason why I am not using create_atoms is because I want to start from the GCMC results. There are other more complex systems with highly saturated binary mixtures that will make it really hard to introduce molecules without overlapping.

Here is my input script:

log             log.lammps append
units           real
atom_style      full
boundary        p p p

pair_style      hybrid/overlay hbond/dreiding/lj 4 6 6.5 90 lj/cut/coul/long 12
bond_style      harmonic
angle_style     harmonic
dihedral_style  charmm
improper_style  umbrella

special_bonds   dreiding
special_bonds   lj/coul 0.0 0.0 1.0
pair_modify     tail yes mix arithmetic
dielectric      1.0
box tilt        large
read_data       data.minimized
include         parm.lammps

kspace_style    pppm 1e-6

neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes

#---------------------#

group           COF     type 1:9
group           CO2     type 10 11

#---------------------#

#molecule co2mol ../co2.mol

timestep 1

thermo_style yaml
thermo 1000

velocity all create 298 12359
run 0
velocity all scale 298

reset_timestep 0

restart 1000 a.tmp.restart b.tmp.restart
dump mydmp all atom 1000 dump.lammpstrj

compute ctCO2 CO2 temp
compute ctCOF COF temp
fix 1 CO2 rigid/nve/small molecule langevin 298 298 $(100*dt) 129840
fix 3 COF nve
fix 4 COF langevin 298 298 $(100*dt) 813247

run 50000

write_data data.nve
write_restart restart.nve

shell rm a.tmp.restart b.tmp.restart

I am starting from a minimized data, but it doesn’t work for the initial (pre-minimization) data as well. In the minimization I used a big potential to try to force the linearity of the molecule.
I also tried running nve without rigid before the final run, but everything just freezes before starting that final run.

Thank you for your attention, every help is welcome!

Márcio Soares

2

If precision of the geometries is the issue at hand, then you can circumvent the error by changing the value for the EPSILON constant in the src/RIGID/rigid_const.h file and recompile LAMMPS.

3

Alternatively, if you are not happy with how close the values of your angles are from 180°, one thing that comes to my mind is the following: you could make a python (or whatever language) code to scan the position of all the atoms in your “data.minimized” data file and identify, for each and all CO2 molecules individually, the given carbon and the two given oxygens. Let’s say one oxygen is point A, the carbon in the middle is point B and the other oxygen is point C. You choose one of these oxygens to be the “reference oxygen” (let’s say) and define a vector u = AB. The “corrected position” of the other oxygen should then be given by a vector AC = 2u. Or even you can just do B + u to find the corrected coordinates of point C.

image

If you are not so happy with the bond lengths either, you could scale the vector u accordingly beforehand to have the proper magnitude.

1 Like
4

Thank you for the suggestion!

I actually thought about doing something like that, but since my initial unreplicated file was just around 40 molecules of CO2, I decided to just convert the PDB file to XYZ and forced the angles at 180º.

Nonetheless, your responses were much apreciated and I might use this scripting idea for some other systems in the future!

Best wishes,
Márcio