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