[gmx-users] Simulation of polarizable Carbon nanotubes
Jashimuddin Ashraf
jashimuddin.ashraf23 at gmail.com
Mon Mar 23 15:59:46 CET 2015
Thanks again for your reply Dr. Lemkul.
I was wondering about the update regarding publishing the codes. Are the
codes ready to be published? Actually I could not get my work done and it
seems I have to wait until the codes/ patch are available.
Thanks in advance.
On Wed, Mar 4, 2015 at 10:13 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> On 3/4/15 11:07 AM, Peter Kroon wrote:
>
>>
>>
>> On 03/04/2015 04:48 PM, Jashimuddin Ashraf wrote:
>>
>>> Thanks for your reply Dr.Peter Kroon and Dr. Lemkul. The virtual_sitesn
>>>
>> Im not a Dr. yet ;)
>>
>>> 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-
>>>
>>> ------------------------------------------------------------
>>> ------------------------------------------------------------
>>> ----------------
>>> .
>>> .
>>> .
>>> 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.
>>>
>> See Justin's answer for an explanation.
>> As for blowing up: if the shell particle moves from it's parent far
>> enough to reach a positively charged particle the attraction will become
>> VERY large, blowing up your system.
>> Check your exclusion list with the original paper. It could be you need
>> to exclude more. For example: CHARMM treats drude/shell particles in
>> such a way that if two atoms are in a 1-4 relationship, their
>> drude/shell particles are as well. According to GROMACS, these would
>> have a 1-6 relationship.
>>
>>
> ...all of which are normal nonbonded relationships. Note that our Drude
> force field does Thole screening for 1-2 and 1-3 interactions (based on the
> parent atoms, not the Drudes, which in a 1-3 interaction of parent atoms
> would actually be 1-5 interaction for Drudes without this special
> classification), but anything beyond that is just "normal."
>
> -Justin
>
>
> Good luck!
>> Peter Kroon
>>
>>>
>>> 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
>
> ==================================================
> --
> 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.
>
More information about the gromacs.org_gmx-users
mailing list