[gmx-users] Simulating PMMA

[Karel de Vries]

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