7
$\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
  • 2
    $\begingroup$ Did you try changing constraints to h-bonds? $\endgroup$ Sep 19 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
    Sep 19 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
    Sep 19 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
    Sep 20 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
    Oct 29 at 18:08
3
$\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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.