[gmx-users] Simulation of polarizable Carbon nanotubes

Thu Jun 25 02:31:24 CEST 2015

(Please don't CC me separately, just post to the list)

On 6/24/15 7:32 PM, S.M. Golam Mortuza wrote:
> Dear Gromacs Users,
>
> I am trying to run a simulation using polarizable force field for carbon
> nanotube (CNT), which has 304 carbon atoms. As a test purpose, I am trying to
> simulate a system containing a pair of KCl ion in an infinite carbon nanotube.
> However, my ultimately goal is to simulate a system comprising water, KCl inside
> the CNT.
>
> I have been trying to mimic the paper
> /http://www.sciencedirect.com/science/article/pii/S0927025607000456/
>
> Also, I have been following the thread on this topic:
> https://www.mail-archive.com/gromacs.org_gmx-users@maillist.sys.kth.se/msg10520.html
>
> Again, as a test purpose, I have put only one dummy atom (atom id 306), which is
> harmonically connected with a shell atom (atom id 305), inside a hexagonal ring
> (atom id: 1 20 38 39 57 58). I have used position restrain for the nanotube and
> the dummy atom. I think I've followed every procedure that Justin and Peter
> suggested in the thread.  Some parts of my itp file of the CNT looks like:
>
> /[ atoms ]
> ;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB
> chargeB      massB
> 1    CJ    1    CNT    C    1    0.25    12.011    ;    qtot    0/
> /2    CJ    1    CNT    C    2    0    12.011    ;    qtot    0
> 3    CJ    1    CNT    C    3    0    12.011    ;    qtot    0/
> /...................................................................................../
> /.....................................................................................
> /
> /305    SL    1    CNT    SL    1    -1.5    1    ;    qtot    0 ; Shell
> 306    VS    1    CNT    VS    300    0    0    ;    qtot    0 ; Dummy/
> /
> /
> /[ polarization ]
> ; virtual_site    shell    functiontype    alpha (in nm^3)
>      306        305    1        0.1296
> /
>
>
> //
> /[bonds]/
> /......../
>
> //
> /[angles]/
> /....../
>
> //
> /[dihedrals]/
> /....../
> /
> /
> /[ virtual_sitesn ]
> ; Site funct from
> 306 1 1 20 38 39 57 58/
>
> /[ exclusions ]
> ; iatom excluded from interaction with i
> 305 1    2    3    4    5    6    7    8    9  ..... 304/   ; I have used
> exclusions from the shell atom with all carbon atoms of the nanotube. I am not
> sure whether I should exclude more exclusions.
>
> Sometimes, the simulation gives error as:
>
> /EM did not converge in 20 iterations/
>
> /Segmentation fault (core dumped)/
>
>
> I tried to play with different parameters, such as polarizablity, exclusion
> list, etc. However, sometimes, I can see following messages in my window, and
> the simulation does not move further.
>
> /starting mdrun 'KCl electrolyte in CNT'
> 100000 steps,    100.0 ps.
> MDStep=    0/ 0 EPot: -4.83364532e+02, rmsF: 2.17e+01
>
> step 0MDStep=    1/ 0 EPot: -4.34285614e+02, rmsF: 2.19e+01
> MDStep=    2/ 0 EPot: -3.03536743e+02, rmsF: 2.22e+01
> MDStep=    3/ 0 EPot: -1.36375397e+02, rmsF: 2.24e+01
> MDStep=    4/ 0 EPot: 1.27241879e+01, rmsF: 2.25e+01
>
> NOTE: Turning on dynamic load balancing
>
> MDStep=    5/ 0 EPot: 1.01205765e+02, rmsF: 2.27e+01
> MDStep=    6/ 0 EPot: 1.15893883e+02, rmsF: 2.28e+01
> MDStep=    7/ 0 EPot: 7.45831146e+01, rmsF: 2.28e+01
> MDStep=    8/ 0 EPot: 1.53622913e+01, rmsF: 2.29e+01
> MDStep=    9/ 0 EPot: -2.12497330e+01, rmsF: 2.29e+01
> MDStep=   10/ 0 EPot: -1.05190954e+01, rmsF: 2.29e+01
> MDStep=   11/ 0 EPot: 4.68694420e+01, rmsF: 2.29e+01
> MDStep=   12/ 0 EPot: 1.29762100e+02, rmsF: 2.29e+01
> MDStep=   13/ 0 EPot: 2.02292053e+02, rmsF: 2.29e+01
> MDStep=   14/ 0 EPot: 2.34559998e+02, rmsF: 2.29e+01
> MDStep=   15/ 0 EPot: 2.15615555e+02, rmsF: 2.29e+01
> MDStep=   16/ 0 EPot: 1.57610504e+02, rmsF: 2.30e+01
> MDStep=   17/ 0 EPot: 8.87061996e+01, rmsF: 2.30e+01
> MDStep=   18/ 0 EPot: 4.16818733e+01, rmsF: 2.30e+01
> MDStep=   19/ 0 EPot: 3.34845123e+01, rmsF: 2.31e+01
> MDStep=   20/ 0 EPot: 6.07384720e+01, rmsF: 2.31e+01
> MDStep=   21/ 0 EPot:          inf, rmsF: 2.30e+01/
>
> I have checked that my shell and dummy atoms remain close to each other.
> However, sometimes the CNT blows up. I would really appreciate if anyone can
> help me with this.
>
> Justin,
> I have seen that your paper on polarizable force field has been published
> already. I am using gromacs version 4.6.1. Can I use the codes for polarizable
> force field in this version. Any sort of help would be highly appreciated.
>
> Please let me know, if you need any further information. Thanks in advance.
>

The only benefit from our code is the anisotropic polarization function.  You 
have defined an isotropic polarization in this model, which is not what the 
paper specifies.  Pull the latest master branch and do a 'git checkout drude' to 
get the code.  Build the manual for the syntax of the new topology directives.

Note that the extended Lagrangian formalism is incompatible with this model 
because the "parent atom" of the shell/Drude is actually a virtual site.  The 
Lagrangian will collapse spectacularly.  Stick with SCF.

De-complicate the system with some model system like benzene or something. 
Don't put in ions or anything yet.  The attachment of a shell to a dummy is a 
bit strange and may need to be considered differently in the code, but that's 
not something I can say with any degree of certainty.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================