REVISED: Implementation of the density tight-binding (PTB) method

[marcelmbn]

This pull request contains the full implementation of the PTB method as described in: J. Chem. Phys. 158, 124111 (2023) The implementation is based on infrastructure and standard routines provided by tblite. Thus, PTB is only available if WITH_TBLITE is true; in other cases, the program is terminated. A key design principle in this implementation was that the code is located in the xtb source code but (except for printouts and similar things) does not rely on xtb infrastructure but on more modern tblite, mctc-lib, multicharge routines. All empirical parameters are included in the source code, so no parameter file is necessary/possible. This includes also the vDZP basis set, as published with the ωB97X-3c method.

A reference stand-alone implementation is provided in grimme-lab/ptb. The results should be similar, however, the reference implementation partially uses single-precision while the present implementation is (currently) based on double precision matrices throughout. That's why the differences can be at the level of 10⁻⁵ (charges) up to 10⁻¹ (polarizabilities).

A typical calculation is invoked via xtb <struc_file> --ptb. Additionally, the run type --alpha is new and refers to the calculation of the polarizability tensor via PTB (not (yet) possible with GFNn-xTB/FF). The mixed vibrational spectra are accessible with xtb <struc_file> --ptb --hess [--alpha].

The procedure that is undergone after the command above is the following:

• Setup of a new calculator type in src/main/setup.f90 and src/ptb/calculator.F90.
• Add basis set, coordination number and EEQ model.
• Enter singlepoint routine in src/ptb/calculator.F90 and add a new tblite wavefunction. Doing it earlier in the code would be possible but restarting PTB is not possible anyway, which is why I decided against it.
• Entering the central twostepscf procedure (src/ptb/twostepscf.F90), which is designed not as a two-fold iterative cycle but as a one-shot procedure containing two diagonalizations. The reason for this design decision is mainly that the similarity between the first and the second iteration is limited – the employed routines are the same but the parameters entering and some adaptions to the variables are different for many contributions.
• twostepscf works only as a wrapper module for collecting the distinct Hamiltonian contribution. Since PTB works differently compared to GFNn-xTB (many variables are re-used several times), all required overlap matrices and coordination numbers are calculated at the beginning of the actual procedure. Then, the approximated effective core potential (ECP), the shell levels, and the EEQ guess charges are calculated.
  After setting up the EEQ guess and the initial Coulomb matrix, the first iteration starts. Within both iterations, H0 is collected and added onto the (constant) V_ECP. After that, the different density-dependent potentials (2nd and 3rd-order electrostatics, Pauli-XC, DFT+U, and (if present) an external electric field are added. Each iteration ends with diagonalizing the full H matrix and calculating the mixed Mulliken-Löwdin charges.
• twostepscf provides a complete wavefunction including the density and atomic/shell charges. Additionally, the dipole and quadrupole moments, all required integrals, and the bond orders are stored.
• The tblite wavefunction type is then copied onto the xtb wavefunction type to be able to use xtb property evaluation and printout routines down the road.
• In a polarizability calculation (--alpha), numgrad_polarizability is called after the twostepscf. It is a "modified" numerical differentiation of the dipole moment with respect to an infinitesimally small external electric field, in which not the standard twostepscf but a single diagonalization-variant in onestepscf is called. The principle in the construction of the effective Hamiltonian is the same, however.
• Within the singlepoint routine, almost no results are printed except for a "Setup" block similar to what GFN-xTB does.
• For printing out properties, a mixture of xtb and special (new) PTB routines is invoked (see ptb_property and main_ptb_json).

A noteworthy addition within this PR is the functionality of mixing GFNn-xTB/FF frequencies with PTB dipole moment derivatives or polarizability derivatives to obtain higher-quality IR and Raman spectra. Currently, this is limited to the combination of GFN2-xTB+PTB as shown in the PTB publication, but principle-wise this could be handled also generally in the future.
Furthermore, the correct Raman activities and absolute differential Raman scattering cross-sections (in m^2/sr) are implemented into src/hessian.F90 and src/freq/io.f90. In this course, also the vibspectrum file was updated.

The implementation comes with tests for the full Hamiltonian including the response approximation and the polarizability itself. It should work with both meson and CMake.

Moreover, some files were reformatted via fprettify and the following settings:
--- formatting changes were deleted from the present PR ---

[awvwgk]

This PR applies formatting changes to many files, I recommend to split applying the formatter fprettify first in a separate and then add the new features.

[marcelmbn]

This PR applies formatting changes to many files, I recommend to split applying the formatter fprettify first in a separate and then add the new features.

Done. Thank you for the advice, this makes life way easier. 👍

[Albkat]

The description of the CLI options is missing (5d3648f).

[Albkat]

Is there a specific reason why you consider fnv for PTB as a special case?
I tested this locally and it seems to work properly.

it is already done