Added initialization of the variables after allocation

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)

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)