diff --git a/ChemistryLib/CreatePhreeqcIO.cpp b/ChemistryLib/CreatePhreeqcIO.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd5ed3a468a88748d848fc05e77f068a4f8b7530
--- /dev/null
+++ b/ChemistryLib/CreatePhreeqcIO.cpp
@@ -0,0 +1,117 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <fstream>
+
+#include "BaseLib/ConfigTreeUtil.h"
+#include "BaseLib/Error.h"
+#include "CreateOutput.h"
+#include "CreatePhreeqcIO.h"
+#include "PhreeqcIO.h"
+#include "PhreeqcIOData/AqueousSolution.h"
+#include "PhreeqcIOData/EquilibriumPhase.h"
+#include "PhreeqcIOData/KineticReactant.h"
+#include "PhreeqcIOData/ReactionRate.h"
+
+namespace ChemistryLib
+{
+std::unique_ptr<PhreeqcIO> createPhreeqcIO(
+ std::size_t const num_nodes,
+ std::vector<std::pair<int, std::string>> const&
+ process_id_to_component_name_map,
+ boost::optional<BaseLib::ConfigTree> const& config,
+ std::string const& output_directory)
+{
+ auto const& num_chemical_systems = num_nodes;
+
+ // database
+ //! \ogs_file_param{prj__chemical_system__database}
+ auto const database = config->getConfigParameter<std::string>("database");
+ auto path_to_database =
+ BaseLib::joinPaths(BaseLib::getProjectDirectory(), database);
+
+ if (BaseLib::IsFileExisting(path_to_database))
+ {
+ INFO("Found the specified thermodynamic database: %s",
+ path_to_database.c_str());
+ }
+ else
+ {
+ OGS_FATAL("Not found the specified thermodynamicdatabase: %s",
+ path_to_database.c_str());
+ }
+
+ // solution
+ auto const aqueous_solution_per_chem_sys = createAqueousSolution(
+ //! \ogs_file_param{prj__chemical_system__solution}
+ config->getConfigSubtree("solution"));
+
+ auto const& components_per_chem_sys =
+ aqueous_solution_per_chem_sys.components;
+ for (auto const& component : components_per_chem_sys)
+ {
+ auto process_id_to_component_name_map_element = std::find_if(
+ process_id_to_component_name_map.begin(),
+ process_id_to_component_name_map.end(),
+ [&component](std::pair<int, std::string> const& map_element) {
+ return map_element.second == component.name;
+ });
+
+ if (process_id_to_component_name_map_element ==
+ process_id_to_component_name_map.end())
+ OGS_FATAL(
+ "Component %s given in <solution>/<components> is not found in "
+ "specified coupled processes (see "
+ "<process>/<process_variables>/<concentration>).",
+ component.name.c_str());
+ }
+ if (components_per_chem_sys.size() + 1 !=
+ process_id_to_component_name_map.size())
+ OGS_FATAL(
+ "The number of components given in <solution>/<components> is not "
+ "in line with the number of transport processes - 1 which stands "
+ "for the transport process of hydrogen.");
+
+ std::vector<AqueousSolution> aqueous_solutions(
+ num_chemical_systems, aqueous_solution_per_chem_sys);
+
+ // equilibrium phases
+ auto const equilibrium_phases_per_chem_sys = createEquilibriumPhases(
+ //! \ogs_file_param{prj__chemical_system__equilibrium_phases}
+ config->getConfigSubtreeOptional("equilibrium_phases"));
+ std::vector<std::vector<EquilibriumPhase>> equilibrium_phases(
+ num_chemical_systems, equilibrium_phases_per_chem_sys);
+
+ // kinetic reactants
+ auto const kinetic_reactants_per_chem_sys = createKineticReactants(
+ //! \ogs_file_param{prj__chemical_system__kinetic_reactants}
+ config->getConfigSubtreeOptional("kinetic_reactants"));
+ std::vector<std::vector<KineticReactant>> kinetic_reactants(
+ num_chemical_systems, kinetic_reactants_per_chem_sys);
+
+ // rates
+ auto reaction_rates = createReactionRates(
+ //! \ogs_file_param{prj__chemical_system__rates}
+ config->getConfigSubtreeOptional("rates"));
+
+ // output
+ auto const project_file_name = BaseLib::joinPaths(
+ output_directory,
+ BaseLib::extractBaseNameWithoutExtension(config->getProjectFileName()));
+ auto output =
+ createOutput(components_per_chem_sys, equilibrium_phases_per_chem_sys,
+ kinetic_reactants_per_chem_sys, project_file_name);
+
+ return std::make_unique<PhreeqcIO>(
+ project_file_name, std::move(path_to_database),
+ std::move(aqueous_solutions), std::move(equilibrium_phases),
+ std::move(kinetic_reactants), std::move(reaction_rates),
+ std::move(output), process_id_to_component_name_map);
+}
+} // namespace ChemistryLib
diff --git a/ChemistryLib/CreatePhreeqcIO.h b/ChemistryLib/CreatePhreeqcIO.h
new file mode 100644
index 0000000000000000000000000000000000000000..57dd86b9513b1cd5e4040d949c74258c71775b3a
--- /dev/null
+++ b/ChemistryLib/CreatePhreeqcIO.h
@@ -0,0 +1,34 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+#include "PhreeqcIO.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace ProcessLib
+{
+class Process;
+}
+
+namespace ChemistryLib
+{
+std::unique_ptr<PhreeqcIO> createPhreeqcIO(
+ std::size_t const num_nodes,
+ std::vector<std::pair<int, std::string>> const&
+ process_id_to_component_name_map,
+ boost::optional<BaseLib::ConfigTree> const& config,
+ std::string const& output_directory);
+} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
new file mode 100644
index 0000000000000000000000000000000000000000..a88a7f58ff1e5707f509c11ab5cd23758647767f
--- /dev/null
+++ b/ChemistryLib/PhreeqcIO.h
@@ -0,0 +1,59 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
+#include "Output.h"
+#include "PhreeqcIOData/AqueousSolution.h"
+#include "PhreeqcIOData/EquilibriumPhase.h"
+#include "PhreeqcIOData/KineticReactant.h"
+#include "PhreeqcIOData/ReactionRate.h"
+
+namespace ChemistryLib
+{
+class PhreeqcIO
+{
+public:
+ PhreeqcIO(std::string const& project_file_name,
+ std::string&& database,
+ std::vector<AqueousSolution>&& aqueous_solutions,
+ std::vector<std::vector<EquilibriumPhase>>&& equilibrium_phases,
+ std::vector<std::vector<KineticReactant>>&& kinetic_reactants,
+ std::vector<ReactionRate>&& reaction_rates,
+ std::unique_ptr<Output>&& output,
+ std::vector<std::pair<int, std::string>> const&
+ process_id_to_component_name_map)
+ : _phreeqc_input_file(project_file_name + "_phreeqc.inp"),
+ _database(std::move(database)),
+ _aqueous_solutions(std::move(aqueous_solutions)),
+ _equilibrium_phases(std::move(equilibrium_phases)),
+ _kinetic_reactants(std::move(kinetic_reactants)),
+ _reaction_rates(std::move(reaction_rates)),
+ _output(std::move(output)),
+ _process_id_to_component_name_map(process_id_to_component_name_map)
+ {
+ }
+
+ std::string const _phreeqc_input_file;
+
+private:
+ std::string const _database;
+ std::vector<AqueousSolution> _aqueous_solutions;
+ std::vector<std::vector<EquilibriumPhase>> _equilibrium_phases;
+ std::vector<std::vector<KineticReactant>> _kinetic_reactants;
+ std::vector<ReactionRate> const _reaction_rates;
+ std::unique_ptr<Output> const _output;
+ std::vector<std::pair<int, std::string>> const&
+ _process_id_to_component_name_map;
+ double _dt;
+};
+} // namespace ChemistryLib