!********************************************************************** !> Module: basis_projection_mod !> !> Description: !> This module provides routines to project molecular orbitals (MO) and !> density matrices (DM) from a primary basis set to an alternative !> (initial) basis set. !> !> The module includes: !> - A C-binding wrapper subroutine (proj_dm_newbas_C) to interface with !> C codes. !> - The main projection routine (proj_dm_newbas) that performs the MO and !> DM projection, including orthogonalization and density matrix computation. !> - A utility function (itoa) to convert integers to character strings. !********************************************************************** module basis_projection_mod implicit none character(len=*), parameter :: module_name = "basis_projection_mod" contains !********************************************************************** !> C-binding wrapper for the MO/DM projection routine. !********************************************************************** subroutine proj_dm_newbas_C(c_handle) bind(C, name="proj_dm_newbas") use c_interop, only: oqp_handle_t, oqp_handle_get_info use types, only: information type(oqp_handle_t) :: c_handle type(information), pointer :: inf inf => oqp_handle_get_info(c_handle) call proj_dm_newbas(inf) end subroutine proj_dm_newbas_C !********************************************************************** !> Main subroutine for MO and DM projection between basis sets. !> !> This routine projects molecular orbitals (MO) and density matrices (DM) !> from a primary basis set to an alternative (initial) basis set. It performs: !> - Overlap matrix computation and normalization, !> - Corresponding orbital projection, !> - Orbital orthogonalization, !> - Density matrix calculation (for both RHF and ROHF/UHF cases), !> !> Input Data: !> - OQP::VEC_MO_A, OQP::DM_A for the alpha !> - OQP::VEC_MO_B, OQP::DM_B for the beta (if applicable). !> !> Output Data: !> - OQP::VEC_MO_A_tmp, OQP::DM_A_tmp for the alpha spin channel. !> - OQP::VEC_MO_B_tmp, OQP::DM_B_tmp for the beta spin channel (if applicable). !> !> @param[in,out] infos Information structure containing basis sets, atomic data, !> molecular properties, and control parameters. !********************************************************************** subroutine proj_dm_newbas(infos) use precision, only: dp use types, only: information use io_constants, only: IW use oqp_tagarray_driver use basis_tools, only: basis_set use mathlib, only: matrix_invsqrt use util, only: measure_time use messages, only: show_message, WITH_ABORT use printing, only: print_module_info use constants, only: tol_int use int1, only: basis_overlap use iso_c_binding, only: c_char use parallel, only: par_env_t use guess, only: get_ab_initio_density, corresponding_orbital_projection use huckel, only: orthogonalize_orbitals use messages, only: show_message, with_abort implicit none character(len=*), parameter :: subroutine_name = "proj_dm_newbas" type(information), target, intent(inout) :: infos integer :: i, j, nbf, nbf2, nbf_alt, nbf2_alt, nat, nact, ndoc, nproj, l0 type(basis_set), pointer :: basis type(basis_set), pointer :: alt_basis character(len=:), allocatable :: basis_file logical :: err integer , parameter :: root = 0 type(par_env_t) :: pe real(kind=dp), allocatable :: sco(:,:) real(kind=dp), contiguous, pointer :: & Smat(:), q(:,:), & dmat_a(:), mo_a(:,:), mo_energy_a(:), & dmat_b(:), mo_b(:,:), mo_energy_b(:) real(kind=dp), contiguous, pointer :: & Smat_alt(:), & dmat_a_alt(:), mo_a_alt(:,:), mo_energy_a_alt(:), & dmat_b_alt(:), mo_b_alt(:,:), mo_energy_b_alt(:) ! tagarray character(len=*), parameter :: tags_alpha(3) = (/ character(len=80) :: & OQP_DM_A, OQP_E_MO_A, OQP_VEC_MO_A /) character(len=*), parameter :: tags_alpha_tmp(3) = (/ character(len=80) :: & "OQP::DM_A_tmp", "OQP::E_MO_A_tmp", "OQP::VEC_MO_A_tmp" /) character(len=*), parameter :: tags_beta(3) = (/ character(len=80) :: & OQP_DM_B, OQP_E_MO_B, OQP_VEC_MO_B /) character(len=*), parameter :: tags_beta_tmp(3) = (/ character(len=80) :: & "OQP::DM_B_tmp", "OQP::E_MO_B_tmp", "OQP::VEC_MO_B_tmp" /) character(len=*), parameter :: tags_general(1) = (/ character(len=80) :: & OQP_SM /) character(len=1,kind=c_char), contiguous, pointer :: basis_filename(:) ! Files open ! 1. XYZ: Read : Geometric data, ATOMS ! 3. LOG: Read Write: Main output file ! open (unit=IW, file=infos%log_filename, position="append") call print_module_info('Basis Projection', "Projecting MOs and DMs from" & // new_line('') // & "initial to primary basis") ! Readings ! load basis set basis => infos%basis alt_basis => infos%alt_basis call pe%init(infos%mpiinfo%comm, infos%mpiinfo%usempi) alt_basis%atoms => infos%atoms basis%atoms => infos%atoms nbf = basis%nbf nbf2 =nbf*(nbf+1)/2 nat = infos%mol_prop%natom nbf_alt = alt_basis%nbf nbf2_alt =nbf_alt*(nbf_alt+1)/2 if (nbf_alt > nbf) then call show_message("Warning: The initial basis set ("//trim(adjustl(itoa(nbf_alt)))//" functions) " // & "exceeds the primary basis ("//trim(adjustl(itoa(nbf)))//"). " // & "Please select a smaller initial basis set.", WITH_ABORT) end if call data_has_tags(infos%dat, tags_general, module_name, subroutine_name, WITH_ABORT) call tagarray_get_data(infos%dat, OQP_SM, smat) allocate(sco(nbf_alt,nbf), & mo_a(nbf,nbf), & mo_b(nbf,nbf), & q(nbf,nbf), & Dmat_a(nbf2), & Dmat_b(nbf2)) ! Load the converged initial-basis orbitals/densities. These are INPUTS -- ! the projection source, read below through the *_alt pointers -- so they ! must be retrieved, NOT reallocated (alloc_or_die would erase them). call data_has_tags(infos%dat, tags_alpha, module_name, subroutine_name, WITH_ABORT) call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a_alt) call tagarray_get_data(infos%dat, OQP_E_MO_A, mo_energy_a_alt) call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a_alt) call data_has_tags(infos%dat, tags_beta, module_name, subroutine_name, WITH_ABORT) call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b_alt) call tagarray_get_data(infos%dat, OQP_E_MO_B, mo_energy_b_alt) call tagarray_get_data(infos%dat, OQP_VEC_MO_B, mo_b_alt) ! allocate alpha_tmp call infos%dat%alloc_or_die("OQP::DM_A_tmp", (/ nbf2 /), dmat_a, description=OQP_DM_A_comment) call infos%dat%alloc_or_die("OQP::E_MO_A_tmp", (/ nbf /), mo_energy_a, description=OQP_E_MO_A_comment) call infos%dat%alloc_or_die("OQP::VEC_MO_A_tmp", (/ nbf, nbf /), mo_a, description=OQP_VEC_MO_A_comment) ! allocate beta_tmp call infos%dat%alloc_or_die("OQP::DM_B_tmp", (/ nbf2 /), dmat_b, description=OQP_DM_B_comment) call infos%dat%alloc_or_die("OQP::E_MO_B_tmp", (/ nbf /), mo_energy_b, description=OQP_E_MO_B_comment) call infos%dat%alloc_or_die("OQP::VEC_MO_B_tmp", (/ nbf, nbf /), mo_b, description=OQP_VEC_MO_B_comment) ! Dmat_b = 0_dp Dmat_a = 0_dp mo_b = 0_dp mo_a = 0_dp mo_energy_a = 0_dp mo_energy_b = 0_dp call basis_overlap(sco, infos%basis, infos%alt_basis, tol=log(10.0d0)*tol_int) do i =1, nbf sco(:,i) = sco(:,i)*basis%bfnrm(i) * alt_basis%bfnrm(:) end do if (infos%control%scftype == 1) then ndoc = infos%mol_prop%nelec/2 nact = 0 else if (infos%control%scftype >= 2) then ndoc = infos%mol_prop%nelec_b nact = infos%mol_prop%nelec_a-infos%mol_prop%nelec_b end if nproj = nbf_alt l0 = nbf call matrix_invsqrt(smat, q, nbf, qrnk=l0) mo_a(1:nbf,1:nbf) = q(1:nbf,1:nbf) call corresponding_orbital_projection(mo_a_alt, sco, mo_a, ndoc, nact, nproj, nbf, nbf_alt, l0) call orthogonalize_orbitals(q, smat, mo_a, nproj, l0, nbf, nbf) if (infos%control%scftype >= 2) then mo_b(1:nbf,1:nbf) = q(1:nbf,1:nbf) call corresponding_orbital_projection(mo_b_alt, sco, mo_b, ndoc, nact, nproj, nbf, nbf_alt, l0) call orthogonalize_orbitals(q, smat, mo_b, nproj, l0, nbf, nbf) else mo_b = mo_a end if ! Calculate Density Matrix if (pe%rank == root) then ! RHF if (infos%control%scftype == 1) then call get_ab_initio_density(Dmat_A, MO_A, infos=infos, basis=basis) ! ROHF/UHF else call get_ab_initio_density(Dmat_A, MO_A, Dmat_B, MO_B, infos, infos%basis) endif endif ! Broadcast MO and density matrices to all processes call pe%bcast(MO_A, nbf*nbf) if (infos%control%scftype >= 2) then call pe%bcast(MO_B, nbf*nbf) endif ! Broadcast the density matrices to all processes if (infos%control%scftype == 1) then call pe%bcast(Dmat_A, nbf2) else call pe%bcast(Dmat_A, nbf2) call pe%bcast(Dmat_B, nbf2) endif call pe%barrier() write(IW, '(/,a,/)') "...... Completed Basis Projection Computation ......" call measure_time(print_total=1, log_unit=iw) close(iw) end subroutine proj_dm_newbas pure function itoa(num) result(str) integer, intent(in) :: num character(len=16) :: str write(str, '(I0)') num end function itoa end module basis_projection_mod