[gmx-users] Parameter suggestion for locating an Iron atom between two carbon ring and preventing system blow up

[Sajjad Kavyani]

Dear experts,

I'm trying to simulate a box of benzene with 25~35 ferrocene (C5H5-Fe-C5H5)
in it.
I have the benzene itp, but in the case of ferrocene I have a problem to
choose suitable parameter.
Each parameters I tried for ferrocene lead the system to blow up.
EVEN in a empty box (6 6 6) with just 20 or 30 ferrocene (pure ferrocene)
my system blows up in just 40000 steps (time step = 0.002)
For pure ferrocene, energy minimization hardly reach forces below 1000 and
just in the case of changing the Fe bonds from 1 to harmonic (type 6) I
could manage to get below zero potential values. But again it hardly
reaches below 1000!
- The procedure before blowing up was:
EM - NVT - NPT (in this step, blows up)
EM - NVT (with posre) - NPT (with posre) - NPT (no posre && in this step.
blows up)
- For t/p coupling the 'berendsen' method was utilized. And constraints was
set to 'none'.
The following is the itp for ferrocene was used.

Is it correct to define 10 bonds to keep Fe in the middle of the two rings?
What should I do to prevent the blow up?


This is atomtype.itp of Fe2+:
[ atomtypes ]

    Fe2+     6      55.84700     0.000 A    0.115816833358 0.0

######
And, this is the ferrocene itp (OPLS ff):


[ moleculetype ]
; molname       nrexcl
FRC                 3


[ atoms ]
;id    type     resnr  residu  a.name  cgnr  charge        mass
  1    opls_714  1    FRC       CP1    1     -0.356        12.01100
  2    opls_715  1    FRC       HP1    2      0.156         1.00800
  3    opls_714  1    FRC       CP2    3     -0.356        12.01100
  4    opls_715  1    FRC       HP2    4      0.156         1.00800
  5    opls_714  1    FRC       CP3    5     -0.356        12.01100
  6    opls_715  1    FRC       HP3    6      0.156         1.00800
  7    opls_714  1    FRC       CP4    7     -0.356        12.01100
  8    opls_715  1    FRC       HP4    8      0.156         1.00800
  9    opls_714  1    FRC       CP5    9     -0.356        12.01100
  10   opls_715  1    FRC       HP5    10     0.156         1.00800
  11   Fe2+      1    FRC       FeS    11     2.000        55.84700
  12   opls_714  1    FRC       CP6    12    -0.356        12.01100
  13   opls_715  1    FRC       HP6    13     0.156         1.00800
  14   opls_714  1    FRC       CP7    14    -0.356        12.01100
  15   opls_715  1    FRC       HP7    15     0.156         1.00800
  16   opls_714  1    FRC       CP8    16    -0.356        12.01100
  17   opls_715  1    FRC       HP8    17     0.156         1.00800
  18   opls_714  1    FRC       CP9    18    -0.356        12.01100
  19   opls_715  1    FRC       HP9    19     0.156         1.00800
  20   opls_714  1    FRC       CP0    20    -0.356        12.01100
  21   opls_715  1    FRC       HP0    21     0.156         1.00800

[ bonds ]
;  i     j     funct   length    force.c.
;;CyclPen1
   1 2      1       0.105     307105.6
   1 3      1       0.138     392459.2
   3 4      1       0.105     307105.6
   3 5      1       0.138     392459.2
   5 6      1       0.105     307105.6
   5 7      1       0.138     392459.2
   7 8      1       0.105     307105.6
   7 9      1       0.138     392459.2
   9 10     1       0.105     307105.6
   9 1      1       0.138     392459.2
;;CyclPen2
   12 13     1       0.105     307105.6
   12 14     1       0.138     392459.2
   14 15     1       0.105     307105.6
   14 16     1       0.138     392459.2
   16 17     1       0.105     307105.6
   16 18     1       0.138     392459.2
   18 19     1       0.105     307105.6
   18 20     1       0.138     392459.2
   20 21     1       0.105     307105.6
   20 12     1       0.138     392459.2
;;Iron with carbons:
   11 1      6       0.202     300500
   11 3      6       0.202     300500
   11 5      6       0.202     300500
   11 7      6       0.202     300500
   11 9      6       0.202     300500
   11 12     6       0.202     300500
   11 14     6       0.202     300500
   11 16     6       0.202     300500
   11 18     6       0.202     300500
   11 20     6       0.202     300500

[ angles ]
; i     j k       funct     angle       force.c.
;;CyclPen1
   1     3 5       2         108         10.0
   3     5 7       2         108         10.0
   5     7 9       2         108         10.0
   7     9 1       2         108         10.0
   9     1 3       2         108         10.0
;;CyclPen2
   12     14 16       2         108         10.0
   14     16 18       2         108         10.0
   16     18 20       2         108         10.0
   18     20 12       2         108         10.0
   20     12 14       2         108         10.0