How to use dispersion correction with TPSSh in Gaussian 16?
The answer is quite simple and already given in the manual:


  You can use this empirical dispersion method with other functionals via the IOps(3/174-178) (S6 should be 1.0).


Parameters and the standalone program can be found can be found on the DFTD3 homepage at Mulliken Center for Theoretical Chemistry. Here are links via The Internet archive to the tables for zero-damping and BJ-damping.

The Overlay3 manual can be found on gaussian.com, or The Internet Archive:


  IOp(3/174)
  
  S6 scale factor in Grimme’s D2/D3/D3BJ dispersion.

0          Default (see subroutine R6DS6).
-1         Set S6 to 0. 
NNNNNNNN   A value of NNNNNNNN/1000000. 

  
  
  
  IOp(3/175)
  
  S8 scale factor in Grimme’s D2/D3/D3BJ dispersion.

0          Default (see subroutine R6DS8).
-1         Set S8 to 0. 
NNNNNNNN   A value of NNNNNNNN/1000000. 

  
  
  
  IOp(3/176)
  
  SR6 scale factor in Grimme’s D2/D3/D3BJ dispersion.

0          Default (see subroutine R6DSR6).
-1         Set SR6 to 0. 
NNNNNNNN   A value of NNNNNNNN/1000000. 

  
  
  
  IOp(3/177)
  
  A1 parameter in Becke-Johnson damping for D3BJ and XDM.

0          Default (see subroutine R6DABJ/XDMABJ).
-1         Set A1 to 0. 
NNNNNNNN   A value of NNNNNN/1000000. 

  
  
  
  IOp(3/178)
  
  A2 parameter in Becke-Johnson damping for D3BJ and XDM.

0          Default (see subroutine R6DABJ/XDMABJ).
-1         Set A2 to 0. 
NNNNNNNN   A value of NNNNNN/1000000 Ang.



For TPSSh the route section should therefore look similar to the following:

#P TPSSh/def2TZVP           ! verbose, method/ basis set
EmpiricalDispersion(GD3BJ)  ! Switch on D3 with Becke-Johnson damping
IOp(3/174=1000000)          ! Set s6 to 1.000
IOp(3/175=2238200)          ! Set s8 to 2.2382
IOp(3/177=452900)           ! Set a1 to 0.4529
IOp(3/178=4655000)          ! Set a2 to 4.6550
opt(MaxCycle=100,VeryTight) ! Run optimisation
int(ultrafinegrid)          ! use finer grid for dfa