Coverage for cosmolayer / cosmosac / interaction_matrices.py: 100%

1""" 

2.. module:: cosmolayer.cosmosac.interaction_matrices 

3 :synopsis: Create interaction matrices for COSMO-SAC calculations. 

4 

5.. functionauthor:: Charlles Abreu <craabreu@gmail.com> 

6""" 

7 

8from collections import defaultdict 

9 

10import numpy as np 

11from numpy.typing import NDArray 

12 

13from .segment_groups import OH, OT, SEGMENT_GROUPS 

14 

15 

16def create_cosmo_sac_2002_matrix( # noqa: PLR0913 

17 temperature: float, # K 

18 *, 

19 min_sigma: float = -0.025, 

20 max_sigma: float = 0.025, 

21 num_points: int = 51, 

22 sigma_hb: float = 0.0084, # e/Ų 

23 alpha_prime: float = 16466.72, # (kcal/mol)/(e/Ų)² 

24 c_hb: float = 85580.0, # (kcal/mol)/(e/Ų)² 

25 gas_constant: float = 0.001987, # kcal/(mol·K) 

26) -> NDArray[np.float64]: 

27 r"""Create an interaction matrix for the COSMO-SAC 2002 model :cite:`Bell2020`. 

28 

29 Computes the pairwise interaction energies between surface segments with given 

30 screening charge densities, ΔW(σ,σ'), divided by the product RT, where R is the 

31 universal gas constant and T is the temperature at which the interaction matrix 

32 is computed. 

33 

34 Parameters 

35 ---------- 

36 temperature : float 

37 The temperature in Kelvin at which the interaction matrix is computed. 

38 

39 Keyword Arguments 

40 ----------------- 

41 min_sigma : float, optional 

42 Minimum screening charge density in e/Ų. Default is -0.025. 

43 max_sigma : float, optional 

44 Maximum screening charge density in e/Ų. Default is 0.025. 

45 num_points : int, optional 

46 Number of discrete points in the sigma grid. Default is 51. 

47 sigma_hb : float, optional 

48 Hydrogen bonding cutoff parameter in e/Ų. Defines the range for 

49 hydrogen bonding interactions. Default is 0.0084 :cite:`Bell2020`. 

50 alpha_prime : float, optional 

51 Misfit energy constant in (kcal/mol)/(e/Ų)². Controls the strength 

52 of electrostatic misfit interactions. Default is 16466.72 :cite:`Bell2020`. 

53 c_hb : float, optional 

54 Hydrogen bonding energy constant in (kcal/mol)/(e/Ų)². Controls the 

55 strength of hydrogen bonding interactions. Default is 85580.0 :cite:`Bell2020`. 

56 gas_constant : float, optional 

57 Universal gas constant in kcal/(mol·K). Default is 0.001987 :cite:`Bell2020`. 

58 

59 Returns 

60 ------- 

61 np.ndarray 

62 Dimensionless interaction energy matrix ΔW(σ,σ') / (RT). 

63 Shape: (num_points, num_points). 

64 """ 

65 

66 grid = np.linspace(min_sigma, max_sigma, num_points) 

67 squared_sum_block = np.add.outer(grid, grid) ** 2 

68 delta = (grid - sigma_hb).clip(min=0) + (grid + sigma_hb).clip(max=0) 

69 hb_block = np.outer(delta, delta).clip(max=0) 

70 energy_matrix = (alpha_prime / 2) * squared_sum_block + c_hb * hb_block 

71 result: NDArray[np.float64] = energy_matrix / (gas_constant * temperature) 

72 return result 

73 

74 

75def create_cosmo_sac_2010_matrices( # noqa: PLR0913 

76 temperature: float, # K 

77 *, 

78 min_sigma: float = -0.025, 

79 max_sigma: float = 0.025, 

80 num_points: int = 51, 

81 a_es: float = 6525.69, # (kcal/mol)/(e/Ų)² 

82 b_es: float = 1.4859e8, # (kcal/mol)K²/(e/Ų)² 

83 c_oh_oh: float = 4013.78, # kcal·Å^4·mol⁻¹·e⁻² 

84 c_ot_ot: float = 932.31, # kcal·Å^4·mol⁻¹·e⁻² 

85 c_oh_ot: float = 3016.43, # kcal·Å^4·mol⁻¹·e⁻² 

86 gas_constant: float = 0.0019872043011606513, # kcal/(mol·K) 

87) -> tuple[NDArray[np.float64], NDArray[np.float64]]: 

