Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

System solvation

Solvation plays a critical role in determining the structural, thermodynamic, and dynamical properties of molecules in solution. Accurately representing the solvent environment is therefore essential for preparing realistic simulation systems and for computing solution-phase properties such as free energies of solvation.

VeloxChem provides the SolvationBuilder class to construct solvated simulation boxes ready for molecular dynamics. Starting from a solute molecule and its force field topology, a periodic simulation box can be filled with solvent molecules — by default SPC/E water — with optional NPT equilibration and charge neutralization. For more complex environments, mixed solvents with user-defined proportions and box dimensions are also supported.

Beyond system preparation, VeloxChem offers the SolvationFepDriver class to compute free energies of solvation via free energy perturbation (FEP). By alchemically decoupling the solute from the solvent across a series of lambda states, the solvation free energy can be rigorously estimated from the resulting thermodynamic cycle.

Building systems

VeloxChem can be used to solvate system in preparation for molecular dynamics simulations. First, create a molecule object for the solute, in our example we choose deprotonated ibuprofen.

import veloxchem as vlx
solute = vlx.Molecule.read_smiles("CC(C)CC1=CC=C(C=C1)C(C)C(=O)[O-]")
solute.show()
Loading...

Second, generate a force field for the solute. By default RESP charges will be computed but, here, we instead make use of semi-empirical partial charges for faster execution.

ff_gen = vlx.MMForceFieldGenerator()
ff_gen.create_topology(solute)

Third, we use of the SolvationBuilder class to create a simulation box with a solute padding 1.0 nm (default) consisting of SPC/E water (default). The solvator will run an NPT equilibration for 5 ps at 300 K.

solvator = vlx.SolvationBuilder()
solvator.solvate(solute, equilibrate=False, neutralize=False)
solvator.show_solvation_box()
Loading...

Mixed solvents

Solvents can be mixed. Molecule objects are created for each type of solvent molecules.

water = vlx.Molecule.read_smiles("O")
propylene_glycol = vlx.Molecule.read_smiles("CC(CO)O")

Let us consider a solvation with a 50/50 mix of number of water and propylene glycol molecules in the solvent.

solvator.custom_solvate(
    solute,
    solvents=[water, propylene_glycol],
    proportion=(50, 50),
    box_size=(25, 25, 25),
)
solvator.show_solvation_box()
Loading...

Add saving files options.

Free energy of solvation

fep_drv = vlx.SolvationFepDriver()

# Settings are here chosen for a quick execution,
# it is recommended to use the default settings
fep_drv.num_em_steps = 500
fep_drv.num_equil_steps = 500
fep_drv.num_steps = 500

# Pre-equilibrate the solvation box
solvator.perform_equilibration()

# Compute solvation free energy
fep_results = fep_drv.compute_solvation_free_energy(ff_gen, solvator)
print(f'Solvation free energy: {fep_results["free_energy"]:.3f} kJ/mol')
Solvation free energy: -387.851 kJ/mol
Classical methods
Molecular dynamics
Classical methods
Conformational search