Okay, here is a Mermaid diagram outlining the general workflow within the src/disp/dftd4.F90 file for calculating the D4 dispersion energy and gradients. This focuses on the main subroutines and their interactions.

graph TD

subgraph Initialization["Model Setup - Usually called externally first"]

InitExt["External Call"] --> NM["newD4Model(...)"]

NM -- "Reads Param Data" --> NM

NM -- "Calculates Reference Polarizabilities" --> NM

NM -- "Calculates Reference C6 Coefficients" --> NM

NM --> DispM["Dispersion Model Data (dispm)"]

end

  

subgraph MainCalculationFlow["Energy/Gradient Calculation"]

Start["Start: d4_gradient_* Routine"] --> InputData{"Input: Mol, dispm, CN, q, Params"}

InputData --> WR["weight_references(...)"]

WR -- "Uses CN, q, zeff, gam, g_a, g_c, wf" --> WR

WR --> F_zeta["zeta/dzeta funcs"]

WR --> F_cngw["cngw/dcngw funcs"]

WR --> WritesWeights["Writes: zetavec, zerovec, derivatives"]

  

InputData --> GAC6_Pair["get_atomic_c6(..., zetavec, ...)"]

DispM --> GAC6_Pair

WritesWeights --> GAC6_Pair

GAC6_Pair --> WritesC6_Pair["Writes: c6_pair, dc6dcn_pair, dc6dq_pair"]

  

InputData --> DG["disp_gradient_* Routine (Pairwise E/Grad)"]

WritesC6_Pair --> DG

DG --> DampingFuncs["Damping Functions"]

DG --> WritesE2["Writes: E_pair, dG_pair, dS_pair, dEdcn_pair, dEdq_pair"]

  

InputData -- "If s9 != 0" --> GAC6_ATM["get_atomic_c6(..., zerovec, ...)"]

DispM --> GAC6_ATM

WritesWeights --> GAC6_ATM

GAC6_ATM --> WritesC6_ATM["Writes: c6_atm, dc6dcn_atm"]

  

InputData -- "If s9 != 0" --> AG["atm_gradient_* Routine (3-Body E/Grad)"]

WritesC6_ATM --> AG

AG --> DAT["deriv_atm_triple func"]

AG --> TS["triple_scale func"]

AG --> WritesE3["Writes: E_atm, dG_atm, dS_atm, dEdcn_atm"]

  

InputData --> Combine{"Combine Results"}

WritesE2 --> Combine

WritesE3 -- "If s9 != 0" --> Combine

Combine -- "Uses dcndr, dqdr (Chain Rule)" --> Combine

Combine --> FinalOutput["Output: Total E, Gradient, Sigma"]

end

  

%% Style adjustments

classDef mySubgraph fill:#f9f,stroke:#333,stroke-width:2px;

class InitExt,Start,FinalOutput fill:#ccf,stroke:#333,stroke-width:2px;

Explanation of the Flow:

  1. Initialization (newD4Model): This sets up the fundamental dispm data structure containing reference polarizabilities (dispm%alpha) and reference C6 coefficients (dispm%c6) based on pre-defined parameters and the chosen charge model. This is typically done once before the main SCF or geometry optimization loop.
  2. Main Entry (d4_gradient_*): These routines (like d4_gradient_neigh or d4_gradient_latp) are the main entry points during a calculation. They take molecular geometry, pre-calculated dispm data, coordination numbers (cn), charges (q), and D4 parameters (par, g_a, g_c, wf) as input.
  3. Weight Calculation (weight_references): Calculates the environment-dependent weights (zetavec) for each reference system of each atom, considering the atom’s current coordination number and charge. It also calculates weights for the zero-charge case (zerovec) and derivatives needed later. Helper functions zeta and cngw are used here.
  4. Atomic C6 Calculation (Pairwise - get_atomic_c6 with zetavec): Uses the calculated weights (zetavec) and the reference C6 values (dispm%c6) to compute the final, environment-dependent C6 coefficient (c6_pair) for each pair of atoms. It also computes the derivatives of these C6 values w.r.t CN (dc6dcn_pair) and charge (dc6dq_pair).
  5. Pairwise Energy/Gradient (disp_gradient_*): Calculates the two-body dispersion energy (E_pair) and its geometric gradient (dG_pair, dS_pair) using the final atomic C6 values (c6_pair), damping functions, and interatomic distances. It also calculates the derivatives of the energy w.r.t. CN (dEdcn_pair) and charge (dEdq_pair).
  6. Atomic C6 Calculation (3-Body - get_atomic_c6 with zerovec): If the three-body term is active (par%s9 /= 0.0_wp), this step calculates the C6 coefficients needed for the ATM term, typically using the zero-charge weights (zerovec).
  7. Three-Body Energy/Gradient (atm_gradient_*): If the three-body term is active, this calculates the Axilrod-Teller-Muto energy (E_atm) and its geometric gradient (dG_atm, dS_atm) using the appropriate C6 values (c6_atm). It also calculates the derivative w.r.t. CN (dEdcn_atm). Helper functions deriv_atm_triple and triple_scale are used.
  8. Combine Results: The main d4_gradient_* routine sums the pairwise and (if applicable) three-body energy contributions. It combines the geometric gradients (dG, dS) with the CN/charge derivatives (dEdcn, dEdq) using the chain rule (multiplying by dcndr, dqdr, etc.) to get the final total gradient and stress tensor.
  9. Output: Returns the total D4 energy, gradient, and stress tensor.

