System solvation

System solvation#

import veloxchem as vlx

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.

solute = vlx.Molecule.read_smiles("CC(C)CC1=CC=C(C=C1)C(C)C(=O)[O-]")

solute.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

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.partial_charges = solute.get_partial_charges(solute.get_charge())

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()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

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()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

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: -302.626 kJ/mol
Contents