1! dftd3 program for computing the dispersion energy and forces from cart 2! and atomic numbers as described in 3! 4! S. Grimme, J. Antony, S. Ehrlich and H. Krieg 5! J. Chem. Phys, 132 (2010), 154104 6! 7! S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem, 32 (2011), 1456 8! (for BJ-damping) 9! 10! Copyright (C) 2009 - 2011 Stefan Grimme, University of Muenster, Germany 11! 12! Repackaging of the original code without any change in the functionality: 13! 14! Copyright (C) 2016, Bálint Aradi 15! 16! This program is free software; you can redistribute it and/or modify 17! it under the terms of the GNU General Public License as published by 18! the Free Software Foundation; either version 1, or (at your option) 19! any later version. 20! 21! This program is distributed in the hope that it will be useful, 22! but WITHOUT ANY WARRANTY; without even the implied warranty of 23! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 24! GNU General Public License for more details. 25! 26! For the GNU General Public License, see <http://www.gnu.org/licenses/> 27! 28 29module dftd3_core 30 use dftd3_sizes 31 use dftd3_common 32 use dftd3_pars 33 implicit none 34 35 36 ! atomic <r^2>/<r^4> values 37 real(wp) r2r4(max_elem) 38 39 ! r2r4 =sqrt(0.5*r2r4(i)*dfloat(i)**0.5 ) with i=elementnumber 40 ! the large number of digits is just to keep the results consistent 41 ! with older versions. They should not imply any higher accuracy than 42 ! the old values 43 data r2r4 / & 44 &2.00734898_wp, 1.56637132_wp, 5.01986934_wp, 3.85379032_wp, 3.64446594_wp,& 45 &3.10492822_wp, 2.71175247_wp, 2.59361680_wp, 2.38825250_wp, 2.21522516_wp,& 46 &6.58585536_wp, 5.46295967_wp, 5.65216669_wp, 4.88284902_wp, 4.29727576_wp,& 47 &4.04108902_wp, 3.72932356_wp, 3.44677275_wp, 7.97762753_wp, 7.07623947_wp,& 48 &6.60844053_wp, 6.28791364_wp, 6.07728703_wp, 5.54643096_wp, 5.80491167_wp,& 49 &5.58415602_wp, 5.41374528_wp, 5.28497229_wp, 5.22592821_wp, 5.09817141_wp,& 50 &6.12149689_wp, 5.54083734_wp, 5.06696878_wp, 4.87005108_wp, 4.59089647_wp,& 51 &4.31176304_wp, 9.55461698_wp, 8.67396077_wp, 7.97210197_wp, 7.43439917_wp,& 52 &6.58711862_wp, 6.19536215_wp, 6.01517290_wp, 5.81623410_wp, 5.65710424_wp,& 53 &5.52640661_wp, 5.44263305_wp, 5.58285373_wp, 7.02081898_wp, 6.46815523_wp,& 54 &5.98089120_wp, 5.81686657_wp, 5.53321815_wp, 5.25477007_wp,11.02204549_wp,& 55 &0.15679528_wp, 9.35167836_wp, 9.06926079_wp, 8.97241155_wp, 8.90092807_wp,& 56 &8.85984840_wp, 8.81736827_wp, 8.79317710_wp, 7.89969626_wp, 8.80588454_wp,& 57 &8.42439218_wp, 8.54289262_wp, 8.47583370_wp, 8.45090888_wp, 8.47339339_wp,& 58 &7.83525634_wp, 8.20702843_wp, 7.70559063_wp, 7.32755997_wp, 7.03887381_wp,& 59 &6.68978720_wp, 6.05450052_wp, 5.88752022_wp, 5.70661499_wp, 5.78450695_wp,& 60 &7.79780729_wp, 7.26443867_wp, 6.78151984_wp, 6.67883169_wp, 6.39024318_wp,& 61 &6.09527958_wp,11.79156076_wp,11.10997644_wp, 9.51377795_wp, 8.67197068_wp,& 62 &8.77140725_wp, 8.65402716_wp, 8.53923501_wp, 8.85024712_wp / 63 64 ! PBE0/def2-QZVP atomic values 65 ! data r2r4 / 66 ! . 8.0589, 3.4698, 29.0974, 14.8517, 11.8799, 7.8715, 5.5588, 67 ! . 4.7566, 3.8025, 3.1036, 26.1552, 17.2304, 17.7210, 12.7442, 68 ! . 9.5361, 8.1652, 6.7463, 5.6004, 29.2012, 22.3934, 19.0598, 69 ! . 16.8590, 15.4023, 12.5589, 13.4788, 12.2309, 11.2809, 10.5569, 70 ! . 10.1428, 9.4907, 13.4606, 10.8544, 8.9386, 8.1350, 7.1251, 71 ! . 6.1971, 30.0162, 24.4103, 20.3537, 17.4780, 13.5528, 11.8451, 72 ! . 11.0355, 10.1997, 9.5414, 9.0061, 8.6417, 8.9975, 14.0834, 73 ! . 11.8333, 10.0179, 9.3844, 8.4110, 7.5152, 32.7622, 27.5708, 74 ! . 23.1671, 21.6003, 20.9615, 20.4562, 20.1010, 19.7475, 19.4828, 75 ! . 15.6013, 19.2362, 17.4717, 17.8321, 17.4237, 17.1954, 17.1631, 76 ! . 14.5716, 15.8758, 13.8989, 12.4834, 11.4421, 10.2671, 8.3549, 77 ! . 7.8496, 7.3278, 7.4820, 13.5124, 11.6554, 10.0959, 9.7340, 78 ! . 8.8584, 8.0125, 29.8135, 26.3157, 19.1885, 15.8542, 16.1305, 79 ! . 15.6161, 15.1226, 16.1576 / 80 81 82 ! scale r4/r2 values of the atoms by sqrt(Z) 83 ! sqrt is also globally close to optimum 84 ! together with the factor 1/2 this yield reasonable 85 ! c8 for he, ne and ar. for larger Z, C8 becomes too large 86 ! which effectively mimics higher R^n terms neglected due 87 ! to stability reasons 88 89 ! r2r4 =sqrt(0.5*r2r4(i)*dfloat(i)**0.5 ) with i=elementnumber 90 ! the large number of digits is just to keep the results consistent 91 ! with older versions. They should not imply any higher accuracy than 92 ! the old values 93 !!data r2r4 / & 94 !! & 2.00734898, 1.56637132, 5.01986934, 3.85379032, 3.64446594, & 95 !! & 3.10492822, 2.71175247, 2.59361680, 2.38825250, 2.21522516, & 96 !! & 6.58585536, 5.46295967, 5.65216669, 4.88284902, 4.29727576, & 97 !! & 4.04108902, 3.72932356, 3.44677275, 7.97762753, 7.07623947, & 98 !! & 6.60844053, 6.28791364, 6.07728703, 5.54643096, 5.80491167, & 99 !! & 5.58415602, 5.41374528, 5.28497229, 5.22592821, 5.09817141, & 100 !! & 6.12149689, 5.54083734, 5.06696878, 4.87005108, 4.59089647, & 101 !! & 4.31176304, 9.55461698, 8.67396077, 7.97210197, 7.43439917, & 102 !! & 6.58711862, 6.19536215, 6.01517290, 5.81623410, 5.65710424, & 103 !! & 5.52640661, 5.44263305, 5.58285373, 7.02081898, 6.46815523, & 104 !! & 5.98089120, 5.81686657, 5.53321815, 5.25477007, 11.02204549, & 105 !! &10.15679528, 9.35167836, 9.06926079, 8.97241155, 8.90092807, & 106 !! & 8.85984840, 8.81736827, 8.79317710, 7.89969626, 8.80588454, & 107 !! & 8.42439218, 8.54289262, 8.47583370, 8.45090888, 8.47339339, & 108 !! & 7.83525634, 8.20702843, 7.70559063, 7.32755997, 7.03887381, & 109 !! & 6.68978720, 6.05450052, 5.88752022, 5.70661499, 5.78450695, & 110 !! & 7.79780729, 7.26443867, 6.78151984, 6.67883169, 6.39024318, & 111 !! & 6.09527958, 11.79156076, 11.10997644, 9.51377795, 8.67197068, & 112 !! & 8.77140725, 8.65402716, 8.53923501, 8.85024712 / 113 114 115 real(wp) rcov(max_elem) 116 117 ! covalent radii 118 ! covalent radii (taken from Pyykko and Atsumi, Chem. Eur. J. 15, 2009, 119 ! values for metals decreased by 10 % 120 ! data rcov/ 121 ! . 0.32, 0.46, 1.20, 0.94, 0.77, 0.75, 0.71, 0.63, 0.64, 0.67 122 ! ., 1.40, 1.25, 1.13, 1.04, 1.10, 1.02, 0.99, 0.96, 1.76, 1.54 123 ! ., 1.33, 1.22, 1.21, 1.10, 1.07, 1.04, 1.00, 0.99, 1.01, 1.09 124 ! ., 1.12, 1.09, 1.15, 1.10, 1.14, 1.17, 1.89, 1.67, 1.47, 1.39 125 ! ., 1.32, 1.24, 1.15, 1.13, 1.13, 1.08, 1.15, 1.23, 1.28, 1.26 126 ! ., 1.26, 1.23, 1.32, 1.31, 2.09, 1.76, 1.62, 1.47, 1.58, 1.57 127 ! ., 1.56, 1.55, 1.51, 1.52, 1.51, 1.50, 1.49, 1.49, 1.48, 1.53 128 ! ., 1.46, 1.37, 1.31, 1.23, 1.18, 1.16, 1.11, 1.12, 1.13, 1.32 129 ! ., 1.30, 1.30, 1.36, 1.31, 1.38, 1.42, 2.01, 1.81, 1.67, 1.58 130 ! ., 1.52, 1.53, 1.54, 1.55 / 131 132 ! these new data are scaled with k2=4./3. and converted a_0 via 133 ! autoang=0.52917726d0 134 !!data rcov/ & 135 !! & 0.80628308, 1.15903197, 3.02356173, 2.36845659, 1.94011865, & 136 !! & 1.88972601, 1.78894056, 1.58736983, 1.61256616, 1.68815527, & 137 !! & 3.52748848, 3.14954334, 2.84718717, 2.62041997, 2.77159820, & 138 !! & 2.57002732, 2.49443835, 2.41884923, 4.43455700, 3.88023730, & 139 !! & 3.35111422, 3.07395437, 3.04875805, 2.77159820, 2.69600923, & 140 !! & 2.62041997, 2.51963467, 2.49443835, 2.54483100, 2.74640188, & 141 !! & 2.82199085, 2.74640188, 2.89757982, 2.77159820, 2.87238349, & 142 !! & 2.94797246, 4.76210950, 4.20778980, 3.70386304, 3.50229216, & 143 !! & 3.32591790, 3.12434702, 2.89757982, 2.84718717, 2.84718717, & 144 !! & 2.72120556, 2.89757982, 3.09915070, 3.22513231, 3.17473967, & 145 !! & 3.17473967, 3.09915070, 3.32591790, 3.30072128, 5.26603625, & 146 !! & 4.43455700, 4.08180818, 3.70386304, 3.98102289, 3.95582657, & 147 !! & 3.93062995, 3.90543362, 3.80464833, 3.82984466, 3.80464833, & 148 !! & 3.77945201, 3.75425569, 3.75425569, 3.72905937, 3.85504098, & 149 !! & 3.67866672, 3.45189952, 3.30072128, 3.09915070, 2.97316878, & 150 !! & 2.92277614, 2.79679452, 2.82199085, 2.84718717, 3.32591790, & 151 !! & 3.27552496, 3.27552496, 3.42670319, 3.30072128, 3.47709584, & 152 !! & 3.57788113, 5.06446567, 4.56053862, 4.20778980, 3.98102289, & 153 !! & 3.82984466, 3.85504098, 3.88023730, 3.90543362 / 154 155 ! these new data are scaled with k2=4./3. and converted a_0 via 156 ! autoang=0.52917726d0 157 data rcov/ & 158 & 0.80628308_wp, 1.15903197_wp, 3.02356173_wp, 2.36845659_wp, 1.94011865_wp, & 159 & 1.88972601_wp, 1.78894056_wp, 1.58736983_wp, 1.61256616_wp, 1.68815527_wp, & 160 & 3.52748848_wp, 3.14954334_wp, 2.84718717_wp, 2.62041997_wp, 2.77159820_wp, & 161 & 2.57002732_wp, 2.49443835_wp, 2.41884923_wp, 4.43455700_wp, 3.88023730_wp, & 162 & 3.35111422_wp, 3.07395437_wp, 3.04875805_wp, 2.77159820_wp, 2.69600923_wp, & 163 & 2.62041997_wp, 2.51963467_wp, 2.49443835_wp, 2.54483100_wp, 2.74640188_wp, & 164 & 2.82199085_wp, 2.74640188_wp, 2.89757982_wp, 2.77159820_wp, 2.87238349_wp, & 165 & 2.94797246_wp, 4.76210950_wp, 4.20778980_wp, 3.70386304_wp, 3.50229216_wp, & 166 & 3.32591790_wp, 3.12434702_wp, 2.89757982_wp, 2.84718717_wp, 2.84718717_wp, & 167 & 2.72120556_wp, 2.89757982_wp, 3.09915070_wp, 3.22513231_wp, 3.17473967_wp, & 168 & 3.17473967_wp, 3.09915070_wp, 3.32591790_wp, 3.30072128_wp, 5.26603625_wp, & 169 & 4.43455700_wp, 4.08180818_wp, 3.70386304_wp, 3.98102289_wp, 3.95582657_wp, & 170 & 3.93062995_wp, 3.90543362_wp, 3.80464833_wp, 3.82984466_wp, 3.80464833_wp, & 171 & 3.77945201_wp, 3.75425569_wp, 3.75425569_wp, 3.72905937_wp, 3.85504098_wp, & 172 & 3.67866672_wp, 3.45189952_wp, 3.30072128_wp, 3.09915070_wp, 2.97316878_wp, & 173 & 2.92277614_wp, 2.79679452_wp, 2.82199085_wp, 2.84718717_wp, 3.32591790_wp, & 174 & 3.27552496_wp, 3.27552496_wp, 3.42670319_wp, 3.30072128_wp, 3.47709584_wp, & 175 & 3.57788113_wp, 5.06446567_wp, 4.56053862_wp, 4.20778980_wp, 3.98102289_wp, & 176 & 3.82984466_wp, 3.85504098_wp, 3.88023730_wp, 3.90543362_wp / 177 178 179contains 180 181 182 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 183 ! set parameters 184 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 185 186 subroutine setfuncpar(func,version,TZ,s6,rs6,s18,rs18,alp) 187 integer version 188 real(wp) s6,rs6,s18,alp,rs18 189 character*(*) func 190 logical TZ 191 ! double hybrid values revised according to procedure in the GMTKN30 pap 192 193 194 if(version.eq.6)then 195 s6 =1.0d0 196 alp =14.0d0 197! BJ damping with parameters from ... 198 select case (func) 199 case ("b2-plyp") 200 rs6 =0.486434 201 s18 =0.672820 202 rs18=3.656466 203 s6 =0.640000 204 case ("b3-lyp") 205 rs6 =0.278672 206 s18 =1.466677 207 rs18=4.606311 208 case ("b97-d") 209 rs6 =0.240184 210 s18 =1.206988 211 rs18=3.864426 212 case ("b-lyp") 213 rs6 =0.448486 214 s18 =1.875007 215 rs18=3.610679 216 case ("b-p") 217 rs6 =0.821850 218 s18 =3.140281 219 rs18=2.728151 220 case ("pbe") 221 rs6 =0.012092 222 s18 =0.358940 223 rs18=5.938951 224 case ("pbe0") 225 rs6 =0.007912 226 s18 =0.528823 227 rs18=6.162326 228 case ("lc-wpbe") 229 rs6 =0.563761 230 s18 =0.906564 231 rs18=3.593680 232 case DEFAULT 233 call stoprun( 'functional name unknown' ) 234 end select 235 endif 236 237 if(version.eq.5)then 238 s6 =1.0d0 239 alp =14.0d0 240! zero damping with parameters from ... 241 select case (func) 242 case ("b2-plyp") 243 rs6 =1.313134 244 s18 =0.717543 245 rs18=0.016035 246 s6 =0.640000 247 case ("b3-lyp") 248 rs6 =1.338153 249 s18 =1.532981 250 rs18=0.013988 251 case ("b97-d") 252 rs6 =1.151808 253 s18 =1.020078 254 rs18=0.035964 255 case ("b-lyp") 256 rs6 =1.279637 257 s18 =1.841686 258 rs18=0.014370 259 case ("b-p") 260 rs6 =1.233460 261 s18 =1.945174 262 rs18=0.000000 263 case ("pbe") 264 rs6 =2.340218 265 s18 =0.000000 266 rs18=0.129434 267 case ("pbe0") 268 rs6 =2.077949 269 s18 =0.000081 270 rs18=0.116755 271 case ("lc-wpbe") 272 rs6 =1.366361 273 s18 =1.280619 274 rs18=0.003160 275 case DEFAULT 276 call stoprun( 'functional name unknown' ) 277 end select 278 endif 279 280 ! DFT-D3 with Becke-Johnson finite-damping, variant 2 with their radii 281 ! SE: Alp is only used in 3-body calculations 282 if (version.eq.4)then 283 s6=1.0d0 284 alp =14.0d0 285 286 select case (func) 287 case ("b-p") 288 rs6 =0.3946 289 s18 =3.2822 290 rs18=4.8516 291 case ("b-lyp") 292 rs6 =0.4298 293 s18 =2.6996 294 rs18=4.2359 295 case ("revpbe") 296 rs6 =0.5238 297 s18 =2.3550 298 rs18=3.5016 299 case ("rpbe") 300 rs6 =0.1820 301 s18 =0.8318 302 rs18=4.0094 303 case ("b97-d") 304 rs6 =0.5545 305 s18 =2.2609 306 rs18=3.2297 307 case ("pbe") 308 rs6 =0.4289 309 s18 =0.7875 310 rs18=4.4407 311 case ("rpw86-pbe") 312 rs6 =0.4613 313 s18 =1.3845 314 rs18=4.5062 315 case ("b3-lyp") 316 rs6 =0.3981 317 s18 =1.9889 318 rs18=4.4211 319 case ("tpss") 320 rs6 =0.4535 321 s18 =1.9435 322 rs18=4.4752 323 case ("hf") 324 rs6 =0.3385 325 s18 =0.9171 326 rs18=2.8830 327 case ("tpss0") 328 rs6 =0.3768 329 s18 =1.2576 330 rs18=4.5865 331 case ("pbe0") 332 rs6 =0.4145 333 s18 =1.2177 334 rs18=4.8593 335 case ("hse06") 336 rs6 =0.383 337 s18 =2.310 338 rs18=5.685 339 case ("revpbe38") 340 rs6 =0.4309 341 s18 =1.4760 342 rs18=3.9446 343 case ("pw6b95") 344 rs6 =0.2076 345 s18 =0.7257 346 rs18=6.3750 347 case ("b2-plyp") 348 rs6 =0.3065 349 s18 =0.9147 350 rs18=5.0570 351 s6=0.64d0 352 case ("dsd-blyp") 353 rs6 =0.0000 354 s18 =0.2130 355 rs18=6.0519 356 s6=0.50d0 357 case ("dsd-blyp-fc") 358 rs6 =0.0009 359 s18 =0.2112 360 rs18=5.9807 361 s6=0.50d0 362 case ("bop") 363 rs6 =0.4870 364 s18 =3.2950 365 rs18=3.5043 366 case ("mpwlyp") 367 rs6 =0.4831 368 s18 =2.0077 369 rs18=4.5323 370 case ("o-lyp") 371 rs6 =0.5299 372 s18 =2.6205 373 rs18=2.8065 374 case ("pbesol") 375 rs6 =0.4466 376 s18 =2.9491 377 rs18=6.1742 378 case ("bpbe") 379 rs6 =0.4567 380 s18 =4.0728 381 rs18=4.3908 382 case ("opbe") 383 rs6 =0.5512 384 s18 =3.3816 385 rs18=2.9444 386 case ("ssb") 387 rs6 =-0.0952 388 s18 =-0.1744 389 rs18=5.2170 390 case ("revssb") 391 rs6 =0.4720 392 s18 =0.4389 393 rs18=4.0986 394 case ("otpss") 395 rs6 =0.4634 396 s18 =2.7495 397 rs18=4.3153 398 case ("b3pw91") 399 rs6 =0.4312 400 s18 =2.8524 401 rs18=4.4693 402 case ("bh-lyp") 403 rs6 =0.2793 404 s18 =1.0354 405 rs18=4.9615 406 case ("revpbe0") 407 rs6 =0.4679 408 s18 =1.7588 409 rs18=3.7619 410 case ("tpssh") 411 rs6 =0.4529 412 s18 =2.2382 413 rs18=4.6550 414 case ("mpw1b95") 415 rs6 =0.1955 416 s18 =1.0508 417 rs18=6.4177 418 case ("pwb6k") 419 rs6 =0.1805 420 s18 =0.9383 421 rs18=7.7627 422 case ("b1b95") 423 rs6 =0.2092 424 s18 =1.4507 425 rs18=5.5545 426 case ("bmk") 427 rs6 =0.1940 428 s18 =2.0860 429 rs18=5.9197 430 case ("cam-b3lyp") 431 rs6 =0.3708 432 s18 =2.0674 433 rs18=5.4743 434 case ("lc-wpbe") 435 rs6 =0.3919 436 s18 =1.8541 437 rs18=5.0897 438 case ("b2gp-plyp") 439 rs6 =0.0000 440 s18 =0.2597 441 rs18=6.3332 442 s6=0.560 443 case ("ptpss") 444 rs6 =0.0000 445 s18 =0.2804 446 rs18=6.5745 447 s6=0.750 448 case ("pwpb95") 449 rs6 =0.0000 450 s18 =0.2904 451 rs18=7.3141 452 s6=0.820 453 ! special HF/DFT with eBSSE correction 454 case ("hf/mixed") 455 rs6 =0.5607 456 s18 =3.9027 457 rs18=4.5622 458 case ("hf/sv") 459 rs6 =0.4249 460 s18 =2.1849 461 rs18=4.2783 462 case ("hf/minis") 463 rs6 =0.1702 464 s18 =0.9841 465 rs18=3.8506 466 case ("b3-lyp/6-31gd") 467 rs6 =0.5014 468 s18 =4.0672 469 rs18=4.8409 470 case ("hcth120") 471 rs6=0.3563 472 s18=1.0821 473 rs18=4.3359 474 ! DFTB3 old, deprecated parameters: 475 ! case ("dftb3") 476 ! rs6=0.7461 477 ! s18=3.209 478 ! rs18=4.1906 479 ! special SCC-DFTB parametrization 480 ! full third order DFTB, self consistent charges, hydrogen pair damping 481 ! exponent 4.2 482 case("dftb3") 483 rs6=0.5719d0 484 s18=0.5883d0 485 rs18=3.6017d0 486 case ("pw1pw") 487 rs6 =0.3807d0 488 s18 =2.3363d0 489 rs18=5.8844d0 490 case ("pwgga") 491 rs6 =0.2211d0 492 s18 =2.6910d0 493 rs18=6.7278d0 494 case ("hsesol") 495 rs6 =0.4650d0 496 s18 =2.9215d0 497 rs18=6.2003d0 498 ! special HF-D3-gCP-SRB/MINIX parametrization 499 case ("hf3c") 500 rs6=0.4171d0 501 s18=0.8777d0 502 rs18=2.9149d0 503 ! special HF-D3-gCP-SRB2/ECP-2G parametrization 504 case ("hf3cv") 505 rs6=0.3063d0 506 s18=0.5022d0 507 rs18=3.9856d0 508 ! special PBEh-D3-gCP/def2-mSVP parametrization 509 case ("pbeh3c", "pbeh-3c") 510 rs6=0.4860d0 511 s18=0.0000d0 512 rs18=4.5000d0 513 514 case DEFAULT 515 call stoprun( 'functional name unknown' ) 516 end select 517 end if 518 519 ! DFT-D3 520 if (version.eq.3)then 521 s6 =1.0d0 522 alp =14.0d0 523 rs18=1.0d0 524 ! default def2-QZVP (almost basis set limit) 525 if (.not.TZ) then 526 select case (func) 527 case ("slater-dirac-exchange") 528 rs6 =0.999 529 s18 =-1.957 530 rs18=0.697 531 case ("b-lyp") 532 rs6=1.094 533 s18=1.682 534 case ("b-p") 535 rs6=1.139 536 s18=1.683 537 case ("b97-d") 538 rs6=0.892 539 s18=0.909 540 case ("revpbe") 541 rs6=0.923 542 s18=1.010 543 case ("pbe") 544 rs6=1.217 545 s18=0.722 546 case ("pbesol") 547 rs6=1.345 548 s18=0.612 549 case ("rpw86-pbe") 550 rs6=1.224 551 s18=0.901 552 case ("rpbe") 553 rs6=0.872 554 s18=0.514 555 case ("tpss") 556 rs6=1.166 557 s18=1.105 558 case ("b3-lyp") 559 rs6=1.261 560 s18=1.703 561 case ("pbe0") 562 rs6=1.287 563 s18=0.928 564 565 case ("hse06") 566 rs6=1.129 567 s18=0.109 568 case ("revpbe38") 569 rs6=1.021 570 s18=0.862 571 case ("pw6b95") 572 rs6=1.532 573 s18=0.862 574 case ("tpss0") 575 rs6=1.252 576 s18=1.242 577 case ("b2-plyp") 578 rs6=1.427 579 s18=1.022 580 s6=0.64 581 case ("pwpb95") 582 rs6=1.557 583 s18=0.705 584 s6=0.82 585 case ("b2gp-plyp") 586 rs6=1.586 587 s18=0.760 588 s6=0.56 589 case ("ptpss") 590 rs6=1.541 591 s18=0.879 592 s6=0.75 593 case ("hf") 594 rs6=1.158 595 s18=1.746 596 case ("mpwlyp") 597 rs6=1.239 598 s18=1.098 599 case ("bpbe") 600 rs6=1.087 601 s18=2.033 602 case ("bh-lyp") 603 rs6=1.370 604 s18=1.442 605 case ("tpssh") 606 rs6=1.223 607 s18=1.219 608 case ("pwb6k") 609 rs6=1.660 610 s18=0.550 611 case ("b1b95") 612 rs6=1.613 613 s18=1.868 614 case ("bop") 615 rs6=0.929 616 s18=1.975 617 case ("o-lyp") 618 rs6=0.806 619 s18=1.764 620 case ("o-pbe") 621 rs6=0.837 622 s18=2.055 623 case ("ssb") 624 rs6=1.215 625 s18=0.663 626 case ("revssb") 627 rs6=1.221 628 s18=0.560 629 case ("otpss") 630 rs6=1.128 631 s18=1.494 632 case ("b3pw91") 633 rs6=1.176 634 s18=1.775 635 case ("revpbe0") 636 rs6=0.949 637 s18=0.792 638 case ("pbe38") 639 rs6=1.333 640 s18=0.998 641 case ("mpw1b95") 642 rs6=1.605 643 s18=1.118 644 case ("mpwb1k") 645 rs6=1.671 646 s18=1.061 647 case ("bmk") 648 rs6=1.931 649 s18=2.168 650 case ("cam-b3lyp") 651 rs6=1.378 652 s18=1.217 653 case ("lc-wpbe") 654 rs6=1.355 655 s18=1.279 656 case ("m05") 657 rs6=1.373 658 s18=0.595 659 case ("m052x") 660 rs6=1.417 661 s18=0.000 662 case ("m06l") 663 rs6=1.581 664 s18=0.000 665 case ("m06") 666 rs6=1.325 667 s18=0.000 668 case ("m062x") 669 rs6=1.619 670 s18=0.000 671 case ("m06hf") 672 rs6=1.446 673 s18=0.000 674 ! DFTB3 (zeta=4.0), old deprecated parameters 675 ! case ("dftb3") 676 ! rs6=1.235 677 ! s18=0.673 678 case ("hcth120") 679 rs6=1.221 680 s18=1.206 681 case DEFAULT 682 call stoprun( 'functional name unknown' ) 683 end select 684 else 685 ! special TZVPP parameter 686 select case (func) 687 case ("b-lyp") 688 rs6=1.243 689 s18=2.022 690 case ("b-p") 691 rs6=1.221 692 s18=1.838 693 case ("b97-d") 694 rs6=0.921 695 s18=0.894 696 case ("revpbe") 697 rs6=0.953 698 s18=0.989 699 case ("pbe") 700 rs6=1.277 701 s18=0.777 702 case ("tpss") 703 rs6=1.213 704 s18=1.176 705 case ("b3-lyp") 706 rs6=1.314 707 s18=1.706 708 case ("pbe0") 709 rs6=1.328 710 s18=0.926 711 case ("pw6b95") 712 rs6=1.562 713 s18=0.821 714 case ("tpss0") 715 rs6=1.282 716 s18=1.250 717 case ("b2-plyp") 718 rs6=1.551 719 s18=1.109 720 s6=0.5 721 case DEFAULT 722 call stoprun( 'functional name unknown (TZ case)' ) 723 end select 724 end if 725 end if 726 ! DFT-D2 727 if (version.eq.2)then 728 rs6=1.1d0 729 s18=0.0d0 730 alp=20.0d0 731 select case (func) 732 case ("b-lyp") 733 s6=1.2 734 case ("b-p") 735 s6=1.05 736 case ("b97-d") 737 s6=1.25 738 case ("revpbe") 739 s6=1.25 740 case ("pbe") 741 s6=0.75 742 case ("tpss") 743 s6=1.0 744 case ("b3-lyp") 745 s6=1.05 746 case ("pbe0") 747 s6=0.6 748 case ("pw6b95") 749 s6=0.5 750 case ("tpss0") 751 s6=0.85 752 case ("b2-plyp") 753 s6=0.55 754 case ("b2gp-plyp") 755 s6=0.4 756 case ("dsd-blyp") 757 s6=0.41 758 alp=60.0d0 759 case DEFAULT 760 call stoprun( 'functional name unknown' ) 761 end select 762 763 end if 764 765 end subroutine setfuncpar 766 767 768 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 769 ! compute energy 770 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 771 772 subroutine edisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & 773 & rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,rthr,cn_thr,& 774 & e6,e8,e10,e12,e63) 775 integer n,iz(*),max_elem,maxc,version,mxc(max_elem) 776 real(wp) xyz(3,*),r0ab(max_elem,max_elem),r2r4(*) 777 real(wp) rs6,rs8,rs10,alp6,alp8,alp10,rcov(max_elem) 778 real(wp) c6ab(max_elem,max_elem,maxc,maxc,3) 779 real(wp) e6, e8, e10, e12, e63 780 logical noabc 781 782 integer iat,jat,kat 783 real(wp) r,r2,r6,r8,tmp,dx,dy,dz,c6,c8,ang,rav 784 real(wp) damp6,damp8,rr,c9,rthr,cn_thr 785 real(wp) cn(n) 786 real(wp) r2ab(n*n),cc6ab(n*n),dmp(n*n),t1,t2,t3,a1,a2,tmp2 787 real(wp) abcthr 788 integer(int64) icomp(n*n) 789 integer ij,ik,jk 790 791 e6 =0 792 e8 =0 793 e10=0 794 e12=0 795 e63=0 796 ! the threebody term uses the same threshold as the 797 abcthr=cn_thr 798 799 ! Becke-Johnson parameters 800 a1=rs6 801 a2=rs8 802 803 804 ! DFT-D2 805 if (version.eq.2)then 806 807 do iat=1,n-1 808 do jat=iat+1,n 809 dx=xyz(1,iat)-xyz(1,jat) 810 dy=xyz(2,iat)-xyz(2,jat) 811 dz=xyz(3,iat)-xyz(3,jat) 812 r2=dx*dx+dy*dy+dz*dz 813 ! if (r2.gt.rthr) cycle 814 r=sqrt(r2) 815 c6=c6ab(iz(jat),iz(iat),1,1,1) 816 damp6=1./(1.+exp(-alp6*(r/(rs6*r0ab(iz(jat),iz(iat)))-1.))) 817 r6=r2**3 818 e6 =e6+c6*damp6/r6 819 end do 820 end do 821 822 else 823 ! DFT-D3 824 call ncoord(n,rcov,iz,xyz,cn,cn_thr) 825 826 icomp=0 827 do iat=1,n-1 828 do jat=iat+1,n 829 dx=xyz(1,iat)-xyz(1,jat) 830 dy=xyz(2,iat)-xyz(2,jat) 831 dz=xyz(3,iat)-xyz(3,jat) 832 r2=dx*dx+dy*dy+dz*dz 833 !THR 834 if (r2.gt.rthr) cycle 835 r =sqrt(r2) 836 tmp2=r0ab(iz(jat),iz(iat)) 837 rr=tmp2/r 838 ! damping 839 if(version.eq.3)then 840 ! DFT-D3 zero-damp 841 tmp=rs6*rr 842 damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) 843 tmp=rs8*rr 844 damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) 845 else 846 ! DFT-D3M zero-damp 847 tmp=(r/(rs6*tmp2))+rs8*tmp2 848 damp6 =1.d0/( 1.d0+6.d0*tmp**(-alp6) ) 849 tmp=(r/tmp2)+rs8*tmp2 850 damp8 =1.d0/( 1.d0+6.d0*tmp**(-alp8) ) 851 endif 852 ! get C6 853 call getc6(maxc,max_elem,c6ab,mxc,iz(iat),iz(jat), & 854 & cn(iat),cn(jat),c6) 855 856 r6=r2**3 857 r8=r6*r2 858 ! r2r4 stored in main as sqrt 859 c8 =3.0d0*c6*r2r4(iz(iat))*r2r4(iz(jat)) 860 861 ! DFT-D3 zero-damp or DFT-D3 M(zero) 862 if((version.eq.3).or.(version.eq.5))then 863 e6=e6+c6*damp6/r6 864 e8=e8+c8*damp8/r8 865 end if 866 ! DFT-D3(BJ) or DFT-D3M(BJ) 867 if((version.eq.4).or.(version.eq.6))then 868 ! use BJ radius 869 tmp=sqrt(c8/c6) 870 e6=e6+ c6/(r6+(a1*tmp+a2)**6) 871 e8=e8+ c8/(r8+(a1*tmp+a2)**8) 872 end if 873 874 ! if (.not.noabc) then 875 if ((.not.noabc).and.(r2.lt.abcthr)) then 876 ij=lin(jat,iat) 877 icomp(ij)=1 878 ! store C6 for C9, calc as sqrt 879 cc6ab(ij)=sqrt(c6) 880 ! store R^2 for abc 881 r2ab(ij)=r2 882 ! store for abc damping 883 dmp(ij)=(1./rr)**(1./3.) 884 !noabc 885 end if 886 end do 887 end do 888 889 if (noabc)return 890 891 ! compute non-additive third-order energy using averaged C6 892 do iat=1,n 893 do jat=1,iat-1 894 ij=lin(jat,iat) 895 if (icomp(ij).eq.1)then 896 do kat=1,jat-1 897 898 ik=lin(kat,iat) 899 jk=lin(kat,jat) 900 if ((icomp(ik).eq.0).or.(icomp(jk).eq.0)) cycle 901 ! damping func product 902 ! tmp=dmp(ik)*dmp(jk)*dmp(ij) 903 rav=(4./3.)/(dmp(ik)*dmp(jk)*dmp(ij)) 904 tmp=1.d0/( 1.d0+6.d0*rav**alp8 ) 905 ! triple C6 coefficient (stored as sqrt) 906 c9=cc6ab(ij)*cc6ab(ik)*cc6ab(jk) 907 ! write(*,*) 'C9', c9 908 ! angular terms with law of cosines 909 910 t1 = (r2ab(ij)+r2ab(jk)-r2ab(ik)) 911 t2 = (r2ab(ij)+r2ab(ik)-r2ab(jk)) 912 t3 = (r2ab(ik)+r2ab(jk)-r2ab(ij)) 913 tmp2=r2ab(ij)*r2ab(jk)*r2ab(ik) 914 ang=(0.375d0*t1*t2*t3/tmp2+1.0d0)/tmp2**1.5d0 915 916 e63=e63-tmp*c9*ang 917 918 end do 919 end if 920 end do 921 end do 922 923 end if 924 925 end subroutine edisp 926 927 928 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 929 ! compute gradient 930 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 931 932 subroutine gdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & 933 & s6,s18,rs6,rs8,rs10,alp6,alp8,alp10,noabc,rthr, & 934 & num,version,echo,g,disp,gnorm,cn_thr,fix) 935 936 integer n,iz(*),max_elem,maxc,version,mxc(max_elem) 937 real(wp) xyz(3,*),r0ab(max_elem,max_elem),r2r4(*) 938 real(wp) c6ab(max_elem,max_elem,maxc,maxc,3) 939 real(wp) g(3,*),s6,s18,rcov(max_elem) 940 real(wp) rs6,rs8,rs10,alp10,alp8,alp6,a1,a2 941 logical noabc,num,echo,fix(n) 942 943 integer iat,jat,i,j,kat 944 real(wp) R0,C6,R42,disp,x1,e6abc 945 real(wp) dx,dy,dz,r2,r,r4,r6,r8,t6,t8,damp1 946 real(wp) damp6,damp8,e6,e8,e10,e12,gnorm,tmp1 947 real(wp) s10,s8,step,dispr,displ,r235,tmp2 948 real(wp) cn(n),gx1,gy1,gz1,gx2,gy2,gz2,rthr,cn_thr 949 950 real(wp) rij(3),r7,r9 951 !d(E)/d(r_ij) derivative wrt. dist. iat-ja 952 real(wp) drij(n*(n+1)/2) 953 real(wp) rcovij 954 !d(C6ij)/d(r_ij) 955 real(wp) expterm 956 !dCN(iat)/d(r_ij) is equal to 957 real(wp) dcn 958 !dCN(jat)/d(r_ij) 959 ! saves (1/r^6*f_dmp + 3*r4r2/r^8*f_dmp) for kat l 960 real(wp) dc6_rest 961 integer linij,linik,linjk 962 ! dE_disp/dCN(iat) in dc6i(iat) 963 real(wp) dc6i(n) 964 ! saves dC6(ij)/dCN(iat) 965 real(wp) dc6ij(n,n) 966 real(wp) dc6iji,dc6ijj 967 logical abccalc(n*(n+1)/2) 968 real(wp) abcthr 969 ! temporary container to store iat gradient 970 real(wp) gtmp(3) 971 real(wp) labc,rabc 972 real(wp) c6abc(n*(n+1)/2) 973 real(wp) r2abc(n*(n+1)/2) 974 real(wp) r3abc(n*(n+1)/2) 975 real(wp) c9,rav,rav3,fdmp,ang,angr9,eabc,dc9,dfdmp,dang 976 real(wp) r2ij,r2jk,r2ik,mijk,imjk,ijmk,rijk3 977 integer kk 978 979 980 dc6i=0.0d0 981 abccalc=.FALSE. 982 abcthr=cn_thr 983 984 985 986 987 !NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN 988 if (num) then 989 if (echo)write(*,*) 'doing numerical gradient O(N^3) ...' 990 991 call edisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & 992 & rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,rthr,cn_thr, & 993 & e6,e8,e10,e12,e6abc) 994 disp=-s6*e6-s18*e8-e6abc 995 996 step=2.d-5 997 998 do i=1,n 999 do j=1,3 1000 xyz(j,i)=xyz(j,i)+step 1001 call edisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & 1002 & rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,rthr,cn_thr, & 1003 & e6,e8,e10,e12,e6abc) 1004 dispr=-s6*e6-s18*e8-e6abc 1005 rabc=e6abc 1006 xyz(j,i)=xyz(j,i)-2*step 1007 call edisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & 1008 & rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,rthr,cn_thr, & 1009 & e6,e8,e10,e12,e6abc) 1010 displ=-s6*e6-s18*e8-e6abc 1011 labc=e6abc 1012 g(j,i)=0.5*(dispr-displ)/step 1013 xyz(j,i)=xyz(j,i)+step 1014 end do 1015 end do 1016 goto 999 1017 end if 1018 1019 ! this is the crucial threshold to reduce the N^3 to an 1020 ! effective N^2. 1021 1022 ! write(*,*)'rthr=',rthr,'rthr2=',rthr2,'rthr3=',rthr3 1023 1024 if (echo)write(*,*) 1025 ! 2222222222222222222222222222222222222222222222222222222222222222222222 1026 if (version.eq.2)then 1027 if (echo)write(*,*) 'doing analytical gradient O(N^2) ...' 1028 disp=0 1029 do iat=1,n-1 1030 do jat=iat+1,n 1031 R0=r0ab(iz(jat),iz(iat))*rs6 1032 dx=(xyz(1,iat)-xyz(1,jat)) 1033 dy=(xyz(2,iat)-xyz(2,jat)) 1034 dz=(xyz(3,iat)-xyz(3,jat)) 1035 r2 =dx*dx+dy*dy+dz*dz 1036 ! if (r2.gt.rthr) cycle 1037 r235=r2**3.5 1038 r =dsqrt(r2) 1039 damp6=exp(-alp6*(r/R0-1.0d0)) 1040 damp1=1.+damp6 1041 c6=c6ab(iz(jat),iz(iat),1,1,1)*s6 1042 tmp1=damp6/(damp1*damp1*r235*R0) 1043 tmp2=6./(damp1*r*r235) 1044 gx1=alp6* dx*tmp1-tmp2*dx 1045 gx2=alp6*(-dx)*tmp1+tmp2*dx 1046 gy1=alp6* dy*tmp1-tmp2*dy 1047 gy2=alp6*(-dy)*tmp1+tmp2*dy 1048 gz1=alp6* dz*tmp1-tmp2*dz 1049 gz2=alp6*(-dz)*tmp1+tmp2*dz 1050 g(1,iat)=g(1,iat)-gx1*c6 1051 g(2,iat)=g(2,iat)-gy1*c6 1052 g(3,iat)=g(3,iat)-gz1*c6 1053 g(1,jat)=g(1,jat)-gx2*c6 1054 g(2,jat)=g(2,jat)-gy2*c6 1055 g(3,jat)=g(3,jat)-gz2*c6 1056 disp=disp+c6*(1./damp1)/r2**3 1057 end do 1058 end do 1059 disp=-disp 1060 1061 1062 goto 999 1063 1064 ! 3333333333333333333333333333333333333333333333333333333333333333333333 1065 ! zero damping 1066 elseif ((version.eq.3).or.(version.eq.5)) then 1067 1068 if (echo)write(*,*) 'doing analytical gradient O(N^2) ...' 1069 call ncoord(n,rcov,iz,xyz,cn,cn_thr) 1070 s8 =s18 1071 s10=s18 1072 1073 disp=0 1074 1075 drij=0.0d0 1076 dc6_rest=0.0d0 1077 kat=0 1078 1079 1080 kk=0 1081 do iat=1,n 1082 do jat=1,iat-1 1083 rij=xyz(:,jat)-xyz(:,iat) 1084 r2=sum(rij*rij) 1085 if (r2.gt.rthr) cycle 1086 linij=lin(iat,jat) 1087 1088 r0=r0ab(iz(jat),iz(iat)) 1089 r42=r2r4(iz(iat))*r2r4(iz(jat)) 1090 ! rcovij=rcov(iz(iat))+rcov(iz(jat)) 1091 ! 1092 ! get_dC6_dCNij calculates the derivative dC6(iat,jat)/dCN(iat) and 1093 ! dC6(iat,jat)/dCN(jat). 1094 ! 1095 call get_dC6_dCNij(maxc,max_elem,c6ab,mxc(iz(iat)), & 1096 & mxc(iz(jat)),cn(iat),cn(jat),iz(iat),iz(jat),iat,jat, & 1097 & c6,dc6iji,dc6ijj) 1098 1099 1100 r=dsqrt(r2) 1101 r6=r2*r2*r2 1102 r7=r6*r 1103 r8=r6*r2 1104 r9=r8*r 1105 1106 ! 1107 ! Calculates damping functions: 1108 if (version.eq.3) then 1109 t6 = (r/(rs6*R0))**(-alp6) 1110 damp6 =1.d0/( 1.d0+6.d0*t6 ) 1111 t8 = (r/(rs8*R0))**(-alp8) 1112 damp8 =1.d0/( 1.d0+6.d0*t8 ) 1113 1114 tmp1=s6*6.d0*damp6*C6/r7 1115 tmp2=s8*6.d0*C6*r42*damp8/r9 1116 ! d(r^(-6))/d(r_ij) 1117 drij(linij)=drij(linij)-tmp1 & 1118 & -4.d0*tmp2 1119 1120 1121 drij(linij)=drij(linij) & 1122 & +tmp1*alp6*t6*damp6 & 1123 & +3.d0*tmp2*alp8*t8*damp8 1124 !d(f_dmp)/d(r_ij) 1125 1126 else ! version.eq.5 1127 t6 = (r/(rs6*R0)+R0*rs8)**(-alp6) 1128 damp6 =1.d0/( 1.d0+6.d0*t6 ) 1129 t8 = (r/R0+R0*rs8)**(-alp8) 1130 damp8 =1.d0/( 1.d0+6.d0*t8 ) 1131 1132 tmp1=s6*6.d0*damp6*C6/r7 1133 tmp2=s8*6.d0*C6*r42*damp8/r9 1134 ! d(r^(-6))/d(r_ij) 1135 drij(linij)=drij(linij)-tmp1 & 1136 & -4.d0*tmp2 1137 1138 1139 drij(linij)=drij(linij) & 1140 & +tmp1*alp6*t6*damp6*r/(r+rs6*R0*R0*rs8) & 1141 & +3.d0*tmp2*alp8*t8*damp8*r/(r+R0*R0*rs8) 1142 !d(f_dmp)/d(r_ij) 1143 1144 end if 1145 1146 1147 if ((.not.noabc).and.(r2.lt.abcthr)) then 1148 ! if (.not.noabc) then 1149 abccalc(linij)=.TRUE. 1150 dc6ij(iat,jat)=dc6iji 1151 dc6ij(jat,iat)=dc6ijj 1152 c6abc(linij)=c6 1153 r2abc(linij)=r2 1154 r3abc(linij)=(r/R0)**(1.0/3.0) 1155 !noabc 1156 end if 1157 1158 dc6_rest=s6/r6*damp6+3.d0*s8*r42/r8*damp8 1159 1160 ! calculate E_disp for sanity check 1161 disp=disp-dc6_rest*c6 1162 ! 1163 ! 1164 ! saving all f_dmp/r6*dC6(ij)/dCN(i) for each atom for later 1165 dc6i(iat)=dc6i(iat)+dc6_rest*dc6iji 1166 dc6i(jat)=dc6i(jat)+dc6_rest*dc6ijj 1167 1168 1169 !jat 1170 end do 1171 !iat 1172 end do 1173 1174 1175 1176 1177 1178 1179 ! end if !version 3+5 1180 1181 ! BJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJ 1182 ! Becke-Johnson finite damping 1183 elseif ((version.eq.4).or.(version.eq.6)) then 1184 a1 =rs6 1185 a2 =rs8 1186 s8 =s18 1187 1188 if (echo)write(*,*) 'doing analytical gradient O(N^2) ...' 1189 disp=0 1190 call ncoord(n,rcov,iz,xyz,cn,cn_thr) 1191 1192 drij=0.0d0 1193 dc6_rest=0.0d0 1194 kat=0 1195 1196 1197 do iat=1,n 1198 do jat=1,iat-1 1199 rij=xyz(:,jat)-xyz(:,iat) 1200 r2=sum(rij*rij) 1201 if (r2.gt.rthr) cycle 1202 1203 linij=lin(iat,jat) 1204 r0=r0ab(iz(jat),iz(iat)) 1205 r42=r2r4(iz(iat))*r2r4(iz(jat)) 1206 ! 1207 ! get_dC6_dCNij calculates the derivative dC6(iat,jat)/dCN(iat) and 1208 ! dC6(iat,jat)/dCN(jat). 1209 ! 1210 call get_dC6_dCNij(maxc,max_elem,c6ab,mxc(iz(iat)), & 1211 & mxc(iz(jat)),cn(iat),cn(jat),iz(iat),iz(jat),iat,jat, & 1212 & c6,dc6iji,dc6ijj) 1213 1214 1215 r=dsqrt(r2) 1216 r4=r2*r2 1217 r6=r4*r2 1218 r8=r6*r2 1219 1220 if ((.not.noabc).and.(r2.lt.abcthr)) then 1221 ! if (.not.noabc) then 1222 abccalc(linij)=.TRUE. 1223 dc6ij(iat,jat)=dc6iji 1224 dc6ij(jat,iat)=dc6ijj 1225 c6abc(linij)=c6 1226 r2abc(linij)=r2 1227 r3abc(linij)=(r/r0)**(1.0/3.0) 1228 !noabc 1229 end if 1230 ! use BJ radius 1231 R0=a1*dsqrt(3.0d0*r42)+a2 1232 1233 t6=(r6+R0**6) 1234 t8=(r8+R0**8) 1235 1236 drij(linij)=drij(linij) & 1237 & -s6*C6*6.0d0*r4*r/(t6*t6) & 1238 & -s8*C6*24.0d0*r42*r6*r/(t8*t8) 1239 1240 dc6_rest=s6/t6+3.d0*s8*r42/t8 1241 ! calculate E_disp for sanity check 1242 disp=disp-dc6_rest*c6 1243 1244 ! saving all (1/r^6...)* dC6/dCN(i) for each atom 1245 dc6i(iat)=dc6i(iat)+dc6_rest*dc6iji 1246 dc6i(jat)=dc6i(jat)+dc6_rest*dc6ijj 1247 1248 !jat 1249 end do 1250 !iat 1251 end do 1252 1253 !version=4 (BJ) 1254 end if 1255 1256 1257 if (.not.noabc)then 1258 1259 if (echo)write(*,*) 'doing analytical gradient O(N^3) ...' 1260 do iat=1,n 1261 do jat=1,iat-1 1262 linij=lin(iat,jat) 1263 if (.NOT.abccalc(linij))cycle 1264 r2ij=r2abc(linij) 1265 do kat=1,jat-1 1266 1267 linik=lin(iat,kat) 1268 linjk=lin(jat,kat) 1269 !cuto 1270 if (.NOT.(abccalc(linjk).AND.abccalc(linik)))cycle 1271 ! calculating the 3body energy: 1272 r2jk=r2abc(linjk) 1273 r2ik=r2abc(linik) 1274 c9=c6abc(linij)*c6abc(linjk)*c6abc(linik) 1275 c9=dsqrt(c9) 1276 rav=r3abc(linij)*r3abc(linjk)*r3abc(linik) 1277 fdmp=1.0d0/(1+6.0d0*(0.75d0*rav)**(-alp8)) 1278 mijk=-r2ij+r2jk+r2ik 1279 imjk= r2ij-r2jk+r2ik 1280 ijmk= r2ij+r2jk-r2ik 1281 rijk3=r2ij*r2jk*r2ik 1282 rav3=rijk3**1.5 1283 ang=0.375d0*ijmk*imjk*mijk/rijk3 1284 angr9=(ang +1.0d0) & 1285 & /rav3 1286 1287 eabc=eabc+c9*angr9*fdmp 1288 !end of 3body energy calculation 1289 1290 !start calculating the derivatives of each part w.r.t. r_ij 1291 1292 r=dsqrt(r2ij) 1293 dfdmp=-2.d0*alp8*(0.75d0*rav)**(-alp8)*fdmp*fdmp 1294 1295 dang=-0.375d0*(r2ij**3+r2ij**2*(r2jk+r2ik) & 1296 & +r2ij*(3.0d0*r2jk**2+2.0*r2jk*r2ik+3.0*r2ik**2) & 1297 & -5.0*(r2jk-r2ik)**2*(r2jk+r2ik)) & 1298 & /(r*rijk3*rav3) 1299 1300 1301 tmp1=dfdmp/r*c9*angr9-dang*c9*fdmp 1302 drij(linij)=drij(linij)+tmp1 1303 1304 !start calculating the derivatives of each part w.r.t. r_jk 1305 r=dsqrt(r2jk) 1306 1307 dang=-0.375d0*(r2jk**3+r2jk**2*(r2ik+r2ij) & 1308 & +r2jk*(3.0d0*r2ik**2+2.0*r2ik*r2ij+3.0*r2ij**2) & 1309 & -5.0*(r2ik-r2ij)**2*(r2ik+r2ij)) & 1310 & /(r*rijk3*rav3) 1311 1312 drij(linjk)=drij(linjk) & 1313 & +dfdmp/r*c9*angr9-dang*c9*fdmp 1314 1315 1316 !start calculating the derivatives of each part w.r.t. r_ik 1317 r=dsqrt(r2abc(linik)) 1318 1319 dang=-0.375d0*(r2ik**3+r2ik**2*(r2jk+r2ij) & 1320 & +r2ik*(3.0d0*r2jk**2+2.0*r2jk*r2ij+3.0*r2ij**2) & 1321 & -5.0*(r2jk-r2ij)**2*(r2jk+r2ij)) & 1322 & /(r*rijk3*rav3) 1323 1324 drij(linik)=drij(linik) & 1325 & +dfdmp/r*c9*angr9-dang*c9*fdmp 1326 1327 1328 1329 ! calculate rest* dc9/dcn(iat) and sum it up for every atom ijk 1330 dc6_rest=angr9*fdmp 1331 1332 dc9=dc6ij(iat,jat)/c6abc(linij)+dc6ij(iat,kat)/c6abc(linik) 1333 dc9=-0.5d0*c9*dc9 1334 dc6i(iat)=dc6i(iat)+dc6_rest*dc9 1335 1336 dc9=dc6ij(jat,iat)/c6abc(linij)+dc6ij(jat,kat)/c6abc(linjk) 1337 dc9=-0.5d0*c9*dc9 1338 dc6i(jat)=dc6i(jat)+dc6_rest*dc9 1339 1340 dc9=dc6ij(kat,iat)/c6abc(linik)+dc6ij(kat,jat)/c6abc(linjk) 1341 dc9=-0.5d0*c9*dc9 1342 dc6i(kat)=dc6i(kat)+dc6_rest*dc9 1343 1344 !kat 1345 end do 1346 !jat 1347 end do 1348 !iat 1349 end do 1350 1351 disp=disp+eabc 1352 !noabc 1353 end if 1354 1355 1356 1357 ! After calculating all derivatives dE/dr_ij w.r.t. distances, 1358 ! the grad w.r.t. the coordinates is calculated dE/dr_ij * dr_ij/dxyz_i 1359 do iat=2,n 1360 gtmp=0.0 1361 do jat=1,iat-1 1362 linij=lin(iat,jat) 1363 rij=xyz(:,jat)-xyz(:,iat) 1364 1365 1366 r2=sum(rij*rij) 1367 r=dsqrt(r2) 1368 if (r2.lt.cn_thr) then 1369 rcovij=rcov(iz(iat))+rcov(iz(jat)) 1370 expterm=exp(-k1*(rcovij/r-1.d0)) 1371 dcn=-k1*rcovij*expterm/ & 1372 & (r*r*(expterm+1.d0)*(expterm+1.d0)) 1373 else 1374 dcn=0.d0 1375 end if 1376 x1=drij(linij)+dcn*(dc6i(iat)+dc6i(jat)) 1377 1378 gtmp=gtmp+x1*rij/r 1379 !g(:,iat)=g(:,iat)+x1*rij/r 1380 g(:,jat)=g(:,jat)-x1*rij/r 1381 1382 !iat 1383 end do 1384 g(:,iat)=g(:,iat)+gtmp 1385 !jat 1386 end do 1387 1388 1389 1390999 continue 1391 gnorm=sum(abs(g(1:3,1:n))) 1392 if (echo)then 1393 write(*,*) 1394 write(*,*)'|G|=',gnorm 1395 end if 1396 1397 !fixed atoms have no gradient 1398 do i=1,n 1399 if (fix(i))g(:,i)=0.0 1400 end do 1401 1402 1403 1404 end subroutine gdisp 1405 1406 1407 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1408 ! The N E W gradC6 routine C 1409 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1410 ! 1411 subroutine get_dC6_dCNij(maxc,max_elem,c6ab,mxci,mxcj,cni,cnj, & 1412 & izi,izj,iat,jat,c6check,dc6i,dc6j) 1413 1414 integer maxc,max_elem 1415 real(wp) c6ab(max_elem,max_elem,maxc,maxc,3) 1416 !mxc(iz(iat)) 1417 integer mxci,mxcj 1418 real(wp) cni,cnj,term 1419 integer iat,jat,izi,izj 1420 real(wp) dc6i,dc6j,c6check 1421 1422 1423 integer a,b 1424 real(wp) zaehler,nenner,dzaehler_i,dnenner_i,dzaehler_j,dnenner_j 1425 real(wp) expterm,cn_refi,cn_refj,c6ref,r 1426 real(wp) c6mem,r_save 1427 1428 1429 1430 c6mem=-1.d99 1431 r_save=9999.0 1432 zaehler=0.0d0 1433 nenner=0.0d0 1434 1435 dzaehler_i=0.d0 1436 dnenner_i=0.d0 1437 dzaehler_j=0.d0 1438 dnenner_j=0.d0 1439 1440 1441 do a=1,mxci 1442 do b=1,mxcj 1443 c6ref=c6ab(izi,izj,a,b,1) 1444 if (c6ref.gt.0) then 1445 ! c6mem=c6ref 1446 cn_refi=c6ab(izi,izj,a,b,2) 1447 cn_refj=c6ab(izi,izj,a,b,3) 1448 r=(cn_refi-cni)*(cn_refi-cni)+(cn_refj-cnj)*(cn_refj-cnj) 1449 if (r.lt.r_save) then 1450 r_save=r 1451 c6mem=c6ref 1452 end if 1453 expterm=exp(k3*r) 1454 zaehler=zaehler+c6ref*expterm 1455 nenner=nenner+expterm 1456 expterm=expterm*2.d0*k3 1457 term=expterm*(cni-cn_refi) 1458 dzaehler_i=dzaehler_i+c6ref*term 1459 dnenner_i =dnenner_i + term 1460 1461 term=expterm*(cnj-cn_refj) 1462 dzaehler_j=dzaehler_j+c6ref*term 1463 dnenner_j =dnenner_j + term 1464 end if 1465 !b 1466 end do 1467 !a 1468 end do 1469 1470 if (nenner.gt.1.0d-99) then 1471 c6check=zaehler/nenner 1472 dc6i=((dzaehler_i*nenner)-(dnenner_i*zaehler)) & 1473 & /(nenner*nenner) 1474 dc6j=((dzaehler_j*nenner)-(dnenner_j*zaehler)) & 1475 & /(nenner*nenner) 1476 else 1477 c6check=c6mem 1478 dc6i=0.0d0 1479 dc6j=0.0d0 1480 end if 1481 end subroutine get_dC6_dCNij 1482 1483 1484 1485 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1486 ! interpolate c6 1487 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1488 1489 subroutine getc6(maxc,max_elem,c6ab,mxc,iat,jat,nci,ncj,c6) 1490 integer maxc,max_elem 1491 integer iat,jat,i,j,mxc(max_elem) 1492 real(wp) nci,ncj,c6,c6mem 1493 real(wp) c6ab(max_elem,max_elem,maxc,maxc,3) 1494 ! the exponential is sensitive to numerics 1495 ! when nci or ncj is much larger than cn1/cn2 1496 real(wp) cn1,cn2,r,rsum,csum,tmp1 1497 real(wp) r_save 1498 1499 c6mem=-1.d+99 1500 rsum=0.0d0 1501 csum=0.0d0 1502 c6 =0.0d0 1503 r_save=1.0d99 1504 do i=1,mxc(iat) 1505 do j=1,mxc(jat) 1506 c6=c6ab(iat,jat,i,j,1) 1507 if (c6.gt.0)then 1508 ! c6mem=c6 1509 cn1=c6ab(iat,jat,i,j,2) 1510 cn2=c6ab(iat,jat,i,j,3) 1511 ! distance 1512 r=(cn1-nci)**2+(cn2-ncj)**2 1513 if (r.lt.r_save) then 1514 r_save=r 1515 c6mem=c6 1516 end if 1517 tmp1=exp(k3*r) 1518 rsum=rsum+tmp1 1519 csum=csum+tmp1*c6 1520 end if 1521 end do 1522 end do 1523 1524 if (rsum.gt.1.0d-99)then 1525 c6=csum/rsum 1526 else 1527 c6=c6mem 1528 end if 1529 1530 end subroutine getc6 1531 1532 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1533 ! compute coordination numbers by adding an inverse damping function 1534 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1535 1536 subroutine ncoord(natoms,rcov,iz,xyz,cn,cn_thr) 1537 integer iz(*),natoms,i 1538 real(wp) xyz(3,*),cn(*),rcov(94) 1539 real(wp) cn_thr 1540 1541 integer iat 1542 real(wp) dx,dy,dz,r,damp,xn,rr,rco,r2 1543 1544 do i=1,natoms 1545 xn=0.0d0 1546 do iat=1,natoms 1547 if (iat.ne.i)then 1548 dx=xyz(1,iat)-xyz(1,i) 1549 dy=xyz(2,iat)-xyz(2,i) 1550 dz=xyz(3,iat)-xyz(3,i) 1551 r2=dx*dx+dy*dy+dz*dz 1552 if (r2.gt.cn_thr) cycle 1553 r=sqrt(r2) 1554 ! covalent distance in Bohr 1555 rco=rcov(iz(i))+rcov(iz(iat)) 1556 rr=rco/r 1557 ! counting function exponential has a better long-range behavior than MH 1558 damp=1.d0/(1.d0+exp(-k1*(rr-1.0d0))) 1559 xn=xn+damp 1560 end if 1561 end do 1562 ! if (iz(i).eq.19) then 1563 ! write(*,*) "Input CN of Kalium" 1564 ! read(*,*),input 1565 ! cn(i)=input 1566 ! else 1567 cn(i)=xn 1568 ! end if 1569 end do 1570 1571 end subroutine ncoord 1572 1573 1574 integer function lin(i1,i2) 1575 integer i1,i2,idum1,idum2 1576 idum1=max(i1,i2) 1577 idum2=min(i1,i2) 1578 lin=idum2+idum1*(idum1-1)/2 1579 return 1580 end function lin 1581 1582 1583 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1584 ! set cut-off radii 1585 ! in parts due to INTEL compiler bug 1586 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1587 1588 subroutine setr0ab(max_elem,autoang,r) 1589 integer max_elem,i,j,k 1590 real(wp) r(max_elem,max_elem),autoang 1591 real(wp) r0ab(4465) 1592 r0ab( 1: 70)=(/ & 1593 & 2.1823, 1.8547, 1.7347, 2.9086, 2.5732, 3.4956, 2.3550 & 1594 &, 2.5095, 2.9802, 3.0982, 2.5141, 2.3917, 2.9977, 2.9484 & 1595 &, 3.2160, 2.4492, 2.2527, 3.1933, 3.0214, 2.9531, 2.9103 & 1596 &, 2.3667, 2.1328, 2.8784, 2.7660, 2.7776, 2.7063, 2.6225 & 1597 &, 2.1768, 2.0625, 2.6395, 2.6648, 2.6482, 2.5697, 2.4846 & 1598 &, 2.4817, 2.0646, 1.9891, 2.5086, 2.6908, 2.6233, 2.4770 & 1599 &, 2.3885, 2.3511, 2.2996, 1.9892, 1.9251, 2.4190, 2.5473 & 1600 &, 2.4994, 2.4091, 2.3176, 2.2571, 2.1946, 2.1374, 2.9898 & 1601 &, 2.6397, 3.6031, 3.1219, 3.7620, 3.2485, 2.9357, 2.7093 & 1602 &, 2.5781, 2.4839, 3.7082, 2.5129, 2.7321, 3.1052, 3.2962 & 1603 &/) 1604 r0ab( 71: 140)=(/ & 1605 & 3.1331, 3.2000, 2.9586, 3.0822, 2.8582, 2.7120, 3.2570 & 1606 &, 3.4839, 2.8766, 2.7427, 3.2776, 3.2363, 3.5929, 3.2826 & 1607 &, 3.0911, 2.9369, 2.9030, 2.7789, 3.3921, 3.3970, 4.0106 & 1608 &, 2.8884, 2.6605, 3.7513, 3.1613, 3.3605, 3.3325, 3.0991 & 1609 &, 2.9297, 2.8674, 2.7571, 3.8129, 3.3266, 3.7105, 3.7917 & 1610 &, 2.8304, 2.5538, 3.3932, 3.1193, 3.1866, 3.1245, 3.0465 & 1611 &, 2.8727, 2.7664, 2.6926, 3.4608, 3.2984, 3.5142, 3.5418 & 1612 &, 3.5017, 2.6190, 2.4797, 3.1331, 3.0540, 3.0651, 2.9879 & 1613 &, 2.9054, 2.8805, 2.7330, 2.6331, 3.2096, 3.5668, 3.3684 & 1614 &, 3.3686, 3.3180, 3.3107, 2.4757, 2.4019, 2.9789, 3.1468 & 1615 &/) 1616 r0ab( 141: 210)=(/ & 1617 & 2.9768, 2.8848, 2.7952, 2.7457, 2.6881, 2.5728, 3.0574 & 1618 &, 3.3264, 3.3562, 3.2529, 3.1916, 3.1523, 3.1046, 2.3725 & 1619 &, 2.3289, 2.8760, 2.9804, 2.9093, 2.8040, 2.7071, 2.6386 & 1620 &, 2.5720, 2.5139, 2.9517, 3.1606, 3.2085, 3.1692, 3.0982 & 1621 &, 3.0352, 2.9730, 2.9148, 3.2147, 2.8315, 3.8724, 3.4621 & 1622 &, 3.8823, 3.3760, 3.0746, 2.8817, 2.7552, 2.6605, 3.9740 & 1623 &, 3.6192, 3.6569, 3.9586, 3.6188, 3.3917, 3.2479, 3.1434 & 1624 &, 4.2411, 2.7597, 3.0588, 3.3474, 3.6214, 3.4353, 3.4729 & 1625 &, 3.2487, 3.3200, 3.0914, 2.9403, 3.4972, 3.7993, 3.6773 & 1626 &, 3.8678, 3.5808, 3.8243, 3.5826, 3.4156, 3.8765, 4.1035 & 1627 &/) 1628 r0ab( 211: 280)=(/ & 1629 & 2.7361, 2.9765, 3.2475, 3.5004, 3.4185, 3.4378, 3.2084 & 1630 &, 3.2787, 3.0604, 2.9187, 3.4037, 3.6759, 3.6586, 3.8327 & 1631 &, 3.5372, 3.7665, 3.5310, 3.3700, 3.7788, 3.9804, 3.8903 & 1632 &, 2.6832, 2.9060, 3.2613, 3.4359, 3.3538, 3.3860, 3.1550 & 1633 &, 3.2300, 3.0133, 2.8736, 3.4024, 3.6142, 3.5979, 3.5295 & 1634 &, 3.4834, 3.7140, 3.4782, 3.3170, 3.7434, 3.9623, 3.8181 & 1635 &, 3.7642, 2.6379, 2.8494, 3.1840, 3.4225, 3.2771, 3.3401 & 1636 &, 3.1072, 3.1885, 2.9714, 2.8319, 3.3315, 3.5979, 3.5256 & 1637 &, 3.4980, 3.4376, 3.6714, 3.4346, 3.2723, 3.6859, 3.8985 & 1638 &, 3.7918, 3.7372, 3.7211, 2.9230, 2.6223, 3.4161, 2.8999 & 1639 &/) 1640 r0ab( 281: 350)=(/ & 1641 & 3.0557, 3.3308, 3.0555, 2.8508, 2.7385, 2.6640, 3.5263 & 1642 &, 3.0277, 3.2990, 3.7721, 3.5017, 3.2751, 3.1368, 3.0435 & 1643 &, 3.7873, 3.2858, 3.2140, 3.1727, 3.2178, 3.4414, 2.5490 & 1644 &, 2.7623, 3.0991, 3.3252, 3.1836, 3.2428, 3.0259, 3.1225 & 1645 &, 2.9032, 2.7621, 3.2490, 3.5110, 3.4429, 3.3845, 3.3574 & 1646 &, 3.6045, 3.3658, 3.2013, 3.6110, 3.8241, 3.7090, 3.6496 & 1647 &, 3.6333, 3.0896, 3.5462, 2.4926, 2.7136, 3.0693, 3.2699 & 1648 &, 3.1272, 3.1893, 2.9658, 3.0972, 2.8778, 2.7358, 3.2206 & 1649 &, 3.4566, 3.3896, 3.3257, 3.2946, 3.5693, 3.3312, 3.1670 & 1650 &, 3.5805, 3.7711, 3.6536, 3.5927, 3.5775, 3.0411, 3.4885 & 1651 &/) 1652 r0ab( 351: 420)=(/ & 1653 & 3.4421, 2.4667, 2.6709, 3.0575, 3.2357, 3.0908, 3.1537 & 1654 &, 2.9235, 3.0669, 2.8476, 2.7054, 3.2064, 3.4519, 3.3593 & 1655 &, 3.2921, 3.2577, 3.2161, 3.2982, 3.1339, 3.5606, 3.7582 & 1656 &, 3.6432, 3.5833, 3.5691, 3.0161, 3.4812, 3.4339, 3.4327 & 1657 &, 2.4515, 2.6338, 3.0511, 3.2229, 3.0630, 3.1265, 2.8909 & 1658 &, 3.0253, 2.8184, 2.6764, 3.1968, 3.4114, 3.3492, 3.2691 & 1659 &, 3.2320, 3.1786, 3.2680, 3.1036, 3.5453, 3.7259, 3.6090 & 1660 &, 3.5473, 3.5327, 3.0018, 3.4413, 3.3907, 3.3593, 3.3462 & 1661 &, 2.4413, 2.6006, 3.0540, 3.1987, 3.0490, 3.1058, 2.8643 & 1662 &, 2.9948, 2.7908, 2.6491, 3.1950, 3.3922, 3.3316, 3.2585 & 1663 &/) 1664 r0ab( 421: 490)=(/ & 1665 & 3.2136, 3.1516, 3.2364, 3.0752, 3.5368, 3.7117, 3.5941 & 1666 &, 3.5313, 3.5164, 2.9962, 3.4225, 3.3699, 3.3370, 3.3234 & 1667 &, 3.3008, 2.4318, 2.5729, 3.0416, 3.1639, 3.0196, 3.0843 & 1668 &, 2.8413, 2.7436, 2.7608, 2.6271, 3.1811, 3.3591, 3.3045 & 1669 &, 3.2349, 3.1942, 3.1291, 3.2111, 3.0534, 3.5189, 3.6809 & 1670 &, 3.5635, 3.5001, 3.4854, 2.9857, 3.3897, 3.3363, 3.3027 & 1671 &, 3.2890, 3.2655, 3.2309, 2.8502, 2.6934, 3.2467, 3.1921 & 1672 &, 3.5663, 3.2541, 3.0571, 2.9048, 2.8657, 2.7438, 3.3547 & 1673 &, 3.3510, 3.9837, 3.6871, 3.4862, 3.3389, 3.2413, 3.1708 & 1674 &, 3.6096, 3.6280, 3.6860, 3.5568, 3.4836, 3.2868, 3.3994 & 1675 &/) 1676 r0ab( 491: 560)=(/ & 1677 & 3.3476, 3.3170, 3.2950, 3.2874, 3.2606, 3.9579, 2.9226 & 1678 &, 2.6838, 3.7867, 3.1732, 3.3872, 3.3643, 3.1267, 2.9541 & 1679 &, 2.8505, 2.7781, 3.8475, 3.3336, 3.7359, 3.8266, 3.5733 & 1680 &, 3.3959, 3.2775, 3.1915, 3.9878, 3.8816, 3.5810, 3.5364 & 1681 &, 3.5060, 3.8097, 3.3925, 3.3348, 3.3019, 3.2796, 3.2662 & 1682 &, 3.2464, 3.7136, 3.8619, 2.9140, 2.6271, 3.4771, 3.1774 & 1683 &, 3.2560, 3.1970, 3.1207, 2.9406, 2.8322, 2.7571, 3.5455 & 1684 &, 3.3514, 3.5837, 3.6177, 3.5816, 3.3902, 3.2604, 3.1652 & 1685 &, 3.7037, 3.6283, 3.5858, 3.5330, 3.4884, 3.5789, 3.4094 & 1686 &, 3.3473, 3.3118, 3.2876, 3.2707, 3.2521, 3.5570, 3.6496 & 1687 &/) 1688 r0ab( 561: 630)=(/ & 1689 & 3.6625, 2.7300, 2.5870, 3.2471, 3.1487, 3.1667, 3.0914 & 1690 &, 3.0107, 2.9812, 2.8300, 2.7284, 3.3259, 3.3182, 3.4707 & 1691 &, 3.4748, 3.4279, 3.4182, 3.2547, 3.1353, 3.5116, 3.9432 & 1692 &, 3.8828, 3.8303, 3.7880, 3.3760, 3.7218, 3.3408, 3.3059 & 1693 &, 3.2698, 3.2446, 3.2229, 3.4422, 3.5023, 3.5009, 3.5268 & 1694 &, 2.6026, 2.5355, 3.1129, 3.2863, 3.1029, 3.0108, 2.9227 & 1695 &, 2.8694, 2.8109, 2.6929, 3.1958, 3.4670, 3.4018, 3.3805 & 1696 &, 3.3218, 3.2815, 3.2346, 3.0994, 3.3937, 3.7266, 3.6697 & 1697 &, 3.6164, 3.5730, 3.2522, 3.5051, 3.4686, 3.4355, 3.4084 & 1698 &, 3.3748, 3.3496, 3.3692, 3.4052, 3.3910, 3.3849, 3.3662 & 1699 &/) 1700 r0ab( 631: 700)=(/ & 1701 & 2.5087, 2.4814, 3.0239, 3.1312, 3.0535, 2.9457, 2.8496 & 1702 &, 2.7780, 2.7828, 2.6532, 3.1063, 3.3143, 3.3549, 3.3120 & 1703 &, 3.2421, 3.1787, 3.1176, 3.0613, 3.3082, 3.5755, 3.5222 & 1704 &, 3.4678, 3.4231, 3.1684, 3.3528, 3.3162, 3.2827, 3.2527 & 1705 &, 3.2308, 3.2029, 3.3173, 3.3343, 3.3092, 3.2795, 3.2452 & 1706 &, 3.2096, 3.2893, 2.8991, 4.0388, 3.6100, 3.9388, 3.4475 & 1707 &, 3.1590, 2.9812, 2.8586, 2.7683, 4.1428, 3.7911, 3.8225 & 1708 &, 4.0372, 3.7059, 3.4935, 3.3529, 3.2492, 4.4352, 4.0826 & 1709 &, 3.9733, 3.9254, 3.8646, 3.9315, 3.7837, 3.7465, 3.7211 & 1710 &, 3.7012, 3.6893, 3.6676, 3.7736, 4.0660, 3.7926, 3.6158 & 1711 &/) 1712 r0ab( 701: 770)=(/ & 1713 & 3.5017, 3.4166, 4.6176, 2.8786, 3.1658, 3.5823, 3.7689 & 1714 &, 3.5762, 3.5789, 3.3552, 3.4004, 3.1722, 3.0212, 3.7241 & 1715 &, 3.9604, 3.8500, 3.9844, 3.7035, 3.9161, 3.6751, 3.5075 & 1716 &, 4.1151, 4.2877, 4.1579, 4.1247, 4.0617, 3.4874, 3.9848 & 1717 &, 3.9280, 3.9079, 3.8751, 3.8604, 3.8277, 3.8002, 3.9981 & 1718 &, 3.7544, 4.0371, 3.8225, 3.6718, 4.3092, 4.4764, 2.8997 & 1719 &, 3.0953, 3.4524, 3.6107, 3.6062, 3.5783, 3.3463, 3.3855 & 1720 &, 3.1746, 3.0381, 3.6019, 3.7938, 3.8697, 3.9781, 3.6877 & 1721 &, 3.8736, 3.6451, 3.4890, 3.9858, 4.1179, 4.0430, 3.9563 & 1722 &, 3.9182, 3.4002, 3.8310, 3.7716, 3.7543, 3.7203, 3.7053 & 1723 &/) 1724 r0ab( 771: 840)=(/ & 1725 & 3.6742, 3.8318, 3.7631, 3.7392, 3.9892, 3.7832, 3.6406 & 1726 &, 4.1701, 4.3016, 4.2196, 2.8535, 3.0167, 3.3978, 3.5363 & 1727 &, 3.5393, 3.5301, 3.2960, 3.3352, 3.1287, 2.9967, 3.6659 & 1728 &, 3.7239, 3.8070, 3.7165, 3.6368, 3.8162, 3.5885, 3.4336 & 1729 &, 3.9829, 4.0529, 3.9584, 3.9025, 3.8607, 3.3673, 3.7658 & 1730 &, 3.7035, 3.6866, 3.6504, 3.6339, 3.6024, 3.7708, 3.7283 & 1731 &, 3.6896, 3.9315, 3.7250, 3.5819, 4.1457, 4.2280, 4.1130 & 1732 &, 4.0597, 3.0905, 2.7998, 3.6448, 3.0739, 3.2996, 3.5262 & 1733 &, 3.2559, 3.0518, 2.9394, 2.8658, 3.7514, 3.2295, 3.5643 & 1734 &, 3.7808, 3.6931, 3.4723, 3.3357, 3.2429, 4.0280, 3.5589 & 1735 &/) 1736 r0ab( 841: 910)=(/ & 1737 & 3.4636, 3.4994, 3.4309, 3.6177, 3.2946, 3.2376, 3.2050 & 1738 &, 3.1847, 3.1715, 3.1599, 3.5555, 3.8111, 3.7693, 3.5718 & 1739 &, 3.4498, 3.3662, 4.1608, 3.7417, 3.6536, 3.6154, 3.8596 & 1740 &, 3.0301, 2.7312, 3.5821, 3.0473, 3.2137, 3.4679, 3.1975 & 1741 &, 2.9969, 2.8847, 2.8110, 3.6931, 3.2076, 3.4943, 3.5956 & 1742 &, 3.6379, 3.4190, 3.2808, 3.1860, 3.9850, 3.5105, 3.4330 & 1743 &, 3.3797, 3.4155, 3.6033, 3.2737, 3.2145, 3.1807, 3.1596 & 1744 &, 3.1461, 3.1337, 3.4812, 3.6251, 3.7152, 3.5201, 3.3966 & 1745 &, 3.3107, 4.1128, 3.6899, 3.6082, 3.5604, 3.7834, 3.7543 & 1746 &, 2.9189, 2.6777, 3.4925, 2.9648, 3.1216, 3.2940, 3.0975 & 1747 &/) 1748 r0ab( 911: 980)=(/ & 1749 & 2.9757, 2.8493, 2.7638, 3.6085, 3.1214, 3.4006, 3.4793 & 1750 &, 3.5147, 3.3806, 3.2356, 3.1335, 3.9144, 3.4183, 3.3369 & 1751 &, 3.2803, 3.2679, 3.4871, 3.1714, 3.1521, 3.1101, 3.0843 & 1752 &, 3.0670, 3.0539, 3.3890, 3.5086, 3.5895, 3.4783, 3.3484 & 1753 &, 3.2559, 4.0422, 3.5967, 3.5113, 3.4576, 3.6594, 3.6313 & 1754 &, 3.5690, 2.8578, 2.6334, 3.4673, 2.9245, 3.0732, 3.2435 & 1755 &, 3.0338, 2.9462, 2.8143, 2.7240, 3.5832, 3.0789, 3.3617 & 1756 &, 3.4246, 3.4505, 3.3443, 3.1964, 3.0913, 3.8921, 3.3713 & 1757 &, 3.2873, 3.2281, 3.2165, 3.4386, 3.1164, 3.1220, 3.0761 & 1758 &, 3.0480, 3.0295, 3.0155, 3.3495, 3.4543, 3.5260, 3.4413 & 1759 &/) 1760 r0ab( 981:1050)=(/ & 1761 & 3.3085, 3.2134, 4.0170, 3.5464, 3.4587, 3.4006, 3.6027 & 1762 &, 3.5730, 3.4945, 3.4623, 2.8240, 2.5960, 3.4635, 2.9032 & 1763 &, 3.0431, 3.2115, 2.9892, 2.9148, 2.7801, 2.6873, 3.5776 & 1764 &, 3.0568, 3.3433, 3.3949, 3.4132, 3.3116, 3.1616, 3.0548 & 1765 &, 3.8859, 3.3719, 3.2917, 3.2345, 3.2274, 3.4171, 3.1293 & 1766 &, 3.0567, 3.0565, 3.0274, 3.0087, 2.9939, 3.3293, 3.4249 & 1767 &, 3.4902, 3.4091, 3.2744, 3.1776, 4.0078, 3.5374, 3.4537 & 1768 &, 3.3956, 3.5747, 3.5430, 3.4522, 3.4160, 3.3975, 2.8004 & 1769 &, 2.5621, 3.4617, 2.9154, 3.0203, 3.1875, 2.9548, 2.8038 & 1770 &, 2.7472, 2.6530, 3.5736, 3.0584, 3.3304, 3.3748, 3.3871 & 1771 &/) 1772 r0ab(1051:1120)=(/ & 1773 & 3.2028, 3.1296, 3.0214, 3.8796, 3.3337, 3.2492, 3.1883 & 1774 &, 3.1802, 3.4050, 3.0756, 3.0478, 3.0322, 3.0323, 3.0163 & 1775 &, 3.0019, 3.3145, 3.4050, 3.4656, 3.3021, 3.2433, 3.1453 & 1776 &, 3.9991, 3.5017, 3.4141, 3.3520, 3.5583, 3.5251, 3.4243 & 1777 &, 3.3851, 3.3662, 3.3525, 2.7846, 2.5324, 3.4652, 2.8759 & 1778 &, 3.0051, 3.1692, 2.9273, 2.7615, 2.7164, 2.6212, 3.5744 & 1779 &, 3.0275, 3.3249, 3.3627, 3.3686, 3.1669, 3.0584, 2.9915 & 1780 &, 3.8773, 3.3099, 3.2231, 3.1600, 3.1520, 3.4023, 3.0426 & 1781 &, 3.0099, 2.9920, 2.9809, 2.9800, 2.9646, 3.3068, 3.3930 & 1782 &, 3.4486, 3.2682, 3.1729, 3.1168, 3.9952, 3.4796, 3.3901 & 1783 &/) 1784 r0ab(1121:1190)=(/ & 1785 & 3.3255, 3.5530, 3.5183, 3.4097, 3.3683, 3.3492, 3.3360 & 1786 &, 3.3308, 2.5424, 2.6601, 3.2555, 3.2807, 3.1384, 3.1737 & 1787 &, 2.9397, 2.8429, 2.8492, 2.7225, 3.3875, 3.4910, 3.4520 & 1788 &, 3.3608, 3.3036, 3.2345, 3.2999, 3.1487, 3.7409, 3.8392 & 1789 &, 3.7148, 3.6439, 3.6182, 3.1753, 3.5210, 3.4639, 3.4265 & 1790 &, 3.4075, 3.3828, 3.3474, 3.4071, 3.3754, 3.3646, 3.3308 & 1791 &, 3.4393, 3.2993, 3.8768, 3.9891, 3.8310, 3.7483, 3.3417 & 1792 &, 3.3019, 3.2250, 3.1832, 3.1578, 3.1564, 3.1224, 3.4620 & 1793 &, 2.9743, 2.8058, 3.4830, 3.3474, 3.6863, 3.3617, 3.1608 & 1794 &, 3.0069, 2.9640, 2.8427, 3.5885, 3.5219, 4.1314, 3.8120 & 1795 &/) 1796 r0ab(1191:1260)=(/ & 1797 & 3.6015, 3.4502, 3.3498, 3.2777, 3.8635, 3.8232, 3.8486 & 1798 &, 3.7215, 3.6487, 3.4724, 3.5627, 3.5087, 3.4757, 3.4517 & 1799 &, 3.4423, 3.4139, 4.1028, 3.8388, 3.6745, 3.5562, 3.4806 & 1800 &, 3.4272, 4.0182, 3.9991, 4.0007, 3.9282, 3.7238, 3.6498 & 1801 &, 3.5605, 3.5211, 3.5009, 3.4859, 3.4785, 3.5621, 4.2623 & 1802 &, 3.0775, 2.8275, 4.0181, 3.3385, 3.5379, 3.5036, 3.2589 & 1803 &, 3.0804, 3.0094, 2.9003, 4.0869, 3.5088, 3.9105, 3.9833 & 1804 &, 3.7176, 3.5323, 3.4102, 3.3227, 4.2702, 4.0888, 3.7560 & 1805 &, 3.7687, 3.6681, 3.6405, 3.5569, 3.4990, 3.4659, 3.4433 & 1806 &, 3.4330, 3.4092, 3.8867, 4.0190, 3.7961, 3.6412, 3.5405 & 1807 &/) 1808 r0ab(1261:1330)=(/ & 1809 & 3.4681, 4.3538, 4.2136, 3.9381, 3.8912, 3.9681, 3.7909 & 1810 &, 3.6774, 3.6262, 3.5999, 3.5823, 3.5727, 3.5419, 4.0245 & 1811 &, 4.1874, 3.0893, 2.7917, 3.7262, 3.3518, 3.4241, 3.5433 & 1812 &, 3.2773, 3.0890, 2.9775, 2.9010, 3.8048, 3.5362, 3.7746 & 1813 &, 3.7911, 3.7511, 3.5495, 3.4149, 3.3177, 4.0129, 3.8370 & 1814 &, 3.7739, 3.7125, 3.7152, 3.7701, 3.5813, 3.5187, 3.4835 & 1815 &, 3.4595, 3.4439, 3.4242, 3.7476, 3.8239, 3.8346, 3.6627 & 1816 &, 3.5479, 3.4639, 4.1026, 3.9733, 3.9292, 3.8667, 3.9513 & 1817 &, 3.8959, 3.7698, 3.7089, 3.6765, 3.6548, 3.6409, 3.5398 & 1818 &, 3.8759, 3.9804, 4.0150, 2.9091, 2.7638, 3.5066, 3.3377 & 1819 &/) 1820 r0ab(1331:1400)=(/ & 1821 & 3.3481, 3.2633, 3.1810, 3.1428, 2.9872, 2.8837, 3.5929 & 1822 &, 3.5183, 3.6729, 3.6596, 3.6082, 3.5927, 3.4224, 3.2997 & 1823 &, 3.8190, 4.1865, 4.1114, 4.0540, 3.6325, 3.5697, 3.5561 & 1824 &, 3.5259, 3.4901, 3.4552, 3.4315, 3.4091, 3.6438, 3.6879 & 1825 &, 3.6832, 3.7043, 3.5557, 3.4466, 3.9203, 4.2919, 4.2196 & 1826 &, 4.1542, 3.7573, 3.7039, 3.6546, 3.6151, 3.5293, 3.4849 & 1827 &, 3.4552, 3.5192, 3.7673, 3.8359, 3.8525, 3.8901, 2.7806 & 1828 &, 2.7209, 3.3812, 3.4958, 3.2913, 3.1888, 3.0990, 3.0394 & 1829 &, 2.9789, 2.8582, 3.4716, 3.6883, 3.6105, 3.5704, 3.5059 & 1830 &, 3.4619, 3.4138, 3.2742, 3.7080, 3.9773, 3.9010, 3.8409 & 1831 &/) 1832 r0ab(1401:1470)=(/ & 1833 & 3.7944, 3.4465, 3.7235, 3.6808, 3.6453, 3.6168, 3.5844 & 1834 &, 3.5576, 3.5772, 3.5959, 3.5768, 3.5678, 3.5486, 3.4228 & 1835 &, 3.8107, 4.0866, 4.0169, 3.9476, 3.6358, 3.5800, 3.5260 & 1836 &, 3.4838, 3.4501, 3.4204, 3.3553, 3.6487, 3.6973, 3.7398 & 1837 &, 3.7405, 3.7459, 3.7380, 2.6848, 2.6740, 3.2925, 3.3386 & 1838 &, 3.2473, 3.1284, 3.0301, 2.9531, 2.9602, 2.8272, 3.3830 & 1839 &, 3.5358, 3.5672, 3.5049, 3.4284, 3.3621, 3.3001, 3.2451 & 1840 &, 3.6209, 3.8299, 3.7543, 3.6920, 3.6436, 3.3598, 3.5701 & 1841 &, 3.5266, 3.4904, 3.4590, 3.4364, 3.4077, 3.5287, 3.5280 & 1842 &, 3.4969, 3.4650, 3.4304, 3.3963, 3.7229, 3.9402, 3.8753 & 1843 &/) 1844 r0ab(1471:1540)=(/ & 1845 & 3.8035, 3.5499, 3.4913, 3.4319, 3.3873, 3.3520, 3.3209 & 1846 &, 3.2948, 3.5052, 3.6465, 3.6696, 3.6577, 3.6388, 3.6142 & 1847 &, 3.5889, 3.3968, 3.0122, 4.2241, 3.7887, 4.0049, 3.5384 & 1848 &, 3.2698, 3.1083, 2.9917, 2.9057, 4.3340, 3.9900, 4.6588 & 1849 &, 4.1278, 3.8125, 3.6189, 3.4851, 3.3859, 4.6531, 4.3134 & 1850 &, 4.2258, 4.1309, 4.0692, 4.0944, 3.9850, 3.9416, 3.9112 & 1851 &, 3.8873, 3.8736, 3.8473, 4.6027, 4.1538, 3.8994, 3.7419 & 1852 &, 3.6356, 3.5548, 4.8353, 4.5413, 4.3891, 4.3416, 4.3243 & 1853 &, 4.2753, 4.2053, 4.1790, 4.1685, 4.1585, 4.1536, 4.0579 & 1854 &, 4.1980, 4.4564, 4.2192, 4.0528, 3.9489, 3.8642, 5.0567 & 1855 &/) 1856 r0ab(1541:1610)=(/ & 1857 & 3.0630, 3.3271, 4.0432, 4.0046, 4.1555, 3.7426, 3.5130 & 1858 &, 3.5174, 3.2884, 3.1378, 4.1894, 4.2321, 4.1725, 4.1833 & 1859 &, 3.8929, 4.0544, 3.8118, 3.6414, 4.6373, 4.6268, 4.4750 & 1860 &, 4.4134, 4.3458, 3.8582, 4.2583, 4.1898, 4.1562, 4.1191 & 1861 &, 4.1069, 4.0639, 4.1257, 4.1974, 3.9532, 4.1794, 3.9660 & 1862 &, 3.8130, 4.8160, 4.8272, 4.6294, 4.5840, 4.0770, 4.0088 & 1863 &, 3.9103, 3.8536, 3.8324, 3.7995, 3.7826, 4.2294, 4.3380 & 1864 &, 4.4352, 4.1933, 4.4580, 4.2554, 4.1072, 5.0454, 5.1814 & 1865 &, 3.0632, 3.2662, 3.6432, 3.8088, 3.7910, 3.7381, 3.5093 & 1866 &, 3.5155, 3.3047, 3.1681, 3.7871, 3.9924, 4.0637, 4.1382 & 1867 &/) 1868 r0ab(1611:1680)=(/ & 1869 & 3.8591, 4.0164, 3.7878, 3.6316, 4.1741, 4.3166, 4.2395 & 1870 &, 4.1831, 4.1107, 3.5857, 4.0270, 3.9676, 3.9463, 3.9150 & 1871 &, 3.9021, 3.8708, 4.0240, 4.1551, 3.9108, 4.1337, 3.9289 & 1872 &, 3.7873, 4.3666, 4.5080, 4.4232, 4.3155, 3.8461, 3.8007 & 1873 &, 3.6991, 3.6447, 3.6308, 3.5959, 3.5749, 4.0359, 4.3124 & 1874 &, 4.3539, 4.1122, 4.3772, 4.1785, 4.0386, 4.7004, 4.8604 & 1875 &, 4.6261, 2.9455, 3.2470, 3.6108, 3.8522, 3.6625, 3.6598 & 1876 &, 3.4411, 3.4660, 3.2415, 3.0944, 3.7514, 4.0397, 3.9231 & 1877 &, 4.0561, 3.7860, 3.9845, 3.7454, 3.5802, 4.1366, 4.3581 & 1878 &, 4.2351, 4.2011, 4.1402, 3.5381, 4.0653, 4.0093, 3.9883 & 1879 &/) 1880 r0ab(1681:1750)=(/ & 1881 & 3.9570, 3.9429, 3.9112, 3.8728, 4.0682, 3.8351, 4.1054 & 1882 &, 3.8928, 3.7445, 4.3415, 4.5497, 4.3833, 4.3122, 3.8051 & 1883 &, 3.7583, 3.6622, 3.6108, 3.5971, 3.5628, 3.5408, 4.0780 & 1884 &, 4.0727, 4.2836, 4.0553, 4.3647, 4.1622, 4.0178, 4.5802 & 1885 &, 4.9125, 4.5861, 4.6201, 2.9244, 3.2241, 3.5848, 3.8293 & 1886 &, 3.6395, 3.6400, 3.4204, 3.4499, 3.2253, 3.0779, 3.7257 & 1887 &, 4.0170, 3.9003, 4.0372, 3.7653, 3.9672, 3.7283, 3.5630 & 1888 &, 4.1092, 4.3347, 4.2117, 4.1793, 4.1179, 3.5139, 4.0426 & 1889 &, 3.9867, 3.9661, 3.9345, 3.9200, 3.8883, 3.8498, 4.0496 & 1890 &, 3.8145, 4.0881, 3.8756, 3.7271, 4.3128, 4.5242, 4.3578 & 1891 &/) 1892 r0ab(1751:1820)=(/ & 1893 & 4.2870, 3.7796, 3.7318, 3.6364, 3.5854, 3.5726, 3.5378 & 1894 &, 3.5155, 4.0527, 4.0478, 4.2630, 4.0322, 4.3449, 4.1421 & 1895 &, 3.9975, 4.5499, 4.8825, 4.5601, 4.5950, 4.5702, 2.9046 & 1896 &, 3.2044, 3.5621, 3.8078, 3.6185, 3.6220, 3.4019, 3.4359 & 1897 &, 3.2110, 3.0635, 3.7037, 3.9958, 3.8792, 4.0194, 3.7460 & 1898 &, 3.9517, 3.7128, 3.5474, 4.0872, 4.3138, 4.1906, 4.1593 & 1899 &, 4.0973, 3.4919, 4.0216, 3.9657, 3.9454, 3.9134, 3.8986 & 1900 &, 3.8669, 3.8289, 4.0323, 3.7954, 4.0725, 3.8598, 3.7113 & 1901 &, 4.2896, 4.5021, 4.3325, 4.2645, 3.7571, 3.7083, 3.6136 & 1902 &, 3.5628, 3.5507, 3.5155, 3.4929, 4.0297, 4.0234, 4.2442 & 1903 &/) 1904 r0ab(1821:1890)=(/ & 1905 & 4.0112, 4.3274, 4.1240, 3.9793, 4.5257, 4.8568, 4.5353 & 1906 &, 4.5733, 4.5485, 4.5271, 2.8878, 3.1890, 3.5412, 3.7908 & 1907 &, 3.5974, 3.6078, 3.3871, 3.4243, 3.1992, 3.0513, 3.6831 & 1908 &, 3.9784, 3.8579, 4.0049, 3.7304, 3.9392, 3.7002, 3.5347 & 1909 &, 4.0657, 4.2955, 4.1705, 4.1424, 4.0800, 3.4717, 4.0043 & 1910 &, 3.9485, 3.9286, 3.8965, 3.8815, 3.8500, 3.8073, 4.0180 & 1911 &, 3.7796, 4.0598, 3.8470, 3.6983, 4.2678, 4.4830, 4.3132 & 1912 &, 4.2444, 3.7370, 3.6876, 3.5935, 3.5428, 3.5314, 3.4958 & 1913 &, 3.4730, 4.0117, 4.0043, 4.2287, 3.9939, 4.3134, 4.1096 & 1914 &, 3.9646, 4.5032, 4.8356, 4.5156, 4.5544, 4.5297, 4.5083 & 1915 &/) 1916 r0ab(1891:1960)=(/ & 1917 & 4.4896, 2.8709, 3.1737, 3.5199, 3.7734, 3.5802, 3.5934 & 1918 &, 3.3724, 3.4128, 3.1877, 3.0396, 3.6624, 3.9608, 3.8397 & 1919 &, 3.9893, 3.7145, 3.9266, 3.6877, 3.5222, 4.0448, 4.2771 & 1920 &, 4.1523, 4.1247, 4.0626, 3.4530, 3.9866, 3.9310, 3.9115 & 1921 &, 3.8792, 3.8641, 3.8326, 3.7892, 4.0025, 3.7636, 4.0471 & 1922 &, 3.8343, 3.6854, 4.2464, 4.4635, 4.2939, 4.2252, 3.7169 & 1923 &, 3.6675, 3.5739, 3.5235, 3.5126, 3.4768, 3.4537, 3.9932 & 1924 &, 3.9854, 4.2123, 3.9765, 4.2992, 4.0951, 3.9500, 4.4811 & 1925 &, 4.8135, 4.4959, 4.5351, 4.5105, 4.4891, 4.4705, 4.4515 & 1926 &, 2.8568, 3.1608, 3.5050, 3.7598, 3.5665, 3.5803, 3.3601 & 1927 &/) 1928 r0ab(1961:2030)=(/ & 1929 & 3.4031, 3.1779, 3.0296, 3.6479, 3.9471, 3.8262, 3.9773 & 1930 &, 3.7015, 3.9162, 3.6771, 3.5115, 4.0306, 4.2634, 4.1385 & 1931 &, 4.1116, 4.0489, 3.4366, 3.9732, 3.9176, 3.8983, 3.8659 & 1932 &, 3.8507, 3.8191, 3.7757, 3.9907, 3.7506, 4.0365, 3.8235 & 1933 &, 3.6745, 4.2314, 4.4490, 4.2792, 4.2105, 3.7003, 3.6510 & 1934 &, 3.5578, 3.5075, 3.4971, 3.4609, 3.4377, 3.9788, 3.9712 & 1935 &, 4.1997, 3.9624, 4.2877, 4.0831, 3.9378, 4.4655, 4.7974 & 1936 &, 4.4813, 4.5209, 4.4964, 4.4750, 4.4565, 4.4375, 4.4234 & 1937 &, 2.6798, 3.0151, 3.2586, 3.5292, 3.5391, 3.4902, 3.2887 & 1938 &, 3.3322, 3.1228, 2.9888, 3.4012, 3.7145, 3.7830, 3.6665 & 1939 &/) 1940 r0ab(2031:2100)=(/ & 1941 & 3.5898, 3.8077, 3.5810, 3.4265, 3.7726, 4.0307, 3.9763 & 1942 &, 3.8890, 3.8489, 3.2706, 3.7595, 3.6984, 3.6772, 3.6428 & 1943 &, 3.6243, 3.5951, 3.7497, 3.6775, 3.6364, 3.9203, 3.7157 & 1944 &, 3.5746, 3.9494, 4.2076, 4.1563, 4.0508, 3.5329, 3.4780 & 1945 &, 3.3731, 3.3126, 3.2846, 3.2426, 3.2135, 3.7491, 3.9006 & 1946 &, 3.8332, 3.8029, 4.1436, 3.9407, 3.7998, 4.1663, 4.5309 & 1947 &, 4.3481, 4.2911, 4.2671, 4.2415, 4.2230, 4.2047, 4.1908 & 1948 &, 4.1243, 2.5189, 2.9703, 3.3063, 3.6235, 3.4517, 3.3989 & 1949 &, 3.2107, 3.2434, 3.0094, 2.8580, 3.4253, 3.8157, 3.7258 & 1950 &, 3.6132, 3.5297, 3.7566, 3.5095, 3.3368, 3.7890, 4.1298 & 1951 &/) 1952 r0ab(2101:2170)=(/ & 1953 & 4.0190, 3.9573, 3.9237, 3.2677, 3.8480, 3.8157, 3.7656 & 1954 &, 3.7317, 3.7126, 3.6814, 3.6793, 3.6218, 3.5788, 3.8763 & 1955 &, 3.6572, 3.5022, 3.9737, 4.3255, 4.1828, 4.1158, 3.5078 & 1956 &, 3.4595, 3.3600, 3.3088, 3.2575, 3.2164, 3.1856, 3.8522 & 1957 &, 3.8665, 3.8075, 3.7772, 4.1391, 3.9296, 3.7772, 4.2134 & 1958 &, 4.7308, 4.3787, 4.3894, 4.3649, 4.3441, 4.3257, 4.3073 & 1959 &, 4.2941, 4.1252, 4.2427, 3.0481, 2.9584, 3.6919, 3.5990 & 1960 &, 3.8881, 3.4209, 3.1606, 3.1938, 2.9975, 2.8646, 3.8138 & 1961 &, 3.7935, 3.7081, 3.9155, 3.5910, 3.4808, 3.4886, 3.3397 & 1962 &, 4.1336, 4.1122, 3.9888, 3.9543, 3.8917, 3.5894, 3.8131 & 1963 &/) 1964 r0ab(2171:2240)=(/ & 1965 & 3.7635, 3.7419, 3.7071, 3.6880, 3.6574, 3.6546, 3.9375 & 1966 &, 3.6579, 3.5870, 3.6361, 3.5039, 4.3149, 4.2978, 4.1321 & 1967 &, 4.1298, 3.8164, 3.7680, 3.7154, 3.6858, 3.6709, 3.6666 & 1968 &, 3.6517, 3.8174, 3.8608, 4.1805, 3.9102, 3.8394, 3.8968 & 1969 &, 3.7673, 4.5274, 4.6682, 4.3344, 4.3639, 4.3384, 4.3162 & 1970 &, 4.2972, 4.2779, 4.2636, 4.0253, 4.1168, 4.1541, 2.8136 & 1971 &, 3.0951, 3.4635, 3.6875, 3.4987, 3.5183, 3.2937, 3.3580 & 1972 &, 3.1325, 2.9832, 3.6078, 3.8757, 3.7616, 3.9222, 3.6370 & 1973 &, 3.8647, 3.6256, 3.4595, 3.9874, 4.1938, 4.0679, 4.0430 & 1974 &, 3.9781, 3.3886, 3.9008, 3.8463, 3.8288, 3.7950, 3.7790 & 1975 &/) 1976 r0ab(2241:2310)=(/ & 1977 & 3.7472, 3.7117, 3.9371, 3.6873, 3.9846, 3.7709, 3.6210 & 1978 &, 4.1812, 4.3750, 4.2044, 4.1340, 3.6459, 3.5929, 3.5036 & 1979 &, 3.4577, 3.4528, 3.4146, 3.3904, 3.9014, 3.9031, 4.1443 & 1980 &, 3.8961, 4.2295, 4.0227, 3.8763, 4.4086, 4.7097, 4.4064 & 1981 &, 4.4488, 4.4243, 4.4029, 4.3842, 4.3655, 4.3514, 4.1162 & 1982 &, 4.2205, 4.1953, 4.2794, 2.8032, 3.0805, 3.4519, 3.6700 & 1983 &, 3.4827, 3.5050, 3.2799, 3.3482, 3.1233, 2.9747, 3.5971 & 1984 &, 3.8586, 3.7461, 3.9100, 3.6228, 3.8535, 3.6147, 3.4490 & 1985 &, 3.9764, 4.1773, 4.0511, 4.0270, 3.9614, 3.3754, 3.8836 & 1986 &, 3.8291, 3.8121, 3.7780, 3.7619, 3.7300, 3.6965, 3.9253 & 1987 &/) 1988 r0ab(2311:2380)=(/ & 1989 & 3.6734, 3.9733, 3.7597, 3.6099, 4.1683, 4.3572, 4.1862 & 1990 &, 4.1153, 3.6312, 3.5772, 3.4881, 3.4429, 3.4395, 3.4009 & 1991 &, 3.3766, 3.8827, 3.8868, 4.1316, 3.8807, 4.2164, 4.0092 & 1992 &, 3.8627, 4.3936, 4.6871, 4.3882, 4.4316, 4.4073, 4.3858 & 1993 &, 4.3672, 4.3485, 4.3344, 4.0984, 4.2036, 4.1791, 4.2622 & 1994 &, 4.2450, 2.7967, 3.0689, 3.4445, 3.6581, 3.4717, 3.4951 & 1995 &, 3.2694, 3.3397, 3.1147, 2.9661, 3.5898, 3.8468, 3.7358 & 1996 &, 3.9014, 3.6129, 3.8443, 3.6054, 3.4396, 3.9683, 4.1656 & 1997 &, 4.0394, 4.0158, 3.9498, 3.3677, 3.8718, 3.8164, 3.8005 & 1998 &, 3.7662, 3.7500, 3.7181, 3.6863, 3.9170, 3.6637, 3.9641 & 1999 &/) 2000 r0ab(2381:2450)=(/ & 2001 & 3.7503, 3.6004, 4.1590, 4.3448, 4.1739, 4.1029, 3.6224 & 2002 &, 3.5677, 3.4785, 3.4314, 3.4313, 3.3923, 3.3680, 3.8698 & 2003 &, 3.8758, 4.1229, 3.8704, 4.2063, 3.9987, 3.8519, 4.3832 & 2004 &, 4.6728, 4.3759, 4.4195, 4.3952, 4.3737, 4.3551, 4.3364 & 2005 &, 4.3223, 4.0861, 4.1911, 4.1676, 4.2501, 4.2329, 4.2208 & 2006 &, 2.7897, 3.0636, 3.4344, 3.6480, 3.4626, 3.4892, 3.2626 & 2007 &, 3.3344, 3.1088, 2.9597, 3.5804, 3.8359, 3.7251, 3.8940 & 2008 &, 3.6047, 3.8375, 3.5990, 3.4329, 3.9597, 4.1542, 4.0278 & 2009 &, 4.0048, 3.9390, 3.3571, 3.8608, 3.8056, 3.7899, 3.7560 & 2010 &, 3.7400, 3.7081, 3.6758, 3.9095, 3.6552, 3.9572, 3.7436 & 2011 &/) 2012 r0ab(2451:2520)=(/ & 2013 & 3.5933, 4.1508, 4.3337, 4.1624, 4.0916, 3.6126, 3.5582 & 2014 &, 3.4684, 3.4212, 3.4207, 3.3829, 3.3586, 3.8604, 3.8658 & 2015 &, 4.1156, 3.8620, 4.1994, 3.9917, 3.8446, 4.3750, 4.6617 & 2016 &, 4.3644, 4.4083, 4.3840, 4.3625, 4.3439, 4.3253, 4.3112 & 2017 &, 4.0745, 4.1807, 4.1578, 4.2390, 4.2218, 4.2097, 4.1986 & 2018 &, 2.8395, 3.0081, 3.3171, 3.4878, 3.5360, 3.5145, 3.2809 & 2019 &, 3.3307, 3.1260, 2.9940, 3.4741, 3.6675, 3.7832, 3.6787 & 2020 &, 3.6156, 3.8041, 3.5813, 3.4301, 3.8480, 3.9849, 3.9314 & 2021 &, 3.8405, 3.8029, 3.2962, 3.7104, 3.6515, 3.6378, 3.6020 & 2022 &, 3.5849, 3.5550, 3.7494, 3.6893, 3.6666, 3.9170, 3.7150 & 2023 &/) 2024 r0ab(2521:2590)=(/ & 2025 & 3.5760, 4.0268, 4.1596, 4.1107, 3.9995, 3.5574, 3.5103 & 2026 &, 3.4163, 3.3655, 3.3677, 3.3243, 3.2975, 3.7071, 3.9047 & 2027 &, 3.8514, 3.8422, 3.8022, 3.9323, 3.7932, 4.2343, 4.4583 & 2028 &, 4.3115, 4.2457, 4.2213, 4.1945, 4.1756, 4.1569, 4.1424 & 2029 &, 4.0620, 4.0494, 3.9953, 4.0694, 4.0516, 4.0396, 4.0280 & 2030 &, 4.0130, 2.9007, 2.9674, 3.8174, 3.5856, 3.6486, 3.5339 & 2031 &, 3.2832, 3.3154, 3.1144, 2.9866, 3.9618, 3.8430, 3.9980 & 2032 &, 3.8134, 3.6652, 3.7985, 3.5756, 3.4207, 4.4061, 4.2817 & 2033 &, 4.1477, 4.0616, 3.9979, 3.6492, 3.8833, 3.8027, 3.7660 & 2034 &, 3.7183, 3.6954, 3.6525, 3.9669, 3.8371, 3.7325, 3.9160 & 2035 &/) 2036 r0ab(2591:2660)=(/ & 2037 & 3.7156, 3.5714, 4.6036, 4.4620, 4.3092, 4.2122, 3.8478 & 2038 &, 3.7572, 3.6597, 3.5969, 3.5575, 3.5386, 3.5153, 3.7818 & 2039 &, 4.1335, 4.0153, 3.9177, 3.8603, 3.9365, 3.7906, 4.7936 & 2040 &, 4.7410, 4.5461, 4.5662, 4.5340, 4.5059, 4.4832, 4.4604 & 2041 &, 4.4429, 4.2346, 4.4204, 4.3119, 4.3450, 4.3193, 4.3035 & 2042 &, 4.2933, 4.1582, 4.2450, 2.8559, 2.9050, 3.8325, 3.5442 & 2043 &, 3.5077, 3.4905, 3.2396, 3.2720, 3.0726, 2.9467, 3.9644 & 2044 &, 3.8050, 3.8981, 3.7762, 3.6216, 3.7531, 3.5297, 3.3742 & 2045 &, 4.3814, 4.2818, 4.1026, 4.0294, 3.9640, 3.6208, 3.8464 & 2046 &, 3.7648, 3.7281, 3.6790, 3.6542, 3.6117, 3.8650, 3.8010 & 2047 &/) 2048 r0ab(2661:2730)=(/ & 2049 & 3.6894, 3.8713, 3.6699, 3.5244, 4.5151, 4.4517, 4.2538 & 2050 &, 4.1483, 3.8641, 3.7244, 3.6243, 3.5589, 3.5172, 3.4973 & 2051 &, 3.4715, 3.7340, 4.0316, 3.9958, 3.8687, 3.8115, 3.8862 & 2052 &, 3.7379, 4.7091, 4.7156, 4.5199, 4.5542, 4.5230, 4.4959 & 2053 &, 4.4750, 4.4529, 4.4361, 4.1774, 4.3774, 4.2963, 4.3406 & 2054 &, 4.3159, 4.3006, 4.2910, 4.1008, 4.1568, 4.0980, 2.8110 & 2055 &, 2.8520, 3.7480, 3.5105, 3.4346, 3.3461, 3.1971, 3.2326 & 2056 &, 3.0329, 2.9070, 3.8823, 3.7928, 3.8264, 3.7006, 3.5797 & 2057 &, 3.7141, 3.4894, 3.3326, 4.3048, 4.2217, 4.0786, 3.9900 & 2058 &, 3.9357, 3.6331, 3.8333, 3.7317, 3.6957, 3.6460, 3.6197 & 2059 &/) 2060 r0ab(2731:2800)=(/ & 2061 & 3.5779, 3.7909, 3.7257, 3.6476, 3.5729, 3.6304, 3.4834 & 2062 &, 4.4368, 4.3921, 4.2207, 4.1133, 3.8067, 3.7421, 3.6140 & 2063 &, 3.5491, 3.5077, 3.4887, 3.4623, 3.6956, 3.9568, 3.8976 & 2064 &, 3.8240, 3.7684, 3.8451, 3.6949, 4.6318, 4.6559, 4.4533 & 2065 &, 4.4956, 4.4641, 4.4366, 4.4155, 4.3936, 4.3764, 4.1302 & 2066 &, 4.3398, 4.2283, 4.2796, 4.2547, 4.2391, 4.2296, 4.0699 & 2067 &, 4.1083, 4.0319, 3.9855, 2.7676, 2.8078, 3.6725, 3.4804 & 2068 &, 3.3775, 3.2411, 3.1581, 3.1983, 2.9973, 2.8705, 3.8070 & 2069 &, 3.7392, 3.7668, 3.6263, 3.5402, 3.6807, 3.4545, 3.2962 & 2070 &, 4.2283, 4.1698, 4.0240, 3.9341, 3.8711, 3.5489, 3.7798 & 2071 &/) 2072 r0ab(2801:2870)=(/ & 2073 & 3.7000, 3.6654, 3.6154, 3.5882, 3.5472, 3.7289, 3.6510 & 2074 &, 3.6078, 3.5355, 3.5963, 3.4480, 4.3587, 4.3390, 4.1635 & 2075 &, 4.0536, 3.7193, 3.6529, 3.5512, 3.4837, 3.4400, 3.4191 & 2076 &, 3.3891, 3.6622, 3.8934, 3.8235, 3.7823, 3.7292, 3.8106 & 2077 &, 3.6589, 4.5535, 4.6013, 4.3961, 4.4423, 4.4109, 4.3835 & 2078 &, 4.3625, 4.3407, 4.3237, 4.0863, 4.2835, 4.1675, 4.2272 & 2079 &, 4.2025, 4.1869, 4.1774, 4.0126, 4.0460, 3.9815, 3.9340 & 2080 &, 3.8955, 2.6912, 2.7604, 3.6037, 3.4194, 3.3094, 3.1710 & 2081 &, 3.0862, 3.1789, 2.9738, 2.8427, 3.7378, 3.6742, 3.6928 & 2082 &, 3.5512, 3.4614, 3.4087, 3.4201, 3.2607, 4.1527, 4.0977 & 2083 &/) 2084 r0ab(2871:2940)=(/ & 2085 & 3.9523, 3.8628, 3.8002, 3.4759, 3.7102, 3.6466, 3.6106 & 2086 &, 3.5580, 3.5282, 3.4878, 3.6547, 3.5763, 3.5289, 3.5086 & 2087 &, 3.5593, 3.4099, 4.2788, 4.2624, 4.0873, 3.9770, 3.6407 & 2088 &, 3.5743, 3.5178, 3.4753, 3.3931, 3.3694, 3.3339, 3.6002 & 2089 &, 3.8164, 3.7478, 3.7028, 3.6952, 3.7669, 3.6137, 4.4698 & 2090 &, 4.5488, 4.3168, 4.3646, 4.3338, 4.3067, 4.2860, 4.2645 & 2091 &, 4.2478, 4.0067, 4.2349, 4.0958, 4.1543, 4.1302, 4.1141 & 2092 &, 4.1048, 3.9410, 3.9595, 3.8941, 3.8465, 3.8089, 3.7490 & 2093 &, 2.7895, 2.5849, 3.6484, 3.0162, 3.1267, 3.2125, 3.0043 & 2094 &, 2.9572, 2.8197, 2.7261, 3.7701, 3.2446, 3.5239, 3.4696 & 2095 &/) 2096 r0ab(2941:3010)=(/ & 2097 & 3.4261, 3.3508, 3.1968, 3.0848, 4.1496, 3.6598, 3.5111 & 2098 &, 3.4199, 3.3809, 3.5382, 3.2572, 3.2100, 3.1917, 3.1519 & 2099 &, 3.1198, 3.1005, 3.5071, 3.5086, 3.5073, 3.4509, 3.3120 & 2100 &, 3.2082, 4.2611, 3.8117, 3.6988, 3.5646, 3.6925, 3.6295 & 2101 &, 3.5383, 3.4910, 3.4625, 3.4233, 3.4007, 3.2329, 3.6723 & 2102 &, 3.6845, 3.6876, 3.6197, 3.4799, 3.3737, 4.4341, 4.0525 & 2103 &, 3.9011, 3.8945, 3.8635, 3.8368, 3.8153, 3.7936, 3.7758 & 2104 &, 3.4944, 3.4873, 3.9040, 3.7110, 3.6922, 3.6799, 3.6724 & 2105 &, 3.5622, 3.6081, 3.5426, 3.4922, 3.4498, 3.3984, 3.4456 & 2106 &, 2.7522, 2.5524, 3.5742, 2.9508, 3.0751, 3.0158, 2.9644 & 2107 &/) 2108 r0ab(3011:3080)=(/ & 2109 & 2.8338, 2.7891, 2.6933, 3.6926, 3.1814, 3.4528, 3.4186 & 2110 &, 3.3836, 3.2213, 3.1626, 3.0507, 4.0548, 3.5312, 3.4244 & 2111 &, 3.3409, 3.2810, 3.4782, 3.1905, 3.1494, 3.1221, 3.1128 & 2112 &, 3.0853, 3.0384, 3.4366, 3.4562, 3.4638, 3.3211, 3.2762 & 2113 &, 3.1730, 4.1632, 3.6825, 3.5822, 3.4870, 3.6325, 3.5740 & 2114 &, 3.4733, 3.4247, 3.3969, 3.3764, 3.3525, 3.1984, 3.5989 & 2115 &, 3.6299, 3.6433, 3.4937, 3.4417, 3.3365, 4.3304, 3.9242 & 2116 &, 3.7793, 3.7623, 3.7327, 3.7071, 3.6860, 3.6650, 3.6476 & 2117 &, 3.3849, 3.3534, 3.8216, 3.5870, 3.5695, 3.5584, 3.5508 & 2118 &, 3.4856, 3.5523, 3.4934, 3.4464, 3.4055, 3.3551, 3.3888 & 2119 &/) 2120 r0ab(3081:3150)=(/ & 2121 & 3.3525, 2.7202, 2.5183, 3.4947, 2.8731, 3.0198, 3.1457 & 2122 &, 2.9276, 2.7826, 2.7574, 2.6606, 3.6090, 3.0581, 3.3747 & 2123 &, 3.3677, 3.3450, 3.1651, 3.1259, 3.0147, 3.9498, 3.3857 & 2124 &, 3.2917, 3.2154, 3.1604, 3.4174, 3.0735, 3.0342, 3.0096 & 2125 &, 3.0136, 2.9855, 2.9680, 3.3604, 3.4037, 3.4243, 3.2633 & 2126 &, 3.1810, 3.1351, 4.0557, 3.5368, 3.4526, 3.3699, 3.5707 & 2127 &, 3.5184, 3.4085, 3.3595, 3.3333, 3.3143, 3.3041, 3.1094 & 2128 &, 3.5193, 3.5745, 3.6025, 3.4338, 3.3448, 3.2952, 4.2158 & 2129 &, 3.7802, 3.6431, 3.6129, 3.5853, 3.5610, 3.5406, 3.5204 & 2130 &, 3.5036, 3.2679, 3.2162, 3.7068, 3.4483, 3.4323, 3.4221 & 2131 &/) 2132 r0ab(3151:3220)=(/ & 2133 & 3.4138, 3.3652, 3.4576, 3.4053, 3.3618, 3.3224, 3.2711 & 2134 &, 3.3326, 3.2950, 3.2564, 2.5315, 2.6104, 3.2734, 3.2299 & 2135 &, 3.1090, 2.9942, 2.9159, 2.8324, 2.8350, 2.7216, 3.3994 & 2136 &, 3.4475, 3.4354, 3.3438, 3.2807, 3.2169, 3.2677, 3.1296 & 2137 &, 3.7493, 3.8075, 3.6846, 3.6104, 3.5577, 3.2052, 3.4803 & 2138 &, 3.4236, 3.3845, 3.3640, 3.3365, 3.3010, 3.3938, 3.3624 & 2139 &, 3.3440, 3.3132, 3.4035, 3.2754, 3.8701, 3.9523, 3.8018 & 2140 &, 3.7149, 3.3673, 3.3199, 3.2483, 3.2069, 3.1793, 3.1558 & 2141 &, 3.1395, 3.4097, 3.5410, 3.5228, 3.5116, 3.4921, 3.4781 & 2142 &, 3.4690, 4.0420, 4.1759, 4.0078, 4.0450, 4.0189, 3.9952 & 2143 &/) 2144 r0ab(3221:3290)=(/ & 2145 & 3.9770, 3.9583, 3.9434, 3.7217, 3.8228, 3.7826, 3.8640 & 2146 &, 3.8446, 3.8314, 3.8225, 3.6817, 3.7068, 3.6555, 3.6159 & 2147 &, 3.5831, 3.5257, 3.2133, 3.1689, 3.1196, 3.3599, 2.9852 & 2148 &, 2.7881, 3.5284, 3.3493, 3.6958, 3.3642, 3.1568, 3.0055 & 2149 &, 2.9558, 2.8393, 3.6287, 3.5283, 4.1511, 3.8259, 3.6066 & 2150 &, 3.4527, 3.3480, 3.2713, 3.9037, 3.8361, 3.8579, 3.7311 & 2151 &, 3.6575, 3.5176, 3.5693, 3.5157, 3.4814, 3.4559, 3.4445 & 2152 &, 3.4160, 4.1231, 3.8543, 3.6816, 3.5602, 3.4798, 3.4208 & 2153 &, 4.0542, 4.0139, 4.0165, 3.9412, 3.7698, 3.6915, 3.6043 & 2154 &, 3.5639, 3.5416, 3.5247, 3.5153, 3.5654, 4.2862, 4.0437 & 2155 &/) 2156 r0ab(3291:3360)=(/ & 2157 & 3.8871, 3.7741, 3.6985, 3.6413, 4.2345, 4.3663, 4.3257 & 2158 &, 4.0869, 4.0612, 4.0364, 4.0170, 3.9978, 3.9834, 3.9137 & 2159 &, 3.8825, 3.8758, 3.9143, 3.8976, 3.8864, 3.8768, 3.9190 & 2160 &, 4.1613, 4.0566, 3.9784, 3.9116, 3.8326, 3.7122, 3.6378 & 2161 &, 3.5576, 3.5457, 4.3127, 3.1160, 2.8482, 4.0739, 3.3599 & 2162 &, 3.5698, 3.5366, 3.2854, 3.1039, 2.9953, 2.9192, 4.1432 & 2163 &, 3.5320, 3.9478, 4.0231, 3.7509, 3.5604, 3.4340, 3.3426 & 2164 &, 4.3328, 3.8288, 3.7822, 3.7909, 3.6907, 3.6864, 3.5793 & 2165 &, 3.5221, 3.4883, 3.4649, 3.4514, 3.4301, 3.9256, 4.0596 & 2166 &, 3.8307, 3.6702, 3.5651, 3.4884, 4.4182, 4.2516, 3.9687 & 2167 &/) 2168 r0ab(3361:3430)=(/ & 2169 & 3.9186, 3.9485, 3.8370, 3.7255, 3.6744, 3.6476, 3.6295 & 2170 &, 3.6193, 3.5659, 4.0663, 4.2309, 4.0183, 3.8680, 3.7672 & 2171 &, 3.6923, 4.5240, 4.4834, 4.1570, 4.3204, 4.2993, 4.2804 & 2172 &, 4.2647, 4.2481, 4.2354, 3.8626, 3.8448, 4.2267, 4.1799 & 2173 &, 4.1670, 3.8738, 3.8643, 3.8796, 4.0575, 4.0354, 3.9365 & 2174 &, 3.8611, 3.7847, 3.7388, 3.6826, 3.6251, 3.5492, 4.0889 & 2175 &, 4.2764, 3.1416, 2.8325, 3.7735, 3.3787, 3.4632, 3.5923 & 2176 &, 3.3214, 3.1285, 3.0147, 2.9366, 3.8527, 3.5602, 3.8131 & 2177 &, 3.8349, 3.7995, 3.5919, 3.4539, 3.3540, 4.0654, 3.8603 & 2178 &, 3.7972, 3.7358, 3.7392, 3.8157, 3.6055, 3.5438, 3.5089 & 2179 &/) 2180 r0ab(3431:3500)=(/ & 2181 & 3.4853, 3.4698, 3.4508, 3.7882, 3.8682, 3.8837, 3.7055 & 2182 &, 3.5870, 3.5000, 4.1573, 4.0005, 3.9568, 3.8936, 3.9990 & 2183 &, 3.9433, 3.8172, 3.7566, 3.7246, 3.7033, 3.6900, 3.5697 & 2184 &, 3.9183, 4.0262, 4.0659, 3.8969, 3.7809, 3.6949, 4.2765 & 2185 &, 4.2312, 4.1401, 4.0815, 4.0580, 4.0369, 4.0194, 4.0017 & 2186 &, 3.9874, 3.8312, 3.8120, 3.9454, 3.9210, 3.9055, 3.8951 & 2187 &, 3.8866, 3.8689, 3.9603, 3.9109, 3.9122, 3.8233, 3.7438 & 2188 &, 3.7436, 3.6981, 3.6555, 3.5452, 3.9327, 4.0658, 4.1175 & 2189 &, 2.9664, 2.8209, 3.5547, 3.3796, 3.3985, 3.3164, 3.2364 & 2190 &, 3.1956, 3.0370, 2.9313, 3.6425, 3.5565, 3.7209, 3.7108 & 2191 &/) 2192 r0ab(3501:3570)=(/ & 2193 & 3.6639, 3.6484, 3.4745, 3.3492, 3.8755, 4.2457, 3.7758 & 2194 &, 3.7161, 3.6693, 3.6155, 3.5941, 3.5643, 3.5292, 3.4950 & 2195 &, 3.4720, 3.4503, 3.6936, 3.7392, 3.7388, 3.7602, 3.6078 & 2196 &, 3.4960, 3.9800, 4.3518, 4.2802, 3.8580, 3.8056, 3.7527 & 2197 &, 3.7019, 3.6615, 3.5768, 3.5330, 3.5038, 3.5639, 3.8192 & 2198 &, 3.8883, 3.9092, 3.9478, 3.7995, 3.6896, 4.1165, 4.5232 & 2199 &, 4.4357, 4.4226, 4.4031, 4.3860, 4.3721, 4.3580, 4.3466 & 2200 &, 4.2036, 4.2037, 3.8867, 4.2895, 4.2766, 4.2662, 4.2598 & 2201 &, 3.8408, 3.9169, 3.8681, 3.8250, 3.7855, 3.7501, 3.6753 & 2202 &, 3.5499, 3.4872, 3.5401, 3.8288, 3.9217, 3.9538, 4.0054 & 2203 &/) 2204 r0ab(3571:3640)=(/ & 2205 & 2.8388, 2.7890, 3.4329, 3.5593, 3.3488, 3.2486, 3.1615 & 2206 &, 3.1000, 3.0394, 2.9165, 3.5267, 3.7479, 3.6650, 3.6263 & 2207 &, 3.5658, 3.5224, 3.4762, 3.3342, 3.7738, 4.0333, 3.9568 & 2208 &, 3.8975, 3.8521, 3.4929, 3.7830, 3.7409, 3.7062, 3.6786 & 2209 &, 3.6471, 3.6208, 3.6337, 3.6519, 3.6363, 3.6278, 3.6110 & 2210 &, 3.4825, 3.8795, 4.1448, 4.0736, 4.0045, 3.6843, 3.6291 & 2211 &, 3.5741, 3.5312, 3.4974, 3.4472, 3.4034, 3.7131, 3.7557 & 2212 &, 3.7966, 3.8005, 3.8068, 3.8015, 3.6747, 4.0222, 4.3207 & 2213 &, 4.2347, 4.2191, 4.1990, 4.1811, 4.1666, 4.1521, 4.1401 & 2214 &, 3.9970, 3.9943, 3.9592, 4.0800, 4.0664, 4.0559, 4.0488 & 2215 &/) 2216 r0ab(3641:3710)=(/ & 2217 & 3.9882, 4.0035, 3.9539, 3.9138, 3.8798, 3.8355, 3.5359 & 2218 &, 3.4954, 3.3962, 3.5339, 3.7595, 3.8250, 3.8408, 3.8600 & 2219 &, 3.8644, 2.7412, 2.7489, 3.3374, 3.3950, 3.3076, 3.1910 & 2220 &, 3.0961, 3.0175, 3.0280, 2.8929, 3.4328, 3.5883, 3.6227 & 2221 &, 3.5616, 3.4894, 3.4241, 3.3641, 3.3120, 3.6815, 3.8789 & 2222 &, 3.8031, 3.7413, 3.6939, 3.4010, 3.6225, 3.5797, 3.5443 & 2223 &, 3.5139, 3.4923, 3.4642, 3.5860, 3.5849, 3.5570, 3.5257 & 2224 &, 3.4936, 3.4628, 3.7874, 3.9916, 3.9249, 3.8530, 3.5932 & 2225 &, 3.5355, 3.4757, 3.4306, 3.3953, 3.3646, 3.3390, 3.5637 & 2226 &, 3.7053, 3.7266, 3.7177, 3.6996, 3.6775, 3.6558, 3.9331 & 2227 &/) 2228 r0ab(3711:3780)=(/ & 2229 & 4.1655, 4.0879, 4.0681, 4.0479, 4.0299, 4.0152, 4.0006 & 2230 &, 3.9883, 3.8500, 3.8359, 3.8249, 3.9269, 3.9133, 3.9025 & 2231 &, 3.8948, 3.8422, 3.8509, 3.7990, 3.7570, 3.7219, 3.6762 & 2232 &, 3.4260, 3.3866, 3.3425, 3.5294, 3.7022, 3.7497, 3.7542 & 2233 &, 3.7494, 3.7370, 3.7216, 3.4155, 3.0522, 4.2541, 3.8218 & 2234 &, 4.0438, 3.5875, 3.3286, 3.1682, 3.0566, 2.9746, 4.3627 & 2235 &, 4.0249, 4.6947, 4.1718, 3.8639, 3.6735, 3.5435, 3.4479 & 2236 &, 4.6806, 4.3485, 4.2668, 4.1690, 4.1061, 4.1245, 4.0206 & 2237 &, 3.9765, 3.9458, 3.9217, 3.9075, 3.8813, 3.9947, 4.1989 & 2238 &, 3.9507, 3.7960, 3.6925, 3.6150, 4.8535, 4.5642, 4.4134 & 2239 &/) 2240 r0ab(3781:3850)=(/ & 2241 & 4.3688, 4.3396, 4.2879, 4.2166, 4.1888, 4.1768, 4.1660 & 2242 &, 4.1608, 4.0745, 4.2289, 4.4863, 4.2513, 4.0897, 3.9876 & 2243 &, 3.9061, 5.0690, 5.0446, 4.6186, 4.6078, 4.5780, 4.5538 & 2244 &, 4.5319, 4.5101, 4.4945, 4.1912, 4.2315, 4.5534, 4.4373 & 2245 &, 4.4224, 4.4120, 4.4040, 4.2634, 4.7770, 4.6890, 4.6107 & 2246 &, 4.5331, 4.4496, 4.4082, 4.3095, 4.2023, 4.0501, 4.2595 & 2247 &, 4.5497, 4.3056, 4.1506, 4.0574, 3.9725, 5.0796, 3.0548 & 2248 &, 3.3206, 3.8132, 3.9720, 3.7675, 3.7351, 3.5167, 3.5274 & 2249 &, 3.3085, 3.1653, 3.9500, 4.1730, 4.0613, 4.1493, 3.8823 & 2250 &, 4.0537, 3.8200, 3.6582, 4.3422, 4.5111, 4.3795, 4.3362 & 2251 &/) 2252 r0ab(3851:3920)=(/ & 2253 & 4.2751, 3.7103, 4.1973, 4.1385, 4.1129, 4.0800, 4.0647 & 2254 &, 4.0308, 4.0096, 4.1619, 3.9360, 4.1766, 3.9705, 3.8262 & 2255 &, 4.5348, 4.7025, 4.5268, 4.5076, 3.9562, 3.9065, 3.8119 & 2256 &, 3.7605, 3.7447, 3.7119, 3.6916, 4.1950, 4.2110, 4.3843 & 2257 &, 4.1631, 4.4427, 4.2463, 4.1054, 4.7693, 5.0649, 4.7365 & 2258 &, 4.7761, 4.7498, 4.7272, 4.7076, 4.6877, 4.6730, 4.4274 & 2259 &, 4.5473, 4.5169, 4.5975, 4.5793, 4.5667, 4.5559, 4.3804 & 2260 &, 4.6920, 4.6731, 4.6142, 4.5600, 4.4801, 4.0149, 3.8856 & 2261 &, 3.7407, 4.1545, 4.2253, 4.4229, 4.1923, 4.5022, 4.3059 & 2262 &, 4.1591, 4.7883, 4.9294, 3.3850, 3.4208, 3.7004, 3.8800 & 2263 &/) 2264 r0ab(3921:3990)=(/ & 2265 & 3.9886, 3.9040, 3.6719, 3.6547, 3.4625, 3.3370, 3.8394 & 2266 &, 4.0335, 4.2373, 4.3023, 4.0306, 4.1408, 3.9297, 3.7857 & 2267 &, 4.1907, 4.3230, 4.2664, 4.2173, 4.1482, 3.6823, 4.0711 & 2268 &, 4.0180, 4.0017, 3.9747, 3.9634, 3.9383, 4.1993, 4.3205 & 2269 &, 4.0821, 4.2547, 4.0659, 3.9359, 4.3952, 4.5176, 4.3888 & 2270 &, 4.3607, 3.9583, 3.9280, 3.8390, 3.7971, 3.7955, 3.7674 & 2271 &, 3.7521, 4.1062, 4.3633, 4.2991, 4.2767, 4.4857, 4.3039 & 2272 &, 4.1762, 4.6197, 4.8654, 4.6633, 4.5878, 4.5640, 4.5422 & 2273 &, 4.5231, 4.5042, 4.4901, 4.3282, 4.3978, 4.3483, 4.4202 & 2274 &, 4.4039, 4.3926, 4.3807, 4.2649, 4.6135, 4.5605, 4.5232 & 2275 &/) 2276 r0ab(3991:4060)=(/ & 2277 & 4.4676, 4.3948, 4.0989, 3.9864, 3.8596, 4.0942, 4.2720 & 2278 &, 4.3270, 4.3022, 4.5410, 4.3576, 4.2235, 4.6545, 4.7447 & 2279 &, 4.7043, 3.0942, 3.2075, 3.5152, 3.6659, 3.8289, 3.7459 & 2280 &, 3.5156, 3.5197, 3.3290, 3.2069, 3.6702, 3.8448, 4.0340 & 2281 &, 3.9509, 3.8585, 3.9894, 3.7787, 3.6365, 4.1425, 4.1618 & 2282 &, 4.0940, 4.0466, 3.9941, 3.5426, 3.8952, 3.8327, 3.8126 & 2283 &, 3.7796, 3.7635, 3.7356, 4.0047, 3.9655, 3.9116, 4.1010 & 2284 &, 3.9102, 3.7800, 4.2964, 4.3330, 4.2622, 4.2254, 3.8195 & 2285 &, 3.7560, 3.6513, 3.5941, 3.5810, 3.5420, 3.5178, 3.8861 & 2286 &, 4.1459, 4.1147, 4.0772, 4.3120, 4.1207, 3.9900, 4.4733 & 2287 &/) 2288 r0ab(4061:4130)=(/ & 2289 & 4.6157, 4.4580, 4.4194, 4.3954, 4.3739, 4.3531, 4.3343 & 2290 &, 4.3196, 4.2140, 4.2339, 4.1738, 4.2458, 4.2278, 4.2158 & 2291 &, 4.2039, 4.1658, 4.3595, 4.2857, 4.2444, 4.1855, 4.1122 & 2292 &, 3.7839, 3.6879, 3.5816, 3.8633, 4.1585, 4.1402, 4.1036 & 2293 &, 4.3694, 4.1735, 4.0368, 4.5095, 4.5538, 4.5240, 4.4252 & 2294 &, 3.0187, 3.1918, 3.5127, 3.6875, 3.7404, 3.6943, 3.4702 & 2295 &, 3.4888, 3.2914, 3.1643, 3.6669, 3.8724, 3.9940, 4.0816 & 2296 &, 3.8054, 3.9661, 3.7492, 3.6024, 4.0428, 4.1951, 4.1466 & 2297 &, 4.0515, 4.0075, 3.5020, 3.9158, 3.8546, 3.8342, 3.8008 & 2298 &, 3.7845, 3.7549, 3.9602, 3.8872, 3.8564, 4.0793, 3.8835 & 2299 &/) 2300 r0ab(4131:4200)=(/ & 2301 & 3.7495, 4.2213, 4.3704, 4.3300, 4.2121, 3.7643, 3.7130 & 2302 &, 3.6144, 3.5599, 3.5474, 3.5093, 3.4853, 3.9075, 4.1115 & 2303 &, 4.0473, 4.0318, 4.2999, 4.1050, 3.9710, 4.4320, 4.6706 & 2304 &, 4.5273, 4.4581, 4.4332, 4.4064, 4.3873, 4.3684, 4.3537 & 2305 &, 4.2728, 4.2549, 4.2032, 4.2794, 4.2613, 4.2491, 4.2375 & 2306 &, 4.2322, 4.3665, 4.3061, 4.2714, 4.2155, 4.1416, 3.7660 & 2307 &, 3.6628, 3.5476, 3.8790, 4.1233, 4.0738, 4.0575, 4.3575 & 2308 &, 4.1586, 4.0183, 4.4593, 4.5927, 4.4865, 4.3813, 4.4594 & 2309 &, 2.9875, 3.1674, 3.4971, 3.6715, 3.7114, 3.6692, 3.4446 & 2310 &, 3.4676, 3.2685, 3.1405, 3.6546, 3.8579, 3.9637, 4.0581 & 2311 &/) 2312 r0ab(4201:4270)=(/ & 2313 & 3.7796, 3.9463, 3.7275, 3.5792, 4.0295, 4.1824, 4.1247 & 2314 &, 4.0357, 3.9926, 3.4827, 3.9007, 3.8392, 3.8191, 3.7851 & 2315 &, 3.7687, 3.7387, 3.9290, 3.8606, 3.8306, 4.0601, 3.8625 & 2316 &, 3.7269, 4.2062, 4.3566, 4.3022, 4.1929, 3.7401, 3.6888 & 2317 &, 3.5900, 3.5350, 3.5226, 3.4838, 3.4594, 3.8888, 4.0813 & 2318 &, 4.0209, 4.0059, 4.2810, 4.0843, 3.9486, 4.4162, 4.6542 & 2319 &, 4.5005, 4.4444, 4.4196, 4.3933, 4.3741, 4.3552, 4.3406 & 2320 &, 4.2484, 4.2413, 4.1907, 4.2656, 4.2474, 4.2352, 4.2236 & 2321 &, 4.2068, 4.3410, 4.2817, 4.2479, 4.1921, 4.1182, 3.7346 & 2322 &, 3.6314, 3.5168, 3.8582, 4.0927, 4.0469, 4.0313, 4.3391 & 2323 &/) 2324 r0ab(4271:4340)=(/ & 2325 & 4.1381, 3.9962, 4.4429, 4.5787, 4.4731, 4.3588, 4.4270 & 2326 &, 4.3957, 2.9659, 3.1442, 3.4795, 3.6503, 3.6814, 3.6476 & 2327 &, 3.4222, 3.4491, 3.2494, 3.1209, 3.6324, 3.8375, 3.9397 & 2328 &, 3.8311, 3.7581, 3.9274, 3.7085, 3.5598, 4.0080, 4.1641 & 2329 &, 4.1057, 4.0158, 3.9726, 3.4667, 3.8802, 3.8188, 3.7989 & 2330 &, 3.7644, 3.7474, 3.7173, 3.9049, 3.8424, 3.8095, 4.0412 & 2331 &, 3.8436, 3.7077, 4.1837, 4.3366, 4.2816, 4.1686, 3.7293 & 2332 &, 3.6709, 3.5700, 3.5153, 3.5039, 3.4684, 3.4437, 3.8663 & 2333 &, 4.0575, 4.0020, 3.9842, 4.2612, 4.0643, 3.9285, 4.3928 & 2334 &, 4.6308, 4.4799, 4.4244, 4.3996, 4.3737, 4.3547, 4.3358 & 2335 &/) 2336 r0ab(4341:4410)=(/ & 2337 & 4.3212, 4.2275, 4.2216, 4.1676, 4.2465, 4.2283, 4.2161 & 2338 &, 4.2045, 4.1841, 4.3135, 4.2562, 4.2226, 4.1667, 4.0932 & 2339 &, 3.7134, 3.6109, 3.4962, 3.8352, 4.0688, 4.0281, 4.0099 & 2340 &, 4.3199, 4.1188, 3.9768, 4.4192, 4.5577, 4.4516, 4.3365 & 2341 &, 4.4058, 4.3745, 4.3539, 2.8763, 3.1294, 3.5598, 3.7465 & 2342 &, 3.5659, 3.5816, 3.3599, 3.4024, 3.1877, 3.0484, 3.7009 & 2343 &, 3.9451, 3.8465, 3.9873, 3.7079, 3.9083, 3.6756, 3.5150 & 2344 &, 4.0829, 4.2780, 4.1511, 4.1260, 4.0571, 3.4865, 3.9744 & 2345 &, 3.9150, 3.8930, 3.8578, 3.8402, 3.8073, 3.7977, 4.0036 & 2346 &, 3.7604, 4.0288, 3.8210, 3.6757, 4.2646, 4.4558, 4.2862 & 2347 &/) 2348 r0ab(4411:4465)=(/ & 2349 & 4.2122, 3.7088, 3.6729, 3.5800, 3.5276, 3.5165, 3.4783 & 2350 &, 3.4539, 3.9553, 3.9818, 4.2040, 3.9604, 4.2718, 4.0689 & 2351 &, 3.9253, 4.4869, 4.7792, 4.4918, 4.5342, 4.5090, 4.4868 & 2352 &, 4.4680, 4.4486, 4.4341, 4.2023, 4.3122, 4.2710, 4.3587 & 2353 &, 4.3407, 4.3281, 4.3174, 4.1499, 4.3940, 4.3895, 4.3260 & 2354 &, 4.2725, 4.1961, 3.7361, 3.6193, 3.4916, 3.9115, 3.9914 & 2355 &, 3.9809, 3.9866, 4.3329, 4.1276, 3.9782, 4.5097, 4.6769 & 2356 &, 4.5158, 4.3291, 4.3609, 4.3462, 4.3265, 4.4341 & 2357 &/) 2358 2359 k=0 2360 do i=1,max_elem 2361 do j=1,i 2362 k=k+1 2363 r(i,j)=r0ab(k)/autoang 2364 r(j,i)=r0ab(k)/autoang 2365 end do 2366 end do 2367 2368 end subroutine setr0ab 2369 2370 2371 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2372 2373 subroutine stoprun(s) 2374 character*(*) s 2375 write(*,*)'program stopped due to: ',s 2376 open(99, file='dscf_problem', action='write', status='replace') 2377 close(99) 2378 !call system('touch dscf_problem') 2379 stop 'must stop!' 2380 end subroutine stoprun 2381 2382 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2383 ! Returns the number of a given element string (h-pu, 1-94) 2384 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2385 2386 elemental subroutine ELEM(KEY1, NAT) 2387 CHARACTER(*), intent(in) :: KEY1 2388 INTEGER, intent(out) :: NAT 2389 2390 character(2), parameter :: ELEMNT(94) = [& 2391 & 'h ','he', & 2392 & 'li','be','b ','c ','n ','o ','f ','ne', & 2393 & 'na','mg','al','si','p ','s ','cl','ar', & 2394 & 'k ','ca','sc','ti','v ','cr','mn','fe','co','ni','cu', & 2395 & 'zn','ga','ge','as','se','br','kr', & 2396 & 'rb','sr','y ','zr','nb','mo','tc','ru','rh','pd','ag', & 2397 & 'cd','in','sn','sb','te','i ','xe', & 2398 & 'cs','ba','la','ce','pr','nd','pm','sm','eu','gd','tb','dy', & 2399 & 'ho','er','tm','yb','lu','hf','ta','w ','re','os','ir','pt', & 2400 & 'au','hg','tl','pb','bi','po','at','rn', & 2401 & 'fr','ra','ac','th','pa','u ','np','pu' ] 2402 2403 CHARACTER(2) :: E 2404 integer :: k, j, n, i 2405 2406 nat=0 2407 e=' ' 2408 k=1 2409 do J=1,len(key1) 2410 if (k.gt.2)exit 2411 N=ICHAR(key1(J:J)) 2412 if (n.ge.ichar('A') .and. n.le.ichar('Z') )then 2413 e(k:k)=char(n+ICHAR('a')-ICHAR('A')) 2414 k=k+1 2415 end if 2416 if (n.ge.ichar('a') .and. n.le.ichar('z') )then 2417 e(k:k)=key1(j:j) 2418 k=k+1 2419 end if 2420 end do 2421 2422 do I=1,94 2423 if (e.eq.elemnt(i))then 2424 NAT=I 2425 RETURN 2426 end if 2427 end do 2428 2429 2430 end subroutine ELEM 2431 2432 2433 2434 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2435 ! B E G I N O F P B C P A R T 2436 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2437 2438 2439 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2440 ! compute coordination numbers by adding an inverse damping function 2441 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2442 2443 subroutine pbcncoord(natoms,rcov,iz,xyz,cn,lat,rep_cn,crit_cn) 2444 integer,intent(in) :: natoms,iz(*) 2445 real(wp),intent(in) :: rcov(94) 2446 2447 integer i,rep_cn(3) 2448 real(wp) xyz(3,*),cn(*),lat(3,3) 2449 2450 integer iat,taux,tauy,tauz 2451 real(wp) dx,dy,dz,r,damp,xn,rr,rco,tau(3) 2452 real(wp), INTENT(IN) :: crit_cn 2453 2454 do i=1,natoms 2455 xn=0.0d0 2456 do iat=1,natoms 2457 do taux=-rep_cn(1),rep_cn(1) 2458 do tauy=-rep_cn(2),rep_cn(2) 2459 do tauz=-rep_cn(3),rep_cn(3) 2460 if (iat.eq.i .and. taux.eq.0 .and. tauy.eq.0 .and.& 2461 & tauz.eq.0) cycle 2462 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 2463 dx=xyz(1,iat)-xyz(1,i)+tau(1) 2464 dy=xyz(2,iat)-xyz(2,i)+tau(2) 2465 dz=xyz(3,iat)-xyz(3,i)+tau(3) 2466 r=(dx*dx+dy*dy+dz*dz) 2467 if (r.gt.crit_cn) cycle 2468 r=sqrt(r) 2469 ! covalent distance in Bohr 2470 rco=rcov(iz(i))+rcov(iz(iat)) 2471 rr=rco/r 2472 ! counting function exponential has a better long-range behavior than MH 2473 damp=1.d0/(1.d0+exp(-k1*(rr-1.0d0))) 2474 xn=xn+damp 2475 ! print '("cn(",I2,I2,"): ",E14.8)',i,iat,damp 2476 2477 end do 2478 end do 2479 end do 2480 end do 2481 cn(i)=xn 2482 end do 2483 2484 end subroutine pbcncoord 2485 2486 2487 2488 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2489 ! compute energy 2490 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 2491 2492 subroutine pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& 2493 & rcov,rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,& 2494 & e6,e8,e10,e12,e63,lat,rthr,rep_vdw,cn_thr,rep_cn) 2495 integer max_elem,maxc 2496 real(wp) r2r4(max_elem),rcov(max_elem) 2497 real(wp) rs6,rs8,rs10,alp6,alp8,alp10 2498 real(wp) rthr,cn_thr,crit_cn 2499 integer rep_vdw(3),rep_cn(3) 2500 integer n,iz(*),version,mxc(max_elem) 2501 ! integer rep_v(3)=rep_vdw!,rep_cn(3) 2502 real(wp) xyz(3,*),r0ab(max_elem,max_elem),lat(3,3) 2503 ! real(wp) rs6,rs8,rs10,alp6,alp8,alp10,rcov(max_elem) 2504 real(wp) c6ab(max_elem,max_elem,maxc,maxc,3) 2505 real(wp) e6, e8, e10, e12, e63 2506 logical noabc 2507 2508 integer iat,jat 2509 real(wp) r,r2,r6,r8,tmp,dx,dy,dz,c6,c8,R0 2510 real(wp) damp6,damp8,rr,r42 2511 real(wp) cn(n),rxyz(3),dxyz(3) 2512 real(wp) cc6ab(n*n),tau(3) 2513 integer ij 2514 integer taux,tauy,tauz,counter 2515 real(wp) a1,a2 2516 real(wp) bj_dmp6,bj_dmp8 2517 2518 2519 e6 =0 2520 e8 =0 2521 e10=0 2522 e12=0 2523 e63=0 2524 tau=(/0.0,0.0,0.0/) 2525 counter=0 2526 crit_cn=cn_thr 2527 ! Becke-Johnson parameters 2528 a1=rs6 2529 a2=rs8 2530 2531 2532 2533 ! DFT-D2 2534 if (version.eq.2)then 2535 2536 2537 do iat=1,n-1 2538 do jat=iat+1,n 2539 c6=c6ab(iz(jat),iz(iat),1,1,1) 2540 do taux=-rep_vdw(1),rep_vdw(1) 2541 do tauy=-rep_vdw(2),rep_vdw(2) 2542 do tauz=-rep_vdw(3),rep_vdw(3) 2543 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 2544 dx=xyz(1,iat)-xyz(1,jat)+tau(1) 2545 dy=xyz(2,iat)-xyz(2,jat)+tau(2) 2546 dz=xyz(3,iat)-xyz(3,jat)+tau(3) 2547 r2=dx*dx+dy*dy+dz*dz 2548 if (r2.gt.rthr) cycle 2549 r=sqrt(r2) 2550 damp6=1./(1.+exp(-alp6*(r/(rs6*r0ab(iz(jat),iz(iat)))-1.))) 2551 r6=r2**3 2552 e6 =e6+c6*damp6/r6 2553 end do 2554 end do 2555 end do 2556 end do 2557 end do 2558 2559 do iat=1,n 2560 jat=iat 2561 c6=c6ab(iz(jat),iz(iat),1,1,1) 2562 do taux=-rep_vdw(1),rep_vdw(1) 2563 do tauy=-rep_vdw(2),rep_vdw(2) 2564 do tauz=-rep_vdw(3),rep_vdw(3) 2565 if (taux.eq.0 .and. tauy.eq.0 .and. tauz.eq.0) cycle 2566 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 2567 dx=tau(1) 2568 dy=tau(2) 2569 dz=tau(3) 2570 r2=dx*dx+dy*dy+dz*dz 2571 if (r2.gt.rthr) cycle 2572 r=sqrt(r2) 2573 damp6=1./(1.+exp(-alp6*(r/(rs6*r0ab(iz(jat),iz(iat)))-1.))) 2574 r6=r2**3 2575 e6 =e6+c6*damp6/r6*0.50d0 2576 end do 2577 end do 2578 end do 2579 end do 2580 2581 2582 2583 else if ((version.eq.3).or.(version.eq.5)) then 2584 ! DFT-D3(zero-damping) 2585 2586 call pbcncoord(n,rcov,iz,xyz,cn,lat,rep_cn,crit_cn) 2587 2588 do iat=1,n-1 2589 do jat=iat+1,n 2590 ! get C6 2591 call getc6(maxc,max_elem,c6ab,mxc,iz(iat),iz(jat),& 2592 & cn(iat),cn(jat),c6) 2593 2594 if (.not.noabc)then 2595 ij=lin(jat,iat) 2596 ! store C6 for C9, calc as sqrt 2597 cc6ab(ij)=sqrt(c6) 2598 end if 2599 do taux=-rep_vdw(1),rep_vdw(1) 2600 do tauy=-rep_vdw(2),rep_vdw(2) 2601 do tauz=-rep_vdw(3),rep_vdw(3) 2602 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 2603 2604 dx=xyz(1,iat)-xyz(1,jat)+tau(1) 2605 dy=xyz(2,iat)-xyz(2,jat)+tau(2) 2606 dz=xyz(3,iat)-xyz(3,jat)+tau(3) 2607 r2=dx*dx+dy*dy+dz*dz 2608 ! cutoff 2609 2610 if (r2.gt.rthr) cycle 2611 r =sqrt(r2) 2612 R0=r0ab(iz(jat),iz(iat)) 2613 rr=R0/r 2614 ! damping 2615 if(version.eq.3)then 2616 ! DFT-D3 zero-damp 2617 tmp=rs6*rr 2618 damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) 2619 tmp=rs8*rr 2620 damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) 2621 else 2622 ! DFT-D3M zero-damp 2623 tmp=(r/(rs6*R0))+rs8*R0 2624 damp6 =1.d0/( 1.d0+6.d0*tmp**(-alp6) ) 2625 tmp=(r/R0)+rs8*R0 2626 damp8 =1.d0/( 1.d0+6.d0*tmp**(-alp8) ) 2627 endif 2628 2629 r6=r2**3 2630 e6 =e6+damp6/r6* c6 2631 ! write(*,*)'e6: ',c6*damp6/r6*autokcal 2632 2633 ! stored in main as sqrt 2634 c8 =3.0d0*r2r4(iz(iat))*r2r4(iz(jat))*c6 2635 r8 =r6*r2 2636 2637 e8 =e8+c8*damp8/r8 2638 2639 end do 2640 end do 2641 end do 2642 end do 2643 end do 2644 2645 do iat=1,n 2646 jat=iat 2647 ! get C6 2648 call getc6(maxc,max_elem,c6ab,mxc,iz(iat),iz(jat),& 2649 & cn(iat),cn(jat),c6) 2650 2651 if (.not.noabc)then 2652 ij=lin(jat,iat) 2653 ! store C6 for C9, calc as sqrt 2654 cc6ab(ij)=sqrt(c6) 2655 end if 2656 do taux=-rep_vdw(1),rep_vdw(1) 2657 do tauy=-rep_vdw(2),rep_vdw(2) 2658 do tauz=-rep_vdw(3),rep_vdw(3) 2659 if (taux.eq.0 .and. tauy.eq.0 .and. tauz.eq.0) cycle 2660 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 2661 2662 dx=tau(1) 2663 dy=tau(2) 2664 dz=tau(3) 2665 r2=dx*dx+dy*dy+dz*dz 2666 ! cutoff 2667 if (r2.gt.rthr) cycle 2668 r =sqrt(r2) 2669 R0=r0ab(iz(jat),iz(iat)) 2670 rr=R0/r 2671 2672 ! damping 2673 if(version.eq.3)then 2674 ! DFT-D3 zero-damp 2675 tmp=rs6*rr 2676 damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) 2677 tmp=rs8*rr 2678 damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) 2679 else 2680 ! DFT-D3M zero-damp 2681 tmp=(r/(rs6*R0))+rs8*R0 2682 damp6 =1.d0/( 1.d0+6.d0*tmp**(-alp6) ) 2683 tmp=(r/R0)+rs8*R0 2684 damp8 =1.d0/( 1.d0+6.d0*tmp**(-alp8) ) 2685 endif 2686 2687 2688 r6=r2**3 2689 2690 e6 =e6+damp6/r6*0.50d0 *C6 2691 2692 ! stored in main as sqrt 2693 c8 =3.0d0*r2r4(iz(iat))*r2r4(iz(jat)) *C6 2694 r8 =r6*r2 2695 2696 e8 =e8+c8*damp8/r8*0.50d0 2697 counter=counter+1 2698 2699 end do 2700 end do 2701 end do 2702 end do 2703 ! write(*,*)'counter(edisp): ',counter 2704 else if((version.eq.4).or.(version.eq.6)) then 2705 2706 2707 ! DFT-D3(BJ-damping) 2708 call pbcncoord(n,rcov,iz,xyz,cn,lat,rep_cn,crit_cn) 2709 2710 do iat=1,n 2711 do jat=iat+1,n 2712 ! get C6 2713 call getc6(maxc,max_elem,c6ab,mxc,iz(iat),iz(jat),& 2714 & cn(iat),cn(jat),c6) 2715 2716 rxyz=xyz(:,iat)-xyz(:,jat) 2717 r42=r2r4(iz(iat))*r2r4(iz(jat)) 2718 bj_dmp6=(a1*dsqrt(3.0d0*r42)+a2)**6 2719 bj_dmp8=(a1*dsqrt(3.0d0*r42)+a2)**8 2720 2721 if (.not.noabc)then 2722 ij=lin(jat,iat) 2723 ! store C6 for C9, calc as sqrt 2724 cc6ab(ij)=sqrt(c6) 2725 end if 2726 do taux=-rep_vdw(1),rep_vdw(1) 2727 do tauy=-rep_vdw(2),rep_vdw(2) 2728 do tauz=-rep_vdw(3),rep_vdw(3) 2729 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 2730 2731 dxyz=rxyz+tau 2732 2733 r2=sum(dxyz*dxyz) 2734 ! cutoff 2735 if (r2.gt.rthr) cycle 2736 r =sqrt(r2) 2737 rr=r0ab(iz(jat),iz(iat))/r 2738 2739 2740 r6=r2**3 2741 2742 e6 =e6+c6/(r6+bj_dmp6) 2743 2744 ! stored in main as sqrt 2745 c8 =3.0d0*c6*r42 2746 r8 =r6*r2 2747 2748 e8 =e8+c8/(r8+bj_dmp8) 2749 2750 counter=counter+1 2751 2752 end do 2753 end do 2754 end do 2755 end do 2756 2757 ! Now the self interaction 2758 jat=iat 2759 ! get C6 2760 call getc6(maxc,max_elem,c6ab,mxc,iz(iat),iz(jat),& 2761 & cn(iat),cn(jat),c6) 2762 r42=r2r4(iz(iat))*r2r4(iz(iat)) 2763 bj_dmp6=(a1*dsqrt(3.0d0*r42)+a2)**6 2764 bj_dmp8=(a1*dsqrt(3.0d0*r42)+a2)**8 2765 2766 if (.not.noabc)then 2767 ij=lin(jat,iat) 2768 ! store C6 for C9, calc as sqrt 2769 cc6ab(ij)=dsqrt(c6) 2770 end if 2771 2772 do taux=-rep_vdw(1),rep_vdw(1) 2773 do tauy=-rep_vdw(2),rep_vdw(2) 2774 do tauz=-rep_vdw(3),rep_vdw(3) 2775 if (taux.eq.0 .and. tauy.eq.0 .and. tauz.eq.0) cycle 2776 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 2777 2778 r2=sum(tau*tau) 2779 ! cutoff 2780 if (r2.gt.rthr) cycle 2781 r =sqrt(r2) 2782 rr=r0ab(iz(jat),iz(iat))/r 2783 2784 2785 r6=r2**3 2786 2787 e6 =e6+c6/(r6+bj_dmp6)*0.50d0 2788 2789 ! stored in main as sqrt 2790 c8 =3.0d0*c6*r42 2791 r8 =r6*r2 2792 2793 e8 =e8+c8/(r8+bj_dmp8)*0.50d0 2794 counter=counter+1 2795 2796 end do 2797 end do 2798 end do 2799 end do 2800 2801 2802 end if 2803 2804 2805 if (noabc)return 2806 2807 ! compute non-additive third-order energy using averaged C6 2808 call pbcthreebody(max_elem,xyz,lat,n,iz,rep_cn,crit_cn,& 2809 & cc6ab,r0ab,e63) 2810 2811 end subroutine pbcedisp 2812 2813 2814 subroutine pbcthreebody(max_elem,xyz,lat,n,iz,repv,cnthr,cc6ab,& 2815 & r0ab,eabc) 2816 integer max_elem 2817 INTEGER :: n,jtaux,jtauy,jtauz,iat,jat,kat 2818 INTEGER :: ktaux,ktauy,ktauz,counter,ij,ik,jk,idum 2819 REAL(WP) :: rij2,rik2,rjk2,c9,rr0ij,rr0ik 2820 REAL(WP) :: rr0jk,geomean,fdamp 2821 REAL(WP) :: r0ij,r0ik,r0jk 2822 REAL(WP),INTENT(OUT)::eabc 2823 REAL(WP) :: tmp1,tmp2,tmp3,tmp4,ang 2824 2825 REAL(WP) ,DIMENSION(3,3),INTENT(IN)::lat 2826 REAL(WP) ,DIMENSION(3,*),INTENT(IN) :: xyz 2827 INTEGER,DIMENSION(*),INTENT(IN)::iz 2828 REAL(WP),DIMENSION(3):: jtau,ktau,ijvec,ikvec,jkvec,dumvec 2829 INTEGER,DIMENSION(3):: repv 2830 REAL(WP),INTENT(IN) ::cnthr 2831 REAL(WP),DIMENSION(n*n),INTENT(IN)::cc6ab 2832 REAL(WP),DIMENSION(max_elem,max_elem),INTENT(IN):: r0ab 2833 REAL(WP),PARAMETER::sr9=0.75d0 2834 REAL(WP),PARAMETER::alp9=-16.0d0 2835 REAL(WP) :: abcthr 2836 INTEGER, DIMENSION(3) :: repmin,repmax 2837 ! REAL(WP) :: time1,time2 2838 2839 counter=0 2840 eabc=0.0d0 2841 abcthr=cnthr 2842 ! abcthr=1.0d99 2843 ! write(*,*)'thr:',(abcthr) 2844 2845 ! call cpu_time(time1) 2846 2847 do iat=3,n 2848 do jat=2,iat-1 2849 ijvec=xyz(:,jat)-xyz(:,iat) 2850 ij=lin(iat,jat) 2851 r0ij=r0ab(iz(iat),iz(jat)) 2852 do kat=1,jat-1 2853 ik=lin(iat,kat) 2854 jk=lin(jat,kat) 2855 ikvec=xyz(:,kat)-xyz(:,iat) 2856 jkvec=xyz(:,kat)-xyz(:,jat) 2857 c9=-1.0d0*(cc6ab(ij)*cc6ab(ik)*cc6ab(jk)) 2858 2859 r0ik=r0ab(iz(iat),iz(kat)) 2860 r0jk=r0ab(iz(jat),iz(kat)) 2861 2862 2863 do jtaux=-repv(1),repv(1) 2864 repmin(1)=max(-repv(1),jtaux-repv(1)) 2865 repmax(1)=min(repv(1),jtaux+repv(1)) 2866 do jtauy=-repv(2),repv(2) 2867 repmin(2)=max(-repv(2),jtauy-repv(2)) 2868 repmax(2)=min(repv(2),jtauy+repv(2)) 2869 do jtauz=-repv(3),repv(3) 2870 repmin(3)=max(-repv(3),jtauz-repv(3)) 2871 repmax(3)=min(repv(3),jtauz+repv(3)) 2872 jtau=jtaux*lat(:,1)+jtauy*lat(:,2)+jtauz*lat(:,3) 2873 dumvec=ijvec+jtau 2874 dumvec=dumvec*dumvec 2875 rij2=SUM(dumvec) 2876 if (rij2.gt.abcthr)cycle 2877 2878 rr0ij=DSQRT(rij2)/r0ij 2879 2880 2881 do ktaux=repmin(1),repmax(1) 2882 do ktauy=repmin(2),repmax(2) 2883 do ktauz=repmin(3),repmax(3) 2884 ktau=ktaux*lat(:,1)+ktauy*lat(:,2)+ktauz*lat(:,3) 2885 dumvec=ikvec+ktau 2886 dumvec=dumvec*dumvec 2887 rik2=SUM(dumvec) 2888 if (rik2.gt.abcthr)cycle 2889 rr0ik=DSQRT(rik2)/r0ik 2890 2891 dumvec=jkvec+ktau-jtau 2892 rjk2=SUM(dumvec*dumvec) 2893 if (rjk2.gt.abcthr)cycle 2894 rr0jk=DSQRT(rjk2)/r0jk 2895 2896 2897 geomean=(rr0ij*rr0ik*rr0jk)**(1.0d0/3.0d0) 2898 ! write(*,*)'geomean:',geomean 2899 fdamp=1./(1.+6.*(sr9*geomean)**alp9) 2900 tmp1 = (rij2+rjk2-rik2) 2901 tmp2 = (rij2+rik2-rjk2) 2902 tmp3 = (rik2+rjk2-rij2) 2903 tmp4=rij2*rjk2*rik2 2904 ang=(0.375d0*tmp1*tmp2*tmp3/tmp4+1.0d0)/tmp4**1.5d0 2905 2906 eabc=eabc+ang*c9*fdamp 2907 2908 end do 2909 end do 2910 end do 2911 2912 end do 2913 end do 2914 end do 2915 2916 end do 2917 end do 2918 end do 2919 2920 do iat=2,n 2921 jat=iat 2922 ij=lin(iat,jat) 2923 ijvec=0.0d0 2924 r0ij=r0ab(iz(iat),iz(jat)) 2925 do kat=1,iat-1 2926 jk=lin(jat,kat) 2927 ik=jk 2928 ikvec=xyz(:,kat)-xyz(:,iat) 2929 jkvec=ikvec 2930 c9=-(cc6ab(ij)*cc6ab(ik)*cc6ab(jk)) 2931 2932 r0ik=r0ab(iz(iat),iz(kat)) 2933 r0jk=r0ab(iz(jat),iz(kat)) 2934 do jtaux=-repv(1),repv(1) 2935 repmin(1)=max(-repv(1),jtaux-repv(1)) 2936 repmax(1)=min(repv(1),jtaux+repv(1)) 2937 do jtauy=-repv(2),repv(2) 2938 repmin(2)=max(-repv(2),jtauy-repv(2)) 2939 repmax(2)=min(repv(2),jtauy+repv(2)) 2940 do jtauz=-repv(3),repv(3) 2941 repmin(3)=max(-repv(3),jtauz-repv(3)) 2942 repmax(3)=min(repv(3),jtauz+repv(3)) 2943 if (jtaux.eq.0 .and. jtauy.eq.0 .and. jtauz.eq.0) cycle 2944 jtau=jtaux*lat(:,1)+jtauy*lat(:,2)+jtauz*lat(:,3) 2945 dumvec=ijvec+jtau 2946 dumvec=dumvec*dumvec 2947 rij2=SUM(dumvec) 2948 if (rij2.gt.abcthr)cycle 2949 2950 rr0ij=DSQRT(rij2)/r0ij 2951 2952 do ktaux=repmin(1),repmax(1) 2953 do ktauy=repmin(2),repmax(2) 2954 do ktauz=repmin(3),repmax(3) 2955 ! every result * 0.5 2956 ktau=ktaux*lat(:,1)+ktauy*lat(:,2)+ktauz*lat(:,3) 2957 dumvec=ikvec+ktau 2958 dumvec=dumvec*dumvec 2959 rik2=SUM(dumvec) 2960 if (rik2.gt.abcthr)cycle 2961 rr0ik=DSQRT(rik2)/r0ik 2962 2963 dumvec=jkvec+ktau-jtau 2964 dumvec=dumvec*dumvec 2965 rjk2=SUM(dumvec) 2966 if (rjk2.gt.abcthr)cycle 2967 rr0jk=DSQRT(rjk2)/r0jk 2968 2969 2970 geomean=(rr0ij*rr0ik*rr0jk)**(1./3.) 2971 fdamp=1./(1.+6.*(sr9*geomean)**alp9) 2972 tmp1 = (rij2+rjk2-rik2) 2973 tmp2 = (rij2+rik2-rjk2) 2974 tmp3 = (rik2+rjk2-rij2) 2975 tmp4=rij2*rjk2*rik2 2976 ang=(0.375d0*tmp1*tmp2*tmp3/tmp4+1.0d0)/tmp4**1.5d0 2977 2978 eabc=eabc+c9*fdamp*ang/2.0 2979 end do 2980 end do 2981 end do 2982 2983 end do 2984 end do 2985 end do 2986 end do 2987 end do 2988 2989 do iat=2,n 2990 do jat=1,iat-1 2991 kat=jat 2992 ij=lin(iat,jat) 2993 jk=lin(jat,kat) 2994 ik=ij 2995 ikvec=xyz(:,kat)-xyz(:,iat) 2996 ijvec=ikvec 2997 jkvec=0.0d0 2998 c9=-(cc6ab(ij)*cc6ab(ik)*cc6ab(jk)) 2999 3000 r0ij=r0ab(iz(iat),iz(jat)) 3001 r0ik=r0ij 3002 r0jk=r0ab(iz(jat),iz(kat)) 3003 3004 do jtaux=-repv(1),repv(1) 3005 repmin(1)=max(-repv(1),jtaux-repv(1)) 3006 repmax(1)=min(repv(1),jtaux+repv(1)) 3007 do jtauy=-repv(2),repv(2) 3008 repmin(2)=max(-repv(2),jtauy-repv(2)) 3009 repmax(2)=min(repv(2),jtauy+repv(2)) 3010 do jtauz=-repv(3),repv(3) 3011 repmin(3)=max(-repv(3),jtauz-repv(3)) 3012 repmax(3)=min(repv(3),jtauz+repv(3)) 3013 jtau=jtaux*lat(:,1)+jtauy*lat(:,2)+jtauz*lat(:,3) 3014 dumvec=ijvec+jtau 3015 dumvec=dumvec*dumvec 3016 rij2=SUM(dumvec) 3017 if (rij2.gt.abcthr)cycle 3018 3019 rr0ij=DSQRT(rij2)/r0ij 3020 3021 do ktaux=repmin(1),repmax(1) 3022 do ktauy=repmin(2),repmax(2) 3023 do ktauz=repmin(3),repmax(3) 3024 ! every result * 0.5 3025 if (jtaux.eq.ktaux .and. jtauy.eq.ktauy& 3026 & .and. jtauz.eq.ktauz) cycle 3027 ktau=ktaux*lat(:,1)+ktauy*lat(:,2)+ktauz*lat(:,3) 3028 dumvec=ikvec+ktau 3029 dumvec=dumvec*dumvec 3030 rik2=SUM(dumvec) 3031 if (rik2.gt.abcthr)cycle 3032 rr0ik=DSQRT(rik2)/r0ik 3033 3034 dumvec=jkvec+ktau-jtau 3035 dumvec=dumvec*dumvec 3036 rjk2=SUM(dumvec) 3037 if (rjk2.gt.abcthr)cycle 3038 rr0jk=DSQRT(rjk2)/r0jk 3039 3040 3041 geomean=(rr0ij*rr0ik*rr0jk)**(1./3.) 3042 fdamp=1./(1.+6.*(sr9*geomean)**alp9) 3043 tmp1 = (rij2+rjk2-rik2) 3044 tmp2 = (rij2+rik2-rjk2) 3045 tmp3 = (rik2+rjk2-rij2) 3046 tmp4=rij2*rjk2*rik2 3047 ang=(0.375d0*tmp1*tmp2*tmp3/tmp4+1.0d0)/tmp4**1.5d0 3048 3049 eabc=eabc+c9*fdamp*ang/2.0 3050 end do 3051 end do 3052 end do 3053 3054 end do 3055 end do 3056 end do 3057 end do 3058 end do 3059 3060 3061 ! And finally the self interaction iat=jat=kat all 3062 3063 idum=0 3064 do iat=1,n 3065 jat=iat 3066 kat=iat 3067 ijvec=0.0d0 3068 ij=lin(iat,iat) 3069 ik=ij 3070 jk=ij 3071 ikvec=ijvec 3072 jkvec=ikvec 3073 c9=-(cc6ab(ij)*cc6ab(ik)*cc6ab(jk)) 3074 3075 r0ij=r0ab(iz(iat),iz(iat)) 3076 r0ik=r0ij 3077 r0jk=r0ij 3078 do jtaux=-repv(1),repv(1) 3079 repmin(1)=max(-repv(1),jtaux-repv(1)) 3080 repmax(1)=min(repv(1),jtaux+repv(1)) 3081 do jtauy=-repv(2),repv(2) 3082 repmin(2)=max(-repv(2),jtauy-repv(2)) 3083 repmax(2)=min(repv(2),jtauy+repv(2)) 3084 do jtauz=-repv(3),repv(3) 3085 repmin(3)=max(-repv(3),jtauz-repv(3)) 3086 repmax(3)=min(repv(3),jtauz+repv(3)) 3087 if (jtaux.eq.0 .and. jtauy.eq.0 .and. jtauz.eq.0) cycle 3088 jtau=jtaux*lat(:,1)+jtauy*lat(:,2)+jtauz*lat(:,3) 3089 dumvec=jtau 3090 dumvec=dumvec*dumvec 3091 rij2=SUM(dumvec) 3092 if (rij2.gt.abcthr)cycle 3093 rr0ij=DSQRT(rij2)/r0ij 3094 3095 do ktaux=repmin(1),repmax(1) 3096 do ktauy=repmin(2),repmax(2) 3097 do ktauz=repmin(3),repmax(3) 3098 if ((ktaux.eq.0) .and.( ktauy.eq.0) .and.( ktauz.eq.0))cycle 3099 if ((ktaux.eq.jtaux) .and. (ktauy.eq.jtauy)& 3100 & .and. (ktauz.eq.jtauz)) cycle 3101 3102 ! every result * 1/6 becaues every triple is counted twice due to the tw 3103 ! 3104 !plus 1/3 becaues every triple is three times in each unitcell 3105 ktau=ktaux*lat(:,1)+ktauy*lat(:,2)+ktauz*lat(:,3) 3106 dumvec=ktau 3107 dumvec=dumvec*dumvec 3108 rik2=SUM(dumvec) 3109 if (rik2.gt.abcthr)cycle 3110 rr0ik=DSQRT(rik2)/r0ik 3111 3112 dumvec=jkvec+ktau-jtau 3113 dumvec=dumvec*dumvec 3114 rjk2=SUM(dumvec) 3115 if (rjk2.gt.abcthr)cycle 3116 rr0jk=DSQRT(rjk2)/r0jk 3117 3118 3119 geomean=(rr0ij*rr0ik*rr0jk)**(1./3.) 3120 fdamp=1./(1.+6.*(sr9*geomean)**alp9) 3121 tmp1 = (rij2+rjk2-rik2) 3122 tmp2 = (rij2+rik2-rjk2) 3123 tmp3 = (rik2+rjk2-rij2) 3124 tmp4=rij2*rjk2*rik2 3125 ang=(0.375d0*tmp1*tmp2*tmp3/tmp4+1.0d0)/tmp4**1.5d0 3126 3127 eabc=eabc+c9*fdamp*ang/6.0d0 3128 3129 end do 3130 end do 3131 end do 3132 end do 3133 end do 3134 end do 3135 3136 3137 end do 3138 3139 end subroutine pbcthreebody 3140 3141 3142 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3143 ! compute gradient 3144 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3145 3146 subroutine pbcgdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& 3147 & rcov,s6,s18,rs6,rs8,rs10,alp6,alp8,alp10,noabc,num,& 3148 & version,g,disp,gnorm,stress,lat,rep_v,rep_cn,& 3149 & crit_vdw,echo,crit_cn) 3150 3151 integer n,iz(*),max_elem,maxc,version,mxc(max_elem) 3152 real(wp) xyz(3,*),r0ab(max_elem,max_elem),r2r4(*) 3153 real(wp) c6ab(max_elem,max_elem,maxc,maxc,3) 3154 real(wp) g(3,*),s6,s18,rcov(max_elem) 3155 real(wp) rs6,rs8,rs10,alp10,alp8,alp6 3156 real(wp) a1,a2 3157 logical noabc,num,echo 3158 ! coversion factors 3159 3160 integer iat,jat,i,j,kat,my,ny,a,b,idum 3161 real(wp) R0,C6,R42,disp,x1,e6abc 3162 real(wp) r2,r,r4,r6,r8,t6,t8,damp1 3163 real(wp) damp6,damp8,damp9,e6,e8,e10,e12,gnorm,tmp1 3164 real(wp) s10,s8,term,step,dispr,displ,r235,tmp2 3165 real(wp) cn(n),rthr 3166 real(wp), DIMENSION(3,3) :: lat,stress,sigma,virialstress,lat_1 3167 integer, DIMENSION(3) :: rep_v,rep_cn 3168 real(wp) crit_vdw,crit_cn 3169 integer taux,tauy,tauz 3170 real(wp), DIMENSION(3) :: tau,dxyz 3171 3172 real(wp) rij(3),r7,r9 3173 real(wp) rcovij 3174 real(wp) expterm 3175 real(wp), allocatable,dimension(:,:,:,:) :: drij 3176 real(wp) dcnn 3177 real(wp) :: dc6_rest 3178 real(wp) vec(3),vec2(3) 3179 real(wp) dc6i(n) 3180 real(wp) dc6ij(n,n) 3181 real(wp) dc6_rest_sum(n*(n+1)/2) 3182 integer linij,linik,linjk 3183 real(wp) abc(3,n) 3184 3185 real(wp) eabc 3186 real(wp) gabc(3,n),glatabc(3,3) 3187 real(wp) sigma_abc(3,3) 3188 real(wp) labc,rabc 3189 real(wp) ,dimension(3) ::ijvec,ikvec,jkvec,jtau,ktau,dumvec 3190 integer jtaux,jtauy,jtauz,ktaux,ktauy,ktauz 3191 real(wp) rij2,rik2,rjk2,c9,c6ij,c6ik,c6jk,geomean,geomean3 3192 real(wp) rr0ij,rr0jk,rr0ik,dc6iji,dc6ijj 3193 real(wp) :: sr9=0.75d0 3194 real(wp), parameter :: alp9=-16.0d0 3195 real(wp),DIMENSION(n*(n+1)) ::c6save 3196 real(wp) abcthr,time1,time2,geomean2,r0av,dc9,dfdmp,dang,ang 3197 integer,dimension(3) ::repv,repmin,repmax 3198 3199 ! R^2 cut-off 3200 rthr=crit_vdw 3201 abcthr=crit_cn 3202 ! write(*,*)'abcthr:', abcthr**(1./1.) 3203 sigma=0.0d0 3204 virialstress=0.0d0 3205 stress=0.0d0 3206 gabc=0.0d0 3207 glatabc=0.0d0 3208 3209 ! testsum=0.0d0 3210 3211 if (echo)write(*,*) 3212 3213 if (num) then 3214 if (echo) & 3215 & write(*,*) 'doing numerical gradient O(N^3) ...' 3216 3217 call pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& 3218 & rcov,rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,& 3219 & e6,e8,e10,e12,e6abc,lat,rthr,rep_v,crit_cn,rep_cn) 3220 3221 3222 disp=-s6*e6-s18*e8-e6abc 3223 3224 step=2.d-5 3225 3226 do i=1,n 3227 do j=1,3 3228 xyz(j,i)=xyz(j,i)+step 3229 call pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& 3230 & rcov,rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,& 3231 & e6,e8,e10,e12,e6abc,lat,rthr,rep_v,crit_cn,rep_cn) 3232 3233 dispr=-s6*e6-s18*e8-e6abc 3234 rabc=e6abc 3235 xyz(j,i)=xyz(j,i)-2*step 3236 call pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& 3237 & rcov,rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,& 3238 & e6,e8,e10,e12,e6abc,lat,rthr,rep_v,crit_cn,rep_cn) 3239 3240 displ=-s6*e6-s18*e8-e6abc 3241 labc=e6abc 3242 gabc(j,i)=0.5*(rabc-labc)/step 3243 g(j,i)=0.5*(dispr-displ)/step 3244 xyz(j,i)=xyz(j,i)+step 3245 end do 3246 end do 3247 if (echo) write(*,*)'Doing numerical stresstensor...' 3248 3249 call xyz_to_abc(xyz,abc,lat,n) 3250 step=2.d-5 3251 if (echo) write(*,*)'step: ',step 3252 do i=1,3 3253 do j=1,3 3254 lat(j,i)=lat(j,i)+step 3255 call abc_to_xyz(abc,xyz,lat,n) 3256 call pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& 3257 & rcov,rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,& 3258 & e6,e8,e10,e12,e6abc,lat,rthr,rep_v,crit_cn,rep_cn) 3259 3260 dispr=-s6*e6-s18*e8-e6abc 3261 labc=e6abc 3262 3263 3264 lat(j,i)=lat(j,i)-2*step 3265 call abc_to_xyz(abc,xyz,lat,n) 3266 call pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& 3267 & rcov,rs6,rs8,rs10,alp6,alp8,alp10,version,noabc,& 3268 & e6,e8,e10,e12,e6abc,lat,rthr,rep_v,crit_cn,rep_cn) 3269 3270 displ=-s6*e6-s18*e8-e6abc 3271 rabc=e6abc 3272 stress(j,i)=(dispr-displ)/(step*2.0) 3273 glatabc(j,i)=(rabc-labc)/(step*2.0d0) 3274 3275 lat(j,i)=lat(j,i)+step 3276 call abc_to_xyz(abc,xyz,lat,n) 3277 3278 end do 3279 end do 3280 3281 sigma=0.0d0 3282 call inv_cell(lat,lat_1) 3283 do a=1,3 3284 do b=1,3 3285 do my=1,3 3286 sigma(a,b)=sigma(a,b)-stress(a,my)*lat(b,my) 3287 end do 3288 end do 3289 end do 3290 3291 goto 999 3292 3293 end if 3294 3295 3296 if (version.eq.2)then 3297 if (echo)write(*,*) 'doing analytical gradient D-old O(N^2) ...' 3298 disp=0 3299 stress=0.0d0 3300 do iat=1,n-1 3301 do jat=iat+1,n 3302 R0=r0ab(iz(jat),iz(iat))*rs6 3303 c6=c6ab(iz(jat),iz(iat),1,1,1)*s6 3304 do taux=-rep_v(1),rep_v(1) 3305 do tauy=-rep_v(2),rep_v(2) 3306 do tauz=-rep_v(3),rep_v(3) 3307 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 3308 dxyz=xyz(:,iat)-xyz(:,jat)+tau 3309 r2 =sum(dxyz*dxyz) 3310 if (r2.gt.rthr) cycle 3311 r235=r2**3.5 3312 r =dsqrt(r2) 3313 damp6=exp(-alp6*(r/R0-1.0d0)) 3314 damp1=1.+damp6 3315 tmp1=damp6/(damp1*damp1*r235*R0) 3316 tmp2=6./(damp1*r*r235) 3317 3318 term=alp6*tmp1-tmp2 3319 g(:,iat)=g(:,iat)-term*dxyz*c6 3320 g(:,jat)=g(:,jat)+term*dxyz*c6 3321 disp=disp+c6*(1./damp1)/r2**3 3322 3323 do ny=1,3 3324 do my=1,3 3325 sigma(my,ny)=sigma(my,ny)+term*dxyz(ny)*dxyz(my)*c6 3326 end do 3327 end do 3328 end do 3329 end do 3330 end do 3331 end do 3332 end do 3333 ! and now the self interaction, only for convenient energy in dispersion 3334 do iat=1,n 3335 jat=iat 3336 R0=r0ab(iz(jat),iz(iat))*rs6 3337 c6=c6ab(iz(jat),iz(iat),1,1,1)*s6 3338 do taux=-rep_v(1),rep_v(1) 3339 do tauy=-rep_v(2),rep_v(2) 3340 do tauz=-rep_v(3),rep_v(3) 3341 if (taux.eq.0 .and. tauy.eq.0 .and. tauz.eq.0) cycle 3342 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 3343 3344 dxyz=tau 3345 ! vec12=(/ dx,dy,dz /) 3346 r2 =sum(dxyz*dxyz) 3347 if (r2.gt.rthr) cycle 3348 r235=r2**3.5 3349 r =dsqrt(r2) 3350 damp6=exp(-alp6*(r/R0-1.0d0)) 3351 damp1=1.+damp6 3352 tmp1=damp6/(damp1*damp1*r235*R0) 3353 tmp2=6./(damp1*r*r235) 3354 disp=disp+(c6*(1./damp1)/r2**3)*0.50d0 3355 term=alp6*tmp1-tmp2 3356 do ny=1,3 3357 do my=1,3 3358 sigma(my,ny)=sigma(my,ny)+term*dxyz(ny)*dxyz(my)*c6*0.5d0 3359 end do 3360 end do 3361 3362 3363 end do 3364 end do 3365 end do 3366 end do 3367 3368 call inv_cell(lat,lat_1) 3369 do a=1,3 3370 do b=1,3 3371 do my=1,3 3372 stress(a,b)=stress(a,b)-sigma(a,my)*lat_1(b,my) 3373 end do 3374 end do 3375 end do 3376 3377 disp=-disp 3378 ! sigma=virialstress 3379 goto 999 3380 end if 3381 3382 if ((version.eq.3).or.(version.eq.5)) then 3383!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3384 ! 3385 ! begin ZERO DAMPING GRADIENT 3386 ! 3387!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3388 3389 if (echo)& 3390 & write(*,*) 'doing analytical gradient O(N^2) ...' 3391 ! precompute for analytical part 3392 call pbcncoord(n,rcov,iz,xyz,cn,lat,rep_cn,crit_cn) 3393 3394 3395 s8 =s18 3396 s10=s18 3397 allocate(drij(-rep_v(3):rep_v(3),-rep_v(2):rep_v(2),& 3398 & -rep_v(1):rep_v(1),n*(n+1)/2)) 3399 3400 disp=0 3401 3402 drij=0.0d0 3403 dc6_rest=0.0d0 3404 dc6_rest_sum=0.0d0 3405 c6save=0.0d0 3406 kat=0 3407 dc6i=0.0d0 3408 3409 3410 do iat=1,n 3411 call get_dC6_dCNij(maxc,max_elem,c6ab,mxc(iz(iat)),& 3412 & mxc(iz(iat)),cn(iat),cn(iat),iz(iat),iz(iat),iat,iat,& 3413 & c6,dc6iji,dc6ijj) 3414 3415 c6save(lin(iat,iat))=c6 3416 dc6ij(iat,iat)=dc6iji 3417 r0=r0ab(iz(iat),iz(iat)) 3418 r42=r2r4(iz(iat))*r2r4(iz(iat)) 3419 rcovij=rcov(iz(iat))+rcov(iz(iat)) 3420 3421 3422 do taux=-rep_v(1),rep_v(1) 3423 do tauy=-rep_v(2),rep_v(2) 3424 do tauz=-rep_v(3),rep_v(3) 3425 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 3426 3427 3428 !first dE/d(tau) saved in drij(i,i,counter) 3429 rij=tau 3430 r2=sum(rij*rij) 3431 ! if (r2.gt.rthr) cycle 3432 3433 if (r2.gt.0.1.and.r2.lt.rthr) then 3434 3435 3436 r=dsqrt(r2) 3437 r6=r2*r2*r2 3438 r7=r6*r 3439 r8=r6*r2 3440 r9=r8*r 3441 3442 ! 3443 ! Calculates damping functions: 3444 if (version.eq.3) then 3445 t6 = (r/(rs6*R0))**(-alp6) 3446 damp6 =1.d0/( 1.d0+6.d0*t6 ) 3447 t8 = (r/(rs8*R0))**(-alp8) 3448 damp8 =1.d0/( 1.d0+6.d0*t8 ) 3449 3450 drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat,& 3451 & iat))& 3452 & +(-s6*(6.0/(r7)*C6*damp6)& 3453 & -s8*(24.0/(r9)*C6*r42*damp8))*0.5d0 3454 3455 3456 drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat,& 3457 & iat))& 3458 & +(s6*C6/r7*6.d0*alp6*t6*damp6*damp6& 3459 & +s8*C6*r42/r9*18.d0*alp8*t8*damp8*damp8)*0.5d0 3460 else !version.eq.5 3461 t6 = (r/(rs6*R0)+R0*rs8)**(-alp6) 3462 damp6 =1.d0/( 1.d0+6.d0*t6 ) 3463 t8 = (r/(R0)+R0*rs8)**(-alp8) 3464 damp8 =1.d0/( 1.d0+6.d0*t8 ) 3465 3466 tmp1=s6*6.d0*damp6*C6/r7 3467 tmp2=s8*6.d0*C6*r42*damp8/r9 3468 drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat, & 3469 & iat)) - (tmp1 +4.d0*tmp2)*0.5d0 ! d(r^(-6))/d(r_ij) 3470 3471 3472 drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat, & 3473 & iat)) +(tmp1*alp6*t6*damp6*r/(r+rs6*R0*R0*rs8) & 3474 & +3.d0*tmp2*alp8*t8*damp8*r/(r+R0*R0*rs8))*0.5d0 !d(f_dmp)/d(r_ij) 3475 endif 3476 ! 3477 ! in dC6_rest all terms BUT C6-term is saved for the kat-loop 3478 ! 3479 dc6_rest=& 3480 & (s6/r6*damp6+3.d0*s8*r42/r8*damp8)*0.50d0 3481 3482 3483 disp=disp-dc6_rest*c6 3484 3485 dc6i(iat)=dc6i(iat)+dc6_rest*(dc6iji+dc6ijj) 3486 ! if (r2.lt.crit_cn) 3487 dc6_rest_sum(lin(iat,iat))=dc6_rest_sum(lin(iat,iat))+dc6_rest 3488 3489 3490 else 3491 drij(tauz,tauy,taux,lin(iat,iat))=0.0d0 3492 end if 3493 3494 3495 end do 3496 end do 3497 end do 3498 3499!!!!!!!!!!!!!!!!!!!!!!!!!! 3500 ! B E G I N jat L O O P 3501!!!!!!!!!!!!!!!!!!!!!!!!!! 3502 do jat=1,iat-1 3503 ! 3504 ! get_dC6_dCNij calculates the derivative dC6(iat,jat)/dCN(iat) and 3505 ! dC6(iat,jat)/dCN(jat). these are saved in dC6ij for the kat loop 3506 ! 3507 call get_dC6_dCNij(maxc,max_elem,c6ab,mxc(iz(iat)),& 3508 & mxc(iz(jat)),cn(iat),cn(jat),iz(iat),iz(jat),iat,jat,& 3509 & c6,dc6iji,dc6ijj) 3510 3511 r0=r0ab(iz(jat),iz(iat)) 3512 r42=r2r4(iz(iat))*r2r4(iz(jat)) 3513 rcovij=rcov(iz(iat))+rcov(iz(jat)) 3514 linij=lin(iat,jat) 3515 3516 dc6ij(iat,jat)=dc6iji 3517 dc6ij(jat,iat)=dc6ijj 3518 c6save(linij)=c6 3519 do taux=-rep_v(1),rep_v(1) 3520 do tauy=-rep_v(2),rep_v(2) 3521 do tauz=-rep_v(3),rep_v(3) 3522 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 3523 3524 3525 rij=xyz(:,jat)-xyz(:,iat)+tau 3526 r2=sum(rij*rij) 3527 if (r2.gt.rthr) cycle 3528 3529 3530 r=dsqrt(r2) 3531 r6=r2*r2*r2 3532 r7=r6*r 3533 r8=r6*r2 3534 r9=r8*r 3535 3536 ! 3537 ! Calculates damping functions: 3538 if (version.eq.3) then 3539 t6 = (r/(rs6*R0))**(-alp6) 3540 damp6 =1.d0/( 1.d0+6.d0*t6 ) 3541 t8 = (r/(rs8*R0))**(-alp8) 3542 damp8 =1.d0/( 1.d0+6.d0*t8 ) 3543 3544 drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,& 3545 & linij)& 3546 & -s6*(6.0/(r7)*C6*damp6)& 3547 & -s8*(24.0/(r9)*C6*r42*damp8) 3548 3549 drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,& 3550 & linij)& 3551 & +s6*C6/r7*6.d0*alp6*t6*damp6*damp6& 3552 & +s8*C6*r42/r9*18.d0*alp8*t8*damp8*damp8 3553 else !version.eq.5 3554 t6 = (r/(rs6*R0)+R0*rs8)**(-alp6) 3555 damp6 =1.d0/( 1.d0+6.d0*t6 ) 3556 t8 = (r/(R0)+R0*rs8)**(-alp8) 3557 damp8 =1.d0/( 1.d0+6.d0*t8 ) 3558 3559 tmp1=s6*6.d0*damp6*C6/r7 3560 tmp2=s8*6.d0*C6*r42*damp8/r9 3561 drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux, & 3562 & linij) - (tmp1 +4.d0*tmp2) ! d(r^(-6))/d(r_ij) 3563 3564 3565 drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,linij) & 3566 & +(tmp1*alp6*t6*damp6*r/(r+rs6*R0*R0*rs8) & 3567 & +3.d0*tmp2*alp8*t8*damp8*r/(r+R0*R0*rs8)) !d(f_dmp)/d(r_ij) 3568 endif 3569 ! 3570 ! in dC6_rest all terms BUT C6-term is saved for the kat-loop 3571 ! 3572 dc6_rest=& 3573 & (s6/r6*damp6+3.d0*s8*r42/r8*damp8) 3574 3575 3576 disp=disp-dc6_rest*c6 3577 3578 dc6i(iat)=dc6i(iat)+dc6_rest*dc6iji 3579 dc6i(jat)=dc6i(jat)+dc6_rest*dc6ijj 3580 ! if (r2.lt.crit_cn) 3581 dc6_rest_sum(linij)=dc6_rest_sum(linij)& 3582 & +dc6_rest 3583 3584 3585 end do 3586 end do 3587 end do 3588 3589 end do 3590 3591 end do 3592 3593 elseif ((version.eq.4).or.(version.eq.6)) then 3594 3595 3596 3597!!!!!!!!!!!!!!!!!!!!!!! 3598 ! NOW THE BJ Gradient ! 3599!!!!!!!!!!!!!!!!!!!!!!! 3600 3601 3602 if (echo) write(*,*) 'doing analytical gradient O(N^2) ...' 3603 call pbcncoord(n,rcov,iz,xyz,cn,lat,rep_cn,crit_cn) 3604 3605 a1 =rs6 3606 a2 =rs8 3607 s8 =s18 3608 3609 allocate(drij(-rep_v(3):rep_v(3),-rep_v(2):rep_v(2),& 3610 & -rep_v(1):rep_v(1),n*(n+1)/2)) 3611 disp=0 3612 drij=0.0d0 3613 dc6_rest=0.0d0 3614 dc6_rest_sum=0.0d0 3615 dc6i(:) = 0.0d0 3616 kat=0 3617 3618 do iat=1,n 3619 call get_dC6_dCNij(maxc,max_elem,c6ab,mxc(iz(iat)),& 3620 & mxc(iz(iat)),cn(iat),cn(iat),iz(iat),iz(iat),iat,iat,& 3621 & c6,dc6iji,dc6ijj) 3622 3623 dc6ij(iat,iat)=dc6iji 3624 c6save(lin(iat,iat))=c6 3625 r42=r2r4(iz(iat))*r2r4(iz(iat)) 3626 rcovij=rcov(iz(iat))+rcov(iz(iat)) 3627 3628 R0=a1*sqrt(3.0d0*r42)+a2 3629 3630 do taux=-rep_v(1),rep_v(1) 3631 do tauy=-rep_v(2),rep_v(2) 3632 do tauz=-rep_v(3),rep_v(3) 3633 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 3634 3635 !first dE/d(tau) saved in drij(i,i,counter) 3636 rij=tau 3637 r2=sum(rij*rij) 3638 ! if (r2.gt.rthr) cycle 3639 3640 ! if (r2.gt.0.1) then 3641 if (r2.gt.0.1.and.r2.lt.rthr) then 3642 ! 3643 ! get_dC6_dCNij calculates the derivative dC6(iat,jat)/dCN(iat) and 3644 ! dC6(iat,jat)/dCN(jat). these are saved in dC6ij for the kat loop 3645 ! 3646 r=dsqrt(r2) 3647 r4=r2*r2 3648 r6=r4*r2 3649 r7=r6*r 3650 r8=r6*r2 3651 r9=r8*r 3652 3653 ! 3654 ! Calculates damping functions: 3655 3656 t6=(r6+R0**6) 3657 t8=(r8+R0**8) 3658 3659 drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat, & 3660 & iat))& 3661 & -s6*C6*6.0d0*r4*r/(t6*t6)*0.5d0& 3662 & -s8*C6*24.0d0*r42*r7/(t8*t8)*0.5d0 3663 3664 3665 ! 3666 ! in dC6_rest all terms BUT C6-term is saved for the kat-loop 3667 ! 3668 dc6_rest=& 3669 & (s6/t6+3.d0*s8*r42/t8)*0.50d0 3670 3671 3672 disp=disp-dc6_rest*c6 3673 3674 dc6i(iat)=dc6i(iat)+dc6_rest*(dc6iji+dc6ijj) 3675 ! if (r2.lt.crit_cn) 3676 dc6_rest_sum(lin(iat,iat))=dc6_rest_sum(lin(iat,iat))+& 3677 & dc6_rest 3678 3679 3680 else 3681 drij(tauz,tauy,taux,lin(iat,iat))=0.0d0 3682 end if 3683 3684 3685 end do 3686 end do 3687 end do 3688 3689!!!!!!!!!!!!!!!!!!!!!!!!!! 3690 ! B E G I N jat L O O P 3691!!!!!!!!!!!!!!!!!!!!!!!!!! 3692 do jat=1,iat-1 3693 ! 3694 ! get_dC6_dCNij calculates the derivative dC6(iat,jat)/dCN(iat) and 3695 ! dC6(iat,jat)/dCN(jat). these are saved in dC6ij for the kat loop 3696 ! 3697 call get_dC6_dCNij(maxc,max_elem,c6ab,mxc(iz(iat)),& 3698 & mxc(iz(jat)),cn(iat),cn(jat),iz(iat),iz(jat),iat,jat,& 3699 & c6,dc6iji,dc6ijj) 3700 3701 r42=r2r4(iz(iat))*r2r4(iz(jat)) 3702 rcovij=rcov(iz(iat))+rcov(iz(jat)) 3703 3704 R0=a1*dsqrt(3.0d0*r42)+a2 3705 3706 linij=lin(iat,jat) 3707 dc6ij(iat,jat)=dc6iji 3708 dc6ij(jat,iat)=dc6ijj 3709 c6save(linij)=c6 3710 do taux=-rep_v(1),rep_v(1) 3711 do tauy=-rep_v(2),rep_v(2) 3712 do tauz=-rep_v(3),rep_v(3) 3713 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 3714 3715 3716 rij=xyz(:,jat)-xyz(:,iat)+tau 3717 r2=sum(rij*rij) 3718 if (r2.gt.rthr) cycle 3719 3720 3721 r=dsqrt(r2) 3722 r4=r2*r2 3723 r6=r4*r2 3724 r7=r6*r 3725 r8=r6*r2 3726 r9=r8*r 3727 3728 ! 3729 ! Calculates damping functions: 3730 t6=(r6+R0**6) 3731 t8=(r8+R0**8) 3732 3733 3734 drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,& 3735 & linij)& 3736 & -s6*C6*6.0d0*r4*r/(t6*t6)& 3737 & -s8*C6*24.0d0*r42*r7/(t8*t8) 3738 3739 ! 3740 ! in dC6_rest all terms BUT C6-term is saved for the kat-loop 3741 ! 3742 dc6_rest=& 3743 & (s6/t6+3.d0*s8*r42/t8) 3744 3745 3746 disp=disp-dc6_rest*c6 3747 3748 dc6i(iat)=dc6i(iat)+dc6_rest*dc6iji 3749 dc6i(jat)=dc6i(jat)+dc6_rest*dc6ijj 3750 ! if (r2.lt.crit_cn) 3751 dc6_rest_sum(lin(iat,jat))=dc6_rest_sum(linij)& 3752 & +dc6_rest 3753 3754 3755 end do 3756 end do 3757 end do 3758 3759 end do 3760 3761 end do 3762 3763 end if 3764 3765 3766 ! 3767!!!!!!!!!!!!!!!!!!!!!!! 3768 !! BEGIN Threebody gradient 3769!!!!!!!!!!!!!!!!!!!!!!! 3770 if (.not.noabc) then 3771 3772 ! write(*,*)'!!!!!!!!!! THREEBODY GRADIENT !!!!!!!!!!' 3773 sr9=0.75d0 3774 eabc=0.0d0 3775 abcthr=crit_cn 3776 repv=rep_cn 3777 ! write(*,*)'thr:',sqrt(abcthr) 3778 3779 call cpu_time(time1) 3780 do iat=3,n 3781 do jat=2,iat-1 3782 linij=lin(iat,jat) 3783 ijvec=xyz(:,jat)-xyz(:,iat) 3784 3785 c6ij=c6save(linij) 3786 do kat=1,jat-1 3787 linik=lin(iat,kat) 3788 linjk=lin(jat,kat) 3789 ikvec=xyz(:,kat)-xyz(:,iat) 3790 jkvec=xyz(:,kat)-xyz(:,jat) 3791 3792 c6ik=c6save(linik) 3793 c6jk=c6save(linjk) 3794 c9=-1.0d0*dsqrt(c6ij*c6ik*c6jk) 3795 3796 do jtaux=-rep_cn(1),rep_cn(1) 3797 repmin(1)=max(-rep_cn(1),jtaux-rep_cn(1)) 3798 repmax(1)=min(rep_cn(1),jtaux+rep_cn(1)) 3799 do jtauy=-rep_cn(2),rep_cn(2) 3800 repmin(2)=max(-rep_cn(2),jtauy-rep_cn(2)) 3801 repmax(2)=min(rep_cn(2),jtauy+rep_cn(2)) 3802 do jtauz=-rep_cn(3),rep_cn(3) 3803 repmin(3)=max(-rep_cn(3),jtauz-rep_cn(3)) 3804 repmax(3)=min(rep_cn(3),jtauz+rep_cn(3)) 3805 jtau=jtaux*lat(:,1)+jtauy*lat(:,2)+jtauz*lat(:,3) 3806 rij2=SUM((ijvec+jtau)*(ijvec+jtau)) 3807 if (rij2.gt.abcthr)cycle 3808 3809 rr0ij=DSQRT(rij2)/r0ab(iz(iat),iz(jat)) 3810 3811 3812 do ktaux=repmin(1),repmax(1) 3813 do ktauy=repmin(2),repmax(2) 3814 do ktauz=repmin(3),repmax(3) 3815 ktau=ktaux*lat(:,1)+ktauy*lat(:,2)+ktauz*lat(:,3) 3816 rik2=SUM((ikvec+ktau)*(ikvec+ktau)) 3817 if (rik2.gt.abcthr)cycle 3818 3819 dumvec=jkvec+ktau-jtau 3820 rjk2=SUM(dumvec*dumvec) 3821 if (rjk2.gt.abcthr)cycle 3822 rr0ik=dsqrt(rik2)/r0ab(iz(iat),iz(kat)) 3823 rr0jk=dsqrt(rjk2)/r0ab(iz(jat),iz(kat)) 3824 geomean2=(rij2*rjk2*rik2) 3825 ! first calculate the three components for the energy calculation fdmp 3826 ! and ang 3827 r0av=(rr0ij*rr0ik*rr0jk)**(1.0d0/3.0d0) 3828 damp9=1./(1.+6.*(sr9*r0av)**alp9) 3829 3830 geomean=dsqrt(geomean2) 3831 geomean3=geomean*geomean2 3832 ang=0.375d0*(rij2+rjk2-rik2)*(rij2-rjk2+rik2)& 3833 & *(-rij2+rjk2+rik2)/(geomean3*geomean2)& 3834 & +1.0d0/(geomean3) 3835 3836 dc6_rest=ang*damp9 3837 eabc=eabc+dc6_rest*c9 3838 ! 3839 !start calculating the gradient components dfdmp, dang and dc9 3840 3841 !dfdmp is the same for all three distances 3842 dfdmp=2.d0*alp9*(0.75d0*r0av)**(alp9)*damp9*damp9 3843 3844 !start calculating the derivatives of each part w.r.t. r_ij 3845 r=dsqrt(rij2) 3846 3847 3848 dang=-0.375d0*(rij2**3+rij2**2*(rjk2+rik2)& 3849 & +rij2*(3.0d0*rjk2**2+2.0*rjk2*rik2+3.0*rik2**2)& 3850 & -5.0*(rjk2-rik2)**2*(rjk2+rik2))& 3851 & /(r*geomean3*geomean2) 3852 3853 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 3854 drij(jtauz,jtauy,jtaux,linij)=& 3855 & drij(jtauz,jtauy,jtaux,linij)-tmp1 3856 3857 !start calculating the derivatives of each part w.r.t. r_ik 3858 3859 r=dsqrt(rik2) 3860 3861 3862 dang=-0.375d0*(rik2**3+rik2**2*(rjk2+rij2)& 3863 & +rik2*(3.0d0*rjk2**2+2.0*rjk2*rij2+3.0*rij2**2)& 3864 & -5.0*(rjk2-rij2)**2*(rjk2+rij2))& 3865 & /(r*geomean3*geomean2) 3866 3867 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 3868 ! tmp1=-dc9 3869 drij(ktauz,ktauy,ktaux,linik)=& 3870 & drij(ktauz,ktauy,ktaux,linik)-tmp1 3871 3872 ! 3873 !start calculating the derivatives of each part w.r.t. r_jk 3874 3875 r=dsqrt(rjk2) 3876 3877 dang=-0.375d0*(rjk2**3+rjk2**2*(rik2+rij2)& 3878 & +rjk2*(3.0d0*rik2**2+2.0*rik2*rij2+3.0*rij2**2)& 3879 & -5.0*(rik2-rij2)**2*(rik2+rij2))& 3880 & /(r*geomean3*geomean2) 3881 3882 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 3883 drij(ktauz-jtauz,ktauy-jtauy,ktaux-jtaux,linjk)=& 3884 & drij(ktauz-jtauz,ktauy-jtauy,ktaux-jtaux,linjk)-tmp1 3885 3886 !calculating the CN derivative dE_disp(ijk)/dCN(i) 3887 3888 dc9=dc6ij(iat,jat)/c6ij+dc6ij(iat,kat)/c6ik 3889 dc9=0.5d0*c9*dc9 3890 dc6i(iat)=dc6i(iat)+dc6_rest*dc9 3891 3892 dc9=dc6ij(jat,iat)/c6ij+dc6ij(jat,kat)/c6jk 3893 dc9=0.5d0*c9*dc9 3894 dc6i(jat)=dc6i(jat)+dc6_rest*dc9 3895 3896 dc9=dc6ij(kat,iat)/c6ik+dc6ij(kat,jat)/c6jk 3897 dc9=0.5d0*c9*dc9 3898 dc6i(kat)=dc6i(kat)+dc6_rest*dc9 3899 3900 3901 end do 3902 end do 3903 end do 3904 end do 3905 end do 3906 end do 3907 end do 3908 end do 3909 end do 3910 3911 ! Now the interaction with jat=iat of the triples iat,iat,kat 3912 do iat=2,n 3913 jat=iat 3914 linij=lin(iat,jat) 3915 ijvec=0.0d0 3916 3917 c6ij=c6save(linij) 3918 do kat=1,iat-1 3919 linjk=lin(jat,kat) 3920 linik=linjk 3921 3922 c6ik=c6save(linik) 3923 c6jk=c6ik 3924 ikvec=xyz(:,kat)-xyz(:,iat) 3925 jkvec=ikvec 3926 c9=-dsqrt(c6ij*c6ik*c6jk) 3927 do jtaux=-repv(1),repv(1) 3928 repmin(1)=max(-repv(1),jtaux-repv(1)) 3929 repmax(1)=min(repv(1),jtaux+repv(1)) 3930 do jtauy=-repv(2),repv(2) 3931 repmin(2)=max(-repv(2),jtauy-repv(2)) 3932 repmax(2)=min(repv(2),jtauy+repv(2)) 3933 do jtauz=-repv(3),repv(3) 3934 repmin(3)=max(-repv(3),jtauz-repv(3)) 3935 repmax(3)=min(repv(3),jtauz+repv(3)) 3936 if (jtaux.eq.0 .and. jtauy.eq.0 .and. jtauz.eq.0) cycle 3937 jtau=jtaux*lat(:,1)+jtauy*lat(:,2)+jtauz*lat(:,3) 3938 dumvec=jtau 3939 rij2=SUM(dumvec*dumvec) 3940 if (rij2.gt.abcthr)cycle 3941 3942 rr0ij=DSQRT(rij2)/r0ab(iz(iat),iz(jat)) 3943 3944 do ktaux=repmin(1),repmax(1) 3945 do ktauy=repmin(2),repmax(2) 3946 do ktauz=repmin(3),repmax(3) 3947 ! every result * 0.5 3948 3949 ktau=ktaux*lat(:,1)+ktauy*lat(:,2)+ktauz*lat(:,3) 3950 dumvec=ikvec+ktau 3951 dumvec=dumvec*dumvec 3952 rik2=SUM(dumvec) 3953 if (rik2.gt.abcthr)cycle 3954 3955 dumvec=jkvec+ktau-jtau 3956 dumvec=dumvec*dumvec 3957 rjk2=SUM(dumvec) 3958 if (rjk2.gt.abcthr)cycle 3959 rr0ik=DSQRT(rik2)/r0ab(iz(iat),iz(kat)) 3960 rr0jk=DSQRT(rjk2)/r0ab(iz(jat),iz(kat)) 3961 3962 3963 geomean2=(rij2*rjk2*rik2) 3964 r0av=(rr0ij*rr0ik*rr0jk)**(1.0d0/3.0d0) 3965 damp9=1./(1.+6.*(sr9*r0av)**alp9) 3966 3967 geomean=dsqrt(geomean2) 3968 geomean3=geomean*geomean2 3969 ang=0.375d0*(rij2+rjk2-rik2)*(rij2-rjk2+rik2)& 3970 & *(-rij2+rjk2+rik2)/(geomean3*geomean2)& 3971 & +1.0d0/(geomean3) 3972 3973 3974 dc6_rest=ang*damp9/2.0d0 3975 eabc=eabc+dc6_rest*c9 3976 3977 ! iat=jat 3978 dfdmp=2.d0*alp9*(0.75d0*r0av)**(alp9)*damp9*damp9 3979 3980 !start calculating the derivatives of each part w.r.t. r_ij 3981 r=dsqrt(rij2) 3982 3983 dang=-0.375d0*(rij2**3+rij2**2*(rjk2+rik2) & 3984 & +rij2*(3.0d0*rjk2**2+2.0*rjk2*rik2+3.0*rik2**2)& 3985 & -5.0*(rjk2-rik2)**2*(rjk2+rik2))& 3986 & /(r*geomean3*geomean2) 3987 3988 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 3989 drij(jtauz,jtauy,jtaux,linij)=& 3990 & drij(jtauz,jtauy,jtaux,linij)-tmp1/2.0 3991 3992 !start calculating the derivatives of each part w.r.t. r_ik 3993 r=dsqrt(rik2) 3994 3995 3996 dang=-0.375d0*(rik2**3+rik2**2*(rjk2+rij2)& 3997 & +rik2*(3.0d0*rjk2**2+2.0*rjk2*rij2+3.0*rij2**2)& 3998 & -5.0*(rjk2-rij2)**2*(rjk2+rij2))& 3999 & /(r*geomean3*geomean2) 4000 4001 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 4002 drij(ktauz,ktauy,ktaux,linik)=& 4003 & drij(ktauz,ktauy,ktaux,linik)-tmp1/2.0 4004 ! 4005 !start calculating the derivatives of each part w.r.t. r_ik 4006 r=dsqrt(rjk2) 4007 4008 dang=-0.375d0*(rjk2**3+rjk2**2*(rik2+rij2)& 4009 & +rjk2*(3.0d0*rik2**2+2.0*rik2*rij2+3.0*rij2**2)& 4010 & -5.0*(rik2-rij2)**2*(rik2+rij2))& 4011 & /(r*geomean3*geomean2) 4012 4013 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 4014 4015 drij(ktauz-jtauz,ktauy-jtauy,ktaux-jtaux,linjk)=& 4016 & drij(ktauz-jtauz,ktauy-jtauy,ktaux-jtaux,linjk)-tmp1/2.0 4017 4018 dc9=dc6ij(iat,jat)/c6ij+dc6ij(iat,kat)/c6ik 4019 dc9=0.5d0*c9*dc9 4020 dc6i(iat)=dc6i(iat)+dc6_rest*dc9 4021 4022 dc9=dc6ij(jat,iat)/c6ij+dc6ij(jat,kat)/c6jk 4023 dc9=0.5d0*c9*dc9 4024 dc6i(jat)=dc6i(jat)+dc6_rest*dc9 4025 4026 dc9=dc6ij(kat,iat)/c6ik+dc6ij(kat,jat)/c6jk 4027 dc9=0.5d0*c9*dc9 4028 dc6i(kat)=dc6i(kat)+dc6_rest*dc9 4029 4030 4031 4032 4033 end do 4034 end do 4035 end do 4036 4037 end do 4038 end do 4039 end do 4040 end do 4041 end do 4042 4043 do iat=2,n 4044 do jat=1,iat-1 4045 kat=jat 4046 linij=lin(iat,jat) 4047 linjk=lin(jat,kat) 4048 linik=linij 4049 4050 c6ij=c6save(linij) 4051 c6ik=c6ij 4052 4053 c6jk=c6save(linjk) 4054 ikvec=xyz(:,kat)-xyz(:,iat) 4055 ijvec=ikvec 4056 jkvec=0.0d0 4057 4058 c9=-1.0d0*dsqrt(c6ij*c6ik*c6jk) 4059 do jtaux=-repv(1),repv(1) 4060 repmin(1)=max(-repv(1),jtaux-repv(1)) 4061 repmax(1)=min(repv(1),jtaux+repv(1)) 4062 do jtauy=-repv(2),repv(2) 4063 repmin(2)=max(-repv(2),jtauy-repv(2)) 4064 repmax(2)=min(repv(2),jtauy+repv(2)) 4065 do jtauz=-repv(3),repv(3) 4066 repmin(3)=max(-repv(3),jtauz-repv(3)) 4067 repmax(3)=min(repv(3),jtauz+repv(3)) 4068 4069 jtau=jtaux*lat(:,1)+jtauy*lat(:,2)+jtauz*lat(:,3) 4070 dumvec=ijvec+jtau 4071 dumvec=dumvec*dumvec 4072 rij2=SUM(dumvec) 4073 if (rij2.gt.abcthr)cycle 4074 4075 rr0ij=SQRT(rij2)/r0ab(iz(iat),iz(jat)) 4076 4077 do ktaux=repmin(1),repmax(1) 4078 do ktauy=repmin(2),repmax(2) 4079 do ktauz=repmin(3),repmax(3) 4080 ! every result * 0.5 4081 if (jtaux.eq.ktaux .and. jtauy.eq.ktauy& 4082 & .and. jtauz.eq.ktauz) cycle 4083 ktau=ktaux*lat(:,1)+ktauy*lat(:,2)+ktauz*lat(:,3) 4084 dumvec=ikvec+ktau 4085 dumvec=dumvec*dumvec 4086 rik2=SUM(dumvec) 4087 if (rik2.gt.abcthr)cycle 4088 rr0ik=SQRT(rik2)/r0ab(iz(iat),iz(kat)) 4089 4090 dumvec=jkvec+ktau-jtau 4091 dumvec=dumvec*dumvec 4092 rjk2=SUM(dumvec) 4093 if (rjk2.gt.abcthr)cycle 4094 rr0jk=SQRT(rjk2)/r0ab(iz(jat),iz(kat)) 4095 4096 ! if (rij*rjk*rik.gt.abcthr)cycle 4097 4098 geomean2=(rij2*rjk2*rik2) 4099 r0av=(rr0ij*rr0ik*rr0jk)**(1.0d0/3.0d0) 4100 damp9=1./(1.+6.d0*(sr9*r0av)**alp9) 4101 4102 geomean=dsqrt(geomean2) 4103 geomean3=geomean*geomean2 4104 ang=0.375d0*(rij2+rjk2-rik2)*(rij2-rjk2+rik2)& 4105 & *(-rij2+rjk2+rik2)/(geomean2*geomean3)& 4106 & +1.0d0/(geomean3) 4107 dc6_rest=ang*damp9/2.0d0 4108 eabc=eabc+dc6_rest*c9 4109 4110 4111 ! jat=kat 4112 dfdmp=2.d0*alp9*(0.75d0*r0av)**(alp9)*damp9*damp9 4113 !start calculating the derivatives of each part w.r.t. r_ij 4114 r=dsqrt(rij2) 4115 4116 dang=-0.375d0*(rij2**3+rij2**2*(rjk2+rik2)& 4117 & +rij2*(3.0d0*rjk2**2+2.0d0*rjk2*rik2+3.0d0*rik2**2)& 4118 & -5.0d0*(rjk2-rik2)**2*(rjk2+rik2))& 4119 & /(r*geomean3*geomean2) 4120 4121 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 4122 drij(jtauz,jtauy,jtaux,linij)=& 4123 & drij(jtauz,jtauy,jtaux,linij)-tmp1/2.0d0 4124 4125 !start calculating the derivatives of each part w.r.t. r_ik 4126 r=dsqrt(rik2) 4127 4128 4129 dang=-0.375d0*(rik2**3+rik2**2*(rjk2+rij2)& 4130 & +rik2*(3.0d0*rjk2**2+2.0*rjk2*rij2+3.0*rij2**2)& 4131 & -5.0*(rjk2-rij2)**2*(rjk2+rij2))& 4132 & /(r*geomean3*geomean2) 4133 4134 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 4135 ! tmp1=-dc9 4136 drij(ktauz,ktauy,ktaux,linik)=& 4137 & drij(ktauz,ktauy,ktaux,linik)-tmp1/2.0d0 4138 ! 4139 !start calculating the derivatives of each part w.r.t. r_jk 4140 r=dsqrt(rjk2) 4141 4142 dang=-0.375d0*(rjk2**3+rjk2**2*(rik2+rij2)& 4143 & +rjk2*(3.0d0*rik2**2+2.0*rik2*rij2+3.0*rij2**2)& 4144 & -5.0d0*(rik2-rij2)**2*(rik2+rij2))& 4145 & /(r*geomean3*geomean2) 4146 4147 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 4148 drij(ktauz-jtauz,ktauy-jtauy,ktaux-jtaux,linjk)=& 4149 & drij(ktauz-jtauz,ktauy-jtauy,ktaux-jtaux,linjk)-tmp1/2.0d0 4150 4151 !calculating the CN derivative dE_disp(ijk)/dCN(i) 4152 4153 dc9=dc6ij(iat,jat)/c6ij+dc6ij(iat,kat)/c6ik 4154 dc9=0.5d0*c9*dc9 4155 dc6i(iat)=dc6i(iat)+dc6_rest*dc9 4156 4157 dc9=dc6ij(jat,iat)/c6ij+dc6ij(jat,kat)/c6jk 4158 dc9=0.5d0*c9*dc9 4159 dc6i(jat)=dc6i(jat)+dc6_rest*dc9 4160 4161 dc9=dc6ij(kat,iat)/c6ik+dc6ij(kat,jat)/c6jk 4162 dc9=0.5d0*c9*dc9 4163 dc6i(kat)=dc6i(kat)+dc6_rest*dc9 4164 4165 4166 4167 4168 end do 4169 end do 4170 end do 4171 4172 end do 4173 end do 4174 end do 4175 end do 4176 end do 4177 4178 4179 ! And finally the self interaction iat=jat=kat all 4180 4181 idum=0 4182 do iat=1,n 4183 jat=iat 4184 kat=iat 4185 ijvec=0.0d0 4186 linij=lin(iat,jat) 4187 linik=lin(iat,kat) 4188 linjk=lin(jat,kat) 4189 ikvec=ijvec 4190 jkvec=ikvec 4191 c6ij=c6save(linij) 4192 c6ik=c6ij 4193 c6jk=c6ij 4194 c9=-(DSQRT(c6ij*c6ij*c6ij)) 4195 4196 do jtaux=-repv(1),repv(1) 4197 repmin(1)=max(-repv(1),jtaux-repv(1)) 4198 repmax(1)=min(repv(1),jtaux+repv(1)) 4199 do jtauy=-repv(2),repv(2) 4200 repmin(2)=max(-repv(2),jtauy-repv(2)) 4201 repmax(2)=min(repv(2),jtauy+repv(2)) 4202 do jtauz=-repv(3),repv(3) 4203 repmin(3)=max(-repv(3),jtauz-repv(3)) 4204 repmax(3)=min(repv(3),jtauz+repv(3)) 4205 if ((jtaux.eq.0) .and.(jtauy.eq.0) .and.(jtauz.eq.0))cycle 4206 jtau=jtaux*lat(:,1)+jtauy*lat(:,2)+jtauz*lat(:,3) 4207 dumvec=jtau 4208 dumvec=dumvec*dumvec 4209 rij2=SUM(dumvec) 4210 if (rij2.gt.abcthr)cycle 4211 rr0ij=SQRT(rij2)/r0ab(iz(iat),iz(jat)) 4212 4213 do ktaux=repmin(1),repmax(1) 4214 do ktauy=repmin(2),repmax(2) 4215 do ktauz=repmin(3),repmax(3) 4216 if ((ktaux.eq.0) .and.( ktauy.eq.0) .and.( ktauz.eq.0))cycle 4217 if ((ktaux.eq.jtaux) .and. (ktauy.eq.jtauy)& 4218 & .and. (ktauz.eq.jtauz)) cycle 4219 4220 ! every result * 1/6 becaues every triple is counted twice due to the tw 4221 ! 4222 !plus 1/3 becaues every triple is three times in each unitcell 4223 ktau=ktaux*lat(:,1)+ktauy*lat(:,2)+ktauz*lat(:,3) 4224 dumvec=ktau 4225 dumvec=dumvec*dumvec 4226 rik2=SUM(dumvec) 4227 if (rik2.gt.abcthr)cycle 4228 rr0ik=SQRT(rik2)/r0ab(iz(iat),iz(kat)) 4229 4230 dumvec=jkvec+ktau-jtau 4231 dumvec=dumvec*dumvec 4232 rjk2=SUM(dumvec) 4233 if (rjk2.gt.abcthr)cycle 4234 rr0jk=SQRT(rjk2)/r0ab(iz(jat),iz(kat)) 4235 4236 geomean2=(rij2*rjk2*rik2) 4237 r0av=(rr0ij*rr0ik*rr0jk)**(1.0d0/3.0d0) 4238 damp9=1./(1.+6.*(sr9*r0av)**alp9) 4239 4240 geomean=dsqrt(geomean2) 4241 geomean3=geomean*geomean2 4242 ang=0.375d0*(rij2+rjk2-rik2)*(rij2-rjk2+rik2)& 4243 & *(-rij2+rjk2+rik2)/(geomean2*geomean3)& 4244 & +1.0d0/(geomean3) 4245 dc6_rest=ang*damp9/6.0d0 4246 eabc=eabc+c9*dc6_rest 4247 4248 ! iat=jat=kat 4249 dfdmp=2.d0*alp9*(0.75d0*r0av)**(alp9)*damp9*damp9 4250 !start calculating the derivatives of each part w.r.t. r_ij 4251 4252 r=dsqrt(rij2) 4253 dang=-0.375d0*(rij2**3+rij2**2*(rjk2+rik2)& 4254 & +rij2*(3.0d0*rjk2**2+2.0*rjk2*rik2+3.0*rik2**2)& 4255 & -5.0*(rjk2-rik2)**2*(rjk2+rik2))& 4256 & /(r*geomean3*geomean2) 4257 4258 4259 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 4260 drij(jtauz,jtauy,jtaux,linij)=& 4261 & drij(jtauz,jtauy,jtaux,linij)-tmp1/6.0d0 4262 4263 !start calculating the derivatives of each part w.r.t. r_ik 4264 4265 r=dsqrt(rik2) 4266 4267 dang=-0.375d0*(rik2**3+rik2**2*(rjk2+rij2)& 4268 & +rik2*(3.0d0*rjk2**2+2.0*rjk2*rij2+3.0*rij2**2)& 4269 & -5.0*(rjk2-rij2)**2*(rjk2+rij2))& 4270 & /(r*geomean3*geomean2) 4271 4272 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 4273 drij(ktauz,ktauy,ktaux,linik)=& 4274 & drij(ktauz,ktauy,ktaux,linik)-tmp1/6.0d0 4275 ! 4276 !start calculating the derivatives of each part w.r.t. r_jk 4277 4278 r=dsqrt(rjk2) 4279 dang=-0.375d0*(rjk2**3+rjk2**2*(rik2+rij2)& 4280 & +rjk2*(3.0d0*rik2**2+2.0*rik2*rij2+3.0*rij2**2)& 4281 & -5.0*(rik2-rij2)**2*(rik2+rij2))& 4282 & /(r*geomean3*geomean2) 4283 4284 tmp1=-dang*c9*damp9+dfdmp/r*c9*ang 4285 drij(ktauz-jtauz,ktauy-jtauy,ktaux-jtaux,linjk)=& 4286 & drij(ktauz-jtauz,ktauy-jtauy,ktaux-jtaux,linjk)-tmp1/6.0d0 4287 4288 4289 !calculating the CN derivative dE_disp(ijk)/dCN(i) 4290 4291 dc9=dc6ij(iat,jat)/c6ij+dc6ij(iat,kat)/c6ik 4292 dc9=0.5d0*c9*dc9 4293 dc6i(iat)=dc6i(iat)+dc6_rest*dc9 4294 4295 dc9=dc6ij(jat,iat)/c6ij+dc6ij(jat,kat)/c6jk 4296 dc9=0.5d0*c9*dc9 4297 dc6i(jat)=dc6i(jat)+dc6_rest*dc9 4298 4299 dc9=dc6ij(kat,iat)/c6ik+dc6ij(kat,jat)/c6jk 4300 dc9=0.5d0*c9*dc9 4301 dc6i(kat)=dc6i(kat)+dc6_rest*dc9 4302 4303 4304 4305 4306 4307 end do 4308 end do 4309 end do 4310 end do 4311 end do 4312 !jtaux 4313 end do 4314 !should exclude tabst 4315 !get type of string, 0=numb 4316 !special case: end of line 4317 !cast string on real and get er 4318 !handle exceptions 4319 !check for integer/real 4320 !if integer, add .0 to string; otherwi 4321 ! Selective dynamics 4322 ! Cartesian or direct 4323 !first line must contain Element Info 4324 !second line contains global scaling f 4325 !the Ang->au conversion is included in 4326 ! reading the lattice constants 4327 ! write(*,'(3F6.2)')lattice(1,i),lattice(2,i),lattice(3,i) 4328 !Ether here are the numbers of each el 4329 ! CONTCAR files have additional Element lin 4330 !tauz 4331 !tauy 4332 !taux 4333 !iat 4334 !i 4335 ! Selective dynamics 4336 ! Cartesian or direct 4337 !first line must contain Element Info 4338 !second line contains global scaling f 4339 ! reading the lattice constants 4340 ! write(*,'(3F6.2)')lattice(1,i),lattice(2,i),lattice(3,i) 4341 !Ether here are the numbers of each el 4342 ! CONTCAR files have additional Element lin 4343 !,r2r4(*) 4344 !,crit_vdw,crit_cn 4345 !BJ-parameter 4346 !taux 4347 !tauy 4348 !tauz 4349 !iat 4350 !tauz 4351 !tauy 4352 !taux 4353 !jat 4354 !iat 4355 !tauz 4356 !tauy 4357 !taux 4358 !iat 4359 !tauz 4360 !tauy 4361 !taux 4362 !jat 4363 !tauz 4364 !tauy 4365 !taux 4366 !iat 4367 !version 4368 !reciprocal radii scaling paramete 4369 !alpha saved with "-" sign 4370 !alp9 is already s 4371 !ktauz 4372 !ktauy 4373 !ktaux 4374 !jtauz 4375 !jtauy 4376 !jtaux 4377 !kat 4378 !jat 4379 !iat 4380 !ktauz 4381 !ktauy 4382 !ktaux 4383 !jtauz 4384 !jtauy 4385 !jtaux 4386 !kat 4387 !iat 4388 ! And now kat=jat, but cycling throug all imagecells without jtau= 4389 ! But this counts only 1/2 4390 !ktauz 4391 !ktauy 4392 !ktaux 4393 !jtauz 4394 !jtauy 4395 !jtaux 4396 !kat 4397 !iat 4398 !If kat and jat are th 4399 !ktauz 4400 !ktauy 4401 !ktaux 4402 !jtauz 4403 !jtauy 4404 !jtaux 4405 !iat 4406 !tauz 4407 !tauy 4408 !taux 4409 !j 4410 !i 4411 !BJ-parameters 4412 ! precalculated dampingterms 4413 !d(C6ij)/d(r_ij) 4414 !d(E)/d(r_ij) der 4415 !dCN(iat)/d(r_ij) 4416 !dCN(jat)/d(r_ij) 4417 !dC6i(iat) saves dE_dsp/dCN(iat) 4418 !dC6(iat,jat)/cCN(iat) in dc6ij(i,j) for ABC- 4419 !threebody gradient 4420 !inverse of 4/3 4421 !jat 4422 !iat 4423 !call edisp...dum1 4424 !call edisp...dum2 4425 !j 4426 !i 4427 !b 4428 !a 4429 !num 4430 !my 4431 !ny 4432 !tauz 4433 !tauy 4434 !taux 4435 !jat 4436 !iat 4437 !my 4438 !ny 4439 !tauz 4440 !tauy 4441 !taux 4442 !iat 4443 !b 4444 !a 4445 !version==2 4446 ! d(r^(-6))/d(tau) 4447 !d(f_dmp)/d(tau) 4448 ! calculate E_disp for sanity check 4449 !r2 < 0.1>rthr 4450 !tauz 4451 !tauy 4452 !taux 4453 ! d(r^(-6))/d(r_ij) 4454 !d(f_dmp)/d(r_ij) 4455 ! calculate E_disp for sanity check 4456 !tauz 4457 !tauy 4458 !taux 4459 !jat 4460 !iat 4461 ! d(1/(r^(6)+R0^6)/d(r) 4462 ! calculate E_disp for sanity check 4463 !r2 < 0.1>rthr 4464 !tauz 4465 !tauy 4466 !taux 4467 ! calculate E_disp for sanity check 4468 !tauz 4469 !tauy 4470 !taux 4471 !jat 4472 !iat 4473 ! version=3 or 4 4474 !alp9 is already sa 4475 !ktauz 4476 !ktauy 4477 !ktauz 4478 !jtauz 4479 !jtauy 4480 !jtaux 4481 !kat 4482 !jat 4483 !iat 4484 !alp9 is already saved 4485 !factor 1/2 for doublecounting 4486 !ktauz 4487 !ktauy 4488 !ktaux 4489 !jtauz 4490 !jtauy 4491 !jtaux 4492 !kat 4493 !iat 4494 ! And now kat=jat, but cycling throug all imagecells without jtau= 4495 ! But this counts only 1/2 4496 !alp9 is already save 4497 !factor 1/2 for doublecounting 4498 !ktauz 4499 !ktauy 4500 !ktaux 4501 !jtauz 4502 !jtauy 4503 !jtaux 4504 !kat 4505 !iat 4506 !if 4507 !If kat and jat are th 4508 !alp9 is already saved 4509 !ktauz 4510 !ktauy 4511 !ktaux 4512 !jtauz 4513 !jtauy 4514 4515 4516 end do 4517 4518 4519 call cpu_time(time2) 4520 4521 ! write(*,*)' eabc(gdisp): ',eabc 4522 ! write(*,'('' time(abc) '',f6.1)')time2-time1 4523 disp=disp-eabc 4524 ! write(*,*)'gdisp:',disp 4525 end if 4526 4527!451 continue 4528 4529 sigma_abc=0.0d0 4530 sigma=0.0d0 4531 4532 ! After calculating all derivatives dE/dr_ij w.r.t. distances, 4533 ! the grad w.r.t. the coordinates is calculated dE/dr_ij * dr_ij/dxyz_i 4534 do iat=2,n 4535 do jat=1,iat-1 4536 linij=lin(iat,jat) 4537 rcovij=rcov(iz(iat))+rcov(iz(jat)) 4538 do taux=-rep_v(1),rep_v(1) 4539 do tauy=-rep_v(2),rep_v(2) 4540 do tauz=-rep_v(3),rep_v(3) 4541 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 4542 4543 rij=xyz(:,jat)-xyz(:,iat)+tau 4544 r2=sum(rij*rij) 4545 if (r2.gt.rthr.or.r2.lt.0.5) cycle 4546 r=dsqrt(r2) 4547 4548 if (r2.lt.crit_cn) then 4549 expterm=exp(-k1*(rcovij/r-1.d0)) 4550 dcnn=-k1*rcovij*expterm/& 4551 & (r2*(expterm+1.d0)*(expterm+1.d0)) 4552 else 4553 dcnn=0.0d0 4554 end if 4555 x1=drij(tauz,tauy,taux,linij)+dcnn*(dc6i(iat)+dc6i(jat)) 4556 vec=x1*rij/r 4557 g(:,iat)=g(:,iat)+vec 4558 g(:,jat)=g(:,jat)-vec 4559 do i=1,3 4560 do j=1,3 4561 sigma(j,i)=sigma(j,i)+vec(j)*rij(i) 4562 end do 4563 end do 4564 4565 4566 4567 end do 4568 end do 4569 end do 4570 end do 4571 end do 4572 4573 do iat=1,n 4574 rcovij=rcov(iz(iat))+rcov(iz(iat)) 4575 do taux=-rep_v(1),rep_v(1) 4576 do tauy=-rep_v(2),rep_v(2) 4577 do tauz=-rep_v(3),rep_v(3) 4578 if (taux.eq.0.and.tauy.eq.0.and.tauz.eq.0) cycle 4579 4580 tau=taux*lat(:,1)+tauy*lat(:,2)+tauz*lat(:,3) 4581 r2=(sum(tau*tau)) 4582 r=dsqrt(r2) 4583 if (r2.lt.crit_cn) then 4584 expterm=exp(-k1*(rcovij/r-1.d0)) 4585 dcnn=-k1*rcovij*expterm/& 4586 & (r2*(expterm+1.d0)*(expterm+1.d0)) 4587 else 4588 dcnn=0.0d0 4589 end if 4590 x1=drij(tauz,tauy,taux,lin(iat,iat))+dcnn*dc6i(iat) 4591 vec=x1*tau/r 4592 vec2(1)=taux 4593 vec2(2)=tauy 4594 vec2(3)=tauz 4595 do i=1,3 4596 do j=1,3 4597 sigma(j,i)=sigma(j,i)+vec(j)*tau(i) 4598 end do 4599 end do 4600 4601 4602 end do 4603 end do 4604 end do 4605 4606 4607 4608 end do 4609 4610 4611 4612 stress=0.0d0 4613 glatabc=0.0d0 4614 call inv_cell(lat,lat_1) 4615 do a=1,3 4616 do b=1,3 4617 do my=1,3 4618 stress(a,b)=stress(a,b)-sigma(a,my)*lat_1(b,my) 4619 end do 4620 end do 4621 end do 4622 4623 4624 4625 ! write(*,*)'drij:',drij(lin(iat,jat),:) 4626 ! write(*,*)'g:',g(1,1:3) 4627 ! write(*,*)'dcn:',sum(dcn(lin(2,1),:)) 4628 4629 4630 4631 deallocate(drij) 4632 4633 4634 4635 4636999 continue 4637!!!!!!!!!!!!!!!!!!!!!!!!!!! 4638 ! 4639 !This is where the D2 gradient and the numerical gradient jump. 4640 ! 4641!!!!!!!!!!!!!!!!!!!!!!!!!! 4642 ! do i=1,n 4643 ! write(*,'(83F17.12)') g(1:3,i) 4644 ! end do 4645 gnorm=sum(abs(g(1:3,1:n))) 4646 if (echo)then 4647 ! write(*,*)'testsum:',testsum*autoev/autoang 4648 write(*,*)'|G(force)| =',gnorm 4649 gnorm=sum(abs(stress(1:3,1:3))) 4650 write(*,*)'|G(stress)|=',gnorm 4651 end if 4652 4653 end subroutine pbcgdisp 4654 4655 4656!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4657 4658 subroutine copyc6(fname,maxc,max_elem,c6ab,maxci, minc6,minc6list,maxc6,& 4659 & maxc6list) 4660 integer maxc,max_elem,maxci(max_elem) 4661 real(wp) c6ab(max_elem,max_elem,maxc,maxc,3) 4662 character*(*) fname 4663 logical minc6,maxc6,minc6list(max_elem),maxc6list(max_elem) 4664 4665 integer iat,jat,i,iadr,jadr,nn,kk 4666 logical special 4667 4668 call init_pars(fname) 4669 c6ab=-1 4670 maxci=0 4671 ! read file 4672 kk=1 4673 !only use values for cn=minimum 4674 if (minc6.or.maxc6) then 4675 do i=1,94 4676 if (minc6list(i))then 4677 c6ab(i,:,1,:,2)=10000000.0 4678 c6ab(:,i,:,1,3)=10000000.0 4679 end if 4680 end do 4681 4682 4683 do nn=1,nlines 4684 special=.false. 4685 iat=int(pars(kk+1)) 4686 jat=int(pars(kk+2)) 4687 call limit(iat,jat,iadr,jadr) 4688 4689 !only CN=minimum for iat 4690 if (minc6list(iat)) then 4691 special=.true. 4692 maxci(iat)=1 4693 maxci(jat)=max(maxci(jat),jadr) 4694 4695 if (pars(kk+3).le.c6ab(iat,jat,1,jadr,2)) then 4696 4697 c6ab(iat,jat,1,jadr,1)=pars(kk) 4698 c6ab(iat,jat,1,jadr,2)=pars(kk+3) 4699 c6ab(iat,jat,1,jadr,3)=pars(kk+4) 4700 4701 c6ab(jat,iat,jadr,1,1)=pars(kk) 4702 c6ab(jat,iat,jadr,1,2)=pars(kk+4) 4703 c6ab(jat,iat,jadr,1,3)=pars(kk+3) 4704 end if 4705 end if 4706 4707 !only CN=minimum for jat 4708 if (minc6list(jat)) then 4709 special=.true. 4710 maxci(iat)=max(maxci(iat),iadr) 4711 maxci(jat)=1 4712 4713 if (pars(kk+4).le.c6ab(iat,jat,iadr,1,3)) then 4714 4715 c6ab(iat,jat,iadr,1,1)=pars(kk) 4716 c6ab(iat,jat,iadr,1,2)=pars(kk+3) 4717 c6ab(iat,jat,iadr,1,3)=pars(kk+4) 4718 4719 c6ab(jat,iat,1,iadr,1)=pars(kk) 4720 c6ab(jat,iat,1,iadr,2)=pars(kk+4) 4721 c6ab(jat,iat,1,iadr,3)=pars(kk+3) 4722 end if 4723 end if 4724 4725 4726 !only CN=minimum for 4727 if (minc6list(iat).and.minc6list(jat)) then 4728 special=.true. 4729 maxci(jat)=1 4730 maxci(iat)=1 4731 4732 if (pars(kk+4).le.c6ab(iat,jat,1,1,3).and.& 4733 & pars(kk+3).le.c6ab(iat,jat,1,1,2)) then 4734 4735 c6ab(iat,jat,1,1,1)=pars(kk) 4736 c6ab(iat,jat,1,1,2)=pars(kk+3) 4737 c6ab(iat,jat,1,1,3)=pars(kk+4) 4738 4739 c6ab(jat,iat,1,1,1)=pars(kk) 4740 c6ab(jat,iat,1,1,2)=pars(kk+4) 4741 c6ab(jat,iat,1,1,3)=pars(kk+3) 4742 end if 4743 end if 4744 4745 4746 4747 !only CN=maximum for iat 4748 if (maxc6list(iat)) then 4749 special=.true. 4750 4751 maxci(iat)=1 4752 maxci(jat)=max(maxci(jat),jadr) 4753 4754 if (pars(kk+3).ge.c6ab(iat,jat,1,jadr,2)) then 4755 4756 c6ab(iat,jat,1,jadr,1)=pars(kk) 4757 c6ab(iat,jat,1,jadr,2)=pars(kk+3) 4758 c6ab(iat,jat,1,jadr,3)=pars(kk+4) 4759 4760 c6ab(jat,iat,jadr,1,1)=pars(kk) 4761 c6ab(jat,iat,jadr,1,2)=pars(kk+4) 4762 c6ab(jat,iat,jadr,1,3)=pars(kk+3) 4763 end if 4764 end if 4765 !only CN=maximum for jat 4766 if (maxc6list(jat)) then 4767 special=.true. 4768 4769 maxci(jat)=1 4770 maxci(iat)=max(maxci(iat),iadr) 4771 4772 if (pars(kk+4).ge.c6ab(iat,jat,iadr,1,3)) then 4773 4774 c6ab(iat,jat,iadr,1,1)=pars(kk) 4775 c6ab(iat,jat,iadr,1,2)=pars(kk+3) 4776 c6ab(iat,jat,iadr,1,3)=pars(kk+4) 4777 4778 c6ab(jat,iat,1,iadr,1)=pars(kk) 4779 c6ab(jat,iat,1,iadr,2)=pars(kk+4) 4780 c6ab(jat,iat,1,iadr,3)=pars(kk+3) 4781 end if 4782 end if 4783 4784 !only CN=maximum for 4785 if (maxc6list(iat).and.maxc6list(jat)) then 4786 special=.true. 4787 maxci(jat)=1 4788 maxci(iat)=1 4789 4790 if (pars(kk+4).ge.c6ab(iat,jat,1,1,3).and.& 4791 & pars(kk+3).ge.c6ab(iat,jat,1,1,2)) then 4792 4793 c6ab(iat,jat,1,1,1)=pars(kk) 4794 c6ab(iat,jat,1,1,2)=pars(kk+3) 4795 c6ab(iat,jat,1,1,3)=pars(kk+4) 4796 4797 c6ab(jat,iat,1,1,1)=pars(kk) 4798 c6ab(jat,iat,1,1,2)=pars(kk+4) 4799 c6ab(jat,iat,1,1,3)=pars(kk+3) 4800 end if 4801 end if 4802 4803 !only CN=minimum for 4804 if (minc6list(iat).and.maxc6list(jat)) then 4805 !and CN=maximum jat 4806 special=.true. 4807 maxci(jat)=1 4808 maxci(iat)=1 4809 4810 if (pars(kk+4).ge.c6ab(iat,jat,1,1,3).and.& 4811 & pars(kk+3).le.c6ab(iat,jat,1,1,2)) then 4812 4813 c6ab(iat,jat,1,1,1)=pars(kk) 4814 c6ab(iat,jat,1,1,2)=pars(kk+3) 4815 c6ab(iat,jat,1,1,3)=pars(kk+4) 4816 4817 c6ab(jat,iat,1,1,1)=pars(kk) 4818 c6ab(jat,iat,1,1,2)=pars(kk+4) 4819 c6ab(jat,iat,1,1,3)=pars(kk+3) 4820 end if 4821 end if 4822 4823 !only CN=maximum for 4824 if (maxc6list(iat).and.minc6list(jat)) then 4825 ! and CN=minimum fo 4826 special=.true. 4827 maxci(jat)=1 4828 maxci(iat)=1 4829 4830 if (pars(kk+4).le.c6ab(iat,jat,1,1,3).and.& 4831 & pars(kk+3).ge.c6ab(iat,jat,1,1,2)) then 4832 4833 c6ab(iat,jat,1,1,1)=pars(kk) 4834 c6ab(iat,jat,1,1,2)=pars(kk+3) 4835 c6ab(iat,jat,1,1,3)=pars(kk+4) 4836 4837 c6ab(jat,iat,1,1,1)=pars(kk) 4838 c6ab(jat,iat,1,1,2)=pars(kk+4) 4839 c6ab(jat,iat,1,1,3)=pars(kk+3) 4840 end if 4841 end if 4842 4843 if (.not.special) then 4844 4845 maxci(iat)=max(maxci(iat),iadr) 4846 maxci(jat)=max(maxci(jat),jadr) 4847 4848 c6ab(iat,jat,iadr,jadr,1)=pars(kk) 4849 c6ab(iat,jat,iadr,jadr,2)=pars(kk+3) 4850 c6ab(iat,jat,iadr,jadr,3)=pars(kk+4) 4851 4852 c6ab(jat,iat,jadr,iadr,1)=pars(kk) 4853 c6ab(jat,iat,jadr,iadr,2)=pars(kk+4) 4854 c6ab(jat,iat,jadr,iadr,3)=pars(kk+3) 4855 end if 4856 kk=(nn*5)+1 4857 end do 4858 4859 4860 4861 !no min/max at all 4862 else 4863 do nn=1,nlines 4864 iat=int(pars(kk+1)) 4865 jat=int(pars(kk+2)) 4866 call limit(iat,jat,iadr,jadr) 4867 maxci(iat)=max(maxci(iat),iadr) 4868 maxci(jat)=max(maxci(jat),jadr) 4869 4870 c6ab(iat,jat,iadr,jadr,1)=pars(kk) 4871 c6ab(iat,jat,iadr,jadr,2)=pars(kk+3) 4872 c6ab(iat,jat,iadr,jadr,3)=pars(kk+4) 4873 4874 c6ab(jat,iat,jadr,iadr,1)=pars(kk) 4875 c6ab(jat,iat,jadr,iadr,2)=pars(kk+4) 4876 c6ab(jat,iat,jadr,iadr,3)=pars(kk+3) 4877 kk=(nn*5)+1 4878 end do 4879 end if 4880 4881 end subroutine copyc6 4882 4883 4884 4885!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4886 4887 subroutine SET_CRITERIA(rthr,lat,tau_max) 4888 4889 REAL(WP) :: r_cutoff,rthr 4890 REAL(WP) :: lat(3,3) 4891 REAL(WP) :: tau_max(3) 4892 REAL(WP) :: norm1(3),norm2(3),norm3(3) 4893 REAL(WP) :: cos10,cos21,cos32 4894 4895 r_cutoff=sqrt(rthr) 4896 ! write(*,*) 'lat',lat 4897 call kreuzprodukt(lat(:,2),lat(:,3),norm1) 4898 call kreuzprodukt(lat(:,3),lat(:,1),norm2) 4899 call kreuzprodukt(lat(:,1),lat(:,2),norm3) 4900 ! write(*,*) 'norm2',norm2 4901 norm1=norm1/VECTORSIZE(norm1) 4902 norm2=norm2/VECTORSIZE(norm2) 4903 norm3=norm3/VECTORSIZE(norm3) 4904 ! write(*,*) 'norm2_',norm2 4905 cos10=SUM(norm1*lat(:,1)) 4906 cos21=SUM(norm2*lat(:,2)) 4907 cos32=SUM(norm3*lat(:,3)) 4908 tau_max(1)=abs(r_cutoff/cos10) 4909 tau_max(2)=abs(r_cutoff/cos21) 4910 tau_max(3)=abs(r_cutoff/cos32) 4911 ! write(*,'(3f8.4)')tau_max(1),tau_max(2),tau_max(3) 4912 end subroutine SET_CRITERIA 4913 4914 4915 subroutine kreuzprodukt(A,B,C) 4916 4917 REAL(WP) :: A(3),B(3) 4918 REAL(WP) :: X,Y,Z 4919 REAL(WP) :: C(3) 4920 4921 X=A(2)*B(3)-B(2)*A(3) 4922 Y=A(3)*B(1)-B(3)*A(1) 4923 Z=A(1)*B(2)-B(1)*A(2) 4924 C=(/X,Y,Z/) 4925 end subroutine kreuzprodukt 4926 4927 function VECTORSIZE(VECT) 4928 4929 REAL(WP) :: VECT(3) 4930 REAL(WP) :: SVECT(3) 4931 REAL(WP) :: VECTORSIZE 4932 4933 SVECT=VECT*VECT 4934 VECTORSIZE=SUM(SVECT) 4935 VECTORSIZE=VECTORSIZE**(0.5) 4936 end function VECTORSIZE 4937 4938 4939!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4940 4941 function determinant(aa) result(det) 4942 real(wp), intent(in) :: aa(:,:) 4943 real(wp) :: det 4944 4945 det = aa(1,1) * (aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3))& 4946 & - aa(1,2) * (aa(2,1) * aa(3,3) - aa(3,1) * aa(2,3))& 4947 & + aa(1,3) * (aa(2,1) * aa(3,2) - aa(3,1) * aa(2,2)) 4948 4949 end function determinant 4950 4951 4952 subroutine inv_cell(x,a) 4953 real(wp), intent(in) :: x(3,3) 4954 real(wp), intent(out) :: a(3,3) 4955 real(wp) det 4956 4957 a = 0.0 4958 det = determinant(x) 4959 ! write(*,*)'Det:',det 4960 a(1,1)=x(2,2)*x(3,3)-x(2,3)*x(3,2) 4961 a(2,1)=x(2,3)*x(3,1)-x(2,1)*x(3,3) 4962 a(3,1)=x(2,1)*x(3,2)-x(2,2)*x(3,1) 4963 a(1,2)=x(1,3)*x(3,2)-x(1,2)*x(3,3) 4964 a(2,2)=x(1,1)*x(3,3)-x(1,3)*x(3,1) 4965 a(3,2)=x(1,2)*x(3,1)-x(1,1)*x(3,2) 4966 a(1,3)=x(1,2)*x(2,3)-x(1,3)*x(2,2) 4967 a(2,3)=x(1,3)*x(2,1)-x(1,1)*x(2,3) 4968 a(3,3)=x(1,1)*x(2,2)-x(1,2)*x(2,1) 4969 a=a/det 4970 end subroutine inv_cell 4971 4972!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4973 4974 subroutine xyz_to_abc(xyz,abc,lat,n) 4975 integer,intent(in) :: n 4976 real(wp), INTENT(in) :: xyz(3,n) 4977 real(wp), intent(in) :: lat(3,3) 4978 real(wp), intent(out) :: abc(3,n) 4979 4980 real(wp) lat_1(3,3) 4981 integer i,j,k 4982 4983 call inv_cell(lat,lat_1) 4984 4985 abc(:,:n)=0.0d0 4986 do i=1,n 4987 do j=1,3 4988 do k=1,3 4989 abc(j,i)=abc(j,i)+lat_1(j,k)*xyz(k,i) 4990 end do 4991 abc(j,i)=dmod(abc(j,i),1.0d0) 4992 end do 4993 end do 4994 4995 end subroutine xyz_to_abc 4996 4997!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4998 4999 subroutine abc_to_xyz(abc,xyz,lat,n) 5000 real(wp), INTENT(in) :: abc(3,*) 5001 real(wp), intent(in) :: lat(3,3) 5002 real(wp), intent(out) :: xyz(3,*) 5003 integer,intent(in) :: n 5004 5005 integer i,j,k 5006 5007 xyz(:,:n)=0.0d0 5008 do i=1,n 5009 do j=1,3 5010 do k=1,3 5011 xyz(j,i)=xyz(j,i)+lat(j,k)*abc(k,i) 5012 end do 5013 end do 5014 end do 5015 5016 end subroutine abc_to_xyz 5017 5018!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 5019 5020end module dftd3_core 5021