[gmx-users] AlF4 - producing LINCS errors

[Matthew Cliff]

Hi,

     I'm very new to MD and I am trying to run MD simulations of 
proteins binding AlF4-. I am running GROMACS 4.5.3
and using OPLSAA forcefield, Tip4p water.
I have created a topology file for AlF4 (see below), and added atom 
types to the forcefield files atomtypes and ffnonbonded (see bottom of 
topology file). The force constants are pretty much fictitious, but I 
think reasonable.

If I run the AlF4- in a box of water with a counterion (Na+) with 
timestep 1fs and constraints=none then nothing untoward happens. 
However, if I run with 2 fs time step and constraints=all-bonds, then it 
goes horribly wrong after 69 steps and says:
"Step 69, time 0.138 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.211433, max 0.232176 (between atoms 2 and 5)
bonds that rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length

Step 70, time 0.14 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.558220, max 0.947733 (between atoms 2 and 4)
bonds that rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
       2      1  142.7    0.1413   0.1214      0.1790
       2      5  166.0    0.1374   0.1126      0.1790
       2      3  142.6    0.1416   0.1203      0.1790
Wrote pdb files with previous and current coordinates"

If I increase the bond strength this takes longer, but still happens. 
Energy minimisation works fine (of course?).

Obviously I would like to do half the calculations. (md.mdp is also at 
the bottom of this e-mail)

Cheers,

Matt Cliff



[ moleculetype ]
; Name nrexcl
ALF      3

[ atoms ]
; charges need to be checked
;   nr      type  resnr resid  atom  cgnr   charge     mass
      1  opls_998     1  ALF      F1     1      -0.8  18.99840   ;Al etc
      2  opls_997     1  ALF      AL     1       2.2  26.98153   ;Al etc
      3  opls_998     1  ALF      F2     1      -0.8  18.99840   ;Al etc
      4  opls_998     1  ALF      F3     1      -0.8  18.99840   ;Al etc
      5  opls_998     1  ALF      F4     1      -0.8  18.99840   ;Al etc

[ bonds ]
; ai  aj  fu    b0, kb, ...
    2   1   1    0.179    87864.0    0.179    87864.0 ;   Al-F
    2   3   1    0.179    87864.0    0.179    87864.0 ;   Al-F
    2   4   1    0.179    87864.0    0.179    87864.0 ;   Al-F
    2   5   1    0.179    87864.0    0.179    87864.0 ;   Al-F

[ angles ]
; ai  aj  ak  fu    c0, c1, ...
    1   2   4   1    90         397.5    90.0       397.5
    1   2   3   1    180         397.5    180.0       397.5
    1   2   5   1    90         397.5    90.0       397.5
    3   2   4   1    90         397.5    90.0       397.5
    3   2   5   1    90         397.5    90.0       397.5
    4   2   5   1    180         397.5    180.0       397.5

[ dihedrals ]
; ai  aj  ak  al  fu    c0, c1, m, ...
    2   3   1   4   2      0.0 1673.6        0.0 1673.6   ; imp making 
it flat
    2   1   4   5   2      0.0 1673.6        0.0 1673.6   ; imp making 
it flat

#ifdef POSRES
[ position_restraints ]
  1 1 1000 1000 1000
  2 1 1000 1000 1000
  3 1 1000 1000 1000
  4 1 1000 1000 1000
  5 1 1000 1000 1000
#endif


In atomtypes.atp:
opls_997    26.98153  ; Al
opls_998    18.99840  ; F

In ffnonbonded.itp:
opls_998   F   9 18.9984     -1.00        A    2.9500e-01  2.22e-01
opls_997  AL  13 26.98153     3.00        A    1.0e-01   3.66e+00 ;

;    Input file
;
title               =  General
cpp                 =  /usr/bin/cpp -traditional
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002    ; ps !
nsteps              =  10000  ; total 2 ns
nstcomm             =  1
nstcalcenergy       =  1
nstxout             =  500
nstvout             =  500
nstfout             =  0
nstxtcout           =  5000
nstlog              =  1000
nstenergy           =  1000
nstlist             =  5
ns_type             =  grid
rlist               =  1.2
rcoulomb            =  1.2
coulombtype         =  pme
vdwtype             =  cut-off
rvdw                =  1.2
optimize_fft        =  yes
DispCorr            =  EnerPres
; V-rescale temperature coupling is on in two groups
Tcoupl              =  v-rescale
tc-grps             =   Other Water_and_ions
tau_t               =   0.1    0.1
ref_t               =   300   300
; Energy monitoring
energygrps          =   Other Water_and_ions
; Isotropic pressure coupling is now on
Pcoupl              =  berendsen
Pcoupltype          =  isotropic
tau_p               =  1.0
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is off (from traj)
gen_vel             =  no
gen_temp            =  300.0
gen_seed            =  173529

-- 
Dr. Matthew J Cliff
Tel: (0114) 2222708
e-mail: m.j.cliff at sheffield.ac.uk
Molecular Biology and Biotechnology
University of Sheffield
Firth Court
Western Bank
Sheffield
S10 2TN
UK