88 r"""Create interaction matrices for the COSMO-SAC 2010 model :cite:`Bell2020`. 

89 

90 Computes the electrostatic and hydrogen bonding parts of pairwise interaction 

91 energies between surface segments with given screening charge densities, ΔW(σ,σ'), 

92 divided by the product RT, where R is the universal gas constant and T is the 

93 temperature at which the interaction matrix is computed. 

94 

95 Parameters 

96 ---------- 

97 temperature : float 

98 The temperature in Kelvin at which the interaction matrix is computed. 

99 

100 Keyword Arguments 

101 ----------------- 

102 min_sigma : float, optional 

103 Minimum screening charge density in e/Ų. Default is -0.025. 

104 max_sigma : float, optional 

105 Maximum screening charge density in e/Ų. Default is 0.025. 

106 num_points : int, optional 

107 Number of discrete points in the sigma grid. Default is 51. 

108 a_es : float, optional 

109 Misfit energy constant in (kcal/mol)/(e/Ų)². Controls the strength 

110 of electrostatic misfit interactions. Default is 6525.69 :cite:`Bell2020`. 

111 b_es : float, optional 

112 Misfit energy constant in (kcal/mol)/(e/Ų)². Controls the strength 

113 of electrostatic misfit interactions. Default is 1.4859e8 :cite:`Bell2020`. 

114 c_oh_oh : float, optional 

115 Hydrogen bonding energy constant in (kcal/mol)/(e/Ų)². Controls the 

116 strength of hydrogen bonding interactions. Default is 4013.78 :cite:`Bell2020`. 

117 c_ot_ot : float, optional 

118 Hydrogen bonding energy constant in (kcal/mol)/(e/Ų)². Controls the 

119 strength of hydrogen bonding interactions. Default is 932.31 :cite:`Bell2020`. 

120 c_oh_ot : float, optional 

121 Hydrogen bonding energy constant in (kcal/mol)/(e/Ų)². Controls the 

122 strength of hydrogen bonding interactions. Default is 3016.43 :cite:`Bell2020`. 

123 gas_constant : float, optional 

124 Universal gas constant in kcal/(mol·K). Default is 0.0019872043011606513 

125 :cite:`Bell2020` 

126 

127 Returns 

128 ------- 

129 tuple[np.ndarray, np.ndarray] 

130 Dimensionless interaction energy matrices ΔW(σ,σ') / (RT). 

131 Shape: (num_points, num_points). 

132 """ 

133 RT = gas_constant * temperature 

134 c_hb: defaultdict[str, defaultdict[str, float]] = defaultdict( 

135 lambda: defaultdict(float) 

136 ) 

137 c_hb[OH][OH] = c_oh_oh 

138 c_hb[OT][OT] = c_ot_ot 

139 c_hb[OH][OT] = c_hb[OT][OH] = c_oh_ot 

140 grid = np.linspace(min_sigma, max_sigma, num_points) 

141 

142 es_block = np.add.outer(grid, grid) ** 2 / RT 

143 hb_block = (np.outer(grid, grid) < 0) * np.subtract.outer(grid, grid) ** 2 / RT 

144 

145 es_matrix = np.block([[es_block] * 3] * 3) 

146 hb_matrix = np.block( 

147 [[c_hb[s][t] * hb_block for t in SEGMENT_GROUPS] for s in SEGMENT_GROUPS] 

148 ) 

149 

150 result_a: NDArray[np.float64] = a_es * es_matrix - hb_matrix 

151 result_b: NDArray[np.float64] = (b_es / temperature**2) * es_matrix 

152 return result_a, result_b