On Thu, Nov 18, 2021 at 7:23 AM ABEL Stephane <Stephane.ABEL.cea.fr> wrote:
> Dear all,
>
> Below a little update of of my testing of AMBER19SB with GROMACS(2021.3).
> I have converted prmto/impcrd files (equivalent to prm7 and rst7) of a
> protein modeled with the AMBER19SB into GROMACS gro/top files using Parmed
> (AT21, fully patched) and ACEPYPE (last version).
> The AMBER files were obtained using tleap (AT21). From the results below
> with the same gro file
>
> I used the ParmEd commands below from Jason Swails
>
> import parmed as pmd
> parm = pmd.load_file("MYTOPFF19SBparm7", "MYTOPFF19SB.rst7")
> parm.save("MYTOPFF19SB.top") # .top saves to GROMACS topology file format
> parm.save("MYTOPFF19SB.gro")
> ```
> And execute it as `python <file.py>`. Does that work for you?
>
> The results with the same mdp of GROMACS at step0
>
> From Parmed (Parmed_Conversion_GROMACS_AMBER19SB.py, VERSION 3.4.1) at
> step0
> Step Time
> 0 0.00000
>
> Energies (kJ/mol)
> Bond Angle Proper Dih. Improper Dih. CMAP Dih.
> 1.52227e+03 7.26001e+03 1.15737e+04 2.23252e+02 3.32921e+03
> LJ-14 Coulomb-14 LJ (SR) Coulomb (SR) Potential
> 1.02351e+04 9.98003e+04 -1.76602e+04 -1.87376e+05 -7.10925e+04
> Pressure (bar)
> 1.77444e-01
>
>
[snip]
> From Sander with coordinates file but not the same
>
> NSTEP = 0 TIME(PS) = 0.000 TEMP(K) = 0.00 PRESS =
> 0.0
> Etot = -15720.7550 EKtot = 0.0000 EPtot =
> -15720.7550
> BOND = 348.4593 ANGLE = 1731.7886 DIHED =
> 2819.0630
> 1-4 NB = 2449.9439 1-4 EEL = 23852.5075 VDWAALS =
> -4733.5065
> EELEC = -42775.8491 EHBOND = 0.0000 RESTRAINT =
> 0.0000
> CMAP = 586.8384
>
>
> SANDER* (AT21) Kcal/mol** (KJ/mol)*** with 1kjoul/mol = 4.1867
> Kcal/mol
> **
> step 0
> BOND = 348.4593* = 1458.89** ---> Similar to GROMACS
> results above
> ANGLE = 1731.7886 = 7250.479** ---> Similar to GROMACS results
> above
> DIHED = 2819.0630 = 11802.57** Similar to GROMACS results above
> (Proper Dih + Improper Dih)
> CMAP = 586.8384 = 3.32921e+03 ---> Different from
> GROMACS obtained with ParmEd (3.32921e+03)
>
Am I missing something? The CMAP terms look identical?
* | mdin Single point
> ; see
> http://ringo.ams.stonybrook.edu/index.php/2012_AMBER_Tutorial_with_Biotin_and_Streptavidin
> &cntrl
> <http://ringo.ams.stonybrook.edu/index.php/2012_AMBER_Tutorial_with_Biotin_and_Streptavidin&cntrl>
> imin=0,
> ntpr=1,
> maxcyc=100,
> ntmin=2,
> ntb=0,
> igb=0,
> cut=999
> /
>
> on can say
>
> 1) The converted file generated by ParmEd is well recognized by GROMACS
> and in particular the CMAP term. IT is not the case from the file obtained
> bt ACEPYPE (since it does not support CMAP for AMBER yet)
> 2) CMAP energy in AMBER and GROMACS are different
Unless I'm misreading your email, the CMAP energies between AMBER and
GROMACS look the same.
> 2) the calculations of the energy for the other bonded and non bonded are
> the same (with rounding error)
> 2) WARNING !!! with the same inpcrd Parmed does not genertae a corrct GRO
> file. The box is to too small
> 3.24900 3.12110 3.31620 ---> end gro file from ParmEd
> 117.23400 182.43200 135.66400 ---> end gro file from ACEPYPE
>
For box dimensions, the acpype numbers you post look absurd -- that's a
1172x1824x1356 angstrom box (compared to the more reasonable 32x31x33
angstrom box that is more prototypical of what tleap will generate when
solvating a system).
3) Note that at this moment ParmEd seems not support mixing 1-4 LJ if the
> prmtop contains atom types from GRLYCAM06 and AMBER19SB ffields. The
> program crashes with the following error
> Traceback (most recent call last):
> File "Parmed_Conversion_GROMACS_AMBER19SB.py", line 3, in <module>
> parm.save("ShuA_BOG_AMBER19SB_GLYCAM06_HMASS_From_Parmed.top") # .top
> saves to GROMACS topology file format
> File
> "/home/stephane_abel/amber20/lib/python3.9/site-packages/parmed/structure.py",
> line 1489, in save
> s = gromacs.GromacsTopologyFile.from_structure(self)
> File
> "/home/stephane_abel/amber20/lib/python3.9/site-packages/parmed/gromacs/gromacstop.py",
> line 1279, in from_structure
> raise GromacsError('Structure has mixed 1-4 scaling which is '
> parmed.exceptions.GromacsError: Structure has mixed 1-4 scaling which is
> not supported by Gromacs
>
This is not a ParmEd limitation - it's a GROMACS one. GROMACS doesn't let
you specify pair-specific 1-4 electrostatic terms (it does for
Lennard-Jones, but not electrostatic). As a result, it's fundamentally
impossible to represent a mixed-scaling system with a GROMACS topology
file. Rather than generate something that is "close" (but still wrong in
definition), ParmEd throws an error. I would urge strong caution if acpype
generates a GROMACS topology file from a mixed force field like this, as
the 1-4 terms are unlikely to match exactly.
If the sugar or non-sugar components to the system are small, then the
difference in 1-4 electrostatic contributions may be small enough to escape
detection (i.e., be confused with round-off/precision errors), but in this
case the treatment is actually wrong (which would be easier to observe by
looking at the computed forces, most likely).
All the best,
Jason
--
Jason M. Swails
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Nov 22 2021 - 06:00:02 PST
