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