input files for testing the energy conservation between amber and gromacs

[xiki-tempula]

I'm interested in converting the amber topology to gromacs topology. Just for a sanity test, I kind of wants to make sure that the converted topology is correct. My plan is to use original amber topology to compute the potential of the system in amber and compare it with the energy generated with gromacs. However, due to the difference in amber and gromacs handle things, I cannot get the energy to agree with each other. I wodner if you can suggest an amber input file and gromacs mdp file which can be used for this purpose? Thank you.

[mrshirts]

Can you print out the energy output you got with InterMol (which, if you configure right, will do the calculation for you with the correctly corresponding input files) There are some things that you can and can't get to agree no matter how you tweak the input files: see https://link.springer.com/article/10.1007/s10822-016-9977-1 for a discussion of these issues.

[xiki-tempula]

I am generating the forcefield for a lipid AR_PS_LAL The amber file is generated with tleap.

tleap.in

source leaprc.protein.ff14SB
source leaprc.lipid17
lipid = sequence {AR PS LAL}
setBox lipid "vdw" 20
saveamberparm lipid AR_PS_LAL.prmtop AR_PS_LAL.inpcrd

Then I run tleap tleap -f tleap.in > tleap.out After that, I run the intermol python convert.py --amb_in AR_PS_LAL.prmtop AR_PS_LAL.inpcrd --gromacs -e

I got

INFO 2019-11-13 16:23:08 Beginning InterMol conversion

INFO 2019-11-13 16:23:08 InterMol is research software. If you make use of InterMol in scientific publications please cite the following reference:
Shirts, M.R., Klein, C., Swails, J.M. et al. J Comput Aided Mol Des (2016). doi:10.1007/s10822-016-9977-1

AR_PS_LAL/AR_PS_LAL.prmtop AR_PS_LAL.inpcrd
INFO 2019-11-13 16:23:08 Evaluating energy of AR_PS_LAL.inpcrd
INFO 2019-11-13 16:23:10 InterMol Conversion Energy Comparison Results

amber input energy file: amber.out
gromacs output energy file: energy.xvg

Energy group summary
=======================================================================
                type       input(amber)  output (gromacs)    diff (gromacs)
-----------------------------------------------------------------------
Not comparable energies: are likely not to be the same
-----------------------------------------------------------------------
        coulomb (LR)         0.00000000      178.46041870      178.46041870
        coulomb (SR)         0.00000000    -1590.01171875    -1590.01171875
       coulomb total      -633.27810640     -633.39700317       -0.11889677
          coulomb-14       778.12776800      778.15429688        0.02652887
            improper         0.00000000        0.00949546        0.00949546
              proper       221.53652400      221.52671719       -0.00980681
            vdw (LR)         0.00000000       -0.12749512       -0.12749512
            vdw (SR)         0.00000000  2591723.50000000  2591723.50000000
           vdw total   2591816.82461040  2591776.12575118      -40.69885922
              vdw-14        52.75312720       52.75324631        0.00011911
-----------------------------------------------------------------------
Comparable energy terms: these should be very close
-----------------------------------------------------------------------
               angle        83.39465120       83.39507294        0.00042174
                bond        44.30270240       44.30338669        0.00068429
              bonded       349.23387760      349.23467228        0.00079468
            dihedral       221.53652400      221.53621265       -0.00031135
           nonbonded   2591183.54650400  2591142.72874801      -40.81775599
           potential   2591532.78038160  2591492.00000000      -40.78038160

---------------- Total Potential Energy Comparison --------------------
Input amber potential energy:       2591532.78038160
Difference in potential energy from amber=>gromacs conversion:       -40.78038160
=======================================================================
INFO 2019-11-13 16:23:10 Finished!

I felt that the energy difference is larger than I would like it to be.

[mrshirts]

That's essentially exactly what you should get. Look over the paper. The nonbondeds are never going to be better than 1 part in 10^5 or so, which is what they are here. There are differences in how different codes evaluate nonbondeds that are not insignificant. Most of these differences will cancel out in thermodynamic properties, though - some differences will not cancel out, but that's a difference between programs, not force fields.

Note that coul 1-4 and vdw 1-4 agree quite well (even though they are on the upper list), which means parameters are being assigned correctly.

[xiki-tempula]

Thank you. Though I'm still a bit worried about the 40 kj/mol difference I agree they might cancel each other out during the actual simulation.

[mrshirts]

Perhaps it's a problem, but it's the problem of the programs, not the new force field parameters!

On Wed, Nov 13, 2019 at 11:50 AM Zhiyi Wu notifications@github.com wrote:

Thank you. Though I'm still a bit worried about the 40 kj/mol difference I agree they might cancel each other out during the actual simulation.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/shirtsgroup/InterMol/issues/361?email_source=notifications&email_token=ABATPVFVTZL7K4JIA3PCA5DQTQV5BA5CNFSM4JM4FFK2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOED6ZMVA#issuecomment-553490004, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABATPVA6LUV64WXOGEPHEITQTQV5BANCNFSM4JM4FFKQ .