[gmx-users] Simulating PMMA

Mark Abraham mark.j.abraham at gmail.com
Fri Sep 30 13:22:29 CEST 2016


Hi,

How did you come to decide to model with OPLS/AA? As part of that process,
you should be considering whether it parameterized a quaternary carbon,
e.g. by inspecting its literature, or reference implementation in BOSS.
Availability of convenient tools is only one of many aspects to consider in
modelling physics...

Mark

On Fri, Sep 30, 2016 at 1:15 PM Karel de Vries <kjcdevries at gmail.com> wrote:

> Hi Justin,
>
> Thanks for your answer; I assumed that TopolGen is not intentionally
> incorrect. What I meant to ask was whether it was intentionally different
> from what the comments in the atomtypes.atp file indicate. I thought that
> perhaps, based on past experience, you know that some choices work better
> than others and included this into your script. I would not want to undo
> this by cooking up my own assignments based on no past experience at all.
>
> Karel
>
> On Fri, Sep 30, 2016 at 12:00 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
> >
> >
> > On 9/30/16 3:35 AM, Karel de Vries wrote:
> >
> >> Hi again,
> >>
> >> I think the problem is in the force field description that I put in the
> >> aminoacids.rtp file. Like I said, I just got the output from Justin's
> >> TopolGen (V1.1) script. Upon comparison of the opls_xxx choices to the
> >> comments in the atomtypes.atp file, a few strange assignments turn up.
> For
> >> example, the Os are taken to be from the C=O in an amide. However, there
> >> is
> >> no nitrogen in PMMA at all. There is a C=O in PMMA, but it's an ester
> >> group. There are a few more doubtful choices, e.g. a C bound to four
> other
> >> Cs is assigned the "alkane CH3" atom.
> >>
> >> I have two questions.
> >> One concerns the TopolGen script: are these seemingly incorrect
> >> assignments
> >> intentional? Or is this just a limitation in the script? There is
> another
> >> TopolGen topology for PMMA online, at
> >> http://gromacs-oplsaa-topologies.wikia.com/wiki/PMMA_3-mer, which has
> >> slightly different but similarly doubtful assignments.
> >>
> >>
> > I would never write a program that is intentionally incorrect.  TopolGen
> > tries to guess the best parameter based on connectivity, and sometimes it
> > guesses wrong, because it is pretty much impossible to design a script
> that
> > will account for the vast chemical space that OPLS-AA encompasses.  I try
> > to make this abundantly clear in the final message printed to the screen
> by
> > TopolGen:
> >
> > "Output topology has been written. An attempt has been made to assign
> > charges and atom types based on existing functional groups, but they may
> > not be correct.  No charge calculations or other parameterization
> > calculations have been done. Guesses have been made for charge groups.
> > Please inspect and correct the topology before using it in any
> > simulations.  The author of the script does NOT guarantee accuracy or
> > usability of any of the content; TopolGen was written as a convenience
> for
> > outputting a skeleton topology, and nothing more."
> >
> > -Justin
> >
> >
> > The second is about the opls-aa force field. Is there any other
> >> documentation on which assignments to make than just the comments in the
> >> atomtypes.atp file?
> >>
> >>
> >> Regarding Chaofu Wu's reply to my original email; thank you very much
> for
> >> your comments. I accidentally sent a template file, not one used for an
> >> actual run. I run gromacs with lower nst* values; when I need
> statistics I
> >> make sure there are about 1000 frames in the energy file. Values I
> adjust
> >> from the parameters I sent are nsteps, nst*, ref_t, ref_p, and for the
> >> initial run also gen_vel / gen_temp. I'm sorry for the confusion.
> >>
> >> Karel de Vries
> >>
> >>
> >> On Wed, Sep 28, 2016 at 11:43 AM, Karel de Vries <kjcdevries at gmail.com>
> >> wrote:
> >>
> >> Hello all,
> >>>
> >>> I'm trying to simulate bulk PMMA, poly(methyl methacrylate) in Gromacs.
> >>> I'm using the opls-aa force field. I added a PMMA unit cell to the
> >>> aminoacids.rtp file (listing below). I'm running a simulation with 30
> >>> molecules of 50 unit cells each. I start by putting straight strands
> into
> >>> the box with gmx insert-molecules.
> >>>
> >>> I think that my input parameters are incorrect. My first reason for
> >>> thinking this is that it is very difficult to equilibrate this system:
> it
> >>> takes a lot of computation time to get the internal pressure to equal
> the
> >>> external one during NPT even at high temperatures. The second reason is
> >>> that I find unrealistic values for the glass transition temperature and
> >>> density of PMMA at room temperature.
> >>>
> >>> I have tried several combinations of NVT and NPT simulations at
> different
> >>> temperatures to get the system to equilibrate. The most successful one
> is
> >>> as follows. I start with an NPT simulation at 1000 bar of isotropic
> >>> pressure at 900 K using Berendsen barostat with tau_p=5 and
> >>> compressibility
> >>> 4.5e-5. Then I slowly decreased the pressure to 1 bar in 8 steps of 10
> ns
> >>> each. Finally I ran a 40 ns NPT equilibration step at 1 bar at the end.
> >>> The
> >>> average pressure tensor over this 40 ns run is as follows:
> >>>    Pressure (bar)
> >>>     9.67017e-01   -4.24140e-01   -1.51426e+00
> >>>    -4.24137e-01   -2.97605e+00   -3.38247e+00
> >>>    -1.51425e+00   -3.38246e+00    2.07151e+00
> >>> I'm not sure how good this is, but it's much better than anything else
> >>> I've done.
> >>>
> >>> I then cool down from 900K to 550K in steps of 50K, 5 ns per
> simulation.
> >>> I
> >>> do not expect the average cooling rate (0.01 K/ps) to be
> unrealistically
> >>> fast: 550K is well above the glass transition temperature of real PMMA
> >>> and
> >>> I sometimes see people in literature use similar rates during annealing
> >>> in
> >>> other polymers. I then cool down from 550K to 300K in steps of 10K, 5
> ns
> >>> per simulation. The average pressure tensor contains values that keep
> >>> getting further away from zero. Moreover, I find a glass transition
> >>> temperature of around 470K, while for real PMMA it should be around
> 400K.
> >>> The final density is 1.01 g/cm3, while it should be around 1.18 g/cm3.
> >>>
> >>> I'm starting to think there is something wrong with my parameters here,
> >>> either in the mdp file or in the force field parameters. Has anybody
> here
> >>> got any experience with PMMA and can they tell me what I'm doing wrong?
> >>>
> >>> Thanks a lot in advance.
> >>> Karel de Vries
> >>>
> >>>
> >>>
> >>>
> >>>
> >>> This is what I added to the .rtp file in the force field. I have
> similar
> >>> definitions for the end chains, which have a H connected to C1 and C2,
> >>> respectively. This was based on output by Justin Lemkul's TopolGen
> >>> script,
> >>> except I took the charges from an external source (don't remember
> where I
> >>> found them, unfortunately).
> >>> [ atoms ]
> >>>      C1    opls_135   -0.376     1
> >>>      C2    opls_135    0.009     1
> >>>      H3    opls_140    0.212     1
> >>>      H4    opls_140    0.155     1
> >>>      C5    opls_145    0.650     2
> >>>      C6    opls_135   -0.288     3
> >>>      H7    opls_140    0.096     3
> >>>      H8    opls_140    0.096     3
> >>>      H9    opls_140    0.096     3
> >>>     O10    opls_236   -0.286     2
> >>>     C11    opls_166   -0.161     2
> >>>     O12    opls_236   -0.563     2
> >>>     H13    opls_140    0.120     2
> >>>     H14    opls_140    0.120     2
> >>>     H15    opls_140    0.120     2
> >>>
> >>> [ bonds ]
> >>>     C1    C2
> >>>     C1    H3
> >>>     C1    H4
> >>>    C11   H13
> >>>    C11   H14
> >>>    C11   H15
> >>>    C11   O10
> >>>     C6    H7
> >>>     C6    H8
> >>>     C6    H9
> >>>     C2    C5
> >>>     C2    C6
> >>>     C5   O10
> >>>     C5   O12
> >>>     C2   +C1
> >>>
> >>>
> >>>
> >>>
> >>> These are typical parameters I've got in the mdp file:
> >>> integrator               = md
> >>> tinit                    = 0
> >>> dt                       = 0.001
> >>> nsteps                   = 1000000 ; 1 ns
> >>>
> >>> emtol                    = 100
> >>> emstep                   = 0.001
> >>> niter                    = 20
> >>> fcstep                   = 0
> >>> nstcgsteep               = 1000
> >>> nbfgscorr                = 10
> >>>
> >>> cutoff-scheme            = Verlet
> >>> nstxout                  = 1000000
> >>> nstvout                  = 1000000
> >>> nstfout                  = 1000000
> >>> nstlog                   = 2000000
> >>> nstenergy                = 10000
> >>> nstxout-compressed       = 10000
> >>> compressed-x-precision   = 1000
> >>> rcoulomb                 = 1.00
> >>> nstlist                  = 45
> >>> ns_type                  = grid
> >>> pbc                      = xyz
> >>> rlist                    = 1.00
> >>> vdw-type                 = cut-off
> >>> vdw-modifier             = Potential-switch
> >>> rvdw-switch              = 0.8
> >>> rvdw                     = 1.00
> >>>
> >>> Tcoupl                   = v-rescale
> >>> tc-grps                  = System
> >>> tau_t                    = 0.2
> >>> ref_t                    = 900
> >>> Pcoupl                   = berendsen
> >>> Pcoupltype               = isotropic
> >>> tau_p                    = 5.0
> >>> compressibility          = 4.5e-5
> >>> ref_p                    = 1.0
> >>>
> >>> gen_vel                  = no
> >>> gen_temp                 = 900
> >>> gen_seed                 = 1000
> >>>
> >>> DispCorr                 = EnerPres
> >>>
> >>>
> >>>
> >>>
> >>>
> > --
> > ==================================================
> >
> > 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.
> >
> --
> 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