From 6f86ed3a2765009ef968c0e2def0e36ce575af4b Mon Sep 17 00:00:00 2001 From: robin-dahl Date: Thu, 27 Mar 2025 13:35:48 +0100 Subject: [PATCH 1/5] Added initialization of the variables after allocation --- src/ddx_core.f90 | 88 ++++++++++++++++++++++++------------------------ 1 file changed, 44 insertions(+), 44 deletions(-) diff --git a/src/ddx_core.f90 b/src/ddx_core.f90 index 32fd8b23..a6e8c4b4 100644 --- a/src/ddx_core.f90 +++ b/src/ddx_core.f90 @@ -451,24 +451,24 @@ subroutine allocate_state(params, constants, state, ddx_error) type(ddx_error_type), intent(inout) :: ddx_error integer :: istatus - allocate(state % psi(constants % nbasis, params % nsph), stat=istatus) + allocate(state % psi(constants % nbasis, params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `psi` allocation failed") return end if - allocate(state % phi_cav(constants % ncav), stat=istatus) + allocate(state % phi_cav(constants % ncav), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phi_cav` allocation failed") return end if - allocate(state % gradphi_cav(3, constants % ncav), stat=istatus) + allocate(state % gradphi_cav(3, constants % ncav), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, & & "allocate_state: `gradphi_cav` allocation failed") return end if allocate(state % q(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `q` " // & & "allocation failed") @@ -478,55 +478,55 @@ subroutine allocate_state(params, constants, state, ddx_error) ! COSMO model if (params % model .eq. 1) then allocate(state % phi_grid(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phi_grid` " // & & "allocation failed") return end if allocate(state % phi(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phi` " // & & "allocation failed") return end if allocate(state % xs(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `xs` " // & & "allocation failed") return end if allocate(state % xs_rel_diff(params % maxiter), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `xs_rel_diff` " // & & "allocation failed") return end if allocate(state % s(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `s` " // & & "allocation failed") return end if allocate(state % s_rel_diff(params % maxiter), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `s_rel_diff` " // & & "allocation failed") return end if allocate(state % sgrid(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `sgrid` " // & & "allocation failed") return end if - allocate(state % zeta(constants % ncav), stat=istatus) + allocate(state % zeta(constants % ncav), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `zeta` " // & & "allocation failed") @@ -535,111 +535,111 @@ subroutine allocate_state(params, constants, state, ddx_error) ! PCM model else if (params % model .eq. 2) then allocate(state % phi_grid(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phi_grid` " // & & "allocation failed") return end if allocate(state % phi(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phi` " // & & "allocation failed") return end if allocate(state % phiinf(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phiinf` " // & & "allocation failed") return end if allocate(state % phieps(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phieps` " // & & "allocation failed") return end if allocate(state % phieps_rel_diff(params % maxiter), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `xs_rel_diff` " // & & "allocation failed") return end if allocate(state % xs(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `xs` " // & & "allocation failed") return end if allocate(state % xs_rel_diff(params % maxiter), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `xs_rel_diff` " // & & "allocation failed") return end if allocate(state % s(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `s` " // & & "allocation failed") return end if allocate(state % s_rel_diff(params % maxiter), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `xs_rel_diff` " // & & "allocation failed") return end if allocate(state % sgrid(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `sgrid` " // & & "allocation failed") return end if allocate(state % y(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `y` " // & & "allocation failed") return end if allocate(state % y_rel_diff(params % maxiter), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `y_rel_diff` " // & & "allocation failed") return end if allocate(state % ygrid(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `ygrid` " // & & "allocation failed") return end if allocate(state % g(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `g` " // & & "allocation failed") return end if allocate(state % qgrid(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `qgrid` " // & & "allocation failed") return end if - allocate(state % zeta(constants % ncav), stat=istatus) + allocate(state % zeta(constants % ncav), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `zeta` " // & & "allocation failed") @@ -648,20 +648,20 @@ subroutine allocate_state(params, constants, state, ddx_error) ! LPB model else if (params % model .eq. 3) then allocate(state % phi_grid(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phi_grid` " // & & "allocation failed") return end if allocate(state % phi(constants % nbasis, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phi` " // & & "allocation failed") return end if - allocate(state % zeta(constants % ncav), stat=istatus) + allocate(state % zeta(constants % ncav), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `zeta` " // & & "allocation failed") @@ -669,7 +669,7 @@ subroutine allocate_state(params, constants, state, ddx_error) return end if allocate(state % x_lpb_rel_diff(params % maxiter), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `x_lpb_rel_diff` " // & & "allocation failed") @@ -677,7 +677,7 @@ subroutine allocate_state(params, constants, state, ddx_error) end if allocate(state % rhs_lpb(constants % nbasis, & & params % nsph, 2), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `rhs_lpb` " // & & "allocation failed") @@ -685,7 +685,7 @@ subroutine allocate_state(params, constants, state, ddx_error) end if allocate(state % rhs_adj_lpb(constants % nbasis, & & params % nsph, 2), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `rhs_adj_lpb` " // & & "allocation failed") @@ -693,69 +693,69 @@ subroutine allocate_state(params, constants, state, ddx_error) end if allocate(state % x_lpb(constants % nbasis, & & params % nsph, 2), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `x_lpb` " // & & "allocation failed") return end if allocate(state % x_adj_lpb(constants % nbasis, & - & params % nsph, 2), stat=istatus) + & params % nsph, 2), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `x_adj_lpb` " // & & "allocation failed") return end if allocate(state % x_adj_lpb_rel_diff(params % maxiter), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, & & "allocate_state: `x_adj_lpb_rel_diff` allocation failed") return end if allocate(state % g_lpb(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `g_lpb` " // & & "allocation failed") return end if allocate(state % f_lpb(params % ngrid, & - & params % nsph), stat=istatus) + & params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `f_lpb` " // & & "allocation failed") return end if - allocate(state % zeta_dip(3, constants % ncav), stat=istatus) + allocate(state % zeta_dip(3, constants % ncav), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `zeta_dip` " // & & "allocation failed") return end if allocate(state % x_adj_re_grid(params % ngrid, params % nsph), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `x_adj_re_grid` " // & & "allocation failed") return end if allocate(state % x_adj_r_grid(params % ngrid, params % nsph), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `x_adj_r_grid` " // & & "allocation failed") return end if allocate(state % x_adj_e_grid(params % ngrid, params % nsph), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `x_adj_e_grid` " // & & "allocation failed") return end if allocate(state % phi_n(params % ngrid, params % nsph), & - & stat=istatus) + & stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `phi_n` " // & & "allocation failed") From de50d817af8cd8029f50c324d5f7e062ae935688 Mon Sep 17 00:00:00 2001 From: robin-dahl Date: Thu, 27 Mar 2025 13:38:33 +0100 Subject: [PATCH 2/5] Added reset of the hsp guess to zero during the LPB guess adjoint routine --- src/ddx_operators.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ddx_operators.f90 b/src/ddx_operators.f90 index dedbf00f..8c26f38d 100644 --- a/src/ddx_operators.f90 +++ b/src/ddx_operators.f90 @@ -863,6 +863,7 @@ subroutine prec_tstarx(params, constants, workspace, x, y, ddx_error) start_time = omp_get_wtime() n_iter = params % maxiter + workspace % hsp_guess = zero call jacobi_diis(params, constants, workspace, constants % inner_tol, & & x(:,:,2), workspace % hsp_guess, n_iter, x_rel_diff, bstarx, & & bx_prec, hnorm, ddx_error) From d9870b40b5b222ba13ced0833894fc55e036b7fb Mon Sep 17 00:00:00 2001 From: Michele Nottoli Date: Fri, 8 Aug 2025 14:41:47 +0200 Subject: [PATCH 3/5] Added checks on the state (on the fortran side). --- src/ddx.f90 | 84 +++++++++++++++++++++++++++++++++++++++++++ src/ddx_core.f90 | 15 ++++++++ src/ddx_cosmo.f90 | 2 ++ src/ddx_lpb.f90 | 5 +++ src/ddx_pcm.f90 | 2 ++ src/pyddx_classes.cpp | 50 +++----------------------- 6 files changed, 113 insertions(+), 45 deletions(-) diff --git a/src/ddx.f90 b/src/ddx.f90 index 4b67ef44..550e5618 100644 --- a/src/ddx.f90 +++ b/src/ddx.f90 @@ -515,6 +515,12 @@ subroutine fill_guess(params, constants, workspace, state, tol, ddx_error) real(dp), intent(in) :: tol type(ddx_error_type), intent(inout) :: ddx_error + if (.not.state % rhs_done) then + call update_error(ddx_error, & + & "In fill_guess, the RHS is not initialized") + return + end if + if (params % model .eq. 1) then call cosmo_guess(params, constants, workspace, state, ddx_error) else if (params % model .eq. 2) then @@ -526,6 +532,11 @@ subroutine fill_guess(params, constants, workspace, state, tol, ddx_error) return end if + if (ddx_error % flag .eq. 0) then + state % guess_done = .true. + state % solved = .false. + end if + end subroutine fill_guess !> Do a guess for the adjoint linear system for the different models @@ -547,6 +558,12 @@ subroutine fill_guess_adjoint(params, constants, workspace, state, tol, ddx_erro real(dp), intent(in) :: tol type(ddx_error_type), intent(inout) :: ddx_error + if (.not.state % adjoint_rhs_done) then + call update_error(ddx_error, & + & "In fill_guess_adjoint, the adjoint RHS is not initialized") + return + end if + if (params % model .eq. 1) then call cosmo_guess_adjoint(params, constants, workspace, state, ddx_error) else if (params % model .eq. 2) then @@ -558,6 +575,11 @@ subroutine fill_guess_adjoint(params, constants, workspace, state, tol, ddx_erro return end if + if (ddx_error % flag .eq. 0) then + state % adjoint_guess_done = .true. + state % adjoint_solved = .false. + end if + end subroutine fill_guess_adjoint !> Solve the primal linear system for the different models @@ -579,6 +601,16 @@ subroutine solve(params, constants, workspace, state, tol, ddx_error) real(dp), intent(in) :: tol type(ddx_error_type), intent(inout) :: ddx_error + if (.not.state % rhs_done) then + call update_error(ddx_error, "In solve, the RHS is not initialized") + return + end if + if (.not.(state % solved .or. state % guess_done)) then + call update_error(ddx_error, & + & "In solve, no guess or previous solution is provided") + return + end if + if (params % model .eq. 1) then call cosmo_solve(params, constants, workspace, state, tol, ddx_error) else if (params % model .eq. 2) then @@ -590,6 +622,10 @@ subroutine solve(params, constants, workspace, state, tol, ddx_error) return end if + if (ddx_error % flag .eq. 0) then + state % solved = .true. + end if + end subroutine solve !> Solve the adjoint linear system for the different models @@ -611,6 +647,17 @@ subroutine solve_adjoint(params, constants, workspace, state, tol, ddx_error) real(dp), intent(in) :: tol type(ddx_error_type), intent(inout) :: ddx_error + if (.not.state % adjoint_rhs_done) then + call update_error(ddx_error, & + & "In solve_adjoint, the RHS is not initialized") + return + end if + if (.not.(state % adjoint_solved .or. state % adjoint_guess_done)) then + call update_error(ddx_error, & + & "In solve_adjoint, no guess or previous solution is provided") + return + end if + if (params % model .eq. 1) then call cosmo_solve_adjoint(params, constants, workspace, state, tol, ddx_error) else if (params % model .eq. 2) then @@ -622,6 +669,10 @@ subroutine solve_adjoint(params, constants, workspace, state, tol, ddx_error) return end if + if (ddx_error % flag .eq. 0) then + state % adjoint_solved = .true. + end if + end subroutine solve_adjoint !> Compute the energy for the different models @@ -646,6 +697,18 @@ subroutine energy(params, constants, workspace, state, solvation_energy, ddx_err ! dummy operation on unused interface arguments if (allocated(workspace % tmp_pot)) continue + if (.not.state % solved) then + call update_error(ddx_error, & + & "In energy, the solution is not available.") + return + end if + if (.not.state % adjoint_rhs_done) then + call update_error(ddx_error, & + & "In energy, the adjoint RHS is not initialized") + return + end if + + if (params % model .eq. 1) then call cosmo_energy(constants, state, solvation_energy, ddx_error) else if (params % model .eq. 2) then @@ -683,6 +746,27 @@ subroutine solvation_force_terms(params, constants, workspace, & real(dp), intent(out) :: force(3, params % nsph) type(ddx_error_type), intent(inout) :: ddx_error + if (.not.state % solved) then + call update_error(ddx_error, & + & "In solvation_force_terms, the solution is not available.") + return + end if + if (.not.state % adjoint_solved) then + call update_error(ddx_error, & + & "In solvation_force_terms, the adjoint solution is not available.") + return + end if + if (.not.state % rhs_done) then + call update_error(ddx_error, & + & "In solvation_force_terms, the RHS is not initialized") + return + end if + if (.not.state % adjoint_rhs_done) then + call update_error(ddx_error, & + & "In solvation_force_terms, the adjoint RHS is not initialized") + return + end if + if (params % model .eq. 1) then call cosmo_solvation_force_terms(params, constants, workspace, & & state, electrostatics % e_cav, force, ddx_error) diff --git a/src/ddx_core.f90 b/src/ddx_core.f90 index a6e8c4b4..0cf74bcf 100644 --- a/src/ddx_core.f90 +++ b/src/ddx_core.f90 @@ -186,6 +186,14 @@ module ddx_core real(dp) :: hsp_time real(dp) :: hsp_adj_time + !! Some flags to see if quantities are properly initialized + logical :: rhs_done + logical :: adjoint_rhs_done + logical :: guess_done + logical :: solved + logical :: adjoint_guess_done + logical :: adjoint_solved + end type ddx_state_type !> Container for the electrostatic properties. Since different methods @@ -451,6 +459,13 @@ subroutine allocate_state(params, constants, state, ddx_error) type(ddx_error_type), intent(inout) :: ddx_error integer :: istatus + state % rhs_done = .false. + state % adjoint_rhs_done = .false. + state % guess_done = .false. + state % solved = .false. + state % adjoint_guess_done = .false. + state % adjoint_solved = .false. + allocate(state % psi(constants % nbasis, params % nsph), stat=istatus, source=0.0_dp) if (istatus .ne. 0) then call update_error(ddx_error, "allocate_state: `psi` allocation failed") diff --git a/src/ddx_cosmo.f90 b/src/ddx_cosmo.f90 index d8047e63..1bb3b19a 100644 --- a/src/ddx_cosmo.f90 +++ b/src/ddx_cosmo.f90 @@ -73,6 +73,8 @@ subroutine cosmo_setup(params, constants, workspace, state, phi_cav, psi, ddx_er state % phi = - state % phi state % phi_cav = phi_cav state % psi = psi + state % rhs_done = .true. + state % adjoint_rhs_done = .true. end subroutine cosmo_setup !> Do a guess for the primal ddCOSMO linear system diff --git a/src/ddx_lpb.f90 b/src/ddx_lpb.f90 index a8e31103..051d21a9 100644 --- a/src/ddx_lpb.f90 +++ b/src/ddx_lpb.f90 @@ -90,6 +90,11 @@ subroutine lpb_setup(params, constants, workspace, state, phi_cav, & & constants % vgrid_nbasis, state % f_lpb, state % rhs_lpb(:,:,2)) state % rhs_lpb(:,:,1) = state % rhs_lpb(:,:,1) + state % rhs_lpb(:,:,2) + if (ddx_error % flag .eq. 0) then + state % rhs_done = .true. + state % adjoint_rhs_done = .true. + end if + end subroutine lpb_setup !> Compute the ddLPB energy diff --git a/src/ddx_pcm.f90 b/src/ddx_pcm.f90 index 199424c2..d84e36a5 100644 --- a/src/ddx_pcm.f90 +++ b/src/ddx_pcm.f90 @@ -70,6 +70,8 @@ subroutine pcm_setup(params, constants, workspace, state, phi_cav, psi, ddx_erro state % phi = - state % phi state % phi_cav = phi_cav state % psi = psi + state % rhs_done = .true. + state % adjoint_rhs_done = .true. end subroutine pcm_setup !> Do a guess for the primal ddPCM linear system diff --git a/src/pyddx_classes.cpp b/src/pyddx_classes.cpp index 513875af..c4816904 100644 --- a/src/pyddx_classes.cpp +++ b/src/pyddx_classes.cpp @@ -359,59 +359,27 @@ class State { } void fill_guess(double tol) { - if (model()->model() == "cosmo") { - ddx_cosmo_guess(model()->holder(), holder(), model()->error()); - } else if (model()->model() == "pcm") { - ddx_pcm_guess(model()->holder(), holder(), model()->error()); - } else if (model()->model() == "lpb") { - ddx_lpb_guess(model()->holder(), holder(), tol, model()->error()); - } else { - throw py::value_error("Model " + model()->model() + " not yet implemented."); - } + ddx_fill_guess(model()->holder(), holder(), tol, model()->error()); throw_if_error(); m_solved = false; } void fill_guess_adjoint(double tol) { - if (model()->model() == "cosmo") { - ddx_cosmo_guess_adjoint(model()->holder(), holder(), model()->error()); - } else if (model()->model() == "pcm") { - ddx_pcm_guess_adjoint(model()->holder(), holder(), model()->error()); - } else if (model()->model() == "lpb") { - ddx_lpb_guess_adjoint(model()->holder(), holder(), tol, model()->error()); - } else { - throw py::value_error("Model " + model()->model() + " not yet implemented."); - } + ddx_fill_guess_adjoint(model()->holder(), holder(), tol, model()->error()); throw_if_error(); m_solved_adjoint = false; } // Solve the forward COSMO / PCM / LPB System. The state is modified in-place. void solve(double tol) { - if (model()->model() == "cosmo") { - ddx_cosmo_solve(model()->holder(), holder(), tol, model()->error()); - } else if (model()->model() == "pcm") { - ddx_pcm_solve(model()->holder(), holder(), tol, model()->error()); - } else if (model()->model() == "lpb") { - ddx_lpb_solve(model()->holder(), holder(), tol, model()->error()); - } else { - throw py::value_error("Model " + model()->model() + " not yet implemented."); - } + ddx_solve(model()->holder(), holder(), tol, model()->error()); throw_if_error(); m_solved = true; } // Solve the adjoint COSMO / PCM / LPB System. The state is modified in-place. void solve_adjoint(double tol) { - if (model()->model() == "cosmo") { - ddx_cosmo_solve_adjoint(model()->holder(), holder(), tol, model()->error()); - } else if (model()->model() == "pcm") { - ddx_pcm_solve_adjoint(model()->holder(), holder(), tol, model()->error()); - } else if (model()->model() == "lpb") { - ddx_lpb_solve_adjoint(model()->holder(), holder(), tol, model()->error()); - } else { - throw py::value_error("Model " + model()->model() + " not yet implemented."); - } + ddx_solve_adjoint(model()->holder(), holder(), tol, model()->error()); throw_if_error(); m_solved_adjoint = true; } @@ -420,15 +388,7 @@ class State { double energy() { check_solved(); double ret = 0.0; - if (model()->model() == "cosmo") { - ret = ddx_cosmo_energy(model()->holder(), holder(), model()->error()); - } else if (model()->model() == "pcm") { - ret = ddx_pcm_energy(model()->holder(), holder(), model()->error()); - } else if (model()->model() == "lpb") { - ret = ddx_lpb_energy(model()->holder(), holder(), model()->error()); - } else { - throw py::value_error("Model " + model()->model() + " not yet implemented."); - } + ret = ddx_energy(model()->holder(), holder(), model()->error()); throw_if_error(); return ret; } From 8d6038d2fcbc39c3391567e868cf9cac05844157 Mon Sep 17 00:00:00 2001 From: Michele Nottoli Date: Fri, 8 Aug 2025 15:01:22 +0200 Subject: [PATCH 4/5] Added two new scratch arrays to be used in prec_tstarx. --- src/ddx_lpb.f90 | 4 ++-- src/ddx_operators.f90 | 9 ++++----- src/ddx_workspace.f90 | 20 ++++++++++++++++++-- 3 files changed, 24 insertions(+), 9 deletions(-) diff --git a/src/ddx_lpb.f90 b/src/ddx_lpb.f90 index 051d21a9..47bd4730 100644 --- a/src/ddx_lpb.f90 +++ b/src/ddx_lpb.f90 @@ -173,8 +173,8 @@ subroutine lpb_guess_adjoint(params, constants, workspace, state, tol, & constants % inner_tol = sqrt(tol) ! guess - workspace % ddcosmo_guess = zero - workspace % hsp_guess = zero + workspace % ddcosmo_adj_guess = zero + workspace % hsp_adj_guess = zero call prec_tstarx(params, constants, workspace, state % rhs_adj_lpb, & & state % x_adj_lpb, ddx_error) diff --git a/src/ddx_operators.f90 b/src/ddx_operators.f90 index 8c26f38d..90f7df33 100644 --- a/src/ddx_operators.f90 +++ b/src/ddx_operators.f90 @@ -851,28 +851,27 @@ subroutine prec_tstarx(params, constants, workspace, x, y, ddx_error) call convert_ddcosmo(params, constants, 1, y(:,:,1)) n_iter = params % maxiter call jacobi_diis(params, constants, workspace, constants % inner_tol, & - & y(:,:,1), workspace % ddcosmo_guess, n_iter, x_rel_diff, lstarx, & + & y(:,:,1), workspace % ddcosmo_adj_guess, n_iter, x_rel_diff, lstarx, & & ldm1x, hnorm, ddx_error) if (ddx_error % flag .ne. 0) then call update_error(ddx_error, 'prec_tstarx: ddCOSMO failed to ' // & & 'converge, exiting') return end if - y(:,:,1) = workspace % ddcosmo_guess + y(:,:,1) = workspace % ddcosmo_adj_guess workspace % s_time = workspace % s_time + omp_get_wtime() - start_time start_time = omp_get_wtime() n_iter = params % maxiter - workspace % hsp_guess = zero call jacobi_diis(params, constants, workspace, constants % inner_tol, & - & x(:,:,2), workspace % hsp_guess, n_iter, x_rel_diff, bstarx, & + & x(:,:,2), workspace % hsp_adj_guess, n_iter, x_rel_diff, bstarx, & & bx_prec, hnorm, ddx_error) if (ddx_error % flag .ne. 0) then call update_error(ddx_error, 'prec_tstarx: HSP failed to ' // & & 'converge, exiting') return end if - y(:,:,2) = workspace % hsp_guess + y(:,:,2) = workspace % hsp_adj_guess workspace % hsp_adj_time = workspace % hsp_adj_time & & + omp_get_wtime() - start_time diff --git a/src/ddx_workspace.f90 b/src/ddx_workspace.f90 index c9df8ac5..0fffced3 100644 --- a/src/ddx_workspace.f90 +++ b/src/ddx_workspace.f90 @@ -97,6 +97,7 @@ module ddx_workspace real(dp), allocatable :: tmp_bmat(:, :) !> ddLPB solutions for the microiterations real(dp), allocatable :: ddcosmo_guess(:,:), hsp_guess(:,:) + real(dp), allocatable :: ddcosmo_adj_guess(:,:), hsp_adj_guess(:,:) !> ddLPB temporary timings real(dp) :: xs_time, s_time, hsp_time, hsp_adj_time end type ddx_workspace_type @@ -267,8 +268,11 @@ subroutine workspace_init(params, constants, workspace, ddx_error) ! Allocations for LPB model if (params % model .eq. 3) then allocate(workspace % tmp_bessel(max(2, params % lmax+1), & - & params % nproc), workspace % ddcosmo_guess(constants % nbasis, params % nsph), & - & workspace % hsp_guess(constants % nbasis, params % nsph), stat=info) + & params % nproc), & + & workspace % ddcosmo_guess(constants % nbasis, params % nsph), & + & workspace % ddcosmo_adj_guess(constants % nbasis, params % nsph), & + & workspace % hsp_guess(constants % nbasis, params % nsph), & + & workspace % hsp_adj_guess(constants % nbasis, params % nsph), stat=info) if (info .ne. 0) then call update_error(ddx_error, "workspace_init: `tmp_bessel` " // & & "allocation failed") @@ -452,6 +456,18 @@ subroutine workspace_free(workspace, ddx_error) call update_error(ddx_error, "`hsp_guess` deallocation failed!") end if end if + if (allocated(workspace % ddcosmo_adj_guess)) then + deallocate(workspace % ddcosmo_adj_guess, stat=istat) + if (istat .ne. 0) then + call update_error(ddx_error, "`ddcosmo_adj_guess` deallocation failed!") + end if + end if + if (allocated(workspace % hsp_adj_guess)) then + deallocate(workspace % hsp_adj_guess, stat=istat) + if (istat .ne. 0) then + call update_error(ddx_error, "`hsp_adj_guess` deallocation failed!") + end if + end if if (allocated(workspace % tmp_rhs)) then deallocate(workspace % tmp_rhs, stat=istat) if (istat .ne. 0) then From f802bdfc45c9ea9019778d958932adb9eee9e12d Mon Sep 17 00:00:00 2001 From: Michele Nottoli Date: Fri, 8 Aug 2025 15:28:30 +0200 Subject: [PATCH 5/5] Fixed an unrelated bug in the C interface (ddx_ddrun). --- src/ddx_cinterface.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ddx_cinterface.f90 b/src/ddx_cinterface.f90 index 723b17e4..9530e537 100644 --- a/src/ddx_cinterface.f90 +++ b/src/ddx_cinterface.f90 @@ -553,7 +553,7 @@ function ddx_ddrun(c_ddx, c_state, c_electrostatics, nbasis, nsph, & c_energy = zero ! setup - do_guess = read_guess.ne.0 + do_guess = read_guess.eq.0 call c_f_pointer(c_ddx, ddx_model) call ddx_setup(c_ddx, c_state, c_electrostatics, nbasis, nsph, psi, c_error) if (ddx_get_error_flag(c_error) .ne. 0) return @@ -573,7 +573,7 @@ function ddx_ddrun(c_ddx, c_state, c_electrostatics, nbasis, nsph, & ! adjoint linear system if (ddx_model%params%force .eq. 1) then if (do_guess) then - call ddx_fill_guess(c_ddx, c_state, tol, c_error) + call ddx_fill_guess_adjoint(c_ddx, c_state, tol, c_error) end if if (ddx_get_error_flag(c_error) .ne. 0) return call ddx_solve_adjoint(c_ddx, c_state, tol, c_error)