Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 7/15/18 9:14 PM, Eric Smoll wrote: Justin, No problem. That makes sense. Can I send my atp, rtp, itp, and gro files used with the pdb2gmx from your drude branch off-list? Yes, please send a full, self-contained example with any relevant commands or information to reproduce whatever problem you're having. -Justin Best, Eric On Fri, Jul 13, 2018 at 8:12 PM, Justin Lemkul wrote: On 7/13/18 8:14 PM, Eric Smoll wrote: Justin, Your pdb2gmx appears to exclude all intermolecular interactions for a molecule with no hydrogen atoms. For instance, for a molecule with "N" atoms indexed from "1" to "N," the first line of the generated exclusions directive has a record with the "1 2 3...N" series. Is this expected? Shouldn't intramolecular nonbonded interactions be permitted at and beyond 1-4 interactions? I don't want to guess based on files I haven't seen. There's not much I can tell you. -Justin Best, Eric On Fri, Jul 13, 2018 at 3:02 PM, Eric Smoll wrote: Hi Justin, Very grateful for the rapid reply and warning. If you suggest that I should tread lightly with polarizable H atoms, I will avoid it altogether. I will alter my troubleshooting plan and focus on using your edited pdb2gmx program (adding to the rtp, of course) to build a core-shell MD topology for my problem molecules. It should be easy to compare with the topology I have prepared with my own tools. They should be identical. Best, Eric On Fri, Jul 13, 2018 at 2:51 PM, Justin Lemkul wrote: On 7/13/18 4:48 PM, Eric Smoll wrote: Hello Justin, Thank you for the guidance. Although more testing is needed, molecules where shells are attached to every atom appear to be working properly. The issue appears to be with models where shells are only attached to heavy atoms and bonds between hydrogen atoms and heavy atoms are constrained. Initial energy minimization tests (emtol = 100) show that bonds between heavy atoms and hydrogen in organic molecules slowly stretch or compress far more than they should while the network of non-hydrogen atoms maintains a sensible geometry. All atom-drude displacements on the heavy atoms converge quickly and are stable. I am still hunting for whatever topology problem is causing this. All exclusions out to 1-4 interactions between atoms and drudes should be properly included (combined action of the moleculetype directive and additional exclusions directives). Thole screening only applies to atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H atoms (no attached shells). It may be difficult to provide any useful guidance without details but troubleshooting suggestions are welcome if you have any. Are there issues associated with adding shells to specific atoms in a molecule? Are simulations that place shells on all atoms (hydrogen and heavy) more stable for some reason? I am going to build a new model where shells are attached to all atoms to see if bonds to hydrogen atoms still slowly compress/stretch during a tight-emtol energy minimization. I have never tested a system like that. Our Drude convention does exactly what you seem to find a problem - Drudes on heavy atoms and not H, with bonds to H constrained. I've never had the issue you're experiencing. Without a full test case of complete inputs, I can't tell you anything. Be forewarned - I know of no efforts by me or anyone else to deal with polarizable H atoms in GROMACS, and that may conflict with constraints, so tread lightly... -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support /Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support /Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Justin, No problem. That makes sense. Can I send my atp, rtp, itp, and gro files used with the pdb2gmx from your drude branch off-list? Best, Eric On Fri, Jul 13, 2018 at 8:12 PM, Justin Lemkul wrote: > > > On 7/13/18 8:14 PM, Eric Smoll wrote: > >> Justin, >> >> Your pdb2gmx appears to exclude all intermolecular interactions for a >> molecule with no hydrogen atoms. For instance, for a molecule with "N" >> atoms indexed from "1" to "N," the first line of the generated exclusions >> directive has a record with the "1 2 3...N" series. >> >> Is this expected? Shouldn't intramolecular nonbonded interactions be >> permitted at and beyond 1-4 interactions? >> > > I don't want to guess based on files I haven't seen. There's not much I > can tell you. > > -Justin > > > Best, >> Eric >> >> On Fri, Jul 13, 2018 at 3:02 PM, Eric Smoll wrote: >> >> Hi Justin, >>> >>> Very grateful for the rapid reply and warning. If you suggest that I >>> should tread lightly with polarizable H atoms, I will avoid it >>> altogether. >>> >>> I will alter my troubleshooting plan and focus on using your edited >>> pdb2gmx program (adding to the rtp, of course) to build a core-shell MD >>> topology for my problem molecules. It should be easy to compare with the >>> topology I have prepared with my own tools. They should be identical. >>> >>> Best, >>> Eric >>> >>> On Fri, Jul 13, 2018 at 2:51 PM, Justin Lemkul wrote: >>> >>> On 7/13/18 4:48 PM, Eric Smoll wrote: Hello Justin, > > Thank you for the guidance. > > Although more testing is needed, molecules where shells are attached to > every atom appear to be working properly. > > The issue appears to be with models where shells are only attached to > heavy > atoms and bonds between hydrogen atoms and heavy atoms are constrained. > Initial energy minimization tests (emtol = 100) show that bonds between > heavy atoms and hydrogen in organic molecules slowly stretch or > compress > far more than they should while the network of non-hydrogen atoms > maintains > a sensible geometry. All atom-drude displacements on the heavy atoms > converge quickly and are stable. > > I am still hunting for whatever topology problem is causing this. All > exclusions out to 1-4 interactions between atoms and drudes should be > properly included (combined action of the moleculetype directive and > additional exclusions directives). Thole screening only applies to > atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H > atoms (no attached shells). It may be difficult to provide any useful > guidance without details but troubleshooting suggestions are welcome if > you > have any. > > Are there issues associated with adding shells to specific atoms in a > molecule? Are simulations that place shells on all atoms (hydrogen and > heavy) more stable for some reason? I am going to build a new model > where > shells are attached to all atoms to see if bonds to hydrogen atoms > still > slowly compress/stretch during a tight-emtol energy minimization. > > I have never tested a system like that. Our Drude convention does exactly what you seem to find a problem - Drudes on heavy atoms and not H, with bonds to H constrained. I've never had the issue you're experiencing. Without a full test case of complete inputs, I can't tell you anything. Be forewarned - I know of no efforts by me or anyone else to deal with polarizable H atoms in GROMACS, and that may conflict with constraints, so tread lightly... -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support /Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. >>> > -- > == > > Justin A. Lemkul, Ph.D. > Assistant Professor > Virginia Tech Department of Biochemistry > > 303 Engel Hall > 340 West Campus Dr. > Blacksburg, VA 24061 > > jalem...@vt.edu | (540) 231-3129 > http://www.thelemkullab.com > > == > > -- > Gromacs Users mailing list > > * Please search the archive at http://www.gromacs.org/Support > /Mailing_Lists/GMX-Users_List before
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 7/13/18 9:42 PM, Eric Smoll wrote: Hi Justin, Using your code to generate the topology for my system still results in the distortion of bonds to hydrogen when constraints=h-bonds and emtol=100 (note that the energy minimization mdp files provided with your paper use an emtol of 1000). However, when I set constraints=none and emtol=100, no distortion of bonds to hydrogen occur. Maybe there is some problem using shells with constraints. I am going to remove all constraints from my system and check if my system is stable with SCF and Lagrangian core-shell molecular dynamics. You're going to have to use a very short time step, but the constraint code is entirely separate from anything dealing with polarization, so I have no clue why the two would be connected. Again, happy to troubleshoot a real system but it's really hard for me to guess at what's going on. -Justin Best, Eric On Fri, Jul 13, 2018 at 6:14 PM, Eric Smoll wrote: Justin, Your pdb2gmx appears to exclude all intermolecular interactions for a molecule with no hydrogen atoms. For instance, for a molecule with "N" atoms indexed from "1" to "N," the first line of the generated exclusions directive has a record with the "1 2 3...N" series. Is this expected? Shouldn't intramolecular nonbonded interactions be permitted at and beyond 1-4 interactions? Best, Eric On Fri, Jul 13, 2018 at 3:02 PM, Eric Smoll wrote: Hi Justin, Very grateful for the rapid reply and warning. If you suggest that I should tread lightly with polarizable H atoms, I will avoid it altogether. I will alter my troubleshooting plan and focus on using your edited pdb2gmx program (adding to the rtp, of course) to build a core-shell MD topology for my problem molecules. It should be easy to compare with the topology I have prepared with my own tools. They should be identical. Best, Eric On Fri, Jul 13, 2018 at 2:51 PM, Justin Lemkul wrote: On 7/13/18 4:48 PM, Eric Smoll wrote: Hello Justin, Thank you for the guidance. Although more testing is needed, molecules where shells are attached to every atom appear to be working properly. The issue appears to be with models where shells are only attached to heavy atoms and bonds between hydrogen atoms and heavy atoms are constrained. Initial energy minimization tests (emtol = 100) show that bonds between heavy atoms and hydrogen in organic molecules slowly stretch or compress far more than they should while the network of non-hydrogen atoms maintains a sensible geometry. All atom-drude displacements on the heavy atoms converge quickly and are stable. I am still hunting for whatever topology problem is causing this. All exclusions out to 1-4 interactions between atoms and drudes should be properly included (combined action of the moleculetype directive and additional exclusions directives). Thole screening only applies to atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H atoms (no attached shells). It may be difficult to provide any useful guidance without details but troubleshooting suggestions are welcome if you have any. Are there issues associated with adding shells to specific atoms in a molecule? Are simulations that place shells on all atoms (hydrogen and heavy) more stable for some reason? I am going to build a new model where shells are attached to all atoms to see if bonds to hydrogen atoms still slowly compress/stretch during a tight-emtol energy minimization. I have never tested a system like that. Our Drude convention does exactly what you seem to find a problem - Drudes on heavy atoms and not H, with bonds to H constrained. I've never had the issue you're experiencing. Without a full test case of complete inputs, I can't tell you anything. Be forewarned - I know of no efforts by me or anyone else to deal with polarizable H atoms in GROMACS, and that may conflict with constraints, so tread lightly... -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support /Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 7/13/18 8:14 PM, Eric Smoll wrote: Justin, Your pdb2gmx appears to exclude all intermolecular interactions for a molecule with no hydrogen atoms. For instance, for a molecule with "N" atoms indexed from "1" to "N," the first line of the generated exclusions directive has a record with the "1 2 3...N" series. Is this expected? Shouldn't intramolecular nonbonded interactions be permitted at and beyond 1-4 interactions? I don't want to guess based on files I haven't seen. There's not much I can tell you. -Justin Best, Eric On Fri, Jul 13, 2018 at 3:02 PM, Eric Smoll wrote: Hi Justin, Very grateful for the rapid reply and warning. If you suggest that I should tread lightly with polarizable H atoms, I will avoid it altogether. I will alter my troubleshooting plan and focus on using your edited pdb2gmx program (adding to the rtp, of course) to build a core-shell MD topology for my problem molecules. It should be easy to compare with the topology I have prepared with my own tools. They should be identical. Best, Eric On Fri, Jul 13, 2018 at 2:51 PM, Justin Lemkul wrote: On 7/13/18 4:48 PM, Eric Smoll wrote: Hello Justin, Thank you for the guidance. Although more testing is needed, molecules where shells are attached to every atom appear to be working properly. The issue appears to be with models where shells are only attached to heavy atoms and bonds between hydrogen atoms and heavy atoms are constrained. Initial energy minimization tests (emtol = 100) show that bonds between heavy atoms and hydrogen in organic molecules slowly stretch or compress far more than they should while the network of non-hydrogen atoms maintains a sensible geometry. All atom-drude displacements on the heavy atoms converge quickly and are stable. I am still hunting for whatever topology problem is causing this. All exclusions out to 1-4 interactions between atoms and drudes should be properly included (combined action of the moleculetype directive and additional exclusions directives). Thole screening only applies to atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H atoms (no attached shells). It may be difficult to provide any useful guidance without details but troubleshooting suggestions are welcome if you have any. Are there issues associated with adding shells to specific atoms in a molecule? Are simulations that place shells on all atoms (hydrogen and heavy) more stable for some reason? I am going to build a new model where shells are attached to all atoms to see if bonds to hydrogen atoms still slowly compress/stretch during a tight-emtol energy minimization. I have never tested a system like that. Our Drude convention does exactly what you seem to find a problem - Drudes on heavy atoms and not H, with bonds to H constrained. I've never had the issue you're experiencing. Without a full test case of complete inputs, I can't tell you anything. Be forewarned - I know of no efforts by me or anyone else to deal with polarizable H atoms in GROMACS, and that may conflict with constraints, so tread lightly... -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support /Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Justin, Your pdb2gmx appears to exclude all intermolecular interactions for a molecule with no hydrogen atoms. For instance, for a molecule with "N" atoms indexed from "1" to "N," the first line of the generated exclusions directive has a record with the "1 2 3...N" series. Is this expected? Shouldn't intramolecular nonbonded interactions be permitted at and beyond 1-4 interactions? Best, Eric On Fri, Jul 13, 2018 at 3:02 PM, Eric Smoll wrote: > Hi Justin, > > Very grateful for the rapid reply and warning. If you suggest that I > should tread lightly with polarizable H atoms, I will avoid it altogether. > > I will alter my troubleshooting plan and focus on using your edited > pdb2gmx program (adding to the rtp, of course) to build a core-shell MD > topology for my problem molecules. It should be easy to compare with the > topology I have prepared with my own tools. They should be identical. > > Best, > Eric > > On Fri, Jul 13, 2018 at 2:51 PM, Justin Lemkul wrote: > >> >> >> On 7/13/18 4:48 PM, Eric Smoll wrote: >> >>> Hello Justin, >>> >>> Thank you for the guidance. >>> >>> Although more testing is needed, molecules where shells are attached to >>> every atom appear to be working properly. >>> >>> The issue appears to be with models where shells are only attached to >>> heavy >>> atoms and bonds between hydrogen atoms and heavy atoms are constrained. >>> Initial energy minimization tests (emtol = 100) show that bonds between >>> heavy atoms and hydrogen in organic molecules slowly stretch or compress >>> far more than they should while the network of non-hydrogen atoms >>> maintains >>> a sensible geometry. All atom-drude displacements on the heavy atoms >>> converge quickly and are stable. >>> >>> I am still hunting for whatever topology problem is causing this. All >>> exclusions out to 1-4 interactions between atoms and drudes should be >>> properly included (combined action of the moleculetype directive and >>> additional exclusions directives). Thole screening only applies to >>> atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H >>> atoms (no attached shells). It may be difficult to provide any useful >>> guidance without details but troubleshooting suggestions are welcome if >>> you >>> have any. >>> >>> Are there issues associated with adding shells to specific atoms in a >>> molecule? Are simulations that place shells on all atoms (hydrogen and >>> heavy) more stable for some reason? I am going to build a new model >>> where >>> shells are attached to all atoms to see if bonds to hydrogen atoms still >>> slowly compress/stretch during a tight-emtol energy minimization. >>> >> >> I have never tested a system like that. Our Drude convention does exactly >> what you seem to find a problem - Drudes on heavy atoms and not H, with >> bonds to H constrained. I've never had the issue you're experiencing. >> Without a full test case of complete inputs, I can't tell you anything. >> >> Be forewarned - I know of no efforts by me or anyone else to deal with >> polarizable H atoms in GROMACS, and that may conflict with constraints, so >> tread lightly... >> >> -Justin >> >> -- >> == >> >> Justin A. Lemkul, Ph.D. >> Assistant Professor >> Virginia Tech Department of Biochemistry >> >> 303 Engel Hall >> 340 West Campus Dr. >> Blacksburg, VA 24061 >> >> jalem...@vt.edu | (540) 231-3129 >> http://www.thelemkullab.com >> >> == >> >> -- >> Gromacs Users mailing list >> >> * Please search the archive at http://www.gromacs.org/Support >> /Mailing_Lists/GMX-Users_List before posting! >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> * For (un)subscribe requests visit >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >> send a mail to gmx-users-requ...@gromacs.org. >> > > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Hi Justin, Very grateful for the rapid reply and warning. If you suggest that I should tread lightly with polarizable H atoms, I will avoid it altogether. I will alter my troubleshooting plan and focus on using your edited pdb2gmx program (adding to the rtp, of course) to build a core-shell MD topology for my problem molecules. It should be easy to compare with the topology I have prepared with my own tools. They should be identical. Best, Eric On Fri, Jul 13, 2018 at 2:51 PM, Justin Lemkul wrote: > > > On 7/13/18 4:48 PM, Eric Smoll wrote: > >> Hello Justin, >> >> Thank you for the guidance. >> >> Although more testing is needed, molecules where shells are attached to >> every atom appear to be working properly. >> >> The issue appears to be with models where shells are only attached to >> heavy >> atoms and bonds between hydrogen atoms and heavy atoms are constrained. >> Initial energy minimization tests (emtol = 100) show that bonds between >> heavy atoms and hydrogen in organic molecules slowly stretch or compress >> far more than they should while the network of non-hydrogen atoms >> maintains >> a sensible geometry. All atom-drude displacements on the heavy atoms >> converge quickly and are stable. >> >> I am still hunting for whatever topology problem is causing this. All >> exclusions out to 1-4 interactions between atoms and drudes should be >> properly included (combined action of the moleculetype directive and >> additional exclusions directives). Thole screening only applies to >> atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H >> atoms (no attached shells). It may be difficult to provide any useful >> guidance without details but troubleshooting suggestions are welcome if >> you >> have any. >> >> Are there issues associated with adding shells to specific atoms in a >> molecule? Are simulations that place shells on all atoms (hydrogen and >> heavy) more stable for some reason? I am going to build a new model where >> shells are attached to all atoms to see if bonds to hydrogen atoms still >> slowly compress/stretch during a tight-emtol energy minimization. >> > > I have never tested a system like that. Our Drude convention does exactly > what you seem to find a problem - Drudes on heavy atoms and not H, with > bonds to H constrained. I've never had the issue you're experiencing. > Without a full test case of complete inputs, I can't tell you anything. > > Be forewarned - I know of no efforts by me or anyone else to deal with > polarizable H atoms in GROMACS, and that may conflict with constraints, so > tread lightly... > > -Justin > > -- > == > > Justin A. Lemkul, Ph.D. > Assistant Professor > Virginia Tech Department of Biochemistry > > 303 Engel Hall > 340 West Campus Dr. > Blacksburg, VA 24061 > > jalem...@vt.edu | (540) 231-3129 > http://www.thelemkullab.com > > == > > -- > Gromacs Users mailing list > > * Please search the archive at http://www.gromacs.org/Support > /Mailing_Lists/GMX-Users_List before posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 7/13/18 4:48 PM, Eric Smoll wrote: Hello Justin, Thank you for the guidance. Although more testing is needed, molecules where shells are attached to every atom appear to be working properly. The issue appears to be with models where shells are only attached to heavy atoms and bonds between hydrogen atoms and heavy atoms are constrained. Initial energy minimization tests (emtol = 100) show that bonds between heavy atoms and hydrogen in organic molecules slowly stretch or compress far more than they should while the network of non-hydrogen atoms maintains a sensible geometry. All atom-drude displacements on the heavy atoms converge quickly and are stable. I am still hunting for whatever topology problem is causing this. All exclusions out to 1-4 interactions between atoms and drudes should be properly included (combined action of the moleculetype directive and additional exclusions directives). Thole screening only applies to atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H atoms (no attached shells). It may be difficult to provide any useful guidance without details but troubleshooting suggestions are welcome if you have any. Are there issues associated with adding shells to specific atoms in a molecule? Are simulations that place shells on all atoms (hydrogen and heavy) more stable for some reason? I am going to build a new model where shells are attached to all atoms to see if bonds to hydrogen atoms still slowly compress/stretch during a tight-emtol energy minimization. I have never tested a system like that. Our Drude convention does exactly what you seem to find a problem - Drudes on heavy atoms and not H, with bonds to H constrained. I've never had the issue you're experiencing. Without a full test case of complete inputs, I can't tell you anything. Be forewarned - I know of no efforts by me or anyone else to deal with polarizable H atoms in GROMACS, and that may conflict with constraints, so tread lightly... -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Hello Justin, Thank you for the guidance. Although more testing is needed, molecules where shells are attached to every atom appear to be working properly. The issue appears to be with models where shells are only attached to heavy atoms and bonds between hydrogen atoms and heavy atoms are constrained. Initial energy minimization tests (emtol = 100) show that bonds between heavy atoms and hydrogen in organic molecules slowly stretch or compress far more than they should while the network of non-hydrogen atoms maintains a sensible geometry. All atom-drude displacements on the heavy atoms converge quickly and are stable. I am still hunting for whatever topology problem is causing this. All exclusions out to 1-4 interactions between atoms and drudes should be properly included (combined action of the moleculetype directive and additional exclusions directives). Thole screening only applies to atom-shell pairs so no screening is needed for 1-2 and 1-3 pairs with H atoms (no attached shells). It may be difficult to provide any useful guidance without details but troubleshooting suggestions are welcome if you have any. Are there issues associated with adding shells to specific atoms in a molecule? Are simulations that place shells on all atoms (hydrogen and heavy) more stable for some reason? I am going to build a new model where shells are attached to all atoms to see if bonds to hydrogen atoms still slowly compress/stretch during a tight-emtol energy minimization. Best, Eric On Mon, Jul 9, 2018 at 6:28 AM, Justin Lemkul wrote: > > > On 7/9/18 1:31 AM, Eric Smoll wrote: > >> Thanks Justin, >> >> That is very helpful. I can run the water example provided in your paper >> with the GROMACS "drude" branch. >> >> I have built my own tools to construct a core-shell MD topology file for >> my >> system. Grompp produces a tpr without complaint but mdrun crashes after >> the first step with a segmentation fault and no other useful errors. I am >> trying to troubleshoot. >> >> Is it required that each drude/shell immediately follow the atom it is >> attached to in the atoms directive and the input coordinate file? For >> convenience, I have attached them to the end of each molecule? >> > > There is no such requirement (you'll note that the water Drude is at the > end of the atom listing in SWM4-NDP, to avoid changes in the SETTLE code). > > In my topology, I manually account for the fact that drude/shell particles >> are attached with bonds which count when evaluating exclusions. >> Specifically, I assume a "3" in the molecules directive will take care of >> all exclusions up to 3 bonds away. To account for the extra bond lengths >> introduced by drude/shell particles, I add manually add certain exclusions >> for particles that should be 3 bonds away according to the atom-atom >> connectivity. Will this work or is there some "gotcha" I don't >> understand. For example, do you invalidate the moleculetype directive by >> providing an exclusions directive? >> > > Sounds right, but without an actual example, I can't tell you anything > really useful. > > -Justin > > > >> Best, >> Eric >> >> On Tue, Jul 3, 2018 at 5:39 PM, Justin Lemkul wrote: >> >> >>> On 7/3/18 6:52 PM, Eric Smoll wrote: >>> >>> Justin, Thanks again for the response and the guidance. Your 2015 JCC paper on the subject was very helpful. Your code uses a different thole_polarization format. Currently, Gromacs uses the following format [ thole_polarization ] ; i j k l function-number thole-a-parameter alpha_ij alpha_kl Your code seems to use [ thole_polarization ] ; i j k l function-number alpha_ij alpha_kl thole-a-parameter-1 thole-a-parameter-2 If I would use 2.6 in the current Gromacs thole_polarization format, what would I use in your thole_polarization format? 1.3 for both atoms to get 1.3 + 1.3 = 2.6. Thole's magic number quickly >>> breaks down for complicated molecules, so our convention is that it is a >>> per-atom value that is summed to give a new value of "a" per pair, e.g. >>> a_i >>> + a_j = a_total. >>> >>> -Justin >>> >>> -- >>> == >>> >>> Justin A. Lemkul, Ph.D. >>> Assistant Professor >>> Virginia Tech Department of Biochemistry >>> >>> 303 Engel Hall >>> 340 West Campus Dr. >>> Blacksburg, VA 24061 >>> >>> jalem...@vt.edu | (540) 231-3129 >>> http://www.thelemkullab.com >>> >>> == >>> >>> -- >>> Gromacs Users mailing list >>> >>> * Please search the archive at http://www.gromacs.org/Support >>> /Mailing_Lists/GMX-Users_List before posting! >>> >>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >>> >>> * For (un)subscribe requests visit >>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >>> send a mail to gmx-users-requ...@gromacs.org. >>> >>>
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 7/9/18 1:31 AM, Eric Smoll wrote: Thanks Justin, That is very helpful. I can run the water example provided in your paper with the GROMACS "drude" branch. I have built my own tools to construct a core-shell MD topology file for my system. Grompp produces a tpr without complaint but mdrun crashes after the first step with a segmentation fault and no other useful errors. I am trying to troubleshoot. Is it required that each drude/shell immediately follow the atom it is attached to in the atoms directive and the input coordinate file? For convenience, I have attached them to the end of each molecule? There is no such requirement (you'll note that the water Drude is at the end of the atom listing in SWM4-NDP, to avoid changes in the SETTLE code). In my topology, I manually account for the fact that drude/shell particles are attached with bonds which count when evaluating exclusions. Specifically, I assume a "3" in the molecules directive will take care of all exclusions up to 3 bonds away. To account for the extra bond lengths introduced by drude/shell particles, I add manually add certain exclusions for particles that should be 3 bonds away according to the atom-atom connectivity. Will this work or is there some "gotcha" I don't understand. For example, do you invalidate the moleculetype directive by providing an exclusions directive? Sounds right, but without an actual example, I can't tell you anything really useful. -Justin Best, Eric On Tue, Jul 3, 2018 at 5:39 PM, Justin Lemkul wrote: On 7/3/18 6:52 PM, Eric Smoll wrote: Justin, Thanks again for the response and the guidance. Your 2015 JCC paper on the subject was very helpful. Your code uses a different thole_polarization format. Currently, Gromacs uses the following format [ thole_polarization ] ; i j k l function-number thole-a-parameter alpha_ij alpha_kl Your code seems to use [ thole_polarization ] ; i j k l function-number alpha_ij alpha_kl thole-a-parameter-1 thole-a-parameter-2 If I would use 2.6 in the current Gromacs thole_polarization format, what would I use in your thole_polarization format? 1.3 for both atoms to get 1.3 + 1.3 = 2.6. Thole's magic number quickly breaks down for complicated molecules, so our convention is that it is a per-atom value that is summed to give a new value of "a" per pair, e.g. a_i + a_j = a_total. -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support /Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 7/3/18 6:52 PM, Eric Smoll wrote: Justin, Thanks again for the response and the guidance. Your 2015 JCC paper on the subject was very helpful. Your code uses a different thole_polarization format. Currently, Gromacs uses the following format [ thole_polarization ] ; i j k l function-number thole-a-parameter alpha_ij alpha_kl Your code seems to use [ thole_polarization ] ; i j k l function-number alpha_ij alpha_kl thole-a-parameter-1 thole-a-parameter-2 If I would use 2.6 in the current Gromacs thole_polarization format, what would I use in your thole_polarization format? 1.3 for both atoms to get 1.3 + 1.3 = 2.6. Thole's magic number quickly breaks down for complicated molecules, so our convention is that it is a per-atom value that is summed to give a new value of "a" per pair, e.g. a_i + a_j = a_total. -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Justin, Thanks again for the response and the guidance. Your 2015 JCC paper on the subject was very helpful. Your code uses a different thole_polarization format. Currently, Gromacs uses the following format [ thole_polarization ] ; i j k l function-number thole-a-parameter alpha_ij alpha_kl Your code seems to use [ thole_polarization ] ; i j k l function-number alpha_ij alpha_kl thole-a-parameter-1 thole-a-parameter-2 If I would use 2.6 in the current Gromacs thole_polarization format, what would I use in your thole_polarization format? Best, Eric On Sun, Jul 1, 2018 at 6:52 PM, Justin Lemkul wrote: > > > On 7/1/18 6:37 PM, Eric Smoll wrote: > >> Justin, >> >> Thank you for the response! I deeply appreciate the effort you put into >> this community. >> >> I am still very confused *but* I learn more with every email. I am now >> certain that I need to add drude coordinates in the gro file. This was >> not >> specified anywhere in the manual. I have also learned that one should not >> completley remove coulomb interactions between nearby intramolecular >> core-shell pairs. Instead, I should leave them in with thole screening. >> You have also clarified that this is only necessary at short range but how >> short is not particularly clear. For an OPLSAA style forcefield, I assume >> it is reasonable to start by thole screening 1-2 and 1-3 coulomb >> interactions, then scaling 1-4 non-bonded interactions by 0.5 without >> thole >> screening, then 1-5 non-bonded interactions and beyond can be a full >> strength. >> > > You can develop a model with whatever convention you like. The advice I'm > giving you is what we do in our Drude force field, which works quite well. > > You mention that you don't really know how the polarization directive works >> since you don't use it. Can I implement drude polarizability with thole >> screening without the polarization directive using your method of just >> defining the core-shell pairs and manually attaching them with bonds? I >> > > Yes. > > like your method since it is less opaque. Do I still need to label the >> shell particles in each pair with an 'S' in the topology (also with zero >> mass and LJ params)? I suspect the shell label is still necessary to >> > > Yes. > >> instruct gromacs to iteratively optimize shell positions at every >> timestep. >> Your method also causes complications since there is no harmonic bond type >> that does not create a connectivity link used to auto-generate exclusions. >> Which brings me to my next cluster of questions >> > > In this case, you *do* want to use my Drude branch, because I coded > pdb2gmx to do precisely that. Please see the .rtp entries in the tarball > from http://mackerell.umaryland.edu/charmm_drude_ff.shtml > > If you specify in the header of the .rtp file that the force field is > polarizable, and specify alpha and Thole for each atom, pdb2gmx will do the > rest, with all proper pairs and exclusions. Perhaps before diving into > this, do something really easy, like ethane or something, for which you > could compute the energy even by hand so you know what the result should be. > > Is there a clean strategy that will prevent me from making a mistake on the >> intramolecular nonbonded interactions when I am using your method (no >> polarization directive)? I am struggling to understand what this strategy >> is. What needs to be excluded and how? There are many pieces here >> > > In our method, the Drudes inherit the same pair and exclusion list as > their parent atom. The Drude-atom bond does not count for the purposes of > determining 1-2, 1-3, and 1-4 interactions. > > -Justin > > > (defaults directive settings, exclusion directive settings, the necessity >> to scale 1-4 all LJ and coulomb interactions, the pairs directive, the >> moleculetype directive, and the thole _polarization directive). >> >> Sorry for the additional questions. I suspect this conversation will drag >> on for a long time. >> >> Best, >> Eric >> >> On Sun, Jul 1, 2018 at 1:27 PM, Justin Lemkul wrote: >> >> >>> On 7/1/18 3:13 PM, Eric Smoll wrote: >>> >>> Thank you, I have compiled the "Drude" branch and am preparing my topology. Unfortunately, the best reference for implementing shell molecular dynamics in Gromacs is an extremely long and confusing 2014 gmx-users thread between you and an extremely confused Gromacs user ( https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users /2014-July/091066.html). I have more questions. I am interested in adding polarization to an OPLS-AA style model. Let's build a model system to add clarity to our communication. Consider a linear, 5 physical-atom (not gromacs-atom) artificial molecule. To describe this artificial molecule in a shell molecular dynamics simulation we need to define 10 gromacs-particles. 5 of these gromacs-particles will be of the
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 7/1/18 6:37 PM, Eric Smoll wrote: Justin, Thank you for the response! I deeply appreciate the effort you put into this community. I am still very confused *but* I learn more with every email. I am now certain that I need to add drude coordinates in the gro file. This was not specified anywhere in the manual. I have also learned that one should not completley remove coulomb interactions between nearby intramolecular core-shell pairs. Instead, I should leave them in with thole screening. You have also clarified that this is only necessary at short range but how short is not particularly clear. For an OPLSAA style forcefield, I assume it is reasonable to start by thole screening 1-2 and 1-3 coulomb interactions, then scaling 1-4 non-bonded interactions by 0.5 without thole screening, then 1-5 non-bonded interactions and beyond can be a full strength. You can develop a model with whatever convention you like. The advice I'm giving you is what we do in our Drude force field, which works quite well. You mention that you don't really know how the polarization directive works since you don't use it. Can I implement drude polarizability with thole screening without the polarization directive using your method of just defining the core-shell pairs and manually attaching them with bonds? I Yes. like your method since it is less opaque. Do I still need to label the shell particles in each pair with an 'S' in the topology (also with zero mass and LJ params)? I suspect the shell label is still necessary to Yes. instruct gromacs to iteratively optimize shell positions at every timestep. Your method also causes complications since there is no harmonic bond type that does not create a connectivity link used to auto-generate exclusions. Which brings me to my next cluster of questions In this case, you *do* want to use my Drude branch, because I coded pdb2gmx to do precisely that. Please see the .rtp entries in the tarball from http://mackerell.umaryland.edu/charmm_drude_ff.shtml If you specify in the header of the .rtp file that the force field is polarizable, and specify alpha and Thole for each atom, pdb2gmx will do the rest, with all proper pairs and exclusions. Perhaps before diving into this, do something really easy, like ethane or something, for which you could compute the energy even by hand so you know what the result should be. Is there a clean strategy that will prevent me from making a mistake on the intramolecular nonbonded interactions when I am using your method (no polarization directive)? I am struggling to understand what this strategy is. What needs to be excluded and how? There are many pieces here In our method, the Drudes inherit the same pair and exclusion list as their parent atom. The Drude-atom bond does not count for the purposes of determining 1-2, 1-3, and 1-4 interactions. -Justin (defaults directive settings, exclusion directive settings, the necessity to scale 1-4 all LJ and coulomb interactions, the pairs directive, the moleculetype directive, and the thole _polarization directive). Sorry for the additional questions. I suspect this conversation will drag on for a long time. Best, Eric On Sun, Jul 1, 2018 at 1:27 PM, Justin Lemkul wrote: On 7/1/18 3:13 PM, Eric Smoll wrote: Thank you, I have compiled the "Drude" branch and am preparing my topology. Unfortunately, the best reference for implementing shell molecular dynamics in Gromacs is an extremely long and confusing 2014 gmx-users thread between you and an extremely confused Gromacs user ( https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users /2014-July/091066.html). I have more questions. I am interested in adding polarization to an OPLS-AA style model. Let's build a model system to add clarity to our communication. Consider a linear, 5 physical-atom (not gromacs-atom) artificial molecule. To describe this artificial molecule in a shell molecular dynamics simulation we need to define 10 gromacs-particles. 5 of these gromacs-particles will be of the gromacs-atom type (marked with an A in the topology) and 5 of these gromacs-particle will be of the gromacs-shell type (marked with an S in the topology). Each gromacs-shell will be "attached" to a corresponding gromacs-atom type as discussed in (1) below. (1) As you clarified earlier, I create a polarization directive specifying I will begin my remarks by saying I do not do any simulations that use a [polarization] directive. In our Drude model, everything is just an atom defined under [atoms] with explicit bonds listed in [bonds]. So I may not have exact answers to everything you want to know. Further, you don't necessarily need my Drude branch to deal with anything you've mentioned here. gromacs-atom/gromacs-shell pairs. For each entry in the polarization directive, I specify the function number as 1 and the polarizability in units of nm^3. This "attaches" the gromacs-shell to the gromacs-atom with
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Justin, Thank you for the response! I deeply appreciate the effort you put into this community. I am still very confused *but* I learn more with every email. I am now certain that I need to add drude coordinates in the gro file. This was not specified anywhere in the manual. I have also learned that one should not completley remove coulomb interactions between nearby intramolecular core-shell pairs. Instead, I should leave them in with thole screening. You have also clarified that this is only necessary at short range but how short is not particularly clear. For an OPLSAA style forcefield, I assume it is reasonable to start by thole screening 1-2 and 1-3 coulomb interactions, then scaling 1-4 non-bonded interactions by 0.5 without thole screening, then 1-5 non-bonded interactions and beyond can be a full strength. You mention that you don't really know how the polarization directive works since you don't use it. Can I implement drude polarizability with thole screening without the polarization directive using your method of just defining the core-shell pairs and manually attaching them with bonds? I like your method since it is less opaque. Do I still need to label the shell particles in each pair with an 'S' in the topology (also with zero mass and LJ params)? I suspect the shell label is still necessary to instruct gromacs to iteratively optimize shell positions at every timestep. Your method also causes complications since there is no harmonic bond type that does not create a connectivity link used to auto-generate exclusions. Which brings me to my next cluster of questions Is there a clean strategy that will prevent me from making a mistake on the intramolecular nonbonded interactions when I am using your method (no polarization directive)? I am struggling to understand what this strategy is. What needs to be excluded and how? There are many pieces here (defaults directive settings, exclusion directive settings, the necessity to scale 1-4 all LJ and coulomb interactions, the pairs directive, the moleculetype directive, and the thole _polarization directive). Sorry for the additional questions. I suspect this conversation will drag on for a long time. Best, Eric On Sun, Jul 1, 2018 at 1:27 PM, Justin Lemkul wrote: > > > On 7/1/18 3:13 PM, Eric Smoll wrote: > >> Thank you, >> >> I have compiled the "Drude" branch and am preparing my topology. >> Unfortunately, the best reference for implementing shell molecular >> dynamics >> in Gromacs is an extremely long and confusing 2014 gmx-users thread >> between >> you and an extremely confused Gromacs user ( >> https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users >> /2014-July/091066.html). >> I have more questions. >> >> I am interested in adding polarization to an OPLS-AA style model. Let's >> build a model system to add clarity to our communication. Consider a >> linear, 5 physical-atom (not gromacs-atom) artificial molecule. To >> describe this artificial molecule in a shell molecular dynamics simulation >> we need to define 10 gromacs-particles. 5 of these gromacs-particles will >> be of the gromacs-atom type (marked with an A in the topology) and 5 of >> these gromacs-particle will be of the gromacs-shell type (marked with an S >> in the topology). Each gromacs-shell will be "attached" to a >> corresponding >> gromacs-atom type as discussed in (1) below. >> >> (1) As you clarified earlier, I create a polarization directive specifying >> > > I will begin my remarks by saying I do not do any simulations that use a > [polarization] directive. In our Drude model, everything is just an atom > defined under [atoms] with explicit bonds listed in [bonds]. So I may not > have exact answers to everything you want to know. Further, you don't > necessarily need my Drude branch to deal with anything you've mentioned > here. > > gromacs-atom/gromacs-shell pairs. For each entry in the polarization >> directive, I specify the function number as 1 and the polarizability in >> units of nm^3. This "attaches" the gromacs-shell to the gromacs-atom with >> a harmonic force. This interaction is not considered a "bonded >> interaction" meaning that the resulting connectivity >> gromacs-shell/gromacs-atom >> connectivity is not considered when excluding interactions at the distance >> specified in the molecule type directive. Only the network defined by the >> gromacs-atom/gromacs-atom connectivity is considered when excluding >> interactions. Assuming the exclusion distance is three bonds and I have >> scaled 1-4 interactions specified (defaults directive settings and pairs >> directive entries), what happens with the gromacs-particle interactions as >> a function of distance? >> 1 bond away: gromacs-particle/gromacs-particle non-bonded interactions >> excluded >> 2 bond away: gromacs-particle/gromacs-particle non-bonded interactions >> excluded >> > > I will note here that indeed nrexcl *should* work this way (e.g. the > Drude/shell
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 7/1/18 3:13 PM, Eric Smoll wrote: Thank you, I have compiled the "Drude" branch and am preparing my topology. Unfortunately, the best reference for implementing shell molecular dynamics in Gromacs is an extremely long and confusing 2014 gmx-users thread between you and an extremely confused Gromacs user ( https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2014-July/091066.html). I have more questions. I am interested in adding polarization to an OPLS-AA style model. Let's build a model system to add clarity to our communication. Consider a linear, 5 physical-atom (not gromacs-atom) artificial molecule. To describe this artificial molecule in a shell molecular dynamics simulation we need to define 10 gromacs-particles. 5 of these gromacs-particles will be of the gromacs-atom type (marked with an A in the topology) and 5 of these gromacs-particle will be of the gromacs-shell type (marked with an S in the topology). Each gromacs-shell will be "attached" to a corresponding gromacs-atom type as discussed in (1) below. (1) As you clarified earlier, I create a polarization directive specifying I will begin my remarks by saying I do not do any simulations that use a [polarization] directive. In our Drude model, everything is just an atom defined under [atoms] with explicit bonds listed in [bonds]. So I may not have exact answers to everything you want to know. Further, you don't necessarily need my Drude branch to deal with anything you've mentioned here. gromacs-atom/gromacs-shell pairs. For each entry in the polarization directive, I specify the function number as 1 and the polarizability in units of nm^3. This "attaches" the gromacs-shell to the gromacs-atom with a harmonic force. This interaction is not considered a "bonded interaction" meaning that the resulting connectivity gromacs-shell/gromacs-atom connectivity is not considered when excluding interactions at the distance specified in the molecule type directive. Only the network defined by the gromacs-atom/gromacs-atom connectivity is considered when excluding interactions. Assuming the exclusion distance is three bonds and I have scaled 1-4 interactions specified (defaults directive settings and pairs directive entries), what happens with the gromacs-particle interactions as a function of distance? 1 bond away: gromacs-particle/gromacs-particle non-bonded interactions excluded 2 bond away: gromacs-particle/gromacs-particle non-bonded interactions excluded I will note here that indeed nrexcl *should* work this way (e.g. the Drude/shell inherits its parent atom's exclusions) but it's actually important in polarizable models to *include* all of these neighboring interactions, and this is what Thole screening is for. It's really important for obtaining a physically realistic molecular polarizability and is one of the significant advantages over additive models. 3 bond away: gromacs-particle/gromacs-particle non-bonded interactions included but scaled by 0.5. I am correct in assuming that gromacs-atom/gromacs-atom, gromacs-atom/gromacs-shell, and gromacs-shell/gromacs-shell coulomb interactions are all scaled down by a factor of 0.5? Presumably, but do check with an easily definable single-point energy. Our Drude model does work this way, but I am not 100% sure what a declaration via [polarization] does. I think that's how it works, too. 4 bond away: gromacs-particle/gromacs-particle non-bonded interactions included (gromacs-atom/gromacs-atom, gromacs-atom/gromacs-shell, and gromacs-shell/gromacs-shell) (2) Do I need to create .gro coordinates for the 10 gromacs-particles or just the 5 gromacs-atom? If I need to declare the coordinates for all 10 gromacs-particles, can I specify each gromacs-shell particle with the same coordinates as the gromacs-atom particle that it is attached to? Yes, you need the coordinates to be present. Initial assignment to the coordinates of the parent atom is fine. The dipoles will adjust as a function of the electric field in the system. (3) Now, according to the manual, thole polarization is intended to screen/damp/diminish coulomb interactions in shell molecular dynamics simulations. The wording suggests that thole polarization only impact gromacs-shell/gromacs-shell interactions (not gromacs-atom/gromacs-atom or gromacs-shell/gromacs-shell interactions). Am I correct? Each entry in No. Thole screening is applied between neighboring dipoles, not neighboring charges. Explicit charge-charge interactions are computed for all possible atom-Drude, atom-atom, and Drude-Drude pairs in the listed interaction and are then screened via the specified constant. the associated directive requires listing 4 numbers. The first two numbers specify a gromacs-atom and its attached gromacs-shell. The last two numbers specify a different gromacs-atom and its attached gromacs-shell. Given my assumption in (1) above, I would only need to include these for all 1-4, 1-5,
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Thank you, I have compiled the "Drude" branch and am preparing my topology. Unfortunately, the best reference for implementing shell molecular dynamics in Gromacs is an extremely long and confusing 2014 gmx-users thread between you and an extremely confused Gromacs user ( https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2014-July/091066.html). I have more questions. I am interested in adding polarization to an OPLS-AA style model. Let's build a model system to add clarity to our communication. Consider a linear, 5 physical-atom (not gromacs-atom) artificial molecule. To describe this artificial molecule in a shell molecular dynamics simulation we need to define 10 gromacs-particles. 5 of these gromacs-particles will be of the gromacs-atom type (marked with an A in the topology) and 5 of these gromacs-particle will be of the gromacs-shell type (marked with an S in the topology). Each gromacs-shell will be "attached" to a corresponding gromacs-atom type as discussed in (1) below. (1) As you clarified earlier, I create a polarization directive specifying gromacs-atom/gromacs-shell pairs. For each entry in the polarization directive, I specify the function number as 1 and the polarizability in units of nm^3. This "attaches" the gromacs-shell to the gromacs-atom with a harmonic force. This interaction is not considered a "bonded interaction" meaning that the resulting connectivity gromacs-shell/gromacs-atom connectivity is not considered when excluding interactions at the distance specified in the molecule type directive. Only the network defined by the gromacs-atom/gromacs-atom connectivity is considered when excluding interactions. Assuming the exclusion distance is three bonds and I have scaled 1-4 interactions specified (defaults directive settings and pairs directive entries), what happens with the gromacs-particle interactions as a function of distance? 1 bond away: gromacs-particle/gromacs-particle non-bonded interactions excluded 2 bond away: gromacs-particle/gromacs-particle non-bonded interactions excluded 3 bond away: gromacs-particle/gromacs-particle non-bonded interactions included but scaled by 0.5. I am correct in assuming that gromacs-atom/gromacs-atom, gromacs-atom/gromacs-shell, and gromacs-shell/gromacs-shell coulomb interactions are all scaled down by a factor of 0.5? 4 bond away: gromacs-particle/gromacs-particle non-bonded interactions included (gromacs-atom/gromacs-atom, gromacs-atom/gromacs-shell, and gromacs-shell/gromacs-shell) (2) Do I need to create .gro coordinates for the 10 gromacs-particles or just the 5 gromacs-atom? If I need to declare the coordinates for all 10 gromacs-particles, can I specify each gromacs-shell particle with the same coordinates as the gromacs-atom particle that it is attached to? (3) Now, according to the manual, thole polarization is intended to screen/damp/diminish coulomb interactions in shell molecular dynamics simulations. The wording suggests that thole polarization only impact gromacs-shell/gromacs-shell interactions (not gromacs-atom/gromacs-atom or gromacs-shell/gromacs-shell interactions). Am I correct? Each entry in the associated directive requires listing 4 numbers. The first two numbers specify a gromacs-atom and its attached gromacs-shell. The last two numbers specify a different gromacs-atom and its attached gromacs-shell. Given my assumption in (1) above, I would only need to include these for all 1-4, 1-5, etc. intramolecular gromacs-atom/gromacs-shell pairs. I assume that 1-4 interactions are still diminished by 0.5. Am I correct? Best, Eric On Tue, Jun 26, 2018 at 8:40 AM, Justin Lemkul wrote: > > > On 6/26/18 10:37 AM, Eric Smoll wrote: > >> Justin, >> >> In a previous thread on core/shell optimization ( >> https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users >> /2014-August/091101.html) >> you state: >> >> "No, I've made wholesale changes to the code in the way that energies and >> forces >> are computed. The current code has bugs. If you want a modified version >> (beta >> software, so caveat emptor), contact me off-list." >> >> Is this still true? Is core/shell optimization broken in GROMACS 2018.1? >> If so, would you mind sharing your Drude implementation? >> > > Check out git master and switch to the drude branch. Don't try to use > domain decomposition, only OpenMP for parallelization. Documentation is > unfortunately somewhat sparse but the code is reasonably commented. > > Note that everything I have done is for extended Lagrangian dynamics, I > haven't tested much with massless shells or use of [polarization] > directives. > > -Justin > > > Best, >> Eric >> >> On Mon, Jun 18, 2018 at 1:14 PM, Justin Lemkul wrote: >> >> >>> On 6/18/18 4:05 PM, Eric Smoll wrote: >>> >>> Justin, Thank you so much for the rapid and clear reply! Sorry to ask for a bit more clarification. The thole_polarization isn't in the manual at all. Is it structured the same way
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 6/26/18 10:37 AM, Eric Smoll wrote: Justin, In a previous thread on core/shell optimization ( https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2014-August/091101.html) you state: "No, I've made wholesale changes to the code in the way that energies and forces are computed. The current code has bugs. If you want a modified version (beta software, so caveat emptor), contact me off-list." Is this still true? Is core/shell optimization broken in GROMACS 2018.1? If so, would you mind sharing your Drude implementation? Check out git master and switch to the drude branch. Don't try to use domain decomposition, only OpenMP for parallelization. Documentation is unfortunately somewhat sparse but the code is reasonably commented. Note that everything I have done is for extended Lagrangian dynamics, I haven't tested much with massless shells or use of [polarization] directives. -Justin Best, Eric On Mon, Jun 18, 2018 at 1:14 PM, Justin Lemkul wrote: On 6/18/18 4:05 PM, Eric Smoll wrote: Justin, Thank you so much for the rapid and clear reply! Sorry to ask for a bit more clarification. The thole_polarization isn't in the manual at all. Is it structured the same way as the [ polarization ] directive in the manual: [ thole_polarization ] ; Atom i j type alpha 1 2 1 0.001 If I want Thole corrections, am I correct in assuming that I should list *all shells* in the system under this thole_polarization directive with (as you pointed out) "i" or "j" as the shell? If "i" is the shell, "j" is the core. If "j" is the shell, "i" is the core. You have to list everything explicitly, including the shielding factor, and between dipole pairs (Thole isn't between an atom and its shell, it's between neighboring dipoles). I honestly don't know what the format is; I completely re-wrote the Thole code for our Drude implementation (still not officially incorporated into a release due to DD issues, but we're close to a fix...) The code for "init_shell_flexcon" was very helpful. Thank you! nstcalcenergy must be set to 1. The code says that domain decomposition is not supported so multi-node MPI calculations are not allowed. I can still use an MPI-enabled GROMACS executable on a single node for shell MD, correct? Thread parallelization is still permitted, correct? Presumably you're limited to OpenMPI, but again I have no idea about this code. I've never actually used it. -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support /Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Justin, In a previous thread on core/shell optimization ( https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2014-August/091101.html) you state: "No, I've made wholesale changes to the code in the way that energies and forces are computed. The current code has bugs. If you want a modified version (beta software, so caveat emptor), contact me off-list." Is this still true? Is core/shell optimization broken in GROMACS 2018.1? If so, would you mind sharing your Drude implementation? Best, Eric On Mon, Jun 18, 2018 at 1:14 PM, Justin Lemkul wrote: > > > On 6/18/18 4:05 PM, Eric Smoll wrote: > >> Justin, >> >> Thank you so much for the rapid and clear reply! Sorry to ask for a bit >> more clarification. >> >> The thole_polarization isn't in the manual at all. Is it structured the >> same way as the [ polarization ] directive in the manual: >> >> [ thole_polarization ] >> ; Atom i j type alpha >> 1 2 1 0.001 >> >> If I want Thole corrections, am I correct in assuming that I should list >> *all shells* in the system under this thole_polarization directive with >> (as >> you pointed out) "i" or "j" as the shell? If "i" is the shell, "j" is the >> core. If "j" is the shell, "i" is the core. >> > > You have to list everything explicitly, including the shielding factor, > and between dipole pairs (Thole isn't between an atom and its shell, it's > between neighboring dipoles). I honestly don't know what the format is; I > completely re-wrote the Thole code for our Drude implementation (still not > officially incorporated into a release due to DD issues, but we're close to > a fix...) > > The code for "init_shell_flexcon" was very helpful. Thank you! >> nstcalcenergy must be set to 1. The code says that domain decomposition >> is not supported so multi-node MPI calculations are not allowed. I can >> still use an MPI-enabled GROMACS executable on a single node for shell MD, >> correct? Thread parallelization is still permitted, correct? >> > > Presumably you're limited to OpenMPI, but again I have no idea about this > code. I've never actually used it. > > > -Justin > > -- > == > > Justin A. Lemkul, Ph.D. > Assistant Professor > Virginia Tech Department of Biochemistry > > 303 Engel Hall > 340 West Campus Dr. > Blacksburg, VA 24061 > > jalem...@vt.edu | (540) 231-3129 > http://www.thelemkullab.com > > == > > -- > Gromacs Users mailing list > > * Please search the archive at http://www.gromacs.org/Support > /Mailing_Lists/GMX-Users_List before posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 6/18/18 4:05 PM, Eric Smoll wrote: Justin, Thank you so much for the rapid and clear reply! Sorry to ask for a bit more clarification. The thole_polarization isn't in the manual at all. Is it structured the same way as the [ polarization ] directive in the manual: [ thole_polarization ] ; Atom i j type alpha 1 2 1 0.001 If I want Thole corrections, am I correct in assuming that I should list *all shells* in the system under this thole_polarization directive with (as you pointed out) "i" or "j" as the shell? If "i" is the shell, "j" is the core. If "j" is the shell, "i" is the core. You have to list everything explicitly, including the shielding factor, and between dipole pairs (Thole isn't between an atom and its shell, it's between neighboring dipoles). I honestly don't know what the format is; I completely re-wrote the Thole code for our Drude implementation (still not officially incorporated into a release due to DD issues, but we're close to a fix...) The code for "init_shell_flexcon" was very helpful. Thank you! nstcalcenergy must be set to 1. The code says that domain decomposition is not supported so multi-node MPI calculations are not allowed. I can still use an MPI-enabled GROMACS executable on a single node for shell MD, correct? Thread parallelization is still permitted, correct? Presumably you're limited to OpenMPI, but again I have no idea about this code. I've never actually used it. -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
Justin, Thank you so much for the rapid and clear reply! Sorry to ask for a bit more clarification. The thole_polarization isn't in the manual at all. Is it structured the same way as the [ polarization ] directive in the manual: [ thole_polarization ] ; Atom i j type alpha 1 2 1 0.001 If I want Thole corrections, am I correct in assuming that I should list *all shells* in the system under this thole_polarization directive with (as you pointed out) "i" or "j" as the shell? If "i" is the shell, "j" is the core. If "j" is the shell, "i" is the core. The code for "init_shell_flexcon" was very helpful. Thank you! nstcalcenergy must be set to 1. The code says that domain decomposition is not supported so multi-node MPI calculations are not allowed. I can still use an MPI-enabled GROMACS executable on a single node for shell MD, correct? Thread parallelization is still permitted, correct? Thanks again, Justin. I hope my questions are clear and easy to answer. Best, Eric On Mon, Jun 18, 2018 at 12:07 PM, Justin Lemkul wrote: > > > On 6/18/18 3:00 PM, Eric Smoll wrote: > >> Hello GROMACS users, >> >> I am looking over the shell (Drude) model for polarization in GROMACS. >> There isn't much information available in the manual (probably because >> this >> feature is rarely used). I was hoping someone knowledgeable about >> polarizable simulations in GROMACS could help answer a few of my >> questions: >> >> (1) How exactly are Thole functions turned on? The 2018.1 manual does not >> specify. I am guessing they are hard-coded into all shell molecular >> dynamics simulations. Am I correct? >> > > They're activated with a [thole_polarization] directive. > > (2) To add a shell to a topology, I suspect I must specify a shell atomtype >> (setting the particle type to "S" for shell) and list the shells in the >> atoms directives. Using the "[ atomtypes ]" format in oplsaa.ff for a >> molecule with a res-name of "ABC", I assume the following will produce two >> shell particles with a charge of +1 and -1. Am I correct? >> >> [ atomtypes ] >> ; atom-type bond-type atomic-number mass charge particle-type sigma >> epsilon >> opls_ Sh 1 0 0 S 0 0 >> ; I am guessing that >> >> [ atoms ] >> ; atom-number atom-type res-number res-name atom-name charge-grp charge >> 1 opls_ 1 ABC Sh1 1 +1 >> 2 opls_ 1 ABC Sh2 2 -1 >> >> > > Looks right. > > (3) I suppose connectivity of each shell to its core is indicated with the >> "[ polarizability ]" directive. I am guessing "i" is the core atom and "j" >> is the shell particle. >> > > It actually doesn't matter. mdrun figures it out on the fly based on > particle type - see init_shell_flexcon in shellfc.cpp. > > (4) Since the code carries out an SCF calculation to relax the position of >> the shells at every timestep, I assume there is no need to specify shell >> masses or to thermostat the shell DOF. Does GROMACS omit shell particles >> from thermostats? The manual does not specify. >> > > Yes. Massless particles are not subjected to dynamical update routines. > > -Justin > > -- > == > > Justin A. Lemkul, Ph.D. > Assistant Professor > Virginia Tech Department of Biochemistry > > 303 Engel Hall > 340 West Campus Dr. > Blacksburg, VA 24061 > > jalem...@vt.edu | (540) 231-3129 > http://www.thelemkullab.com > > == > > -- > Gromacs Users mailing list > > * Please search the archive at http://www.gromacs.org/Support > /Mailing_Lists/GMX-Users_List before posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Shell (Drude) model for polarization in GROMACS
On 6/18/18 3:00 PM, Eric Smoll wrote: Hello GROMACS users, I am looking over the shell (Drude) model for polarization in GROMACS. There isn't much information available in the manual (probably because this feature is rarely used). I was hoping someone knowledgeable about polarizable simulations in GROMACS could help answer a few of my questions: (1) How exactly are Thole functions turned on? The 2018.1 manual does not specify. I am guessing they are hard-coded into all shell molecular dynamics simulations. Am I correct? They're activated with a [thole_polarization] directive. (2) To add a shell to a topology, I suspect I must specify a shell atomtype (setting the particle type to "S" for shell) and list the shells in the atoms directives. Using the "[ atomtypes ]" format in oplsaa.ff for a molecule with a res-name of "ABC", I assume the following will produce two shell particles with a charge of +1 and -1. Am I correct? [ atomtypes ] ; atom-type bond-type atomic-number mass charge particle-type sigma epsilon opls_ Sh 1 0 0 S 0 0 ; I am guessing that [ atoms ] ; atom-number atom-type res-number res-name atom-name charge-grp charge 1 opls_ 1 ABC Sh1 1 +1 2 opls_ 1 ABC Sh2 2 -1 Looks right. (3) I suppose connectivity of each shell to its core is indicated with the "[ polarizability ]" directive. I am guessing "i" is the core atom and "j" is the shell particle. It actually doesn't matter. mdrun figures it out on the fly based on particle type - see init_shell_flexcon in shellfc.cpp. (4) Since the code carries out an SCF calculation to relax the position of the shells at every timestep, I assume there is no need to specify shell masses or to thermostat the shell DOF. Does GROMACS omit shell particles from thermostats? The manual does not specify. Yes. Massless particles are not subjected to dynamical update routines. -Justin -- == Justin A. Lemkul, Ph.D. Assistant Professor Virginia Tech Department of Biochemistry 303 Engel Hall 340 West Campus Dr. Blacksburg, VA 24061 jalem...@vt.edu | (540) 231-3129 http://www.thelemkullab.com == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
