1C> \ingroup nwxc 2C> @{ 3C> 4C> \file nwxcm_c_pw91lda.F 5C> The nwxcm_c_pw91lda functional 6C> 7C> @} 8C> 9C> \ingroup nwxc_priv 10C> @{ 11C> 12C> \brief Evaluate the nwxcm_c_pw91lda functional [1] 13C> 14C> \f{eqnarray*}{ 15C> {\it t_1} &=& \rho_\beta+\rho_\alpha\\\\ 16C> {\it t_2} &=& {\it t_1}^{0.3333333333333333}\\\\ 17C> {\it t_3} &=& {{1}\over{{\it t_2}}}\\\\ 18C> {\it t_4} &=& 0.1325688999052018\,{\it t_3}+1.0\\\\ 19C> {\it t_5} &=& \sqrt{{\it t_2}}\\\\ 20C> {\it t_6} &=& {{1}\over{{\it t_5}}}\\\\ 21C> {\it t_7} &=& \log \left({{1.269642451250142\,{\it t_5}} 22C> \over{0.7876233178997433\,{\it t_6}\,\left(0.7876233178997433 23C> \,{\it t_6}\,\left(0.02414199311453321\,{\it t_6} 24C> +0.10186556948\right)+0.22308199064\right)+0.47231125998}} 25C> +1.0\right)\\\\ 26C> {\it t_8} &=& \rho_\alpha-\rho_\beta\\\\ 27C> {\it t_9} &=& {{1}\over{{\it t_1}}}\\\\ 28C> {\it t_{10}} &=& 0.06901399211255826\,{\it t_3}+1.0\\\\ 29C> {\it t_{11}} &=& \log \left({{1.269642451250142\,{ 30C> \it t_5}}\over{0.7876233178997433\,{\it t_6}\, 31C> \left(0.7876233178997433\,{\it t_6}\,\left(0.01321299881039884 32C> \,{\it t_6}+0.029729725188\right)+0.12236585478\right) 33C> +0.3497952466}}+1.0\right)\\\\ 34C> {\it t_{12}} &=& \rho_s^{0.3333333333333333}\\\\ 35C> {\it t_{13}} &=& \sqrt{{\it t_{12}}}\\\\ 36C> {\it t_{14}} &=& {{1}\over{{\it t_{13}}}}\\\\ 37C> {\it t_{15}} &=& {{1}\over{{\it t_{12}}}}\\\\ 38C> {\it t_{16}} &=& \log \left({{1.269642451250142\,{ 39C> \it t_{13}}}\over{0.7876233178997433\,\left(0.7876233178997433 40C> \,\left(0.02414199311453321\,{\it t_{14}}+0.10186556948\right) 41C> \,{\it t_{14}}+0.22308199064\right)\,{\it t_{14}} 42C> +0.47231125998}}+1.0\right)\\\\ 43C> {\it t_{17}} &=& 0.1325688999052018\,{\it t_{15}}+1.0\\\\ 44C> {\it t_{18}} &=& \log \left({{1.269642451250142\,{ 45C> \it t_{13}}}\over{0.7876233178997433\,\left(0.7876233178997433 46C> \,\left(0.01321299881039884\,{\it t_{14}} 47C> +0.029729725188\right)\,{\it t_{14}}+0.12236585478\right)\,{ 48C> \it t_{14}}+0.3497952466}}+1.0\right)\\\\ 49C> {\it t_{19}} &=& 0.06901399211255826\,{\it t_{15}}+1.0\\\\ 50C> f &=& 1.0\,{\it t_1}\,\left(0.5848223622634648\, 51C> \left(1.923661050931536\,\left({\it t_8}\,{\it t_9} 52C> +1.0\right)^{{{4}\over{3}}}+1.923661050931536\,\left(1.0-{ 53C> \it t_8}\,{\it t_9}\right)^{{{4}\over{3}}} 54C> -3.847322101863072\right)\,\left({{{\it t_8}^4\, 55C> \left(1.709920934161365\,\left(0.0621814\,{\it t_4}\,{\it t_7} 56C> -0.0310907\,\left(0.1274696188700087\,{\it t_3}+1.0\right) 57C> \,\log \left({{1.269642451250142\,{\it t_5}} 58C> \over{0.7876233178997433\,{\it t_6}\,\left(0.7876233178997433 59C> \,{\it t_6}\,\left(0.01530901310039024\,{\it t_6} 60C> +0.10465751434\right)+0.19269083139\right)+0.43896648423}} 61C> +1.0\right)\right)-0.0337738\,{\it t_{10}}\,{ 62C> \it t_{11}}\right)}\over{{\it t_1}^4}}+0.0337738\,{\it t_{10}} 63C> \,{\it t_{11}}\right)-0.0621814\,{\it t_4}\,{\it t_7}\right)\\\\ 64C> g &=& 0\\\\ 65C> G &=& 1.0\,\left(0.5848223622634643\,\left(0.0337738\,{ 66C> \it t_{18}}\,{\it t_{19}}+1.0\,\left(1.709920934161365\, 67C> \left(0.0621814\,{\it t_{16}}\,{\it t_{17}}-0.0310907 68C> \,\log \left({{1.269642451250142\,{\it t_{13}}} 69C> \over{0.7876233178997433\,\left(0.7876233178997433\, 70C> \left(0.01530901310039024\,{\it t_{14}}+0.10465751434\right) 71C> \,{\it t_{14}}+0.19269083139\right)\,{\it t_{14}} 72C> +0.43896648423}}+1.0\right)\,\left(0.1274696188700087\,{ 73C> \it t_{15}}+1.0\right)\right)-0.0337738\,{\it t_{18}}\,{ 74C> \it t_{19}}\right)\right)-0.0621814\,{\it t_{16}}\,{ 75C> \it t_{17}}\right)\,\rho_s\\\\ 76C> \f} 77C> 78C> Code generated with Maxima 5.34.0 [2] 79C> driven by autoxc [3]. 80C> 81C> ### References ### 82C> 83C> [1] JP Perdew, Y Wang, Phys.Rev.B 45, 13244 (1992) , DOI: 84C> <a href="https://doi.org/10.1103/PhysRevB.45.13244 "> 85C> 10.1103/PhysRevB.45.13244 </a> 86C> 87C> [2] Maxima, a computer algebra system, 88C> <a href="http://maxima.sourceforge.net/"> 89C> http://maxima.sourceforge.net/</a> 90C> 91C> [3] autoxc, revision 27097 2015-05-08 92C> 93 subroutine nwxcm_c_pw91lda(param,tol_rho,ipol,nq,wght, 94 +rho,fnc,Amat) 95c $Id: $ 96#ifdef NWXC_QUAD_PREC 97 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 98 integer, parameter :: rk=selected_real_kind(30) 99#else 100 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 101 integer, parameter :: rk=selected_real_kind(15) 102#endif 103#include "nwxc_param.fh" 104 double precision param(*) !< [Input] Parameters of functional 105 double precision tol_rho !< [Input] The lower limit on the density 106 integer ipol !< [Input] The number of spin channels 107 integer nq !< [Input] The number of points 108 double precision wght !< [Input] The weight of the functional 109 double precision rho(nq,*) !< [Input] The density 110 double precision fnc(nq) !< [Output] The value of the functional 111c 112c Sampling Matrices for the XC Kernel 113c 114 double precision Amat(nq,*) !< [Output] The derivative wrt rho 115 integer iq 116 double precision tmp 117 double precision rhoa,rhob 118 double precision gammaaa,gammaab,gammabb 119 double precision taua,taub 120 double precision nwxcm_heaviside 121 external nwxcm_heaviside 122CDIR$ NOVECTOR 123 do iq = 1, nq 124 if (ipol.eq.1) then 125 rhoa = 0.5d0*rho(iq,R_T) 126 if (rhoa.gt.tol_rho) then 127 t1 = rhoa**3.333333333333333d-1 128 t2 = t1**5.0d-1 129 t3 = 1/t2 130 t4 = 1/t1 131 t5 = 2.1508070719090538d-2*t3+1.0186556948d-1 132 t6 = 7.016926042943222d-1*t3*t5+2.2308199064d-1 133 t7 = 7.016926042943222d-1*t3*t6+4.7231125998d-1 134 t8 = 1/t7 135 t9 = 1.425125466450768d+0*t2*t8+1.0d+0 136 t10 = log(t9) 137 t11 = 1.0522000558389213d-1*t4+1.0d+0 138 t12 = 1/t2**3 139 t13 = 1/rhoa**6.666666666666667d-1 140 fnc(iq) = fnc(iq)-1.243628d-1*rhoa*log(1.4251254664507676d+0 141 1 *t2/(7.016926042943223d-1*t3*(7.016926042943223d-1*(2.150 142 2 807071909054d-2*t3+1.0186556948d-1)*t3+2.2308199064d-1)+4 143 3 .7231125998d-1)+1.0d+0)*(1.0522000558389215d-1*t4+1.0d+0) 144 4 *wght 145 Amat(iq,D1_RA) = 2.0d+0*rhoa*(1.090454542535705d-3*t10/rhoa* 146 1 *1.3333333333333333d+0-6.21814d-2*t11*(1.1876045553756398 147 2 d-1*t13*t3*t8-1.425125466450768d+0*t2*(7.016926042943222d 148 3 -1*t3*(-5.847438369119352d-2*t12*t13*t5-1.257671179685424 149 4 2d-3/rhoa**1.3333333333333336d+0)-5.847438369119352d-2*t1 150 5 2*t13*t6)/t7**2)/t9)*wght-6.21814d-2*t10*t11*wght+Amat(iq 151 6 ,D1_RA) 152 endif ! rhoa.gt.tol_rho 153 else ! ipol.eq.1 154 rhoa = rho(iq,R_A) 155 rhob = rho(iq,R_B) 156 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 157 t1 = rhob+rhoa 158 t2 = t1**3.333333333333333d-1 159 t3 = 1/t2 160 t4 = 1.325688999052018d-1*t3+1.0d+0 161 t5 = t2**5.0d-1 162 t6 = 1/t5 163 t7 = 2.4141993114533214d-2*t6+1.0186556948d-1 164 t8 = 7.876233178997433d-1*t6*t7+2.2308199064d-1 165 t9 = 7.876233178997433d-1*t6*t8+4.7231125998d-1 166 t10 = 1/t9 167 t11 = 1.269642451250142d+0*t10*t5+1.0d+0 168 t12 = log(t11) 169 t13 = rhoa-rhob 170 t14 = 1/t1 171 t15 = 1.0d+0-t13*t14 172 t16 = t13*t14+1.0d+0 173 t17 = 1.923661050931536d+0*t16**1.3333333333333333d+0+1.9236 174 1 61050931536d+0*t15**1.3333333333333333d+0-3.8473221018630 175 2 72d+0 176 t18 = 6.901399211255826d-2*t3+1.0d+0 177 t19 = 1.3212998810398843d-2*t6+2.9729725188d-2 178 t20 = 7.876233178997433d-1*t19*t6+1.2236585478d-1 179 t21 = 7.876233178997433d-1*t20*t6+3.497952466d-1 180 t22 = 1/t21 181 t23 = 1.269642451250142d+0*t22*t5+1.0d+0 182 t24 = log(t23) 183 t25 = t13**4 184 t26 = 1/t1**4 185 t27 = 1.2746961887000874d-1*t3+1.0d+0 186 t28 = 1.530901310039024d-2*t6+1.0465751434d-1 187 t29 = 7.876233178997433d-1*t28*t6+1.9269083139d-1 188 t30 = 7.876233178997433d-1*t29*t6+4.3896648423d-1 189 t31 = 1/t30 190 t32 = 1.269642451250142d+0*t31*t5+1.0d+0 191 t33 = log(t32) 192 t34 = 1.709920934161365d+0*(6.21814d-2*t12*t4-3.10907d-2*t27 193 1 *t33)-3.37738d-2*t18*t24 194 t35 = t25*t26*t34+3.37738d-2*t18*t24 195 t36 = 5.848223622634648d-1*t17*t35-6.21814d-2*t12*t4 196 t37 = 1/t1**1.3333333333333336d+0 197 t38 = 1/t1**6.666666666666667d-1 198 t39 = 1/t5**3 199 t40 = 2.11607075208357d-1*t10*t38*t6-1.269642451250142d+0*t5 200 1 *(7.876233178997433d-1*t6*(-1.3127055298329054d-1*t38*t39 201 2 *t7-3.169132786263567d-3*t37)-1.3127055298329054d-1*t38*t 202 3 39*t8)/t9**2 203 t41 = 1/t11 204 t42 = -6.21814d-2*t4*t40*t41 205 t43 = 1/t1**1.3333333333333333d+0 206 t44 = 2.747773264188438d-3*t12*t43 207 t45 = 2.11607075208357d-1*t22*t38*t6-1.269642451250142d+0*t5 208 1 *(7.876233178997433d-1*(-1.3127055298329054d-1*t19*t38*t3 209 2 9-1.7344776604086162d-3*t37)*t6-1.3127055298329054d-1*t20 210 3 *t38*t39)/t21**2 211 t46 = 1/t23 212 t47 = 3.37738d-2*t18*t45*t46 213 t48 = -7.769549222703733d-4*t24*t43 214 t49 = t25*t26*(1.709920934161365d+0*(-3.10907d-2*t27*(2.1160 215 1 7075208357d-1*t31*t38*t6-1.269642451250142d+0*t5*(7.87623 216 2 3178997433d-1*(-1.3127055298329054d-1*t28*t38*t39-2.00962 217 3 26153166658d-3*t37)*t6-1.3127055298329054d-1*t29*t38*t39) 218 4 /t30**2)/t32+1.3210398931339265d-3*t33*t43-2.747773264188 219 5 438d-3*t12*t43+6.21814d-2*t4*t40*t41)-3.37738d-2*t18*t45* 220 6 t46+7.769549222703733d-4*t24*t43) 221 t50 = -4*t25*t34/t1**5 222 t51 = t13**3 223 t52 = 1/t1**2 224 t53 = t13*t52 225 t54 = -t14 226 t55 = t15**3.333333333333333d-1 227 t56 = -t13*t52 228 t57 = t16**3.333333333333333d-1 229 t58 = 1.0d+0*t36*wght 230 fnc(iq) = 1.0d+0*t1*t36*wght+fnc(iq) 231 Amat(iq,D1_RA) = 1.0d+0*t1*(5.848223622634648d-1*t35*(2.5648 232 1 81401242048d+0*(t56+t14)*t57+2.564881401242048d+0*(t54+t5 233 2 3)*t55)+5.848223622634648d-1*t17*(4*t26*t34*t51+t50+t49+t 234 3 48+t47)+t44+t42)*wght+t58+Amat(iq,D1_RA) 235 Amat(iq,D1_RB) = 1.0d+0*t1*(5.848223622634648d-1*t35*(2.5648 236 1 81401242048d+0*(t56+t54)*t57+2.564881401242048d+0*(t53+t1 237 2 4)*t55)+5.848223622634648d-1*t17*(-4*t26*t34*t51+t50+t49+ 238 3 t48+t47)+t44+t42)*wght+t58+Amat(iq,D1_RB) 239 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 240 t1 = rhoa**3.333333333333333d-1 241 t2 = t1**5.0d-1 242 t3 = 1/t2 243 t4 = 1.530901310039024d-2*t3+1.0465751434d-1 244 t5 = 7.876233178997433d-1*t3*t4+1.9269083139d-1 245 t6 = 7.876233178997433d-1*t3*t5+4.3896648423d-1 246 t7 = 1/t6 247 t8 = 1.269642451250142d+0*t2*t7+1.0d+0 248 t9 = log(t8) 249 t10 = 1/t1 250 t11 = 1.2746961887000874d-1*t10+1.0d+0 251 t12 = 2.4141993114533214d-2*t3+1.0186556948d-1 252 t13 = 7.876233178997433d-1*t12*t3+2.2308199064d-1 253 t14 = 7.876233178997433d-1*t13*t3+4.7231125998d-1 254 t15 = 1/t14 255 t16 = 1.269642451250142d+0*t15*t2+1.0d+0 256 t17 = log(t16) 257 t18 = 1.325688999052018d-1*t10+1.0d+0 258 t19 = 1.3212998810398843d-2*t3+2.9729725188d-2 259 t20 = 7.876233178997433d-1*t19*t3+1.2236585478d-1 260 t21 = 7.876233178997433d-1*t20*t3+3.497952466d-1 261 t22 = 1/t21 262 t23 = 1.269642451250142d+0*t2*t22+1.0d+0 263 t24 = log(t23) 264 t25 = 6.901399211255826d-2*t10+1.0d+0 265 t26 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(6. 266 1 21814d-2*t17*t18-3.10907d-2*t11*t9)-3.37738d-2*t24*t25)+3 267 2 .37738d-2*t24*t25)-6.21814d-2*t17*t18 268 t27 = 1/rhoa**1.3333333333333333d+0 269 t28 = 1/rhoa**1.3333333333333336d+0 270 t29 = 1/t2**3 271 t30 = 1/rhoa**6.666666666666667d-1 272 t31 = 1/t16 273 t32 = 2.11607075208357d-1*t15*t3*t30-1.269642451250142d+0*t2 274 1 *(7.876233178997433d-1*t3*(-1.3127055298329054d-1*t12*t29 275 2 *t30-3.169132786263567d-3*t28)-1.3127055298329054d-1*t13* 276 3 t29*t30)/t14**2 277 t33 = 1/t23 278 t34 = 2.11607075208357d-1*t22*t3*t30-1.269642451250142d+0*t2 279 1 *(7.876233178997433d-1*t3*(-1.3127055298329054d-1*t19*t29 280 2 *t30-1.7344776604086162d-3*t28)-1.3127055298329054d-1*t20 281 3 *t29*t30)/t21**2 282 fnc(iq) = 1.0d+0*rhoa*t26*wght+fnc(iq) 283 Amat(iq,D1_RA) = 1.0d+0*rhoa*(5.848223622634643d-1*(1.0d+0*( 284 1 1.709920934161365d+0*(1.3210398931339265d-3*t27*t9-3.1090 285 2 7d-2*t11*(2.11607075208357d-1*t3*t30*t7-1.269642451250142 286 3 d+0*t2*(7.876233178997433d-1*t3*(-1.3127055298329054d-1*t 287 4 29*t30*t4-2.0096226153166658d-3*t28)-1.3127055298329054d- 288 5 1*t29*t30*t5)/t6**2)/t8+6.21814d-2*t18*t31*t32-2.74777326 289 6 4188438d-3*t17*t27)-3.37738d-2*t25*t33*t34+7.769549222703 290 7 733d-4*t24*t27)+3.37738d-2*t25*t33*t34-7.769549222703733d 291 8 -4*t24*t27)-6.21814d-2*t18*t31*t32+2.747773264188438d-3*t 292 9 17*t27)*wght+1.0d+0*t26*wght+Amat(iq,D1_RA) 293 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 294 t1 = rhob**3.333333333333333d-1 295 t2 = t1**5.0d-1 296 t3 = 1/t2 297 t4 = 1.530901310039024d-2*t3+1.0465751434d-1 298 t5 = 7.876233178997433d-1*t3*t4+1.9269083139d-1 299 t6 = 7.876233178997433d-1*t3*t5+4.3896648423d-1 300 t7 = 1/t6 301 t8 = 1.269642451250142d+0*t2*t7+1.0d+0 302 t9 = log(t8) 303 t10 = 1/t1 304 t11 = 1.2746961887000874d-1*t10+1.0d+0 305 t12 = 2.4141993114533214d-2*t3+1.0186556948d-1 306 t13 = 7.876233178997433d-1*t12*t3+2.2308199064d-1 307 t14 = 7.876233178997433d-1*t13*t3+4.7231125998d-1 308 t15 = 1/t14 309 t16 = 1.269642451250142d+0*t15*t2+1.0d+0 310 t17 = log(t16) 311 t18 = 1.325688999052018d-1*t10+1.0d+0 312 t19 = 1.3212998810398843d-2*t3+2.9729725188d-2 313 t20 = 7.876233178997433d-1*t19*t3+1.2236585478d-1 314 t21 = 7.876233178997433d-1*t20*t3+3.497952466d-1 315 t22 = 1/t21 316 t23 = 1.269642451250142d+0*t2*t22+1.0d+0 317 t24 = log(t23) 318 t25 = 6.901399211255826d-2*t10+1.0d+0 319 t26 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(6. 320 1 21814d-2*t17*t18-3.10907d-2*t11*t9)-3.37738d-2*t24*t25)+3 321 2 .37738d-2*t24*t25)-6.21814d-2*t17*t18 322 t27 = 1/rhob**1.3333333333333333d+0 323 t28 = 1/rhob**1.3333333333333336d+0 324 t29 = 1/t2**3 325 t30 = 1/rhob**6.666666666666667d-1 326 t31 = 1/t16 327 t32 = 2.11607075208357d-1*t15*t3*t30-1.269642451250142d+0*t2 328 1 *(7.876233178997433d-1*t3*(-1.3127055298329054d-1*t12*t29 329 2 *t30-3.169132786263567d-3*t28)-1.3127055298329054d-1*t13* 330 3 t29*t30)/t14**2 331 t33 = 1/t23 332 t34 = 2.11607075208357d-1*t22*t3*t30-1.269642451250142d+0*t2 333 1 *(7.876233178997433d-1*t3*(-1.3127055298329054d-1*t19*t29 334 2 *t30-1.7344776604086162d-3*t28)-1.3127055298329054d-1*t20 335 3 *t29*t30)/t21**2 336 fnc(iq) = 1.0d+0*rhob*t26*wght+fnc(iq) 337 Amat(iq,D1_RB) = 1.0d+0*rhob*(5.848223622634643d-1*(1.0d+0*( 338 1 1.709920934161365d+0*(1.3210398931339265d-3*t27*t9-3.1090 339 2 7d-2*t11*(2.11607075208357d-1*t3*t30*t7-1.269642451250142 340 3 d+0*t2*(7.876233178997433d-1*t3*(-1.3127055298329054d-1*t 341 4 29*t30*t4-2.0096226153166658d-3*t28)-1.3127055298329054d- 342 5 1*t29*t30*t5)/t6**2)/t8+6.21814d-2*t18*t31*t32-2.74777326 343 6 4188438d-3*t17*t27)-3.37738d-2*t25*t33*t34+7.769549222703 344 7 733d-4*t24*t27)+3.37738d-2*t25*t33*t34-7.769549222703733d 345 8 -4*t24*t27)-6.21814d-2*t18*t31*t32+2.747773264188438d-3*t 346 9 17*t27)*wght+1.0d+0*t26*wght+Amat(iq,D1_RB) 347 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 348 endif ! ipol.eq.1 349 enddo ! iq 350 end 351C> 352C> \brief Evaluate the nwxcm_c_pw91lda functional [1] 353C> 354C> \f{eqnarray*}{ 355C> {\it t_1} &=& \rho_\beta+\rho_\alpha\\\\ 356C> {\it t_2} &=& {\it t_1}^{0.3333333333333333}\\\\ 357C> {\it t_3} &=& {{1}\over{{\it t_2}}}\\\\ 358C> {\it t_4} &=& 0.1325688999052018\,{\it t_3}+1.0\\\\ 359C> {\it t_5} &=& \sqrt{{\it t_2}}\\\\ 360C> {\it t_6} &=& {{1}\over{{\it t_5}}}\\\\ 361C> {\it t_7} &=& \log \left({{1.269642451250142\,{\it t_5}} 362C> \over{0.7876233178997433\,{\it t_6}\,\left(0.7876233178997433 363C> \,{\it t_6}\,\left(0.02414199311453321\,{\it t_6} 364C> +0.10186556948\right)+0.22308199064\right)+0.47231125998}} 365C> +1.0\right)\\\\ 366C> {\it t_8} &=& \rho_\alpha-\rho_\beta\\\\ 367C> {\it t_9} &=& {{1}\over{{\it t_1}}}\\\\ 368C> {\it t_{10}} &=& 0.06901399211255826\,{\it t_3}+1.0\\\\ 369C> {\it t_{11}} &=& \log \left({{1.269642451250142\,{ 370C> \it t_5}}\over{0.7876233178997433\,{\it t_6}\, 371C> \left(0.7876233178997433\,{\it t_6}\,\left(0.01321299881039884 372C> \,{\it t_6}+0.029729725188\right)+0.12236585478\right) 373C> +0.3497952466}}+1.0\right)\\\\ 374C> {\it t_{12}} &=& \rho_s^{0.3333333333333333}\\\\ 375C> {\it t_{13}} &=& \sqrt{{\it t_{12}}}\\\\ 376C> {\it t_{14}} &=& {{1}\over{{\it t_{13}}}}\\\\ 377C> {\it t_{15}} &=& {{1}\over{{\it t_{12}}}}\\\\ 378C> {\it t_{16}} &=& \log \left({{1.269642451250142\,{ 379C> \it t_{13}}}\over{0.7876233178997433\,\left(0.7876233178997433 380C> \,\left(0.02414199311453321\,{\it t_{14}}+0.10186556948\right) 381C> \,{\it t_{14}}+0.22308199064\right)\,{\it t_{14}} 382C> +0.47231125998}}+1.0\right)\\\\ 383C> {\it t_{17}} &=& 0.1325688999052018\,{\it t_{15}}+1.0\\\\ 384C> {\it t_{18}} &=& \log \left({{1.269642451250142\,{ 385C> \it t_{13}}}\over{0.7876233178997433\,\left(0.7876233178997433 386C> \,\left(0.01321299881039884\,{\it t_{14}} 387C> +0.029729725188\right)\,{\it t_{14}}+0.12236585478\right)\,{ 388C> \it t_{14}}+0.3497952466}}+1.0\right)\\\\ 389C> {\it t_{19}} &=& 0.06901399211255826\,{\it t_{15}}+1.0\\\\ 390C> f &=& 1.0\,{\it t_1}\,\left(0.5848223622634648\, 391C> \left(1.923661050931536\,\left({\it t_8}\,{\it t_9} 392C> +1.0\right)^{{{4}\over{3}}}+1.923661050931536\,\left(1.0-{ 393C> \it t_8}\,{\it t_9}\right)^{{{4}\over{3}}} 394C> -3.847322101863072\right)\,\left({{{\it t_8}^4\, 395C> \left(1.709920934161365\,\left(0.0621814\,{\it t_4}\,{\it t_7} 396C> -0.0310907\,\left(0.1274696188700087\,{\it t_3}+1.0\right) 397C> \,\log \left({{1.269642451250142\,{\it t_5}} 398C> \over{0.7876233178997433\,{\it t_6}\,\left(0.7876233178997433 399C> \,{\it t_6}\,\left(0.01530901310039024\,{\it t_6} 400C> +0.10465751434\right)+0.19269083139\right)+0.43896648423}} 401C> +1.0\right)\right)-0.0337738\,{\it t_{10}}\,{ 402C> \it t_{11}}\right)}\over{{\it t_1}^4}}+0.0337738\,{\it t_{10}} 403C> \,{\it t_{11}}\right)-0.0621814\,{\it t_4}\,{\it t_7}\right)\\\\ 404C> g &=& 0\\\\ 405C> G &=& 1.0\,\left(0.5848223622634643\,\left(0.0337738\,{ 406C> \it t_{18}}\,{\it t_{19}}+1.0\,\left(1.709920934161365\, 407C> \left(0.0621814\,{\it t_{16}}\,{\it t_{17}}-0.0310907 408C> \,\log \left({{1.269642451250142\,{\it t_{13}}} 409C> \over{0.7876233178997433\,\left(0.7876233178997433\, 410C> \left(0.01530901310039024\,{\it t_{14}}+0.10465751434\right) 411C> \,{\it t_{14}}+0.19269083139\right)\,{\it t_{14}} 412C> +0.43896648423}}+1.0\right)\,\left(0.1274696188700087\,{ 413C> \it t_{15}}+1.0\right)\right)-0.0337738\,{\it t_{18}}\,{ 414C> \it t_{19}}\right)\right)-0.0621814\,{\it t_{16}}\,{ 415C> \it t_{17}}\right)\,\rho_s\\\\ 416C> \f} 417C> 418C> Code generated with Maxima 5.34.0 [2] 419C> driven by autoxc [3]. 420C> 421C> ### References ### 422C> 423C> [1] JP Perdew, Y Wang, Phys.Rev.B 45, 13244 (1992) , DOI: 424C> <a href="https://doi.org/10.1103/PhysRevB.45.13244 "> 425C> 10.1103/PhysRevB.45.13244 </a> 426C> 427C> [2] Maxima, a computer algebra system, 428C> <a href="http://maxima.sourceforge.net/"> 429C> http://maxima.sourceforge.net/</a> 430C> 431C> [3] autoxc, revision 27097 2015-05-08 432C> 433 subroutine nwxcm_c_pw91lda_d2(param,tol_rho,ipol,nq,wght, 434 +rho,fnc,Amat,Amat2) 435c $Id: $ 436#ifdef NWXC_QUAD_PREC 437 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 438 integer, parameter :: rk=selected_real_kind(30) 439#else 440 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 441 integer, parameter :: rk=selected_real_kind(15) 442#endif 443#include "nwxc_param.fh" 444 double precision param(*) !< [Input] Parameters of functional 445 double precision tol_rho !< [Input] The lower limit on the density 446 integer ipol !< [Input] The number of spin channels 447 integer nq !< [Input] The number of points 448 double precision wght !< [Input] The weight of the functional 449 double precision rho(nq,*) !< [Input] The density 450 double precision fnc(nq) !< [Output] The value of the functional 451c 452c Sampling Matrices for the XC Kernel 453c 454 double precision Amat(nq,*) !< [Output] The derivative wrt rho 455c 456c Sampling Matrices for the XC Kernel 457c 458 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 459c 460 integer iq 461 double precision tmp 462 double precision rhoa,rhob 463 double precision gammaaa,gammaab,gammabb 464 double precision taua,taub 465 double precision nwxcm_heaviside 466 external nwxcm_heaviside 467CDIR$ NOVECTOR 468 do iq = 1, nq 469 if (ipol.eq.1) then 470 rhoa = 0.5d0*rho(iq,R_T) 471 if (rhoa.gt.tol_rho) then 472 t1 = rhoa**3.333333333333333d-1 473 t2 = t1**5.0d-1 474 t3 = 1/t2 475 t4 = 1/t1 476 t5 = 2.1508070719090538d-2*t3+1.0186556948d-1 477 t6 = 7.016926042943222d-1*t3*t5+2.2308199064d-1 478 t7 = 7.016926042943222d-1*t3*t6+4.7231125998d-1 479 t8 = 1/t7 480 t9 = 1.425125466450768d+0*t2*t8+1.0d+0 481 t10 = log(t9) 482 t11 = 1.0522000558389213d-1*t4+1.0d+0 483 t12 = 1/rhoa**1.3333333333333333d+0 484 t13 = 1/t9 485 t14 = 1/t7**2 486 t15 = 1/rhoa**1.3333333333333336d+0 487 t16 = 1/t2**3 488 t17 = 1/rhoa**6.666666666666667d-1 489 t18 = -5.847438369119352d-2*t16*t17*t5-1.2576711796854242d-3 490 1 *t15 491 t19 = 7.016926042943222d-1*t18*t3-5.847438369119352d-2*t16*t 492 1 17*t6 493 t20 = 1.1876045553756398d-1*t17*t3*t8-1.425125466450768d+0*t 494 1 14*t19*t2 495 t21 = 1.090454542535705d-3*t10*t12-6.21814d-2*t11*t13*t20 496 t22 = 2.0d+0*t21*wght 497 t23 = -7.269696950238034d-4*t10/rhoa**2.333333333333333d+0 498 t24 = log(1.425125466450768d+0*t2/(7.016926042943222d-1*t3*( 499 1 7.016926042943222d-1*(1.1771443702974158d-2*t3+2.97297251 500 2 88d-2)*t3+1.2236585478d-1)+3.497952466d-1)+1.0d+0) 501 t25 = 5.477644184000001d-2*t4+1.0d+0 502 t26 = 1/rhoa**2 503 t27 = 2.18090908507141d-3*t12*t13*t20 504 t28 = 6.21814d-2*t11*t20**2/t9**2 505 t29 = 1/rhoa**1.6666666666666669d+0 506 t30 = 1/t2**5 507 t31 = -6.21814d-2*t11*t13*(-3.9586818512521327d-2*t29*t3*t8- 508 1 9.896704628130328d-3*t15*t16*t8+2.850250932901536d+0*t19* 509 2 *2*t2/t7**3-1.425125466450768d+0*t14*t2*(1.46185959227983 510 3 75d-2*t15*t30*t6+1.949146123039784d-2*t16*t29*t6+7.016926 511 4 042943222d-1*t3*(1.4618595922798375d-2*t15*t30*t5+1.94914 512 5 6123039784d-2*t16*t29*t5+9.432533847640683d-4/rhoa**2.333 513 6 3333333333334d+0)-1.1694876738238703d-1*t16*t17*t18)-2.37 514 7 52091107512796d-1*t14*t17*t19*t3) 515 fnc(iq) = fnc(iq)-1.243628d-1*rhoa*log(1.4251254664507676d+0 516 1 *t2/(7.016926042943223d-1*t3*(7.016926042943223d-1*(2.150 517 2 807071909054d-2*t3+1.0186556948d-1)*t3+2.2308199064d-1)+4 518 3 .7231125998d-1)+1.0d+0)*(1.0522000558389215d-1*t4+1.0d+0) 519 4 *wght 520 Amat(iq,D1_RA) = 2.0d+0*rhoa*t21*wght-6.21814d-2*t10*t11*wgh 521 1 t+Amat(iq,D1_RA) 522 Amat2(iq,D2_RA_RA) = 2.0d+0*rhoa*(t31+t28+t27+8.443450000000 523 1 001d-3*t24*t25*t26+t23)*wght+t22+Amat2(iq,D2_RA_RA) 524 Amat2(iq,D2_RA_RB) = 2.0d+0*rhoa*(t31+t28+t27-8.44345d-3*t24 525 1 *t25*t26+t23)*wght+t22+Amat2(iq,D2_RA_RB) 526 endif ! rhoa.gt.tol_rho 527 else ! ipol.eq.1 528 rhoa = rho(iq,R_A) 529 rhob = rho(iq,R_B) 530 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 531 t1 = rhob+rhoa 532 t2 = t1**3.333333333333333d-1 533 t3 = 1/t2 534 t4 = 1.325688999052018d-1*t3+1.0d+0 535 t5 = t2**5.0d-1 536 t6 = 1/t5 537 t7 = 2.4141993114533214d-2*t6+1.0186556948d-1 538 t8 = 7.876233178997433d-1*t6*t7+2.2308199064d-1 539 t9 = 7.876233178997433d-1*t6*t8+4.7231125998d-1 540 t10 = 1/t9 541 t11 = 1.269642451250142d+0*t10*t5+1.0d+0 542 t12 = log(t11) 543 t13 = rhoa-rhob 544 t14 = 1/t1 545 t15 = 1.0d+0-t13*t14 546 t16 = t13*t14+1.0d+0 547 t17 = 1.923661050931536d+0*t16**1.3333333333333333d+0+1.9236 548 1 61050931536d+0*t15**1.3333333333333333d+0-3.8473221018630 549 2 72d+0 550 t18 = 6.901399211255826d-2*t3+1.0d+0 551 t19 = 1.3212998810398843d-2*t6+2.9729725188d-2 552 t20 = 7.876233178997433d-1*t19*t6+1.2236585478d-1 553 t21 = 7.876233178997433d-1*t20*t6+3.497952466d-1 554 t22 = 1/t21 555 t23 = 1.269642451250142d+0*t22*t5+1.0d+0 556 t24 = log(t23) 557 t25 = t13**4 558 t26 = 1/t1**4 559 t27 = 1.2746961887000874d-1*t3+1.0d+0 560 t28 = 1.530901310039024d-2*t6+1.0465751434d-1 561 t29 = 7.876233178997433d-1*t28*t6+1.9269083139d-1 562 t30 = 7.876233178997433d-1*t29*t6+4.3896648423d-1 563 t31 = 1/t30 564 t32 = 1.269642451250142d+0*t31*t5+1.0d+0 565 t33 = log(t32) 566 t34 = 1.709920934161365d+0*(6.21814d-2*t12*t4-3.10907d-2*t27 567 1 *t33)-3.37738d-2*t18*t24 568 t35 = t25*t26*t34+3.37738d-2*t18*t24 569 t36 = 5.848223622634648d-1*t17*t35-6.21814d-2*t12*t4 570 t37 = 1/t1**1.3333333333333336d+0 571 t38 = 1/t1**6.666666666666667d-1 572 t39 = 1/t5**3 573 t40 = -1.3127055298329054d-1*t38*t39*t7-3.169132786263567d-3 574 1 *t37 575 t41 = 7.876233178997433d-1*t40*t6-1.3127055298329054d-1*t38* 576 1 t39*t8 577 t42 = 1/t9**2 578 t43 = 2.11607075208357d-1*t10*t38*t6-1.269642451250142d+0*t4 579 1 1*t42*t5 580 t44 = 1/t11 581 t45 = -6.21814d-2*t4*t43*t44 582 t46 = 1/t1**1.3333333333333333d+0 583 t47 = 2.747773264188438d-3*t12*t46 584 t48 = -1.3127055298329054d-1*t19*t38*t39-1.7344776604086162d 585 1 -3*t37 586 t49 = 7.876233178997433d-1*t48*t6-1.3127055298329054d-1*t20* 587 1 t38*t39 588 t50 = 1/t21**2 589 t51 = 2.11607075208357d-1*t22*t38*t6-1.269642451250142d+0*t4 590 1 9*t5*t50 591 t52 = 1/t23 592 t53 = 3.37738d-2*t18*t51*t52 593 t54 = -7.769549222703733d-4*t24*t46 594 t55 = -1.3127055298329054d-1*t28*t38*t39-2.0096226153166658d 595 1 -3*t37 596 t56 = 7.876233178997433d-1*t55*t6-1.3127055298329054d-1*t29* 597 1 t38*t39 598 t57 = 1/t30**2 599 t58 = 2.11607075208357d-1*t31*t38*t6-1.269642451250142d+0*t5 600 1 *t56*t57 601 t59 = 1/t32 602 t60 = 1.709920934161365d+0*(-3.10907d-2*t27*t58*t59+1.321039 603 1 8931339265d-3*t33*t46-2.747773264188438d-3*t12*t46+6.2181 604 2 4d-2*t4*t43*t44)-3.37738d-2*t18*t51*t52+7.769549222703733 605 3 d-4*t24*t46 606 t61 = t25*t26*t60 607 t62 = 1/t1**5 608 t63 = -4*t25*t34*t62 609 t64 = t13**3 610 t65 = 4*t26*t34*t64+t63+t61+t54+t53 611 t66 = 1/t1**2 612 t67 = t13*t66 613 t68 = -t14 614 t69 = t68+t67 615 t70 = t15**3.333333333333333d-1 616 t71 = -t13*t66 617 t72 = t71+t14 618 t73 = t16**3.333333333333333d-1 619 t74 = 2.564881401242048d+0*t72*t73+2.564881401242048d+0*t69* 620 1 t70 621 t75 = 5.848223622634648d-1*t35*t74+5.848223622634648d-1*t17* 622 1 t65+t47+t45 623 t76 = 1.0d+0*t36*wght 624 t77 = -4*t26*t34*t64+t63+t61+t54+t53 625 t78 = t67+t14 626 t79 = t71+t68 627 t80 = 2.564881401242048d+0*t73*t79+2.564881401242048d+0*t70* 628 1 t78 629 t81 = 5.848223622634648d-1*t35*t80+5.848223622634648d-1*t17* 630 1 t77+t47+t45 631 t82 = t43**2 632 t83 = 1/t11**2 633 t84 = 6.21814d-2*t4*t82*t83 634 t85 = 1/t1**2.3333333333333334d+0 635 t86 = 1/t5**5 636 t87 = 1/t1**1.6666666666666669d+0 637 t88 = 2.539284902500284d+0*t41**2*t5/t9**3-1.269642451250142 638 1 d+0*t42*t5*(7.876233178997433d-1*t6*(8.751370198886037d-2 639 2 *t39*t7*t87+6.563527649164527d-2*t37*t7*t86+4.75369917939 640 3 5351d-3*t85)+8.751370198886037d-2*t39*t8*t87+6.5635276491 641 4 64527d-2*t37*t8*t86-2.625411059665811d-1*t38*t39*t40)-1.4 642 5 107138347223802d-1*t10*t6*t87-4.23214150416714d-1*t38*t41 643 6 *t42*t6-3.52678458680595d-2*t10*t37*t39 644 t89 = -6.21814d-2*t4*t44*t88 645 t90 = 5.495546528376876d-3*t43*t44*t46 646 t91 = 1/t1**2.333333333333333d+0 647 t92 = -3.663697685584584d-3*t12*t91 648 t93 = t51**2 649 t94 = 1/t23**2 650 t95 = -3.37738d-2*t18*t93*t94 651 t96 = -1.269642451250142d+0*t5*t50*(7.876233178997433d-1*t6* 652 1 (8.751370198886037d-2*t19*t39*t87+6.563527649164527d-2*t1 653 2 9*t37*t86+2.601716490612924d-3*t85)+8.751370198886037d-2* 654 3 t20*t39*t87+6.563527649164527d-2*t20*t37*t86-2.6254110596 655 4 65811d-1*t38*t39*t48)-1.4107138347223802d-1*t22*t6*t87-4. 656 5 23214150416714d-1*t38*t49*t50*t6+2.539284902500284d+0*t49 657 6 **2*t5/t21**3-3.52678458680595d-2*t22*t37*t39 658 t97 = 3.37738d-2*t18*t52*t96 659 t98 = -1.5539098445407465d-3*t46*t51*t52 660 t99 = 1.0359398963604977d-3*t24*t91 661 t100 = t25*t26*(-3.37738d-2*t18*t52*t96+3.37738d-2*t18*t93*t 662 1 94+1.709920934161365d+0*(-1.7613865241785687d-3*t33*t91+3 663 2 .663697685584584d-3*t12*t91+6.21814d-2*t4*t44*t88-3.10907 664 3 d-2*t27*t59*(-1.269642451250142d+0*t5*t57*(7.876233178997 665 4 433d-1*t6*(8.751370198886037d-2*t28*t39*t87+6.56352764916 666 5 4527d-2*t28*t37*t86+3.0144339229749983d-3*t85)+8.75137019 667 6 8886037d-2*t29*t39*t87+6.563527649164527d-2*t29*t37*t86-2 668 7 .625411059665811d-1*t38*t39*t55)-1.4107138347223802d-1*t3 669 8 1*t6*t87-4.23214150416714d-1*t38*t56*t57*t6+2.53928490250 670 9 0284d+0*t5*t56**2/t30**3-3.52678458680595d-2*t31*t37*t39) 671 : -6.21814d-2*t4*t82*t83+2.642079786267853d-3*t46*t58*t59+3 672 ; .10907d-2*t27*t58**2/t32**2-5.495546528376876d-3*t43*t44* 673 < t46)-1.0359398963604977d-3*t24*t91+1.5539098445407465d-3* 674 = t46*t51*t52) 675 t101 = -8*t25*t60*t62 676 t102 = 20*t25*t34/t1**6 677 t103 = t13**2 678 t104 = 12*t103*t26*t34 679 t105 = 1/t15**6.666666666666666d-1 680 t106 = 1/t1**3 681 t107 = -2*t106*t13 682 t108 = 2*t66 683 t109 = 1/t16**6.666666666666666d-1 684 t110 = 2*t106*t13 685 t111 = -2*t66 686 fnc(iq) = 1.0d+0*t1*t36*wght+fnc(iq) 687 Amat(iq,D1_RA) = 1.0d+0*t1*t75*wght+t76+Amat(iq,D1_RA) 688 Amat(iq,D1_RB) = 1.0d+0*t1*t81*wght+t76+Amat(iq,D1_RB) 689 Amat2(iq,D2_RA_RA) = 1.0d+0*t1*(5.848223622634648d-1*t17*(t9 690 1 9+t98+t97+t95-32*t34*t62*t64+8*t26*t60*t64+t104+t102+t101 691 2 +t100)+t92+t90+t89+t84+1.1696447245269297d+0*t65*t74+5.84 692 3 8223622634648d-1*t35*(2.564881401242048d+0*(t111+t110)*t7 693 4 3+8.549604670806825d-1*t109*t72**2+2.564881401242048d+0*( 694 5 t108+t107)*t70+8.549604670806825d-1*t105*t69**2))*wght+2. 695 6 0d+0*t75*wght+Amat2(iq,D2_RA_RA) 696 Amat2(iq,D2_RA_RB) = 1.0d+0*t1*(5.848223622634648d-1*t17*(t9 697 1 9+t98+t97+t95-12*t103*t26*t34+t102+t101+t100)+t92+t90+t89 698 2 +t84+5.848223622634648d-1*t65*t80+5.848223622634648d-1*t3 699 3 5*(8.549604670806825d-1*t109*t72*t79+8.549604670806825d-1 700 4 *t105*t69*t78+5.129762802484096d+0*t106*t13*t73-5.1297628 701 5 02484096d+0*t106*t13*t70)+5.848223622634648d-1*t74*t77)*w 702 6 ght+1.0d+0*t81*wght+1.0d+0*t75*wght+Amat2(iq,D2_RA_RB) 703 Amat2(iq,D2_RB_RB) = 1.0d+0*t1*(5.848223622634648d-1*t17*(t9 704 1 9+t98+t97+t95+32*t34*t62*t64-8*t26*t60*t64+t104+t102+t101 705 2 +t100)+t92+t90+t89+t84+1.1696447245269297d+0*t77*t80+5.84 706 3 8223622634648d-1*t35*(8.549604670806825d-1*t109*t79**2+8. 707 4 549604670806825d-1*t105*t78**2+2.564881401242048d+0*(t110 708 5 +t108)*t73+2.564881401242048d+0*(t111+t107)*t70))*wght+2. 709 6 0d+0*t81*wght+Amat2(iq,D2_RB_RB) 710 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 711 t1 = rhoa**3.333333333333333d-1 712 t2 = t1**5.0d-1 713 t3 = 1/t2 714 t4 = 1.530901310039024d-2*t3+1.0465751434d-1 715 t5 = 7.876233178997433d-1*t3*t4+1.9269083139d-1 716 t6 = 7.876233178997433d-1*t3*t5+4.3896648423d-1 717 t7 = 1/t6 718 t8 = 1.269642451250142d+0*t2*t7+1.0d+0 719 t9 = log(t8) 720 t10 = 1/t1 721 t11 = 1.2746961887000874d-1*t10+1.0d+0 722 t12 = 2.4141993114533214d-2*t3+1.0186556948d-1 723 t13 = 7.876233178997433d-1*t12*t3+2.2308199064d-1 724 t14 = 7.876233178997433d-1*t13*t3+4.7231125998d-1 725 t15 = 1/t14 726 t16 = 1.269642451250142d+0*t15*t2+1.0d+0 727 t17 = log(t16) 728 t18 = 1.325688999052018d-1*t10+1.0d+0 729 t19 = 1.3212998810398843d-2*t3+2.9729725188d-2 730 t20 = 7.876233178997433d-1*t19*t3+1.2236585478d-1 731 t21 = 7.876233178997433d-1*t20*t3+3.497952466d-1 732 t22 = 1/t21 733 t23 = 1.269642451250142d+0*t2*t22+1.0d+0 734 t24 = log(t23) 735 t25 = 6.901399211255826d-2*t10+1.0d+0 736 t26 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(6. 737 1 21814d-2*t17*t18-3.10907d-2*t11*t9)-3.37738d-2*t24*t25)+3 738 2 .37738d-2*t24*t25)-6.21814d-2*t17*t18 739 t27 = 1/rhoa**1.3333333333333333d+0 740 t28 = 1/t8 741 t29 = 1/t6**2 742 t30 = 1/rhoa**1.3333333333333336d+0 743 t31 = 1/t2**3 744 t32 = 1/rhoa**6.666666666666667d-1 745 t33 = -1.3127055298329054d-1*t31*t32*t4-2.0096226153166658d- 746 1 3*t30 747 t34 = 7.876233178997433d-1*t3*t33-1.3127055298329054d-1*t31* 748 1 t32*t5 749 t35 = 2.11607075208357d-1*t3*t32*t7-1.269642451250142d+0*t2* 750 1 t29*t34 751 t36 = 1/t16 752 t37 = 1/t14**2 753 t38 = -1.3127055298329054d-1*t12*t31*t32-3.169132786263567d- 754 1 3*t30 755 t39 = 7.876233178997433d-1*t3*t38-1.3127055298329054d-1*t13* 756 1 t31*t32 757 t40 = 2.11607075208357d-1*t15*t3*t32-1.269642451250142d+0*t2 758 1 *t37*t39 759 t41 = 1/t23 760 t42 = 1/t21**2 761 t43 = -1.3127055298329054d-1*t19*t31*t32-1.7344776604086162d 762 1 -3*t30 763 t44 = 7.876233178997433d-1*t3*t43-1.3127055298329054d-1*t20* 764 1 t31*t32 765 t45 = 2.11607075208357d-1*t22*t3*t32-1.269642451250142d+0*t2 766 1 *t42*t44 767 t46 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(1. 768 1 3210398931339265d-3*t27*t9+6.21814d-2*t18*t36*t40-3.10907 769 2 d-2*t11*t28*t35-2.747773264188438d-3*t17*t27)-3.37738d-2* 770 3 t25*t41*t45+7.769549222703733d-4*t24*t27)+3.37738d-2*t25* 771 4 t41*t45-7.769549222703733d-4*t24*t27)-6.21814d-2*t18*t36* 772 5 t40+2.747773264188438d-3*t17*t27 773 t47 = 1/rhoa**2.333333333333333d+0 774 t48 = 1/rhoa**2.3333333333333334d+0 775 t49 = 1/rhoa**1.6666666666666669d+0 776 t50 = 1/t2**5 777 t51 = 1/t16**2 778 t52 = t40**2 779 t53 = -1.269642451250142d+0*t2*t37*(7.876233178997433d-1*t3* 780 1 (6.563527649164527d-2*t12*t30*t50+8.751370198886037d-2*t1 781 2 2*t31*t49+4.753699179395351d-3*t48)+6.563527649164527d-2* 782 3 t13*t30*t50+8.751370198886037d-2*t13*t31*t49-2.6254110596 783 4 65811d-1*t31*t32*t38)-1.4107138347223802d-1*t15*t3*t49+2. 784 5 539284902500284d+0*t2*t39**2/t14**3-4.23214150416714d-1*t 785 6 3*t32*t37*t39-3.52678458680595d-2*t15*t30*t31 786 t54 = 1/t23**2 787 t55 = t45**2 788 t56 = -1.269642451250142d+0*t2*t42*(7.876233178997433d-1*t3* 789 1 (6.563527649164527d-2*t19*t30*t50+8.751370198886037d-2*t1 790 2 9*t31*t49+2.601716490612924d-3*t48)+6.563527649164527d-2* 791 3 t20*t30*t50+8.751370198886037d-2*t20*t31*t49-2.6254110596 792 4 65811d-1*t31*t32*t43)-1.4107138347223802d-1*t22*t3*t49+2. 793 5 539284902500284d+0*t2*t44**2/t21**3-4.23214150416714d-1*t 794 6 3*t32*t42*t44-3.52678458680595d-2*t22*t30*t31 795 fnc(iq) = 1.0d+0*rhoa*t26*wght+fnc(iq) 796 Amat(iq,D1_RA) = 1.0d+0*rhoa*t46*wght+1.0d+0*t26*wght+Amat(i 797 1 q,D1_RA) 798 Amat2(iq,D2_RA_RA) = 1.0d+0*rhoa*(5.848223622634643d-1*(1.0d 799 1 +0*(1.709920934161365d+0*(-1.7613865241785687d-3*t47*t9+3 800 2 .10907d-2*t11*t35**2/t8**2-3.10907d-2*t11*t28*(-1.4107138 801 3 347223802d-1*t3*t49*t7-3.52678458680595d-2*t30*t31*t7+2.5 802 4 39284902500284d+0*t2*t34**2/t6**3-1.269642451250142d+0*t2 803 5 *t29*(7.876233178997433d-1*t3*(6.563527649164527d-2*t30*t 804 6 4*t50+8.751370198886037d-2*t31*t4*t49+3.0144339229749983d 805 7 -3*t48)+6.563527649164527d-2*t30*t5*t50+8.751370198886037 806 8 d-2*t31*t49*t5-2.625411059665811d-1*t31*t32*t33)-4.232141 807 9 50416714d-1*t29*t3*t32*t34)+6.21814d-2*t18*t36*t53-6.2181 808 : 4d-2*t18*t51*t52+3.663697685584584d-3*t17*t47-5.495546528 809 ; 376876d-3*t27*t36*t40+2.642079786267853d-3*t27*t28*t35)-3 810 < .37738d-2*t25*t41*t56+3.37738d-2*t25*t54*t55-1.0359398963 811 = 604977d-3*t24*t47+1.5539098445407465d-3*t27*t41*t45)+3.37 812 > 738d-2*t25*t41*t56-3.37738d-2*t25*t54*t55+1.0359398963604 813 ? 977d-3*t24*t47-1.5539098445407465d-3*t27*t41*t45)-6.21814 814 @ d-2*t18*t36*t53+6.21814d-2*t18*t51*t52-3.663697685584584d 815 1 -3*t17*t47+5.495546528376876d-3*t27*t36*t40)*wght+2.0d+0* 816 2 t46*wght+Amat2(iq,D2_RA_RA) 817 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 818 t1 = rhob**3.333333333333333d-1 819 t2 = t1**5.0d-1 820 t3 = 1/t2 821 t4 = 1.530901310039024d-2*t3+1.0465751434d-1 822 t5 = 7.876233178997433d-1*t3*t4+1.9269083139d-1 823 t6 = 7.876233178997433d-1*t3*t5+4.3896648423d-1 824 t7 = 1/t6 825 t8 = 1.269642451250142d+0*t2*t7+1.0d+0 826 t9 = log(t8) 827 t10 = 1/t1 828 t11 = 1.2746961887000874d-1*t10+1.0d+0 829 t12 = 2.4141993114533214d-2*t3+1.0186556948d-1 830 t13 = 7.876233178997433d-1*t12*t3+2.2308199064d-1 831 t14 = 7.876233178997433d-1*t13*t3+4.7231125998d-1 832 t15 = 1/t14 833 t16 = 1.269642451250142d+0*t15*t2+1.0d+0 834 t17 = log(t16) 835 t18 = 1.325688999052018d-1*t10+1.0d+0 836 t19 = 1.3212998810398843d-2*t3+2.9729725188d-2 837 t20 = 7.876233178997433d-1*t19*t3+1.2236585478d-1 838 t21 = 7.876233178997433d-1*t20*t3+3.497952466d-1 839 t22 = 1/t21 840 t23 = 1.269642451250142d+0*t2*t22+1.0d+0 841 t24 = log(t23) 842 t25 = 6.901399211255826d-2*t10+1.0d+0 843 t26 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(6. 844 1 21814d-2*t17*t18-3.10907d-2*t11*t9)-3.37738d-2*t24*t25)+3 845 2 .37738d-2*t24*t25)-6.21814d-2*t17*t18 846 t27 = 1/rhob**1.3333333333333333d+0 847 t28 = 1/t8 848 t29 = 1/t6**2 849 t30 = 1/rhob**1.3333333333333336d+0 850 t31 = 1/t2**3 851 t32 = 1/rhob**6.666666666666667d-1 852 t33 = -1.3127055298329054d-1*t31*t32*t4-2.0096226153166658d- 853 1 3*t30 854 t34 = 7.876233178997433d-1*t3*t33-1.3127055298329054d-1*t31* 855 1 t32*t5 856 t35 = 2.11607075208357d-1*t3*t32*t7-1.269642451250142d+0*t2* 857 1 t29*t34 858 t36 = 1/t16 859 t37 = 1/t14**2 860 t38 = -1.3127055298329054d-1*t12*t31*t32-3.169132786263567d- 861 1 3*t30 862 t39 = 7.876233178997433d-1*t3*t38-1.3127055298329054d-1*t13* 863 1 t31*t32 864 t40 = 2.11607075208357d-1*t15*t3*t32-1.269642451250142d+0*t2 865 1 *t37*t39 866 t41 = 1/t23 867 t42 = 1/t21**2 868 t43 = -1.3127055298329054d-1*t19*t31*t32-1.7344776604086162d 869 1 -3*t30 870 t44 = 7.876233178997433d-1*t3*t43-1.3127055298329054d-1*t20* 871 1 t31*t32 872 t45 = 2.11607075208357d-1*t22*t3*t32-1.269642451250142d+0*t2 873 1 *t42*t44 874 t46 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(1. 875 1 3210398931339265d-3*t27*t9+6.21814d-2*t18*t36*t40-3.10907 876 2 d-2*t11*t28*t35-2.747773264188438d-3*t17*t27)-3.37738d-2* 877 3 t25*t41*t45+7.769549222703733d-4*t24*t27)+3.37738d-2*t25* 878 4 t41*t45-7.769549222703733d-4*t24*t27)-6.21814d-2*t18*t36* 879 5 t40+2.747773264188438d-3*t17*t27 880 t47 = 1/rhob**2.333333333333333d+0 881 t48 = 1/rhob**2.3333333333333334d+0 882 t49 = 1/rhob**1.6666666666666669d+0 883 t50 = 1/t2**5 884 t51 = 1/t16**2 885 t52 = t40**2 886 t53 = -1.269642451250142d+0*t2*t37*(7.876233178997433d-1*t3* 887 1 (6.563527649164527d-2*t12*t30*t50+8.751370198886037d-2*t1 888 2 2*t31*t49+4.753699179395351d-3*t48)+6.563527649164527d-2* 889 3 t13*t30*t50+8.751370198886037d-2*t13*t31*t49-2.6254110596 890 4 65811d-1*t31*t32*t38)-1.4107138347223802d-1*t15*t3*t49+2. 891 5 539284902500284d+0*t2*t39**2/t14**3-4.23214150416714d-1*t 892 6 3*t32*t37*t39-3.52678458680595d-2*t15*t30*t31 893 t54 = 1/t23**2 894 t55 = t45**2 895 t56 = -1.269642451250142d+0*t2*t42*(7.876233178997433d-1*t3* 896 1 (6.563527649164527d-2*t19*t30*t50+8.751370198886037d-2*t1 897 2 9*t31*t49+2.601716490612924d-3*t48)+6.563527649164527d-2* 898 3 t20*t30*t50+8.751370198886037d-2*t20*t31*t49-2.6254110596 899 4 65811d-1*t31*t32*t43)-1.4107138347223802d-1*t22*t3*t49+2. 900 5 539284902500284d+0*t2*t44**2/t21**3-4.23214150416714d-1*t 901 6 3*t32*t42*t44-3.52678458680595d-2*t22*t30*t31 902 fnc(iq) = 1.0d+0*rhob*t26*wght+fnc(iq) 903 Amat(iq,D1_RB) = 1.0d+0*rhob*t46*wght+1.0d+0*t26*wght+Amat(i 904 1 q,D1_RB) 905 Amat2(iq,D2_RB_RB) = 1.0d+0*rhob*(5.848223622634643d-1*(1.0d 906 1 +0*(1.709920934161365d+0*(-1.7613865241785687d-3*t47*t9+3 907 2 .10907d-2*t11*t35**2/t8**2-3.10907d-2*t11*t28*(-1.4107138 908 3 347223802d-1*t3*t49*t7-3.52678458680595d-2*t30*t31*t7+2.5 909 4 39284902500284d+0*t2*t34**2/t6**3-1.269642451250142d+0*t2 910 5 *t29*(7.876233178997433d-1*t3*(6.563527649164527d-2*t30*t 911 6 4*t50+8.751370198886037d-2*t31*t4*t49+3.0144339229749983d 912 7 -3*t48)+6.563527649164527d-2*t30*t5*t50+8.751370198886037 913 8 d-2*t31*t49*t5-2.625411059665811d-1*t31*t32*t33)-4.232141 914 9 50416714d-1*t29*t3*t32*t34)+6.21814d-2*t18*t36*t53-6.2181 915 : 4d-2*t18*t51*t52+3.663697685584584d-3*t17*t47-5.495546528 916 ; 376876d-3*t27*t36*t40+2.642079786267853d-3*t27*t28*t35)-3 917 < .37738d-2*t25*t41*t56+3.37738d-2*t25*t54*t55-1.0359398963 918 = 604977d-3*t24*t47+1.5539098445407465d-3*t27*t41*t45)+3.37 919 > 738d-2*t25*t41*t56-3.37738d-2*t25*t54*t55+1.0359398963604 920 ? 977d-3*t24*t47-1.5539098445407465d-3*t27*t41*t45)-6.21814 921 @ d-2*t18*t36*t53+6.21814d-2*t18*t51*t52-3.663697685584584d 922 1 -3*t17*t47+5.495546528376876d-3*t27*t36*t40)*wght+2.0d+0* 923 2 t46*wght+Amat2(iq,D2_RB_RB) 924 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 925 endif ! ipol.eq.1 926 enddo ! iq 927 end 928C> 929C> \brief Evaluate the nwxcm_c_pw91lda functional [1] 930C> 931C> \f{eqnarray*}{ 932C> {\it t_1} &=& \rho_\beta+\rho_\alpha\\\\ 933C> {\it t_2} &=& {\it t_1}^{0.3333333333333333}\\\\ 934C> {\it t_3} &=& {{1}\over{{\it t_2}}}\\\\ 935C> {\it t_4} &=& 0.1325688999052018\,{\it t_3}+1.0\\\\ 936C> {\it t_5} &=& \sqrt{{\it t_2}}\\\\ 937C> {\it t_6} &=& {{1}\over{{\it t_5}}}\\\\ 938C> {\it t_7} &=& \log \left({{1.269642451250142\,{\it t_5}} 939C> \over{0.7876233178997433\,{\it t_6}\,\left(0.7876233178997433 940C> \,{\it t_6}\,\left(0.02414199311453321\,{\it t_6} 941C> +0.10186556948\right)+0.22308199064\right)+0.47231125998}} 942C> +1.0\right)\\\\ 943C> {\it t_8} &=& \rho_\alpha-\rho_\beta\\\\ 944C> {\it t_9} &=& {{1}\over{{\it t_1}}}\\\\ 945C> {\it t_{10}} &=& 0.06901399211255826\,{\it t_3}+1.0\\\\ 946C> {\it t_{11}} &=& \log \left({{1.269642451250142\,{ 947C> \it t_5}}\over{0.7876233178997433\,{\it t_6}\, 948C> \left(0.7876233178997433\,{\it t_6}\,\left(0.01321299881039884 949C> \,{\it t_6}+0.029729725188\right)+0.12236585478\right) 950C> +0.3497952466}}+1.0\right)\\\\ 951C> {\it t_{12}} &=& \rho_s^{0.3333333333333333}\\\\ 952C> {\it t_{13}} &=& \sqrt{{\it t_{12}}}\\\\ 953C> {\it t_{14}} &=& {{1}\over{{\it t_{13}}}}\\\\ 954C> {\it t_{15}} &=& {{1}\over{{\it t_{12}}}}\\\\ 955C> {\it t_{16}} &=& \log \left({{1.269642451250142\,{ 956C> \it t_{13}}}\over{0.7876233178997433\,\left(0.7876233178997433 957C> \,\left(0.02414199311453321\,{\it t_{14}}+0.10186556948\right) 958C> \,{\it t_{14}}+0.22308199064\right)\,{\it t_{14}} 959C> +0.47231125998}}+1.0\right)\\\\ 960C> {\it t_{17}} &=& 0.1325688999052018\,{\it t_{15}}+1.0\\\\ 961C> {\it t_{18}} &=& \log \left({{1.269642451250142\,{ 962C> \it t_{13}}}\over{0.7876233178997433\,\left(0.7876233178997433 963C> \,\left(0.01321299881039884\,{\it t_{14}} 964C> +0.029729725188\right)\,{\it t_{14}}+0.12236585478\right)\,{ 965C> \it t_{14}}+0.3497952466}}+1.0\right)\\\\ 966C> {\it t_{19}} &=& 0.06901399211255826\,{\it t_{15}}+1.0\\\\ 967C> f &=& 1.0\,{\it t_1}\,\left(0.5848223622634648\, 968C> \left(1.923661050931536\,\left({\it t_8}\,{\it t_9} 969C> +1.0\right)^{{{4}\over{3}}}+1.923661050931536\,\left(1.0-{ 970C> \it t_8}\,{\it t_9}\right)^{{{4}\over{3}}} 971C> -3.847322101863072\right)\,\left({{{\it t_8}^4\, 972C> \left(1.709920934161365\,\left(0.0621814\,{\it t_4}\,{\it t_7} 973C> -0.0310907\,\left(0.1274696188700087\,{\it t_3}+1.0\right) 974C> \,\log \left({{1.269642451250142\,{\it t_5}} 975C> \over{0.7876233178997433\,{\it t_6}\,\left(0.7876233178997433 976C> \,{\it t_6}\,\left(0.01530901310039024\,{\it t_6} 977C> +0.10465751434\right)+0.19269083139\right)+0.43896648423}} 978C> +1.0\right)\right)-0.0337738\,{\it t_{10}}\,{ 979C> \it t_{11}}\right)}\over{{\it t_1}^4}}+0.0337738\,{\it t_{10}} 980C> \,{\it t_{11}}\right)-0.0621814\,{\it t_4}\,{\it t_7}\right)\\\\ 981C> g &=& 0\\\\ 982C> G &=& 1.0\,\left(0.5848223622634643\,\left(0.0337738\,{ 983C> \it t_{18}}\,{\it t_{19}}+1.0\,\left(1.709920934161365\, 984C> \left(0.0621814\,{\it t_{16}}\,{\it t_{17}}-0.0310907 985C> \,\log \left({{1.269642451250142\,{\it t_{13}}} 986C> \over{0.7876233178997433\,\left(0.7876233178997433\, 987C> \left(0.01530901310039024\,{\it t_{14}}+0.10465751434\right) 988C> \,{\it t_{14}}+0.19269083139\right)\,{\it t_{14}} 989C> +0.43896648423}}+1.0\right)\,\left(0.1274696188700087\,{ 990C> \it t_{15}}+1.0\right)\right)-0.0337738\,{\it t_{18}}\,{ 991C> \it t_{19}}\right)\right)-0.0621814\,{\it t_{16}}\,{ 992C> \it t_{17}}\right)\,\rho_s\\\\ 993C> \f} 994C> 995C> Code generated with Maxima 5.34.0 [2] 996C> driven by autoxc [3]. 997C> 998C> ### References ### 999C> 1000C> [1] JP Perdew, Y Wang, Phys.Rev.B 45, 13244 (1992) , DOI: 1001C> <a href="https://doi.org/10.1103/PhysRevB.45.13244 "> 1002C> 10.1103/PhysRevB.45.13244 </a> 1003C> 1004C> [2] Maxima, a computer algebra system, 1005C> <a href="http://maxima.sourceforge.net/"> 1006C> http://maxima.sourceforge.net/</a> 1007C> 1008C> [3] autoxc, revision 27097 2015-05-08 1009C> 1010 subroutine nwxcm_c_pw91lda_d3(param,tol_rho,ipol,nq,wght, 1011 +rho,fnc,Amat,Amat2,Amat3) 1012c $Id: $ 1013#ifdef NWXC_QUAD_PREC 1014 implicit real(kind=selected_real_kind(30))(a-h,o-z),integer(i-n) 1015 integer, parameter :: rk=selected_real_kind(30) 1016#else 1017 implicit real(kind=selected_real_kind(15))(a-h,o-z),integer(i-n) 1018 integer, parameter :: rk=selected_real_kind(15) 1019#endif 1020#include "nwxc_param.fh" 1021 double precision param(*) !< [Input] Parameters of functional 1022 double precision tol_rho !< [Input] The lower limit on the density 1023 integer ipol !< [Input] The number of spin channels 1024 integer nq !< [Input] The number of points 1025 double precision wght !< [Input] The weight of the functional 1026 double precision rho(nq,*) !< [Input] The density 1027 double precision fnc(nq) !< [Output] The value of the functional 1028c 1029c Sampling Matrices for the XC Kernel 1030c 1031 double precision Amat(nq,*) !< [Output] The derivative wrt rho 1032c 1033c Sampling Matrices for the XC Kernel 1034c 1035 double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 1036c 1037c Sampling Matrices for the XC Kernel 1038c 1039 double precision Amat3(nq,*) !< [Output] The 3rd derivative wrt rho 1040c 1041 integer iq 1042 double precision tmp 1043 double precision rhoa,rhob 1044 double precision gammaaa,gammaab,gammabb 1045 double precision taua,taub 1046 double precision nwxcm_heaviside 1047 external nwxcm_heaviside 1048CDIR$ NOVECTOR 1049 do iq = 1, nq 1050 if (ipol.eq.1) then 1051 rhoa = 0.5d0*rho(iq,R_T) 1052 if (rhoa.gt.tol_rho) then 1053 t1 = rhoa**3.333333333333333d-1 1054 t2 = t1**5.0d-1 1055 t3 = 1/t2 1056 t4 = 1/t1 1057 t5 = 2.1508070719090538d-2*t3+1.0186556948d-1 1058 t6 = 7.016926042943222d-1*t3*t5+2.2308199064d-1 1059 t7 = 7.016926042943222d-1*t3*t6+4.7231125998d-1 1060 t8 = 1/t7 1061 t9 = 1.425125466450768d+0*t2*t8+1.0d+0 1062 t10 = log(t9) 1063 t11 = 1.0522000558389213d-1*t4+1.0d+0 1064 t12 = 1/rhoa**1.3333333333333333d+0 1065 t13 = 1/t9 1066 t14 = 1/t7**2 1067 t15 = 1/rhoa**1.3333333333333336d+0 1068 t16 = 1/t2**3 1069 t17 = 1/rhoa**6.666666666666667d-1 1070 t18 = -5.847438369119352d-2*t16*t17*t5-1.2576711796854242d-3 1071 1 *t15 1072 t19 = 7.016926042943222d-1*t18*t3-5.847438369119352d-2*t16*t 1073 1 17*t6 1074 t20 = 1.1876045553756398d-1*t17*t3*t8-1.425125466450768d+0*t 1075 1 14*t19*t2 1076 t21 = 1.090454542535705d-3*t10*t12-6.21814d-2*t11*t13*t20 1077 t22 = 2.0d+0*t21*wght 1078 t23 = 1/rhoa**2.333333333333333d+0 1079 t24 = -7.269696950238034d-4*t10*t23 1080 t25 = 1.1771443702974158d-2*t3+2.9729725188d-2 1081 t26 = 7.016926042943222d-1*t25*t3+1.2236585478d-1 1082 t27 = 7.016926042943222d-1*t26*t3+3.497952466d-1 1083 t28 = 1/t27 1084 t29 = 1.425125466450768d+0*t2*t28+1.0d+0 1085 t30 = log(t29) 1086 t31 = 5.477644184000001d-2*t4+1.0d+0 1087 t32 = 1/rhoa**2 1088 t33 = 2.18090908507141d-3*t12*t13*t20 1089 t34 = 1/t9**2 1090 t35 = t20**2 1091 t36 = 6.21814d-2*t11*t34*t35 1092 t37 = 1/t7**3 1093 t38 = t19**2 1094 t39 = 1/rhoa**2.3333333333333334d+0 1095 t40 = 1/rhoa**1.6666666666666669d+0 1096 t41 = 1/t2**5 1097 t42 = 1.4618595922798375d-2*t15*t41*t5+1.949146123039784d-2* 1098 1 t16*t40*t5+9.432533847640683d-4*t39 1099 t43 = 1.4618595922798375d-2*t15*t41*t6+1.949146123039784d-2* 1100 1 t16*t40*t6+7.016926042943222d-1*t3*t42-1.1694876738238703 1101 2 d-1*t16*t17*t18 1102 t44 = -3.9586818512521327d-2*t3*t40*t8-9.896704628130328d-3* 1103 1 t15*t16*t8-1.425125466450768d+0*t14*t2*t43+2.850250932901 1104 2 536d+0*t2*t37*t38-2.3752091107512796d-1*t14*t17*t19*t3 1105 t45 = -6.21814d-2*t11*t13*t44 1106 t46 = t45+t36+t33+8.443450000000001d-3*t30*t31*t32+t24 1107 t47 = t45+t36+t33-8.44345d-3*t30*t31*t32+t24 1108 t48 = 8.48131310861104d-4*t10/rhoa**3.333333333333333d+0 1109 t49 = 1/rhoa**3 1110 t50 = -2.1809090850714105d-3*t13*t20*t23 1111 t51 = 3.37738d-2*(1.1876045553756398d-1*t17*t28*t3-1.4251254 1112 1 66450768d+0*t2*(7.016926042943222d-1*(-5.847438369119352d 1113 2 -2*t16*t17*t25-6.883279156869946d-4*t15)*t3-5.84743836911 1114 3 9352d-2*t16*t17*t26)/t27**2)*t31/t29-3.083347652359653d-4 1115 4 *t12*t30 1116 t52 = -3.2713636276071156d-3*t12*t34*t35 1117 t53 = 3.2713636276071156d-3*t12*t13*t44 1118 t54 = -1.243628d-1*t11*t20**3/t9**3 1119 t55 = 1.865442d-1*t11*t20*t34*t44 1120 t56 = 1/rhoa**2.666666666666667d+0 1121 t57 = 1/t2**7 1122 t58 = 1/rhoa**2.0d+0 1123 t59 = -6.21814d-2*t11*t13*(2.474176157032582d-3*t41*t58*t8+3 1124 1 .29890154271011d-2*t3*t56*t8+9.89670462813033d-3*t16*t39* 1125 2 t8-8.550752798704606d+0*t19**3*t2/t7**4-1.425125466450768 1126 3 d+0*t14*t2*(-6.091081634499322d-3*t57*t58*t6-1.6242884358 1127 4 664864d-2*t16*t56*t6-1.4618595922798375d-2*t39*t41*t6+7.0 1128 5 16926042943222d-1*t3*(-6.091081634499322d-3*t5*t57*t58-1. 1129 6 6242884358664864d-2*t16*t5*t56-1.4618595922798375d-2*t39* 1130 7 t41*t5-1.161599075681677d-3/rhoa**3.3333333333333337d+0)- 1131 8 1.7542315107358056d-1*t16*t17*t42+4.385578776839513d-2*t1 1132 9 5*t18*t41+5.847438369119352d-2*t16*t18*t40)+8.55075279870 1133 : 4606d+0*t19*t2*t37*t43-3.5628136661269194d-1*t14*t17*t3*t 1134 ; 43+1.1876045553756398d-1*t14*t19*t3*t40+7.125627332253839 1135 < d-1*t17*t3*t37*t38+2.969011388439098d-2*t14*t15*t16*t19) 1136 fnc(iq) = fnc(iq)-1.243628d-1*rhoa*log(1.4251254664507676d+0 1137 1 *t2/(7.016926042943223d-1*t3*(7.016926042943223d-1*(2.150 1138 2 807071909054d-2*t3+1.0186556948d-1)*t3+2.2308199064d-1)+4 1139 3 .7231125998d-1)+1.0d+0)*(1.0522000558389215d-1*t4+1.0d+0) 1140 4 *wght 1141 Amat(iq,D1_RA) = 2.0d+0*rhoa*t21*wght-6.21814d-2*t10*t11*wgh 1142 1 t+Amat(iq,D1_RA) 1143 Amat2(iq,D2_RA_RA) = 2.0d+0*rhoa*t46*wght+t22+Amat2(iq,D2_RA 1144 1 _RA) 1145 Amat2(iq,D2_RA_RB) = 2.0d+0*rhoa*t47*wght+t22+Amat2(iq,D2_RA 1146 1 _RB) 1147 Amat3(iq,D3_RA_RA_RA) = 2.0d+0*rhoa*(t59+t55+t54+t53+t52+7.5 1148 1 00000000000002d-1*t32*t51+t50-2.533035d-2*t30*t31*t49+t48 1149 2 )*wght+3.0d+0*t46*wght+Amat3(iq,D3_RA_RA_RA) 1150 Amat3(iq,D3_RA_RA_RB) = 2.0d+0*rhoa*(t59+t55+t54+t53+t52-2.4 1151 1 999999999999994d-1*t32*t51+t50+8.44345d-3*t30*t31*t49+t48 1152 2 )*wght+2.0d+0*t47*wght+1.0d+0*t46*wght+Amat3(iq,D3_RA_RA_ 1153 3 RB) 1154 endif ! rhoa.gt.tol_rho 1155 else ! ipol.eq.1 1156 rhoa = rho(iq,R_A) 1157 rhob = rho(iq,R_B) 1158 if (rhoa.gt.tol_rho.and.rhob.gt.tol_rho) then 1159 t1 = rhob+rhoa 1160 t2 = t1**3.333333333333333d-1 1161 t3 = 1/t2 1162 t4 = 1.325688999052018d-1*t3+1.0d+0 1163 t5 = t2**5.0d-1 1164 t6 = 1/t5 1165 t7 = 2.4141993114533214d-2*t6+1.0186556948d-1 1166 t8 = 7.876233178997433d-1*t6*t7+2.2308199064d-1 1167 t9 = 7.876233178997433d-1*t6*t8+4.7231125998d-1 1168 t10 = 1/t9 1169 t11 = 1.269642451250142d+0*t10*t5+1.0d+0 1170 t12 = log(t11) 1171 t13 = rhoa-rhob 1172 t14 = 1/t1 1173 t15 = 1.0d+0-t13*t14 1174 t16 = t13*t14+1.0d+0 1175 t17 = 1.923661050931536d+0*t16**1.3333333333333333d+0+1.9236 1176 1 61050931536d+0*t15**1.3333333333333333d+0-3.8473221018630 1177 2 72d+0 1178 t18 = 6.901399211255826d-2*t3+1.0d+0 1179 t19 = 1.3212998810398843d-2*t6+2.9729725188d-2 1180 t20 = 7.876233178997433d-1*t19*t6+1.2236585478d-1 1181 t21 = 7.876233178997433d-1*t20*t6+3.497952466d-1 1182 t22 = 1/t21 1183 t23 = 1.269642451250142d+0*t22*t5+1.0d+0 1184 t24 = log(t23) 1185 t25 = t13**4 1186 t26 = 1/t1**4 1187 t27 = 1.2746961887000874d-1*t3+1.0d+0 1188 t28 = 1.530901310039024d-2*t6+1.0465751434d-1 1189 t29 = 7.876233178997433d-1*t28*t6+1.9269083139d-1 1190 t30 = 7.876233178997433d-1*t29*t6+4.3896648423d-1 1191 t31 = 1/t30 1192 t32 = 1.269642451250142d+0*t31*t5+1.0d+0 1193 t33 = log(t32) 1194 t34 = 1.709920934161365d+0*(6.21814d-2*t12*t4-3.10907d-2*t27 1195 1 *t33)-3.37738d-2*t18*t24 1196 t35 = t25*t26*t34+3.37738d-2*t18*t24 1197 t36 = 5.848223622634648d-1*t17*t35-6.21814d-2*t12*t4 1198 t37 = 1/t1**1.3333333333333336d+0 1199 t38 = 1/t1**6.666666666666667d-1 1200 t39 = 1/t5**3 1201 t40 = -1.3127055298329054d-1*t38*t39*t7-3.169132786263567d-3 1202 1 *t37 1203 t41 = 7.876233178997433d-1*t40*t6-1.3127055298329054d-1*t38* 1204 1 t39*t8 1205 t42 = 1/t9**2 1206 t43 = 2.11607075208357d-1*t10*t38*t6-1.269642451250142d+0*t4 1207 1 1*t42*t5 1208 t44 = 1/t11 1209 t45 = -6.21814d-2*t4*t43*t44 1210 t46 = 1/t1**1.3333333333333333d+0 1211 t47 = 2.747773264188438d-3*t12*t46 1212 t48 = -1.3127055298329054d-1*t19*t38*t39-1.7344776604086162d 1213 1 -3*t37 1214 t49 = 7.876233178997433d-1*t48*t6-1.3127055298329054d-1*t20* 1215 1 t38*t39 1216 t50 = 1/t21**2 1217 t51 = 2.11607075208357d-1*t22*t38*t6-1.269642451250142d+0*t4 1218 1 9*t5*t50 1219 t52 = 1/t23 1220 t53 = 3.37738d-2*t18*t51*t52 1221 t54 = -7.769549222703733d-4*t24*t46 1222 t55 = -1.3127055298329054d-1*t28*t38*t39-2.0096226153166658d 1223 1 -3*t37 1224 t56 = 7.876233178997433d-1*t55*t6-1.3127055298329054d-1*t29* 1225 1 t38*t39 1226 t57 = 1/t30**2 1227 t58 = 2.11607075208357d-1*t31*t38*t6-1.269642451250142d+0*t5 1228 1 *t56*t57 1229 t59 = 1/t32 1230 t60 = 1.709920934161365d+0*(-3.10907d-2*t27*t58*t59+1.321039 1231 1 8931339265d-3*t33*t46-2.747773264188438d-3*t12*t46+6.2181 1232 2 4d-2*t4*t43*t44)-3.37738d-2*t18*t51*t52+7.769549222703733 1233 3 d-4*t24*t46 1234 t61 = t25*t26*t60 1235 t62 = 1/t1**5 1236 t63 = -4*t25*t34*t62 1237 t64 = t13**3 1238 t65 = 4*t26*t34*t64+t63+t61+t54+t53 1239 t66 = 1/t1**2 1240 t67 = t13*t66 1241 t68 = -t14 1242 t69 = t68+t67 1243 t70 = t15**3.333333333333333d-1 1244 t71 = -t13*t66 1245 t72 = t71+t14 1246 t73 = t16**3.333333333333333d-1 1247 t74 = 2.564881401242048d+0*t72*t73+2.564881401242048d+0*t69* 1248 1 t70 1249 t75 = 5.848223622634648d-1*t35*t74+5.848223622634648d-1*t17* 1250 1 t65+t47+t45 1251 t76 = 1.0d+0*t36*wght 1252 t77 = -4*t26*t34*t64+t63+t61+t54+t53 1253 t78 = t67+t14 1254 t79 = t71+t68 1255 t80 = 2.564881401242048d+0*t73*t79+2.564881401242048d+0*t70* 1256 1 t78 1257 t81 = 5.848223622634648d-1*t35*t80+5.848223622634648d-1*t17* 1258 1 t77+t47+t45 1259 t82 = t43**2 1260 t83 = 1/t11**2 1261 t84 = 6.21814d-2*t4*t82*t83 1262 t85 = t41**2 1263 t86 = 1/t9**3 1264 t87 = 1/t1**2.3333333333333334d+0 1265 t88 = 1/t5**5 1266 t89 = 1/t1**1.6666666666666669d+0 1267 t90 = 8.751370198886037d-2*t39*t7*t89+6.563527649164527d-2*t 1268 1 37*t7*t88+4.753699179395351d-3*t87 1269 t91 = 7.876233178997433d-1*t6*t90+8.751370198886037d-2*t39*t 1270 1 8*t89+6.563527649164527d-2*t37*t8*t88-2.625411059665811d- 1271 2 1*t38*t39*t40 1272 t92 = -1.269642451250142d+0*t42*t5*t91-1.4107138347223802d-1 1273 1 *t10*t6*t89+2.539284902500284d+0*t5*t85*t86-4.23214150416 1274 2 714d-1*t38*t41*t42*t6-3.52678458680595d-2*t10*t37*t39 1275 t93 = -6.21814d-2*t4*t44*t92 1276 t94 = 5.495546528376876d-3*t43*t44*t46 1277 t95 = 1/t1**2.333333333333333d+0 1278 t96 = -3.663697685584584d-3*t12*t95 1279 t97 = t51**2 1280 t98 = 1/t23**2 1281 t99 = -3.37738d-2*t18*t97*t98 1282 t100 = t49**2 1283 t101 = 1/t21**3 1284 t102 = 8.751370198886037d-2*t19*t39*t89+6.563527649164527d-2 1285 1 *t19*t37*t88+2.601716490612924d-3*t87 1286 t103 = 8.751370198886037d-2*t20*t39*t89+6.563527649164527d-2 1287 1 *t20*t37*t88+7.876233178997433d-1*t102*t6-2.6254110596658 1288 2 11d-1*t38*t39*t48 1289 t104 = -1.4107138347223802d-1*t22*t6*t89-4.23214150416714d-1 1290 1 *t38*t49*t50*t6-1.269642451250142d+0*t103*t5*t50+2.539284 1291 2 902500284d+0*t100*t101*t5-3.52678458680595d-2*t22*t37*t39 1292 t105 = 3.37738d-2*t104*t18*t52 1293 t106 = -1.5539098445407465d-3*t46*t51*t52 1294 t107 = 1.0359398963604977d-3*t24*t95 1295 t108 = t58**2 1296 t109 = 1/t32**2 1297 t110 = t56**2 1298 t111 = 1/t30**3 1299 t112 = 8.751370198886037d-2*t28*t39*t89+6.563527649164527d-2 1300 1 *t28*t37*t88+3.0144339229749983d-3*t87 1301 t113 = 8.751370198886037d-2*t29*t39*t89+6.563527649164527d-2 1302 1 *t29*t37*t88+7.876233178997433d-1*t112*t6-2.6254110596658 1303 2 11d-1*t38*t39*t55 1304 t114 = -1.4107138347223802d-1*t31*t6*t89-4.23214150416714d-1 1305 1 *t38*t56*t57*t6-1.269642451250142d+0*t113*t5*t57+2.539284 1306 2 902500284d+0*t110*t111*t5-3.52678458680595d-2*t31*t37*t39 1307 t115 = 3.37738d-2*t18*t97*t98+1.709920934161365d+0*(-1.76138 1308 1 65241785687d-3*t33*t95+3.663697685584584d-3*t12*t95+6.218 1309 2 14d-2*t4*t44*t92-6.21814d-2*t4*t82*t83+2.642079786267853d 1310 3 -3*t46*t58*t59-3.10907d-2*t114*t27*t59-5.495546528376876d 1311 4 -3*t43*t44*t46+3.10907d-2*t108*t109*t27)-1.03593989636049 1312 5 77d-3*t24*t95+1.5539098445407465d-3*t46*t51*t52-3.37738d- 1313 6 2*t104*t18*t52 1314 t116 = t115*t25*t26 1315 t117 = -8*t25*t60*t62 1316 t118 = 1/t1**6 1317 t119 = 20*t118*t25*t34 1318 t120 = t13**2 1319 t121 = 12*t120*t26*t34 1320 t122 = t99-32*t34*t62*t64+8*t26*t60*t64+t121+t119+t117+t116+ 1321 1 t107+t106+t105 1322 t123 = t69**2 1323 t124 = 1/t15**6.666666666666666d-1 1324 t125 = 1/t1**3 1325 t126 = -2*t125*t13 1326 t127 = 2*t66 1327 t128 = t127+t126 1328 t129 = t72**2 1329 t130 = 1/t16**6.666666666666666d-1 1330 t131 = 2*t125*t13 1331 t132 = -2*t66 1332 t133 = t132+t131 1333 t134 = 2.564881401242048d+0*t133*t73+2.564881401242048d+0*t1 1334 1 28*t70+8.549604670806825d-1*t129*t130+8.549604670806825d- 1335 2 1*t123*t124 1336 t135 = t96+t94+t93+t84+1.1696447245269297d+0*t65*t74+5.84822 1337 1 3622634648d-1*t134*t35+5.848223622634648d-1*t122*t17 1338 t136 = t99-12*t120*t26*t34+t119+t117+t116+t107+t106+t105 1339 t137 = 8.549604670806825d-1*t130*t72*t79+8.549604670806825d- 1340 1 1*t124*t69*t78+5.129762802484096d+0*t125*t13*t73-5.129762 1341 2 802484096d+0*t125*t13*t70 1342 t138 = t96+t94+t93+t84+5.848223622634648d-1*t65*t80+5.848223 1343 1 622634648d-1*t74*t77+5.848223622634648d-1*t137*t35+5.8482 1344 2 23622634648d-1*t136*t17 1345 t139 = t99+32*t34*t62*t64-8*t26*t60*t64+t121+t119+t117+t116+ 1346 1 t107+t106+t105 1347 t140 = t78**2 1348 t141 = t132+t126 1349 t142 = t79**2 1350 t143 = t131+t127 1351 t144 = 2.564881401242048d+0*t143*t73+2.564881401242048d+0*t1 1352 1 41*t70+8.549604670806825d-1*t130*t142+8.549604670806825d- 1353 2 1*t124*t140 1354 t145 = t96+t94+t93+t84+1.1696447245269297d+0*t77*t80+5.84822 1355 1 3622634648d-1*t144*t35+5.848223622634648d-1*t139*t17 1356 t146 = t43**3 1357 t147 = 1/t11**3 1358 t148 = -1.243628d-1*t146*t147*t4 1359 t149 = 1.865442d-1*t4*t43*t83*t92 1360 t150 = -8.243319792565315d-3*t46*t82*t83 1361 t151 = 1/t1**3.3333333333333337d+0 1362 t152 = 1/t1**2.0d+0 1363 t153 = 1/t5**7 1364 t154 = 1/t1**2.666666666666667d+0 1365 t155 = 7.617854707500852d+0*t41*t5*t86*t91-6.34821225625071d 1366 1 -1*t38*t42*t6*t91-1.269642451250142d+0*t42*t5*(-3.9381165 1367 2 894987163d-1*t38*t39*t90+2.625411059665811d-1*t39*t40*t89 1368 3 +7.876233178997433d-1*t6*(-1.3127055298329054d-1*t7*t87*t 1369 4 88-1.4585616998143394d-1*t154*t39*t7-5.469606374303773d-2 1370 5 *t152*t153*t7-1.1708185015918181d-2*t151)-1.3127055298329 1371 6 054d-1*t8*t87*t88+1.9690582947493582d-1*t37*t40*t88-1.458 1372 7 5616998143394d-1*t154*t39*t8-5.469606374303773d-2*t152*t1 1373 8 53*t8)-7.617854707500852d+0*t41**3*t5/t9**4+4.23214150416 1374 9 71406d-1*t41*t42*t6*t89+1.763392293402975d-2*t10*t152*t88 1375 : +7.053569173611901d-2*t10*t39*t87+1.269642451250142d+0*t3 1376 ; 8*t6*t85*t86+2.3511897245373004d-1*t10*t154*t6+1.05803537 1377 < 60417849d-1*t37*t39*t41*t42 1378 t156 = -6.21814d-2*t155*t4*t44 1379 t157 = 8.243319792565315d-3*t44*t46*t92 1380 t158 = -1.0991093056753751d-2*t43*t44*t95 1381 t159 = 1/t1**3.333333333333333d+0 1382 t160 = 8.548627933030694d-3*t12*t159 1383 t161 = 1/t15**1.6666666666666669d+0 1384 t162 = 6*t13*t26 1385 t163 = -6*t125 1386 t164 = 1/t16**1.6666666666666669d+0 1387 t165 = -6*t13*t26 1388 t166 = 6*t125 1389 t167 = t51**3 1390 t168 = 1/t23**3 1391 t169 = 6.75476d-2*t167*t168*t18 1392 t170 = -1.013214d-1*t104*t18*t51*t98 1393 t171 = 2.33086476681112d-3*t46*t97*t98 1394 t172 = -1.269642451250142d+0*t5*t50*(2.625411059665811d-1*t3 1395 1 9*t48*t89+7.876233178997433d-1*t6*(-1.3127055298329054d-1 1396 2 *t19*t87*t88-1.4585616998143394d-1*t154*t19*t39-5.4696063 1397 3 74303773d-2*t152*t153*t19-6.407931356509611d-3*t151)-1.31 1398 4 27055298329054d-1*t20*t87*t88+1.9690582947493582d-1*t37*t 1399 5 48*t88-3.9381165894987163d-1*t102*t38*t39-1.4585616998143 1400 6 394d-1*t154*t20*t39-5.469606374303773d-2*t152*t153*t20)+4 1401 7 .2321415041671406d-1*t49*t50*t6*t89+1.763392293402975d-2* 1402 8 t152*t22*t88+7.053569173611901d-2*t22*t39*t87-6.348212256 1403 9 25071d-1*t103*t38*t50*t6+1.269642451250142d+0*t100*t101*t 1404 : 38*t6+2.3511897245373004d-1*t154*t22*t6+1.058035376041784 1405 ; 9d-1*t37*t39*t49*t50-7.617854707500852d+0*t49**3*t5/t21** 1406 < 4+7.617854707500852d+0*t101*t103*t49*t5 1407 t173 = 3.37738d-2*t172*t18*t52 1408 t174 = -2.33086476681112d-3*t104*t46*t52 1409 t175 = 3.107819689081493d-3*t51*t52*t95 1410 t176 = -2.4171930915078277d-3*t159*t24 1411 t177 = t25*t26*(-2.33086476681112d-3*t46*t97*t98+1.013214d-1 1412 1 *t104*t18*t51*t98+1.709920934161365d+0*(-5.28415957253570 1413 2 6d-3*t58*t59*t95+1.0991093056753751d-2*t43*t44*t95-1.8654 1414 3 42d-1*t4*t43*t83*t92-8.243319792565315d-3*t44*t46*t92-3.1 1415 4 0907d-2*t27*t59*(-1.269642451250142d+0*t5*t57*(2.62541105 1416 5 9665811d-1*t39*t55*t89+7.876233178997433d-1*t6*(-1.312705 1417 6 5298329054d-1*t28*t87*t88-1.4585616998143394d-1*t154*t28* 1418 7 t39-5.469606374303773d-2*t152*t153*t28-7.424439106586571d 1419 8 -3*t151)-1.3127055298329054d-1*t29*t87*t88+1.969058294749 1420 9 3582d-1*t37*t55*t88-3.9381165894987163d-1*t112*t38*t39-1. 1421 : 4585616998143394d-1*t154*t29*t39-5.469606374303773d-2*t15 1422 ; 2*t153*t29)+4.2321415041671406d-1*t56*t57*t6*t89+1.763392 1423 < 293402975d-2*t152*t31*t88+7.053569173611901d-2*t31*t39*t8 1424 = 7-6.34821225625071d-1*t113*t38*t57*t6+1.269642451250142d+ 1425 > 0*t110*t111*t38*t6+2.3511897245373004d-1*t154*t31*t6+1.05 1426 ? 80353760417849d-1*t37*t39*t56*t57-7.617854707500852d+0*t5 1427 @ *t56**3/t30**4+7.617854707500852d+0*t111*t113*t5*t56)+8.2 1428 1 43319792565315d-3*t46*t82*t83+3.96311967940178d-3*t114*t4 1429 2 6*t59-6.21814d-2*t27*t58**3/t32**3+9.327209999999999d-2*t 1430 3 109*t114*t27*t58-3.96311967940178d-3*t108*t109*t46+6.2181 1431 4 4d-2*t155*t4*t44+1.243628d-1*t146*t147*t4+4.1099018897499 1432 5 934d-3*t159*t33-8.548627933030694d-3*t12*t159)-3.10781968 1433 6 9081493d-3*t51*t52*t95+2.33086476681112d-3*t104*t46*t52-3 1434 7 .37738d-2*t172*t18*t52+2.4171930915078277d-3*t159*t24-6.7 1435 8 5476d-2*t167*t168*t18) 1436 t178 = -12*t115*t25*t62 1437 t179 = 60*t118*t25*t60 1438 t180 = 36*t120*t26*t60 1439 t181 = -120*t25*t34/t1**7 1440 t182 = -144*t120*t34*t62 1441 t183 = 24*t13*t26*t34 1442 t184 = 2.0d+0*t138*wght 1443 t185 = -12*t120*t26*t60 1444 t186 = 48*t120*t34*t62 1445 t187 = -24*t13*t26*t34 1446 fnc(iq) = 1.0d+0*t1*t36*wght+fnc(iq) 1447 Amat(iq,D1_RA) = 1.0d+0*t1*t75*wght+t76+Amat(iq,D1_RA) 1448 Amat(iq,D1_RB) = 1.0d+0*t1*t81*wght+t76+Amat(iq,D1_RB) 1449 Amat2(iq,D2_RA_RA) = 2.0d+0*t75*wght+1.0d+0*t1*t135*wght+Ama 1450 1 t2(iq,D2_RA_RA) 1451 Amat2(iq,D2_RA_RB) = 1.0d+0*t81*wght+1.0d+0*t75*wght+1.0d+0* 1452 1 t1*t138*wght+Amat2(iq,D2_RA_RB) 1453 Amat2(iq,D2_RB_RB) = 2.0d+0*t81*wght+1.0d+0*t1*t145*wght+Ama 1454 1 t2(iq,D2_RB_RB) 1455 Amat3(iq,D3_RA_RA_RA) = 1.0d+0*t1*(1.7544670867903944d+0*t12 1456 1 2*t74+5.848223622634648d-1*t35*(2.564881401242048d+0*(t16 1457 2 6+t165)*t73-5.69973644720455d-1*t164*t72**3+2.56488140124 1458 3 20473d+0*t130*t133*t72+2.564881401242048d+0*(t163+t162)*t 1459 4 70-5.69973644720455d-1*t161*t69**3+2.5648814012420473d+0* 1460 5 t124*t128*t69)+1.7544670867903944d+0*t134*t65+5.848223622 1461 6 634648d-1*t17*(-96*t60*t62*t64+240*t118*t34*t64+12*t115*t 1462 7 26*t64+t183+t182+t181+t180+t179+t178+t177+t176+t175+t174+ 1463 8 t173+t171+t170+t169)+t160+t158+t157+t156+t150+t149+t148)* 1464 9 wght+3.0d+0*t135*wght+Amat3(iq,D3_RA_RA_RA) 1465 Amat3(iq,D3_RA_RA_RB) = 1.0d+0*t1*(5.848223622634648d-1*t122 1466 1 *t80+5.848223622634648d-1*t35*(-5.69973644720455d-1*t129* 1467 2 t164*t79+8.549604670806825d-1*t130*t133*t79-5.69973644720 1468 3 455d-1*t123*t161*t78+8.549604670806825d-1*t124*t128*t78+2 1469 4 .564881401242048d+0*(t165+2*t125)*t73+3.41984186832273d+0 1470 5 *t125*t13*t130*t72+2.564881401242048d+0*(t162-2*t125)*t70 1471 6 -3.41984186832273d+0*t124*t125*t13*t69)+5.848223622634648 1472 7 d-1*t134*t77+1.1696447245269297d+0*t136*t74+1.16964472452 1473 8 69297d+0*t137*t65+5.848223622634648d-1*t17*(-32*t60*t62*t 1474 9 64+80*t118*t34*t64+4*t115*t26*t64+t187+t186+t185+t181+t17 1475 : 9+t178+t177+t176+t175+t174+t173+t171+t170+t169)+t160+t158 1476 ; +t157+t156+t150+t149+t148)*wght+1.0d+0*t135*wght+t184+Ama 1477 < t3(iq,D3_RA_RA_RB) 1478 Amat3(iq,D3_RA_RB_RB) = 1.0d+0*t1*(1.1696447245269297d+0*t13 1479 1 6*t80+5.848223622634648d-1*t35*(3.41984186832273d+0*t125* 1480 2 t13*t130*t79-3.41984186832273d+0*t124*t125*t13*t78-1.5389 1481 3 288407452287d+1*t13*t26*t73-5.129762802484096d+0*t125*t73 1482 4 -5.69973644720455d-1*t142*t164*t72+8.549604670806825d-1*t 1483 5 130*t143*t72+1.5389288407452287d+1*t13*t26*t70+5.12976280 1484 6 2484096d+0*t125*t70-5.69973644720455d-1*t140*t161*t69+8.5 1485 7 49604670806825d-1*t124*t141*t69)+1.1696447245269297d+0*t1 1486 8 37*t77+5.848223622634648d-1*t139*t74+5.848223622634648d-1 1487 9 *t144*t65+5.848223622634648d-1*t17*(32*t60*t62*t64-80*t11 1488 : 8*t34*t64-4*t115*t26*t64+t186+t185+t183+t181+t179+t178+t1 1489 ; 77+t176+t175+t174+t173+t171+t170+t169)+t160+t158+t157+t15 1490 < 6+t150+t149+t148)*wght+1.0d+0*t145*wght+t184+Amat3(iq,D3_ 1491 = RA_RB_RB) 1492 Amat3(iq,D3_RB_RB_RB) = 1.0d+0*t1*(1.7544670867903944d+0*t13 1493 1 9*t80+5.848223622634648d-1*t35*(-5.69973644720455d-1*t164 1494 2 *t79**3+2.5648814012420473d+0*t130*t143*t79-5.69973644720 1495 3 455d-1*t161*t78**3+2.5648814012420473d+0*t124*t141*t78+2. 1496 4 564881401242048d+0*(t165+t163)*t73+2.564881401242048d+0*( 1497 5 t166+t162)*t70)+1.7544670867903944d+0*t144*t77+5.84822362 1498 6 2634648d-1*t17*(96*t60*t62*t64-240*t118*t34*t64-12*t115*t 1499 7 26*t64+t187+t182+t181+t180+t179+t178+t177+t176+t175+t174+ 1500 8 t173+t171+t170+t169)+t160+t158+t157+t156+t150+t149+t148)* 1501 9 wght+3.0d+0*t145*wght+Amat3(iq,D3_RB_RB_RB) 1502 elseif (rhoa.gt.tol_rho.and.rhob.le.tol_rho) then 1503 t1 = rhoa**3.333333333333333d-1 1504 t2 = t1**5.0d-1 1505 t3 = 1/t2 1506 t4 = 1.530901310039024d-2*t3+1.0465751434d-1 1507 t5 = 7.876233178997433d-1*t3*t4+1.9269083139d-1 1508 t6 = 7.876233178997433d-1*t3*t5+4.3896648423d-1 1509 t7 = 1/t6 1510 t8 = 1.269642451250142d+0*t2*t7+1.0d+0 1511 t9 = log(t8) 1512 t10 = 1/t1 1513 t11 = 1.2746961887000874d-1*t10+1.0d+0 1514 t12 = 2.4141993114533214d-2*t3+1.0186556948d-1 1515 t13 = 7.876233178997433d-1*t12*t3+2.2308199064d-1 1516 t14 = 7.876233178997433d-1*t13*t3+4.7231125998d-1 1517 t15 = 1/t14 1518 t16 = 1.269642451250142d+0*t15*t2+1.0d+0 1519 t17 = log(t16) 1520 t18 = 1.325688999052018d-1*t10+1.0d+0 1521 t19 = 1.3212998810398843d-2*t3+2.9729725188d-2 1522 t20 = 7.876233178997433d-1*t19*t3+1.2236585478d-1 1523 t21 = 7.876233178997433d-1*t20*t3+3.497952466d-1 1524 t22 = 1/t21 1525 t23 = 1.269642451250142d+0*t2*t22+1.0d+0 1526 t24 = log(t23) 1527 t25 = 6.901399211255826d-2*t10+1.0d+0 1528 t26 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(6. 1529 1 21814d-2*t17*t18-3.10907d-2*t11*t9)-3.37738d-2*t24*t25)+3 1530 2 .37738d-2*t24*t25)-6.21814d-2*t17*t18 1531 t27 = 1/rhoa**1.3333333333333333d+0 1532 t28 = 1/t8 1533 t29 = 1/t6**2 1534 t30 = 1/rhoa**1.3333333333333336d+0 1535 t31 = 1/t2**3 1536 t32 = 1/rhoa**6.666666666666667d-1 1537 t33 = -1.3127055298329054d-1*t31*t32*t4-2.0096226153166658d- 1538 1 3*t30 1539 t34 = 7.876233178997433d-1*t3*t33-1.3127055298329054d-1*t31* 1540 1 t32*t5 1541 t35 = 2.11607075208357d-1*t3*t32*t7-1.269642451250142d+0*t2* 1542 1 t29*t34 1543 t36 = 1/t16 1544 t37 = 1/t14**2 1545 t38 = -1.3127055298329054d-1*t12*t31*t32-3.169132786263567d- 1546 1 3*t30 1547 t39 = 7.876233178997433d-1*t3*t38-1.3127055298329054d-1*t13* 1548 1 t31*t32 1549 t40 = 2.11607075208357d-1*t15*t3*t32-1.269642451250142d+0*t2 1550 1 *t37*t39 1551 t41 = 1/t23 1552 t42 = 1/t21**2 1553 t43 = -1.3127055298329054d-1*t19*t31*t32-1.7344776604086162d 1554 1 -3*t30 1555 t44 = 7.876233178997433d-1*t3*t43-1.3127055298329054d-1*t20* 1556 1 t31*t32 1557 t45 = 2.11607075208357d-1*t22*t3*t32-1.269642451250142d+0*t2 1558 1 *t42*t44 1559 t46 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(1. 1560 1 3210398931339265d-3*t27*t9+6.21814d-2*t18*t36*t40-3.10907 1561 2 d-2*t11*t28*t35-2.747773264188438d-3*t17*t27)-3.37738d-2* 1562 3 t25*t41*t45+7.769549222703733d-4*t24*t27)+3.37738d-2*t25* 1563 4 t41*t45-7.769549222703733d-4*t24*t27)-6.21814d-2*t18*t36* 1564 5 t40+2.747773264188438d-3*t17*t27 1565 t47 = 1/rhoa**2.333333333333333d+0 1566 t48 = 1/t8**2 1567 t49 = t35**2 1568 t50 = 1/t6**3 1569 t51 = t34**2 1570 t52 = 1/rhoa**2.3333333333333334d+0 1571 t53 = 1/rhoa**1.6666666666666669d+0 1572 t54 = 1/t2**5 1573 t55 = 6.563527649164527d-2*t30*t4*t54+8.751370198886037d-2*t 1574 1 31*t4*t53+3.0144339229749983d-3*t52 1575 t56 = 7.876233178997433d-1*t3*t55+6.563527649164527d-2*t30*t 1576 1 5*t54+8.751370198886037d-2*t31*t5*t53-2.625411059665811d- 1577 2 1*t31*t32*t33 1578 t57 = -1.4107138347223802d-1*t3*t53*t7-3.52678458680595d-2*t 1579 1 30*t31*t7-1.269642451250142d+0*t2*t29*t56+2.5392849025002 1580 2 84d+0*t2*t50*t51-4.23214150416714d-1*t29*t3*t32*t34 1581 t58 = 1/t16**2 1582 t59 = t40**2 1583 t60 = 1/t14**3 1584 t61 = t39**2 1585 t62 = 6.563527649164527d-2*t12*t30*t54+8.751370198886037d-2* 1586 1 t12*t31*t53+4.753699179395351d-3*t52 1587 t63 = 7.876233178997433d-1*t3*t62+6.563527649164527d-2*t13*t 1588 1 30*t54+8.751370198886037d-2*t13*t31*t53-2.625411059665811 1589 2 d-1*t31*t32*t38 1590 t64 = -1.269642451250142d+0*t2*t37*t63+2.539284902500284d+0* 1591 1 t2*t60*t61-1.4107138347223802d-1*t15*t3*t53-4.23214150416 1592 2 714d-1*t3*t32*t37*t39-3.52678458680595d-2*t15*t30*t31 1593 t65 = 1/t23**2 1594 t66 = t45**2 1595 t67 = 1/t21**3 1596 t68 = t44**2 1597 t69 = 6.563527649164527d-2*t19*t30*t54+8.751370198886037d-2* 1598 1 t19*t31*t53+2.601716490612924d-3*t52 1599 t70 = 7.876233178997433d-1*t3*t69+6.563527649164527d-2*t20*t 1600 1 30*t54+8.751370198886037d-2*t20*t31*t53-2.625411059665811 1601 2 d-1*t31*t32*t43 1602 t71 = -1.269642451250142d+0*t2*t42*t70+2.539284902500284d+0* 1603 1 t2*t67*t68-1.4107138347223802d-1*t22*t3*t53-4.23214150416 1604 2 714d-1*t3*t32*t42*t44-3.52678458680595d-2*t22*t30*t31 1605 t72 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(-1 1606 1 .7613865241785687d-3*t47*t9+6.21814d-2*t18*t36*t64-6.2181 1607 2 4d-2*t18*t58*t59-3.10907d-2*t11*t28*t57+3.10907d-2*t11*t4 1608 3 8*t49+3.663697685584584d-3*t17*t47-5.495546528376876d-3*t 1609 4 27*t36*t40+2.642079786267853d-3*t27*t28*t35)-3.37738d-2*t 1610 5 25*t41*t71+3.37738d-2*t25*t65*t66-1.0359398963604977d-3*t 1611 6 24*t47+1.5539098445407465d-3*t27*t41*t45)+3.37738d-2*t25* 1612 7 t41*t71-3.37738d-2*t25*t65*t66+1.0359398963604977d-3*t24* 1613 8 t47-1.5539098445407465d-3*t27*t41*t45)-6.21814d-2*t18*t36 1614 9 *t64+6.21814d-2*t18*t58*t59-3.663697685584584d-3*t17*t47+ 1615 : 5.495546528376876d-3*t27*t36*t40 1616 t73 = 1/rhoa**3.333333333333333d+0 1617 t74 = 1/rhoa**3.3333333333333337d+0 1618 t75 = 1/rhoa**2.666666666666667d+0 1619 t76 = 1/t2**7 1620 t77 = 1/rhoa**2.0d+0 1621 t78 = 1/t16**3 1622 t79 = t40**3 1623 t80 = -1.269642451250142d+0*t2*t37*(7.876233178997433d-1*t3* 1624 1 (-5.469606374303773d-2*t12*t76*t77-1.4585616998143394d-1* 1625 2 t12*t31*t75-1.1708185015918181d-2*t74-1.3127055298329054d 1626 3 -1*t12*t52*t54)-5.469606374303773d-2*t13*t76*t77-1.458561 1627 4 6998143394d-1*t13*t31*t75-3.9381165894987163d-1*t31*t32*t 1628 5 62-1.3127055298329054d-1*t13*t52*t54+1.9690582947493582d- 1629 6 1*t30*t38*t54+2.625411059665811d-1*t31*t38*t53)+1.7633922 1630 7 93402975d-2*t15*t54*t77+2.3511897245373004d-1*t15*t3*t75+ 1631 8 7.617854707500852d+0*t2*t39*t60*t63-6.34821225625071d-1*t 1632 9 3*t32*t37*t63+1.269642451250142d+0*t3*t32*t60*t61+4.23214 1633 : 15041671406d-1*t3*t37*t39*t53+7.053569173611901d-2*t15*t3 1634 ; 1*t52-7.617854707500852d+0*t2*t39**3/t14**4+1.05803537604 1635 < 17849d-1*t30*t31*t37*t39 1636 t81 = 1/t23**3 1637 t82 = t45**3 1638 t83 = -1.269642451250142d+0*t2*t42*(7.876233178997433d-1*t3* 1639 1 (-5.469606374303773d-2*t19*t76*t77-1.4585616998143394d-1* 1640 2 t19*t31*t75-6.407931356509611d-3*t74-1.3127055298329054d- 1641 3 1*t19*t52*t54)-5.469606374303773d-2*t20*t76*t77-1.4585616 1642 4 998143394d-1*t20*t31*t75-3.9381165894987163d-1*t31*t32*t6 1643 5 9-1.3127055298329054d-1*t20*t52*t54+1.9690582947493582d-1 1644 6 *t30*t43*t54+2.625411059665811d-1*t31*t43*t53)+1.76339229 1645 7 3402975d-2*t22*t54*t77+2.3511897245373004d-1*t22*t3*t75+7 1646 8 .617854707500852d+0*t2*t44*t67*t70-6.34821225625071d-1*t3 1647 9 *t32*t42*t70+1.269642451250142d+0*t3*t32*t67*t68+4.232141 1648 : 5041671406d-1*t3*t42*t44*t53+7.053569173611901d-2*t22*t31 1649 ; *t52-7.617854707500852d+0*t2*t44**3/t21**4+1.058035376041 1650 < 7849d-1*t30*t31*t42*t44 1651 fnc(iq) = 1.0d+0*rhoa*t26*wght+fnc(iq) 1652 Amat(iq,D1_RA) = 1.0d+0*rhoa*t46*wght+1.0d+0*t26*wght+Amat(i 1653 1 q,D1_RA) 1654 Amat2(iq,D2_RA_RA) = 1.0d+0*rhoa*t72*wght+2.0d+0*t46*wght+Am 1655 1 at2(iq,D2_RA_RA) 1656 Amat3(iq,D3_RA_RA_RA) = 1.0d+0*rhoa*(5.848223622634643d-1*(1 1657 1 .0d+0*(1.709920934161365d+0*(4.1099018897499934d-3*t73*t9 1658 2 +6.21814d-2*t18*t36*t80-6.21814d-2*t11*t35**3/t8**3+1.243 1659 3 628d-1*t18*t78*t79-3.10907d-2*t11*t28*(-1.269642451250142 1660 4 d+0*t2*t29*(7.876233178997433d-1*t3*(-5.469606374303773d- 1661 5 2*t4*t76*t77-1.4585616998143394d-1*t31*t4*t75-7.424439106 1662 6 586571d-3*t74-1.3127055298329054d-1*t4*t52*t54)-5.4696063 1663 7 74303773d-2*t5*t76*t77-1.4585616998143394d-1*t31*t5*t75-3 1664 8 .9381165894987163d-1*t31*t32*t55-1.3127055298329054d-1*t5 1665 9 *t52*t54+1.9690582947493582d-1*t30*t33*t54+2.625411059665 1666 : 811d-1*t31*t33*t53)+1.763392293402975d-2*t54*t7*t77+2.351 1667 ; 1897245373004d-1*t3*t7*t75+7.053569173611901d-2*t31*t52*t 1668 < 7-7.617854707500852d+0*t2*t34**3/t6**4+7.617854707500852d 1669 = +0*t2*t34*t50*t56-6.34821225625071d-1*t29*t3*t32*t56+4.23 1670 > 21415041671406d-1*t29*t3*t34*t53+1.269642451250142d+0*t3* 1671 ? t32*t50*t51+1.0580353760417849d-1*t29*t30*t31*t34)-8.5486 1672 @ 27933030694d-3*t17*t73-1.865442d-1*t18*t40*t58*t64-8.2433 1673 1 19792565315d-3*t27*t36*t64+8.243319792565315d-3*t27*t58*t 1674 2 59+9.327209999999999d-2*t11*t35*t48*t57+3.96311967940178d 1675 3 -3*t27*t28*t57-3.96311967940178d-3*t27*t48*t49+1.09910930 1676 4 56753751d-2*t36*t40*t47-5.284159572535706d-3*t28*t35*t47) 1677 5 -3.37738d-2*t25*t41*t83-6.75476d-2*t25*t81*t82+2.41719309 1678 6 15078277d-3*t24*t73+1.013214d-1*t25*t45*t65*t71+2.3308647 1679 7 6681112d-3*t27*t41*t71-2.33086476681112d-3*t27*t65*t66-3. 1680 8 107819689081493d-3*t41*t45*t47)+3.37738d-2*t25*t41*t83+6. 1681 9 75476d-2*t25*t81*t82-2.4171930915078277d-3*t24*t73-1.0132 1682 : 14d-1*t25*t45*t65*t71-2.33086476681112d-3*t27*t41*t71+2.3 1683 ; 3086476681112d-3*t27*t65*t66+3.107819689081493d-3*t41*t45 1684 < *t47)-6.21814d-2*t18*t36*t80-1.243628d-1*t18*t78*t79+8.54 1685 = 8627933030694d-3*t17*t73+1.865442d-1*t18*t40*t58*t64+8.24 1686 > 3319792565315d-3*t27*t36*t64-8.243319792565315d-3*t27*t58 1687 ? *t59-1.0991093056753751d-2*t36*t40*t47)*wght+3.0d+0*t72*w 1688 @ ght+Amat3(iq,D3_RA_RA_RA) 1689 elseif (rhoa.le.tol_rho.and.rhob.gt.tol_rho) then 1690 t1 = rhob**3.333333333333333d-1 1691 t2 = t1**5.0d-1 1692 t3 = 1/t2 1693 t4 = 1.530901310039024d-2*t3+1.0465751434d-1 1694 t5 = 7.876233178997433d-1*t3*t4+1.9269083139d-1 1695 t6 = 7.876233178997433d-1*t3*t5+4.3896648423d-1 1696 t7 = 1/t6 1697 t8 = 1.269642451250142d+0*t2*t7+1.0d+0 1698 t9 = log(t8) 1699 t10 = 1/t1 1700 t11 = 1.2746961887000874d-1*t10+1.0d+0 1701 t12 = 2.4141993114533214d-2*t3+1.0186556948d-1 1702 t13 = 7.876233178997433d-1*t12*t3+2.2308199064d-1 1703 t14 = 7.876233178997433d-1*t13*t3+4.7231125998d-1 1704 t15 = 1/t14 1705 t16 = 1.269642451250142d+0*t15*t2+1.0d+0 1706 t17 = log(t16) 1707 t18 = 1.325688999052018d-1*t10+1.0d+0 1708 t19 = 1.3212998810398843d-2*t3+2.9729725188d-2 1709 t20 = 7.876233178997433d-1*t19*t3+1.2236585478d-1 1710 t21 = 7.876233178997433d-1*t20*t3+3.497952466d-1 1711 t22 = 1/t21 1712 t23 = 1.269642451250142d+0*t2*t22+1.0d+0 1713 t24 = log(t23) 1714 t25 = 6.901399211255826d-2*t10+1.0d+0 1715 t26 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(6. 1716 1 21814d-2*t17*t18-3.10907d-2*t11*t9)-3.37738d-2*t24*t25)+3 1717 2 .37738d-2*t24*t25)-6.21814d-2*t17*t18 1718 t27 = 1/rhob**1.3333333333333333d+0 1719 t28 = 1/t8 1720 t29 = 1/t6**2 1721 t30 = 1/rhob**1.3333333333333336d+0 1722 t31 = 1/t2**3 1723 t32 = 1/rhob**6.666666666666667d-1 1724 t33 = -1.3127055298329054d-1*t31*t32*t4-2.0096226153166658d- 1725 1 3*t30 1726 t34 = 7.876233178997433d-1*t3*t33-1.3127055298329054d-1*t31* 1727 1 t32*t5 1728 t35 = 2.11607075208357d-1*t3*t32*t7-1.269642451250142d+0*t2* 1729 1 t29*t34 1730 t36 = 1/t16 1731 t37 = 1/t14**2 1732 t38 = -1.3127055298329054d-1*t12*t31*t32-3.169132786263567d- 1733 1 3*t30 1734 t39 = 7.876233178997433d-1*t3*t38-1.3127055298329054d-1*t13* 1735 1 t31*t32 1736 t40 = 2.11607075208357d-1*t15*t3*t32-1.269642451250142d+0*t2 1737 1 *t37*t39 1738 t41 = 1/t23 1739 t42 = 1/t21**2 1740 t43 = -1.3127055298329054d-1*t19*t31*t32-1.7344776604086162d 1741 1 -3*t30 1742 t44 = 7.876233178997433d-1*t3*t43-1.3127055298329054d-1*t20* 1743 1 t31*t32 1744 t45 = 2.11607075208357d-1*t22*t3*t32-1.269642451250142d+0*t2 1745 1 *t42*t44 1746 t46 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(1. 1747 1 3210398931339265d-3*t27*t9+6.21814d-2*t18*t36*t40-3.10907 1748 2 d-2*t11*t28*t35-2.747773264188438d-3*t17*t27)-3.37738d-2* 1749 3 t25*t41*t45+7.769549222703733d-4*t24*t27)+3.37738d-2*t25* 1750 4 t41*t45-7.769549222703733d-4*t24*t27)-6.21814d-2*t18*t36* 1751 5 t40+2.747773264188438d-3*t17*t27 1752 t47 = 1/rhob**2.333333333333333d+0 1753 t48 = 1/t8**2 1754 t49 = t35**2 1755 t50 = 1/t6**3 1756 t51 = t34**2 1757 t52 = 1/rhob**2.3333333333333334d+0 1758 t53 = 1/rhob**1.6666666666666669d+0 1759 t54 = 1/t2**5 1760 t55 = 6.563527649164527d-2*t30*t4*t54+8.751370198886037d-2*t 1761 1 31*t4*t53+3.0144339229749983d-3*t52 1762 t56 = 7.876233178997433d-1*t3*t55+6.563527649164527d-2*t30*t 1763 1 5*t54+8.751370198886037d-2*t31*t5*t53-2.625411059665811d- 1764 2 1*t31*t32*t33 1765 t57 = -1.4107138347223802d-1*t3*t53*t7-3.52678458680595d-2*t 1766 1 30*t31*t7-1.269642451250142d+0*t2*t29*t56+2.5392849025002 1767 2 84d+0*t2*t50*t51-4.23214150416714d-1*t29*t3*t32*t34 1768 t58 = 1/t16**2 1769 t59 = t40**2 1770 t60 = 1/t14**3 1771 t61 = t39**2 1772 t62 = 6.563527649164527d-2*t12*t30*t54+8.751370198886037d-2* 1773 1 t12*t31*t53+4.753699179395351d-3*t52 1774 t63 = 7.876233178997433d-1*t3*t62+6.563527649164527d-2*t13*t 1775 1 30*t54+8.751370198886037d-2*t13*t31*t53-2.625411059665811 1776 2 d-1*t31*t32*t38 1777 t64 = -1.269642451250142d+0*t2*t37*t63+2.539284902500284d+0* 1778 1 t2*t60*t61-1.4107138347223802d-1*t15*t3*t53-4.23214150416 1779 2 714d-1*t3*t32*t37*t39-3.52678458680595d-2*t15*t30*t31 1780 t65 = 1/t23**2 1781 t66 = t45**2 1782 t67 = 1/t21**3 1783 t68 = t44**2 1784 t69 = 6.563527649164527d-2*t19*t30*t54+8.751370198886037d-2* 1785 1 t19*t31*t53+2.601716490612924d-3*t52 1786 t70 = 7.876233178997433d-1*t3*t69+6.563527649164527d-2*t20*t 1787 1 30*t54+8.751370198886037d-2*t20*t31*t53-2.625411059665811 1788 2 d-1*t31*t32*t43 1789 t71 = -1.269642451250142d+0*t2*t42*t70+2.539284902500284d+0* 1790 1 t2*t67*t68-1.4107138347223802d-1*t22*t3*t53-4.23214150416 1791 2 714d-1*t3*t32*t42*t44-3.52678458680595d-2*t22*t30*t31 1792 t72 = 5.848223622634643d-1*(1.0d+0*(1.709920934161365d+0*(-1 1793 1 .7613865241785687d-3*t47*t9+6.21814d-2*t18*t36*t64-6.2181 1794 2 4d-2*t18*t58*t59-3.10907d-2*t11*t28*t57+3.10907d-2*t11*t4 1795 3 8*t49+3.663697685584584d-3*t17*t47-5.495546528376876d-3*t 1796 4 27*t36*t40+2.642079786267853d-3*t27*t28*t35)-3.37738d-2*t 1797 5 25*t41*t71+3.37738d-2*t25*t65*t66-1.0359398963604977d-3*t 1798 6 24*t47+1.5539098445407465d-3*t27*t41*t45)+3.37738d-2*t25* 1799 7 t41*t71-3.37738d-2*t25*t65*t66+1.0359398963604977d-3*t24* 1800 8 t47-1.5539098445407465d-3*t27*t41*t45)-6.21814d-2*t18*t36 1801 9 *t64+6.21814d-2*t18*t58*t59-3.663697685584584d-3*t17*t47+ 1802 : 5.495546528376876d-3*t27*t36*t40 1803 t73 = 1/rhob**3.333333333333333d+0 1804 t74 = 1/rhob**3.3333333333333337d+0 1805 t75 = 1/rhob**2.666666666666667d+0 1806 t76 = 1/t2**7 1807 t77 = 1/rhob**2.0d+0 1808 t78 = 1/t16**3 1809 t79 = t40**3 1810 t80 = -1.269642451250142d+0*t2*t37*(7.876233178997433d-1*t3* 1811 1 (-5.469606374303773d-2*t12*t76*t77-1.4585616998143394d-1* 1812 2 t12*t31*t75-1.1708185015918181d-2*t74-1.3127055298329054d 1813 3 -1*t12*t52*t54)-5.469606374303773d-2*t13*t76*t77-1.458561 1814 4 6998143394d-1*t13*t31*t75-3.9381165894987163d-1*t31*t32*t 1815 5 62-1.3127055298329054d-1*t13*t52*t54+1.9690582947493582d- 1816 6 1*t30*t38*t54+2.625411059665811d-1*t31*t38*t53)+1.7633922 1817 7 93402975d-2*t15*t54*t77+2.3511897245373004d-1*t15*t3*t75+ 1818 8 7.617854707500852d+0*t2*t39*t60*t63-6.34821225625071d-1*t 1819 9 3*t32*t37*t63+1.269642451250142d+0*t3*t32*t60*t61+4.23214 1820 : 15041671406d-1*t3*t37*t39*t53+7.053569173611901d-2*t15*t3 1821 ; 1*t52-7.617854707500852d+0*t2*t39**3/t14**4+1.05803537604 1822 < 17849d-1*t30*t31*t37*t39 1823 t81 = 1/t23**3 1824 t82 = t45**3 1825 t83 = -1.269642451250142d+0*t2*t42*(7.876233178997433d-1*t3* 1826 1 (-5.469606374303773d-2*t19*t76*t77-1.4585616998143394d-1* 1827 2 t19*t31*t75-6.407931356509611d-3*t74-1.3127055298329054d- 1828 3 1*t19*t52*t54)-5.469606374303773d-2*t20*t76*t77-1.4585616 1829 4 998143394d-1*t20*t31*t75-3.9381165894987163d-1*t31*t32*t6 1830 5 9-1.3127055298329054d-1*t20*t52*t54+1.9690582947493582d-1 1831 6 *t30*t43*t54+2.625411059665811d-1*t31*t43*t53)+1.76339229 1832 7 3402975d-2*t22*t54*t77+2.3511897245373004d-1*t22*t3*t75+7 1833 8 .617854707500852d+0*t2*t44*t67*t70-6.34821225625071d-1*t3 1834 9 *t32*t42*t70+1.269642451250142d+0*t3*t32*t67*t68+4.232141 1835 : 5041671406d-1*t3*t42*t44*t53+7.053569173611901d-2*t22*t31 1836 ; *t52-7.617854707500852d+0*t2*t44**3/t21**4+1.058035376041 1837 < 7849d-1*t30*t31*t42*t44 1838 fnc(iq) = 1.0d+0*rhob*t26*wght+fnc(iq) 1839 Amat(iq,D1_RB) = 1.0d+0*rhob*t46*wght+1.0d+0*t26*wght+Amat(i 1840 1 q,D1_RB) 1841 Amat2(iq,D2_RB_RB) = 1.0d+0*rhob*t72*wght+2.0d+0*t46*wght+Am 1842 1 at2(iq,D2_RB_RB) 1843 Amat3(iq,D3_RB_RB_RB) = 1.0d+0*rhob*(5.848223622634643d-1*(1 1844 1 .0d+0*(1.709920934161365d+0*(4.1099018897499934d-3*t73*t9 1845 2 +6.21814d-2*t18*t36*t80-6.21814d-2*t11*t35**3/t8**3+1.243 1846 3 628d-1*t18*t78*t79-3.10907d-2*t11*t28*(-1.269642451250142 1847 4 d+0*t2*t29*(7.876233178997433d-1*t3*(-5.469606374303773d- 1848 5 2*t4*t76*t77-1.4585616998143394d-1*t31*t4*t75-7.424439106 1849 6 586571d-3*t74-1.3127055298329054d-1*t4*t52*t54)-5.4696063 1850 7 74303773d-2*t5*t76*t77-1.4585616998143394d-1*t31*t5*t75-3 1851 8 .9381165894987163d-1*t31*t32*t55-1.3127055298329054d-1*t5 1852 9 *t52*t54+1.9690582947493582d-1*t30*t33*t54+2.625411059665 1853 : 811d-1*t31*t33*t53)+1.763392293402975d-2*t54*t7*t77+2.351 1854 ; 1897245373004d-1*t3*t7*t75+7.053569173611901d-2*t31*t52*t 1855 < 7-7.617854707500852d+0*t2*t34**3/t6**4+7.617854707500852d 1856 = +0*t2*t34*t50*t56-6.34821225625071d-1*t29*t3*t32*t56+4.23 1857 > 21415041671406d-1*t29*t3*t34*t53+1.269642451250142d+0*t3* 1858 ? t32*t50*t51+1.0580353760417849d-1*t29*t30*t31*t34)-8.5486 1859 @ 27933030694d-3*t17*t73-1.865442d-1*t18*t40*t58*t64-8.2433 1860 1 19792565315d-3*t27*t36*t64+8.243319792565315d-3*t27*t58*t 1861 2 59+9.327209999999999d-2*t11*t35*t48*t57+3.96311967940178d 1862 3 -3*t27*t28*t57-3.96311967940178d-3*t27*t48*t49+1.09910930 1863 4 56753751d-2*t36*t40*t47-5.284159572535706d-3*t28*t35*t47) 1864 5 -3.37738d-2*t25*t41*t83-6.75476d-2*t25*t81*t82+2.41719309 1865 6 15078277d-3*t24*t73+1.013214d-1*t25*t45*t65*t71+2.3308647 1866 7 6681112d-3*t27*t41*t71-2.33086476681112d-3*t27*t65*t66-3. 1867 8 107819689081493d-3*t41*t45*t47)+3.37738d-2*t25*t41*t83+6. 1868 9 75476d-2*t25*t81*t82-2.4171930915078277d-3*t24*t73-1.0132 1869 : 14d-1*t25*t45*t65*t71-2.33086476681112d-3*t27*t41*t71+2.3 1870 ; 3086476681112d-3*t27*t65*t66+3.107819689081493d-3*t41*t45 1871 < *t47)-6.21814d-2*t18*t36*t80-1.243628d-1*t18*t78*t79+8.54 1872 = 8627933030694d-3*t17*t73+1.865442d-1*t18*t40*t58*t64+8.24 1873 > 3319792565315d-3*t27*t36*t64-8.243319792565315d-3*t27*t58 1874 ? *t59-1.0991093056753751d-2*t36*t40*t47)*wght+3.0d+0*t72*w 1875 @ ght+Amat3(iq,D3_RB_RB_RB) 1876 endif ! rhoa.gt.tol_rho.and.rhob.gt.tol_rho 1877 endif ! ipol.eq.1 1878 enddo ! iq 1879 end 1880C> @} 1881