[gmx-users] Simulating PMMA
Karel de Vries
kjcdevries at gmail.com
Sun Oct 2 20:06:07 CEST 2016
Mark,
I understand your concern.
I need to have the hydrogens for my final answer, so an all-atom force
field is definitely inevitable. Justification for choosing OPLS over the
other all-atom force fields is a little more tricky, but I have found
literature where they used OPLS-aa for similar purposes so I figured it
should work. As for the presence of a quaternary carbon, I expect that to
be the "alkane C", opls_139. Given that opls_135 through 138 are "alkane
CHn", this sounds like the most reasonable interpretation of opls_139.
I will try to find better documentation for the exact meaning of these atom
types though, thanks.
On Fri, Sep 30, 2016 at 1:49 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> On 9/30/16 7:14 AM, Karel de Vries 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.
>>
>>
> Nope, it's just based on the normal principal of trying to construct
> additive force fields from building blocks. Take what's known, assign a
> best match, otherwise default to some stock parameter, then clearly warn
> the user that the script can be dangerous :)
>
> -Justin
>
>
> 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.
>>>
>>>
> --
> ==================================================
>
> 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