import veloxchem as vlxBuilding 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()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.partial_charges = solute.get_partial_charges(solute.get_charge())
ff_gen.create_topology(solute)* Info * Sum of partial charges is not a whole number.
* Info * Compensating by removing 1.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
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) VeloxChem Solvation Builder
=============================
* Info * Solvating the solute with cspce molecules
* Info * Padding: 1.0 nm
* Info * The box size is: 2.78 x 2.78 x 2.78 nm^3
* Info * The volume of the solute is: 0.23 nm^3
* Info * The volume available for the solvent is: 21.35 nm^3
* Info * The experimental density of cspce is: 1000 kg/m^3
* Info * Solvated system with 711 solvent molecules out of 711 requested
* Info * Time to solvate the system: 0.37 s
* Info * The density of the solvent after packing is: 995 kg/m^3
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),
)* Info * Solvated system with 70 solvent molecules out of 70 requested
* 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 70 solvent molecules out of 70 requested
solvator.show_solvation_box()Loading...
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)* Info * The box size after equilibration is: 2.31 x 2.31 x 2.31 nm^3
* Info * The density of the solvent after equilibration is: 172 kg/m^3
* 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
VeloxChem/OpenMM Solvation Free Energy Calculation
==================================================
Simulation parameters:
----------------------
Temperature: 298.15 K
Pressure: 1 atm
Timestep: 2.0 fs
Non-Bonded Cutoff: 1.0 nm
Thermodynamic Ensemble: NPT
Alchemical Parameters:
----------------------
Number of equilibration steps: 500
Number of steps per lambda simulation: 500
Number of snapshots per lambda simulation: 500
Simulation time per lambda: 0.00 ns
Lambdas in Stage 1: [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
Lambdas in Stage 2: [1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.15, 0.1, 0.05, 0.03, 0.0]
Lambdas in Stage 3: [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
Lambdas in Stage 4: [1.0, 0.5, 0.0]
* Info * Starting solvated simulation (Stage 1)...
* Info * Generating systems for stage 1
* Info * Running lambda = 0.0, stage = 1...
* Info * Time spent in energy minimization: 2.59 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 479 with Neff_max = 15.76
* Info * Frames saved = 21 | No subsampling (g=1.40 <= 20)
* Info * Lambda = 0.0 completed. Elapsed time: 1.57 s
* Info * Performance: 2.29 ns/hour
* Info * Running lambda = 0.2, stage = 1...
* Info * Time spent in energy minimization: 2.77 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 341 with Neff_max = 24.24
* Info * Frames saved = 159 | No subsampling (g=6.60 <= 20)
* Info * Lambda = 0.2 completed. Elapsed time: 1.44 s
* Info * Performance: 2.50 ns/hour
* Info * Running lambda = 0.4, stage = 1...
* Info * Time spent in energy minimization: 2.77 s
* Info * Time spent in pre-equilibration: 0.35 s
* Info * Equilibration detected at frame 437 with Neff_max = 30.72
* Info * Frames saved = 63 | No subsampling (g=2.08 <= 20)
* Info * Lambda = 0.4 completed. Elapsed time: 1.49 s
* Info * Performance: 2.41 ns/hour
* Info * Running lambda = 0.6, stage = 1...
* Info * Time spent in energy minimization: 2.65 s
* Info * Time spent in pre-equilibration: 0.38 s
* Info * Equilibration detected at frame 471 with Neff_max = 25.02
* Info * Frames saved = 29 | No subsampling (g=1.20 <= 20)
* Info * Lambda = 0.6 completed. Elapsed time: 1.39 s
* Info * Performance: 2.59 ns/hour
* Info * Running lambda = 0.8, stage = 1...
* Info * Time spent in energy minimization: 2.73 s
* Info * Time spent in pre-equilibration: 0.37 s
* Info * Equilibration detected at frame 413 with Neff_max = 32.13
* Info * Frames saved = 87 | No subsampling (g=2.74 <= 20)
* Info * Lambda = 0.8 completed. Elapsed time: 1.49 s
* Info * Performance: 2.42 ns/hour
* Info * Running lambda = 1.0, stage = 1...
* Info * Time spent in energy minimization: 2.66 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 489 with Neff_max = 12.00
* Info * Frames saved = 11 | No subsampling (g=1.00 <= 20)
* Info * Lambda = 1.0 completed. Elapsed time: 1.51 s
* Info * Performance: 2.39 ns/hour
* Info * Total simulation time for stage 1: 8.89 s. Total simulated time: 0.01 ns.
* Info * Overall performance for stage 1: 2.43 ns/hour
* Info * Recalculating energies for stage 1...
* Info * Recalculating energies for Forcefield 0...
* Info * Forcefield 0 energy recalculation took 0.92 seconds.
* Info * Recalculating energies for Forcefield 1...
* Info * Forcefield 1 energy recalculation took 0.92 seconds.
* Info * Recalculating energies for Forcefield 2...
* Info * Forcefield 2 energy recalculation took 0.92 seconds.
* Info * Recalculating energies for Forcefield 3...
* Info * Forcefield 3 energy recalculation took 0.92 seconds.
* Info * Recalculating energies for Forcefield 4...
* Info * Forcefield 4 energy recalculation took 0.93 seconds.
* Info * Recalculating energies for Forcefield 5...
* Info * Forcefield 5 energy recalculation took 0.92 seconds.
* Info * Energy recalculation for stage 1 took 5.55 seconds.
* Info * Calculating the free energy with MBAR for stage 1...
******* JAX 64-bit mode is now on! *******
* JAX is now set to 64-bit mode! *
* This MAY cause problems with other *
* uses of JAX in the same code. *
******************************************
Free energy for stage 1: -379.7875 +/- 0.6907 kJ/mol
* Info * Removing GSC potential (Stage 2)...
* Info * Generating systems for stage 2
* Info * Running lambda = 1.0, stage = 2...
* Info * Time spent in energy minimization: 2.67 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 411 with Neff_max = 33.32
* Info * Frames saved = 89 | No subsampling (g=2.70 <= 20)
* Info * Lambda = 1.0 completed. Elapsed time: 1.48 s
* Info * Performance: 2.44 ns/hour
* Info * Running lambda = 0.8, stage = 2...
* Info * Time spent in energy minimization: 2.86 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 433 with Neff_max = 18.22
* Info * Frames saved = 67 | No subsampling (g=3.73 <= 20)
* Info * Lambda = 0.8 completed. Elapsed time: 1.52 s
* Info * Performance: 2.37 ns/hour
* Info * Running lambda = 0.6, stage = 2...
* Info * Time spent in energy minimization: 2.90 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 479 with Neff_max = 22.00
* Info * Frames saved = 21 | No subsampling (g=1.00 <= 20)
* Info * Lambda = 0.6 completed. Elapsed time: 1.47 s
* Info * Performance: 2.45 ns/hour
* Info * Running lambda = 0.4, stage = 2...
* Info * Time spent in energy minimization: 2.55 s
* Info * Time spent in pre-equilibration: 0.35 s
* Info * Equilibration detected at frame 361 with Neff_max = 10.79
* Info * Frames saved = 139 | No subsampling (g=12.98 <= 20)
* Info * Lambda = 0.4 completed. Elapsed time: 1.49 s
* Info * Performance: 2.41 ns/hour
* Info * Running lambda = 0.3, stage = 2...
* Info * Time spent in energy minimization: 2.37 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 448 with Neff_max = 14.56
* Info * Frames saved = 52 | No subsampling (g=3.64 <= 20)
* Info * Lambda = 0.3 completed. Elapsed time: 1.41 s
* Info * Performance: 2.55 ns/hour
* Info * Running lambda = 0.2, stage = 2...
* Info * Time spent in energy minimization: 2.60 s
* Info * Time spent in pre-equilibration: 0.38 s
* Info * Equilibration detected at frame 458 with Neff_max = 24.14
* Info * Frames saved = 42 | No subsampling (g=1.78 <= 20)
* Info * Lambda = 0.2 completed. Elapsed time: 1.40 s
* Info * Performance: 2.57 ns/hour
* Info * Running lambda = 0.15, stage = 2...
* Info * Time spent in energy minimization: 2.93 s
* Info * Time spent in pre-equilibration: 0.37 s
* Info * Equilibration detected at frame 424 with Neff_max = 60.09
* Info * Frames saved = 76 | No subsampling (g=1.28 <= 20)
* Info * Lambda = 0.15 completed. Elapsed time: 1.54 s
* Info * Performance: 2.33 ns/hour
* Info * Running lambda = 0.1, stage = 2...
* Info * Time spent in energy minimization: 2.99 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 469 with Neff_max = 18.76
* Info * Frames saved = 31 | No subsampling (g=1.71 <= 20)
* Info * Lambda = 0.1 completed. Elapsed time: 1.31 s
* Info * Performance: 2.75 ns/hour
* Info * Running lambda = 0.05, stage = 2...
* Info * Time spent in energy minimization: 2.98 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 423 with Neff_max = 14.37
* Info * Frames saved = 77 | No subsampling (g=5.43 <= 20)
* Info * Lambda = 0.05 completed. Elapsed time: 1.43 s
* Info * Performance: 2.51 ns/hour
* Info * Running lambda = 0.03, stage = 2...
* Info * Time spent in energy minimization: 2.36 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 441 with Neff_max = 14.37
* Info * Frames saved = 59 | No subsampling (g=4.18 <= 20)
* Info * Lambda = 0.03 completed. Elapsed time: 1.42 s
* Info * Performance: 2.53 ns/hour
* Info * Running lambda = 0.0, stage = 2...
* Info * Time spent in energy minimization: 2.69 s
* Info * Time spent in pre-equilibration: 0.36 s
* Info * Equilibration detected at frame 423 with Neff_max = 9.07
* Info * Frames saved = 77 | No subsampling (g=8.60 <= 20)
* Info * Lambda = 0.0 completed. Elapsed time: 1.45 s
* Info * Performance: 2.48 ns/hour
* Info * Total simulation time for stage 2: 15.93 s. Total simulated time: 0.01 ns.
* Info * Overall performance for stage 2: 2.49 ns/hour
* Info * Recalculating energies for stage 2...
* Info * Recalculating energies for Forcefield 0...
* Info * Forcefield 0 energy recalculation took 1.90 seconds.
* Info * Recalculating energies for Forcefield 1...
* Info * Forcefield 1 energy recalculation took 1.77 seconds.
* Info * Recalculating energies for Forcefield 2...
* Info * Forcefield 2 energy recalculation took 1.78 seconds.
* Info * Recalculating energies for Forcefield 3...
* Info * Forcefield 3 energy recalculation took 1.77 seconds.
* Info * Recalculating energies for Forcefield 4...
* Info * Forcefield 4 energy recalculation took 1.77 seconds.
* Info * Recalculating energies for Forcefield 5...
* Info * Forcefield 5 energy recalculation took 1.81 seconds.
* Info * Recalculating energies for Forcefield 6...
* Info * Forcefield 6 energy recalculation took 1.77 seconds.
* Info * Recalculating energies for Forcefield 7...
* Info * Forcefield 7 energy recalculation took 1.88 seconds.
* Info * Recalculating energies for Forcefield 8...
* Info * Forcefield 8 energy recalculation took 1.79 seconds.
* Info * Recalculating energies for Forcefield 9...
* Info * Forcefield 9 energy recalculation took 1.77 seconds.
* Info * Recalculating energies for Forcefield 10...
* Info * Forcefield 10 energy recalculation took 1.77 seconds.
* Info * Energy recalculation for stage 2 took 19.80 seconds.
* Info * Calculating the free energy with MBAR for stage 2...
Free energy for stage 2: 63.1101 +/- 0.2640 kJ/mol
* Info * Starting vacuum simulation (Stage 3)...
* Info * Generating systems for stage 3
* Info * Running lambda = 0.0, stage = 3...
* Info * Time spent in energy minimization: 0.19 s
* Info * Time spent in pre-equilibration: 0.02 s
* Info * Equilibration detected at frame 268 with Neff_max = 70.63
* Info * Frames saved = 232 | No subsampling (g=3.30 <= 20)
* Info * Lambda = 0.0 completed. Elapsed time: 0.71 s
* Info * Performance: 5.08 ns/hour
* Info * Running lambda = 0.2, stage = 3...
* Info * Time spent in energy minimization: 0.17 s
* Info * Time spent in pre-equilibration: 0.02 s
* Info * Equilibration detected at frame 146 with Neff_max = 102.31
* Info * Frames saved = 354 | No subsampling (g=3.47 <= 20)
* Info * Lambda = 0.2 completed. Elapsed time: 0.70 s
* Info * Performance: 5.12 ns/hour
* Info * Running lambda = 0.4, stage = 3...
* Info * Time spent in energy minimization: 0.19 s
* Info * Time spent in pre-equilibration: 0.01 s
* Info * Equilibration detected at frame 143 with Neff_max = 112.68
* Info * Frames saved = 357 | No subsampling (g=3.18 <= 20)
* Info * Lambda = 0.4 completed. Elapsed time: 0.70 s
* Info * Performance: 5.17 ns/hour
* Info * Running lambda = 0.6, stage = 3...
* Info * Time spent in energy minimization: 0.16 s
* Info * Time spent in pre-equilibration: 0.02 s
* Info * Equilibration detected at frame 316 with Neff_max = 84.78
* Info * Frames saved = 184 | No subsampling (g=2.18 <= 20)
* Info * Lambda = 0.6 completed. Elapsed time: 0.70 s
* Info * Performance: 5.12 ns/hour
* Info * Running lambda = 0.8, stage = 3...
* Info * Time spent in energy minimization: 0.16 s
* Info * Time spent in pre-equilibration: 0.02 s
* Info * Equilibration detected at frame 187 with Neff_max = 92.82
* Info * Frames saved = 313 | No subsampling (g=3.38 <= 20)
* Info * Lambda = 0.8 completed. Elapsed time: 0.70 s
* Info * Performance: 5.13 ns/hour
* Info * Running lambda = 1.0, stage = 3...
* Info * Time spent in energy minimization: 0.17 s
* Info * Time spent in pre-equilibration: 0.01 s
* Info * Equilibration detected at frame 94 with Neff_max = 118.66
* Info * Frames saved = 406 | No subsampling (g=3.43 <= 20)
* Info * Lambda = 1.0 completed. Elapsed time: 0.69 s
* Info * Performance: 5.21 ns/hour
* Info * Total simulation time for stage 3: 4.20 s. Total simulated time: 0.01 ns.
* Info * Overall performance for stage 3: 5.14 ns/hour
* Info * Recalculating energies for stage 3...
* Info * Recalculating energies for Forcefield 0...
* Info * Forcefield 0 energy recalculation took 1.79 seconds.
* Info * Recalculating energies for Forcefield 1...
* Info * Forcefield 1 energy recalculation took 1.78 seconds.
* Info * Recalculating energies for Forcefield 2...
* Info * Forcefield 2 energy recalculation took 1.76 seconds.
* Info * Recalculating energies for Forcefield 3...
* Info * Forcefield 3 energy recalculation took 1.75 seconds.
* Info * Recalculating energies for Forcefield 4...
* Info * Forcefield 4 energy recalculation took 1.80 seconds.
* Info * Recalculating energies for Forcefield 5...
* Info * Forcefield 5 energy recalculation took 1.76 seconds.
* Info * Energy recalculation for stage 3 took 10.63 seconds.
* Info * Calculating the free energy with MBAR for stage 3...
Free energy for stage 3: -15.9119 +/- 0.0806 kJ/mol
* Info * Removing GSC potential in vacuum (Stage 4)...
* Info * Generating systems for stage 4
* Info * Running lambda = 1.0, stage = 4...
* Info * Time spent in energy minimization: 0.17 s
* Info * Time spent in pre-equilibration: 0.01 s
* Info * Equilibration detected at frame 215 with Neff_max = 83.68
* Info * Frames saved = 285 | No subsampling (g=3.42 <= 20)
* Info * Lambda = 1.0 completed. Elapsed time: 0.70 s
* Info * Performance: 5.14 ns/hour
* Info * Running lambda = 0.5, stage = 4...
* Info * Time spent in energy minimization: 0.19 s
* Info * Time spent in pre-equilibration: 0.01 s
* Info * Equilibration detected at frame 225 with Neff_max = 94.33
* Info * Frames saved = 275 | No subsampling (g=2.93 <= 20)
* Info * Lambda = 0.5 completed. Elapsed time: 0.74 s
* Info * Performance: 4.87 ns/hour
* Info * Running lambda = 0.0, stage = 4...
* Info * Time spent in energy minimization: 1.16 s
* Info * Time spent in pre-equilibration: 0.02 s
* Info * Equilibration detected at frame 72 with Neff_max = 102.08
* Info * Frames saved = 428 | No subsampling (g=4.20 <= 20)
* Info * Lambda = 0.0 completed. Elapsed time: 0.69 s
* Info * Performance: 5.21 ns/hour
* Info * Total simulation time for stage 4: 2.13 s. Total simulated time: 0.00 ns.
* Info * Overall performance for stage 4: 5.07 ns/hour
* Info * Recalculating energies for stage 4...
* Info * Recalculating energies for Forcefield 0...
* Info * Forcefield 0 energy recalculation took 0.98 seconds.
* Info * Recalculating energies for Forcefield 1...
* Info * Forcefield 1 energy recalculation took 0.95 seconds.
* Info * Recalculating energies for Forcefield 2...
* Info * Forcefield 2 energy recalculation took 0.98 seconds.
* Info * Energy recalculation for stage 4 took 2.91 seconds.
* Info * Calculating the free energy with MBAR for stage 4...
Free energy for stage 4: 9.4980 +/- 0.1298 kJ/mol
Final free energy: -310.2635 kJ/mol
print(f'Solvation free energy: {fep_results["free_energy"]:.3f} kJ/mol')Solvation free energy: -310.264 kJ/mol