[gmx-users] Simulating PMMA
Karel de Vries
kjcdevries at gmail.com
Fri Sep 30 13:14:34 CEST 2016
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.
>
More information about the gromacs.org_gmx-users
mailing list