Feature Contribution: ABF Scan for MM PES in QForce

[mattiafelice-palermo]

Hello Selim,

I hope everything's well!

I've been addressing the issue with the spiky MM PES we previously discussed, I eventually ended up adding the ABF scan feature we discussed about in QForce. My preliminary tests show that it is working as expected However, there are several enhancements I'd like to implement (see mattiafelice-palermo#1 ).

If you believe this feature would be beneficial for QForce, and if it does not conflict with the updates you're currently developing, I'd be happy to create a pull request for further review and integration. In that case, let me know if there's a specific branch you'd prefer for this contribution, or if I should direct it towards master.

Thanks,

Mattia

[selimsami]

HI Mattia,

Thanks for the update. Curious to give it a try!
Is there a short explanation on how it works?

I made a branch called abf_scan. Could you perhaps add it there for some testing?

[mattiafelice-palermo]

Hello,

sure, you can find more information in this tutorial by Luca Muccioli.

Thank you for creating the branch, I will send a pull request. Keep in mind that, compared to the other methods available in QForce, this is slower, but has (should have) the advantage that it works well with the problematic dihedrals we have spoken about.

Mattia

[mattiafelice-palermo]

Hello Selim,

since in the other issue you mentioned releasing a new version of Q-Force, I was wondering if this could fit into the new release.

I have been playing with this in my local machine, and the improvement is definitelty remarkable. For example, take a look at the dihedral profiles for Agarose:

With current Q-Force gromacs implementation:

With ABF:

The two dihedrals that are still not on spot I don't think are due to the ABF procedure per se. I had some issues with Q-Force "not allocating enough multiplicity" for some moieties of other molecules (I will open an issue as soon as I get the time), that was preventing obtaning a good dihedral. Once this was taken into account manually, the dihedral profiles were then well reproduced.

After some tweaking with the colvars module, now the performance of this method is comparable to the standard gromacs one. Also, I managed to simplify the procedure a lot after pushing GROMACS developer to fix a 3 years old bug that was preventing the use of andersen-massive thermostat, and therefore using larger timesteps keeping by constraining the H bonds. This should result in a much cleaner implementation.

If you think this could be a good addition for the new release, then instead of working on implementing this in the current one, I might directly contribute to the new one.

Let me know!

[mattiafelice-palermo]

P.s.: also, I don't know if you have been working into having some sort of "native support" for torsiondrive in Q-Force v2. But since the developers of torsiondrive seem not to be maintaining anymore the repository, I have forked a version where I have implemented support for ORCA, as well as a graceful management of cases where some dihedrals SCF or geometry optimization cycles fail:
https://github.com/mattiafelice-palermo/torsiondrive-on-steroids/tree/master

It would be nice to capitalize on these various efforts and make them available into a new release of Q-Force.

[selimsami]

Hi Mattia,

I am really sorry for never getting around to checking the ABF module.
It's been really hectic, plus I was not (and still am not) sure how ABF fits with all the changes we're planning in 2.0.

Your results look really impressive though, seemingly taking care of one of the main issues I was struggling to solve with Q-Force: unreliability of the dihedral scans.

We are going somewhat on a different direction with Q-Force though, namely we're planning to remove MM relaxations completely and instead fitting the energies and forces on the QM optimized scan structures.

Together with this, also planning other changes that might improve the performance of the dihedrals: 1) adding more bonded coupling terms, 2) fitting the Hessian and dihedrals together, and 3) (hopefully) fitting the equilibrium bonds lengths and angles.

As you can see, with all these changes, I am not sure if ABF will be needed and if it will, how will it fit to the new version.

So my suggestion to you: Wait until the initial 2.0 release (sometime in August), try the same molecule again, and if your issues are not fixed, then adjust your ABF to 2.0 and we'll be more than happy to include this in the final 2.0 release in that case.

I hope that clarifies things a bit. I'll also talk with @MFSJMenger about a native support for torsiondrive, we've also had issues and discussions about this before.

[mattiafelice-palermo]

Hello,

thanks for the update, looks very interesting. Indeed, I have also been studying in these months how these Class I force fields perform when compared to the PES obtained from DFT calculations for large samples of conformations. I also concluded that, without cross terms like bond-bond and bond-angle (the only ones implemented in gromacs, unless one wants to migrate to LAMMPS), a good reproduction of the QM PES is not possible, exception made for very simple molecules.

Following this, I started working as well on minimizing the PES difference by fitting to QM energies, which seems to align with your approach. But I didn't try on dihedrals. The impression I got when experimenting with this is that it can lead to unphysical results, especially if not enough degrees of freedom (in terms of available bonded interactions terms) are available.

So I'm really curious if the inclusion of forces and cross terms in the minimization procedure can actually deliver good results. I will then wait for your new release and see whether this ABF implementation fits with your new approach.

Looking forward to it :)

[selimsami]

Yeah GROMACS is proving very limiting - we are adding OpenMM in 2.0, which has the flexibility of using any random functional form with string equations. Also adding initial MChem support for future...

You can see for some initial small molecules how powerful the cross terms can be to almost exactly reproduce the QM PES:
https://pubs.acs.org/doi/full/10.1021/acs.jpclett.4c00587
However, the challenge is going to be how they will work together with dihedrals and non-bonded interactions.

We are also adding the option to Q-Force to insert additional distorted structures (Wigner, xTB, whatever else). Curious how it's all gonna come together!

[mattiafelice-palermo]

Very interesting work, I will take my time to read through it! Worst case for dihedrals, at least there's a plan B that seems to work reasonably well for quite complex molecules ;)

And thanks for the pointer to OpenMM, indeed the possibility of using string equations is super interesting!

Talk to you soon :)