Acellera / htmd

HTMD: Programming Environment for Molecular Discovery
https://software.acellera.com/docs/latest/htmd/index.html
Other
255 stars 58 forks source link

Parameterize: treatment of Amber impropers #198

Closed katkovs closed 7 years ago

katkovs commented 7 years ago

Hello HTMD team,

Another question about parameterize in combination with GAFF, this time regarding improper dihedrals. I noticed that in couple of molecules, antechamber assigns certain impropres, while htmd does not. I am attaching the input mol2 file, the output mol2 and frcmod files using antechamber or htmd, and the log file of the parameterize run. Additionally, in the htmd frcmod file there are a lot of parameters which are not used. E.g. in the htmd produced mol2 file there are no za, zb, cf, zc and ce atomtypes, but they are there in the htmd produced frcmod file, including lots of bonded terms involving them. Maybe the impropers are hidden somewhere there, but I can't know this.

So the questions are: 1) Are the impropers not written on puprose? 2) Are these additional terms in the htmd produced frcmod file harmless (or maybe necessary)?

Thanks, Kalina

AFC.zip

j3mdamas commented 7 years ago

About the "additional" terms, I may be mistaken here, and @mj-harvey can probably answer you more accurately, but the frcmod is just an addition to the underlying GAFF, while our frcmod is totally self-contained (all parameters are there).

About torsions, I think you may be right. @mj-harvey, is there any particular reason for not writing the impropers? In the code (https://github.com/Acellera/htmd/blob/de615f12639594e82a38e20f5cb324e0a8e74e04/htmd/parameterization/ff.py#L204-L205), it just prints the header.

mj-harvey commented 7 years ago

If we aren't emitting the impropers terms in the output frcmod that's definitely a bug!

j3mdamas commented 7 years ago

Do you want me to fix that bug?

mj-harvey commented 7 years ago

@j3mdamas Yes, please. @katkovs could you attach an example of a compound with GAFF impropers please?

katkovs commented 7 years ago

The AFC.zip I attached in the first comment contains such an example. Input mol2 file, as well as the outputs from htmd and antechamber are there. The lactate ion I attached in issue #197 should work too, it should have one improper. Let me know if you need anything more.

mj-harvey commented 7 years ago

Super, thanks.

j3mdamas commented 7 years ago

Thanks for reporting. I'll let you know ASAP when it's fixed in master.

João

katkovs commented 7 years ago

Thanks!

j3mdamas commented 7 years ago

@katkovs, this is going to take a little longer. The fact is that impropers aren't supported yet in the frcmod format, and there was a bug on erroring function for that. I've fixed the error on that, but the implementation of impropers on frcmod will have to wait a bit.

Commit: ff034f6a4faa1a6439fb44c7051116551e2b6bf5

João

mj-harvey commented 7 years ago

@katkovs Amber impropers have the same functional form as dihedrals, so we treat them together. I've committed a fix that will correctly separate them when writing the frcmod.

katkovs commented 7 years ago

Perfect! thanks

katkovs commented 7 years ago

Hello HTMD team,

Writing of impropers seems to still be an issue. I am attaching the resulting frcmod, mol2, and prmtop files using antechamber or htmd on the same molecule as before.

AFC2.zip

There seem to be two issues:

  1. Not all imporpers are written out in the prmtop file (compare antechamber.prmtop where there are 9 impropers, vs. htmd.prmtop where there are only 3).
  2. The force constants for the impropers that do get written seem to get lost on prmtop generation with tleap. (In the frcmod file all impropers have force constant of 1.1, while in the prmtop file, the dihedral force constant at the corresponding index 29 is 0.)

I am not familiar with all the details of formatting of Amber files, so maybe I'm missing something. Could these terms be accounted for somewhere else?...

On another note, the parameters for the bonds and angles involving the fluorine atom are slightly different. (Compare bonds f-cf in antechamber.frcmod with zg-zf in htmd.frcmod; angles f-cf-ce in antechamber.frcmod and zg-zf-ze in htmd.frcmod; and angles f-cf-c in antechamber.frcmod and zg-zf-zi in htmd.frcmod.) Are other bonded parameters allowed to change in the dihedral fitting procedure?...

FYI, I am using the most recent htmd version (cloned it from here today).
Let me know if I can provide any more information or files.

Thanks, Kalina

j3mdamas commented 7 years ago

@mj-harvey can you check this?

mj-harvey commented 7 years ago

@katkovs can you double-check that you are using the antechamber, tleap etc that are part of the ambermini package installed as an HTMD dependency.

katkovs commented 7 years ago

Yes, all ambertools-related programs are from the ambermini package installed with conda, version 15.0.4 as far as I see with "conda list". There are no other AmberTools installed on the machine where I ran both htmd and pure antechamber calculations.

mj-harvey commented 7 years ago

@katkovs I can see there are some important errors in frcmod HTMD has written out, I'll investigate and update here later.

mj-harvey commented 7 years ago

@katkovs problem fixed in HEAD. Two problems: 1) Frcmods were being written with the 1-4 ee scaling term inverted 2) fitted dihedrals were being fit with that term == 1.0, rather than Amber's 1/1.2. Not all MM codes can accept multiple values of that term, even when it's speficied per dihedral in the prmtop (Namd 2.10 and earlier, ACEMD 3212 and earlier)

