8
$\begingroup$

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.

$\endgroup$
7
  • 3
    $\begingroup$ Did you try changing constraints to h-bonds? $\endgroup$ Commented Sep 19, 2021 at 6:27
  • $\begingroup$ 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.) $\endgroup$
    – dwhswenson
    Commented Sep 19, 2021 at 7:19
  • $\begingroup$ yes, I will be trying changing to h-bonds. But my question is, what if i want to have all-bonds to be constrained? $\endgroup$
    – megamence
    Commented Sep 19, 2021 at 23:41
  • 2
    $\begingroup$ ...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? $\endgroup$
    – S R Maiti
    Commented Sep 20, 2021 at 13:25
  • 1
    $\begingroup$ 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. $\endgroup$
    – Tyberius
    Commented Oct 29, 2021 at 18:08

1 Answer 1

4
$\begingroup$

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.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .