Sire icon indicating copy to clipboard operation
Sire copied to clipboard

Remove water topology naming restriction required by SOMD

Open lohedges opened this issue 6 years ago • 0 comments

Currently SOMD requires an exact AMBER water topology naming convention. This is down to the way that it applies constraints, e.g. from openmmfrenergyst.cpp:

  if (flag_constraint == NONE)
            {
                //JM 10/16 If constraint water flag is on and if molecule is a water molecule then apply constraint
                if (flag_constraint_water and molfirstresname == ResName("WAT"))
                    system_openmm->addConstraint(idx0, idx1, r0 * OpenMM::NmPerAngstrom);
                else
                    bondStretch_openmm->addBond(idx0, idx1, r0 * OpenMM::NmPerAngstrom, k * 2.0 * OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm);

              //cout << "\nBOND ADDED TO "<< atom0.toStdString() << " AND " << atom1.toStdString() << "\n";
            }
            else if (flag_constraint == ALLBONDS || flag_constraint == HANGLES)
            {
                system_openmm->addConstraint(idx0, idx1, r0 * OpenMM::NmPerAngstrom);
                //cout << "\nALLBONDS or HANGLES ADDED BOND CONSTRAINT TO " << atom0.toStdString() << " AND " << atom1.toStdString() << "\n";
            }
            else if (flag_constraint == HBONDS)
            {

                if ((atom0[6] == 'H') || (atom1[6] == 'H'))
                {
                    system_openmm->addConstraint(idx0, idx1, r0 * OpenMM::NmPerAngstrom);
                    //cout << "\nHBONDS ADDED BOND CONSTRAINT TO " << atom0.toStdString() << " AND " << atom1.toStdString() << "\n";
                }
                else
                {
                    bondStretch_openmm->addBond(idx0, idx1, r0 * OpenMM::NmPerAngstrom, k * 2.0 * OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm);
                    //cout << "\nHBONDS ADDED BOND TO " << atom0.toStdString() << " AND " << atom1.toStdString() << "\n";
                }
            }

Note that the water molecules require ResName("WAT") and the atoms are identified by a character in the string representing their name, e.g. if ((atom0[6] == 'H') || (atom1[6] == 'H')) (not sure why the element property isn't checked so that any name can be used).

This restriction causes problems for BioSimSpace, see, e.g. here, since we need to know the exact water topology in order to convert it to the correct format required by SOMD. Removing this naming restriction should hopefully be easy since there is now built in functionality within Sire to be able to search for water molecules within a system, e.g. system.search("water") returns a Sire.Mol._Mol.SelectResult object containing all of the water molecules in the system.

lohedges avatar Aug 13 '19 08:08 lohedges