mj-harvey commented 7 years ago

Hi,

Pretty sure that the frcmod is being generated correctly now.

Re the bond and angle terms sometimes being different, we take whatever Amber's parmchk2 emits. That output seems to be geometry sensitive.

katkovs commented 7 years ago

Hello,

thanks for the fix and information.

Regarding the differences in the bonds and angles, yes, you are right, it was a parmchk thing. I used parmchk to create antechamber.frcmod, and not parmchk2 as HTMD does. So, there is no problem here.

Thanks for fixing the issues with the fudge factors. As I am not in depth familiar with Amber formatting, I didn't even notice they were written wrongly.

However, both issues with the writing of the improper dihedrals that I mentioned in my previous comment are still there. Moreover, I noticed that the force constants of some of the proper dihedrals that were not optimized by HTMD have also changed (and not due to different version of parmchk). It almost seems as some of the parameters for the impropers are "leaking" into the propers, in cases where a group of four atom types is found as a proper and improper in the molecule. I'm just speculating though. It also doesn't explain why for example, the improper involving the fluorine atom is not there.

I upload now the new mol.frcmod, strucutre.prmtop, and mol.mol2 created by HTMD. I also upload the output of tleap, since it already gives some warnings about the impropers. Finally, I also upload the HTMD-produced and the "pure antechamber" topologies in Gromacs format, since it is much easier to read and identify differences than the wretched Amber topology files.

AFC3.zip

mj-harvey commented 7 years ago

Ok, I hadn't properly checked the tleap output, I see the issue now. HEAD now writes frcmods that tleap can eat.

katkovs commented 7 years ago

Ok, great! So now only the issue with the missing impropers remains.

j3mdamas commented 7 years ago

We've had this issue with the impropers in the past, but I think it was with charmm. But we don't do any impropers in parameterize, so it's either not being generated by the typer, or it's being lost in the way.

On Feb 10, 2017 2:07 PM, "katkovs" notifications@github.com wrote:

Ok, great! So now only the issue with the missing impropers remains.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/Acellera/htmd/issues/198#issuecomment-278937842, or mute the thread https://github.com/notifications/unsubscribe-auth/AKTzQk9erDe_mPeawYGm4QukJE47wJphks5rbGD_gaJpZM4LACZ2 .

mj-harvey commented 7 years ago

I understand the cause of the problem now. It will take a little while to fix, unfortunately.

katkovs commented 7 years ago

OK, thanks. Just one question: Were the impropers not taken properly into account also in the fitting procedure, or are they simply not written out in the end? I'm asking cause if it's the latter, then I can just simply add them manually in my topology, as a workaround until the bug is fixed.

mj-harvey commented 7 years ago

@katkovs, I cant comment with certainty just yet. Will have an update in a few days.

j3mdamas commented 7 years ago

@mj-harvey, for reference of the discussion we had today, here is the issue with the improper problems on CHARMM: #74

j3mdamas commented 7 years ago

This has to be bumped into 1.8.0.

mj-harvey commented 7 years ago

Should be fixed now in HEAD, @j3mdamas could you confirm please?

j3mdamas commented 7 years ago

It's on my pipeline.

katkovs commented 7 years ago

Hello HTMD team,

Is there any update on this? Or approximate time frame when it is expected to be fixed?

Thanks, Kalina

j3mdamas commented 7 years ago

Hi Kalina,

This issue in particular is tied with other issues in parameterize. A definite solution is going to be available at maximum in 1 month, but probably earlier. If earlier, I'll post here.

Thanks for reporting and for helping us improve the software João

katkovs commented 7 years ago

OK, thanks for the quick reply!

j3mdamas commented 7 years ago

@raimis, this would be a good check of the things we were talking before (validation of parameterize output against gaff output)

raimis commented 7 years ago

Hi @katkovs

We just have released HTMD 1.9.5, which fixes several improper issues.

I tried with your structure:

parameterize -m ACF3/mol.mol2 --no-min --no-esp --no-torsions

mol.frcmod file contains 5 improper, as expected from the structure:

...
IMPR
cf-za-zb-za     1.100000 180.000000 2.000000
zb-ce-cf-zc     1.100000 180.000000 2.000000
ca-cf-ce-ha     1.100000 180.000000 2.000000
ca-ca-ca-ce     1.100000 180.000000 2.000000
ca-ca-ca-ha     1.100000 180.000000 2.000000
...

And in structure.prmtop their are applied 9 times.

Could you verify it solves the issue, please.

raimis commented 7 years ago

@katkovs?

katkovs commented 7 years ago

Hello htmd team,

Sorry for the late reply.

Unfortunately, the project i needed the parameters for is on hold, and I won't be able to confirm if the problem was solved.

As far as I'm concerned, you may close this issue.

Thank you for the efforts.

Best regards Kalina

On Aug 21, 2017 10:57 AM, "Raimondas Galvelis" notifications@github.com wrote:

@katkovs https://github.com/katkovs?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/Acellera/htmd/issues/198#issuecomment-323686848, or mute the thread https://github.com/notifications/unsubscribe-auth/AWaV-7K21Mz0vUy_h17iDFBqKPhDuJWGks5saUaBgaJpZM4LACZ2 .

raimis commented 7 years ago

@j3mdamas, let's close this!