Dear Carlos,
In the Gromacs topologies written by ParmEd, most parameters are
assigned by atom number rather than type - at least in my test cases.
The only exception are the VdW parameters, but they are easy to copy
along with the atom type.
This makes it much easier to change the topologies for CMAP, since
renaming the atom types does not mess with other parameters. After some
testing, I came up with a short Python script, which I will paste at the
end of this E-mail.
I ran a very short simulation of BPTI in OPC water, and I was able to
reproduce the CMAP energy very accurately using gmx -rerun (MAE = 0.16
kcal/mol). The other parameters are not entirely accurate, but I believe
that this is due to coordinate precision rather than force field
problems. E.g., the bond terms have a MAE of 2.3 kcal/mol, in a system
of 15000 particles (I ran my tests using ntf=1 for easier comparison).
I agree that we should get Jason's opinion on this. I will post the same
information on the Github issue mentioned before, to see what he has to
say about it.
Thanks for your help.
Best,
Franz
(ff19_to_gromacs.py)
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import parmed as pmd
parm = pmd.load_file("topology.parm7", "coord.rst7")
known_types = {}
used_names = set(a.atom_type.name for a in parm.atoms)
for cmap in parm.cmaps:
if cmap.type not in known_types:
atom_type = cmap.atom3.atom_type
i = 0
while(f"XC{i}") in used_names:
i += 1
new_typename = f"XC{i}"
used_names.add(new_typename)
new_type = pmd.AtomType(new_typename, atom_type.number,
atom_type.mass, atom_type.atomic_number, new_typename, atom_type.charge)
new_type.epsilon = atom_type.epsilon
new_type.rmin = atom_type.rmin
known_types[cmap.type] = new_type
cmap.atom3.atom_type = known_types[cmap.type]
cmap.atom3.type = known_types[cmap.type].name
parm.save("converted.top") # .top saves to GROMACS topology file format
parm.save("converted.gro")
On 21.06.23 18:07, Carlos Simmerling wrote:
> I think it's helpful to consider separately the way the cmaps and
> other parameters are stored in the library files used as input for
> building the topology with leap, vs the way they are stored in the
> final topology. The library files are more general, and mostly based
> on atom types (except the cmaps). inside the topology file, however,
> things are mostly atom number based, with pointers into lists of
> parameter sets for bonds, angles and so on, including the cmaps. the
> way the cmaps are stored in the topology is very similar to the way
> the bonds are stored - a set of atom numbers, followed by a pointer
> into the list of cmap parameters.
>
> %FLAG CMAP_INDEX
> %FORMAT(6I8)
> 5 7 9 22 24 9
>
>
> the way I see it, the different options for adding multiple cmaps
> would have been to (1) use residue names, or (2) as Jason suggested,
> creating different atom types for each amino acid. the reason that we
> didn't use atom types is that it would have involved duplicating (in
> the parameter library files) all of the bond/angle/dihedral terms
> involving the Calpha atoms for each of the different amino acids, as
> well as their vdw info. If any changes to these parameters are needed
> in the future, each copy would need to be updated, which can be
> error-prone.
>
> Once the topology file is built, I am not sure if the atom types are
> used other than for user selection of atoms by their atom type. All of
> the parameters are based on atom numbers at that point. (see
> https://ambermd.org/FileFormats.php#topology). At this point, one can
> change the atom types at will and it should not (I think) affect any
> of the Amber energy calculations - since the force field parameters
> are already encoded by atom numbers, there is no need to make
> duplicates of any of the force field terms. So it worked for you to
> change atom types after the topology is built, but it would be far
> more cumbersome to split the Calpha atom types before the topology is
> built. However, I believe that current Amber allows only 2 characters
> for each atom type, so your example can't be done inside the Amber
> topology file.
>
> It seems to me that the cmap conversion into a gromacs
> topology should proceed similarly to how all of the other terms are
> handled. But I don't know how the parmed conversion works - for
> example, does it convert it so gromacs can select bond parameters
> based on atom type as well, or does it use the atom numbers that are
> used in the Amber topology? Are cmaps the only thing that uses atom
> types instead of numbers? If gromacs can only used cmaps that are
> selected by atom types, then a script or parmed modification to change
> the atom types is definitely the easiest approach, though I honestly
> don't see why cmaps need atom types defined but the other force field
> terms do not. In any case, one can easily (probably with a script)
> look at the atom numbers in the cmap list, and for each cmap type,
> alter the atom type for a specific atom number before converting to
> the gromacs format. Then gromacs should be able to know which cmap to
> use for each residue. This is essentially what you did by hand, but
> should not be difficult inside parmed since the info is all there.
> This is far easier than starting out with major duplication of the
> non-cmap parameters in the Amber force field library files.
>
> It could also help to involve Jason in these discussions since I may
> have missed something important here.
>
> carlos
>
> On Wed, Jun 21, 2023 at 11:25 AM Franz Waibl via AMBER
> <amber.ambermd.org> wrote:
>
> Hi,
>
> I am also no expert on Gromacs topologies, but it seems like the [
> cmap
> ] entry only defines which atoms should have CMAP terms at all, while
> the [ cmaptypes ] entry defines a mapping from atom types to CMAP
> type.
> For example, one definition in my .top file looks like this:
>
> [ cmaptypes ]
>
> C N XC C N 1 24 24\
> (followed by 576 numbers)
>
> To adapt my example topology, I copied the XC atom type to
> generate XC1
> and XC2, and added the corresponding entries. I did not need any
> further
> changes.
>
> Best,
> Franz
>
> On 21.06.23 16:52, Carlos Simmerling wrote:
> > Hi Franz,
> > thanks that's very helpful. Amber uses the residue name inside Leap
> > while it is building the topology, but inside the
> final topology file
> > itself the cmaps are specified using atom number and cmap id
> (pointing
> > to a table of cmaps that are stored in the topology file
> itself), so
> > no residue name is needed. I'm not sure why Jason suggested that
> the
> > prmtop uses the residue name (unless I misunderstood his comment on
> > the github issue). It sounds like converting what is in the Amber
> > prmtop to something needed for gromacs should be
> straightforward. As I
> > said I'm not familiar myself with how gromacs implements cmaps
> or how
> > they are stored. Does it assign cmaps based on atom type?
> > carlos
> >
> > On Wed, Jun 21, 2023 at 10:29 AM Franz Waibl via AMBER
> > <amber.ambermd.org> wrote:
> >
> > Dear Carlos, dear Amber users,
> >
> > thanks for the reply. I am aware that Amber is free for
> academic use
> > (thanks for that, by the way!), but it would still be useful
> to be
> > able
> > to run ff19SB simulations on Gromacs, e.g., to use
> algorithms that
> > are
> > not implemented in Amber.
> >
> > Gromacs does support multiple CMAPs. If I understand
> correctly, the
> > problem is that ff19SB CMAP terms depend on the residue name in
> > addition
> > to the atom type. To convert the topology correctly, one would
> > have to
> > define separate atom types for the CA atom of each amino acid
> > (XC1, XC2,
> > ...). This also seems to be what CHARMM-GUI is doing (see
> the Github
> > issue mentioned before).
> > In fact, I ran a very short test simulation in Amber and
> compared the
> > results to energy values from gmx mdrun -rerun. After manually
> > tweaking
> > the topology as described above, I was able to obtain
> consistent CMAP
> > energies between Amber and Gromacs.
> >
> > However, it is very tedious to do this manually for a larger
> > system, so
> > a ParmEd-based solution would be highly appreciated.
> >
> > Best regards,
> > Franz
> >
> >
> > On 21.06.23 15:14, Carlos Simmerling wrote:
> > > My understanding is that the issue is on the gromacs end, not
> > parmed.
> > > I'm not a gromacs user but you could ask them if multiple
> cmaps are
> > > supported. If yes, then we can work to update the conversion
> > process.
> > > Alternatively, Amber is free and has a good library of
> tutorials.
> > > Carlos
> > >
> > > On Wed, Jun 21, 2023, 8:59 AM Franz Waibl via AMBER
> > > <amber.ambermd.org> wrote:
> > >
> > > Dear Amber users,
> > >
> > > I would like to run simulations with the Amber-ff19SB
> > forcefield on
> > > Gromacs. I have seen an Issue on Github
> > > (https://github.com/ParmEd/ParmEd/issues/1292) stating
> that
> > ParmEd
> > > does
> > > not handle the CMAP term conversion correctly. I also ran
> > some tests
> > > with the ParmEd version in Ambertools23 installed via
> conda,
> > and it
> > > seems that the same problem occurs: Only one set of CMAP
> > > parameters is
> > > written to the Gromacs topology, rather than different
> > parameters for
> > > each amino acid.
> > >
> > > Do you know if this problem can be avoided somehow in
> > ParmEd? If no,
> > > what is the recommended way to obtain ff19SB
> topologies for
> > Gromacs?
> > >
> > > Thanks for your help.
> > >
> > > Best regards,
> > > Franz
> > >
> > > --
> > > ETH Zürich
> > > Franz Waibl
> > > Lab. für Physikalische Chemie
> > > HCI G 227
> > > Vladimir-Prelog-Weg 1-5/10
> > > 8093 Zürich
> > >
> > > Telefon +41 44 632 55 04
> > > franz.waibl.phys.chem.ethz.ch
> > >
> > >
> > > _______________________________________________
> > > AMBER mailing list
> > > AMBER.ambermd.org
> > > http://lists.ambermd.org/mailman/listinfo/amber
> > >
> >
> > --
> > ETH Zürich
> > Franz Waibl
> > Lab. für Physikalische Chemie
> > HCI G 227
> > Vladimir-Prelog-Weg 1-5/10
> > 8093 Zürich
> >
> > Telefon +41 44 632 55 04
> > franz.waibl.phys.chem.ethz.ch
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> >
>
> --
> ETH Zürich
> Franz Waibl
> Lab. für Physikalische Chemie
> HCI G 227
> Vladimir-Prelog-Weg 1-5/10
> 8093 Zürich
>
> Telefon +41 44 632 55 04
> franz.waibl.phys.chem.ethz.ch
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
ETH Zürich
Franz Waibl
Lab. für Physikalische Chemie
HCI G 227
Vladimir-Prelog-Weg 1-5/10
8093 Zürich
Telefon +41 44 632 55 04
franz.waibl.phys.chem.ethz.ch
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jun 22 2023 - 01:00:03 PDT
