fix: correct usage of nonbond energy for close contacts

[DaniBodor]

Previously, a distinct set of parameters was used to calculate Van der Waals energy, depending on whether atoms are part of the same residue or not. Electrostatic energy was calculated in the same way for any atom pair.

In the current PR, this is changed such that:

• VDW and electrostatic energy are set to 0 for any atom pairs that are within 3.6 A distance of each other (as a proxy for covalent pairs and 1-3 pairs)
• VDW energy uses low-energy parameters to calculate energy for atom pairs between 3.6 A and 4.2 A distance (proxy for 1-4 pairs)
  See feat: redefine potential energy terms between atom pairs, depending on bond distance #357 for more details

Tests were adjusted and added according to above rules.

Also, nomenclature has been improved for:

• nomenclature for vanderwaals and electrostatic energy has been unified
• nomenclature for inter_ and intra_ parameters adjusted to _main and _14 parameters in the parser. Not sure I am totally happy with this nomenclature, so open to suggestions.

[DaniBodor]

@cbaakman , do you think we need to create a test to check that the correct parameters are used?

Do you have an idea how to do this in a simple way?

[cbaakman]

@cbaakman , do you think we need to create a test to check that the correct parameters are used?

Do you have an idea how to do this in a simple way?

We could at least check the output of the function _intra_partners, use a well known pdb structure with well known binding partners.

[cbaakman]

I see that numpy is still used here. I wonder whether we could use torch instead.

[DaniBodor]

@cbaakman , do you think we need to create a test to check that the correct parameters are used?
Do you have an idea how to do this in a simple way?

We could at least check the output of the function _intra_partners, use a well known pdb structure with well known binding partners.

I actually since noticed that in PR #348 , you added tests for exactly this (although the function used there is not part of the PR, maybe you forgot to push that file?). I am guessing you looked at the structure by eye and found residues at different covalent bonding distances to test individually. I can indeed use those same residues to test the output of _intra_partners as you suggest.

It would be even nicer to test whether the output of _get_vdw_energy is correct, but that would be slightly more complex. I'll take a look whether I can make that work as well.

[DaniBodor]

I see that numpy is still used here. I wonder whether we could use torch instead.

What would the advantage be? Speed?

[cbaakman]

What would the advantage be? Speed?

Yes, we could use the GPU. As I understand, numpy is always CPU.

[DaniBodor]

@cbaakman , can you please re-review this, as it is now a completely different PR than it was the first time you looked at it.

[DaniBodor]

In case you review this while the builds are still failing, please ignore that. It is a conflict between newest version of pytorch and pytorch geometric, and is being dealt with in a separate PR (#386)

[LilySnow]

For the test, we should compare with real force field calculations.

For example, for l LYS A 3 C TYR B 72 CZ 8.363386
charge1= 0.500; charge2= 0.265; dist = 8.363386

the CNS calculated Eelec energy using OPLS force field is:

0.1 * (0.50.265/8.363386)(1.0 - (8.3633862/8.52))**2 * 332.0636
0.0005348836008092561

Also, for the attached pdb file, the inter-chain energy should be: Evdw = -33.1338 with hydrogens (or -29.729 without hydrogens), Eelec = -8.87.

1A2K_3002.pdb.zip

I think we should move to openMM soon to avoid all these...

[DaniBodor]

For the test, we should compare with real force field calculations.

For example, for l LYS A 3 C TYR B 72 CZ 8.363386 charge1= 0.500; charge2= 0.265; dist = 8.363386

the CNS calculated Eelec energy using OPLS force field is:

0.1 * (0.5_0.265/8.363386)_(1.0 - (8.3633862/8.52))**2 * 332.0636
0.0005348836008092561

Also, for the attached pdb file, the inter-chain energy should be: Evdw = -33.1338 with hydrogens (or -29.729 without hydrogens), Eelec = -8.87.

1A2K_3002.pdb.zip

I think we should move to openMM soon to avoid all these...

Given that we are planning to move to openMM within a reasonable amount of time, do you think it's ok to leave the tests as they are for now, rather than spend time on doing it as you suggest here?

[cbaakman]

Nicely done!