Solvation Builder#
import veloxchem as vlx
Create a molecule object for deprotonated ibuprofen and generate a forcefield with semiempirical partial charges to use in the following solvationfep calculations.
ibuprof = vlx.Molecule.read_smiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)[O-]')
ibuprof.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
ff_gen_solute = vlx.MMForceFieldGenerator()
ff_gen_solute.partial_charges = ibuprof.get_partial_charges(
ibuprof.get_charge())
ff_gen_solute.create_topology(ibuprof)
ff_gen_solute.write_gromacs_files('ibuprof','MOL')
* Info * Sum of partial charges is not a whole number.
* Info * Compensating by removing 2.000e-06 from the largest charge.
* Info * Using GAFF (v2.11) parameters.
Reference: J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, J. Comput. Chem. 2004,
25, 1157-1174.
* Info * Updated bond length 13-15 (c -o ) to 0.139 nm
* Info * Updated bond angle 14-13-15 (o -c -o ) to 119.379 deg
Load the SolvationBuilder and create a box of padding 1.0 nm (default). Solvate the molecule with water (SPC/E water model). After solvating the box, the SolvationBuilder will run a NPT equilibration of 5 ps at 300K. Once the equilibration is finished, gromacs files for the equilibrated system is written.
solvator = vlx.SolvationBuilder()
solvator.solvate(ibuprof, equilibrate=True, neutralize=False)
solvator.system_molecule.show()
solvator.write_gromacs_files(solute_ff=ff_gen_solute)
VeloxChem System Builder
==========================
* Info * Solvating the solute with spce molecules
* Info * Padding: 1.0 nm
* Info * NPT Equilibration of the box requested
* Info * The volume of the solute is: 0.23 nm^3
* Info * The box size is: 2.79 x 2.79 x 2.79 nm^3
* Info * The volume available for the solvent is: 21.57 nm^3
* Info * The experimental density of spce is: 1000 kg/m^3
* Info * Warning: 90 attempts have been made to insert the solvent molecule
* Info * Consider reducing the target density
* Info * Warning: 90 attempts have been made to insert the solvent molecule
* Info * Consider reducing the target density
* Info * Solvated system with 718 solvent molecules out of 718 requested
* Info * Time to solvate the system: 0.34 s
* Info * The density of the solvent after packing is: 995 kg/m^3
* Info * The box size after equilibration is: 2.79 x 2.79 x 2.79 nm^3
* Info * The density of the solvent after equilibration is: 1004 kg/m^3
* Info * Equilibrating the system
* Info * Duration: 5.0 ps
* Info * Temperature: 300 K
* Info * Pressure: 1 bar
* Info * Elapsed time to equilibrate the system: 1.31 s
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
* Info * Generating the ForceField for the solvent
* Info * solute.itp file written
* Info * solvent_1.itp file written
* Info * system.top file written
* Info * system.gro file written
Free Energy of Solvation#
The free energy of solvation for ibuprofen in a mix of water and propylene glycol is computed with the SolvationFepDriver using the gromacs files generated from the SolvationBuilder as input. The functon requires a gro and topology file for both the solvated system and the molecule in vacuum.
solvationfep = vlx.SolvationFepDriver()
# Optional, saving trajectory from each lambda simulation in xtc format
#solvationfep.save_trajectory_xtc = True
# The number of steps here has been chosen for a quick execution in the
# Notebook, we recommend using 500 000 steps for production runs (1 ns)
# per Lambda.
# This calculation is commentted out to avoid long execution times in the notebook.
#solvationfep.num_steps = 10000
#delta_f, final_energy =solvationfep.compute_solvation_from_gromacs_files(
# 'system.gro','system.top','ibuprof.gro','ibuprof.top')
Custom solvate#
Following is an example of how to solvate ibuprofen in a mixed solvent of equal amounts of water and propylene glycol. A VeloxChem molecule object is created for solvent molecules. A list of the solvent molecules, the quantity of each molecule and the box dimension (Å) is specified in the function. Gromacs files for the system is written.
solvator = vlx.SolvationBuilder()
propylene_glycol = vlx.Molecule.read_smiles('CC(CO)O')
water = vlx.Molecule.read_smiles('O')
solvator.custom_solvate(ibuprof,
solvents=[water, propylene_glycol],
quantities=[200,200],
box=(40,40,40))
solvator.system_molecule.show()
solvator.write_gromacs_files(solute_ff=ff_gen_solute)
* Info * Solvated system with 200 solvent molecules out of 200 requested
* Info * Solvated system with 200 solvent molecules out of 200 requested
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
* Info * Generating the ForceField for the solvent
* Info * Generating the ForceField for the solvent
* Info * solute.itp file written
* Info * solvent_1.itp file written
* Info * solvent_2.itp file written
* Info * system.top file written
* Info * system.gro file written