Calculation of Self Solvation Free Energy

Do Imports

[1]:
import os
from rdkit import Chem
import pygromos
from pygromos.files.forcefield.openff.openff import OpenFF
from pygromos.files.gromos_system.gromos_system import Gromos_System
from pygromos.simulations.approaches.solvation_free_energy_calculation.solvation_free_energy import Solvation_free_energy_calculation
from pygromos.simulations.hpc_queuing.submission_systems.local import LOCAL
Warning: Unable to load toolkit 'OpenEye Toolkit'. The Open Force Field Toolkit does not require the OpenEye Toolkits, and can use RDKit/AmberTools instead. However, if you have a valid license for the OpenEye Toolkits, consider installing them for faster performance and additional file format support: https://docs.eyesopen.com/toolkits/python/quickstart-python/linuxosx.html OpenEye offers free Toolkit licenses for academics: https://www.eyesopen.com/academic-licensing

Choose Molecule to run calculation for

[2]:
smiles = "c1ccccc1"
workfolder = project_dir = os.path.abspath("example_files/solvation_free_energy")

create the gromos_system from a smile and get the number of atoms

[3]:
groSys = Gromos_System(work_folder=workfolder, system_name="test", in_smiles=smiles, forcefield=OpenFF(), auto_convert=True)
number_of_atoms = groSys.mol.GetNumAtoms()
print("Number of atoms:", number_of_atoms)
Found off: /home/mlehner/PyGromosTools/pygromos/data/ff/SMIRNOFF/openff_unconstrained-2.0.0.offxml
Found off: /home/mlehner/PyGromosTools/pygromos/data/ff/SMIRNOFF/openff_unconstrained-2.0.0.offxml
Number of atoms: 12
[4]:
subSys = LOCAL() # use the local submission system (for cluster use LSF)
[5]:
n_points = 5 # Number of Lambda points to calculate (typically 21)

Create The Solvation free energy calculation system

[6]:
sf = Solvation_free_energy_calculation(input_system=groSys, # Gromos_System, SMILES (str) or rdkit Mol
                                       work_folder=workfolder, # Folder to do calculations in
                                       system_name="test", # Name of the system (does not need to be smiles but convenient)
                                       forcefield=OpenFF(), # Force field to use
                                       density=789, # density of the liquid in kg/L
                                       num_molecules=512, # number of molecules used for the calculation
                                       num_atoms=number_of_atoms, # number of atoms in one molecule
                                       subsystem=subSys, # Subsystem to use for calculation local or lsf
                                       amberscaling=False, # Whether to use amberscaling (for openforcefield recommended)
                                       n_points=n_points) # Number of Lambda points to calculate (typically 21)
Found off: /home/mlehner/PyGromosTools/pygromos/data/ff/SMIRNOFF/openff_unconstrained-2.0.0.offxml
File imd is empty , can not be written!
/home/mlehner/PyGromosTools/pygromos/simulations/approaches/solvation_free_energy_calculation/solvation_free_energy.py:146: UserWarning: Folder does already exist
  warnings.warn("Folder does already exist")
/home/mlehner/PyGromosTools/pygromos/files/gromos_system/gromos_system.py:869: UserWarning: Did not change file path as its only promised None
  warnings.warn("Did not change file path as its only promised " + str(file_obj.path))

Create Liquid

[7]:
sf.create_liq()
/home/mlehner/PyGromosTools/pygromos/simulations/approaches/solvation_free_energy_calculation/solvation_free_energy.py:225: UserWarning: Folder does already exist
  warnings.warn("Folder does already exist")

Minimize Liquid

[8]:
emin_sys, jobID = sf.minimize_liq(in_gromos_simulation_system=sf.groSys_liq,prev_JobID=-1)
################################################################################

test_emin
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/test_emin/analysis/data/test_emin.cnf
 GOING TO SKIPT THIS SUBMISSION!

Change the number of cores for longer runs

[9]:
sf.subsystem.nomp = 6

Equilibrate System

[10]:
eq_sys, jobID = sf.eq_liq(in_gromos_simulation_system=emin_sys,prev_JobID=jobID)
################################################################################

eqtest_eq_1
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/eqtest_eq_1/analysis/data/eqtest_eq_1.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

eqtest_eq_2
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/eqtest_eq_2/analysis/data/eqtest_eq_2.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

