Is there a particular reason why NWChem does not support QM/MM dynamics with all atoms?

[hjjvandam]

NWChem has QM/MM capabilities that allow one to optimize the positions of all atoms (QM, Link atoms, MM Solute, and Solvent). When running a dynamics simulation with the QM/MM module the time evolution only of region 1 (typically the QM or QMLink region) is evaluated. All other atoms are kept fixed. Given that the code can calculate forces on all atoms, is there a particular reason for restricting the dynamics to just the QM atoms?

[ebylaska]

Probably because it’s expensive and one needs to implement a far to local field expansion to make the one electron operator cost more reasonable. I’m sure I explained this to you many years ago. Knowing myself probably multiple times. This is also described in a recent NWChemex writeUp to interface it with exaalt. Note the aimd/mm method in the plane wave code runs dynamics for both the QM and mm atoms very fast.

---

From: Hubertus van Dam notifications@github.com Sent: Saturday, October 31, 2020 9:14:01 PM To: nwchemgit/nwchem nwchem@noreply.github.com Cc: Subscribed subscribed@noreply.github.com Subject: [nwchemgit/nwchem] Is there a particular reason why NWChem does not support QM/MM dynamics with all atoms? (#264)

NWChem has QM/MM capabilities that allow one to optimize the positions of all atoms (QM, Link atoms, MM Solute, and Solvent). When running a dynamics simulation with the QM/MM module the time evolution only of region 1 (typically the QM or QMLink region) is evaluated. All other atoms are kept fixed. Given that the code can calculate forces on all atoms, is there a particular reason for restricting the dynamics to just the QM atoms?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/nwchemgit/nwchem/issues/264, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AATEVFOUZXDWXHSTDQDIKJ3SNTOATANCNFSM4TGI4SKQ.

[hjjvandam]

Ok, that makes sense. I was wondering if doing the forces per region meant that some terms were missing. So I should try the planewave AIMD/MM functionality instead. I will have a go at that.

[hjjvandam]

So I am trying to turn the qmmm_opt0 QA test into an AIMD/MM calculation as I am interested in doing QM/MM AIMD on biomolecular systems. The classical MD capabilities in NWChem obviously recognize the usual protein residues without problem. I can get NWChem to produce a PDB file using "prepare" containing a solvated version of the protein. For the AIMD calculation, however, I have to explicitly specify the force field parameters. This requires that the atom naming convention I use allows me to uniquely define the force field parameters, but I should also allow the PSPW module to recognize what the element of the atom is. The atom names generated by the classical MD module may start with a number (e.g. 2HB for residue ALA) which would not be recognized as a Hydrogen atom by the PSPW I think. Other tricky atoms are alpha-Carbons which are labeled "CA" by the classical MD module, but they are Carbons and not Calcium atoms. So I could translate the alpha-Carbon of Alanine (CA ALA) to something like C_A_ALA, so the PSPW module recognizes it as a Carbon but also maintaining the information that this is a particular Carbon atom that is part of an Alanine residue. I could subsequently look the force field parameters up from the Amber force fields as implemented in the classical MD module and insert them into the input. Is this naming convention sensible? I.e. rename 2HB of Alanine to H_B2_ALA, and CA of Alanine to C_A_ALA, and then specify force field parameters accordingly? I guess this is flexible enough to avoid clashes with more general applications that involve, for example, inorganic structures such as a variety of materials (including zeolites, oxides, etc.). Would you agree?

[Straatsma]

Hi Huub,

The atom names in the MD module have the following convention:

The first two characters are the element name, except when a hydrogen in cases where more than one hydrogen is bonded to the same atom. An alpha carbon would start with a space: “ CA”, while a Calcium atom would not: “Ca “ or “CA “. A hydrogen is “ H” if only hydrogen to an atom, “2H”, “3H”, “4H” as needed for groups with multiple hydrogens to the same atom.

This is the official PDB naming convention (which is strictly followed in NWChem, while not in other MD codes).

Tjerk

Dr. T. P. Straatsma

Distinguished Research Scientist, National Center for Computational Sciences, ORNL Adjunct Professor, Department of Chemistry, University of Alabama

Oak Ridge National Laboratory P.O.Box 2008, MS 6008 Oak Ridge, TN 37831-6373

Bldg 5600, A325 Tel. (865) 241-5864 Adm. (865) 574-0972 Fax. (865) 241-8484

From: Hubertus van Dam notifications@github.com Sent: Tuesday, November 03, 2020 23:55 To: nwchemgit/nwchem nwchem@noreply.github.com Cc: Subscribed subscribed@noreply.github.com Subject: [EXTERNAL] Re: [nwchemgit/nwchem] Is there a particular reason why NWChem does not support QM/MM dynamics with all atoms? (#264)

So I am trying to turn the qmmm_opt0 QA test into an AIMD/MM calculation as I am interested in doing QM/MM AIMD on biomolecular systems. The classical MD capabilities in NWChem obviously recognize the usual protein residues without problem. I can get NWChem to produce a PDB file using "prepare" containing a solvated version of the protein. For the AIMD calculation, however, I have to explicitly specify the force field parameters. This requires that the atom naming convention I use allows me to uniquely define the force field parameters, but I should also allow the PSPW module to recognize what the element of the atom is. The atom names generated by the classical MD module may start with a number (e.g. 2HB for residue ALA) which would not be recognized as a Hydrogen atom by the PSPW I think. Other tricky atoms are alpha-Carbons which are labeled "CA" by the classical MD module, but they are Carbons and not Calcium atoms. So I could translate the alpha-Carbon of Alanine (CA ALA) to something like C_A_ALA, so the PSPW module recognizes it as a Carbon but also maintaining the information that this is a particular Carbon atom that is part of an Alanine residue. I could subsequently look the force field parameters up from the Amber force fields as implemented in the classical MD module and insert them into the input. Is this naming convention sensible? I.e. rename 2HB of Alanine to H_B2_ALA, and CA of Alanine to C_A_ALA, and then specify force field parameters accordingly? I guess this is flexible enough to avoid clashes with more general applications that involve, for example, inorganic structures such as a variety of materials (including zeolites, oxides, etc.). Would you agree?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/nwchemgit/nwchem/issues/264#issuecomment-721511279, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEDIDLCCMDWPKNXWHRII62DSODNDDANCNFSM4TGI4SKQ.

[hjjvandam]

Upon looking closer it is a bit different. It seems like the QMMM block in the PSPW code expects a specification of all the force field parameters. The atoms are identified by their rank in the geometry. So basically what needs to happen is that in the geometry we need to have the element for each atom and its coordinates which is easily extracted from a PDB file. For the force field parameters we need to basically extract them from the topology file that can be generated with the prepare module. Doing this is not difficult, just tedious. I have started putting a script together that can do this. It will take a PDB file and produce an input file for the PSPW QMMM capability. Obviously this script will be useful only for biomolecular systems.