[gmx-users] Simulating PMMA
Justin Lemkul
jalemkul at vt.edu
Fri Sep 30 13:49:47 CEST 2016
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
==================================================
More information about the gromacs.org_gmx-users
mailing list