Reading/using virtual sites automatically

[SofiaBariami]

Hello, I am trying to modify the code in Sire to automatically read and recognise virtual sites. So far I am doing this manually, by inserting the parameters of the virtual sites in the code of openmmfrenergyst.cpp. What we would like to do is generate a MORPH.pert file with BSS that has a specific section for the virtual sites like the following and modify Sire source code to understand the information.

version 1
molecule MOL
    virtual-site
        name           X4
        initial_type   X4
        final_type     X4
        initial_LJ     3.39680 0.10060
        final_LJ       3.37680 0.09821
        initial_charge -0.39020
        final_charge   -0.23020
    endvirtual-site
endmolecule

Please correct me if I'm wrong, but from my understanding, we should create a new perturbations template at corelib/src/libs/SireIO/perturbationslibrary.cpp similar to the one that is used for atoms. Then we need to define the parameters of the virtual sites :

 // virtual sites template: 
        if ( words[0] == "virtual-site")
        {
            vsname = "";
            chargeinit_vs = 0 * mod_electron;
            chargefinal_vs = 0 * mod_electron;
            sigmainit_vs = 0 * angstrom;
            sigmafinal_vs = 0 * angstrom;
            epsiloninit_vs = 0 * kcal_per_mol;
            epsilonfinal_vs = 0 * kcal_per_mol;
            inVsite = true;
            continue;
        }

and finally set the parameters:

new_templates[current].setInitCharge_vs(name, chargeinit_vs);
new_templates[current].setFinalCharge_vs(name, chargefinal_vs);
LJParameter ljinit_vs = LJParameter(sigmainit_vs, epsiloninit_vs);
new_templates[current].setInitLJ_vs(name, ljinit_vs);
LJParameter ljfinal_vs = LJParameter(sigmafinal, epsilonfinal);
new_templates[current].setFinalLJ_vs(name, ljfinal_vs);

We also need to set the v-sites properties names and property maps. For example:

PropertyName charge_property_vs = PropertyName("charge");
editmol.setProperty( initial_charge_property_vs, editmol.property(charge_property_vs) );

// virtual sites property map
  PropertyMap charge_map_vs;
  charge_map_vs.set("initial_charge_vs", initial_charge_property_vs);
  charge_map_vs.set("charge_vs", charge_property_vs);
  charge_map_vs.set("final_charge_vs", final_charge_property_vs);
  perturbations.append( ChargePerturbation(charge_map_vs) );

void PerturbationsTemplate::setInitCharge_vs(const QString &atomname, const SireUnits::Dimension::Charge &atomcharge)
{
  initcharges_vs.insert(atomname, atomcharge);
}

void PerturbationsTemplate::setFinalCharge_vs(const QString &atomname, const SireUnits::Dimension::Charge &atomcharge)
{
  finalcharges_vs.insert(atomname, atomcharge);
}

After this is done, we need to modify the source code at openmmfrenergyst.cpp in order for Sire to understand the parameters as charges and LJ. For example, the corresponding part of the code for atoms would be the following:

if (molecule.hasProperty("perturbations"))
        {
            if (Debug)
                qDebug() << "Molecule Perturbed number = " << i;
            AtomCharges atomcharges_start = molecule.property("initial_charge").asA<AtomCharges>();
            AtomCharges atomcharges_final = molecule.property("final_charge").asA<AtomCharges>();

            start_charges = atomcharges_start.toVector();
            final_charges = atomcharges_final.toVector();

            AtomLJs atomvdws_start = molecule.property("initial_LJ").asA<AtomLJs>();
            AtomLJs atomvdws_final = molecule.property("final_LJ").asA<AtomLJs>();
            start_LJs = atomvdws_start.toVector();
            final_LJs = atomvdws_final.toVector();
        }

For virtual sites, we cannot assign the parameters as “AtomCharges” or “AtomLJs” because these properties are just used to describe atoms. My question is if we should define new “V-siteCharges” and “V-siteLJs” and a class V-siteProperty to include the new parameters? To do that do we also need to make two new files “SireMM/v-sitesljs.h/cpp” and “SireMM/v-sitescharges.h/cpp” that define the property, similar to the existing ones for atoms and also define a new v-siteProperty similar to AtomProperty<SireMM::LJParameter> for example? If so how can we do that? I'd really appreciate any help on that issue. Thank you very much!

[chryswoods]

You could do that. But, what are the reasons not to use AtomCharges or AtomLJs? In the Monte Carlo code a virtual site is just a dummy atom (zero charge and zero LJ parameters). Would it be easier to adapt the perturbation code to allow virtual sites to use zero charge and LJ parameters, rather than add a new v-siteProperty?

[jmichel80]

hi @SofiaBariami the v-sites do not map cleanly on atom properties because so far we are treating v-sites as a molecules properties that provides a list of virtual sites to apply to a molecule. This was easier to code via a property dictionary attached to a molecule to get going. @chryswoods suggestion could work if the code reading in v-sites is updated to edit Sire molecules to add new atoms. This means changing a bit the code you have been writing, but would minimise additional changes needed in Sire. Could you review this alternative route, post a ''walkthrough'' of the changes required and ask @jmichel80 @lohedges @chryswoods to comment on the solution ?

[SofiaBariami]

Thank you @chryswoods and @jmichel80 for the suggestions. Here are some thoughts about the implementation:

First of all, if we set the virtual site as a zero-mass atom, does this mean that we do not need to call setVirtualSite() at all? We can generate pdb files that also contain the coordinates of the virtual sites. This can be used as an alternative to getting the local coordinates for the LocalCoordinatesSite. Since the vSite is bonded to the parent atom, we probably won't have an issue defining its position during an MD.

Charges and LJ parameters can still be read from the property dictionary and added to nonbonded_openmm and custom_forcefield with addParticle. I just need to append the list of the atoms and add a particle with zero mass. Looking at the OpenMM integrators documentation here, we read that "CustomIntegrator ignores any particle whose mass is 0. It is skipped when doing per-DOF computations, and is not included when computing sums over degrees of freedom". However, at line 3091 of openmmfrenergyst.cpp, we see that there is a check for particles with zero mass. Do you think that this is going to be an issue?