How to extend a force field? (Don't use inheritance)

[scmartin]

I have a force field which is an extension of another force field. The new force field is in it's own .lt file. Both force fields use 'Data XXX By Type', but I believe I found a bug. As an example if my new forcefield lt file is something like

import "forcefield1.lt" 

FF2 inherits FF1 {
write_once('Data Masses') {
...
}

write_once('Data Bonds By Type'){
@bond:ab @atom:a @atom:b
@bond:bc @atom:b @atom:c
}

write_once('Data Angles By Type') {
@angle:A @atom:a @atom:b @atom:c @bond:* @bond:*
...
}

}

then moltemplate fails to find the angles. After messing around with it for a bit, I found that it seems to be checking in the first forcefield file. The output from Moltemplate for generating 3-body interactions by atom/bond type looks similar to

checking @/angle:FF2/A type requirements:
['@atom:FF2/a', '@atom:FF2/b', '@atom:FF2/c', '@bond:FF1/ab '@bond:FF1/bc']

If I remove the @bond:* variables from the data angles by type lines, it works correctly.

[hothello]

Can you provide a complete input deck that reproduces this issue?

As a general comment: it only makes sense to use the bond terms constraints if you need to distinguish a specific triplet using something other than the atom types. But I agree that the behaviour of nbody_by_type.py with the @bond:* @bond:* terms added should match that without them.

[jewettaij]

See the end of this discussion for a workaround ("workaround 2")

I posted several long comments.

Feel free to skip to the end for my final recommendation. (Don't use inheritance. See "workaround 2".)

It's actually not a bug. But your question reveals that I need to add more documentation. I'll do that today.

[jewettaij]

I will take a look why the "@bond:*" gets mapped to "@/bond:FF1/*" not "@/bond:FF2/*"

This happens because when a force field inherits from another force-field there is a fundamental ambiguity regarding how wildcards ("*" or "?") are interpreted:

• Do you want to match with all the possible bond types defined in both FF1 and FF2?
• Or just the new ones in FF2? (I don't think that's what you want.)

Workaround 1 (not recommended)

You can clarify this using:

    @angle:A_B_C @atom:a @atom:b @atom:c @bond:/FF*/* @bond:/FF*/*

But you will still probably have issues.

The same problem occurs for "@Atoms" as well as "@bonds" (and any other "@" variable).
If you want your "Angle By Type" section to be able to refer to atom types defined in either FF1 or FF2, use:

    @angle:A_B_X @atom:/FF*/a @atom:/FF*/b @atom:/FF*/* @bond:/FF*/* @bond:/FF*/*

And if you defined a new atom type in FF2 and you want to match it using *, you must clarify this the same way:

    @angle:A_B_X @atom:/FF*/a @atom:/FF*/b @atom:/FF*/* @bond:/FF*/* @bond:/FF*/*

This method is ugly and treacherous. (It would be easy to accidentally cause moltemplate to think that you are trying to define a new atom type of type "@atom:a" in FF2, even if "@atom:a" already exists in FF1.)

I have a better solution below. ("Workaround 2")

[jewettaij]

Workaround 2 (recommended)

I suggest you implement your forcefield extension this way instead:

# --- new file "ff2.lt" ---

import "ff1.lt"

FF1 {
  # ADD YOUR NEW ATOM TYPES AND BONDED INTERACTIONS HERE
  # Note: This will not overwrite FF1.
  #        (Instead, it appends to the existing FF1 defined earlier.)
  write_once("Data Angles By Type") {
    @angle:A_B_C @atom:a @atom:b @atom:c @bond:* @bond:*
  }
}

Molecules that use the new force field version would import the new "ff2.lt" file.

import "ff2.lt"

MoleculeV2 inherits FF1 {
   # Here you can refer to @atoms, @bonds, etc... defined in either file:
   # - "ff1.lt"
   # - "ff2.lt"
}

That is how the LOPLS extension was implemented.
(See the force-field file at loplsaa.lt and example that demonstrates how to use it: alkane_chain.)

Thanks for posting this question, by the way.
It is not obvious at all that this is the correct solution.

[jewettaij]

Documentation updated in a few different places. Hopefully people will notice it in the future.
I don't want to alienate people like you who are developing new force fields.

(And if you want to share the new force-field you've created, please do!)

I will leave this question open in case other people run into this problem.