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