Fragments with possibly incorrect charge?

[mattiafelice-palermo]

Hello Salim,

I'm trying to parametrize the monomer of the agarose molecule:

Q-Force identifies 8 dihedrals and generates the corresponding fragments. Upon closer inspection, some of the fragments have a +1 charge and multiplicity 2:

* xyz   1   2
  C    -5.484870    -0.530386     1.555744
  C    -4.328913     0.278134     2.112156
  O    -5.062165    -1.817042     1.137077
  H    -5.884522    -0.026790     0.674256
  H    -6.288293    -0.598733     2.297159
  O    -3.859606    -0.438482     3.275576
  C    -4.675721     1.704887     2.508366
  H    -3.529407     0.298694     1.366312
  H    -4.624796    -2.232333     1.890247
[...]

These fragments have been capped with hydrogens as expected where parts of the molecule have been removed. I am not sure whether this is a bug or the intended behaviour, but I thought it was worth reporting it. I expected for neutral molecules the fragment total charge to be 0 and in singlet state, in order to be "closer" to the original molecule. But maybe there's further, more nuanced reasoning behind how Q-Force assigns the charge and spin to the fragments.

In case they can be useful, here are the initial structure (init.xyz) of the agarose monomer I have used, and the input files generated for each fragment: agarose.zip

Please let me know if you need further information and thank you for your insights!

Mattia

[selimsami]

Hi Mattia, That seems weird, I would not expect any charged fragments for this molecule. Could you please send me the Hessian output files as well so I can re-run qforce to generate the fragments? Thanks!

[selimsami]

Also, your fragments don't seem to be much smaller than your molecule, so you might as well close the fragmentation. But still would be interesting to see why the charged fragments appear.

Additionally, I can already predict that a standard relaxed dihedral scan won't go well for your central dihedrals with so many intramolecular hydrogen bonding.

We have almost finished some improved dihedral scans with xtb-torsiondrive (which is already implemented and you can already try this) with a single-point QM calculation at the end matching your Hessian QM method (this part is almost done). I imagine this should improve/solve your future problems... But I'm getting ahead of myself here.

[mattiafelice-palermo]

Thanks for the insight! Here you can find the hessian file: agarose_monomer.zip

Indeed I also did not expect it to work. In fact I noticed the strange charge/multiplicity while writing a script to convert the ORCA input file into torsiondrive input ;)

I tried the xtb implementation in torsiondrive, which I guess for organic molecules would work very well! I wonder if for more exotic molecules (e.g. this I am working on) the geometries will be good enough for recalculating the energies at a higher level of theory.

An idea for the future, since you're already implementing the xtb-torsiondrive, could be to exploit the CREST tool also from the Grimme group to generate conformers to be used as multiple starting point of the wavefront scan. But I digress :)

Thanks again and let me know if you need other information regarding the agarose case.

[selimsami]

@mattiafelice-palermo, I need all the files generated by the QM software, not only .hess.

We have had a discussion about using CREST, we might indeed do that in the future.

You are right that xtb will not be a good starting point in some cases. We also plan to add the option to combine torsiondrive with the other QM software in the future.

[mattiafelice-palermo]

Sure, here are all the files from the QM calculation performed with ORCA (splitted in two zips due to file size upload limit): agarose_1.zip agarose_2.zip

Thanks!

[selimsami]

Hi @mattiafelice-palermo,

Found your issue (well, at least your solution).

In ORCA, the default charge method seems to be for some reason ESP, which is quite unreliable in my opinion! Switching this to Hirshfeld/CM5, which is the default for other software, fixes your issue and all fragments seem to be neutral.

You can do this by:

[qm::software(orca)] charge_method = cm5

I'll also shortly update the code. ORCA module was added by a user, and is mostly not used by the developers (but seems very popular amongst users!), so some of the defaults are a bit different... We might wanna eventually brings these in unison.

[selimsami]

Update pushed. If this fixes your issue, feel free to close it.

[mattiafelice-palermo]

Thank you very much for troubleshooting this! I will check this out for sure :)

Mattia