submodule(mqc_mbe_fragment_distribution_scheme) mqc_unfragmented_workflow implicit none contains module subroutine unfragmented_calculation(sys_geom, config, result_out, json_data) !! Run unfragmented calculation on the entire system (nlevel=0) !! This is a simple single-process calculation without MPI distribution !! If result_out is present, returns result instead of writing JSON and destroying it !! If json_data is present, populates it for centralized JSON output use mqc_error, only: error_t use mqc_vibrational_analysis, only: compute_vibrational_frequencies, & compute_vibrational_analysis, print_vibrational_analysis use mqc_thermochemistry, only: thermochemistry_result_t, compute_thermochemistry use mqc_json_output_types, only: json_output_data_t, OUTPUT_MODE_UNFRAGMENTED type(system_geometry_t), intent(in) :: sys_geom type(driver_config_t), intent(in) :: config !! Driver configuration (includes method_config, calc_type, etc.) type(calculation_result_t), intent(out), optional :: result_out type(json_output_data_t), intent(out), optional :: json_data type(calculation_result_t) :: result integer :: total_atoms type(physical_fragment_t) :: full_system type(error_t) :: error integer :: i total_atoms = sys_geom%total_atoms call logger%info("============================================") call logger%info("Running unfragmented calculation") call logger%info(" Total atoms: "//to_char(total_atoms)) call logger%info("============================================") ! Build the full system as a single fragment ! For overlapping fragments, we use the full system directly (not concatenating fragments) full_system%n_atoms = total_atoms full_system%n_caps = 0 allocate (full_system%element_numbers(total_atoms)) allocate (full_system%coordinates(3, total_atoms)) ! Copy all atoms from system geometry full_system%element_numbers = sys_geom%element_numbers full_system%coordinates = sys_geom%coordinates ! Set charge and multiplicity from system full_system%charge = sys_geom%charge full_system%multiplicity = sys_geom%multiplicity call full_system%compute_nelec() ! Validate geometry (check for spatially overlapping atoms) call check_duplicate_atoms(full_system, error) if (error%has_error()) then call logger%error(error%get_full_trace()) error stop "Overlapping atoms in unfragmented system" end if ! Process the full system call do_fragment_work(0_int64, result, config%method_config, phys_frag=full_system, calc_type=config%calc_type) ! Check for calculation errors if (result%has_error) then call logger%error("Unfragmented calculation failed: "//result%error%get_message()) if (present(result_out)) then result_out = result return else error stop "Unfragmented calculation failed" end if end if call logger%info("============================================") call logger%info("Unfragmented calculation completed") block character(len=2048) :: result_line ! Large buffer for Hessian matrix rows integer :: current_log_level, iatom, i, j real(dp) :: hess_norm write (result_line, '(a,f25.15)') " Final energy: ", result%energy%total() call logger%info(trim(result_line)) if (result%has_dipole) then write (result_line, '(a,3f15.8)') " Dipole (e*Bohr): ", result%dipole call logger%info(trim(result_line)) write (result_line, '(a,f15.8)') " Dipole magnitude (Debye): ", norm2(result%dipole)*2.541746_dp call logger%info(trim(result_line)) end if if (result%has_gradient) then write (result_line, '(a,f25.15)') " Gradient norm: ", sqrt(sum(result%gradient**2)) call logger%info(trim(result_line)) ! Print full gradient if verbose and system is small call logger%configuration(level=current_log_level) if (current_log_level >= verbose_level .and. total_atoms < 100) then call logger%info(" ") call logger%info("Gradient (Hartree/Bohr):") do iatom = 1, total_atoms write (result_line, '(a,i5,a,3f20.12)') " Atom ", iatom, ": ", & result%gradient(1, iatom), result%gradient(2, iatom), result%gradient(3, iatom) call logger%info(trim(result_line)) end do call logger%info(" ") end if end if if (result%has_hessian) then ! Compute Frobenius norm of Hessian hess_norm = sqrt(sum(result%hessian**2)) write (result_line, '(a,f25.15)') " Hessian Frobenius norm: ", hess_norm call logger%info(trim(result_line)) ! Print full Hessian if verbose and system is small call logger%configuration(level=current_log_level) if (current_log_level >= verbose_level .and. total_atoms < 20) then call logger%info(" ") call logger%info("Hessian matrix (Hartree/Bohr^2):") do i = 1, 3*total_atoms write (result_line, '(a,i5,a,999f15.8)') " Row ", i, ": ", (result%hessian(i, j), j=1, 3*total_atoms) call logger%info(trim(result_line)) end do call logger%info(" ") end if ! Compute and print vibrational analysis block real(dp), allocatable :: frequencies(:), eigenvalues(:), projected_hessian(:, :) real(dp), allocatable :: reduced_masses(:), force_constants(:) real(dp), allocatable :: cart_disp(:, :), fc_mdyne(:), ir_intensities(:) integer :: ii, jj ! First get projected Hessian for verbose output call logger%info(" Computing vibrational analysis (projecting trans/rot modes)...") call compute_vibrational_frequencies(result%hessian, sys_geom%element_numbers, frequencies, eigenvalues, & coordinates=sys_geom%coordinates, project_trans_rot=.true., & projected_hessian_out=projected_hessian) ! Print projected mass-weighted Hessian if verbose and small system if (current_log_level >= verbose_level .and. total_atoms < 20) then if (allocated(projected_hessian)) then call logger%info(" ") call logger%info("Mass-weighted Hessian after trans/rot projection (a.u.):") do ii = 1, 3*total_atoms write (result_line, '(a,i5,a,999f15.8)') " Row ", ii, ": ", & (projected_hessian(ii, jj), jj=1, 3*total_atoms) call logger%info(trim(result_line)) end do call logger%info(" ") end if end if ! Compute full vibrational analysis and print (with IR intensities if available) if (result%has_dipole_derivatives) then call compute_vibrational_analysis(result%hessian, sys_geom%element_numbers, frequencies, & reduced_masses, force_constants, cart_disp, & coordinates=sys_geom%coordinates, & project_trans_rot=.true., & force_constants_mdyne=fc_mdyne, & dipole_derivatives=result%dipole_derivatives, & ir_intensities=ir_intensities) else call compute_vibrational_analysis(result%hessian, sys_geom%element_numbers, frequencies, & reduced_masses, force_constants, cart_disp, & coordinates=sys_geom%coordinates, & project_trans_rot=.true., & force_constants_mdyne=fc_mdyne) end if if (allocated(frequencies)) then ! Compute thermochemistry for JSON output block type(thermochemistry_result_t) :: thermo_result integer :: n_modes, n_at n_at = size(sys_geom%element_numbers) n_modes = size(frequencies) call compute_thermochemistry(sys_geom%coordinates, sys_geom%element_numbers, & frequencies, n_at, n_modes, thermo_result, & temperature=config%hessian%temperature, pressure=config%hessian%pressure) ! Print vibrational analysis to log if (allocated(ir_intensities)) then call print_vibrational_analysis(frequencies, reduced_masses, force_constants, & cart_disp, sys_geom%element_numbers, & force_constants_mdyne=fc_mdyne, & ir_intensities=ir_intensities, & coordinates=sys_geom%coordinates, & electronic_energy=result%energy%total(), & temperature=config%hessian%temperature, pressure=config%hessian%pressure) else call print_vibrational_analysis(frequencies, reduced_masses, force_constants, & cart_disp, sys_geom%element_numbers, & force_constants_mdyne=fc_mdyne, & coordinates=sys_geom%coordinates, & electronic_energy=result%energy%total(), & temperature=config%hessian%temperature, pressure=config%hessian%pressure) end if ! Populate json_data if present (for centralized JSON output) if (present(json_data)) then json_data%output_mode = OUTPUT_MODE_UNFRAGMENTED json_data%total_energy = result%energy%total() json_data%has_energy = result%has_energy json_data%has_vibrational = .true. ! Copy vibrational data allocate (json_data%frequencies(n_modes)) allocate (json_data%reduced_masses(n_modes)) allocate (json_data%force_constants(n_modes)) json_data%frequencies = frequencies json_data%reduced_masses = reduced_masses json_data%force_constants = fc_mdyne json_data%thermo = thermo_result if (allocated(ir_intensities)) then allocate (json_data%ir_intensities(n_modes)) json_data%ir_intensities = ir_intensities json_data%has_ir_intensities = .true. end if ! Copy dipole if available if (result%has_dipole) then allocate (json_data%dipole(3)) json_data%dipole = result%dipole json_data%has_dipole = .true. end if ! Copy gradient if available if (result%has_gradient) then allocate (json_data%gradient(3, total_atoms)) json_data%gradient = result%gradient json_data%has_gradient = .true. end if ! Copy hessian if available if (result%has_hessian) then allocate (json_data%hessian(3*total_atoms, 3*total_atoms)) json_data%hessian = result%hessian json_data%has_hessian = .true. end if end if if (allocated(ir_intensities)) deallocate (ir_intensities) end block deallocate (frequencies, reduced_masses, force_constants, cart_disp, fc_mdyne) end if if (allocated(eigenvalues)) deallocate (eigenvalues) if (allocated(projected_hessian)) deallocate (projected_hessian) end block end if end block call logger%info("============================================") ! Return result to caller or handle json_data if (present(result_out)) then ! Transfer result to output (for dynamics/optimization) result_out = result else ! Populate json_data for non-Hessian case if present ! (Hessian case already handled above in the vibrational block) if (present(json_data) .and. .not. result%has_hessian) then json_data%output_mode = OUTPUT_MODE_UNFRAGMENTED json_data%total_energy = result%energy%total() json_data%has_energy = result%has_energy if (result%has_dipole) then allocate (json_data%dipole(3)) json_data%dipole = result%dipole json_data%has_dipole = .true. end if if (result%has_gradient) then allocate (json_data%gradient(3, total_atoms)) json_data%gradient = result%gradient json_data%has_gradient = .true. end if end if call result%destroy() end if end subroutine unfragmented_calculation end submodule mqc_unfragmented_workflow