Using a ready made topology file instead of the Hessain calculation

[xiki-tempula]

I'm interested in supplying a Gromacs topology file for the small molecule.
I wonder if the Qforce could load the topology, identify the soft torsions and then run QM calculations to refit the dihedral.
Then read the refitted QM data and spit out a new topology file with the refitted dihedrals.

[selimsami]

Definitely possible but would require decent amount of coding.

Writing an ITP reader is a complex task, would be ideal to outsource this to existing packages. If perhaps MDAnalysis could be used to easily read topology files, then that would simplify the problem.

Then one would write the read FF parameters to the Q-Force terms object and skip the Hessian fitting stage and connect directly to the dihedral fitting.

If you are interested in tackling this, I could give more pointers.

[xiki-tempula]

So my thought is that parmed could be used, depending on the detailed implementation, this could get the reading and writing topology done in one go. I contributed code to MDAnalysis but I felt the current MDAnalysis is too bulky.

Definitely possible but would require decent amount of coding.

I will need to think about this when I got the TorsionDrive thing done and then see if I have time for it.

[selimsami]

Parmed might indeed be a better option, MDAnalysis is indeed too bulky.

[xiki-tempula]

@selimsami If we go this route, do we still need the NBO?

[selimsami]

If you want Q-Force to determine which dihedrals are flexible, you need at least WBO analysis with, say xTB. If Q-Force is provided which dihedrals are flexible (also possible since anyway it's reading a FF file), then you could even do away with the WBO analysis. NBO I think is unnecessary, I will also phase it out eventually. Bond orders seem to be sufficient.

[xiki-tempula]

@selimsami My boss wants me to do this, so I will start working on this now.

So I would imagine that I could provide the topology information
Bond, angle, dihedrals, charge, LJ, mixing rules and fudge LJ/Q.

and the information of NBO. I wonder would you mind point me toward where should I supply this information? Thanks.

[xiki-tempula]

@selimsami Sorry, I wonder how do you like the way that I write do string for the orca interface? I find it a bit difficult to understand the code without docstring. I wonder if you mind that I write the docstring in the format that I used in the orca parser? Thanks.

[selimsami]

I think I would make a new function similar to run_qforce / run_hessian_fitting_for_external that is in main.py.

We have some plans to add checkpoints to Q-Force, so it knows which stage it is and can read data. In that case we might restructure things a bit. But all this would happen in main.py anyway, so no big changes for this PR.

Docstrings in orca interface are perfectly fine! Defintely feel free to use them. That is my bad, I should also add to existing code as well.

[xiki-tempula]

@selimsami Thanks. I wonder if there is an example of run_hessian_fitting_for_external? So I could use IDE to track what is the input and what is the expected output.

[xiki-tempula]

Ok, it is the same as run_qforce. I will have a look, Thanks.

[selimsami]

Maybe this helps.

qforce_interface.zip

If i'm not mistaken, you need to add a new line for "ext_h_cap".

[xiki-tempula]

@selimsami Thanks for the help and I'm reading the code now.

I might be naive but I think we don't really need the LP? The only place that I think it is used is in https://github.com/selimsami/qforce/blob/master/qforce/molecule/topology.py#L44 as an identifier.
I might miss something as I haven't read the code after this step yet.

If this is the case, we could swap all the LP with the pre atom total bond order? At the end of day, if LP is the same, the total bond order is the same.

[selimsami]

If this is the case, we could swap all the LP with the pre atom total bond order? At the end of day, if LP is the same, the total bond order is the same.

Yes, that's exactly my plan for depreciating the LPs :)

[xiki-tempula]

@selimsami Before you depreciate the LPs, do you think it would be fine if I just pass the integer total bond order to the lone_e? So I could do some testing. Thanks.