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 Molecular and Gamma-point periodic CC4S export from restricted RI-HF.
10 : !> RI_FLAVOR MO stores a weighted AO-pair triangle. The full vertex is
11 : !> C^T (T + T^T) C, dressed with the symmetric Coulomb inverse square root.
12 : !> The SCF tensors and orbitals are read only. MO transformations are distributed
13 : !> and bounded output pieces are gathered to the writing rank.
14 : ! **************************************************************************************************
15 : MODULE cc4s_export
16 : USE ISO_FORTRAN_ENV, ONLY: file_storage_size
17 : USE cell_types, ONLY: cell_type,&
18 : plane_distance,&
19 : use_perd_none,&
20 : use_perd_xyz
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
23 : dbcsr_release,&
24 : dbcsr_type
25 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26 : cp_fm_to_dbcsr_row_template
27 : USE cp_fm_diag, ONLY: cp_fm_geeig_canon
28 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
29 : cp_fm_struct_release,&
30 : cp_fm_struct_type
31 : USE cp_fm_types, ONLY: cp_fm_create,&
32 : cp_fm_release,&
33 : cp_fm_to_fm,&
34 : cp_fm_type
35 : USE cp_log_handling, ONLY: cp_get_default_logger,&
36 : cp_logger_get_default_io_unit,&
37 : cp_logger_type
38 : USE cp_output_handling, ONLY: cp_p_file,&
39 : cp_print_key_finished_output,&
40 : cp_print_key_generate_filename,&
41 : cp_print_key_should_output,&
42 : cp_print_key_unit_nr,&
43 : cp_printkey_is_on
44 : USE dbt_api, ONLY: &
45 : dbt_contract, dbt_copy, dbt_copy_matrix_to_tensor, dbt_create, dbt_destroy, dbt_get_block, &
46 : dbt_get_info, dbt_iterator_blocks_left, dbt_iterator_next_block, dbt_iterator_start, &
47 : dbt_iterator_stop, dbt_iterator_type, dbt_nblks_total, dbt_pgrid_create, &
48 : dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
49 : USE hfx_types, ONLY: hfx_type
50 : USE ieee_arithmetic, ONLY: ieee_is_finite
51 : USE input_constants, ONLY: do_potential_coulomb,&
52 : do_potential_truncated,&
53 : sic_none,&
54 : xc_none
55 : USE input_cp2k_hfx, ONLY: ri_mo
56 : USE input_section_types, ONLY: section_vals_get,&
57 : section_vals_get_subs_vals,&
58 : section_vals_type,&
59 : section_vals_val_get
60 : USE kinds, ONLY: default_path_length,&
61 : dp,&
62 : dp_size,&
63 : int_8
64 : USE machine, ONLY: m_walltime
65 : USE message_passing, ONLY: mp_para_env_type
66 : USE qs_environment_types, ONLY: get_qs_env,&
67 : qs_environment_type
68 : USE qs_mo_types, ONLY: get_mo_set,&
69 : has_uniform_occupation,&
70 : mo_set_type
71 : USE qs_tensors_types, ONLY: create_2c_tensor,&
72 : create_3c_tensor
73 : USE scf_control_types, ONLY: scf_control_type
74 : #include "./base/base_uses.f90"
75 :
76 : IMPLICIT NONE
77 : PRIVATE
78 : PUBLIC :: write_cc4s
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cc4s_export'
80 : ! MPI statistics count bytes in default integers. DBT/DBM also use default-integer
81 : ! counts and grow buffers by up to 30%; leave a factor of two for data storage.
82 : INTEGER(int_8), PARAMETER :: output_count_limit = INT(HUGE(0), int_8)/dp_size
83 : INTEGER(int_8), PARAMETER :: backend_data_limit = output_count_limit/2
84 : ! Up to four 64-bit indices per block, with allocation/communication headroom.
85 : INTEGER(int_8), PARAMETER :: backend_block_limit = INT(HUGE(0), int_8)/64
86 : REAL(dp), PARAMETER :: elements_per_mib = 1024.0_dp**2/REAL(dp_size, dp)
87 :
88 : CONTAINS
89 :
90 : ! **************************************************************************************************
91 : !> \brief Check the requested export, canonicalize private Fock/overlap copies and write tensors.
92 : !> \param qs_env QS environment after SCF
93 : !> \param converged Whether the SCF actually converged (including the outer SCF)
94 : ! **************************************************************************************************
95 24709 : SUBROUTINE write_cc4s(qs_env, converged)
96 : TYPE(qs_environment_type), POINTER :: qs_env
97 : LOGICAL, INTENT(IN) :: converged
98 :
99 : CHARACTER(len=*), PARAMETER :: routineN = 'write_cc4s'
100 : REAL(dp), PARAMETER :: hf_fraction_tolerance = 1.0e-14_dp
101 :
102 : CHARACTER(len=8) :: method
103 : CHARACTER(len=default_path_length) :: filename
104 : INTEGER :: frozen, functional, handle, homo, max_iter, method_id, nao, nmo, nri, nstored, &
105 : output_unit, poisson_periodic, stage_handle
106 : INTEGER(KIND=int_8) :: buffer_elements, total_elements
107 : INTEGER, DIMENSION(2) :: dims2
108 : INTEGER, DIMENSION(3) :: dims3
109 : LOGICAL :: admm_section_present, do_kpoints, &
110 : enabled, periodic
111 : REAL(dp) :: amp_conv, cutoff_radius, energy_conv, export_start, fermi_energy, max_cutoff, &
112 : max_memory_mb, memory_mb, output_buffer_mib, vertex_norm, write_start
113 24709 : REAL(dp), ALLOCATABLE :: eigenvalues(:)
114 : REAL(dp), DIMENSION(5) :: elapsed
115 : TYPE(cell_type), POINTER :: cell
116 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
117 : TYPE(cp_fm_type) :: fm_active, fm_coeff, fm_fock, &
118 : fm_overlap, fm_work
119 : TYPE(cp_fm_type), POINTER :: mo_coeff
120 : TYPE(cp_logger_type), POINTER :: logger
121 24709 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
122 : TYPE(dbcsr_type) :: coefficient_matrix
123 172963 : TYPE(dbt_type) :: coefficients
124 : TYPE(dft_control_type), POINTER :: dft_control
125 24709 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
126 24709 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
127 : TYPE(mp_para_env_type), POINTER :: para_env
128 : TYPE(scf_control_type), POINTER :: scf_control
129 : TYPE(section_vals_type), POINTER :: admm_section, input, print_section
130 :
131 49418 : logger => cp_get_default_logger()
132 24709 : CALL get_qs_env(qs_env, input=input)
133 24709 : print_section => section_vals_get_subs_vals(input, "DFT%PRINT%CC4S")
134 24709 : enabled = cp_printkey_is_on(logger%iter_info, print_section)
135 24709 : IF (.NOT. enabled) RETURN
136 6 : IF (.NOT. qs_env%single_point_run) THEN
137 0 : CALL cp_abort(__LOCATION__, "CC4S: export requires a single-point RUN_TYPE ENERGY or ENERGY_FORCE")
138 : END IF
139 6 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, print_section), cp_p_file)) RETURN
140 6 : CALL section_vals_val_get(print_section, "FILENAME", c_val=filename)
141 6 : IF (filename(1:1) == '=' .OR. filename == '__STD_OUT__') THEN
142 : CALL cp_abort(__LOCATION__, &
143 0 : "CC4S: FILENAME must be a file prefix, not an exact name (=name) or __STD_OUT__")
144 : END IF
145 6 : CALL timeset(routineN, handle)
146 6 : export_start = m_walltime()
147 6 : elapsed = 0.0_dp
148 6 : output_unit = cp_logger_get_default_io_unit(logger)
149 :
150 6 : CALL section_vals_val_get(print_section, "METHOD", i_val=method_id)
151 6 : CALL section_vals_val_get(print_section, "N_FROZEN_CORE", i_val=frozen)
152 6 : CALL section_vals_val_get(print_section, "MAX_ITER", i_val=max_iter)
153 6 : CALL section_vals_val_get(print_section, "ENERGY_CONVERGENCE", r_val=energy_conv)
154 6 : CALL section_vals_val_get(print_section, "AMPLITUDES_CONVERGENCE", r_val=amp_conv)
155 6 : CALL section_vals_val_get(print_section, "MAX_MEMORY", r_val=max_memory_mb)
156 6 : CALL section_vals_val_get(print_section, "OUTPUT_BUFFER_SIZE", r_val=output_buffer_mib)
157 6 : IF (.NOT. ieee_is_finite(output_buffer_mib)) THEN
158 0 : CPABORT("CC4S: OUTPUT_BUFFER_SIZE must be finite")
159 : END IF
160 6 : IF (output_buffer_mib < 1.0_dp/elements_per_mib .OR. &
161 : output_buffer_mib > REAL(HUGE(1), dp)/elements_per_mib) THEN
162 : CALL cp_abort(__LOCATION__, &
163 0 : "CC4S: OUTPUT_BUFFER_SIZE must hold at least one Real64 value and fit the MPI integer count")
164 : END IF
165 6 : buffer_elements = FLOOR(output_buffer_mib*elements_per_mib, KIND=int_8)
166 0 : SELECT CASE (method_id)
167 : CASE (1)
168 0 : method = 'MP2'
169 : CASE (2)
170 0 : method = 'CCSD'
171 : CASE (3)
172 6 : method = 'CCSD_T'
173 : CASE DEFAULT
174 6 : CPABORT("CC4S: unsupported method")
175 : END SELECT
176 6 : IF (.NOT. ieee_is_finite(max_memory_mb)) CPABORT("CC4S: MAX_MEMORY must be finite")
177 6 : IF (max_iter < 1 .OR. energy_conv <= 0.0_dp .OR. amp_conv <= 0.0_dp .OR. max_memory_mb <= 0.0_dp) THEN
178 : CALL cp_abort(__LOCATION__, &
179 0 : "CC4S: iteration, convergence and memory limits must be positive")
180 : END IF
181 : CALL get_qs_env(qs_env, x_data=x_data, mos=mos, para_env=para_env, dft_control=dft_control, &
182 : cell=cell, do_kpoints=do_kpoints, scf_control=scf_control, &
183 6 : matrix_ks=matrix_ks, matrix_s=matrix_s)
184 6 : IF (dft_control%qs_control%cdft) THEN
185 0 : CPABORT("CC4S: constrained DFT is unsupported")
186 : END IF
187 6 : IF (.NOT. converged) THEN
188 0 : CPABORT("CC4S: export requires a converged SCF")
189 : END IF
190 6 : CALL section_vals_val_get(input, "DFT%POISSON%PERIODIC", i_val=poisson_periodic)
191 6 : IF (do_kpoints) THEN
192 0 : CPABORT("CC4S: k-point sampling is unsupported")
193 : END IF
194 16 : periodic = ALL(cell%perd == 1) .AND. poisson_periodic == use_perd_xyz
195 18 : IF (.NOT. periodic .AND. &
196 : (ANY(cell%perd /= 0) .OR. poisson_periodic /= use_perd_none)) THEN
197 : CALL cp_abort(__LOCATION__, &
198 0 : "CC4S: CELL and POISSON must both use PERIODIC NONE or both PERIODIC XYZ")
199 : END IF
200 6 : IF (periodic .AND. dft_control%charge /= 0) THEN
201 0 : CPABORT("CC4S: periodic cells must be neutral")
202 : END IF
203 6 : IF (dft_control%nspins /= 1 .OR. dft_control%uks .OR. dft_control%roks .OR. dft_control%smear) THEN
204 : CALL cp_abort(__LOCATION__, &
205 0 : "CC4S: a restricted closed-shell reference without smearing is required")
206 : END IF
207 6 : admm_section => section_vals_get_subs_vals(input, "DFT%AUXILIARY_DENSITY_MATRIX_METHOD")
208 6 : CALL section_vals_get(admm_section, explicit=admm_section_present)
209 6 : IF (admm_section_present) THEN
210 0 : CPABORT("CC4S: AUXILIARY_DENSITY_MATRIX_METHOD is unsupported")
211 : END IF
212 6 : IF (dft_control%dft_plus_u .OR. dft_control%sic_method_id /= sic_none) THEN
213 0 : CPABORT("CC4S: DFT+U and self-interaction corrections are unsupported")
214 : END IF
215 6 : CALL section_vals_val_get(input, "DFT%XC%XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=functional)
216 6 : IF (functional /= xc_none) THEN
217 : CALL cp_abort(__LOCATION__, &
218 0 : "CC4S: canonical HF orbitals required; use XC_FUNCTIONAL NONE")
219 : END IF
220 6 : IF (.NOT. ASSOCIATED(x_data)) THEN
221 0 : CPABORT("CC4S: enable HF with RI_FLAVOR MO")
222 : END IF
223 6 : IF (SIZE(x_data, 1) /= 1) THEN
224 : CALL cp_abort(__LOCATION__, &
225 0 : "CC4S: exactly one HF section is allowed")
226 : END IF
227 6 : IF (ABS(x_data(1, 1)%general_parameter%fraction - 1.0_dp) > hf_fraction_tolerance) THEN
228 0 : CPABORT("CC4S: HF FRACTION must be 1")
229 : END IF
230 6 : IF (.NOT. x_data(1, 1)%do_hfx_ri) THEN
231 0 : CPABORT("CC4S: enable HF/RI with RI_FLAVOR MO")
232 : END IF
233 6 : IF (.NOT. ASSOCIATED(x_data(1, 1)%ri_data)) THEN
234 0 : CPABORT("CC4S: RI-HFX data are unavailable")
235 : END IF
236 6 : IF (x_data(1, 1)%ri_data%flavor /= ri_mo) THEN
237 0 : CPABORT("CC4S: HF/RI requires RI_FLAVOR MO")
238 : END IF
239 6 : cutoff_radius = 0.0_dp
240 6 : IF (periodic) THEN
241 2 : IF (.NOT. x_data(1, 1)%ri_data%same_op .OR. &
242 : x_data(1, 1)%ri_data%hfx_pot%potential_type /= do_potential_truncated) THEN
243 : CALL cp_abort(__LOCATION__, &
244 0 : "CC4S: periodic export requires TRUNCATED HF interaction and the same RI metric; use RI_METRIC HFX")
245 : END IF
246 2 : cutoff_radius = x_data(1, 1)%ri_data%hfx_pot%cutoff_radius
247 : max_cutoff = 0.5_dp*MIN(plane_distance(1, 0, 0, cell), &
248 2 : plane_distance(0, 1, 0, cell), plane_distance(0, 0, 1, cell))
249 2 : IF (cutoff_radius <= 0.0_dp .OR. cutoff_radius >= max_cutoff) THEN
250 : CALL cp_abort(__LOCATION__, &
251 0 : "CC4S: positive HF CUTOFF_RADIUS must be smaller than half the minimal cell plane spacing")
252 : END IF
253 : ELSE
254 4 : IF (.NOT. x_data(1, 1)%ri_data%same_op .OR. &
255 : x_data(1, 1)%ri_data%hfx_pot%potential_type /= do_potential_coulomb) THEN
256 : CALL cp_abort(__LOCATION__, &
257 0 : "CC4S: molecular HF interaction and RI metric must both be COULOMB")
258 : END IF
259 : END IF
260 :
261 6 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff, homo=homo, nmo=nstored, nao=nao)
262 6 : IF (homo < 1 .OR. homo > nstored) THEN
263 0 : CPABORT("CC4S: no occupied orbitals")
264 : END IF
265 6 : IF (.NOT. has_uniform_occupation(mos(1), first_mo=1, last_mo=homo, occupation=2.0_dp) .OR. &
266 : .NOT. has_uniform_occupation(mos(1), first_mo=homo + 1, last_mo=nstored, occupation=0.0_dp)) THEN
267 : CALL cp_abort(__LOCATION__, &
268 0 : "CC4S: occupations must be 2 below the HOMO and 0 above")
269 : END IF
270 6 : IF (frozen < 0 .OR. frozen >= homo) THEN
271 : CALL cp_abort(__LOCATION__, &
272 0 : "CC4S: N_FROZEN_CORE must leave at least one occupied orbital")
273 : END IF
274 6 : CALL dbt_get_info(x_data(1, 1)%ri_data%t_3c_int_ctr_1(1, 1), nfull_total=dims3)
275 6 : nri = dims3(2)
276 6 : CPASSERT(dims3(1) == nao .AND. dims3(3) == nao .AND. nri > 0)
277 : ! Check quadratic storage before canonicalization; the full AO vertex is never gathered.
278 18 : IF (limited_product([INT(nao, int_8), INT(nao, int_8)], backend_data_limit) > backend_data_limit) THEN
279 0 : CPABORT("CC4S: canonical matrices exceed backend count limits")
280 : END IF
281 18 : IF (limited_product([INT(nri, int_8), INT(nri, int_8)], backend_data_limit) > backend_data_limit) THEN
282 0 : CPABORT("CC4S: RI metric exceeds backend count limits")
283 : END IF
284 : ! Canonicalization precedes batch selection; reject an insufficient budget before allocation.
285 6 : memory_mb = REAL(dp_size, dp)*6.0_dp*REAL(nao, dp)**2/1.0e6_dp
286 6 : IF (memory_mb > max_memory_mb) THEN
287 0 : CPABORT("CC4S: canonicalization exceeds MAX_MEMORY [MB]")
288 : END IF
289 6 : elapsed = 0.0_dp
290 :
291 : ! Always diagonalize private copies, including for OT or a complete stored MO set.
292 : ! Drop overlap null modes using the same cutoff as SCF; do not export placeholder energies.
293 6 : NULLIFY (fm_struct)
294 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, &
295 6 : template_fmstruct=mo_coeff%matrix_struct)
296 6 : CALL cp_fm_create(fm_fock, fm_struct)
297 6 : CALL cp_fm_create(fm_overlap, fm_struct)
298 6 : CALL cp_fm_create(fm_coeff, fm_struct)
299 6 : CALL cp_fm_create(fm_work, fm_struct)
300 6 : CALL cp_fm_struct_release(fm_struct)
301 6 : CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, fm_fock)
302 6 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, fm_overlap)
303 18 : ALLOCATE (eigenvalues(nao))
304 : CALL cp_fm_geeig_canon(fm_fock, fm_overlap, fm_coeff, eigenvalues, fm_work, &
305 6 : scf_control%eps_eigval, nmo_retained=nmo)
306 6 : IF (nmo <= homo) THEN
307 : CALL cp_abort(__LOCATION__, &
308 0 : "CC4S: no virtual orbitals after removing overlap null modes")
309 : END IF
310 6 : IF (eigenvalues(homo + 1) <= eigenvalues(homo)) THEN
311 : CALL cp_abort(__LOCATION__, &
312 0 : "CC4S: a positive HOMO-LUMO gap is required")
313 : END IF
314 6 : fermi_energy = 0.5_dp*(eigenvalues(homo) + eigenvalues(homo + 1))
315 6 : nmo = nmo - frozen
316 6 : total_elements = vertex_size(nri, nmo)
317 6 : buffer_elements = MIN(buffer_elements, total_elements, output_count_limit)
318 : ! The initial coefficient copy can have a different block partition. Bound its
319 : ! metadata by the number of elements before entering DBCSR/DBT conversion.
320 18 : IF (limited_product([INT(nao, int_8), INT(nmo, int_8)], backend_block_limit) > backend_block_limit) THEN
321 0 : CPABORT("CC4S: coefficients exceed backend count limits")
322 : END IF
323 6 : IF (output_unit > 0) THEN
324 3 : WRITE (output_unit, '(/,T3,A)') "CC4S| Reference export"
325 3 : WRITE (output_unit, '(T3,A,T60,A21)') "CC4S| Reference:", "Restricted RI-HF"
326 3 : IF (periodic) THEN
327 1 : WRITE (output_unit, '(T3,A,T60,A21)') "CC4S| Boundary conditions:", "Periodic Gamma"
328 1 : WRITE (output_unit, '(T3,A,T60,A21)') "CC4S| Coulomb interaction:", "TRUNCATED"
329 1 : WRITE (output_unit, '(T3,A,T69,F12.6)') "CC4S| Coulomb cutoff radius [bohr]:", cutoff_radius
330 : ELSE
331 2 : WRITE (output_unit, '(T3,A,T60,A21)') "CC4S| Boundary conditions:", "Isolated"
332 2 : WRITE (output_unit, '(T3,A,T60,A21)') "CC4S| Coulomb interaction:", "COULOMB"
333 : END IF
334 3 : WRITE (output_unit, '(T3,A)') "CC4S|"
335 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Number of auxiliary basis functions:", nri
336 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Number of exported orbitals:", nmo
337 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Number of exported occupied orbitals:", homo - frozen
338 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Number of exported virtual orbitals:", nmo - homo + frozen
339 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Number of frozen occupied orbitals:", frozen
340 : WRITE (output_unit, '(T3,A,T63,F18.6)') &
341 3 : "CC4S| Expected Coulomb vertex size [MiB]:", REAL(total_elements, dp)/elements_per_mib
342 3 : WRITE (output_unit, '(T3,A)') "CC4S|"
343 : WRITE (output_unit, '(T3,A,T69,F12.6)') &
344 3 : "CC4S| Actual output buffer size per MPI rank [MiB]:", REAL(buffer_elements, dp)/elements_per_mib
345 : END IF
346 6 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nmo, template_fmstruct=mo_coeff%matrix_struct)
347 6 : CALL cp_fm_create(fm_active, fm_struct, name="CC4S active orbitals")
348 6 : CALL cp_fm_struct_release(fm_struct)
349 6 : CALL cp_fm_to_fm(fm_coeff, fm_active, ncol=nmo, source_start=frozen + 1)
350 6 : CALL cp_fm_to_dbcsr_row_template(coefficient_matrix, fm_active, matrix_s(1)%matrix)
351 6 : CALL dbt_create(coefficient_matrix, coefficients, name="CC4S coefficients")
352 6 : CALL dbt_copy_matrix_to_tensor(coefficient_matrix, coefficients)
353 6 : CALL dbcsr_release(coefficient_matrix)
354 6 : CALL cp_fm_release(fm_active)
355 6 : CALL cp_fm_release(fm_fock)
356 6 : CALL cp_fm_release(fm_overlap)
357 6 : CALL cp_fm_release(fm_coeff)
358 6 : CALL cp_fm_release(fm_work)
359 :
360 6 : CALL dbt_get_info(x_data(1, 1)%ri_data%t_2c_int(1, 1), nfull_total=dims2)
361 18 : CPASSERT(ALL(dims2 == nri))
362 : ! Transform and export one orbital range at a time; never retain the complete vertex.
363 : CALL export_vertex_batches(x_data(1, 1)%ri_data%t_3c_int_ctr_1(1, 1), &
364 : x_data(1, 1)%ri_data%t_2c_int(1, 1), coefficients, para_env, &
365 : logger, print_section, nmo, max_memory_mb, buffer_elements, &
366 : periodic, cutoff_radius, cell%hmat, output_unit, vertex_norm, elapsed, &
367 6 : homo - frozen, eigenvalues(frozen + 1:frozen + nmo), fermi_energy)
368 6 : CALL dbt_destroy(coefficients)
369 6 : CALL timeset("cc4s_write_files", stage_handle)
370 6 : write_start = m_walltime()
371 6 : IF (para_env%is_source()) THEN
372 3 : CALL write_cc4s_in(logger, print_section, method, homo - frozen, max_iter, energy_conv, amp_conv)
373 : END IF
374 6 : elapsed(5) = elapsed(5) + m_walltime() - write_start
375 6 : CALL timestop(stage_handle)
376 6 : IF (para_env%is_source()) THEN
377 3 : IF (output_unit > 0) THEN
378 3 : WRITE (output_unit, '(/,T3,A)') "CC4S| Export completed"
379 3 : WRITE (output_unit, '(T3,A,T56,ES25.16)') "CC4S| Vertex norm squared [Ha]:", vertex_norm
380 : WRITE (output_unit, '(T3,A,T56,ES25.16)') &
381 56 : "CC4S| Sum of exported orbital energies [Ha]:", SUM(eigenvalues(frozen + 1:frozen + nmo))
382 3 : WRITE (output_unit, '(T3,A)') "CC4S|"
383 3 : filename = cp_print_key_generate_filename(logger, print_section, 'EigenEnergies', '.yaml', .FALSE.)
384 3 : WRITE (output_unit, '(T3,A,T41,A)') "CC4S| Orbital energies metadata:", &
385 6 : ADJUSTR(filename(:MAX(40, LEN_TRIM(filename))))
386 3 : filename = cp_print_key_generate_filename(logger, print_section, 'EigenEnergies', '.elements', .FALSE.)
387 3 : WRITE (output_unit, '(T3,A,T41,A)') "CC4S| Orbital energies data:", &
388 6 : ADJUSTR(filename(:MAX(40, LEN_TRIM(filename))))
389 3 : filename = cp_print_key_generate_filename(logger, print_section, 'CoulombVertex', '.yaml', .FALSE.)
390 3 : WRITE (output_unit, '(T3,A,T41,A)') "CC4S| Coulomb vertex metadata:", &
391 6 : ADJUSTR(filename(:MAX(40, LEN_TRIM(filename))))
392 3 : filename = cp_print_key_generate_filename(logger, print_section, 'CoulombVertex', '.elements', .FALSE.)
393 3 : WRITE (output_unit, '(T3,A,T41,A)') "CC4S| Coulomb vertex file:", &
394 6 : ADJUSTR(filename(:MAX(40, LEN_TRIM(filename))))
395 3 : CALL print_vertex_file_size(output_unit, filename)
396 : END IF
397 : END IF
398 6 : DEALLOCATE (eigenvalues)
399 6 : CALL para_env%sync()
400 6 : elapsed(1) = m_walltime() - export_start
401 : ! Report the longest rank-local elapsed times; only the source rank writes files.
402 6 : CALL para_env%max(elapsed)
403 6 : IF (output_unit > 0) THEN
404 3 : WRITE (output_unit, '(T3,A,T69,F12.3)') "CC4S| Total export time [s]:", elapsed(1)
405 3 : WRITE (output_unit, '(T3,A,T69,F12.3)') "CC4S| Orbital transformation time [s]:", elapsed(2)
406 3 : WRITE (output_unit, '(T3,A,T69,F12.3)') "CC4S| Metric contraction time [s]:", elapsed(3)
407 3 : WRITE (output_unit, '(T3,A,T69,F12.3)') "CC4S| Output gathering time [s]:", elapsed(4)
408 3 : WRITE (output_unit, '(T3,A,T69,F12.3)') "CC4S| File writing time [s]:", elapsed(5)
409 : END IF
410 6 : IF (para_env%is_source()) THEN
411 3 : IF (output_unit > 0) THEN
412 3 : WRITE (output_unit, '(/,T3,A)') "CC4S| Generated CC4S input"
413 3 : filename = cp_print_key_generate_filename(logger, print_section, extension='.in', my_local=.FALSE.)
414 3 : WRITE (output_unit, '(T3,A,T41,A)') "CC4S| Generated input file:", &
415 6 : ADJUSTR(filename(:MAX(40, LEN_TRIM(filename))))
416 3 : IF (method == 'CCSD_T') THEN
417 3 : WRITE (output_unit, '(T3,A,T73,A8)') "CC4S| Requested method:", 'CCSD(T)'
418 : ELSE
419 0 : WRITE (output_unit, '(T3,A,T73,A8)') "CC4S| Requested method:", TRIM(method)
420 : END IF
421 3 : IF (method /= 'MP2') THEN
422 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Maximum CCSD iterations:", max_iter
423 3 : WRITE (output_unit, '(T3,A,T73,ES8.1)') "CC4S| CCSD energy convergence [Ha]:", energy_conv
424 3 : WRITE (output_unit, '(T3,A,T73,ES8.1)') "CC4S| CCSD amplitudes convergence:", amp_conv
425 : END IF
426 : END IF
427 : END IF
428 6 : CALL timestop(handle)
429 74187 : END SUBROUTINE write_cc4s
430 :
431 : ! **************************************************************************************************
432 : !> \brief Report the actual size of the closed binary vertex file, in MiB, on the writing rank.
433 : !> \param output_unit Log output unit
434 : !> \param filename Export filename generated by the print key
435 : ! **************************************************************************************************
436 3 : SUBROUTINE print_vertex_file_size(output_unit, filename)
437 : INTEGER, INTENT(IN) :: output_unit
438 : CHARACTER(len=*), INTENT(IN) :: filename
439 :
440 : INTEGER :: io_status
441 : INTEGER(int_8) :: file_size
442 : REAL(dp) :: size_mib
443 :
444 3 : file_size = -1_int_8
445 3 : size_mib = -1.0_dp
446 3 : INQUIRE (FILE=TRIM(filename), SIZE=file_size, IOSTAT=io_status)
447 3 : IF (io_status == 0) THEN
448 3 : IF (file_size >= 0_int_8) THEN
449 3 : size_mib = REAL(file_size, dp)*REAL(file_storage_size, dp)/(8.0_dp*1024.0_dp**2)
450 : END IF
451 : END IF
452 3 : IF (size_mib >= 0.0_dp) THEN
453 3 : WRITE (output_unit, '(T3,A,T63,F18.6)') "CC4S| Written Coulomb vertex size [MiB]:", size_mib
454 : ELSE
455 0 : WRITE (output_unit, '(T3,A,T63,A18)') "CC4S| Written Coulomb vertex size [MiB]:", "unavailable"
456 : END IF
457 3 : END SUBROUTINE print_vertex_file_size
458 :
459 : ! **************************************************************************************************
460 : !> \brief Transform and stream orbital batches, using private auxiliary slices of the SCF triangle.
461 : !> \param integrals Live SCF tensor, ordered (AO, RI, AO)
462 : !> \param metric Live SCF inverse square root of the RI metric
463 : !> \param coefficients Distributed active canonical orbitals
464 : !> \param para_env MPI parallel environment
465 : !> \param logger Logger for output selection and naming
466 : !> \param print_section CC4S print key
467 : !> \param nmo Number of active orbitals
468 : !> \param max_memory_mb Estimated additional workspace limit per rank, in decimal MB
469 : !> \param buffer_elements Maximum number of Real64 values in each output piece
470 : !> \param periodic Whether the cell is three-dimensionally periodic at Gamma
471 : !> \param cutoff_radius Truncated Coulomb cutoff in bohr
472 : !> \param hmat Cell vectors in bohr
473 : !> \param output_unit Log output unit
474 : !> \param vertex_norm Squared vertex norm on the writing rank
475 : !> \param elapsed Accumulated wall times: total, transformation, metric, gathering, writing
476 : !> \param homo Number of active occupied orbitals
477 : !> \param eigenvalues Active canonical orbital energies
478 : !> \param fermi_energy Energy inside the HOMO-LUMO gap
479 : ! **************************************************************************************************
480 6 : SUBROUTINE export_vertex_batches(integrals, metric, coefficients, para_env, logger, print_section, &
481 : nmo, max_memory_mb, buffer_elements, periodic, cutoff_radius, hmat, &
482 6 : output_unit, vertex_norm, elapsed, homo, eigenvalues, fermi_energy)
483 : TYPE(dbt_type), INTENT(INOUT) :: integrals, metric, coefficients
484 : TYPE(mp_para_env_type), POINTER :: para_env
485 : TYPE(cp_logger_type), POINTER :: logger
486 : TYPE(section_vals_type), POINTER :: print_section
487 : INTEGER, INTENT(IN) :: nmo
488 : REAL(dp), INTENT(IN) :: max_memory_mb
489 : INTEGER(KIND=int_8), INTENT(IN) :: buffer_elements
490 : LOGICAL, INTENT(IN) :: periodic
491 : REAL(dp), INTENT(IN) :: cutoff_radius, hmat(3, 3)
492 : INTEGER, INTENT(IN) :: output_unit
493 : REAL(dp), INTENT(OUT) :: vertex_norm
494 : REAL(dp), INTENT(INOUT) :: elapsed(5)
495 : INTEGER, INTENT(IN) :: homo
496 : REAL(dp), INTENT(IN) :: eigenvalues(nmo), fermi_energy
497 :
498 : INTEGER :: block_width, bounds2(2, 2), bounds3(2, 3), handle, im, ir, iunit, mb, &
499 : metric_blocks(2), mfirst, mlast, nao, nri, pdims(3), pdims2(2), rb, rfirst, rlast
500 6 : INTEGER, ALLOCATABLE :: ao(:), d1(:), d2(:), d3(:), mo(:), ri(:)
501 : REAL(dp) :: memory_mb, pressure, start, trial_m, &
502 : trial_r
503 6 : REAL(dp), ALLOCATABLE :: values(:)
504 30 : TYPE(dbt_pgrid_type) :: grid, grid2
505 150 : TYPE(dbt_type) :: ao_triangle, coeff, coeff_batch, gamma, &
506 114 : half, metric_copy, triangle
507 :
508 30 : ALLOCATE (ao(dbt_nblks_total(integrals, 1)), ri(dbt_nblks_total(integrals, 2)))
509 6 : CALL dbt_get_info(integrals, blk_size_1=ao, blk_size_2=ri)
510 28 : nao = SUM(ao)
511 96 : nri = SUM(ri)
512 : ! Smaller orbital blocks may be necessary for one output stripe to fit the backend.
513 6 : block_width = 16
514 : DO WHILE (limited_product([INT(nri, int_8), INT(nmo, int_8), INT(block_width, int_8)], &
515 24 : backend_data_limit) > backend_data_limit .AND. block_width > 1)
516 0 : block_width = block_width/2
517 : END DO
518 18 : ALLOCATE (mo(1 + (nmo - 1)/block_width))
519 16 : mo = block_width
520 6 : mo(SIZE(mo)) = nmo - (SIZE(mo) - 1)*block_width
521 6 : CALL dbt_get_info(metric, nblks_total=metric_blocks)
522 : ! Reblocking can refine both original and requested auxiliary partitions.
523 24 : IF (limited_product([INT(para_env%num_pe, int_8), INT(metric_blocks(1), int_8) + SIZE(ri, KIND=int_8), &
524 : INT(metric_blocks(2), int_8) + SIZE(ri, KIND=int_8)], &
525 6 : backend_block_limit) > backend_block_limit) THEN
526 0 : CPABORT("CC4S: RI metric metadata exceeds backend limits")
527 : END IF
528 6 : mb = SIZE(mo)
529 6 : rb = SIZE(ri)
530 : ! Bound the busiest rank without assuming balanced distribution. Include private operands,
531 : ! simultaneously live results, and conservative redistribution/communication copies.
532 6 : pressure = batch_pressure(ao, ri, mo, mb, rb, para_env%num_pe, buffer_elements, max_memory_mb)
533 12 : DO WHILE (pressure > 1.0_dp .AND. (mb > 1 .OR. rb > 1))
534 6 : trial_m = HUGE(0.0_dp)
535 6 : trial_r = HUGE(0.0_dp)
536 6 : IF (mb > 1) trial_m = batch_pressure(ao, ri, mo, mb/2 + MODULO(mb, 2), rb, para_env%num_pe, &
537 4 : buffer_elements, max_memory_mb)
538 6 : IF (rb > 1) trial_r = batch_pressure(ao, ri, mo, mb, rb/2 + MODULO(rb, 2), para_env%num_pe, &
539 6 : buffer_elements, max_memory_mb)
540 6 : IF (trial_m <= trial_r) THEN
541 2 : mb = mb/2 + MODULO(mb, 2)
542 : ELSE
543 4 : rb = rb/2 + MODULO(rb, 2)
544 : END IF
545 12 : pressure = MIN(trial_m, trial_r)
546 : END DO
547 6 : memory_mb = batch_workspace(nao, nri, nmo, mo, ri, mb, rb, buffer_elements)
548 6 : IF (.NOT. batch_counts_fit(ao, ri, mo, mb, rb, para_env%num_pe)) THEN
549 0 : CPABORT("CC4S: minimum batch exceeds backend count limits")
550 : END IF
551 6 : IF (output_unit > 0) THEN
552 : WRITE (output_unit, '(T3,A,T69,F12.3)') &
553 3 : "CC4S| Estimated export workspace per MPI rank [MB]:", memory_mb
554 : WRITE (output_unit, '(T3,A,T69,F12.3)') &
555 3 : "CC4S| Export workspace budget per MPI rank [MB]:", max_memory_mb
556 : END IF
557 6 : IF (memory_mb > max_memory_mb) THEN
558 : CALL cp_abort(__LOCATION__, &
559 0 : "CC4S: minimum transformation batch exceeds MAX_MEMORY; increase the limit or reduce OUTPUT_BUFFER_SIZE")
560 : END IF
561 6 : IF (output_unit > 0) THEN
562 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Maximum orbitals per transformation batch:", batch_width(mo, mb)
563 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Maximum auxiliary functions per batch:", batch_width(ri, rb)
564 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Number of orbital batches:", 1 + (SIZE(mo) - 1)/mb
565 3 : WRITE (output_unit, '(T3,A,T69,I12)') "CC4S| Number of auxiliary batches:", 1 + (SIZE(ri) - 1)/rb
566 : END IF
567 6 : pdims = 0
568 6 : pdims2 = 0
569 24 : CALL dbt_pgrid_create(para_env, pdims, grid, tensor_dims=[SIZE(mo), SIZE(ri), SIZE(ao)])
570 18 : CALL dbt_pgrid_create(para_env, pdims2, grid2, tensor_dims=[SIZE(ao), SIZE(mo)])
571 6 : CALL create_2c_tensor(coeff, d1, d2, grid2, ao, mo, name="CC4S active coefficients")
572 6 : DEALLOCATE (d1, d2)
573 6 : CALL dbt_copy(coefficients, coeff)
574 : ! The metric and integrals can have different auxiliary block partitions.
575 6 : CALL create_2c_tensor(metric_copy, d1, d2, grid2, ri, ri, name="CC4S inverse metric")
576 6 : DEALLOCATE (d1, d2)
577 6 : CALL dbt_copy(metric, metric_copy)
578 18 : ALLOCATE (values(buffer_elements))
579 6 : vertex_norm = 0.0_dp
580 6 : iunit = -1
581 6 : start = m_walltime()
582 6 : IF (para_env%is_source()) THEN
583 3 : CALL write_eigen_energies_files(logger, print_section, nmo, homo, eigenvalues, fermi_energy)
584 3 : CALL open_coulomb_vertex_file(logger, print_section, nri, nmo, periodic, cutoff_radius, hmat, iunit)
585 : END IF
586 6 : elapsed(5) = elapsed(5) + m_walltime() - start
587 14 : DO im = 1, SIZE(mo), mb
588 8 : CALL dbt_create(coeff, coeff_batch, name="CC4S coefficient batch")
589 : CALL create_3c_tensor(gamma, d1, d2, d3, grid, ri, mo, mo, [1, 2], [3], &
590 8 : name="CC4S Coulomb vertex batch")
591 8 : DEALLOCATE (d1, d2, d3)
592 10 : mfirst = SUM(mo(:im - 1)) + 1
593 20 : mlast = SUM(mo(:im + MIN(mb - 1, SIZE(mo) - im)))
594 24 : bounds2(:, 1) = [1, nao]
595 24 : bounds2(:, 2) = [mfirst, mlast]
596 8 : CALL dbt_copy(coeff, coeff_batch, bounds=bounds2)
597 8 : DO ir = 1, SIZE(ri), rb
598 164 : rfirst = SUM(ri(:ir - 1)) + 1
599 296 : rlast = SUM(ri(:ir + MIN(rb - 1, SIZE(ri) - ir)))
600 20 : start = m_walltime()
601 20 : CALL timeset("cc4s_transform", handle)
602 60 : bounds3(:, 1) = [1, nao]
603 60 : bounds3(:, 2) = [rfirst, rlast]
604 60 : bounds3(:, 3) = [1, nao]
605 20 : CALL dbt_create(integrals, ao_triangle)
606 : CALL create_3c_tensor(half, d1, d2, d3, grid, mo, ri, ao, [1], [2, 3], &
607 20 : name="CC4S half transform")
608 20 : DEALLOCATE (d1, d2, d3)
609 : ! Copy only this auxiliary slice. Contractions never receive live SCF tensors.
610 20 : CALL dbt_copy(integrals, ao_triangle, bounds=bounds3)
611 : CALL dbt_contract(1.0_dp, coeff_batch, ao_triangle, 0.0_dp, half, &
612 : contract_1=[1], notcontract_1=[2], contract_2=[1], notcontract_2=[2, 3], &
613 20 : map_1=[1], map_2=[2, 3], filter_eps=0.0_dp)
614 : ! Restore the weighted AO triangle before transforming the second orbital index.
615 : ! Sum, not average: diagonal atom blocks carry weight 1/2.
616 : CALL dbt_contract(1.0_dp, coeff_batch, ao_triangle, 1.0_dp, half, &
617 : contract_1=[1], notcontract_1=[2], contract_2=[3], notcontract_2=[2, 1], &
618 20 : map_1=[1], map_2=[2, 3], filter_eps=0.0_dp)
619 20 : CALL dbt_destroy(ao_triangle)
620 : CALL create_3c_tensor(triangle, d1, d2, d3, grid, ri, mo, mo, [1, 2], [3], &
621 20 : name="CC4S MO batch")
622 20 : DEALLOCATE (d1, d2, d3)
623 : CALL dbt_contract(1.0_dp, half, coeff, 0.0_dp, triangle, &
624 : contract_1=[3], notcontract_1=[1, 2], contract_2=[1], notcontract_2=[2], &
625 20 : map_1=[3, 1], map_2=[2], filter_eps=0.0_dp)
626 20 : CALL dbt_destroy(half)
627 20 : CALL timestop(handle)
628 20 : elapsed(2) = elapsed(2) + m_walltime() - start
629 20 : start = m_walltime()
630 20 : CALL timeset("cc4s_metric", handle)
631 : CALL dbt_contract(1.0_dp, metric_copy, triangle, 1.0_dp, gamma, &
632 : contract_1=[1], notcontract_1=[2], contract_2=[1], notcontract_2=[2, 3], &
633 20 : map_1=[1], map_2=[2, 3], filter_eps=0.0_dp)
634 20 : CALL dbt_destroy(triangle)
635 20 : CALL timestop(handle)
636 60 : elapsed(3) = elapsed(3) + m_walltime() - start
637 : END DO
638 8 : CALL write_coulomb_vertex_batch(para_env, nri, nmo, gamma, iunit, mfirst, mlast, values, vertex_norm, elapsed)
639 8 : CALL dbt_destroy(gamma)
640 14 : CALL dbt_destroy(coeff_batch)
641 : END DO
642 6 : start = m_walltime()
643 6 : IF (para_env%is_source()) CALL cp_print_key_finished_output(iunit, logger, print_section)
644 6 : elapsed(5) = elapsed(5) + m_walltime() - start
645 6 : DEALLOCATE (values)
646 6 : CALL dbt_destroy(metric_copy)
647 6 : CALL dbt_destroy(coeff)
648 6 : CALL dbt_pgrid_destroy(grid)
649 6 : CALL dbt_pgrid_destroy(grid2)
650 6 : DEALLOCATE (ao, ri, mo)
651 6 : END SUBROUTINE export_vertex_batches
652 :
653 : ! **************************************************************************************************
654 : !> \brief Positive product, saturated at limit+1, without overflowing intermediate products.
655 : !> \param factors Nonnegative factors
656 : !> \param limit Maximum supported value, strictly below HUGE(int_8)
657 : !> \return Product or limit+1 if the limit would be exceeded
658 : ! **************************************************************************************************
659 234 : PURE FUNCTION limited_product(factors, limit) RESULT(value)
660 : INTEGER(int_8), INTENT(IN) :: factors(:), limit
661 : INTEGER(int_8) :: value
662 :
663 : INTEGER :: i
664 :
665 234 : value = 1_int_8
666 1116 : DO i = 1, SIZE(factors)
667 1116 : IF (factors(i) == 0_int_8) THEN
668 234 : value = 0_int_8
669 : RETURN
670 : END IF
671 : END DO
672 1116 : DO i = 1, SIZE(factors)
673 882 : IF (value > limit/factors(i)) THEN
674 0 : value = limit + 1_int_8
675 : RETURN
676 : END IF
677 1116 : value = value*factors(i)
678 : END DO
679 : END FUNCTION limited_product
680 :
681 : ! **************************************************************************************************
682 : !> \brief Validate the complete stream extent before constructing offsets or opening files.
683 : !> \param nri Number of auxiliary functions
684 : !> \param nmo Number of active orbitals
685 : !> \return Number of exported Real64 elements
686 : ! **************************************************************************************************
687 6 : FUNCTION vertex_size(nri, nmo) RESULT(elements)
688 : INTEGER, INTENT(IN) :: nri, nmo
689 : INTEGER(int_8) :: elements
690 :
691 : INTEGER(int_8), PARAMETER :: file_limit = HUGE(0_int_8)/dp_size
692 :
693 6 : CPASSERT(nri > 0 .AND. nmo > 0)
694 24 : elements = limited_product([INT(nri, int_8), INT(nmo, int_8), INT(nmo, int_8)], file_limit)
695 6 : IF (elements > file_limit) CPABORT("CC4S: vertex exceeds stream position range")
696 6 : END FUNCTION vertex_size
697 :
698 : ! **************************************************************************************************
699 : !> \brief Intersect a block segment with an output piece before narrowing to local indices.
700 : !> \param start Zero-based global offset of the block segment
701 : !> \param length Number of values in the block segment
702 : !> \param offset Zero-based global offset of the output piece
703 : !> \param count Number of values in the output piece
704 : !> \param lo First overlapping block index, or 1 for an empty intersection
705 : !> \param hi Last overlapping block index, or 0 for an empty intersection
706 : ! **************************************************************************************************
707 15432183 : PURE SUBROUTINE segment_overlap(start, length, offset, count, lo, hi)
708 : INTEGER(int_8), INTENT(IN) :: start
709 : INTEGER, INTENT(IN) :: length
710 : INTEGER(int_8), INTENT(IN) :: offset, count
711 : INTEGER, INTENT(OUT) :: lo, hi
712 :
713 : INTEGER(int_8) :: first, last
714 :
715 15432183 : first = MAX(start, offset)
716 15432183 : last = MIN(start + INT(length, int_8), offset + count)
717 15432183 : lo = 1
718 15432183 : hi = 0
719 15432183 : IF (first >= last) RETURN
720 14863 : lo = INT(first - start) + 1
721 14863 : hi = INT(last - start)
722 : END SUBROUTINE segment_overlap
723 :
724 : ! **************************************************************************************************
725 : !> \brief Bound backend data and metadata independently of the user memory budget.
726 : !> \param ao AO block sizes
727 : !> \param ri Auxiliary block sizes
728 : !> \param mo MO block sizes
729 : !> \param mb MO blocks in a batch
730 : !> \param rb Auxiliary blocks in a batch
731 : !> \param nranks Maximum replication multiplicity
732 : !> \param pressure Maximum ratio of a count bound to its supported limit
733 : !> \return Whether every intermediate is safely representable, even on one rank
734 : ! **************************************************************************************************
735 22 : FUNCTION batch_counts_fit(ao, ri, mo, mb, rb, nranks, pressure) RESULT(fits)
736 : INTEGER, INTENT(IN) :: ao(:), ri(:), mo(:), mb, rb, nranks
737 : REAL(dp), INTENT(OUT), OPTIONAL :: pressure
738 : LOGICAL :: fits
739 :
740 : INTEGER(int_8), PARAMETER :: product_limit = HUGE(0_int_8)/dp_size
741 :
742 : INTEGER(int_8) :: a, b, block_count, data_count, m, na, &
743 : nb, nm, nq, nr, q, r
744 :
745 116 : a = SUM(INT(ao, int_8))
746 412 : r = SUM(INT(ri, int_8))
747 62 : m = SUM(INT(mo, int_8))
748 22 : b = INT(batch_width(mo, mb), int_8)
749 22 : q = INT(batch_width(ri, rb), int_8)
750 22 : na = SIZE(ao, KIND=int_8)
751 22 : nr = SIZE(ri, KIND=int_8)
752 22 : nm = SIZE(mo, KIND=int_8)
753 22 : nb = INT(mb, int_8)
754 22 : nq = INT(rb, int_8)
755 : data_count = MAX(limited_product([a, a, q], product_limit), &
756 : limited_product([b, q, a], product_limit), &
757 : limited_product([q, m, b], product_limit), &
758 286 : limited_product([r, m, b], product_limit))
759 : ! Allow refinement in three dimensions and replicated partial-result metadata
760 : ! received from every rank. Payload buffers in TAS replication are per peer.
761 : block_count = MAX(limited_product([8_int_8, INT(nranks, int_8), na, nm], product_limit), &
762 : limited_product([8_int_8, INT(nranks, int_8), na, na, nq], product_limit), &
763 : limited_product([8_int_8, INT(nranks, int_8), nb, nq, na], product_limit), &
764 : limited_product([8_int_8, INT(nranks, int_8), nq, nm, nb], product_limit), &
765 550 : limited_product([8_int_8, INT(nranks, int_8), nr, nm, nb], product_limit))
766 22 : fits = data_count <= backend_data_limit .AND. block_count <= backend_block_limit
767 22 : IF (PRESENT(pressure)) pressure = MAX(REAL(data_count, dp)/REAL(backend_data_limit, dp), &
768 16 : REAL(block_count, dp)/REAL(backend_block_limit, dp))
769 22 : END FUNCTION batch_counts_fit
770 :
771 : ! **************************************************************************************************
772 : !> \brief Score a candidate batch against both workspace and integer-count limits.
773 : !> \param ao AO block sizes
774 : !> \param ri Auxiliary block sizes
775 : !> \param mo MO block sizes
776 : !> \param mb MO blocks in a batch
777 : !> \param rb Auxiliary blocks in a batch
778 : !> \param nranks Maximum replication multiplicity
779 : !> \param buffer_elements Output buffer length
780 : !> \param max_memory_mb Additional workspace budget
781 : !> \return Score no larger than one for an admissible batch
782 : ! **************************************************************************************************
783 16 : FUNCTION batch_pressure(ao, ri, mo, mb, rb, nranks, buffer_elements, max_memory_mb) RESULT(pressure)
784 : INTEGER, INTENT(IN) :: ao(:), ri(:), mo(:), mb, rb, nranks
785 : INTEGER(int_8), INTENT(IN) :: buffer_elements
786 : REAL(dp), INTENT(IN) :: max_memory_mb
787 : REAL(dp) :: pressure
788 :
789 : REAL(dp) :: count_pressure
790 :
791 418 : pressure = batch_workspace(SUM(ao), SUM(ri), SUM(mo), mo, ri, mb, rb, buffer_elements)/max_memory_mb
792 16 : IF (.NOT. batch_counts_fit(ao, ri, mo, mb, rb, nranks, count_pressure)) pressure = MAX(pressure, count_pressure)
793 16 : END FUNCTION batch_pressure
794 :
795 : ! **************************************************************************************************
796 : !> \brief Largest number of functions in a group of consecutive tensor blocks.
797 : !> \param sizes Block sizes
798 : !> \param blocks Number of blocks per batch
799 : !> \return Maximum batch width
800 : ! **************************************************************************************************
801 94 : FUNCTION batch_width(sizes, blocks) RESULT(width)
802 : INTEGER, INTENT(IN) :: sizes(:), blocks
803 : INTEGER :: width
804 :
805 : INTEGER :: i
806 :
807 94 : width = 0
808 94 : DO i = 1, SIZE(sizes), blocks
809 1068 : width = MAX(width, SUM(sizes(i:i + MIN(blocks - 1, SIZE(sizes) - i))))
810 : END DO
811 94 : END FUNCTION batch_width
812 :
813 : ! **************************************************************************************************
814 : !> \brief Conservative busiest-rank workspace estimate, without a balanced-distribution assumption.
815 : !> Four copies allow for private operands, contraction results and redistribution buffers.
816 : !> Excludes existing SCF storage, allocator overhead and CC4S. This is not an RSS limit.
817 : !> \param nao Number of AO functions
818 : !> \param nri Number of auxiliary functions
819 : !> \param nmo Number of active orbitals
820 : !> \param mo Orbital block sizes
821 : !> \param ri Auxiliary block sizes
822 : !> \param mb Orbital blocks per batch
823 : !> \param rb Auxiliary blocks per batch
824 : !> \param buffer_elements Reusable output buffer length
825 : !> \return Estimated additional memory per rank in decimal MB
826 : ! **************************************************************************************************
827 22 : FUNCTION batch_workspace(nao, nri, nmo, mo, ri, mb, rb, buffer_elements) RESULT(memory_mb)
828 : INTEGER, INTENT(IN) :: nao, nri, nmo, mo(:), ri(:), mb, rb
829 : INTEGER(KIND=int_8), INTENT(IN) :: buffer_elements
830 : REAL(dp) :: memory_mb
831 :
832 : REAL(dp) :: a, b, fixed, m, q, r, work
833 :
834 22 : a = REAL(nao, dp)
835 22 : r = REAL(nri, dp)
836 22 : m = REAL(nmo, dp)
837 22 : b = REAL(batch_width(mo, mb), dp)
838 22 : q = REAL(batch_width(ri, rb), dp)
839 22 : fixed = 4.0_dp*(r*r + a*m) + 2.0_dp*REAL(buffer_elements, dp)
840 22 : work = 4.0_dp*(a*a*q + b*q*a + q*m*b + r*m*b + a*b)
841 22 : memory_mb = REAL(dp_size, dp)*MAX(6.0_dp*a*a, fixed + work)/1.0e6_dp
842 22 : END FUNCTION batch_workspace
843 :
844 : ! **************************************************************************************************
845 : !> \brief Write EigenEnergies.yaml and EigenEnergies.elements (CC4S format, Hartree units).
846 : !> \param logger Logger for output selection and naming
847 : !> \param print_section CC4S print key
848 : !> \param nMO Number of active canonical orbitals
849 : !> \param homo Number of active occupied orbitals
850 : !> \param eigenvalues Active orbital energies in Hartree
851 : !> \param fermi_energy Energy inside the HOMO-LUMO gap
852 : ! **************************************************************************************************
853 3 : SUBROUTINE write_eigen_energies_files(logger, print_section, nMO, homo, eigenvalues, fermi_energy)
854 : TYPE(cp_logger_type), POINTER :: logger
855 : TYPE(section_vals_type), POINTER :: print_section
856 : INTEGER, INTENT(IN) :: nMO, homo
857 : REAL(KIND=dp), INTENT(IN) :: eigenvalues(nMO), fermi_energy
858 :
859 : INTEGER :: i, iunit
860 :
861 : iunit = cp_print_key_unit_nr(logger, print_section, extension='.yaml', middle_name='EigenEnergies', &
862 3 : file_status='REPLACE', file_position='REWIND')
863 3 : WRITE (iunit, '(A)') 'version: 100'
864 3 : WRITE (iunit, '(A)') 'type: Tensor'
865 3 : WRITE (iunit, '(A)') 'scalarType: Real64'
866 3 : WRITE (iunit, '(A)') 'dimensions:'
867 3 : WRITE (iunit, '(A,I0)') '- length: ', nMO
868 3 : WRITE (iunit, '(A)') ' type: State'
869 3 : WRITE (iunit, '(A)') 'elements:'
870 3 : WRITE (iunit, '(A)') ' type: TextFile'
871 3 : WRITE (iunit, '(A)') 'unit: 1.0'
872 3 : WRITE (iunit, '(A)') 'metaData:'
873 3 : WRITE (iunit, '(A,ES25.16)') ' fermiEnergy: ', fermi_energy
874 3 : WRITE (iunit, '(A,I0)') ' homoIndex: ', homo
875 3 : WRITE (iunit, '(A)') ' energies:'
876 56 : DO i = 1, nMO
877 56 : WRITE (iunit, '(A,ES25.16)') ' - ', eigenvalues(i)
878 : END DO
879 3 : CALL cp_print_key_finished_output(iunit, logger, print_section)
880 :
881 : iunit = cp_print_key_unit_nr(logger, print_section, extension='.elements', middle_name='EigenEnergies', &
882 3 : file_status='REPLACE', file_position='REWIND')
883 56 : DO i = 1, nMO
884 56 : WRITE (iunit, '(ES25.16)') eigenvalues(i)
885 : END DO
886 3 : CALL cp_print_key_finished_output(iunit, logger, print_section)
887 :
888 3 : END SUBROUTINE write_eigen_energies_files
889 :
890 : ! **************************************************************************************************
891 : !> \brief Gather and write one orbital range in bounded pieces, preserving Fortran element order.
892 : !> \param para_env MPI parallel environment
893 : !> \param nRI Number of auxiliary functions
894 : !> \param nMO Number of active canonical orbitals
895 : !> \param gamma Distributed vertex, with data only in the selected third-index range
896 : !> \param iunit Open stream unit on the writing rank
897 : !> \param first_mo First orbital in this batch
898 : !> \param last_mo Last orbital in this batch
899 : !> \param values Reusable output buffer
900 : !> \param vertex_norm Accumulated squared norm on the writing rank
901 : !> \param elapsed Accumulated wall times
902 : ! **************************************************************************************************
903 8 : SUBROUTINE write_coulomb_vertex_batch(para_env, nRI, nMO, gamma, iunit, first_mo, last_mo, &
904 8 : values, vertex_norm, elapsed)
905 : TYPE(mp_para_env_type), POINTER :: para_env
906 : INTEGER, INTENT(IN) :: nRI, nMO
907 : TYPE(dbt_type), INTENT(INOUT) :: gamma
908 : INTEGER, INTENT(IN) :: iunit, first_mo, last_mo
909 : REAL(dp), INTENT(INOUT) :: values(:), vertex_norm, elapsed(5)
910 :
911 : INTEGER :: column, handle, hi, ind(3), j, k, lo, &
912 : sizes(3)
913 : INTEGER(KIND=int_8) :: count, first, last, mirror, nvalues, &
914 : offset, row, segment, start, used, &
915 : within
916 8 : INTEGER, ALLOCATABLE :: off1(:), off2(:), off3(:)
917 : LOGICAL :: found
918 : REAL(dp) :: start_time
919 8 : REAL(dp), ALLOCATABLE :: block(:, :, :)
920 : TYPE(dbt_iterator_type) :: iter
921 :
922 16 : start_time = m_walltime()
923 8 : CALL timeset("cc4s_gather", handle)
924 8 : CPASSERT(SIZE(values, KIND=int_8) <= output_count_limit)
925 0 : ALLOCATE (off1(dbt_nblks_total(gamma, 1)), off2(dbt_nblks_total(gamma, 2)), &
926 56 : off3(dbt_nblks_total(gamma, 3)))
927 8 : CALL dbt_get_info(gamma, blk_offset_1=off1, blk_offset_2=off2, blk_offset_3=off3)
928 114 : DO column = first_mo, last_mo
929 106 : offset = INT(nRI, int_8)*INT(nMO, int_8)*INT(column - 1, int_8)
930 106 : nvalues = offset + INT(nRI, int_8)*INT(column, int_8)
931 : ! Gather one triangle only. Write its transpose from identical values for exact symmetry.
932 5548 : DO WHILE (offset < nvalues)
933 5434 : count = MIN(SIZE(values, KIND=int_8), nvalues - offset)
934 3287490 : values = 0.0_dp
935 : ! Blocks own disjoint elements; threads fill disjoint parts of the bounded buffer.
936 : !$OMP PARALLEL DEFAULT(NONE) &
937 : !$OMP SHARED(gamma, off1, off2, off3, nRI, nMO, offset, count, values) &
938 5434 : !$OMP PRIVATE(iter, ind, sizes, block, found, j, k, first, last, start, lo, hi)
939 : CALL dbt_iterator_start(iter, gamma)
940 : DO WHILE (dbt_iterator_blocks_left(iter))
941 : CALL dbt_iterator_next_block(iter, ind, blk_size=sizes)
942 : first = INT(off1(ind(1)) - 1, int_8) + INT(nRI, int_8)* &
943 : (INT(off2(ind(2)) - 1, int_8) + INT(nMO, int_8)*INT(off3(ind(3)) - 1, int_8))
944 : last = first + INT(sizes(1), int_8) + INT(nRI, int_8)* &
945 : (INT(sizes(2) - 1, int_8) + INT(nMO, int_8)*INT(sizes(3) - 1, int_8))
946 : IF (first >= offset + count .OR. last <= offset) CYCLE
947 : CALL dbt_get_block(gamma, ind, block, found)
948 : CPASSERT(found)
949 : DO k = 1, sizes(3)
950 : DO j = 1, sizes(2)
951 : start = first + INT(nRI, int_8)*(INT(j - 1, int_8) + INT(nMO, int_8)*INT(k - 1, int_8))
952 : CALL segment_overlap(start, sizes(1), offset, count, lo, hi)
953 : IF (lo > hi) CYCLE
954 : values(start + lo - offset:start + hi - offset) = block(lo:hi, j, k)
955 : END DO
956 : END DO
957 : DEALLOCATE (block)
958 : END DO
959 : CALL dbt_iterator_stop(iter)
960 : !$OMP END PARALLEL
961 285034 : CALL para_env%sum(values(:count), para_env%source)
962 5434 : CALL timestop(handle)
963 5434 : elapsed(4) = elapsed(4) + m_walltime() - start_time
964 5434 : start_time = m_walltime()
965 5434 : CALL timeset("cc4s_write_files", handle)
966 5434 : IF (para_env%is_source()) THEN
967 72617 : CALL write_coulomb_vertex_block(iunit, offset, values(:count))
968 72617 : vertex_norm = vertex_norm + SUM(values(:count)**2)
969 2717 : used = 0_int_8
970 5985 : DO WHILE (used < count)
971 3268 : within = MODULO(offset + used, INT(nRI, int_8)*INT(nMO, int_8))
972 3268 : row = within/INT(nRI, int_8) + 1_int_8
973 3268 : within = MODULO(within, INT(nRI, int_8))
974 3268 : segment = MIN(count - used, INT(nRI, int_8) - within)
975 3268 : IF (row < column) THEN
976 3002 : mirror = within + INT(nRI, int_8)*(INT(column - 1, int_8) + INT(nMO, int_8)*(row - 1_int_8))
977 67234 : CALL write_coulomb_vertex_block(iunit, mirror, values(used + 1:used + segment))
978 67234 : vertex_norm = vertex_norm + SUM(values(used + 1:used + segment)**2)
979 : END IF
980 3268 : used = used + segment
981 : END DO
982 : END IF
983 5434 : CALL timestop(handle)
984 5434 : elapsed(5) = elapsed(5) + m_walltime() - start_time
985 5434 : start_time = m_walltime()
986 5434 : CALL timeset("cc4s_gather", handle)
987 5434 : offset = offset + count
988 : END DO
989 : END DO
990 8 : DEALLOCATE (off1, off2, off3)
991 8 : CALL timestop(handle)
992 8 : elapsed(4) = elapsed(4) + m_walltime() - start_time
993 16 : END SUBROUTINE write_coulomb_vertex_batch
994 :
995 : ! **************************************************************************************************
996 : !> \brief Write the Coulomb vertex header and open its raw binary elements file.
997 : !> \param logger Logger for output selection and naming
998 : !> \param print_section CC4S print key
999 : !> \param nRI Number of auxiliary functions
1000 : !> \param nMO Number of active canonical orbitals
1001 : !> \param periodic Whether the cell is three-dimensionally periodic at Gamma
1002 : !> \param cutoff_radius Truncated Coulomb cutoff in bohr (zero for molecular Coulomb)
1003 : !> \param hmat Cell vectors as columns, in bohr
1004 : !> \param iunit Open stream unit for vertex elements
1005 : ! **************************************************************************************************
1006 3 : SUBROUTINE open_coulomb_vertex_file(logger, print_section, nRI, nMO, periodic, cutoff_radius, hmat, iunit)
1007 : TYPE(cp_logger_type), POINTER :: logger
1008 : TYPE(section_vals_type), POINTER :: print_section
1009 : INTEGER, INTENT(IN) :: nRI, nMO
1010 : LOGICAL, INTENT(IN) :: periodic
1011 : REAL(KIND=dp), INTENT(IN) :: cutoff_radius, hmat(3, 3)
1012 : INTEGER, INTENT(OUT) :: iunit
1013 :
1014 : INTEGER :: i
1015 :
1016 : iunit = cp_print_key_unit_nr(logger, print_section, extension='.yaml', middle_name='CoulombVertex', &
1017 3 : file_status='REPLACE', file_position='REWIND')
1018 3 : WRITE (iunit, '(A)') 'version: 100'
1019 3 : WRITE (iunit, '(A)') 'type: Tensor'
1020 3 : WRITE (iunit, '(A)') 'scalarType: Real64'
1021 3 : WRITE (iunit, '(A)') 'dimensions:'
1022 3 : WRITE (iunit, '(A,I0)') '- length: ', nRI
1023 3 : WRITE (iunit, '(A)') ' type: AuxiliaryField'
1024 3 : WRITE (iunit, '(A,I0)') '- length: ', nMO
1025 3 : WRITE (iunit, '(A)') ' type: State'
1026 3 : WRITE (iunit, '(A,I0)') '- length: ', nMO
1027 3 : WRITE (iunit, '(A)') ' type: State'
1028 3 : WRITE (iunit, '(A)') 'elements:'
1029 3 : WRITE (iunit, '(A)') ' type: IeeeBinaryFile'
1030 3 : WRITE (iunit, '(A)') 'unit: 1.0'
1031 3 : WRITE (iunit, '(A)') 'metaData:'
1032 3 : IF (periodic) THEN
1033 1 : WRITE (iunit, '(A)') ' boundaryConditions: "periodic-3d-gamma"'
1034 1 : WRITE (iunit, '(A)') ' interactionPotential: "TRUNCATED"'
1035 1 : WRITE (iunit, '(A,ES25.16)') ' cutoffRadiusBohr: ', cutoff_radius
1036 1 : WRITE (iunit, '(A)') ' madelungCorrection: 0.0'
1037 1 : WRITE (iunit, '(A)') ' latticeVectorsBohr:'
1038 4 : DO i = 1, 3
1039 : WRITE (iunit, '(A,ES25.16,A,ES25.16,A,ES25.16,A)') &
1040 4 : ' - [', hmat(1, i), ', ', hmat(2, i), ', ', hmat(3, i), ']'
1041 : END DO
1042 : ELSE
1043 2 : WRITE (iunit, '(A)') ' boundaryConditions: "molecular"'
1044 2 : WRITE (iunit, '(A)') ' interactionPotential: "COULOMB"'
1045 : END IF
1046 3 : CALL cp_print_key_finished_output(iunit, logger, print_section)
1047 :
1048 : iunit = cp_print_key_unit_nr(logger, print_section, extension='.elements', middle_name='CoulombVertex', &
1049 : file_status='REPLACE', file_position='REWIND', &
1050 3 : file_form='UNFORMATTED', file_access='STREAM')
1051 3 : END SUBROUTINE open_coulomb_vertex_file
1052 :
1053 : ! **************************************************************************************************
1054 : !> \brief Write a contiguous vertex block at an explicit position in the raw elements file.
1055 : !> \param iunit Open stream unit for vertex elements
1056 : !> \param offset Zero-based element offset in auxiliary-orbital-orbital storage order
1057 : !> \param values Contiguous block of vertex elements
1058 : ! **************************************************************************************************
1059 5719 : SUBROUTINE write_coulomb_vertex_block(iunit, offset, values)
1060 : INTEGER, INTENT(IN) :: iunit
1061 : INTEGER(KIND=int_8), INTENT(IN) :: offset
1062 : REAL(KIND=dp), CONTIGUOUS, INTENT(IN) :: values(:)
1063 :
1064 : INTEGER :: io_status
1065 : INTEGER(KIND=int_8) :: position
1066 :
1067 : ! Use CP2K's Real64 byte size, independent of direct-access record-length units.
1068 5719 : CPASSERT(offset >= 0_int_8)
1069 5719 : CPASSERT(SIZE(values, KIND=int_8) > 0_int_8)
1070 5719 : CPASSERT(SIZE(values, KIND=int_8) <= HUGE(position)/INT(dp_size, int_8) - offset)
1071 5719 : position = offset*INT(dp_size, int_8) + 1_int_8
1072 5719 : WRITE (iunit, POS=position, IOSTAT=io_status) values
1073 5719 : IF (io_status /= 0) THEN
1074 0 : CPABORT("CC4S: failed writing CoulombVertex.elements")
1075 : END IF
1076 5719 : END SUBROUTINE write_coulomb_vertex_block
1077 :
1078 : ! **************************************************************************************************
1079 : !> \brief Write the CC4S algorithm sequence and solver settings to cc4s.in.
1080 : !> \param logger Logger for output selection and naming
1081 : !> \param print_section CC4S print key
1082 : !> \param method Requested correlation method: MP2, CCSD or CCSD_T
1083 : !> \param nocc Number of exported occupied spatial orbitals
1084 : !> \param max_iter Maximum coupled-cluster iterations
1085 : !> \param energy_conv CC4S energy convergence threshold
1086 : !> \param amp_conv CC4S amplitude convergence threshold
1087 : ! **************************************************************************************************
1088 3 : SUBROUTINE write_cc4s_in(logger, print_section, method, nocc, max_iter, energy_conv, amp_conv)
1089 : TYPE(cp_logger_type), POINTER :: logger
1090 : TYPE(section_vals_type), POINTER :: print_section
1091 : CHARACTER(LEN=*), INTENT(IN) :: method
1092 : INTEGER, INTENT(IN) :: nocc, max_iter
1093 : REAL(dp), INTENT(IN) :: energy_conv, amp_conv
1094 :
1095 : CHARACTER(len=default_path_length) :: filename
1096 : INTEGER :: iunit
1097 :
1098 : iunit = cp_print_key_unit_nr(logger, print_section, extension='.in', &
1099 3 : file_status='REPLACE', file_position='REWIND')
1100 :
1101 : ! --- common header: read tensors, slice, build integrals ---
1102 3 : filename = cp_print_key_generate_filename(logger, print_section, 'CoulombVertex', '.yaml', .FALSE.)
1103 3 : CALL write_tensor_read(iunit, filename, 'CoulombVertex')
1104 3 : filename = cp_print_key_generate_filename(logger, print_section, 'EigenEnergies', '.yaml', .FALSE.)
1105 3 : CALL write_tensor_read(iunit, filename, 'EigenEnergies')
1106 3 : WRITE (iunit, '(A)') '- name: DefineHolesAndParticles'
1107 3 : WRITE (iunit, '(A)') ' in:'
1108 3 : WRITE (iunit, '(A)') ' eigenEnergies: EigenEnergies'
1109 3 : WRITE (iunit, '(A)') ' out:'
1110 3 : WRITE (iunit, '(A)') ' slicedEigenEnergies: EigenEnergies'
1111 3 : WRITE (iunit, '(A)') '- name: SliceOperator'
1112 3 : WRITE (iunit, '(A)') ' in:'
1113 3 : WRITE (iunit, '(A)') ' slicedEigenEnergies: EigenEnergies'
1114 3 : WRITE (iunit, '(A)') ' operator: CoulombVertex'
1115 3 : WRITE (iunit, '(A)') ' out:'
1116 3 : WRITE (iunit, '(A)') ' slicedOperator: CoulombVertex'
1117 3 : WRITE (iunit, '(A)') '- name: VertexCoulombIntegrals'
1118 3 : WRITE (iunit, '(A)') ' in:'
1119 3 : WRITE (iunit, '(A)') ' slicedCoulombVertex: CoulombVertex'
1120 3 : WRITE (iunit, '(A)') ' out:'
1121 3 : WRITE (iunit, '(A)') ' coulombIntegrals: CoulombIntegrals'
1122 :
1123 : ! --- MP2 (always included) ---
1124 3 : WRITE (iunit, '(A)') '- name: SecondOrderPerturbationTheory'
1125 3 : WRITE (iunit, '(A)') ' in:'
1126 3 : WRITE (iunit, '(A)') ' coulombIntegrals: CoulombIntegrals'
1127 3 : WRITE (iunit, '(A)') ' slicedEigenEnergies: EigenEnergies'
1128 3 : WRITE (iunit, '(A)') ' out:'
1129 3 : WRITE (iunit, '(A)') ' energy: Mp2Energy'
1130 :
1131 : ! --- CoupledCluster (CCSD, CCSD_T) ---
1132 6 : SELECT CASE (TRIM(method))
1133 : CASE ('MP2')
1134 : ! The MP2 driver is complete above.
1135 : CASE ('CCSD', 'CCSD_T')
1136 3 : CALL write_coupled_cluster_step(iunit, 'Ccsd', 'Amplitudes', nocc, max_iter, energy_conv, amp_conv)
1137 : CASE DEFAULT
1138 3 : CPABORT("CC4S: unsupported method")
1139 : END SELECT
1140 :
1141 : ! --- PerturbativeTriples (CCSD_T only) ---
1142 3 : IF (TRIM(method) == 'CCSD_T') THEN
1143 3 : WRITE (iunit, '(A)') '- name: PerturbativeTriples'
1144 3 : WRITE (iunit, '(A)') ' in:'
1145 3 : WRITE (iunit, '(A)') ' slicedEigenEnergies: EigenEnergies'
1146 3 : WRITE (iunit, '(A)') ' amplitudes: Amplitudes'
1147 3 : WRITE (iunit, '(A)') ' coulombIntegrals: CoulombIntegrals'
1148 3 : WRITE (iunit, '(A)') ' out: {}'
1149 : END IF
1150 :
1151 3 : CALL cp_print_key_finished_output(iunit, logger, print_section)
1152 :
1153 3 : END SUBROUTINE write_cc4s_in
1154 :
1155 : ! **************************************************************************************************
1156 : !> \brief Write a single CoupledCluster step to the cc4s.in file.
1157 : !> \param iunit Output file unit
1158 : !> \param cc_method CC4S coupled-cluster method
1159 : !> \param amp_key CC4S output key for amplitudes
1160 : !> \param nocc Number of exported occupied spatial orbitals
1161 : !> \param max_iter Maximum coupled-cluster iterations
1162 : !> \param energy_conv CC4S energy convergence threshold
1163 : !> \param amp_conv CC4S amplitude convergence threshold
1164 : ! **************************************************************************************************
1165 3 : SUBROUTINE write_coupled_cluster_step(iunit, cc_method, amp_key, nocc, max_iter, energy_conv, amp_conv)
1166 : INTEGER, INTENT(IN) :: iunit
1167 : CHARACTER(LEN=*), INTENT(IN) :: cc_method, amp_key
1168 : INTEGER, INTENT(IN) :: nocc, max_iter
1169 : REAL(dp), INTENT(IN) :: energy_conv, amp_conv
1170 :
1171 3 : WRITE (iunit, '(A)') '- name: CoupledCluster'
1172 3 : WRITE (iunit, '(A)') ' in:'
1173 3 : WRITE (iunit, '(A)') ' method: '//TRIM(cc_method)
1174 : ! Match CC4S v1.1 defaults using the occupied space after frozen-core removal.
1175 3 : WRITE (iunit, '(A,I0)') ' integralsSliceSize: ', nocc
1176 3 : WRITE (iunit, '(A)') ' slicedEigenEnergies: EigenEnergies'
1177 3 : WRITE (iunit, '(A)') ' coulombIntegrals: CoulombIntegrals'
1178 3 : WRITE (iunit, '(A)') ' slicedCoulombVertex: CoulombVertex'
1179 3 : WRITE (iunit, '(A,I0)') ' maxIterations: ', max_iter
1180 3 : WRITE (iunit, '(A,ES25.16)') ' energyConvergence: ', energy_conv
1181 3 : WRITE (iunit, '(A,ES25.16)') ' amplitudesConvergence: ', amp_conv
1182 3 : WRITE (iunit, '(A)') ' mixer:'
1183 3 : WRITE (iunit, '(A)') ' type: DiisMixer'
1184 3 : WRITE (iunit, '(A)') ' maxResidua: 4'
1185 3 : WRITE (iunit, '(A)') ' out:'
1186 3 : WRITE (iunit, '(A)') ' amplitudes: '//TRIM(amp_key)
1187 3 : WRITE (iunit, '(A)') ' energy: CoupledClusterEnergy'
1188 :
1189 3 : END SUBROUTINE write_coupled_cluster_step
1190 :
1191 : ! **************************************************************************************************
1192 : !> \brief Write a tensor reference relative to the generated CC4S input's directory.
1193 : !> \param iunit CC4S input unit
1194 : !> \param filename Generated tensor header filename
1195 : !> \param destination CC4S tensor name
1196 : ! **************************************************************************************************
1197 6 : SUBROUTINE write_tensor_read(iunit, filename, destination)
1198 : INTEGER, INTENT(IN) :: iunit
1199 : CHARACTER(len=*), INTENT(IN) :: filename, destination
1200 :
1201 : INTEGER :: i
1202 :
1203 6 : WRITE (iunit, '(A)') '- name: Read'
1204 6 : WRITE (iunit, '(A)') ' in:'
1205 6 : WRITE (iunit, '(A)', ADVANCE='NO') " fileName: '"
1206 : ! A quoted YAML scalar preserves spaces and punctuation; apostrophes are doubled.
1207 224 : DO i = INDEX(filename, '/', BACK=.TRUE.) + 1, LEN_TRIM(filename)
1208 218 : WRITE (iunit, '(A)', ADVANCE='NO') filename(i:i)
1209 224 : IF (filename(i:i) == "'") WRITE (iunit, '(A)', ADVANCE='NO') "'"
1210 : END DO
1211 6 : WRITE (iunit, '(A)') "'"
1212 6 : WRITE (iunit, '(A)') ' out:'
1213 6 : WRITE (iunit, '(A)') ' destination: '//destination
1214 6 : END SUBROUTINE write_tensor_read
1215 :
1216 : END MODULE cc4s_export
|
