[gmx-users] COM motion removal, LINCS error, PBC artifact?

[Leonardo Sepulveda Durán]

Hello!!!

I saw my protein was diffusing from the center of the waterbox to the
box edge, and then stayed there. That was rare because I used nstcomm
= 1 and comm_mode  = linear. Anyway, that simulations with diffusion
seemed to be ok, there was no warnings during the 5 ns simulation.

I added a new line to avoid COM movement of the solute

comm_grps = Protein SOL

As I read in the forum. It actually solved the problem of COM motion
(protein still rotates), Nevertheless the situation becomes
problematic because the following error becomes to iterate:

Step 600421, time 1200.84 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 0.001487 (between atoms 238 and 240) rms 0.000065
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
    238    239   31.5    0.1000   0.1000      0.1000

and after a while, this warning was added:

Step 600436  Warning: pressure scaling more than 1%, mu: 1.01362 1.01362 1.01362

after that appeared this other kind of warning:

Step 600460, time 1200.92 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 9.876248 (between atoms 575 and 576) rms 0.588436
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
    238    239   90.0    1.6668   1.0861      0.1000
    238    240   89.9    2.5654   0.9383      0.1000
    574    575   34.0    0.1396   0.1394      0.1390
    575    576   90.1    0.1501   1.1855      0.1090
    575    579   30.8    0.1395   0.1409      0.1390
    577    581   42.5    0.1396   0.1348      0.1390
    581    582   64.9    0.1095   0.0949      0.1090
    581    583   42.1    0.1396   0.1288      0.1390
Constraint error in algorithm Lincs at step 600460

and finaly, the output iterate in this kind of warning:

Constraint error in algorithm Lincs at step 600464
Correcting invalid box:
old box (3x3):
   old box[    0]={-3.66371e+04,  0.00000e+00, -0.00000e+00}
   old box[    1]={ 0.00000e+00, -3.89088e+04, -0.00000e+00}
   old box[    2]={ 0.00000e+00,  0.00000e+00, -4.15310e+04}
new box (3x3):
   new box[    0]={-3.66371e+04,  0.00000e+00, -0.00000e+00}
   new box[    1]={ 0.00000e+00, -3.89088e+04, -0.00000e+00}
   new box[    2]={ 0.00000e+00,  3.89088e+04, -4.15310e+04}

As the only changue I made was the comm_grps line, I guess the
centering procedure would be giving problems, or maybe another
changues have to be made to mdp file before using it. My mdp is the
following:

title                   = Simulated annealing
cpp                     = /lib/cpp
constraints             = all-bonds     ;  if dt = 0.002 ps
integrator              = md
;
;       molecular dynamics
;
dt                      = 0.002         ; ps !
nsteps                  = 2500000       ; 5 ns
comm_mode               = linear
nstcomm                 = 1
comm_grps               = Protein SOL
nstxout                 = 1000          ; collect data coordinates every 2 ps
nstvout                 = 1000          ; write velocities to output trajectory
nstfout                 = 0             ; write forces
nstenergy               = 10            ; write energies
;
;       non bonding setup
;
nstlist                 = 5             ; update neighbor list GROMACS
= 10, GROMOS96 = 5
ns_type                 = grid
rlist                   = 0.8
coulombtype             = cut-off
vdwtype                 = cut-off
rcoulomb                = 0.8           ; nm!
rvdw                    = 1.4           ; nm! GROMACS >= 0.9, GROMOS96 >= 1.4
;
Tcoupl                  = berendsen
tau_t                   = 0.1     0.1
tc-grps                 = Protein SOL
ref_t                   = 498     498
;
;       annealing options
;
;zero_temp_time         = 0             ; a ese tiempo se tendra 0 ªK
annealing               = single single
annealing_npoints       = 2      2
annealing_time          = 0 50   0 50
annealing_temp          = 1 498  1 498
;
;       Pressure coupling 
;
Pcoupl                  = berendsen
tau_p                   = 0.5
compressibility         = 4.5e-5        ; compresibilidad del H2O a
300K y 1 atm (en bar-1);
ref_p                   = 1.0
;
;       initial velocities
;
gen_vel                 = yes
gen_temp                = 1
gen_seed                = 173529

I dont know what is happening. CAn someone give some advice???

Ahhh... and another question. In non COM removal dynamics the protein
diffuses to the boundary, then begins to jump beetween box faces. But
it mantains there, it do not goes again to the box center, it seems to
be trapped for some periodic boubdary artifact. So I am not sure if
that kind of simulations are flawed. Or they are ok???. And if they
are Ok, there is some way to recenter protein in the box in that
simulations using some program after they are finished??

Thanks in advance


Thanks 

Leonardo