Virtual Sites Implementation

[SofiaBariami]

I am trying to create a Virtual Site (VS) in Sire with parameters that are stored in a "virtual-sites" property of a molecule. I am testing my implementation by running single point energy calculations. However, the energy that I get when running the calculation with the molecule containing VS is the same as the one that I get without the VS. That means that Sire does not recognise the implementation of the VS. To implement the VS, I followed the same process that was followed for the Restraints and TIP4P, using the OpenMM function LocalCoordinatesSite, instead of the ThreeParticleAverageSite used for TIP4P. The files that I changed are the the openmmfrenergyst.* and the openmmmdintegrator.*.
The main part of the code that I added is the following:

    OpenMM::VirtualSite * vsite = NULL;

    if (VirtualSite_flag == true)
    {

           int nVSites;
           int vsIndex;
           int atom1;
           int atom2;
           int atom3;
           double p1;
           double p2;
           double p3;
           double wo1;
           double wo2;
           double wo3;
           double wx1;
           double wx2;
           double wx3;
           double wy1;
           double wy2;
           double wy3;
           double charge;
           double sigma;
           double epsilon;

        int openmmindex;
        vsite = new OpenMM::LocalCoordinatesSite(atom1, atom2, atom3, OpenMM::Vec3(wo1, wo2, wo3), OpenMM::Vec3(wx1, wx2, wx3), OpenMM::Vec3(wy1, wy2, wy3), OpenMM::Vec3(p1,p2,p3));
        nonbond_openmm->addParticle(charge, sigma * OpenMM::NmPerAngstrom, epsilon * OpenMM::KJPerKcal);
        system_openmm->setVirtualSite(openmmindex, vsite);
    }

        if (VirtualSite_flag == true)
        {
            bool hasVirtualSites = molecule.hasProperty("virtual-sites");

            if (hasVirtualSites)
            {
                Properties virtualSites = molecule.property("virtual-sites").asA<Properties>();

                int nvirtualsites = virtualSites.property(QString("nvirtualsites")).asA<VariantProperty>().toInt();

                if (Debug)
                    qDebug() << "nvirtualsites = " << nvirtualsites;
                for (int i = 0; i < nvirtualsites; i++)
                {

                    int nVSites = virtualSites.property(QString("nvirtualsites(%1)").arg(i)).asA<VariantProperty>().toInt();
                    int vsIndex = virtualSites.property(QString("index(%1)").arg(i)).asA<VariantProperty>().toInt();
                    int atom1 = virtualSites.property(QString("atom1(%1)").arg(i)).asA<VariantProperty>().toInt();
                    int atom2 = virtualSites.property(QString("atom2(%1)").arg(i)).asA<VariantProperty>().toInt();
                    int atom3 = virtualSites.property(QString("atom3(%1)").arg(i)).asA<VariantProperty>().toInt();
                    double p1 = virtualSites.property(QString("p1(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double p2 = virtualSites.property(QString("p2(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double p3 = virtualSites.property(QString("p3(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wo1 = virtualSites.property(QString("wo1(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wo2 = virtualSites.property(QString("wo2(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wo3 = virtualSites.property(QString("wo3(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wx1 = virtualSites.property(QString("wx1(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wx2 = virtualSites.property(QString("wx2(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wx3 = virtualSites.property(QString("wx3(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wy1 = virtualSites.property(QString("wy1(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wy2 = virtualSites.property(QString("wy2(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double wy3 = virtualSites.property(QString("wy3(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double charge = virtualSites.property(QString("charge(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double sigma = virtualSites.property(QString("sigma(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    double epsilon = virtualSites.property(QString("epsilon(%1)").arg(i)).asA<VariantProperty>().toDouble();
                    //QString name = virtualSites.property(QString("name(%1)").arg(i)).asA<VariantProperty>().toString();
                    //QString type = virtualSites.property(QString("type(%1)").arg(i)).asA<VariantProperty>().toString();
                    int openmmindex = AtomNumToOpenMMIndex[vsIndex];

                    if (Debug)
                    {
                        qDebug() << "atom1 " << atom1 << " p1 " << p1 << " wx1 " << wx1 << " wy1 " << wy1 << " sigma " << sigma << " epsilon " << epsilon;
                    }

                    int vSitedim = 20;
                    std::vector<double> vSites_params(vSitedim);

                    vSites_params[0] = nVSites;
                    vSites_params[1] = vsIndex;
                    vSites_params[2] = atom1;
                    vSites_params[3] = atom2;
                    vSites_params[4] = atom3;
                    vSites_params[5] = p1;
                    vSites_params[6] = p2;
                    vSites_params[7] = p3;
                    vSites_params[8] = wo1;
                    vSites_params[9] = wo2;
                    vSites_params[10] = wo3;
                    vSites_params[11] = wx1;
                    vSites_params[12] = wx2;
                    vSites_params[13] = wx3;
                    vSites_params[14] = wy1;
                    vSites_params[15] = wy2;
                    vSites_params[16] = wy3;
                    vSites_params[17] = charge;
                    vSites_params[18] = sigma * OpenMM::NmPerAngstrom;
                    vSites_params[19] = epsilon * (OpenMM::KJPerKcal);
                    //vSites_params[20] = name;
                    //vSites_params[21] = type;

                  OpenMM::LocalCoordinatesSite * vsite = new OpenMM::LocalCoordinatesSite(atom1, atom2, atom3, OpenMM::Vec3(wo1, wo2, wo3), OpenMM::Vec3(wx1, wx2, wx3), OpenMM::Vec3(wy1, wy2, wy3), OpenMM::Vec3(p1,p2,p3) );
                nonbond_openmm->addParticle(charge, sigma, epsilon);
                system_openmm->setVirtualSite(openmmindex, vsite);
                        // virtualSites_openmm->addParticle(vSites_params);
                }
            }
        }//end of virtual sites flag

(You can find all the files at the qube_combRules-feature_xml branch)
My question is if this is the correct way to set the VS:

nonbond_openmm->addParticle(charge, sigma * OpenMM::NmPerAngstrom, epsilon * OpenMM::KJPerKcal);
system_openmm->setVirtualSite(openmmindex, vsite);

Also, the code compiles, but there are warnings like the following:

/users/sofia/Software/Sire_xml/Sire/corelib/src/libs/SireMove/openmmfrenergyst.cpp:1312:210: warning: 'atom1' may be used uninitialized in this function -Wmaybe-uninitialized]
 LocalCoordinatesSite(atom1, atom2, atom3, OpenMM::Vec3(wo1, wo2, wo3), OpenMM::Vec3(wx1, wx2, wx3), OpenMM::Vec3(wy1, wy2, wy3), OpenMM::Vec3(p1,p2,p3));
  
/users/sofia/Software/Sire_xml/Sire/corelib/src/libs/SireMove/openmmfrenergyst.cpp:1286:36: warning: unused variable 'vsite' [-Wunused-variable]
     OpenMM::LocalCoordinatesSite * vsite = NULL;
                                    ^~~~~

Should I worry about them?

Finally, the script that I am using, to get the energies, along with the input files can be found here and here.
I would gladly appreciate some comments/ thoughts.
Thanks a lot!

[jmichel80]

Hi @SofiaBariami as I told you on Monday Sire doesn't support virtual sites for energy calculations so system.energy() is not going to work. If you refer to the notes from our meeting on Monday I told to use integrator.getPotentialEnergy() to get the single point energy calculated by OpenMM. This means you need to create an OpenMM Move to initialise the integrator.

So your script will need code like

moves = setupMoves(system, ranseed, gpu.val)
integrator = move0.integrator()
integrator.getPotentialEnergy() 

See how this is done in https://github.com/michellab/Sire/blob/devel/wrapper/Tools/OpenMMMD.py
for inspiration

It looks like your C++ warning are because you declared variables inside a if statement but you assign values in another block, the compiler is unable to tell if this could lead to runtime errors. This is best avoided by rethinking where and when you declare your variables. Focus on this once this code is actually executed by your python scripts.

[jmichel80]

hi @SofiaBariami have you been able to progress with this issue ? keep me posted.

[SofiaBariami]

Hi Julien,
Yes, I managed to do it. The problem is that I am now getting a Segmentation Fault every time I am trying to calculate energies of molecules with virtual sites. The script works fine with molecules without any VS. I am looking at this problem right now

[SofiaBariami]

The problem was that I had also written code for openmmmdintegrator.cpp, that wasn't needed. After I commended it out of the file, the implementation of the VS worked. The energies are different from the ones reported from the Cole group, so I am now working to fix that.

[jmichel80]

Great ! Check out the thread here openmm/openmm#2045 , it could be with how you build excluded pairs for non-bonded interactions Julien

…

-------------------------------------------------------------- Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/ ------------------------------------------------------------- On Fri, Nov 1, 2019 at 3:49 PM Sofia Bariami <[email protected]<mailto:[email protected]>> wrote: The problem was that I had also written code for openmmmdintegrator.cpp, that wasn't needed. After I commended it out of the file, the implementation of the VS worked. The energies are different from the ones reported from the Cole group, so I am now working to fix that. — You are receiving this because you commented. Reply to this email directly, view it on GitHub<#299?email_source=notifications&email_token=ACZN3ZFYUEPITDMLNMVCBOTQRRFX7A5CNFSM4JEXQT52YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEC3KAEA#issuecomment-548839440>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACZN3ZCWPKYSPJ6T6YCVGQLQRRFX7ANCNFSM4JEXQT5Q>. The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.