# Set minimal log verbosity
setVerbosity(MinimalLog)

# Load the species contained in the molecule database
database = CosmoRSSpeciesDatabase()
acetic_acid = database.exportSpecies('acetic acid')
water = database.exportSpecies('water')
hydronium = database.exportSpecies('hydronium')

# Calculate the CosmoRealSpecies for acetate
# Define the acetate molecule
elements = [Carbon, Hydrogen, Hydrogen, Hydrogen, Carbon, Oxygen, Oxygen]
cartesian_coordinates = [[-0.006526715882, -0.001732896708,  0.010141495754],
                         [ 1.085957422037, -0.002262777684, -0.015712581148],
                         [-0.347930089302, -0.002482824972, -1.028451790718],
                         [-0.354378068774, -0.911147451933,  0.501258230599],
                         [-0.51999582117 ,  1.265644780621,  0.735999158168],
                         [-0.172396533587,  2.381199982022,  0.2424156233  ],
                         [-1.246063552721,  1.095251237566,  1.762408676291]]*Angstrom
configuration = MoleculeConfiguration(
    elements=elements,
    cartesian_coordinates=cartesian_coordinates
)

# Set the DFT-COSMO calculator
basis_set = [
    BasisGGAPseudoDojo.Hydrogen_High,
    BasisGGAPseudoDojo.Carbon_High,
    BasisGGAPseudoDojo.Oxygen_High,
]
poisson_solver = ParallelConjugateGradientSolver(
    boundary_conditions=[
        [MultipoleBoundaryCondition(), MultipoleBoundaryCondition()],
        [MultipoleBoundaryCondition(), MultipoleBoundaryCondition()],
        [MultipoleBoundaryCondition(), MultipoleBoundaryCondition()]
    ]
)
solvation_parameters = CosmoSolvationParameters()
calculator = LCAOCalculator(
    basis_set=basis_set,
    charge=-1.0,
    poisson_solver=poisson_solver,
    solvation_parameters=solvation_parameters
)
configuration.setCalculator(calculator)

# Optimize the geometry
optimized_configuration = OptimizeGeometry(
    configuration=configuration,
    trajectory_filename=None
)

# Calculate the solvation energy
solvation_energy = SolvationEnergy(
    configuration=optimized_configuration
)

# Create the species for acetate
cosmo_rs = CosmoRS(
    configuration=optimized_configuration,
    solvation_energy=solvation_energy
)
acetate = cosmo_rs.evaluate()

# Calculate the pKa of acetic acid
acetic_acid_pka = Acidity(
    reactant_species=[acetic_acid, water],
    product_species=[acetate, hydronium],
    solvent_components=water,
    temperature=298*Kelvin,
    parameters=CosmoRSParameters()
)
pka = acetic_acid_pka.acidity()
nlprint(f'The calculated pKa of acetic acid is {pka}')