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