Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 10 : 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 10 : INTEGER, DIMENSION(:, :), POINTER :: conf
85 : LOGICAL :: explicit, should_stop
86 10 : 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 10 : NULLIFY (bsse_section, n_frags, Em, conf)
92 20 : logger => cp_get_default_logger()
93 10 : root_section => force_env%root_section
94 10 : bsse_section => section_vals_get_subs_vals(force_env%force_env_section, "BSSE")
95 10 : n_frags => section_vals_get_subs_vals(bsse_section, "FRAGMENT")
96 10 : CALL section_vals_get(n_frags, n_repetition=Num_of_Frag)
97 :
98 : ! Number of configurations
99 10 : num_of_conf = 0
100 34 : DO k = 1, Num_of_frag
101 34 : num_of_conf = num_of_conf + FACT(Num_of_frag)/(FACT(k)*FACT(Num_of_frag - k))
102 : END DO
103 40 : ALLOCATE (conf(num_of_conf, Num_of_frag))
104 30 : ALLOCATE (Em(num_of_conf))
105 10 : CALL gen_Nbody_conf(Num_of_frag, conf)
106 :
107 10 : should_stop = .FALSE.
108 10 : istart = 0
109 10 : fragment_energies_section => section_vals_get_subs_vals(bsse_section, "FRAGMENT_ENERGIES")
110 10 : CALL section_vals_get(fragment_energies_section, explicit=explicit)
111 10 : 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 10 : CALL cp_add_iter_level(logger%iter_info, "BSSE")
120 10 : CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=istart)
121 :
122 : ! Evaluating the energy of the N-body cluster terms
123 52 : DO i = istart + 1, num_of_conf
124 42 : 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 42 : root_section, globenv, should_stop)
127 42 : IF (should_stop) EXIT
128 :
129 : ! If no signal was received in the inner loop let's check also at this stage
130 42 : CALL external_control(should_stop, "BSSE", globenv=globenv)
131 42 : 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 42 : i_rep_val=i)
137 94 : CALL write_bsse_restart(bsse_section, root_section)
138 : END DO
139 10 : IF (.NOT. should_stop) CALL dump_bsse_results(conf, Em, num_of_frag, bsse_section)
140 10 : CALL cp_rm_iter_level(logger%iter_info, "BSSE")
141 10 : DEALLOCATE (Em)
142 10 : DEALLOCATE (conf)
143 :
144 30 : 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 42 : 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 42 : INTEGER, DIMENSION(:, :), POINTER :: conf_loc
170 : REAL(KIND=dp) :: my_energy
171 42 : REAL(KIND=dp), DIMENSION(:), POINTER :: Em_loc
172 :
173 42 : NULLIFY (conf_loc, Em_loc)
174 42 : should_stop = .FALSE.
175 : ! Count the number of subconfiguration to evaluate..
176 154 : Num_of_sub_frag = COUNT(conf == 1)
177 42 : Num_of_sub_conf = 0
178 42 : IF (Num_of_sub_frag == 1) THEN
179 20 : CALL eval_bsse_energy_low(force_env, conf, conf, n_frags, root_section, globenv, Em)
180 : ELSE
181 22 : my_energy = 0.0_dp
182 70 : DO k = 1, Num_of_sub_frag
183 : Num_of_sub_conf = Num_of_sub_conf + &
184 70 : FACT(Num_of_sub_frag)/(FACT(k)*FACT(Num_of_sub_frag - k))
185 : END DO
186 88 : ALLOCATE (conf_loc(Num_of_sub_conf, Num_of_sub_frag))
187 66 : ALLOCATE (Em_loc(Num_of_sub_conf))
188 104 : Em_loc = 0.0_dp
189 22 : CALL gen_Nbody_conf(Num_of_sub_frag, conf_loc)
190 22 : CALL make_plan_conf(conf, conf_loc)
191 104 : DO i = 1, Num_of_sub_conf
192 : CALL eval_bsse_energy_low(force_env, conf, conf_loc(i, :), n_frags, &
193 82 : root_section, globenv, Em_loc(i))
194 82 : CALL external_control(should_stop, "BSSE", globenv=globenv)
195 104 : IF (should_stop) EXIT
196 : END DO
197 : ! Energy
198 82 : k = COUNT(conf == 1)
199 104 : DO i = 1, Num_of_sub_conf
200 310 : j = COUNT(conf_loc(i, :) == 1)
201 104 : my_energy = my_energy + (-1.0_dp)**(k + j)*Em_loc(i)
202 : END DO
203 22 : Em = my_energy
204 22 : DEALLOCATE (Em_loc)
205 22 : DEALLOCATE (conf_loc)
206 : END IF
207 :
208 42 : 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 102 : 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 102 : 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 102 : 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 102 : CALL section_vals_get(n_frags, n_repetition=num_of_frag)
247 102 : CPASSERT(SIZE(conf) == num_of_frag)
248 102 : NULLIFY (subsys, particles, para_env, cell, atom_index, atom_type, tmplist, &
249 102 : force_env_section)
250 102 : CALL force_env_get(force_env, force_env_section=force_env_section)
251 102 : CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
252 102 : bsse_section => section_vals_get_subs_vals(force_env_section, "BSSE")
253 102 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
254 102 : dft_section => section_vals_get_subs_vals(force_env_section, "DFT")
255 :
256 306 : ALLOCATE (my_conf(SIZE(conf)))
257 382 : my_conf = conf
258 : CALL force_env_get(force_env=force_env, subsys=subsys, para_env=para_env, &
259 102 : cell=cell)
260 102 : CALL cp_subsys_get(subsys, particles=particles)
261 102 : isize = 0
262 102 : ALLOCATE (atom_index(isize))
263 382 : DO i = 1, num_of_frag
264 382 : IF (conf(i) == 1) THEN
265 : !
266 : ! Get the list of atoms creating the present fragment
267 : !
268 212 : old_size = isize
269 212 : CALL section_vals_val_get(n_frags, "LIST", i_rep_section=i, n_rep_val=n_rep)
270 212 : IF (n_rep /= 0) THEN
271 508 : DO ir = 1, n_rep
272 296 : CALL section_vals_val_get(n_frags, "LIST", i_rep_section=i, i_rep_val=ir, i_vals=tmplist)
273 296 : CALL reallocate(atom_index, 1, isize + SIZE(tmplist))
274 1848 : atom_index(isize + 1:isize + SIZE(tmplist)) = tmplist
275 508 : isize = SIZE(atom_index)
276 : END DO
277 : END IF
278 212 : my_conf(i) = isize - old_size
279 212 : 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 102 : dft_section)
284 : !
285 : ! Get names and modify the ghost ones
286 : !
287 306 : ALLOCATE (atom_type(isize))
288 730 : DO j = 1, isize
289 628 : my_targ = atom_index(j)
290 1068 : DO k = 1, SIZE(particles%els)
291 1068 : CALL get_atomic_kind(particles%els(k)%atomic_kind, atom_list=atom_list, name=name)
292 3422 : IF (ANY(atom_list == my_targ)) EXIT
293 : END DO
294 730 : atom_type(j) = name
295 : END DO
296 382 : DO i = 1, SIZE(conf_loc)
297 382 : IF (my_conf(i) /= 0 .AND. conf_loc(i) == 0) THEN
298 562 : DO j = SUM(my_conf(1:i - 1)) + 1, SUM(my_conf(1:i))
299 286 : 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 102 : present_charge, present_multpl)
305 : !
306 : ! Let's start setting up environments and calculations
307 : !
308 102 : energy = 0.0_dp
309 102 : IF (method_name_id == do_qs) THEN
310 : 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 102 : 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 102 : force_env_section=force_env_section, subsys_section=subsys_section)
319 :
320 102 : ALLOCATE (qs_env)
321 102 : 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 102 : use_motion_section=.FALSE.)
325 102 : CALL cp_subsys_release(subsys_loc)
326 :
327 : !
328 : ! Evaluate Energy
329 : !
330 102 : CALL qs_energies(qs_env)
331 102 : CALL get_qs_env(qs_env, energy=qs_energy)
332 102 : energy = qs_energy%total
333 102 : CALL qs_env_release(qs_env)
334 102 : DEALLOCATE (qs_env)
335 : END BLOCK
336 : ELSE
337 0 : CPABORT("")
338 : END IF
339 102 : DEALLOCATE (atom_index)
340 102 : DEALLOCATE (atom_type)
341 102 : DEALLOCATE (my_conf)
342 :
343 204 : 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 102 : 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 102 : CHARACTER(len=20*SIZE(conf)) :: conf_loc_s, conf_s
369 : TYPE(cp_logger_type), POINTER :: logger
370 :
371 102 : NULLIFY (logger)
372 102 : logger => cp_get_default_logger()
373 : iw = cp_print_key_unit_nr(logger, bsse_section, "PRINT%PROGRAM_RUN_INFO", &
374 102 : extension=".log")
375 102 : IF (iw > 0) THEN
376 51 : WRITE (conf_s, fmt="(1000I0)", iostat=istat) conf;
377 51 : IF (istat .NE. 0) conf_s = "exceeded"
378 51 : CALL compress(conf_s, full=.TRUE.)
379 51 : WRITE (conf_loc_s, fmt="(1000I0)", iostat=istat) conf_loc;
380 51 : IF (istat .NE. 0) conf_loc_s = "exceeded"
381 51 : CALL compress(conf_loc_s, full=.TRUE.)
382 :
383 51 : WRITE (UNIT=iw, FMT="(/,T2,A)") REPEAT("-", 79)
384 51 : WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-"
385 : WRITE (UNIT=iw, FMT="(T2,A,T5,A,T30,A,T55,A,T80,A)") &
386 51 : "-", "BSSE CALCULATION", "FRAGMENT CONF: "//TRIM(conf_s), "FRAGMENT SUBCONF: "//TRIM(conf_loc_s), "-"
387 51 : WRITE (UNIT=iw, FMT="(T2,A,T30,A,I6,T55,A,I6,T80,A)") "-", "CHARGE =", present_charge, "MULTIPLICITY =", &
388 102 : present_multpl, "-"
389 51 : WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-"
390 51 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,T60,A,T80,A)") "-", "ATOM INDEX", "ATOM NAME", "-"
391 51 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,T60,A,T80,A)") "-", "----------", "---------", "-"
392 365 : DO i = 1, SIZE(atom_index)
393 365 : WRITE (UNIT=iw, FMT="(T2,A,T20,I6,T61,A,T80,A)") "-", atom_index(i), TRIM(atom_type(i)), "-"
394 : END DO
395 51 : WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79)
396 :
397 : CALL cp_print_key_finished_output(iw, logger, bsse_section, &
398 51 : "PRINT%PROGRAM_RUN_INFO")
399 :
400 : END IF
401 102 : 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 102 : 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 102 : INTEGER, DIMENSION(:), POINTER :: glb_conf, sub_conf
423 : LOGICAL :: explicit
424 : TYPE(section_vals_type), POINTER :: configurations
425 :
426 102 : present_charge = 0
427 102 : present_multpl = 0
428 102 : NULLIFY (configurations, glb_conf, sub_conf)
429 : ! Loop over all configurations to pick up the right one
430 204 : configurations => section_vals_get_subs_vals(bsse_section, "CONFIGURATION")
431 102 : CALL section_vals_get(configurations, explicit=explicit, n_repetition=nconf)
432 102 : 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 102 : CALL section_vals_val_set(dft_section, "CHARGE", i_val=present_charge)
454 102 : CALL section_vals_val_set(dft_section, "MULTIPLICITY", i_val=present_multpl)
455 102 : 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 10 : 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 10 : NULLIFY (logger)
477 10 : logger => cp_get_default_logger()
478 : iw = cp_print_key_unit_nr(logger, bsse_section, "PRINT%PROGRAM_RUN_INFO", &
479 10 : extension=".log")
480 :
481 10 : IF (iw > 0) THEN
482 5 : WRITE (UNIT=iw, FMT="(/,T2,A)") REPEAT("-", 79)
483 5 : WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-"
484 : WRITE (UNIT=iw, FMT="(T2,A,T36,A,T80,A)") &
485 5 : "-", "BSSE RESULTS", "-"
486 5 : WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-"
487 28 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,F16.6,T80,A)") "-", "CP-corrected Total energy:", SUM(Em), "-"
488 5 : WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-"
489 28 : DO i = 1, SIZE(conf, 1)
490 23 : IF (i .GT. 1) THEN
491 114 : IF (SUM(conf(i - 1, :)) == 1 .AND. SUM(conf(i, :)) /= 1) THEN
492 5 : WRITE (UNIT=iw, FMT="(T2,A,T80,A)") "-", "-"
493 : END IF
494 : END IF
495 88 : WRITE (UNIT=iw, FMT="(T2,A,T24,I3,A,F16.6,T80,A)") "-", SUM(conf(i, :)), "-body contribution:", Em(i), "-"
496 : END DO
497 16 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,F16.6,T80,A)") "-", "BSSE-free interaction energy:", SUM(Em(Num_of_frag + 1:)), "-"
498 5 : WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79)
499 : END IF
500 :
501 : CALL cp_print_key_finished_output(iw, logger, bsse_section, &
502 10 : "PRINT%PROGRAM_RUN_INFO")
503 :
504 10 : 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 32 : 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 32 : my_ind = 0
521 : !
522 : ! Set up the N-body configurations
523 : !
524 416 : conf = 0
525 104 : DO k = 1, Num_of_frag
526 104 : CALL build_Nbody_conf(1, Num_of_frag, conf, k, my_ind)
527 : END DO
528 32 : 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 192 : 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 192 : kloc = k - 1
547 192 : my_ind0 = my_ind
548 192 : IF (kloc /= 0) THEN
549 184 : DO i = ldown, lup
550 120 : CALL build_Nbody_conf(i + 1, lup, conf, kloc, my_ind)
551 184 : conf(my_ind0 + 1:my_ind, i) = 1
552 64 : my_ind0 = my_ind
553 : END DO
554 : ELSE
555 256 : DO i = ldown, lup
556 128 : my_ind = my_ind + 1
557 256 : conf(my_ind, i) = 1
558 : END DO
559 : END IF
560 192 : END SUBROUTINE build_Nbody_conf
561 :
562 : ! **************************************************************************************************
563 : !> \brief ...
564 : !> \param num ...
565 : !> \return ...
566 : ! **************************************************************************************************
567 368 : RECURSIVE FUNCTION FACT(num) RESULT(my_fact)
568 : INTEGER, INTENT(IN) :: num
569 : INTEGER :: my_fact
570 :
571 368 : IF (num <= 1) THEN
572 : my_fact = 1
573 : ELSE
574 152 : my_fact = num*FACT(num - 1)
575 : END IF
576 368 : END FUNCTION FACT
577 :
578 : ! **************************************************************************************************
579 : !> \brief ...
580 : !> \param main_conf ...
581 : !> \param conf ...
582 : ! **************************************************************************************************
583 22 : 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 88 : ALLOCATE (tmp_conf(SIZE(conf, 1), SIZE(main_conf)))
591 310 : tmp_conf = 0
592 82 : ind = 0
593 82 : DO i = 1, SIZE(main_conf)
594 82 : IF (main_conf(i) /= 0) THEN
595 48 : ind = ind + 1
596 480 : tmp_conf(:, i) = conf(:, ind)
597 : END IF
598 : END DO
599 22 : DEALLOCATE (conf)
600 88 : ALLOCATE (conf(SIZE(tmp_conf, 1), SIZE(tmp_conf, 2)))
601 620 : conf = tmp_conf
602 22 : DEALLOCATE (tmp_conf)
603 :
604 22 : 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 42 : 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 42 : logger => cp_get_default_logger()
622 : ires = cp_print_key_unit_nr(logger, bsse_section, "PRINT%RESTART", &
623 42 : extension=".restart", do_backup=.FALSE., file_position="REWIND")
624 :
625 42 : IF (ires > 0) THEN
626 21 : CALL write_restart_header(ires)
627 21 : 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 42 : "PRINT%RESTART")
632 :
633 42 : END SUBROUTINE write_bsse_restart
634 :
635 : END MODULE bsse
|