eqtest_eq_3
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/eqtest_eq_3/analysis/data/eqtest_eq_3.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

eqtest_eq_4
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/eqtest_eq_4/analysis/data/eqtest_eq_4.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

eqtest_eq_5
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/eqtest_eq_5/analysis/data/eqtest_eq_5.cnf
 GOING TO SKIPT THIS SUBMISSION!

Do TI calculation

[11]:
ti_sys, jobID = sf.ti_liq(in_gromos_simulation_system=eq_sys,
                          prev_JobID=jobID,
                          n_points=n_points)
/home/mlehner/PyGromosTools/pygromos/simulations/approaches/solvation_free_energy_calculation/solvation_free_energy.py:481: UserWarning: Folder does already exist
  warnings.warn("Folder does already exist")
################################################################################

test_L0.0_0
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.0_0/test_L0.0_0/analysis/data/test_L0.0_0.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.0_1
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.0_1/test_L0.0_1/analysis/data/test_L0.0_1.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.0_2
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.0_2/test_L0.0_2/analysis/data/test_L0.0_2.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.25_0
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.25_0/test_L0.25_0/analysis/data/test_L0.25_0.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.25_1
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.25_1/test_L0.25_1/analysis/data/test_L0.25_1.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.25_2
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.25_2/test_L0.25_2/analysis/data/test_L0.25_2.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.5_0
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.5_0/test_L0.5_0/analysis/data/test_L0.5_0.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.5_1
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.5_1/test_L0.5_1/analysis/data/test_L0.5_1.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.5_2
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.5_2/test_L0.5_2/analysis/data/test_L0.5_2.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.75_0
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.75_0/test_L0.75_0/analysis/data/test_L0.75_0.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.75_1
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.75_1/test_L0.75_1/analysis/data/test_L0.75_1.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L0.75_2
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L0.75_2/test_L0.75_2/analysis/data/test_L0.75_2.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L1.0_0
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L1.0_0/test_L1.0_0/analysis/data/test_L1.0_0.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L1.0_1
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L1.0_1/test_L1.0_1/analysis/data/test_L1.0_1.cnf
 GOING TO SKIPT THIS SUBMISSION!
################################################################################

test_L1.0_2
################################################################################

============================================================
FOUND RESULT: /home/mlehner/PyGromosTools/docs/sphinx_project/Examples/example_files/solvation_free_energy/test_liq/ti/test_L1.0_2/test_L1.0_2/analysis/data/test_L1.0_2.cnf
 GOING TO SKIPT THIS SUBMISSION!

Read out results

[12]:
sf.calculate_solvation_free_energy(n_points=n_points)
/home/mlehner/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:296: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block2_values] [items->Index(['lambda', 'totals', 'baths', 'bonded', 'nonbonded', 'special', 'eds',
       'precalclam'],
      dtype='object')]

  self.database.to_hdf(
/home/mlehner/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:296: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block2_values] [items->Index(['lambda', 'totals', 'baths', 'bonded', 'nonbonded', 'special', 'eds',
       'precalclam'],
      dtype='object')]

  self.database.to_hdf(
/home/mlehner/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:296: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block2_values] [items->Index(['lambda', 'totals', 'baths', 'bonded', 'nonbonded', 'special', 'eds',
       'precalclam'],
      dtype='object')]

  self.database.to_hdf(
/home/mlehner/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:296: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block2_values] [items->Index(['lambda', 'totals', 'baths', 'bonded', 'nonbonded', 'special', 'eds',
       'precalclam'],
      dtype='object')]

  self.database.to_hdf(
   Lambda        avg       rmsd        ee
0    0.00  28.975699   6.600739  0.634101
1    0.25 -17.135306  15.295040  1.130627
2    0.50  32.891360  24.085484  0.814267
3    0.75  55.630164  45.259040  4.219731
4    1.00  33.014152  19.132328  1.016960
Solvation Free Energy: -23.47933332028344 pm -2.056752317486494
/home/mlehner/PyGromosTools/pygromos/files/trajectory/_general_trajectory.py:296: PerformanceWarning:
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block2_values] [items->Index(['lambda', 'totals', 'baths', 'bonded', 'nonbonded', 'special', 'eds',
       'precalclam'],
      dtype='object')]

  self.database.to_hdf(
[12]:
(-23.47933332028344, -2.056752317486494)