# Very slow performance for an OPLS-AA methanol simulation

I am running a methanol simulation with a topology I obtained from Tinker for OPLS-AA, at 298 K and 1 bar. This is my .itp file:

[ atomtypes ]
; atm_kind   symbol      mass    charge      ptype   sigma   epsilon
MTH0      C01    12.0110     0.000       A   0.35000   0.27614
MTH1      O02    15.9990     0.000       A   0.31200   0.71128
MTH2      H03     1.0080     0.000       A   0.25000   0.12552
MTH3      H04     1.0080     0.000       A   0.25000   0.12552
MTH4      H05     1.0080     0.000       A   0.25000   0.12552
MTH5      H06     1.0080     0.000       A   0.00000   0.00000

[ moleculetype ]
; residue    nrexcl
MTH      3

[ atoms ]
; srno   atm_kind    resnr   res     sym     cgnr    charge      mass
1         MTH0      1    MTH      C01      1     0.145    12.0110
2         MTH1      1    MTH      O02      1    -0.683    15.9990
3         MTH2      1    MTH      H03      1     0.040     1.0080
4         MTH3      1    MTH      H04      1     0.040     1.0080
5         MTH4      1    MTH      H05      1     0.040     1.0080
6         MTH5      1    MTH      H06      1     0.418     1.0080

[ bonds ]
; p1     p2      funct   r0      ks
1     2     1      0.1410 267776.000
1     3     1      0.1090 284512.000
1     4     1      0.1090 284512.000
1     5     1      0.1090 284512.000
2     6     1      0.0945 462750.400

[ angles ]
; p1     p2      p3      funct   theta0      kb
2     1     3     1    109.500    292.880
2     1     4     1    109.500    292.880
2     1     5     1    109.500    292.880
3     1     4     1    107.800    276.144
3     1     5     1    107.800    276.144
4     1     5     1    107.800    276.144
1     2     6     1    108.500    460.240

[ dihedrals ]
; PROPER DIHEDRAL ANGLES
; p1     p2      p3      p4      funct   c0      c1      c2      c3      c4      c5
3    1    2    6      3       0.736    2.209    0.000   -2.946    0.000    0.000
4    1    2    6      3       0.736    2.209    0.000   -2.946    0.000    0.000
5    1    2    6      3       0.736    2.209    0.000   -2.946    0.000    0.000


When I run the npt equilibration step before my production run, I get the following note:

NOTE 1 [file unknown]:
You are using constraints on all bonds, whereas the forcefield has been
parametrized only with constraints involving hydrogen atoms. We suggest
using constraints = h-bonds instead, this will also improve performance.


These are my constraints in my .mdp file:

continuation            = yes       ; Restarting after NPT
constraint_algorithm    = lincs     ; holonomic constraints
constraints             = all-bonds   ; bonds involving H are constrained



Furthermore, my performance is pretty slow. I have 3000 methanol molecules, so 18000 atoms, and I am taking upwards of 24 hours to get 50 ns. What could be the cause of such a slow down? Methanol is a small, well-studied molecule and I want to benchmark TINKER before I move forward.

Any advice you have would be appreciated.

• Did you try changing constraints to h-bonds? Sep 19 at 6:27
• What hardware are you running on? I once spent a couple days debugging performance on the "national flagship supercomputer" of a certain European country before I realized the performance issues were because the computer in question was old and due for replacement. (Also, I agree with @TristanMaxson that you should use h-bonds instead of all-bonds -- the comment in your mdp doesn't match the statement.) Sep 19 at 7:19
• yes, I will be trying changing to h-bonds. But my question is, what if i want to have all-bonds to be constrained? Sep 19 at 23:41
• ...Also, with 18000 atoms, that kind of speed doesn't sound unreasonable to me. It depends on how many processors you are running on. I don't remeber where I read this, but gromacs can become slow if the order of atoms in the [atoms] section is different. In my gromacs installation, the order of methanol atoms in OPLS-AA is different from the one you write, so maybe you can try using the default one and see if it makes the simulation faster? Sep 20 at 13:25
• Were any of the suggestions in the comments helpful? This gromacs page mentions that certain forcefields (including OPLS/aa) were parameterized with constrained h-bonds and shouldn't be used with all-atom constraints.
– Tyberius
Oct 29 at 18:08

After running some tests, the best way to overcome the slow performance was to simply use the constraints no h-bonds. There was no significant difference in densities and diffusion coefficient when constraint was changed to all-bonds.