[gmx-users] Simulation of polarizable Carbon nanotubes

[Jashimuddin Ashraf]

Dear users,

I want to perform a molecular dynamic simulation of polarizable carbon
nanotubes. I intend to implement this paper-

http://www.sciencedirect.com/science/article/pii/S0927025607000456

I digged up the manual but could not find much help from it. I went through
some mails in the gromacs user maillist, studied some .itp files and
learned some elementary stuffs regarding the addition of virtual sites and
shell atoms.I understand that I have to add both shell and virtual sites in
this case.

Now, before jumping right into a big nanotube molecule, I was trying to
perform a simulation with a single benzene ring with a virtual site placed
at the center and a shell attached to the virtual site.

In my forcefield.itp file have the virtual site and the shell declared like
this-

-------------------------------------------------------------------------------------------------------------------------------------
[ atomtypes ]
; name       mass      charge    ptype      sigma      eps
  CJ1   1     12.01100     0.25   A      3.40000e-01  3.60100e-01
  VS    1     0        0      D       0         0
  SL     1     0          -1.5      S      0        0

[ bondtypes ]
`; i    j     func  b0          kb        gamma
  CJ1    CJ1      3    0.1418   478.9000        21.867
  VS     SL       1    0.06       2409
-------------------------------------------------------------------------------------------------------------------------------------

Here, I have considered a bond between the virtual site and the shell (the
paper mentions something like it but does not provide with the bond
length). Is it a mistake?

and inside my topol.top file, I have-


-------------------------------------------------------------------------------------------------------------------------------------
[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass
typeB    chargeB      massB
     1        CJ1      1      C      C      1        0.25     12.011   ;
qtot 0.25
     2        CJ1      1      C      C      1        0.25     12.011   ;
qtot 0.5
     3        CJ1      1      C      C      1        0.25     12.011   ;
qtot 0.75
     4        CJ1      1      C      C      1        0.25     12.011   ;
qtot 1.00
     5        CJ1      1      C      C      1        0.25     12.011   ;
qtot 1.25
     6        CJ1      1      C      C      1        0.25     12.011   ;
qtot 1.50
     7         VS      1      C     VS      1        0          0      ;
qtot 1.50
     8         SL      1      C      S      1       -1.5          0    ;
qtot 0


[ polarization ]
; virtual_site    shell    functiontype    alpha (in nm^3)
 7        8    1        0.1


[ bonds ]
;  ai    aj funct            c0            c1            c2            c3
    1     2     3
    1     5     3
    2     4     3
    3     4     3
    3     6     3
    5     6     3
    7     8     1

[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1
c2            c3            c4            c5
    5     1     2     4     3
    2     1     5     6     3
    1     2     4     3     3
    6     3     4     2     3
    4     3     6     5     3
    1     5     6     3     3

[ virtual_sites3 ]
; detailed calculation not shown here
; Dummy from            funct    a        b
 7     4     5     6    1     0.5        0

-------------------------------------------------------------------------------------------------------------------------------------

I ran an energy minimization for an emtol of 100 but it gives me a result
like this-

Steepest Descents converged to machine precision in 141 steps,
but did not reach the requested Fmax < 100.
Potential Energy  = -6.1810545e+06
Maximum force     =  9.1656689e+12 on atom 4
Norm of force     =  4.5828344e+12

-------------------------------------------------------------------------------------------------------------------------------------
If I run a production MD with this, the simulation blows up with this-
MDStep=   40/18 EPot:          nan, rmsF:    nan
Warning: Only triclinic boxes with the first vector parallel to the x-axis
and the second vector in the xy-plane are supported.
         Box (3x3):
            Box[    0]={         nan,          nan,          nan}
            Box[    1]={         nan,          nan,          nan}
            Box[    2]={         nan,          nan,          nan}
         Can not fix pbc.
MDStep=   40/19 EPot:          nan, rmsF:    nan
step 40: EM did not converge in 20 iterations, RMS force nan

-------------------------------------------------------------------------------------------------------------------------------------
Is something wrong with my system itself? or is there anything wrong with
my methods?

I have been stuck at this for a very long time now and anything- any
comment or hint would be very much helpful for me.

Thanks in advance,
Jashimuddin Ashraf