This diagram provides a high-level overview of how data flows through the key subroutines in dftd4.F90 to arrive at the final dispersion correction.

D4 dispersion

Okay, let’s break down how the force (which is the negative of the energy gradient, ( $\mathbf{F} = -\nabla E )$) is calculated in src/disp/dftd4.F90.

The total D4 dispersion energy depends on the atomic coordinates ((\mathbf{R})) in multiple ways:

  1. Directly: Through the distances ((r_{ij})) in the (r^{-n}) terms and the damping functions.
  2. Indirectly: Through the coordination numbers ((CN)) and atomic charges ((q)), which themselves depend on the geometry ((CN(\mathbf{R}), q(\mathbf{R}))).

Therefore, the gradient calculation uses the chain rule:

$$
\nabla_{\mathbf{R}} E_{D4} = \underbrace{\frac{\partial E_{D4}}{\partial \mathbf{R}}}{\text{Direct Term}} + \sum_k \underbrace{\frac{\partial E{D4}}{\partial CN_k}}{\text{dEdcn}} \underbrace{\nabla{\mathbf{R}} CN_k}{\text{dcndr}} + \sum_k \underbrace{\frac{\partial E{D4}}{\partial q_k}}{\text{dEdq}} \underbrace{\nabla{\mathbf{R}} q_k}_{\text{dqdr}}
$$

Here’s how the code calculates these terms:

  1. Direct Geometric Gradient ($\frac{\partial E_{D4}}{\partial \mathbf{R}}$):

    • This is calculated within the lower-level gradient routines: disp_gradient_neigh/disp_gradient_latp for pairwise terms and atm_gradient_neigh/atm_gradient_latp (via deriv_atm_triple) for three-body terms.
    • They differentiate the damped (r^{-n}) expressions with respect to the coordinates.
    • For example, in disp_gradient_neigh (line 1422), ddisp holds the derivative of the damped terms w.r.t (r^2), and dG = -c6(iat, jat)*ddisp*rij calculates the contribution to the gradient vector between atoms iat and jat.
    • These direct contributions (dG) are accumulated directly into the gradient array (e.g., lines 1429-1430).

  2. Energy Derivatives w.r.t CN and Charge (($\frac{\partial E_{D4}}{\partial CN_k}) and (\frac{\partial E_{D4}}{\partial q_k}$)):

    • These terms represent how sensitive the dispersion energy is to changes in coordination number or charge at a fixed geometry.
    • They are calculated within the disp_gradient_neigh/disp_gradient_latp routines using the previously computed dc6dcn and dc6dq.
    • Example from disp_gradient_neigh (lines 1420-1421, 1425-1426):
      1
      2
      3
      4
      5
      ! dEdcn(k) accumulates dE/dCN_k
      dEdcn(iat) = dEdcn(iat) - dc6dcn(iat, jat) * disp
      ! dEdq(k) accumulates dE/dq_k
      dEdq(iat) = dEdq(iat) - dc6dq(iat, jat) * disp
      ! ... symmetric terms for jat ...
    • Similarly, the three-body gradient routines (atm_gradient_*) accumulate contributions to dEdcn (lines 1614-1616, via dCN from deriv_atm_triple).

  3. Chain Rule Application (Indirect Terms):

    • The geometric derivatives of CN (( \nabla_{\mathbf{R}} CN_k )) and charge (( \nabla_{\mathbf{R}} q_k )) are provided as input arrays (dcndr and dqdr, respectively).
    • The final application of the chain rule (summing the products) happens in the main wrapper routines (d4_full_gradient_neigh, d4_gradient_neigh, etc.).
    • Matrix-vector multiplications (mctc_gemv) perform this summation efficiently. Example from d4_full_gradient_neigh (lines 1235-1240):
      1
      2
      3
      4
      5
      6
      7
      ! gradient = gradient + dcndr^T * dEdcn
      call mctc_gemv(dcndr, dEdcn, gradient, beta=1.0_wp)
      ! gradient = gradient + dqdr^T * dEdq  (if dqdr is present)
      if (present(dqdr)) then
         call mctc_gemv(dqdr, dEdq, gradient, beta=1.0_wp)
      end if
      ! ... similar terms for stress tensor sigma using dcndL, dqdL ...
    • The beta=1.0_wp argument means the result is added to the existing gradient array (which already contains the direct geometric derivative from step 1).

In summary:

The final gradient array accumulated in the main d4_*_gradient_* routines contains the total gradient of the D4 energy (($\nabla_{\mathbf{R}} E_{D4})$). It’s computed by summing:

  • The direct derivatives from differentiating distance-dependent terms.
  • The indirect contributions calculated via the chain rule, using the energy sensitivities (dEdcn, dEdq) multiplied by the geometric derivatives of CN and charge (dcndr, dqdr).

The force on each atom is then simply the negative of the corresponding elements in this final gradient array.

subgraph one
a1-->a2
end

$$
\textbf{math}test
$$