LINCS Warning with linear CO2 virtual site

GROMACS version: 2020.1-Ubuntu-2020.1-1
GROMACS modification: No
Here post your question:
So I am doing a simulation of protein solvated in gas - CO2 RCSB PDB - 5WSK: Structure of Ribulose-1,5-bisphosphate carboxylase/oxygenase from wheat. I took topology and PDB of this gas from Virtual Sites. I took the PDB model from PDB and then removed water and all non-standard residues.

Then I ran these gmx commands using the files provided. (I used OPLS/AA force field)

gmx pdb2gmx -f rubisco_mfix.pdb -o 1.gro -p topol.top -water none
gmx editconf -f 1.gro -o box.gro -c -d 3.0 -bt dodecahedron
gmx insert-molecules -f box.gro -ci co2.pdb -nmol 40 -try 20 -o boxco2.gro

Then I included topol.itp (which is itp of co2 with ff included) in topology (topol.top) and I added CO2 molecules to the system.
I didn’t neutralize this system, because I think that you don’t need to neutralize gas solvated systems (If I have to, how can I do it?) (With all gmx grommp I was getting a Warning about this, so I used -maxwarn 1)
Then I ran these gmx commands (mdp files included) :

gmx grompp -f minimize.mdp -c water_ions.gro -p topol.top -o energymin.tpr
gmx mdrun -v -s energymin.tpr -deffnm em
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt

All of these runs were successful. Then I attempted NPT run and it returned this error :

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt
Step 0, time 0 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 1.241085, max 56.916912 (between atoms 70025 and 70026)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
  69980  69981   90.0    0.2132  11.0419      0.2132
  69985  69986   90.0    0.2132   6.3848      0.2132
  69990  69991   93.1    0.2132  10.8002      0.2132
  69995  69996   99.1    0.2132   9.3370      0.2132
  70000  70001  100.8    0.2132   8.7460      0.2132
  70010  70011  113.2    0.2132   7.8905      0.2132
  70015  70016   90.0    0.2132  11.0835      0.2132
  70020  70021   90.0    0.2132   8.9943      0.2132
  70025  70026   90.0    0.2132  12.3463      0.2132
  70030  70031   90.0    0.2132   9.1664      0.2132
  70035  70036   90.0    0.2132   5.8325      0.2132
  70040  70041   90.0    0.2132   3.0451      0.2132
  70045  70046   90.0    0.2132   7.3950      0.2132
  70050  70051   90.0    0.2132   4.0567      0.2132
  70055  70056   90.0    0.2132  11.7571      0.2132
  70060  70061   90.0    0.2132   7.8807      0.2132
  70070  70071   90.0    0.2132   9.5708      0.2132
  70075  70076   90.0    0.2132   6.0452      0.2132
  70080  70081   90.0    0.2132   7.4494      0.2132
  70085  70086   90.0    0.2132   2.7008      0.2132
  70090  70091   90.0    0.2132   5.1878      0.2132
  70095  70096   94.6    0.2132  10.1661      0.2132
  70100  70101   90.0    0.2132   5.2383      0.2132
  70105  70106  114.0    0.2132   7.9965      0.2132
  70110  70111   90.0    0.2132   2.7225      0.2132
  70115  70116  100.1    0.2132  10.9836      0.2132
  70120  70121   90.0    0.2132   4.2884      0.2132
  70125  70126   90.0    0.2132   5.1484      0.2132
  70130  70131   90.0    0.2132   7.7973      0.2132
  70135  70136   90.0    0.2132   9.1673      0.2132
  70140  70141   90.0    0.2132   7.0320      0.2132
  70145  70146   90.0    0.2132   6.3205      0.2132
  70150  70151   90.0    0.2132   7.4730      0.2132
  70155  70156   90.0    0.2132  11.0906      0.2132
  70160  70161   90.0    0.2132   2.3694      0.2132
  70165  70166   90.0    0.2132  10.8323      0.2132
  70170  70171   90.0    0.2132   9.2923      0.2132
  70175  70176   90.0    0.2132   8.7614      0.2132
  70175  70176   90.0    0.2132   8.7614      0.2132
  70175  70176   90.0    0.2132   8.7614      0.2132
  70175  70176   90.0    0.2132   8.7614      0.2132

Atoms given in this message are atoms of my CO2 Virtual Sites (checked in .gro file)Then it produced other LINCS errors in further steps.
My question is - how do I solve this problem? I am pretty sure that problem lies in my CO2 virtual sites model.

; minimize.mdp
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = PME
rvdw            = 1.0
rlist           = 1.0
rcoulomb        = 1.0
fourierspacing  = 0.12
pme-order       = 4
ewald-rtol      = 1e-5

;nvt.mdp
define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 50000     ; 2 * 50000 = 100 ps
dt                      = 0.001     ; 2 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Nonbonded settings 
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
tau_t                   = 0.1     0.1           ; time constant, in ps
ref_t                   = 300     300           ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl                  = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 300       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed

;npt.mdp 
define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 50000     ; 2 * 50000 = 100 ps
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = yes       ; Restarting after NVT 
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Nonbonded settings 
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
tau_t                   = 0.1     0.1           ; time constant, in ps
ref_t                   = 300     300           ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
pcoupltype              = isotropic             ; uniform scaling of box vectors
tau_p                   = 2.0                   ; time constant, in ps
ref_p                   = 1.0                   ; reference pressure, in bar
compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
refcoord_scaling        = com
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = no        ; Velocity generation is off 


;topol.itp
#include "oplsaa.ff/forcefield.itp"

; Here we introduce a special atom type to be used for each mass center
[ atomtypes ]
; name  bond_type    mass    charge   ptype          sigma      epsilon
  MCO   MCO          0.000   0.000    A              0.000      0.000

[ moleculetype ]
; name  nrexcl
CO2     2

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB
     1   opls_272      1    CO2     O1      1    -0.350      0.0000 
     2   opls_271      1    CO2      C      1     0.700      0.0000
     3   opls_272      1    CO2     O2      1    -0.350      0.0000
     4        MCO      1    CO2     M1      1     0.000     22.0049     
     5        MCO      1    CO2     M2      1     0.000     22.0049

[ constraints ]
; There are no bonds in this system
; Instead, we fix the distance between the mass centers such that
; the virtual sites can be reconstructed
    4   5   1   0.213173

[ virtual_sites2 ]
; the M--O distance is 0.125 - 0.1065865 = 0.0184135
; the M--M distance is 0.213173
; therefore, the fraction of the distance along the M--M length is (0.213173+0.0184135)/0.213173 = 1.0851116
; thus placing the virtual O sites beyond the M--M distance 
; site  ai  aj  funct   a
     1   4   5      1   1.0851116   ; relative to mass center 4, extends beyond mass center 5 
     2   4   5      1   0.5000      ; right in the middle
     3   5   4      1   1.0851116   ; as in the case of site 1