1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Module to perform a counterpoise correction (BSSE) 8!> \par History 9!> 6.2005 created [tlaino] 10!> \author Teodoro Laino 11! ************************************************************************************************** 12MODULE bsse 13 USE atomic_kind_types, ONLY: get_atomic_kind 14 USE cell_types, ONLY: cell_type 15 USE cp2k_info, ONLY: write_restart_header 16 USE cp_external_control, ONLY: external_control 17 USE cp_log_handling, ONLY: cp_get_default_logger,& 18 cp_logger_type 19 USE cp_output_handling, ONLY: cp_add_iter_level,& 20 cp_iterate,& 21 cp_print_key_finished_output,& 22 cp_print_key_unit_nr,& 23 cp_rm_iter_level 24 USE cp_para_types, ONLY: cp_para_env_type 25 USE cp_subsys_methods, ONLY: create_small_subsys 26 USE cp_subsys_types, ONLY: cp_subsys_get,& 27 cp_subsys_release,& 28 cp_subsys_type 29 USE force_env_types, ONLY: force_env_get,& 30 force_env_type 31 USE global_types, ONLY: global_environment_type 32 USE input_constants, ONLY: do_qs 33 USE input_section_types, ONLY: section_vals_get,& 34 section_vals_get_subs_vals,& 35 section_vals_type,& 36 section_vals_val_get,& 37 section_vals_val_set,& 38 section_vals_write 39 USE kinds, ONLY: default_string_length,& 40 dp 41 USE memory_utilities, ONLY: reallocate 42 USE particle_list_types, ONLY: particle_list_type 43 USE qs_energy, ONLY: qs_energies 44 USE qs_energy_types, ONLY: qs_energy_type 45 USE qs_environment, ONLY: qs_init 46 USE qs_environment_types, ONLY: get_qs_env,& 47 qs_env_create,& 48 qs_env_release,& 49 qs_environment_type 50 USE string_utilities, ONLY: compress 51#include "./base/base_uses.f90" 52 53 IMPLICIT NONE 54 PRIVATE 55 56 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 57 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bsse' 58 59 PUBLIC :: do_bsse_calculation 60 61CONTAINS 62 63! ************************************************************************************************** 64!> \brief Perform an COUNTERPOISE CORRECTION (BSSE) 65!> For a 2-body system the correction scheme can be represented as: 66!> 67!> E_{AB}^{2} = E_{AB}(AB) - E_A(AB) - E_B(AB) [BSSE-corrected interaction energy] 68!> E_{AB}^{2,uncorr} = E_{AB}(AB) - E_A(A) - E_B(B) 69!> E_{AB}^{CP} = E_{AB}(AB) + [ E_A(A) - E_A(AB) ] + [ E_B(B) - E_B(AB) ] 70!> [CP-corrected total energy of AB] 71!> \param force_env ... 72!> \param globenv ... 73!> \par History 74!> 06.2005 created [tlaino] 75!> \author Teodoro Laino 76! ************************************************************************************************** 77 SUBROUTINE do_bsse_calculation(force_env, globenv) 78 TYPE(force_env_type), POINTER :: force_env 79 TYPE(global_environment_type), POINTER :: globenv 80 81 CHARACTER(LEN=*), PARAMETER :: routineN = 'do_bsse_calculation', & 82 routineP = moduleN//':'//routineN 83 84 INTEGER :: i, istart, k, num_of_conf, Num_of_Frag 85 INTEGER, DIMENSION(:, :), POINTER :: conf 86 LOGICAL :: explicit, should_stop 87 REAL(KIND=dp), DIMENSION(:), POINTER :: Em 88 TYPE(cp_logger_type), POINTER :: logger 89 TYPE(section_vals_type), POINTER :: bsse_section, fragment_energies_section, & 90 n_frags, root_section 91 92 NULLIFY (bsse_section, n_frags, Em, conf) 93 logger => cp_get_default_logger() 94 root_section => force_env%root_section 95 bsse_section => section_vals_get_subs_vals(force_env%force_env_section, "BSSE") 96 n_frags => section_vals_get_subs_vals(bsse_section, "FRAGMENT") 97 CALL section_vals_get(n_frags, n_repetition=Num_of_Frag) 98 99 ! Number of configurations 100 num_of_conf = 0 101 DO k = 1, Num_of_frag 102 num_of_conf = num_of_conf + FACT(Num_of_frag)/(FACT(k)*FACT(Num_of_frag - k)) 103 END DO 104 ALLOCATE (conf(num_of_conf, Num_of_frag)) 105 ALLOCATE (Em(num_of_conf)) 106 CALL gen_Nbody_conf(Num_of_frag, conf) 107 108 should_stop = .FALSE. 109 istart = 0 110 fragment_energies_section => section_vals_get_subs_vals(bsse_section, "FRAGMENT_ENERGIES") 111 CALL section_vals_get(fragment_energies_section, explicit=explicit) 112 IF (explicit) THEN 113 CALL section_vals_val_get(fragment_energies_section, "_DEFAULT_KEYWORD_", n_rep_val=istart) 114 DO i = 1, istart 115 CALL section_vals_val_get(fragment_energies_section, "_DEFAULT_KEYWORD_", r_val=Em(i), & 116 i_rep_val=i) 117 END DO 118 END IF 119 ! Setup the iteration level for BSSE 120 CALL cp_add_iter_level(logger%iter_info, "BSSE") 121 CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=istart) 122 123 ! Evaluating the energy of the N-body cluster terms 124 DO i = istart + 1, num_of_conf 125 CALL cp_iterate(logger%iter_info, last=(i == num_of_conf), iter_nr=i) 126 CALL eval_bsse_energy(conf(i, :), Em(i), force_env, n_frags, & 127 root_section, globenv, should_stop) 128 IF (should_stop) EXIT 129 130 ! If no signal was received in the inner loop let's check also at this stage 131 CALL external_control(should_stop, "BSSE", globenv=globenv) 132 IF (should_stop) EXIT 133 134 ! Dump Restart info only if the calculation of the energy of a configuration 135 ! ended nicely.. 136 CALL section_vals_val_set(fragment_energies_section, "_DEFAULT_KEYWORD_", r_val=Em(i), & 137 i_rep_val=i) 138 CALL write_bsse_restart(bsse_section, root_section) 139 END DO 140 IF (.NOT. should_stop) CALL dump_bsse_results(conf, Em, num_of_frag, bsse_section) 141 CALL cp_rm_iter_level(logger%iter_info, "BSSE") 142 DEALLOCATE (Em) 143 DEALLOCATE (conf) 144 145 END SUBROUTINE do_bsse_calculation 146 147! ************************************************************************************************** 148!> \brief Evaluate the N-body energy contribution to the BSSE evaluation 149!> \param conf ... 150!> \param Em ... 151!> \param force_env ... 152!> \param n_frags ... 153!> \param root_section ... 154!> \param globenv ... 155!> \param should_stop ... 156!> \par History 157!> 07.2005 created [tlaino] 158!> \author Teodoro Laino 159! ************************************************************************************************** 160 SUBROUTINE eval_bsse_energy(conf, Em, force_env, n_frags, root_section, & 161 globenv, should_stop) 162 INTEGER, DIMENSION(:), INTENT(IN) :: conf 163 REAL(KIND=dp), INTENT(OUT) :: Em 164 TYPE(force_env_type), POINTER :: force_env 165 TYPE(section_vals_type), POINTER :: n_frags, root_section 166 TYPE(global_environment_type), POINTER :: globenv 167 LOGICAL, INTENT(OUT) :: should_stop 168 169 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_bsse_energy', & 170 routineP = moduleN//':'//routineN 171 172 INTEGER :: i, j, k, Num_of_sub_conf, Num_of_sub_frag 173 INTEGER, DIMENSION(:, :), POINTER :: conf_loc 174 REAL(KIND=dp) :: my_energy 175 REAL(KIND=dp), DIMENSION(:), POINTER :: Em_loc 176 177 NULLIFY (conf_loc, Em_loc) 178 should_stop = .FALSE. 179 ! Count the number of subconfiguration to evaluate.. 180 Num_of_sub_frag = COUNT(conf == 1) 181 Num_of_sub_conf = 0 182 IF (Num_of_sub_frag == 1) THEN 183 CALL eval_bsse_energy_low(force_env, conf, conf, n_frags, root_section, globenv, Em) 184 ELSE 185 my_energy = 0.0_dp 186 DO k = 1, Num_of_sub_frag 187 Num_of_sub_conf = Num_of_sub_conf + & 188 FACT(Num_of_sub_frag)/(FACT(k)*FACT(Num_of_sub_frag - k)) 189 END DO 190 ALLOCATE (conf_loc(Num_of_sub_conf, Num_of_sub_frag)) 191 ALLOCATE (Em_loc(Num_of_sub_conf)) 192 Em_loc = 0.0_dp 193 CALL gen_Nbody_conf(Num_of_sub_frag, conf_loc) 194 CALL make_plan_conf(conf, conf_loc) 195 DO i = 1, Num_of_sub_conf 196 CALL eval_bsse_energy_low(force_env, conf, conf_loc(i, :), n_frags, & 197 root_section, globenv, Em_loc(i)) 198 CALL external_control(should_stop, "BSSE", globenv=globenv) 199 IF (should_stop) EXIT 200 END DO 201 ! Energy 202 k = COUNT(conf == 1) 203 DO i = 1, Num_of_sub_conf 204 j = COUNT(conf_loc(i, :) == 1) 205 my_energy = my_energy + (-1.0_dp)**(k + j)*Em_loc(i) 206 END DO 207 Em = my_energy 208 DEALLOCATE (Em_loc) 209 DEALLOCATE (conf_loc) 210 END IF 211 212 END SUBROUTINE eval_bsse_energy 213 214! ************************************************************************************************** 215!> \brief Evaluate the N-body energy contribution to the BSSE evaluation 216!> \param force_env ... 217!> \param conf ... 218!> \param conf_loc ... 219!> \param n_frags ... 220!> \param root_section ... 221!> \param globenv ... 222!> \param energy ... 223!> \par History 224!> 07.2005 created [tlaino] 225!> 2014/09/17 made atom list to be read from repeated occurance of LIST [LTong] 226!> \author Teodoro Laino 227! ************************************************************************************************** 228 SUBROUTINE eval_bsse_energy_low(force_env, conf, conf_loc, n_frags, & 229 root_section, globenv, energy) 230 TYPE(force_env_type), POINTER :: force_env 231 INTEGER, DIMENSION(:), INTENT(IN) :: conf, conf_loc 232 TYPE(section_vals_type), POINTER :: n_frags, root_section 233 TYPE(global_environment_type), POINTER :: globenv 234 REAL(KIND=dp), INTENT(OUT) :: energy 235 236 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_bsse_energy_low', & 237 routineP = moduleN//':'//routineN 238 239 CHARACTER(LEN=default_string_length) :: name 240 CHARACTER(len=default_string_length), & 241 DIMENSION(:), POINTER :: atom_type 242 INTEGER :: i, ir, isize, j, k, method_name_id, & 243 my_targ, n_rep, num_of_frag, old_size, & 244 present_charge, present_multpl 245 INTEGER, DIMENSION(:), POINTER :: atom_index, atom_list, my_conf, tmplist 246 TYPE(cell_type), POINTER :: cell 247 TYPE(cp_para_env_type), POINTER :: para_env 248 TYPE(cp_subsys_type), POINTER :: subsys, subsys_loc 249 TYPE(particle_list_type), POINTER :: particles 250 TYPE(qs_energy_type), POINTER :: qs_energy 251 TYPE(qs_environment_type), POINTER :: qs_env 252 TYPE(section_vals_type), POINTER :: bsse_section, dft_section, & 253 force_env_section, subsys_section 254 255 CALL section_vals_get(n_frags, n_repetition=num_of_frag) 256 CPASSERT(SIZE(conf) == num_of_frag) 257 NULLIFY (subsys_loc, subsys, particles, para_env, cell, atom_index, atom_type, tmplist, & 258 force_env_section) 259 CALL force_env_get(force_env, force_env_section=force_env_section) 260 CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id) 261 bsse_section => section_vals_get_subs_vals(force_env_section, "BSSE") 262 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS") 263 dft_section => section_vals_get_subs_vals(force_env_section, "DFT") 264 265 ALLOCATE (my_conf(SIZE(conf))) 266 my_conf = conf 267 CALL force_env_get(force_env=force_env, subsys=subsys, para_env=para_env, & 268 cell=cell) 269 CALL cp_subsys_get(subsys, particles=particles) 270 isize = 0 271 ALLOCATE (atom_index(isize)) 272 DO i = 1, num_of_frag 273 IF (conf(i) == 1) THEN 274 ! 275 ! Get the list of atoms creating the present fragment 276 ! 277 old_size = isize 278 CALL section_vals_val_get(n_frags, "LIST", i_rep_section=i, n_rep_val=n_rep) 279 IF (n_rep /= 0) THEN 280 DO ir = 1, n_rep 281 CALL section_vals_val_get(n_frags, "LIST", i_rep_section=i, i_rep_val=ir, i_vals=tmplist) 282 CALL reallocate(atom_index, 1, isize + SIZE(tmplist)) 283 atom_index(isize + 1:isize + SIZE(tmplist)) = tmplist 284 isize = SIZE(atom_index) 285 END DO 286 END IF 287 my_conf(i) = isize - old_size 288 CPASSERT(conf(i) /= 0) 289 END IF 290 END DO 291 CALL conf_info_setup(present_charge, present_multpl, conf, conf_loc, bsse_section, & 292 dft_section) 293 ! 294 ! Get names and modify the ghost ones 295 ! 296 ALLOCATE (atom_type(isize)) 297 DO j = 1, isize 298 my_targ = atom_index(j) 299 DO k = 1, SIZE(particles%els) 300 CALL get_atomic_kind(particles%els(k)%atomic_kind, atom_list=atom_list, name=name) 301 IF (ANY(atom_list == my_targ)) EXIT 302 END DO 303 atom_type(j) = name 304 END DO 305 DO i = 1, SIZE(conf_loc) 306 IF (my_conf(i) /= 0 .AND. conf_loc(i) == 0) THEN 307 DO j = SUM(my_conf(1:i - 1)) + 1, SUM(my_conf(1:i)) 308 atom_type(j) = TRIM(atom_type(j))//"_ghost" 309 END DO 310 END IF 311 END DO 312 CALL dump_bsse_info(atom_index, atom_type, conf, conf_loc, bsse_section, & 313 present_charge, present_multpl) 314 ! 315 ! Let's start setting up environments and calculations 316 ! 317 energy = 0.0_dp 318 IF (method_name_id == do_qs) THEN 319 CALL create_small_subsys(subsys_loc, big_subsys=subsys, & 320 small_para_env=para_env, small_cell=cell, sub_atom_index=atom_index, & 321 sub_atom_kind_name=atom_type, para_env=para_env, & 322 force_env_section=force_env_section, subsys_section=subsys_section) 323 324 CALL qs_env_create(qs_env, globenv) 325 CALL qs_init(qs_env, para_env, root_section, globenv=globenv, cp_subsys=subsys_loc, & 326 force_env_section=force_env_section, subsys_section=subsys_section, & 327 use_motion_section=.FALSE.) 328 CALL cp_subsys_release(subsys_loc) 329 330 ! 331 ! Evaluate Energy 332 ! 333 CALL qs_energies(qs_env) 334 CALL get_qs_env(qs_env, energy=qs_energy) 335 energy = qs_energy%total 336 CALL qs_env_release(qs_env) 337 ELSE 338 CPABORT("") 339 END IF 340 DEALLOCATE (atom_index) 341 DEALLOCATE (atom_type) 342 DEALLOCATE (my_conf) 343 344 END SUBROUTINE eval_bsse_energy_low 345 346! ************************************************************************************************** 347!> \brief Dumps bsse information (configuration fragment) 348!> \param atom_index ... 349!> \param atom_type ... 350!> \param conf ... 351!> \param conf_loc ... 352!> \param bsse_section ... 353!> \param present_charge ... 354!> \param present_multpl ... 355!> \par History 356!> 07.2005 created [tlaino] 357!> \author Teodoro Laino 358! ************************************************************************************************** 359 SUBROUTINE dump_bsse_info(atom_index, atom_type, conf, conf_loc, bsse_section, & 360 present_charge, present_multpl) 361 INTEGER, DIMENSION(:), POINTER :: atom_index 362 CHARACTER(len=default_string_length), & 363 DIMENSION(:), POINTER :: atom_type 364 INTEGER, DIMENSION(:), INTENT(IN) :: conf, conf_loc 365 TYPE(section_vals_type), POINTER :: bsse_section 366 INTEGER, INTENT(IN) :: present_charge, present_multpl 367 368 CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_bsse_info', routineP = moduleN//':'//routineN 369 370 INTEGER :: i, istat, iw 371 CHARACTER(len=20*SIZE(conf)) :: conf_loc_s, conf_s 372 TYPE(cp_logger_type), POINTER :: logger 373 374 NULLIFY (logger) 375 logger => cp_get_default_logger() 376 iw = cp_print_key_unit_nr(logger, bsse_section, "PRINT%PROGRAM_RUN_INFO", & 377 extension=".log") 378 IF (iw > 0) THEN 379 WRITE (conf_s, fmt="(1000I0)", iostat=istat) conf; 380 IF (istat .NE. 0) conf_s = "exceeded" 381 CALL compress(conf_s, full=.TRUE.) 382 WRITE (conf_loc_s, fmt="(1000I0)", iostat=istat) conf_loc; 383 IF (istat .NE. 0) conf_loc_s = "exceeded" 384 CALL compress(conf_loc_s, full=.TRUE.) 385 386 WRITE (UNIT=iw, FMT="(/,T2,A)") REPEAT("-", 79) 387 WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-" 388 WRITE (UNIT=iw, FMT="(T2,A,T5,A,T30,A,T55,A,T80,A)") & 389 "-", "BSSE CALCULATION", "FRAGMENT CONF: "//TRIM(conf_s), "FRAGMENT SUBCONF: "//TRIM(conf_loc_s), "-" 390 WRITE (UNIT=iw, FMT="(T2,A,T30,A,I6,T55,A,I6,T80,A)") "-", "CHARGE =", present_charge, "MULTIPLICITY =", & 391 present_multpl, "-" 392 WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-" 393 WRITE (UNIT=iw, FMT="(T2,A,T20,A,T60,A,T80,A)") "-", "ATOM INDEX", "ATOM NAME", "-" 394 WRITE (UNIT=iw, FMT="(T2,A,T20,A,T60,A,T80,A)") "-", "----------", "---------", "-" 395 DO i = 1, SIZE(atom_index) 396 WRITE (UNIT=iw, FMT="(T2,A,T20,I6,T61,A,T80,A)") "-", atom_index(i), TRIM(atom_type(i)), "-" 397 END DO 398 WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79) 399 400 CALL cp_print_key_finished_output(iw, logger, bsse_section, & 401 "PRINT%PROGRAM_RUN_INFO") 402 403 END IF 404 END SUBROUTINE dump_bsse_info 405 406! ************************************************************************************************** 407!> \brief Read modified parameters for configurations 408!> \param present_charge ... 409!> \param present_multpl ... 410!> \param conf ... 411!> \param conf_loc ... 412!> \param bsse_section ... 413!> \param dft_section ... 414!> \par History 415!> 09.2007 created [tlaino] 416!> \author Teodoro Laino - University of Zurich 417! ************************************************************************************************** 418 SUBROUTINE conf_info_setup(present_charge, present_multpl, conf, conf_loc, & 419 bsse_section, dft_section) 420 INTEGER, INTENT(OUT) :: present_charge, present_multpl 421 INTEGER, DIMENSION(:), INTENT(IN) :: conf, conf_loc 422 TYPE(section_vals_type), POINTER :: bsse_section, dft_section 423 424 CHARACTER(LEN=*), PARAMETER :: routineN = 'conf_info_setup', & 425 routineP = moduleN//':'//routineN 426 427 INTEGER :: i, nconf 428 INTEGER, DIMENSION(:), POINTER :: glb_conf, sub_conf 429 LOGICAL :: explicit 430 TYPE(section_vals_type), POINTER :: configurations 431 432 present_charge = 0 433 present_multpl = 0 434 NULLIFY (configurations, glb_conf, sub_conf) 435 ! Loop over all configurations to pick up the right one 436 configurations => section_vals_get_subs_vals(bsse_section, "CONFIGURATION") 437 CALL section_vals_get(configurations, explicit=explicit, n_repetition=nconf) 438 IF (explicit) THEN 439 DO i = 1, nconf 440 CALL section_vals_val_get(configurations, "GLB_CONF", i_rep_section=i, i_vals=glb_conf) 441 IF (SIZE(glb_conf) /= SIZE(conf)) & 442 CALL cp_abort(__LOCATION__, & 443 "GLB_CONF requires a binary description of the configuration. Number of integer "// & 444 "different from the number of fragments defined!") 445 CALL section_vals_val_get(configurations, "SUB_CONF", i_rep_section=i, i_vals=sub_conf) 446 IF (SIZE(sub_conf) /= SIZE(conf)) & 447 CALL cp_abort(__LOCATION__, & 448 "SUB_CONF requires a binary description of the configuration. Number of integer "// & 449 "different from the number of fragments defined!") 450 IF (ALL(conf == glb_conf) .AND. ALL(conf_loc == sub_conf)) THEN 451 CALL section_vals_val_get(configurations, "CHARGE", i_rep_section=i, & 452 i_val=present_charge) 453 CALL section_vals_val_get(configurations, "MULTIPLICITY", i_rep_section=i, & 454 i_val=present_multpl) 455 END IF 456 END DO 457 END IF 458 ! Setup parameter for this configuration 459 CALL section_vals_val_set(dft_section, "CHARGE", i_val=present_charge) 460 CALL section_vals_val_set(dft_section, "MULTIPLICITY", i_val=present_multpl) 461 END SUBROUTINE conf_info_setup 462 463! ************************************************************************************************** 464!> \brief Dumps results 465!> \param conf ... 466!> \param Em ... 467!> \param num_of_frag ... 468!> \param bsse_section ... 469!> \par History 470!> 09.2007 created [tlaino] 471!> \author Teodoro Laino - University of Zurich 472! ************************************************************************************************** 473 SUBROUTINE dump_bsse_results(conf, Em, num_of_frag, bsse_section) 474 INTEGER, DIMENSION(:, :), INTENT(IN) :: conf 475 REAL(KIND=dp), DIMENSION(:), POINTER :: Em 476 INTEGER, INTENT(IN) :: num_of_frag 477 TYPE(section_vals_type), POINTER :: bsse_section 478 479 CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_bsse_results', & 480 routineP = moduleN//':'//routineN 481 482 INTEGER :: i, iw 483 TYPE(cp_logger_type), POINTER :: logger 484 485 NULLIFY (logger) 486 logger => cp_get_default_logger() 487 iw = cp_print_key_unit_nr(logger, bsse_section, "PRINT%PROGRAM_RUN_INFO", & 488 extension=".log") 489 490 IF (iw > 0) THEN 491 WRITE (UNIT=iw, FMT="(/,T2,A)") REPEAT("-", 79) 492 WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-" 493 WRITE (UNIT=iw, FMT="(T2,A,T36,A,T80,A)") & 494 "-", "BSSE RESULTS", "-" 495 WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-" 496 WRITE (UNIT=iw, FMT="(T2,A,T20,A,F16.6,T80,A)") "-", "CP-corrected Total energy:", SUM(Em), "-" 497 WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-" 498 DO i = 1, SIZE(conf, 1) 499 IF (i .GT. 1) THEN 500 IF (SUM(conf(i - 1, :)) == 1 .AND. SUM(conf(i, :)) /= 1) THEN 501 WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-" 502 END IF 503 END IF 504 WRITE (UNIT=iw, FMT="(T2,A,T24,I3,A,F16.6,T80,A)") "-", SUM(conf(i, :)), "-body contribution:", Em(i), "-" 505 END DO 506 WRITE (UNIT=iw, FMT="(T2,A,T20,A,F16.6,T80,A)") "-", "BSSE-free interaction energy:", SUM(Em(Num_of_frag + 1:)), "-" 507 WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79) 508 END IF 509 510 CALL cp_print_key_finished_output(iw, logger, bsse_section, & 511 "PRINT%PROGRAM_RUN_INFO") 512 513 END SUBROUTINE dump_bsse_results 514 515! ************************************************************************************************** 516!> \brief generate the N-body configuration for the N-body BSSE evaluation 517!> \param Num_of_frag ... 518!> \param conf ... 519!> \par History 520!> 07.2005 created [tlaino] 521!> \author Teodoro Laino 522! ************************************************************************************************** 523 SUBROUTINE gen_Nbody_conf(Num_of_frag, conf) 524 INTEGER, INTENT(IN) :: Num_of_frag 525 INTEGER, DIMENSION(:, :), POINTER :: conf 526 527 INTEGER :: k, my_ind 528 529 my_ind = 0 530 ! 531 ! Set up the N-body configurations 532 ! 533 conf = 0 534 DO k = 1, Num_of_frag 535 CALL build_Nbody_conf(1, Num_of_frag, conf, k, my_ind) 536 END DO 537 END SUBROUTINE gen_Nbody_conf 538 539! ************************************************************************************************** 540!> \brief ... 541!> \param ldown ... 542!> \param lup ... 543!> \param conf ... 544!> \param k ... 545!> \param my_ind ... 546! ************************************************************************************************** 547 RECURSIVE SUBROUTINE build_Nbody_conf(ldown, lup, conf, k, my_ind) 548 INTEGER, INTENT(IN) :: ldown, lup 549 INTEGER, DIMENSION(:, :), POINTER :: conf 550 INTEGER, INTENT(IN) :: k 551 INTEGER, INTENT(INOUT) :: my_ind 552 553 INTEGER :: i, kloc, my_ind0 554 555 kloc = k - 1 556 my_ind0 = my_ind 557 IF (kloc /= 0) THEN 558 DO i = ldown, lup 559 CALL build_Nbody_conf(i + 1, lup, conf, kloc, my_ind) 560 conf(my_ind0 + 1:my_ind, i) = 1 561 my_ind0 = my_ind 562 END DO 563 ELSE 564 DO i = ldown, lup 565 my_ind = my_ind + 1 566 conf(my_ind, i) = 1 567 END DO 568 END IF 569 END SUBROUTINE build_Nbody_conf 570 571! ************************************************************************************************** 572!> \brief ... 573!> \param num ... 574!> \return ... 575! ************************************************************************************************** 576 RECURSIVE FUNCTION FACT(num) RESULT(my_fact) 577 INTEGER, INTENT(IN) :: num 578 INTEGER :: my_fact 579 580 IF (num <= 1) THEN 581 my_fact = 1 582 ELSE 583 my_fact = num*FACT(num - 1) 584 END IF 585 END FUNCTION FACT 586 587! ************************************************************************************************** 588!> \brief ... 589!> \param main_conf ... 590!> \param conf ... 591! ************************************************************************************************** 592 SUBROUTINE make_plan_conf(main_conf, conf) 593 INTEGER, DIMENSION(:), INTENT(IN) :: main_conf 594 INTEGER, DIMENSION(:, :), POINTER :: conf 595 596 INTEGER :: i, ind 597 INTEGER, DIMENSION(:, :), POINTER :: tmp_conf 598 599 ALLOCATE (tmp_conf(SIZE(conf, 1), SIZE(main_conf))) 600 tmp_conf = 0 601 ind = 0 602 DO i = 1, SIZE(main_conf) 603 IF (main_conf(i) /= 0) THEN 604 ind = ind + 1 605 tmp_conf(:, i) = conf(:, ind) 606 END IF 607 END DO 608 DEALLOCATE (conf) 609 ALLOCATE (conf(SIZE(tmp_conf, 1), SIZE(tmp_conf, 2))) 610 conf = tmp_conf 611 DEALLOCATE (tmp_conf) 612 613 END SUBROUTINE make_plan_conf 614 615! ************************************************************************************************** 616!> \brief Writes restart for BSSE calculations 617!> \param bsse_section ... 618!> \param root_section ... 619!> \par History 620!> 01.2008 created [tlaino] 621!> \author Teodoro Laino 622! ************************************************************************************************** 623 SUBROUTINE write_bsse_restart(bsse_section, root_section) 624 625 TYPE(section_vals_type), POINTER :: bsse_section, root_section 626 627 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_bsse_restart', & 628 routineP = moduleN//':'//routineN 629 630 INTEGER :: ires 631 TYPE(cp_logger_type), POINTER :: logger 632 633 logger => cp_get_default_logger() 634 ires = cp_print_key_unit_nr(logger, bsse_section, "PRINT%RESTART", & 635 extension=".restart", do_backup=.FALSE., file_position="REWIND") 636 637 IF (ires > 0) THEN 638 CALL write_restart_header(ires) 639 CALL section_vals_write(root_section, unit_nr=ires, hide_root=.TRUE.) 640 ENDIF 641 642 CALL cp_print_key_finished_output(ires, logger, bsse_section, & 643 "PRINT%RESTART") 644 645 END SUBROUTINE write_bsse_restart 646 647END MODULE bsse 648