{"action":"closed","issue":{"url":"https://api.github.com/repos/choderalab/yank/issues/585","repository_url":"https://api.github.com/repos/choderalab/yank","labels_url":"https://api.github.com/repos/choderalab/yank/issues/585/labels{/name}","comments_url":"https://api.github.com/repos/choderalab/yank/issues/585/comments","events_url":"https://api.github.com/repos/choderalab/yank/issues/585/events","html_url":"https://github.com/choderalab/yank/issues/585","id":191809060,"number":585,"title":"Fixes for heterogeneous dispersion correction ","user":{"login":"jchodera","id":3656088,"avatar_url":"https://avatars.githubusercontent.com/u/3656088?v=3","gravatar_id":"","url":"https://api.github.com/users/jchodera","html_url":"https://github.com/jchodera","followers_url":"https://api.github.com/users/jchodera/followers","following_url":"https://api.github.com/users/jchodera/following{/other_user}","gists_url":"https://api.github.com/users/jchodera/gists{/gist_id}","starred_url":"https://api.github.com/users/jchodera/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/jchodera/subscriptions","organizations_url":"https://api.github.com/users/jchodera/orgs","repos_url":"https://api.github.com/users/jchodera/repos","events_url":"https://api.github.com/users/jchodera/events{/privacy}","received_events_url":"https://api.github.com/users/jchodera/received_events","type":"User","site_admin":false},"labels":[{"id":62138045,"url":"https://api.github.com/repos/choderalab/yank/labels/enhancement","name":"enhancement","color":"84b6eb","default":true}],"state":"closed","locked":false,"assignee":null,"assignees":[],"milestone":{"url":"https://api.github.com/repos/choderalab/yank/milestones/1","html_url":"https://github.com/choderalab/yank/milestones/yank%201.0","labels_url":"https://api.github.com/repos/choderalab/yank/milestones/1/labels","id":477995,"number":1,"title":"yank 1.0","description":null,"creator":{"login":"jchodera","id":3656088,"avatar_url":"https://avatars.githubusercontent.com/u/3656088?v=3","gravatar_id":"","url":"https://api.github.com/users/jchodera","html_url":"https://github.com/jchodera","followers_url":"https://api.github.com/users/jchodera/followers","following_url":"https://api.github.com/users/jchodera/following{/other_user}","gists_url":"https://api.github.com/users/jchodera/gists{/gist_id}","starred_url":"https://api.github.com/users/jchodera/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/jchodera/subscriptions","organizations_url":"https://api.github.com/users/jchodera/orgs","repos_url":"https://api.github.com/users/jchodera/repos","events_url":"https://api.github.com/users/jchodera/events{/privacy}","received_events_url":"https://api.github.com/users/jchodera/received_events","type":"User","site_admin":false},"open_issues":33,"closed_issues":104,"state":"open","created_at":"2013-11-07T21:09:29Z","updated_at":"2016-11-30T01:56:50Z","due_on":null,"closed_at":null},"comments":3,"created_at":"2016-11-26T14:39:38Z","updated_at":"2016-11-30T01:56:50Z","closed_at":"2016-11-30T01:56:50Z","body":"We've noticed some surprises in implementing the heterogeneous dispersion correction scheme from [our paper with @mrshirts and @davidlmobley](http://www.choderalab.org/publications/2014/4/26/accurate-and-efficient-corrections-for-missing-dispersion-interactions-in-molecular-simulations).\r\n\r\nIn validating YANK on binding free energy calculations, we've been working to reproduce some of the [T4 lysozyme L99A results from @davidlmobley](http://www.choderalab.org/publications/2014/4/26/predicting-absolute-ligand-binding-free-energies-to-a-simple-model-site). We perform the simulation with PME and an isotropic dispersion correction at 300 K and 1 atm, using a cutoff of 10A. For the complex leg of the alchemical cycle, we enlarge the cutoff to 16A at both lambda=1 (fully interacting ligand) and lambda=0 (noninteracting ligand), and estimate the free energy for this leg between the enlarged endpoints using MBAR.\r\n\r\nSurprisingly, this method seems to add a great deal to the statistical error, which seems to differ from the findings of our earlier work. Enlarging the cutoffs at the endpoints for ethylbenzene binding to T4 lysozyme adds ~1.2 kT to the statistical error in a 1.5 ns/replica simulation, and dominates the statistical error in the computed binding free energy. \r\n\r\nTo study this problem in greater detail, I ran a replica-exchange simulation of Src kinase in explicit solvent, using different nonbonded cutoffs: `[9, 10, 11, ..., 24] A` (which set both the real-space and LJ cutoffs), both with and without Lennard-Jones switching. You can see the statistical error from a 100 ps simulation by reading across from the corresponding row. For example, reweighting from 9A to any other cutoff gives ~0.2 kT error (for a short simulation of 100 ps/replica) using the default OpenMM PME cutoff of 5e-4.\r\n```\r\n2016-11-26 09:19:34,109: dDeltaf_ij:\r\n2016-11-26 09:19:34,109:    0.000   0.204   0.193   0.207   0.198   0.198   0.200   0.190   0.206   0.203   0.200   0.205   0.191   0.206   0.206   0.204\r\n2016-11-26 09:19:34,109:    0.204   0.000   0.198   0.199   0.207   0.192   0.196   0.196   0.202   0.196   0.195   0.200   0.187   0.200   0.205   0.214\r\n2016-11-26 09:19:34,109:    0.193   0.198   0.000   0.205   0.204   0.198   0.196   0.192   0.206   0.190   0.196   0.200   0.193   0.202   0.209   0.210\r\n2016-11-26 09:19:34,109:    0.207   0.199   0.205   0.000   0.211   0.192   0.204   0.196   0.209   0.197   0.204   0.209   0.193   0.205   0.208   0.215\r\n2016-11-26 09:19:34,110:    0.198   0.207   0.204   0.211   0.000   0.202   0.192   0.202   0.205   0.200   0.199   0.209   0.195   0.204   0.204   0.215\r\n2016-11-26 09:19:34,110:    0.198   0.192   0.198   0.192   0.202   0.000   0.197   0.179   0.207   0.187   0.195   0.195   0.181   0.191   0.191   0.202\r\n2016-11-26 09:19:34,110:    0.200   0.196   0.196   0.204   0.192   0.197   0.000   0.195   0.206   0.195   0.193   0.204   0.190   0.199   0.204   0.211\r\n2016-11-26 09:19:34,110:    0.190   0.196   0.192   0.196   0.202   0.179   0.195   0.000   0.209   0.188   0.198   0.192   0.183   0.191   0.200   0.201\r\n2016-11-26 09:19:34,110:    0.206   0.202   0.206   0.209   0.205   0.207   0.206   0.209   0.000   0.211   0.204   0.210   0.195   0.206   0.212   0.211\r\n2016-11-26 09:19:34,110:    0.203   0.196   0.190   0.197   0.200   0.187   0.195   0.188   0.211   0.000   0.198   0.185   0.191   0.190   0.203   0.208\r\n2016-11-26 09:19:34,111:    0.200   0.195   0.196   0.204   0.199   0.195   0.193   0.198   0.204   0.198   0.000   0.207   0.180   0.205   0.197   0.210\r\n2016-11-26 09:19:34,111:    0.205   0.200   0.200   0.209   0.209   0.195   0.204   0.192   0.210   0.185   0.207   0.000   0.201   0.197   0.211   0.206\r\n2016-11-26 09:19:34,111:    0.191   0.187   0.193   0.193   0.195   0.181   0.190   0.183   0.195   0.191   0.180   0.201   0.000   0.194   0.190   0.202\r\n2016-11-26 09:19:34,111:    0.206   0.200   0.202   0.205   0.204   0.191   0.199   0.191   0.206   0.190   0.205   0.197   0.194   0.000   0.214   0.203\r\n2016-11-26 09:19:34,111:    0.206   0.205   0.209   0.208   0.204   0.191   0.204   0.200   0.212   0.203   0.197   0.211   0.190   0.214   0.000   0.223\r\n2016-11-26 09:19:34,112:    0.204   0.214   0.210   0.215   0.215   0.202   0.211   0.201   0.211   0.208   0.210   0.206   0.202   0.203   0.223   0.000\r\n```\r\nReducing the PME tolerance to just 1e-4 helps significantly, reducing the error from 9A >> 24A by a factor of 3:\r\n```\r\n2016-11-26 09:20:28,283: dDeltaf_ij:\r\n2016-11-26 09:20:28,284:    0.000   0.052   0.057   0.062   0.064   0.069   0.069   0.072   0.074   0.073   0.071   0.075   0.071   0.076   0.075   0.075\r\n2016-11-26 09:20:28,285:    0.052   0.000   0.040   0.039   0.046   0.048   0.050   0.052   0.053   0.052   0.050   0.055   0.052   0.055   0.056   0.054\r\n2016-11-26 09:20:28,285:    0.057   0.040   0.000   0.036   0.031   0.040   0.036   0.040   0.042   0.040   0.037   0.041   0.039   0.044   0.044   0.043\r\n2016-11-26 09:20:28,286:    0.062   0.039   0.036   0.000   0.035   0.033   0.039   0.037   0.042   0.037   0.037   0.041   0.039   0.041   0.043   0.040\r\n2016-11-26 09:20:28,286:    0.064   0.046   0.031   0.035   0.000   0.036   0.033   0.039   0.037   0.036   0.035   0.036   0.036   0.037   0.040   0.037\r\n2016-11-26 09:20:28,287:    0.069   0.048   0.040   0.033   0.036   0.000   0.040   0.035   0.041   0.035   0.036   0.039   0.037   0.039   0.041   0.039\r\n2016-11-26 09:20:28,288:    0.069   0.050   0.036   0.039   0.033   0.040   0.000   0.040   0.035   0.039   0.034   0.038   0.036   0.038   0.037   0.038\r\n2016-11-26 09:20:28,288:    0.072   0.052   0.040   0.037   0.039   0.035   0.040   0.000   0.042   0.035   0.037   0.038   0.037   0.039   0.041   0.040\r\n2016-11-26 09:20:28,289:    0.074   0.053   0.042   0.042   0.037   0.041   0.035   0.042   0.000   0.044   0.035   0.042   0.037   0.040   0.038   0.039\r\n2016-11-26 09:20:28,290:    0.073   0.052   0.040   0.037   0.036   0.035   0.039   0.035   0.044   0.000   0.038   0.035   0.039   0.035   0.040   0.037\r\n2016-11-26 09:20:28,290:    0.071   0.050   0.037   0.037   0.035   0.036   0.034   0.037   0.035   0.038   0.000   0.041   0.033   0.039   0.037   0.036\r\n2016-11-26 09:20:28,291:    0.075   0.055   0.041   0.041   0.036   0.039   0.038   0.038   0.042   0.035   0.041   0.000   0.041   0.034   0.043   0.039\r\n2016-11-26 09:20:28,293:    0.071   0.052   0.039   0.039   0.036   0.037   0.036   0.037   0.037   0.039   0.033   0.041   0.000   0.042   0.036   0.040\r\n2016-11-26 09:20:28,293:    0.076   0.055   0.044   0.041   0.037   0.039   0.038   0.039   0.040   0.035   0.039   0.034   0.042   0.000   0.044   0.038\r\n2016-11-26 09:20:28,294:    0.075   0.056   0.044   0.043   0.040   0.041   0.037   0.041   0.038   0.040   0.037   0.043   0.036   0.044   0.000   0.043\r\n2016-11-26 09:20:28,295:    0.075   0.054   0.043   0.040   0.037   0.039   0.038   0.040   0.039   0.037   0.036   0.039   0.040   0.038   0.043   0.000\r\n```\r\nThis comes with a performance loss, however---for Src, ns/day goes from 43 to 27 ns/day on a (slow) GTX-680 for a 9A cutoff. The optimally performant cutoff for 1e-4 tolerance is actually 10-12A, so if we go to a 12A cutoff, we actually see the error is 5x smaller than before.\r\n\r\nIn summary, I recommend we take the following steps:\r\n* Expose the PME tolerance as an option (likely in the `solvent` section?)\r\n* Set the default PME tolerance in YANK to 1e-4 or smaller\r\n* Retool the examples to use a cutoff of 10-12A (perhaps 12A)\r\n\r\nIt may also be useful to expose Lennard-Jones switching and the switch width (where the switching distance is `cutoff - switch_width`) as options while we're at it, with the recommendation that the switch width be 1 A and the switch be on by default for better behavior."}}
{"actor":{"display_login":"jchodera"}}