[gmx-users] Simulation of polarizable Carbon nanotubes
Justin Lemkul
jalemkul at vt.edu
Wed Mar 4 16:59:34 CET 2015
On 3/4/15 10:48 AM, Jashimuddin Ashraf wrote:
> Thanks for your reply Dr.Peter Kroon and Dr. Lemkul. The virtual_sitesn
> directive is very helpful in declaring the virtual site positions.
>
> Now, the simulation runs fine in vacuum. But as soon as I solvate the
> system in water, the production MD fails with-
>
What water model are you using?
> ----------------------------------------------------------------------------------------------------------------------------------------
> .
> .
> .
> MDStep= 91/ 9 EPot: -6.87315918e+02, rmsF: 1.75e+03
> MDStep= 91/10 EPot: -7.87661133e+02, rmsF: 1.44e+03
> MDStep= 91/11 EPot: -8.65109253e+02, rmsF: 1.22e+03
> MDStep= 91/12 EPot: -9.25257446e+02, rmsF: 1.06e+03
> MDStep= 91/13 EPot: -9.72217041e+02, rmsF: 9.53e+02
> MDStep= 91/14 EPot: -1.00904297e+03, rmsF: 8.72e+02
> MDStep= 91/15 EPot: -1.03802905e+03, rmsF: 8.10e+02
> MDStep= 91/16 EPot: -1.06091382e+03, rmsF: 7.61e+02
> MDStep= 91/17 EPot: -1.07902722e+03, rmsF: 7.23e+02
> MDStep= 91/18 EPot: -1.09339343e+03, rmsF: 6.93e+02
> MDStep= 91/19 EPot: -1.10480664e+03, rmsF: 6.70e+02
> step 91: EM did not converge in 20 iterations, RMS force 579.612
> MDStep= 92/ 0 EPot: -1.14934949e+03, rmsF: 7.17e+02
> Segmentation fault (core dumped)
> ----------------------------------------------------------------------------------------------------------------------------------------
>
> I found that at lower value of the polarizability (alpha_zz), the
> simulation goes fine (in water). The simulation also goes well if I
> decrease the charges of the carbon atoms and the shell. But at the values
> of alpha (0.1296 nm^3) and charges of the carbon atoms (0.25e) as the paper
> mentioned, the simulation ends with this error.
>
> In your message, you were kind to reply with-
>
> "If your simulation blows up, try analysing the distance between the
> shell and it's parent; it should give you a hint about what's going wrong."
>
> In my case, the distance between the shell and it's parent goes higher with
> the increase of the value of alpha. But I cannot understand what is causing
> the system to blow up. It would be very much helpful if you could help me
> with this error.
>
This makes sense. Alpha is the atomic polarizability, directly proportional to
the shell charge and inversely proportional to the force constant of the bond.
If you increase alpha with a fixed charge, the bond force constant is lower and
therefore you get a larger displacement.
-Justin
> thanks in advance.
>
> On Mon, Feb 23, 2015 at 6:39 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 2/23/15 3:47 AM, Peter Kroon wrote:
>>
>>>
>>> On 02/22/2015 05:33 PM, Jashimuddin Ashraf wrote:
>>>
>>>> Thanks for your reply Dr. Lemkul
>>>>
>>>> I changed the bond length to zero in my topology file, forcefield file
>>>> and
>>>> in the n2t file. I also added the [ exclusion] section-
>>>>
>>> If memory serves, you don't need a bond between a shell particle and
>>> it's "parent"; the polarization directive takes care of it's position.
>>> If you want, you can define a non-interacting connection (bond type 5)
>>> between them.
>>>
>>
>> This is correct. The [polarization] directive is just another bonded
>> interaction. The value of the force constant is back-calculated from the
>> atomic polarizability, alpha.
>>
>>
>>>> [ exclusions ]
>>>> ; iatom excluded from interaction with i
>>>> 1 2 3 4 5 6 7 8
>>>> 2 1 3 4 5 6 7 8
>>>> 3 1 2 4 5 6 7 8
>>>> 4 1 2 3 5 6 7 8
>>>> 5 1 2 3 4 6 7 8
>>>> 6 1 2 3 4 5 7 8
>>>> 7 1 2 3 4 5 6
>>>> 8 1 2 3 4 5 6
>>>>
>>>>
>>>> in my topology. But the same error keeps showing up. What did I possibly
>>>> do
>>>> wrong?
>>>>
>>> AFAIK this will exclude all atoms from eachother. Just the last line
>>> should be enough if you just want to exclude the shell-carbon interactions
>>>
>>>>
>>>> In your mail, you were kind to reply with-
>>>>
>>>> "I only looked at the paper briefly, but it seems they are working with a
>>>> model that makes use of anisotropic polarization. In GROMACS, this is
>>>> currently only available for water, so the model would not be supported."
>>>>
>>>> Are you suggesting that, this model is something we should not work with
>>>> in
>>>> GROMACS right now? (I am really sorry, I could not understand this part
>>>> properly)
>>>>
>>>
>> From my quick read of the paper (focusing only a narrow section of the
>> methods), the authors list alpha_zz only, implying that the shell is only
>> polarizable along the normal to the ring. Therefore the other elements of
>> the polarization tensor are zero.
>>
>>
>>>> Also, it would be very nice, if you could give us an idea about the time
>>>> required to publish the codes and adding them in the new versions of
>>>> GROMACS.
>>>>
>>>>
>> That depends on reviewers :)
>>
>> The paper should be submitted this week, after which I will be working to
>> submit the code to the master development branch. The code will not be
>> included in 5.1 (I missed the deadline, due to some debugging that took
>> longer than expected), but the code will be available as a patch to the
>> master branch in the coming weeks.
>>
>> -Justin
>>
>>
>> Thanks in advance,
>>>>
>>> Lastly, for defining the virtual site you could look at the
>>> virtual_sitesn directives and define a center of geometry or a center of
>>> mass. Im not sure if it has any advantages over a virtual site as it is
>>> defined now though.
>>> If your simulation blows up, try analysing the distance between the
>>> shell and it's parent; it should give you a hint about what's going wrong.
>>>
>>> Good luck!
>>>
>>> Regards,
>>> Peter Kroon
>>>
>>>>
>>>> On Sat, Feb 21, 2015 at 12:20 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>>
>>>>
>>>>> On 2/20/15 12:57 PM, Jashimuddin Ashraf wrote:
>>>>>
>>>>> 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?
>>>>>>
>>>>>>
>>>>>> The equilibrium bond length should be zero (i.e. no induced
>>>>> polarization).
>>>>>
>>>>> 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
>>>>>>
>>>>>>
>>>>>> I only looked at the paper briefly, but it seems they are working
>>>>> with a
>>>>> model that makes use of anisotropic polarization. In GROMACS, this is
>>>>> currently only available for water, so the model would not be supported.
>>>>> It will be soon (I know I've been saying that for a while, but our paper
>>>>> regarding Drude simulations in GROMACS is just about done, after which I
>>>>> can provide the code).
>>>>>
>>>>>
>>>>> [ 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
>>>>>>
>>>>>>
>>>>>> You're missing a critical element here; the paper says that the shell
>>>>> does
>>>>> not interact with the carbon atoms of the ring, so you need to define
>>>>> exclusions manually.
>>>>>
>>>>> -Justin
>>>>>
>>>>>
>>>>> ------------------------------------------------------------
>>>>>
>>>>>> ------------------------------------------------------------
>>>>>> -------------
>>>>>>
>>>>>> 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
>>>>>>
>>>>>>
>>>>>> --
>>>>> ==================================================
>>>>>
>>>>> 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
>>>>>
>>>>> ==================================================
>>>>> --
>>>>> Gromacs Users mailing list
>>>>>
>>>>> * Please search the archive at http://www.gromacs.org/
>>>>> Support/Mailing_Lists/GMX-Users_List before posting!
>>>>>
>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>>
>>>>> * For (un)subscribe requests visit
>>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>>>> send a mail to gmx-users-request at gromacs.org.
>>>>>
>>>>>
>>>
>>>
>>>
>>>
>> --
>> ==================================================
>>
>> 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
>>
>> ==================================================
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at http://www.gromacs.org/
>> Support/Mailing_Lists/GMX-Users_List before posting!
>>
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-request at gromacs.org.
>>
--
==================================================
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
==================================================
More information about the gromacs.org_gmx-users
mailing list