1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Debug the derivatives of the the rotational matrices 8!> 9!> \param sepi ... 10!> \param sepj ... 11!> \param rjiv ... 12!> \param ij_matrix ... 13!> \param do_invert ... 14!> \date 04.2008 [tlaino] 15!> \author Teodoro Laino [tlaino] - University of Zurich 16! ************************************************************************************************** 17SUBROUTINE check_rotmat_der(sepi, sepj, rjiv, ij_matrix, do_invert) 18 19 USE kinds, ONLY: dp 20 USE semi_empirical_int_utils, ONLY: rotmat 21 USE semi_empirical_types, ONLY: rotmat_create, & 22 rotmat_release, & 23 rotmat_type, & 24 semi_empirical_type, & 25 se_int_control_type 26#include "./base/base_uses.f90" 27 IMPLICIT NONE 28 TYPE(semi_empirical_type), POINTER :: sepi, sepj 29 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rjiv 30 TYPE(rotmat_type), POINTER :: ij_matrix 31 LOGICAL, INTENT(IN) :: do_invert 32 33 CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', & 34 routineN = 'check_rotmat_der', routineP = moduleN//':'//routineN 35 36 LOGICAL :: check_value 37 REAL(KIND=dp) :: dx, r, r0(3), x(3) 38 TYPE(rotmat_type), POINTER :: matrix, matrix_m, matrix_n, & 39 matrix_p 40 INTEGER :: imap(3), i, j, k, l 41 42 NULLIFY (matrix_m, matrix_p, matrix_n, matrix) 43 CALL rotmat_create(matrix_p) 44 CALL rotmat_create(matrix_m) 45 CALL rotmat_create(matrix_n) 46 dx = 1.0E-6_dp 47 imap(1) = 1 48 imap(2) = 2 49 imap(3) = 3 50 IF (do_invert) THEN 51 imap(1) = 3 52 imap(2) = 2 53 imap(3) = 1 54 END IF 55 ! Check derivatives: analytical VS numerical 56 WRITE (*, *) "DEBUG::"//routineP 57 DO j = 1, 3 58 x = 0.0_dp 59 x(imap(j)) = dx 60 DO i = 1, 2 61 IF (i == 1) matrix => matrix_p 62 IF (i == 2) matrix => matrix_m 63 r0 = rjiv + (-1.0_dp)**(i - 1)*x 64 r = SQRT(DOT_PRODUCT(r0, r0)) 65 CALL rotmat(sepi, sepj, r0, r, matrix, do_derivatives=.FALSE., debug_invert=do_invert) 66 END DO 67 ! SP 68 matrix_n%sp_d(j, :, :) = (matrix_p%sp - matrix_m%sp)/(2.0_dp*dx) 69 DO i = 1, 3 70 DO k = 1, 3 71 IF (.NOT. check_value(matrix_n%sp_d(j, k, i), ij_matrix%sp_d(j, k, i), dx, 0.1_dp)) THEN 72 WRITE (*, *) "ERROR for SP rotation matrix derivative SP(j,k,i), j,k,i::", j, k, i 73 CPABORT("") 74 END IF 75 END DO 76 END DO 77 ! PP 78 matrix_n%pp_d(j, :, :, :) = (matrix_p%pp - matrix_m%pp)/(2.0_dp*dx) 79 DO i = 1, 3 80 DO k = 1, 3 81 DO l = 1, 6 82 IF (.NOT. check_value(matrix_n%pp_d(j, l, k, i), ij_matrix%pp_d(j, l, k, i), dx, 0.1_dp)) THEN 83 WRITE (*, *) "ERROR for PP rotation matrix derivative PP(j,l,k,i), j,l,k,i::", j, l, k, i 84 CPABORT("") 85 END IF 86 END DO 87 END DO 88 END DO 89 ! d-orbitals debug 90 IF (sepi%dorb .OR. sepj%dorb) THEN 91 ! SD 92 matrix_n%sd_d(j, :, :) = (matrix_p%sd - matrix_m%sd)/(2.0_dp*dx) 93 DO i = 1, 5 94 DO k = 1, 5 95 IF (.NOT. check_value(matrix_n%sd_d(j, k, i), ij_matrix%sd_d(j, k, i), dx, 0.1_dp)) THEN 96 WRITE (*, *) "ERROR for SD rotation matrix derivative SD(j,k,i), j,k,i::", j, k, i 97 CPABORT("") 98 END IF 99 END DO 100 END DO 101 ! DP 102 matrix_n%pd_d(j, :, :, :) = (matrix_p%pd - matrix_m%pd)/(2.0_dp*dx) 103 DO i = 1, 3 104 DO k = 1, 5 105 DO l = 1, 15 106 IF (.NOT. check_value(matrix_n%pd_d(j, l, k, i), ij_matrix%pd_d(j, l, k, i), dx, 0.1_dp)) THEN 107 WRITE (*, *) "ERROR for DP rotation matrix derivative DP(j,l,k,i), j,l,k,i::", j, l, k, i 108 CPABORT("") 109 END IF 110 END DO 111 END DO 112 END DO 113 ! DD 114 matrix_n%dd_d(j, :, :, :) = (matrix_p%dd - matrix_m%dd)/(2.0_dp*dx) 115 DO i = 1, 5 116 DO k = 1, 5 117 DO l = 1, 15 118 IF (.NOT. check_value(matrix_n%dd_d(j, l, k, i), ij_matrix%dd_d(j, l, k, i), dx, 0.1_dp)) THEN 119 WRITE (*, *) "ERROR for DD rotation matrix derivative DD(j,l,k,i), j,l,k,i::", j, l, k, i 120 CPABORT("") 121 END IF 122 END DO 123 END DO 124 END DO 125 END IF 126 END DO 127 CALL rotmat_release(matrix_p) 128 CALL rotmat_release(matrix_m) 129 CALL rotmat_release(matrix_n) 130END SUBROUTINE check_rotmat_der 131 132! ************************************************************************************************** 133!> \brief Check Numerical Vs Analytical 134!> \param sepi ... 135!> \param sepj ... 136!> \param rijv ... 137!> \param se_int_control ... 138!> \param se_taper ... 139!> \param invert ... 140!> \param ii ... 141!> \param kk ... 142!> \param v_d ... 143!> \par History 144!> 04.2008 created [tlaino] 145!> \author Teodoro Laino - Zurich University 146!> \note 147!> Debug routine 148! ************************************************************************************************** 149SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, v_d) 150 151 USE kinds, ONLY: dp 152 USE semi_empirical_int_utils, ONLY: rotmat 153 USE semi_empirical_int_arrays, ONLY: indexb 154 USE semi_empirical_int_num, ONLY: terep_num 155 USE semi_empirical_types, ONLY: semi_empirical_type, & 156 rotmat_type, & 157 rotmat_create, & 158 rotmat_release, & 159 se_int_control_type, & 160 se_taper_type 161 USE semi_empirical_int_utils, ONLY: rot_2el_2c_first 162#include "./base/base_uses.f90" 163 IMPLICIT NONE 164 TYPE(semi_empirical_type), POINTER :: sepi, sepj 165 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv 166 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 167 LOGICAL, INTENT(IN) :: invert 168 TYPE(se_taper_type), POINTER :: se_taper 169 INTEGER, INTENT(IN) :: ii, kk 170 REAL(KIND=dp), DIMENSION(3, 45, 45), & 171 INTENT(IN) :: v_d 172 173 CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', & 174 routineN = 'rot_2el_2c_first', routineP = moduleN//':'//routineN 175 176 REAL(KIND=dp), DIMENSION(491) :: rep 177 LOGICAL, DIMENSION(45, 45) :: logv 178 REAL(KIND=dp), DIMENSION(45, 45) :: v_n, v_p, v_m 179 180 LOGICAL :: check_value 181 REAL(KIND=dp) :: dx, r, r0(3), x(3) 182 TYPE(rotmat_type), POINTER :: matrix 183 INTEGER :: imap(3), i, j, k, limkl 184 185 NULLIFY (matrix) 186 dx = 1.0E-6_dp 187 imap(1) = 1 188 imap(2) = 2 189 imap(3) = 3 190 IF (invert) THEN 191 imap(1) = 3 192 imap(2) = 2 193 imap(3) = 1 194 END IF 195 limkl = indexb(kk, kk) 196 ! Check derivatives: analytical VS numerical 197 WRITE (*, *) "DEBUG::"//routineP 198 DO j = 1, 3 199 x = 0.0_dp 200 x(imap(j)) = dx 201 DO i = 1, 2 202 r0 = rijv + (-1.0_dp)**(i - 1)*x 203 r = SQRT(DOT_PRODUCT(r0, r0)) 204 205 CALL rotmat_create(matrix) 206 CALL rotmat(sepi, sepj, r0, r, matrix, do_derivatives=.FALSE., debug_invert=invert) 207 208 ! Compute integrals in diatomic frame 209 CALL terep_num(sepi, sepj, r, rep, se_taper=se_taper, se_int_control=se_int_control) 210 211 IF (i == 1) THEN 212 CALL rot_2el_2c_first(sepi, sepj, r0, se_int_control, se_taper, invert, ii, kk, rep, logv, matrix, & 213 v_p, lgrad=.FALSE.) 214 END IF 215 IF (i == 2) THEN 216 CALL rot_2el_2c_first(sepi, sepj, r0, se_int_control, se_taper, invert, ii, kk, rep, logv, matrix, & 217 v_m, lgrad=.FALSE.) 218 END IF 219 CALL rotmat_release(matrix) 220 END DO 221 ! Check numerical VS analytical 222 DO i = 1, 45 223 DO k = 1, limkl 224 ! Compute the numerical derivative 225 v_n(i, k) = (v_p(i, k) - v_m(i, k))/(2.0_dp*dx) 226 IF (.NOT. check_value(v_d(j, i, k), v_n(i, k), dx, 0.1_dp)) THEN 227 WRITE (*, *) "ERROR for rot_2el_2c_first derivative V_D(j,i,k), j,i,k::", j, i, k 228 CPABORT("") 229 END IF 230 END DO 231 END DO 232 END DO 233END SUBROUTINE rot_2el_2c_first_debug 234 235! ************************************************************************************************** 236!> \brief Check Numerical Vs Analytical ssss 237!> \param sepi ... 238!> \param sepj ... 239!> \param r ... 240!> \param dssss ... 241!> \param itype ... 242!> \param se_int_control ... 243!> \param se_taper ... 244!> \par History 245!> 04.2008 created [tlaino] 246!> \author Teodoro Laino - Zurich University 247!> \note 248!> Debug routine 249! ************************************************************************************************** 250SUBROUTINE check_dssss_nucint_ana(sepi, sepj, r, dssss, itype, se_int_control, se_taper) 251 252 USE kinds, ONLY: dp 253 USE semi_empirical_int_num, ONLY: ssss_nucint_num 254 USE semi_empirical_types, ONLY: semi_empirical_type, & 255 se_int_control_type, & 256 se_taper_type 257#include "./base/base_uses.f90" 258 IMPLICIT NONE 259 TYPE(semi_empirical_type), POINTER :: sepi, sepj 260 REAL(dp), INTENT(IN) :: r 261 REAL(dp), INTENT(IN) :: dssss 262 INTEGER, INTENT(IN) :: itype 263 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 264 TYPE(se_taper_type), POINTER :: se_taper 265 266 CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', & 267 routineN = 'check_dssss_nucint_ana', routineP = moduleN//':'//routineN 268 269 REAL(dp) :: delta, nssss, od, rn, ssssm, ssssp 270 LOGICAL :: check_value 271 272 delta = 1.0E-8_dp 273 od = 0.5_dp/delta 274 rn = r + delta 275 CALL ssss_nucint_num(sepi, sepj, rn, ssssp, itype, se_taper, se_int_control) 276 rn = r - delta 277 CALL ssss_nucint_num(sepi, sepj, rn, ssssm, itype, se_taper, se_int_control) 278 nssss = od*(ssssp - ssssm) 279 ! check 280 WRITE (*, *) "DEBUG::"//routineP 281 IF (.NOT. check_value(nssss, dssss, delta, 0.1_dp)) THEN 282 WRITE (*, *) "ERROR for SSSS derivative SSSS" 283 CPABORT("") 284 END IF 285END SUBROUTINE check_dssss_nucint_ana 286 287! ************************************************************************************************** 288!> \brief Check Numerical Vs Analytical core 289!> \param sepi ... 290!> \param sepj ... 291!> \param r ... 292!> \param dcore ... 293!> \param itype ... 294!> \param se_int_control ... 295!> \param se_taper ... 296!> \par History 297!> 04.2008 created [tlaino] 298!> \author Teodoro Laino - Zurich University 299!> \note 300!> Debug routine 301! ************************************************************************************************** 302SUBROUTINE check_dcore_nucint_ana(sepi, sepj, r, dcore, itype, se_int_control, se_taper) 303 304 USE kinds, ONLY: dp 305 USE semi_empirical_int_num, ONLY: core_nucint_num 306 USE semi_empirical_types, ONLY: semi_empirical_type, & 307 se_int_control_type, & 308 se_taper_type 309#include "./base/base_uses.f90" 310 IMPLICIT NONE 311 TYPE(semi_empirical_type), POINTER :: sepi, sepj 312 REAL(dp), INTENT(IN) :: r 313 REAL(dp), DIMENSION(10, 2), INTENT(IN) :: dcore 314 INTEGER, INTENT(IN) :: itype 315 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 316 TYPE(se_taper_type), POINTER :: se_taper 317 318 CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', & 319 routineN = 'check_dcore_nucint_ana', routineP = moduleN//':'//routineN 320 321 INTEGER :: i, j 322 REAL(dp) :: delta, od, rn 323 REAL(dp), DIMENSION(10, 2) :: corem, corep, ncore 324 LOGICAL :: check_value 325 326 delta = 1.0E-8_dp 327 od = 0.5_dp/delta 328 rn = r + delta 329 CALL core_nucint_num(sepi, sepj, rn, corep, itype, se_taper, se_int_control) 330 rn = r - delta 331 CALL core_nucint_num(sepi, sepj, rn, corem, itype, se_taper, se_int_control) 332 ncore = od*(corep - corem) 333 ! check 334 WRITE (*, *) "DEBUG::"//routineP 335 DO i = 1, 2 336 DO j = 1, 10 337 IF (.NOT. check_value(ncore(j, i), dcore(j, i), delta, 0.1_dp)) THEN 338 WRITE (*, *) "ERROR for CORE derivative CORE(j,i), j,i::", j, i 339 CPABORT("") 340 END IF 341 END DO 342 END DO 343END SUBROUTINE check_dcore_nucint_ana 344 345! ************************************************************************************************** 346!> \brief Low level comparison function between numberical and analytical values 347!> \param num ... 348!> \param ana ... 349!> \param minval ... 350!> \param thrs ... 351!> \return ... 352!> \par History 353!> 04.2008 created [tlaino] 354!> \author Teodoro Laino - Zurich University 355!> \note 356!> Debug routine 357! ************************************************************************************************** 358FUNCTION check_value(num, ana, minval, thrs) RESULT(passed) 359 USE kinds, ONLY: dp 360 IMPLICIT NONE 361 REAL(KIND=dp) :: num, ana, thrs, minval 362 LOGICAL :: passed 363 364 passed = .TRUE. 365 366 IF ((ABS(num) < minval) .AND. (ABS(ana) > minval)) THEN 367 WRITE (*, *) "WARNING ---> ", num, ana, thrs 368 ELSE IF ((ABS(num) > minval) .AND. (ABS(ana) < minval)) THEN 369 WRITE (*, *) "WARNING ---> ", num, ana, thrs 370 ELSE IF ((ABS(num) < minval) .AND. (ABS(ana) < minval)) THEN 371 ! skip.. 372 RETURN 373 END IF 374 IF (ABS((num - ana)/num*100._dp) > thrs) THEN 375 WRITE (*, *) ABS(num - ana)/num*100._dp, thrs 376 passed = .FALSE. 377 END IF 378 IF (.NOT. passed) THEN 379 WRITE (*, *) "ANALYTICAL ::", ana 380 WRITE (*, *) "NUMERICAL ::", num 381 END IF 382END FUNCTION check_value 383 384! ************************************************************************************************** 385!> \brief Check Numerical Vs Analytical 386!> \param sepi ... 387!> \param sepj ... 388!> \param rijv ... 389!> \param itype ... 390!> \param se_int_control ... 391!> \param se_taper ... 392!> \param e1b ... 393!> \param e2a ... 394!> \param de1b ... 395!> \param de2a ... 396!> \par History 397!> 04.2008 created [tlaino] 398!> \author Teodoro Laino - Zurich University 399!> \note 400!> Debug routine 401! ************************************************************************************************** 402SUBROUTINE check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a) 403 404 USE kinds, ONLY: dp 405 USE semi_empirical_int_num, ONLY: rotnuc_num, & 406 drotnuc_num 407 USE semi_empirical_types, ONLY: semi_empirical_type, & 408 se_int_control_type, & 409 se_taper_type 410#include "./base/base_uses.f90" 411 IMPLICIT NONE 412 TYPE(semi_empirical_type), POINTER :: sepi, sepj 413 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 414 INTEGER, INTENT(IN) :: itype 415 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 416 TYPE(se_taper_type), POINTER :: se_taper 417 REAL(dp), DIMENSION(45), INTENT(IN), & 418 OPTIONAL :: e1b, e2a 419 REAL(dp), DIMENSION(3, 45), & 420 INTENT(IN), OPTIONAL :: de1b, de2a 421 422 CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', & 423 routineN = 'check_drotnuc_ana', routineP = moduleN//':'//routineN 424 425 INTEGER :: i, j 426 LOGICAL :: l_de1b, l_de2a, l_e1b, l_e2a, & 427 lgrad, check_value 428 REAL(dp) :: delta 429 REAL(KIND=dp), DIMENSION(45) :: e1b2, e2a2 430 REAL(KIND=dp), DIMENSION(3, 45) :: de1b2, de2a2 431 432 l_e1b = PRESENT(e1b) 433 l_e2a = PRESENT(e2a) 434 l_de1b = PRESENT(de1b) 435 l_de2a = PRESENT(de2a) 436 lgrad = l_de1b .OR. l_de2a 437 delta = 1.0E-5_dp 438 ! Check value of integrals 439 WRITE (*, *) "DEBUG::"//routineP 440 CALL rotnuc_num(sepi, sepj, rijv, e1b2, e2a2, itype, se_int_control, se_taper=se_taper) 441 IF (l_e1b) THEN 442 DO j = 1, 45 443 IF (.NOT. check_value(e1b2(j), e1b(j), delta, 0.1_dp)) THEN 444 WRITE (*, *) "ERROR for E1B value E1B(j), j::", j 445 CPABORT("") 446 END IF 447 END DO 448 END IF 449 IF (l_e2a) THEN 450 DO J = 1, 45 451 IF (.NOT. check_value(e2a2(j), e2a(j), delta, 0.1_dp)) THEN 452 WRITE (*, *) "ERROR for E2A value E2A(j), j::", j 453 CPABORT("") 454 END IF 455 END DO 456 END IF 457 458 ! Check derivatives 459 IF (lgrad) THEN 460 CALL drotnuc_num(sepi, sepj, rijv, de1b2, de2a2, itype, delta=delta, & 461 se_int_control=se_int_control, se_taper=se_taper) 462 IF (l_de1b) THEN 463 DO i = 1, 3 464 DO j = 1, 45 465 ! Additional check on the value of the integral before checking derivatives 466 IF (ABS(e1b2(j)) > delta) THEN 467 IF (.NOT. check_value(de1b2(i, j), de1b(i, j), delta, 0.1_dp)) THEN 468 WRITE (*, *) "ERROR for derivative of E1B. DE1B(i,j), i,j::", i, j 469 CPABORT("") 470 END IF 471 END IF 472 END DO 473 END DO 474 END IF 475 IF (l_de2a) THEN 476 DO i = 1, 3 477 DO j = 1, 45 478 ! Additional check on the value of the integral before checking derivatives 479 IF (ABS(e2a2(j)) > delta) THEN 480 IF (.NOT. check_value(de2a2(i, j), de2a(i, j), delta, 0.1_dp)) THEN 481 WRITE (*, *) "ERROR for derivative of E2A. DE2A(i,j), i,j::", i, j 482 CPABORT("") 483 END IF 484 END IF 485 END DO 486 END DO 487 END IF 488 END IF 489END SUBROUTINE check_drotnuc_ana 490 491! ************************************************************************************************** 492!> \brief Check Numerical Vs Analytical CORE-CORE 493!> \param sepi ... 494!> \param sepj ... 495!> \param rijv ... 496!> \param itype ... 497!> \param se_int_control ... 498!> \param se_taper ... 499!> \param enuc ... 500!> \param denuc ... 501!> \par History 502!> 04.2007 created [tlaino] 503!> \author Teodoro Laino - Zurich University 504!> \note 505!> Debug routine 506! ************************************************************************************************** 507SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, enuc, denuc) 508 509 USE kinds, ONLY: dp 510 USE semi_empirical_int_num, ONLY: corecore_num, & 511 dcorecore_num 512 USE semi_empirical_types, ONLY: semi_empirical_type, & 513 se_int_control_type, & 514 se_taper_type 515#include "./base/base_uses.f90" 516 IMPLICIT NONE 517 TYPE(semi_empirical_type), POINTER :: sepi, sepj 518 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 519 INTEGER, INTENT(IN) :: itype 520 REAL(dp), INTENT(IN), OPTIONAL :: enuc 521 REAL(dp), DIMENSION(3), INTENT(IN), & 522 OPTIONAL :: denuc 523 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 524 TYPE(se_taper_type), POINTER :: se_taper 525 526 CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', & 527 routineN = 'check_dcorecore_ana', routineP = moduleN//':'//routineN 528 529 LOGICAL :: check_value 530 INTEGER :: j 531 REAL(dp) :: enuc_num, delta 532 REAL(dp), DIMENSION(3) :: denuc_num 533 534 WRITE (*, *) "DEBUG::"//routineP 535 delta = 1.0E-7_dp 536 ! check 537 IF (PRESENT(enuc)) THEN 538 CALL corecore_num(sepi, sepj, rijv, enuc_num, itype, se_int_control, se_taper) 539 IF (.NOT. check_value(enuc, enuc_num, delta, 0.001_dp)) THEN 540 WRITE (*, *) "ERROR for CORE-CORE energy value (numerical different from analytical)!!" 541 CPABORT("") 542 END IF 543 END IF 544 IF (PRESENT(denuc)) THEN 545 CALL dcorecore_num(sepi, sepj, rijv, denuc_num, itype, delta, se_int_control, se_taper) 546 DO j = 1, 3 547 IF (.NOT. check_value(denuc(j), denuc_num(j), delta, 0.001_dp)) THEN 548 WRITE (*, *) "ERROR for CORE-CORE energy derivative value (numerical different from analytical). DENUC(j), j::", j 549 CPABORT("") 550 END IF 551 END DO 552 END IF 553END SUBROUTINE check_dcorecore_ana 554 555! ************************************************************************************************** 556!> \brief Check Numerical Vs Analytical DTEREP_ANA 557!> \param sepi ... 558!> \param sepj ... 559!> \param r ... 560!> \param ri ... 561!> \param dri ... 562!> \param se_int_control ... 563!> \param se_taper ... 564!> \param lgrad ... 565!> \par History 566!> 04.2007 created [tlaino] 567!> \author Teodoro Laino - Zurich University 568!> \note 569!> Debug routine 570! ************************************************************************************************** 571SUBROUTINE check_dterep_ana(sepi, sepj, r, ri, dri, se_int_control, se_taper, lgrad) 572 573 USE kinds, ONLY: dp 574 USE semi_empirical_int_num, ONLY: terep_num 575 USE semi_empirical_types, ONLY: semi_empirical_type, & 576 se_int_control_type, & 577 se_taper_type 578#include "./base/base_uses.f90" 579 IMPLICIT NONE 580 TYPE(semi_empirical_type), POINTER :: sepi, sepj 581 REAL(dp), INTENT(IN) :: r 582 REAL(dp), DIMENSION(491), INTENT(IN) :: ri, dri 583 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 584 LOGICAL, INTENT(IN) :: lgrad 585 TYPE(se_taper_type), POINTER :: se_taper 586 587 CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', & 588 routineN = 'check_dterep_ana', routineP = moduleN//':'//routineN 589 590 LOGICAL :: check_value 591 INTEGER :: j, i 592 REAL(dp) :: delta, od, rn 593 REAL(dp), DIMENSION(491) :: nri, ri0, rim, rip 594 595 delta = 1.0E-8_dp 596 od = 0.5_dp/delta 597 rn = r 598 CALL terep_num(sepi, sepj, rn, ri0, se_taper, se_int_control) 599 IF (lgrad) THEN 600 rn = r + delta 601 CALL terep_num(sepi, sepj, rn, rip, se_taper, se_int_control) 602 rn = r - delta 603 CALL terep_num(sepi, sepj, rn, rim, se_taper, se_int_control) 604 nri = od*(rip - rim) 605 END IF 606 ! check 607 WRITE (*, *) "DEBUG::"//routineP 608 DO j = 1, 491 609 IF (ABS(ri(j) - ri0(j)) > EPSILON(0.0_dp)) THEN 610 WRITE (*, *) "Error in value of the integral RI", j, ri(j), ri0(j) 611 CPABORT("") 612 END IF 613 IF (lgrad) THEN 614 IF (.NOT. check_value(dri(j), nri(j), delta*100.0_dp, 0.1_dp)) THEN 615 WRITE (*, *) "ERROR for derivative of RI integral, RI(j), j::", j 616 WRITE (*, *) "FULL SET OF INTEGRALS: INDX ANAL NUM DIFF" 617 DO i = 1, 491 618 WRITE (*, '(I5,3F15.9)') i, dri(i), nri(i), nri(i) - dri(i) 619 END DO 620 CPABORT("") 621 END IF 622 END IF 623 END DO 624END SUBROUTINE check_dterep_ana 625 626! ************************************************************************************************** 627!> \brief Check Numerical Vs Analytical ROTINT_ANA 628!> \param sepi ... 629!> \param sepj ... 630!> \param rijv ... 631!> \param w ... 632!> \param dw ... 633!> \param se_int_control ... 634!> \param se_taper ... 635!> \par History 636!> 04.2008 created [tlaino] 637!> \author Teodoro Laino - Zurich University 638!> \note 639!> Debug routine 640! ************************************************************************************************** 641SUBROUTINE check_rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper) 642 643 USE kinds, ONLY: dp 644 USE semi_empirical_int_num, ONLY: rotint_num, & 645 drotint_num 646 USE semi_empirical_types, ONLY: semi_empirical_type, & 647 se_int_control_type, & 648 se_taper_type 649#include "./base/base_uses.f90" 650 IMPLICIT NONE 651 TYPE(semi_empirical_type), POINTER :: sepi, sepj 652 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 653 REAL(dp), DIMENSION(2025), INTENT(IN), & 654 OPTIONAL :: w 655 REAL(dp), DIMENSION(3, 2025), & 656 INTENT(IN), OPTIONAL :: dw 657 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 658 TYPE(se_taper_type), POINTER :: se_taper 659 660 CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', & 661 routineN = 'rotint_ana', routineP = moduleN//':'//routineN 662 663 REAL(dp), DIMENSION(2025) :: w2 664 REAL(dp), DIMENSION(3, 2025) :: dw2 665 LOGICAL :: check_value 666 INTEGER :: j, i 667 REAL(KIND=dp) :: delta 668 669 delta = 1.0E-6_dp 670 WRITE (*, *) "DEBUG::"//routineP 671 IF (PRESENT(w)) THEN 672 w2 = 0.0_dp 673 CALL rotint_num(sepi, sepj, rijv, w2, se_int_control, se_taper=se_taper) 674 DO j = 1, 2025 675 IF (.NOT. check_value(w(j), w2(j), delta, 0.1_dp)) THEN 676 WRITE (*, *) "ERROR for integral value W(j), j::", j 677 CPABORT("") 678 END IF 679 END DO 680 END IF 681 IF (PRESENT(dw)) THEN 682 ! Numerical derivatives are obviosly a big problem.. 683 ! First of all let's decide if the value we get for delta is compatible 684 ! with a reasonable value of the integral.. (compatible if the value of the 685 ! integral is greater than 1.0E-6) 686 CALL drotint_num(sepi, sepj, rijv, dw2, delta=delta, se_int_control=se_int_control, se_taper=se_taper) 687 CALL rotint_num(sepi, sepj, rijv, w2, se_int_control=se_int_control, se_taper=se_taper) 688 DO i = 1, 3 689 DO j = 1, 2025 690 IF ((ABS(w2(j)) > delta) .AND. (ABS(dw2(i, j)) > delta*10)) THEN 691 IF (.NOT. check_value(dw(i, j), dw2(i, j), delta, 0.1_dp)) THEN 692 WRITE (*, *) "ERROR for derivative of the integral value W(j). DW(i,j) i,j::", i, j 693 CPABORT("") 694 END IF 695 END IF 696 END DO 697 END DO 698 END IF 699END SUBROUTINE check_rotint_ana 700