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