{"action":"created","issue":{"url":"https://api.github.com/repos/choderalab/openmm/issues/16","repository_url":"https://api.github.com/repos/choderalab/openmm","labels_url":"https://api.github.com/repos/choderalab/openmm/issues/16/labels{/name}","comments_url":"https://api.github.com/repos/choderalab/openmm/issues/16/comments","events_url":"https://api.github.com/repos/choderalab/openmm/issues/16/events","html_url":"https://github.com/choderalab/openmm/issues/16","id":144349641,"number":16,"title":"Proposal for new improper-matching algorithm","user":{"login":"rafwiewiora","id":12088757,"avatar_url":"https://avatars.githubusercontent.com/u/12088757?v=3","gravatar_id":"","url":"https://api.github.com/users/rafwiewiora","html_url":"https://github.com/rafwiewiora","followers_url":"https://api.github.com/users/rafwiewiora/followers","following_url":"https://api.github.com/users/rafwiewiora/following{/other_user}","gists_url":"https://api.github.com/users/rafwiewiora/gists{/gist_id}","starred_url":"https://api.github.com/users/rafwiewiora/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/rafwiewiora/subscriptions","organizations_url":"https://api.github.com/users/rafwiewiora/orgs","repos_url":"https://api.github.com/users/rafwiewiora/repos","events_url":"https://api.github.com/users/rafwiewiora/events{/privacy}","received_events_url":"https://api.github.com/users/rafwiewiora/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":28,"created_at":"2016-03-29T18:34:35Z","updated_at":"2016-03-30T19:41:39Z","closed_at":null,"body":"@swails @jchodera @peastman \r\n\r\nI have looked into the problems with the impropers and come up with a revised version of the matching code. Have a glance at the existing code: https://github.com/pandegroup/openmm/blob/master/wrappers/python/simtk/openmm/app/forcefield.py#L1360-L1390\r\n\r\nWe had:\r\n* first has to match (central atom) -- this becomes the 3rd atom in the torsion added to System\r\n* a permutation of the other 3 then matches the order as written in the ffxml --> ParmEd ensures atom types are alphabetically-ordered --> I observed the permutations here to lead to problems when atom 1 and 4 or 2 and 4 had the same type. The code as is takes the 4th atom through as whatever was matched in the first matching permutation here - while it appears tLeap applies sorting by the index of the atom in such a situation\r\n* then we had the bizarre stuff - problems here were tLeap appears to ensure that if you have a no-wildcard improper torsion, you should match the exact order as written in the parameters you loaded (assuming the parameter file had those alphabetized, so it's rather the alphabetically ordered torsion has to be adhered to if no wildcards. The code here didn't make that distinction and would do sorting only based on the elements of the atoms involved. Now I confirmed the rule on line 1385 applies to wildcard containing impropers in tLeap but I failed to find any justification for line 1387. Perhaps authors if present here could explain? \r\n\r\nAll these observations I made by running a piece of dev code I wrote. It loads in a leaprc (I was using leaprc.ff14SB) to Parmed, turns off all propers and all impropers but one (by setting all those other phi_k 's to zero) and writes out a .dat with these new parameters. These are then loaded back in using a leaprc hooked up to the modified .dat, converted to ffxml and we do energy comparisons on ala3/villin/DNA/RNA and log out the energies, the improper as written in the .dat, and the sequences as assigned from the prmtop and from ffxml-System with element, type and index information. Some rules were immediately evident:\r\n* treat torsions using wild cards and not using wild cards separately\r\n* if no wildcards - ordering of atom types must match the ordering in the ffxml (but first atom becomes 3rd atom of course), however if there are >1 same atom types within such matched torsion - order the identical atom types by indexes of the atoms\r\n* if wildcards - atoms of the same elements are ordered by indexes of the atoms - this applies to 1-4 and 2-4, atoms 1 and 2 are ordered by index regardless of what they are - at least one of these will be from a wildcard, so makes sense to treat 4 differently to 1/2\r\n\r\nThis gives the new proposal code:\r\n``` python\r\nif type1 in types1:\r\n    for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):\r\n        if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:\r\n            a1 = torsion[t2[1]]\r\n            a2 = torsion[t3[1]]\r\n            a4 = torsion[t4[1]]\r\n            e1 = data.atoms[a1].element\r\n            e2 = data.atoms[a2].element\r\n            e4 = data.atoms[a4].element\r\n            if not hasWildcard:\r\n                if t2[0] == t4[0] and a1 > a4:\r\n                    (a1, a4) = (a4, a1)\r\n                if t3[0] == t4[0] and a2 > a4:\r\n                    (a2, a4) = (a4, a2)\r\n                if t2[0] == t3[0] and a1 > a2:\r\n                    (a1, a2) = (a2, a1)\r\n            elif hasWildcard:\r\n                if e1 == e4 and a1 > a4:\r\n                    (a1, a4) = (a4, a1)\r\n                if e2 == e4 and a2 > a4:\r\n                    (a2, a4) = (a4, a2)\r\n                if a1 > a2:\r\n                    (a1, a2) = (a2, a1)    \r\n            match = (a1, a2, torsion[0], a4, tordef)\r\n            break\r\n```\r\n\r\nTwo issues here remain:\r\n* this appeared to work perfectly for leaprc.ff14SB in my dev code, but when I went back to not iterating over every improper but checking the energies of everything at once - one improper wouldn't match anymore. (I also wrote a new validation function that checks the atom index sequences for all impropers AMBER vs OpenMM). This is because there is an RNA torsion that has a `C5` atom type and a `CB`, and `CB` comes before `C5`in the `parm10.dat` but the alphabetical ordering in Leap doesn't seem to correct it. (Leap just uses a `strcmp` so C5 should come before CB, no idea...). In my dev code that went through Parmed first, which enforced the right alphabetization (and so it does for ffxmls)\r\n\r\n* some impropers for some other systems do not match also\r\n\r\nBUT\r\n\r\nthe relative errors have decreased a lot. All are now below the tolerance for all other forces of `1e-5` except in ffAM1, ffPM3, ff96 and one residue in phospho-ff + 14SB. You can see this by looking at the last commit diff (compare second red line change with second green line): https://github.com/choderalab/openmm/pull/15/commits/cc63dd80d83091b16bebd5a59451fb6db880d3b8\r\n\r\nThe ones that are above the 1e-5 are below the 1e-2 tolerance - the highest is around 3e-3. So this is a big improvement. I've spent since the morning in LEaP code trying to justify this with their code and my patience has run out for now - perhaps I'll try again later. But from your previous discussions (https://github.com/pandegroup/openmm/issues/220 - still open btw, perhaps good to close now as this will overtake) it sounds like this is not even internally consistent in AMBER - in this case I suggest settling on:\r\n* there will be some relative error but most important thing is these torsions keep things flat\r\n* choose matching scheme based on what gives the lowest relative error in energies - my issue here is I came up with this new code just by looking at the results of LeAP not through understanding their code. Hence maybe I introduced some big issue that will come up only in some systems I don't test here and the older code would have dealt better? On the other hand the test systems here are fairly comprehensive so should cover most stuff. What I would like to do is incorporate some code that will loop over all residues in an ffxml and build systems in LeAP for those to compare energies using absolutely all parameters. I'll work on that when I deal with other priorities - that would be the ultimate validation. I worry a bit about removing this mass rule that I said I saw no evidence for: https://github.com/pandegroup/openmm/blob/master/wrappers/python/simtk/openmm/app/forcefield.py#L1387 (anyone remember the origin of this?)\r\n\r\nAlso I can output an 'improper log' that lists all differences between AMBER and OpenMM assignments found in testing (in atom indexes, this is using the new algorithm):\r\nhttps://gist.github.com/rafwiewiora/22bcaa5beb517642957b75ea07e82bf1\r\n\r\nWhat do you think?"},"comment":{"url":"https://api.github.com/repos/choderalab/openmm/issues/comments/203602429","html_url":"https://github.com/choderalab/openmm/issues/16#issuecomment-203602429","issue_url":"https://api.github.com/repos/choderalab/openmm/issues/16","id":203602429,"user":{"login":"swails","id":779022,"avatar_url":"https://avatars.githubusercontent.com/u/779022?v=3","gravatar_id":"","url":"https://api.github.com/users/swails","html_url":"https://github.com/swails","followers_url":"https://api.github.com/users/swails/followers","following_url":"https://api.github.com/users/swails/following{/other_user}","gists_url":"https://api.github.com/users/swails/gists{/gist_id}","starred_url":"https://api.github.com/users/swails/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/swails/subscriptions","organizations_url":"https://api.github.com/users/swails/orgs","repos_url":"https://api.github.com/users/swails/repos","events_url":"https://api.github.com/users/swails/events{/privacy}","received_events_url":"https://api.github.com/users/swails/received_events","type":"User","site_admin":false},"created_at":"2016-03-30T19:41:39Z","updated_at":"2016-03-30T19:41:39Z","body":"I personally don't think the impropers are a big deal -- they're there to get planarity correct and avoid bad things.  The force constants are not carefully picked (if they're too weak, make them stronger), and many other programs assign them differently than LEaP and claim to implement the same FF (see GROMACS).\r\n\r\nHowever, there are a sizeable number of people who will dismiss any calculation that implements any term (including impropers) differently than Amber because, well, *we just don't know*.  I'm of the opinion that the approximation made in improper ordering is easily lost in the noise of the FF as a whole.  But all the reviewer has to say is \"prove it\" and I have no solid reply until either a) I do a comprehensive study with a wide array of molecules to demonstrate (and even that only weakens an outright dismissal -- it doesn't invalidate it) or b) do a comprehensive comparison (out to convergence) with both approaches for my system and prove the answers are the same.\r\n\r\nSo I think the most sensible thing to do is to train OpenMM to assign periodic impropers the way Amber does and fix the ambiguity in Amber that results in the aforementioned inconsistencies.  As a note, the only time I'm aware of ambiguous ordering coming out of tleap is when 2+ atom types have the same name in the improper (and neither is the central atom).  Every other time the improper is well-defined.  But I don't know the exact rules well enough."}}