You may have a forcefield expression that contains a free term that is not dependent on a forcefield component.

[julien]

Hi Chris,

I am working on a script to process a MD trajectory using Sire. To avoid a large performance hit, I am not recreating a system from scratch for every frame in a DCD trajectory, rather I update the system coordinates.

The script appears to work as I would like to, and gives me the same energies as an unoptimised old script.  However I now get this curious error message for every frame, apart from the 1st frame.

"You may have a forcefield expression that contains a free term (E_{gridwaterff}^{CLJ} + E_{gridwater:otherff}^{CLJ}) that is not dependent on a forcefield component. This is because each term must have dimensions of energy, and the addition of a constant dimensionless value is not dimensionally correct. Please check that the remainder is zero."
"You may have a forcefield expression that contains a free term (E_{gridwaterff}^{CLJ} + E_{gridwater:otherff}^{CLJ}) that is not dependent on a forcefield component. This is because each term must have dimensions of energy, and the addition of a constant dimensionless value is not dimensionally correct. Please check that the remainder is zero."
"You may have a forcefield expression that contains a free term (E_{gridwater:otherff}^{CLJ}) that is not dependent on a forcefield component. This is because each term must have dimensions of energy, and the addition of a constant dimensionless value is not dimensionally correct. Please check that the remainder is zero."
"You may have a forcefield expression that contains a free term (E_{gridwater:otherff}^{CLJ}) that is not dependent on a forcefield component. This is because each term must have dimensions of energy, and the addition of a constant dimensionless value is not dimensionally correct. Please check that the remainder is zero."

This debug statement comes from void FFSymbolExpression::expandInTermsOf(const QSet<Symbol> &ffsymbols) in forcefields.cpp

The subroutine that initialises each snapshot is a bit complex, but it roughly follows that structure

* Go to the desired DCD frame. Buffer atom coordinates of each molecule.
* Change the coordinate of each molecule. Also assign the molecule to a group new_gridwater (water molecule within a 3D grid) or group new_other (everything else)
* Then in details

    system.update(changedmols)
    system.removeAllForceFields()
    system.remove(MGName("gridwater"))
    system.remove(MGName("other"))

    system.add(new_gridwater)
    system.add(new_other)

(the new groups have the same name as the old ones that were in the system "gridwater" and "other")

* Then call a setupForceFields() subroutine where new forcefields involving the groups gridwater and other are created using the new molecule groups.

* The first two debug messages are triggered by the call to system.removeAllForceFields().
* The other two debug messages are triggered by the code below in the subroutine setupForceFields()

    for forcefield in forcefields:
        system.add(forcefield)

Now my questions :-)

- Am I doing something bad?

- Is there a way to cleanly remove groups/forcefields from a system without triggering these messages?

- If not, would you suggest that I convert the debug statement into a user warning, and catch the warning in my script, since I still seem to
get the desired behaviour?

Best wishes,


Julien

[Christopher Woods]

  Hi Julien,

The debug message occurs when the code sees that an energy expression contains a constant term, e.g. 

E_total = E_clj + lambda

lambda in this case is the constant term, and is dimensionally incorrect (as it does not have units of energy). I've not made this an error as there may be some (admittedly strange) cases where the user may want to add a constant to the total energy. In your case the error occurs because, while removing and adding forcefields, you are temporarily changing the energy expression to have constant terms (E_{gridwaterff}^{CLJ} and E_{gridwater:otherff}^{CLJ}). When the system doesn't contain the "gridwaterff" or "gridwater:other" forcefields, Sire temporarily creates dummy symbols for them so that the energy expression will still be valid. These dummy symbols are constants, rather than energy symbols, and so they trigger the debug warning (once while removing the forcefields, and once while adding them). The messages are harmless in your case, so will not affect the results. You can comment out the warning statement in the code if the messages are annoying.

(1) You are not doing anything bad. The error comes because Sire cannot yet update forcefields in a system, and so you have to remove and then re-add them. Sire could be modified so that you can update forcefields like you can update molecules and molecule groups, which may be slightly faster than deleting and re-adding forcefields.

(2) You can cleanly remove forcefields from a system only when the energy expressions in the system don't refer to energy components from those forcefields. The simplest way to do this is to remove all energy expressions, then remove forcefields, and then add back energy expressions that don't involve those forcefields.

(3) Yes, perhaps this error should be a user warning, as anyone who knows what they are doing should be able to catch the warning if they truly want to have constant terms in their energy expressions.

  Cheers,

  Christopher