53#include "./base/base_uses.f90"
58 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .true.
59 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'bsse'
83 INTEGER :: i, istart, k, num_of_conf, num_of_frag
84 INTEGER,
DIMENSION(:, :),
POINTER :: conf
85 LOGICAL :: explicit, should_stop
86 REAL(kind=
dp),
DIMENSION(:),
POINTER :: em
91 NULLIFY (bsse_section, n_frags, em, conf)
93 root_section => force_env%root_section
100 DO k = 1, num_of_frag
101 num_of_conf = num_of_conf + fact(num_of_frag)/(fact(k)*fact(num_of_frag - k))
103 ALLOCATE (conf(num_of_conf, num_of_frag))
104 ALLOCATE (em(num_of_conf))
105 CALL gen_nbody_conf(num_of_frag, conf)
107 should_stop = .false.
120 CALL cp_iterate(logger%iter_info, last=.false., iter_nr=istart)
123 DO i = istart + 1, num_of_conf
124 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 root_section, globenv, should_stop)
127 IF (should_stop)
EXIT
131 IF (should_stop)
EXIT
137 CALL write_bsse_restart(bsse_section, root_section)
139 IF (.NOT. should_stop)
CALL dump_bsse_results(conf, em, num_of_frag, bsse_section)
159 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
166 LOGICAL,
INTENT(OUT) :: should_stop
168 INTEGER :: i, j, k, num_of_sub_conf, num_of_sub_frag
169 INTEGER,
DIMENSION(:, :),
POINTER :: conf_loc
170 REAL(kind=
dp) :: my_energy
171 REAL(kind=
dp),
DIMENSION(:),
POINTER :: em_loc
173 NULLIFY (conf_loc, em_loc)
174 should_stop = .false.
176 num_of_sub_frag = count(conf == 1)
178 IF (num_of_sub_frag == 1)
THEN
179 CALL eval_bsse_energy_low(force_env, conf, conf, n_frags, root_section, globenv, em)
182 DO k = 1, num_of_sub_frag
183 num_of_sub_conf = num_of_sub_conf + &
184 fact(num_of_sub_frag)/(fact(k)*fact(num_of_sub_frag - k))
186 ALLOCATE (conf_loc(num_of_sub_conf, num_of_sub_frag))
187 ALLOCATE (em_loc(num_of_sub_conf))
189 CALL gen_nbody_conf(num_of_sub_frag, conf_loc)
190 CALL make_plan_conf(conf, conf_loc)
191 DO i = 1, num_of_sub_conf
192 CALL eval_bsse_energy_low(force_env, conf, conf_loc(i, :), n_frags, &
193 root_section, globenv, em_loc(i))
195 IF (should_stop)
EXIT
199 DO i = 1, num_of_sub_conf
200 j = count(conf_loc(i, :) == 1)
201 my_energy = my_energy + (-1.0_dp)**(k + j)*em_loc(i)
205 DEALLOCATE (conf_loc)
208 END SUBROUTINE eval_bsse_energy
224 SUBROUTINE eval_bsse_energy_low(force_env, conf, conf_loc, n_frags, &
225 root_section, globenv, energy)
227 INTEGER,
DIMENSION(:),
INTENT(IN) :: conf, conf_loc
230 REAL(kind=
dp),
INTENT(OUT) :: energy
232 CHARACTER(LEN=default_string_length) :: name
233 CHARACTER(len=default_string_length), &
234 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 INTEGER,
DIMENSION(:),
POINTER :: atom_index, atom_list, my_conf, tmplist
244 force_env_section, subsys_section
247 cpassert(
SIZE(conf) == num_of_frag)
248 NULLIFY (subsys, particles, para_env, cell, atom_index, atom_type, tmplist, &
250 CALL force_env_get(force_env, force_env_section=force_env_section)
256 ALLOCATE (my_conf(
SIZE(conf)))
258 CALL force_env_get(force_env=force_env, subsys=subsys, para_env=para_env, &
262 ALLOCATE (atom_index(isize))
263 DO i = 1, num_of_frag
264 IF (conf(i) == 1)
THEN
273 CALL reallocate(atom_index, 1, isize +
SIZE(tmplist))
274 atom_index(isize + 1:isize +
SIZE(tmplist)) = tmplist
275 isize =
SIZE(atom_index)
278 my_conf(i) = isize - old_size
279 cpassert(conf(i) /= 0)
282 CALL conf_info_setup(present_charge, present_multpl, conf, conf_loc, bsse_section, &
287 ALLOCATE (atom_type(isize))
289 my_targ = atom_index(j)
290 DO k = 1,
SIZE(particles%els)
291 CALL get_atomic_kind(particles%els(k)%atomic_kind, atom_list=atom_list, name=name)
292 IF (any(atom_list == my_targ))
EXIT
296 DO i = 1,
SIZE(conf_loc)
297 IF (my_conf(i) /= 0 .AND. conf_loc(i) == 0)
THEN
298 DO j = sum(my_conf(1:i - 1)) + 1, sum(my_conf(1:i))
299 atom_type(j) = trim(atom_type(j))//
"_ghost"
303 CALL dump_bsse_info(atom_index, atom_type, conf, conf_loc, bsse_section, &
304 present_charge, present_multpl)
309 IF (method_name_id ==
do_qs)
THEN
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 force_env_section=force_env_section, subsys_section=subsys_section)
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 use_motion_section=.false.)
339 DEALLOCATE (atom_index)
340 DEALLOCATE (atom_type)
343 END SUBROUTINE eval_bsse_energy_low
358 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
367 INTEGER :: i, istat, iw
368 CHARACTER(len=20*SIZE(conf)) :: conf_loc_s, conf_s
369 TYPE(cp_logger_type),
POINTER :: logger
372 logger => cp_get_default_logger()
373 iw = cp_print_key_unit_nr(logger, bsse_section,
"PRINT%PROGRAM_RUN_INFO", &
376 WRITE (conf_s, fmt=
"(1000I0)", iostat=istat) conf;
377 IF (istat .NE. 0) conf_s =
"exceeded"
378 CALL compress(conf_s, full=.true.)
379 WRITE (conf_loc_s, fmt=
"(1000I0)", iostat=istat) conf_loc;
380 IF (istat .NE. 0) conf_loc_s =
"exceeded"
381 CALL compress(conf_loc_s, full=.true.)
383 WRITE (unit=iw, fmt=
"(/,T2,A)") repeat(
"-", 79)
384 WRITE (unit=iw, fmt=
"(T2,A,T80,A)")
"-",
"-"
385 WRITE (unit=iw, fmt=
"(T2,A,T5,A,T30,A,T55,A,T80,A)") &
386 "-",
"BSSE CALCULATION",
"FRAGMENT CONF: "//trim(conf_s),
"FRAGMENT SUBCONF: "//trim(conf_loc_s),
"-"
387 WRITE (unit=iw, fmt=
"(T2,A,T30,A,I6,T55,A,I6,T80,A)")
"-",
"CHARGE =", present_charge,
"MULTIPLICITY =", &
389 WRITE (unit=iw, fmt=
"(T2,A,T80,A)")
"-",
"-"
390 WRITE (unit=iw, fmt=
"(T2,A,T20,A,T60,A,T80,A)")
"-",
"ATOM INDEX",
"ATOM NAME",
"-"
391 WRITE (unit=iw, fmt=
"(T2,A,T20,A,T60,A,T80,A)")
"-",
"----------",
"---------",
"-"
392 DO i = 1,
SIZE(atom_index)
393 WRITE (unit=iw, fmt=
"(T2,A,T20,I6,T61,A,T80,A)")
"-", atom_index(i), trim(atom_type(i)),
"-"
395 WRITE (unit=iw, fmt=
"(T2,A)") repeat(
"-", 79)
397 CALL cp_print_key_finished_output(iw, logger, bsse_section, &
398 "PRINT%PROGRAM_RUN_INFO")
401 END SUBROUTINE dump_bsse_info
415 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
422 INTEGER,
DIMENSION(:),
POINTER :: glb_conf, sub_conf
424 TYPE(section_vals_type),
POINTER :: configurations
428 NULLIFY (configurations, glb_conf, sub_conf)
430 configurations => section_vals_get_subs_vals(bsse_section,
"CONFIGURATION")
431 CALL section_vals_get(configurations, explicit=explicit, n_repetition=nconf)
434 CALL section_vals_val_get(configurations,
"GLB_CONF", i_rep_section=i, i_vals=glb_conf)
435 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 "different from the number of fragments defined!")
439 CALL section_vals_val_get(configurations,
"SUB_CONF", i_rep_section=i, i_vals=sub_conf)
440 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 "different from the number of fragments defined!")
444 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 i_val=present_charge)
447 CALL section_vals_val_get(configurations,
"MULTIPLICITY", i_rep_section=i, &
448 i_val=present_multpl)
453 CALL section_vals_val_set(dft_section,
"CHARGE", i_val=present_charge)
454 CALL section_vals_val_set(dft_section,
"MULTIPLICITY", i_val=present_multpl)
455 END SUBROUTINE conf_info_setup
467 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
474 TYPE(cp_logger_type),
POINTER :: logger
477 logger => cp_get_default_logger()
478 iw = cp_print_key_unit_nr(logger, bsse_section,
"PRINT%PROGRAM_RUN_INFO", &
482 WRITE (unit=iw, fmt=
"(/,T2,A)") repeat(
"-", 79)
483 WRITE (unit=iw, fmt=
"(T2,A,T80,A)")
"-",
"-"
484 WRITE (unit=iw, fmt=
"(T2,A,T36,A,T80,A)") &
485 "-",
"BSSE RESULTS",
"-"
486 WRITE (unit=iw, fmt=
"(T2,A,T80,A)")
"-",
"-"
487 WRITE (unit=iw, fmt=
"(T2,A,T20,A,F16.6,T80,A)")
"-",
"CP-corrected Total energy:", sum(em),
"-"
488 WRITE (unit=iw, fmt=
"(T2,A,T80,A)")
"-",
"-"
489 DO i = 1,
SIZE(conf, 1)
491 IF (sum(conf(i - 1, :)) == 1 .AND. sum(conf(i, :)) /= 1)
THEN
492 WRITE (unit=iw, fmt=
"(T2,A,T80,A)")
"-",
"-"
495 WRITE (unit=iw, fmt=
"(T2,A,T24,I3,A,F16.6,T80,A)")
"-", sum(conf(i, :)),
"-body contribution:", em(i),
"-"
497 WRITE (unit=iw, fmt=
"(T2,A,T20,A,F16.6,T80,A)")
"-",
"BSSE-free interaction energy:", sum(em(num_of_frag + 1:)),
"-"
498 WRITE (unit=iw, fmt=
"(T2,A)") repeat(
"-", 79)
501 CALL cp_print_key_finished_output(iw, logger, bsse_section, &
502 "PRINT%PROGRAM_RUN_INFO")
504 END SUBROUTINE dump_bsse_results
514 SUBROUTINE gen_nbody_conf(Num_of_frag, conf)
515 INTEGER,
INTENT(IN) :: num_of_frag
516 INTEGER,
DIMENSION(:, :),
POINTER :: conf
525 DO k = 1, num_of_frag
526 CALL build_nbody_conf(1, num_of_frag, conf, k, my_ind)
528 END SUBROUTINE gen_nbody_conf
538 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
544 INTEGER :: i, kloc, my_ind0
550 CALL build_nbody_conf(i + 1, lup, conf, kloc, my_ind)
551 conf(my_ind0 + 1:my_ind, i) = 1
560 END SUBROUTINE build_nbody_conf
567 RECURSIVE FUNCTION fact(num)
RESULT(my_fact)
568 INTEGER,
INTENT(IN) :: num
574 my_fact = num*fact(num - 1)
583 SUBROUTINE make_plan_conf(main_conf, conf)
584 INTEGER,
DIMENSION(:),
INTENT(IN) :: main_conf
585 INTEGER,
DIMENSION(:, :),
POINTER :: conf
588 INTEGER,
DIMENSION(:, :),
POINTER :: tmp_conf
590 ALLOCATE (tmp_conf(
SIZE(conf, 1),
SIZE(main_conf)))
593 DO i = 1,
SIZE(main_conf)
594 IF (main_conf(i) /= 0)
THEN
596 tmp_conf(:, i) = conf(:, ind)
600 ALLOCATE (conf(
SIZE(tmp_conf, 1),
SIZE(tmp_conf, 2)))
602 DEALLOCATE (tmp_conf)
604 END SUBROUTINE make_plan_conf
614 SUBROUTINE write_bsse_restart(bsse_section, root_section)
616 TYPE(section_vals_type),
POINTER :: bsse_section, root_section
619 TYPE(cp_logger_type),
POINTER :: logger
621 logger => cp_get_default_logger()
622 ires = cp_print_key_unit_nr(logger, bsse_section,
"PRINT%RESTART", &
623 extension=
".restart", do_backup=.false., file_position=
"REWIND")
626 CALL write_restart_header(ires)
627 CALL section_vals_write(root_section, unit_nr=ires, hide_root=.true.)
630 CALL cp_print_key_finished_output(ires, logger, bsse_section, &
633 END SUBROUTINE write_bsse_restart
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Module to perform a counterpoise correction (BSSE)
subroutine, public do_bsse_calculation(force_env, globenv)
Perform an COUNTERPOISE CORRECTION (BSSE) For a 2-body system the correction scheme can be represente...
Handles all functions related to the CELL.
some minimal info about CP2K, including its version and license
subroutine, public write_restart_header(iunit)
Writes the header for the restart file.
Routines to handle the external control of CP2K.
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
subroutine, public cp_rm_iter_level(iteration_info, level_name, n_rlevel_att)
Removes an iteration level.
subroutine, public cp_add_iter_level(iteration_info, level_name, n_rlevel_new)
Adds an iteration level.
Initialize a small environment for a particular calculation.
subroutine, public create_small_subsys(small_subsys, big_subsys, small_cell, small_para_env, sub_atom_index, sub_atom_kind_name, para_env, force_env_section, subsys_section, ignore_outside_box)
updates the molecule information of the given subsys
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_release(subsys)
releases a subsys (see doc/ReferenceCounting.html)
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Define type storing the global information of a run. Keep the amount of stored data small....
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Utility routines for the memory handling.
Interface to the message passing library MPI.
represent a simple array based list of the given type
Perform a QUICKSTEP wavefunction optimization (single point)
subroutine, public qs_energies(qs_env, consistent_energies, calc_forces)
Driver routine for QUICKSTEP single point wavefunction optimization.
subroutine, public qs_env_release(qs_env)
releases the given qs_env (see doc/ReferenceCounting.html)
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
subroutine, public qs_env_create(qs_env, globenv)
allocates and intitializes a qs_env
subroutine, public qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, cell, cell_ref, qmmm, qmmm_env_qm, force_env_section, subsys_section, use_motion_section, silent)
Read the input and the database files for the setup of the QUICKSTEP environment.
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Type defining parameters related to the simulation cell.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represents a system: atoms, molecules, their pos,vel,...
wrapper to abstract the force evaluation of the various methods
contains the initially parsed file and the initial parallel environment
stores all the informations relevant to an mpi environment
represent a list of objects