LCOV - code coverage report
Current view: top level - src - bsse.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.7 % 238 235
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 11 11

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

Generated by: LCOV version 2.0-1