forlilab / waterkit

Tool to predict water molecules placement and energy in ligand binding sites

Waterkit

Tool to predict water molecules placement and energy in ligand binding sites

Eberhardt, J., Forli, S. WaterKit: Thermodynamic Profiling of Protein Hydration Sites, 2023, Journal of Chemical Theory and Computation

Prerequisites

You need, at a minimum (requirements):

• Python (>=3.7)
• OpenBabel
• Numpy
• Scipy
• Pandas
• tqdm (progress bar)
• AutoDock Vina and AutoGrid (for generating maps)
• AmberTools (protein preparation and gist calculations)
• OpenMM (minimization)
• ParmED (files conversion)
• gridData (read dx files from GIST)
• Sphinx (documentation)
• Sphinx_rtd_theme (documentation)

Installation

I highly recommand you to install the Anaconda distribution (https://www.continuum.io/downloads) if you want a clean python environnment with nearly all the prerequisites already installed. To install everything properly, you just have to do this:

$ conda create -n waterkit -c conda-forge python=3 mkl numpy scipy pandas \
    openbabel parmed ambertools openmm netcdf4 griddataformats tqdm \
    sphinx sphinx_rtd_theme
$ conda activate waterkit
$ pip install vina

We can now install the WaterKit package

$ git clone https://github.com/forlilab/waterkit
$ cd waterkit
$ pip install -e .

Finally we will need to compile a special version of autogrid. This version supports up to 56 different atom types. You will need it to run WaterKit.

$ cd autodocksuite-4.2.6-src/autogrid
$ autoreconf -i
$ mkdir x86_64Linux2 # for x86_64architecture
$ cd x86_64Linux2
$ ../configure
$ make
$ make install # optional

Quick tutorial

Receptor preparation

Conversion to PDBQT using AmberTools and wk_prepare_receptor.py script. If you added hydrogen atoms with REDUCE, you can keep them by using the --keep_hydrogen argument. The protonation states for histidine residues (HIE, HID, HIP) will then be automatically assigned based the presence of hydrogen atoms on NE2 or ND1. Histidine residues already set HIE, HID or HIP will remain unchanged. Per default, the protonation state of histidine residues is set to HIE.

$ wk_prepare_receptor.py -i protein.pdb -o protein_prepared --pdb --amber_pdbqt

The following protein coordinate files will be generated:protein_prepared.pdb and protein_prepared_amber.pdbqt. The PDBQT file will be used by WaterKit and the PDB file will be used to create the trajectory file at the end.

Sample water molecule positions with WaterKit

Run WaterKit

$ mkdir traj
# Generate 10.000 frames using 16 cpus (tip3p, 300 K, 3 hydration layers)
# X Y Z define the center and SX SY SZ the size (in Angstrom) of the box
# If it has issue locating autogrid, specify the path with the argument --autogrid_exec_path
$ run_waterkit.py -i protein_prepared_amber.pdbqt -c X Y Z -s SX SY SZ -n 10000 -j 16 -o traj
# Create ensemble trajectory
$ wk_make_trajectory.py -r protein_prepared.pdb -w traj -o protein
# Minimize each conformation (100 steps, 2.5 kcal/mol/A**2 restraints on heavy atoms, CUDA)
$ wk_minimize_trajectory.py -p protein_system.prmtop -t protein_system.nc -o protein_min.nc

Run Grid Inhomogeneous Solvation Theory (GIST)

1. Create input file for cpptraj

# gist.inp
parm protein_system.prmtop
trajin protein.nc
# X Y Z define the center and GX GY GZ the size (in gridpoints) of the box
# Example: if SX = SY = SZ = 24 Angstrom and gridspacn = 0.5, then GX = GY = GZ = 48
gist gridspacn 0.5 gridcntr X Y Z griddim GX GY GZ
go
quit

Usually you would choose the same parameters as the AutoGrid maps (npts and gridcenter). Unlike AutoGrid, the default the grid spacing in GIST is 0.5 A, so you will have to choose box dimension accordingly to match the Autogrid maps dimensions. More informations on GIST are available here: https://amber-md.github.io/cpptraj/CPPTRAJ.xhtml#magicparlabel-4672

2. Run GIST

$ cpptraj -i gist.inp

Identification of hydration sites

Hydration sites are identified based on the oxygen density map (gO) in an iterative way, by selecting the voxel with the highest density, then the second highest and so on, whil keeping a minimum distance of 2.5 A between them. Energies for each of those identified hydration sites are then computed by adding all the surrounding voxels and Gaussian weighted by their distances from the hydration site.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
from gridData import Grid
from waterkit.analysis import HydrationSites
from waterkit.analysis import blur_map

gO = Grid("gist-gO.dx")
esw = Grid('gist-Esw-dens.dx')
eww = Grid('gist-Eww-dens.dx')
tst = Grid('gist-dTStrans-dens.dx')
tso = Grid('gist-dTSorient-dens.dx')
dg = (esw + 2 * eww) - (tst + tso)

# Identification of hydration site positions using gO
hs = HydrationSites(gridsize=0.5, water_radius=1.4, min_water_distance=2.5, min_density=2.0)
hydration_sites = hs.find(gO) # can pass "gist-gO.dx" directly also

# Get Gaussian smoothed energy for hydration sites only
dg_energy = hs.hydration_sites_energy(dg, water_radius=1.4)
hs.export_to_pdb("hydration_sites_dG_smoothed.pdb", hydration_sites, dg_energy)

# ... or get the whole Gaussian smoothed map
map_smooth = blur_map(dg, radius=1.4)
map_smooth.export("gist-dG-dens_smoothed.dx")

????

PROFIT!!!

waterkit's Issues

Placing water with predicted .dx files only?

Hi, I would like to use waterkit to perform explicit water placement with predicted .dx files (from other algorithms) without performing any MD simulations.
Can waterkit do that?

Issue creating pdbqt file with ligand bound system

Hi,

I am trying to generate the waterkit receptor using the updated wk_prepare_receptor.py script.

I have generated my Amber lib and frcmod files for my ligand using the following procedure:

antechamber -i ligand.pdb -fi pdb -o ligand.prepi -fo prepi -c bcc -nc 1 -at gaff2
parmchk2 -i ligand.prepi -f prepi -o ligand.frcmod -s gaff2 

Then to generate the lib file I have run:

tleap
loadamberprep ligand.prepi
saveoff UNK ligand.lib 

However, when running the following command the script fails to fun:

wk_prepare_receptor.py -i protein_ligand.pdb --lib ligand.lib --frcmod ligand.frcmod -o protein_prepared_ligand_test --pdb --amber_pdbqt

The error I obtain is as follows:

Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.
INFO:WaterKit receptor preparation:Amber lib parameter files for residue(s): UNK
INFO:WaterKit receptor preparation:Removed all hydrogen atoms
INFO:WaterKit receptor preparation:Removed all water molecules
WARNING:WaterKit receptor preparation:Histidine protonation will be automatically assigned to HIE: HIS - 11, HIS - 31, HIS - 73, HIS - 173, HIS - 183, HIS - 187, HIS - 192
ERROR:WaterKit receptor preparation:Could not generate topology/coordinates files with tleap
Traceback (most recent call last):
  File "/Users/michaelcarter/opt/miniconda3/envs/waterkit_corr/lib/python3.9/site-packages/pdb4amber/utils.py", line 10, in easy_call
    output = subprocess.check_output(
  File "/Users/michaelcarter/opt/miniconda3/envs/waterkit_corr/lib/python3.9/subprocess.py", line 424, in check_output
    return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
  File "/Users/michaelcarter/opt/miniconda3/envs/waterkit_corr/lib/python3.9/subprocess.py", line 528, in run
    raise CalledProcessError(retcode, process.args,
subprocess.CalledProcessError: Command 'tleap -s -f leap.template.in > leap.template.out' returned non-zero exit status 31.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/michaelcarter/DD_tools/waterkit/scripts/wk_prepare_receptor.py", line 761, in prepare
    easy_call('tleap -s -f %s > %s' % (tleap_input, tleap_output), shell=True)
  File "/Users/michaelcarter/opt/miniconda3/envs/waterkit_corr/lib/python3.9/site-packages/pdb4amber/utils.py", line 14, in easy_call
    raise RuntimeError(e.output.decode())
RuntimeError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/michaelcarter/opt/miniconda3/envs/waterkit_corr/bin/wk_prepare_receptor.py", line 7, in <module>
    exec(compile(f.read(), __file__, 'exec'))
  File "/Users/michaelcarter/DD_tools/waterkit/scripts/wk_prepare_receptor.py", line 862, in <module>
    main()
  File "/Users/michaelcarter/DD_tools/waterkit/scripts/wk_prepare_receptor.py", line 846, in main
    pr.prepare(pdb_filename, lib_files, frcmod_files, clean)
  File "/Users/michaelcarter/DD_tools/waterkit/scripts/wk_prepare_receptor.py", line 765, in prepare
    raise RuntimeError(error_msg)
RuntimeError: Could not generate topology/coordinates files with tleap

I think this error stems from the lib file as I can run the command fine without it. Do you have any suggestions on how to get this working?

Receptor preparation

Hello,
I want to ask what is the difference between a pdb file processed by wk_prepare_recepter.py and a pdb file processed directly by the pdb4amber command?

failed steps in receptor preparation; documentation on receptor prep?

tleap has failed to generate the topology/coords for my input pdb:

python /media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py -i MDH-258_molrep7_refmac3-coot-0.pdb -o MDH-258-xtal_prepared --pdb --amber_pdbqt
WARNING:WaterKit receptor preparation:Found residue(s) with missing heavy atoms: [<Residue LYS[32]; chain=A; segid=A--->, 4], [<Residue ASN[34]; chain=A; segid=A--->, 3], [<Residue LYS[38]; chain=A; segid=A--->, 3], [<Residue LYS[46]; chain=A; segid=A--->, 2], [<Residue LYS[78]; chain=A; segid=A--->, 4], [<Residue LYS[87]; chain=A; segid=A--->, 3], [<Residue LYS[120]; chain=A; segid=A--->, 2], [<Residue LYS[159]; chain=A; segid=A--->, 3], [<Residue LYS[162]; chain=A; segid=A--->, 4], [<Residue ARG[181]; chain=A; segid=A--->, 6], [<Residue GLN[204]; chain=A; segid=A--->, 4], [<Residue LYS[207]; chain=A; segid=A--->, 1], [<Residue ASP[211]; chain=A; segid=A--->, 3], [<Residue ASP[212]; chain=A; segid=A--->, 3], [<Residue ASP[213]; chain=A; segid=A--->, 3], [<Residue ASN[214]; chain=A; segid=A--->, 3], [<Residue LYS[227]; chain=A; segid=A--->, 4], [<Residue LYS[239]; chain=A; segid=A--->, 2], [<Residue SER[244]; chain=A; segid=A--->, 1], [<Residue LYS[247]; chain=A; segid=A--->, 4], [<Residue LYS[259]; chain=A; segid=A--->, 3], [<Residue LYS[272]; chain=A; segid=A--->, 2], [<Residue GLN[289]; chain=A; segid=A--->, 3], [<Residue GLN[298]; chain=A; segid=A--->, 3], [<Residue GLU[322]; chain=A; segid=A--->, 3], [<Residue LYS[331]; chain=A; segid=A--->, 4]
INFO:WaterKit receptor preparation:Removed all hydrogen atoms
INFO:WaterKit receptor preparation:Removed all water molecules
INFO:WaterKit receptor preparation:Removed all non-standard Amber residues: LIG
INFO:WaterKit receptor preparation:Histidine protonation states were automatically set to: HIE - 44, HIE - 86, HIE - 128, HIE - 170, HIE - 178, HIE - 255, HIE - 300, HIE - 310
INFO:WaterKit receptor preparation:Lysine protonation states were automatically set to: LYN - 32, LYN - 38, LYN - 46, LYN - 52, LYN - 67, LYN - 70, LYN - 78, LYN - 81, LYN - 87, LYN - 109, LYN - 112, LYN - 120, LYN - 123, LYN - 126, LYN - 159, LYN - 162, LYN - 165, LYN - 207, LYN - 221, LYN - 227, LYN - 239, LYN - 247, LYN - 250, LYN - 256, LYN - 259, LYN - 272, LYN - 291, LYN - 296, LYN - 325, LYN - 328, LYN - 331
INFO:WaterKit receptor preparation:Cysteine protonation states were automatically set to: CYM - 121, CYM - 134, CYM - 135, CYM - 163, CYM - 195, CYM - 205, CYM - 248, CYM - 261, CYM - 309, CYM - 334
ERROR:WaterKit receptor preparation:Could not generate topology/coordinates files with tleap
Traceback (most recent call last):
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/pdb4amber/utils.py", line 10, in easy_call
    output = subprocess.check_output(
             ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/subprocess.py", line 466, in check_output
    return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/subprocess.py", line 571, in run
    raise CalledProcessError(retcode, process.args,
subprocess.CalledProcessError: Command 'tleap -s -f leap.template.in > leap.template.out' returned non-zero exit status 31.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 892, in prepare
    easy_call('tleap -s -f %s > %s' % (tleap_input, tleap_output), shell=True)
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/pdb4amber/utils.py", line 14, in easy_call
    raise RuntimeError(e.output.decode())
RuntimeError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 997, in <module>
    main()
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 981, in main
    pr.prepare(pdb_filename, lib_files, frcmod_files, clean)
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 896, in prepare
    raise RuntimeError(error_msg)
RuntimeError: Could not generate topology/coordinates files with tleap

I was able to get waterkit to work on a test case on a crystal structure of rhodopsin, but this was only after I 1) ran the protein preparation wizard in Maestro and 2) exported the pdb in PyMOL. For some reason parmed isnt detecting the pdb produced by Maestro. Below is the error message I got from waterkit after putting it through Maestro (reproduced on two different protein structures):

 python /media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py -i 3pqr_prepared.pd b -o 3pqr_wk_prepared --pdb --amber_pdbqt
Traceback (most recent call last):
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 997, in <module>
    main()
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 981, in main
    pr.prepare(pdb_filename, lib_files, frcmod_files, clean)
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 728, in prepare
    receptor = pmd.load_file(pdb_filename)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/parmed/formats/registry.py", line 180, in load_file
    raise FormatNotFound('Could not identify file format')
parmed.exceptions.FormatNotFound: Could not identify file format

On my protein that I am trying to directly work on, I put it through Maestro and PyMOL, but got a different error message:

python /media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py -i  mdh-258-pymol.pdb -o mdh-258-xtal_prepared --pdb --amber_pdbqt
WARNING:WaterKit receptor preparation:Found residue(s) with missing heavy atoms: [<Residue LYS[32]; chain=A>, 4], [<Residue PRO[33]; chain=A>, 7], [< Residue ASN[34]; chain=A>, 4], [<Residue TYR[35]; chain=A>, 12], [<Residue ALA[36]; chain=A>, 5], [<Residue LEU[37]; chain=A>, 8], [<Residue LYS[38];  chain=A>, 4], [<Residue PHE[39]; chain=A>, 11], [<Residue THR[40]; chain=A>, 7], [<Residue LEU[41]; chain=A>, 8], [<Residue ALA[42]; chain=A>, 5], [ <Residue GLY[43]; chain=A>, 4], [<Residue HIS[44]; chain=A>, 10], [<Residue THR[45]; chain=A>, 7], [<Residue LYS[46]; chain=A>, 4], [<Residue ALA[47] ; chain=A>, 5], [<Residue VAL[48]; chain=A>, 7], [<Residue SER[49]; chain=A>, 6], [<Residue SER[50]; chain=A>, 6], [<Residue VAL[51]; chain=A>, 7], [ <Residue LYS[52]; chain=A>, 9], [<Residue PHE[53]; chain=A>, 11], [<Residue SER[54]; chain=A>, 6], [<Residue PRO[55]; chain=A>, 7], [<Residue ASN[56] ; chain=A>, 8], [<Residue GLY[57]; chain=A>, 4], [<Residue GLU[58]; chain=A>, 9], [<Residue TRP[59]; chain=A>, 14], [<Residue LEU[60]; chain=A>, 8],  [<Residue ALA[61]; chain=A>, 5], [<Residue SER[62]; chain=A>, 6], [<Residue SER[63]; chain=A>, 6], [<Residue SER[64]; chain=A>, 6], [<Residue ALA[65] ; chain=A>, 5], [<Residue ASP[66]; chain=A>, 8], [<Residue LYS[67]; chain=A>, 9], [<Residue LEU[68]; chain=A>, 8], [<Residue ILE[69]; chain=A>, 8], [ <Residue LYS[70]; chain=A>, 9], [<Residue ILE[71]; chain=A>, 8], [<Residue TRP[72]; chain=A>, 14], [<Residue GLY[73]; chain=A>, 4], [<Residue ALA[74] ; chain=A>, 5], [<Residue TYR[75]; chain=A>, 12], [<Residue ASP[76]; chain=A>, 8], [<Residue GLY[77]; chain=A>, 4], [<Residue LYS[78]; chain=A>, 4],  [<Residue PHE[79]; chain=A>, 11], [<Residue GLU[80]; chain=A>, 9], [<Residue LYS[81]; chain=A>, 9], [<Residue THR[82]; chain=A>, 7], [<Residue ILE[83 ]; chain=A>, 8], [<Residue SER[84]; chain=A>, 6], [<Residue GLY[85]; chain=A>, 4], [<Residue HIS[86]; chain=A>, 10], [<Residue LYS[87]; chain=A>, 4],  [<Residue LEU[88]; chain=A>, 8], [<Residue GLY[89]; chain=A>, 4], [<Residue ILE[90]; chain=A>, 8], [<Residue SER[91]; chain=A>, 6], [<Residue ASP[92 ]; chain=A>, 8], [<Residue VAL[93]; chain=A>, 7], [<Residue ALA[94]; chain=A>, 5], [<Residue TRP[95]; chain=A>, 14], [<Residue SER[96]; chain=A>, 6],  [<Residue SER[97]; chain=A>, 6], [<Residue ASP[98]; chain=A>, 8], [<Residue SER[99]; chain=A>, 6], [<Residue ASN[100]; chain=A>, 8], [<Residue LEU[1 01]; chain=A>, 8], [<Residue LEU[102]; chain=A>, 8], [<Residue VAL[103]; chain=A>, 7], [<Residue SER[104]; chain=A>, 6], [<Residue ALA[105]; chain=A> , 5], [<Residue SER[106]; chain=A>, 6], [<Residue ASP[107]; chain=A>, 8], [<Residue ASP[108]; chain=A>, 8], [<Residue LYS[109]; chain=A>, 9], [<Resid ue THR[110]; chain=A>, 7], [<Residue LEU[111]; chain=A>, 8], [<Residue LYS[112]; chain=A>, 9], [<Residue ILE[113]; chain=A>, 8], [<Residue TRP[114];  chain=A>, 14], [<Residue ASP[115]; chain=A>, 8], [<Residue VAL[116]; chain=A>, 7], [<Residue SER[117]; chain=A>, 6], [<Residue SER[118]; chain=A>, 6] , [<Residue GLY[119]; chain=A>, 4], [<Residue LYS[120]; chain=A>, 4], [<Residue CYS[121]; chain=A>, 6], [<Residue LEU[122]; chain=A>, 8], [<Residue L YS[123]; chain=A>, 9], [<Residue THR[124]; chain=A>, 7], [<Residue LEU[125]; chain=A>, 8], [<Residue LYS[126]; chain=A>, 9], [<Residue GLY[127]; chai n=A>, 4], [<Residue HIS[128]; chain=A>, 10], [<Residue SER[129]; chain=A>, 6], [<Residue ASN[130]; chain=A>, 8], [<Residue TYR[131]; chain=A>, 12], [ <Residue VAL[132]; chain=A>, 7], [<Residue PHE[133]; chain=A>, 11], [<Residue CYS[134]; chain=A>, 6], [<Residue CYS[135]; chain=A>, 6], [<Residue ASN [136]; chain=A>, 8], [<Residue PHE[137]; chain=A>, 11], [<Residue ASN[138]; chain=A>, 8], [<Residue PRO[139]; chain=A>, 7], [<Residue GLN[140]; chain =A>, 9], [<Residue SER[141]; chain=A>, 6], [<Residue ASN[142]; chain=A>, 8], [<Residue LEU[143]; chain=A>, 8], [<Residue ILE[144]; chain=A>, 8], [<Re sidue VAL[145]; chain=A>, 7], [<Residue SER[146]; chain=A>, 6], [<Residue GLY[147]; chain=A>, 4], [<Residue SER[148]; chain=A>, 6], [<Residue PHE[149 ]; chain=A>, 11], [<Residue ASP[150]; chain=A>, 8], [<Residue GLU[151]; chain=A>, 9], [<Residue SER[152]; chain=A>, 6], [<Residue VAL[153]; chain=A>,  7], [<Residue ARG[154]; chain=A>, 11], [<Residue ILE[155]; chain=A>, 8], [<Residue TRP[156]; chain=A>, 14], [<Residue ASP[157]; chain=A>, 8], [<Resi due VAL[158]; chain=A>, 7], [<Residue LYS[159]; chain=A>, 4], [<Residue THR[160]; chain=A>, 7], [<Residue GLY[161]; chain=A>, 4], [<Residue LYS[162];  chain=A>, 4], [<Residue CYS[163]; chain=A>, 6], [<Residue LEU[164]; chain=A>, 8], [<Residue LYS[165]; chain=A>, 9], [<Residue THR[166]; chain=A>, 7] , [<Residue LEU[167]; chain=A>, 8], [<Residue PRO[168]; chain=A>, 7], [<Residue ALA[169]; chain=A>, 5], [<Residue HIS[170]; chain=A>, 10], [<Residue  SER[171]; chain=A>, 6], [<Residue ASP[172]; chain=A>, 8], [<Residue PRO[173]; chain=A>, 7], [<Residue VAL[174]; chain=A>, 7], [<Residue SER[175]; cha in=A>, 6], [<Residue ALA[176]; chain=A>, 5], [<Residue VAL[177]; chain=A>, 7], [<Residue HIS[178]; chain=A>, 10], [<Residue PHE[179]; chain=A>, 11],  [<Residue ASN[180]; chain=A>, 8], [<Residue ARG[181]; chain=A>, 4], [<Residue ASP[182]; chain=A>, 8], [<Residue GLY[183]; chain=A>, 4], [<Residue SER [184]; chain=A>, 6], [<Residue LEU[185]; chain=A>, 8], [<Residue ILE[186]; chain=A>, 8], [<Residue VAL[187]; chain=A>, 7], [<Residue SER[188]; chain= A>, 6], [<Residue SER[189]; chain=A>, 6], [<Residue SER[190]; chain=A>, 6], [<Residue TYR[191]; chain=A>, 12], [<Residue ASP[192]; chain=A>, 8], [<Re sidue GLY[193]; chain=A>, 4], [<Residue LEU[194]; chain=A>, 8], [<Residue CYS[195]; chain=A>, 6], [<Residue ARG[196]; chain=A>, 11], [<Residue ILE[19 7]; chain=A>, 8], [<Residue TRP[198]; chain=A>, 14], [<Residue ASP[199]; chain=A>, 8], [<Residue THR[200]; chain=A>, 7], [<Residue ALA[201]; chain=A> , 5], [<Residue SER[202]; chain=A>, 6], [<Residue GLY[203]; chain=A>, 4], [<Residue GLN[204]; chain=A>, 4], [<Residue CYS[205]; chain=A>, 6], [<Resid ue LEU[206]; chain=A>, 8], [<Residue LYS[207]; chain=A>, 4], [<Residue THR[208]; chain=A>, 7], [<Residue LEU[209]; chain=A>, 8], [<Residue ILE[210];  chain=A>, 8], [<Residue ASP[211]; chain=A>, 4], [<Residue ASP[212]; chain=A>, 4], [<Residue ASP[213]; chain=A>, 4], [<Residue ASN[214]; chain=A>, 4],  [<Residue PRO[215]; chain=A>, 7], [<Residue PRO[216]; chain=A>, 7], [<Residue VAL[217]; chain=A>, 7], [<Residue SER[218]; chain=A>, 6], [<Residue PH E[219]; chain=A>, 11], [<Residue VAL[220]; chain=A>, 7], [<Residue LYS[221]; chain=A>, 9], [<Residue PHE[222]; chain=A>, 11], [<Residue SER[223]; cha in=A>, 6], [<Residue PRO[224]; chain=A>, 7], [<Residue ASN[225]; chain=A>, 8], [<Residue GLY[226]; chain=A>, 4], [<Residue LYS[227]; chain=A>, 4], [< Residue TYR[228]; chain=A>, 12], [<Residue ILE[229]; chain=A>, 8], [<Residue LEU[230]; chain=A>, 8], [<Residue ALA[231]; chain=A>, 5], [<Residue ALA[ 232]; chain=A>, 5], [<Residue THR[233]; chain=A>, 7], [<Residue LEU[234]; chain=A>, 8], [<Residue ASP[235]; chain=A>, 8], [<Residue ASN[236]; chain=A >, 8], [<Residue THR[237]; chain=A>, 7], [<Residue LEU[238]; chain=A>, 8], [<Residue LYS[239]; chain=A>, 4], [<Residue LEU[240]; chain=A>, 8], [<Resi due TRP[241]; chain=A>, 14], [<Residue ASP[242]; chain=A>, 8], [<Residue TYR[243]; chain=A>, 12], [<Residue SER[244]; chain=A>, 4], [<Residue SER[245 ]; chain=A>, 6], [<Residue GLY[246]; chain=A>, 4], [<Residue LYS[247]; chain=A>, 4], [<Residue CYS[248]; chain=A>, 6], [<Residue LEU[249]; chain=A>,  8], [<Residue LYS[250]; chain=A>, 9], [<Residue THR[251]; chain=A>, 7], [<Residue TYR[252]; chain=A>, 12], [<Residue THR[253]; chain=A>, 7], [<Residu e GLY[254]; chain=A>, 4], [<Residue HIS[255]; chain=A>, 10], [<Residue LYS[256]; chain=A>, 9], [<Residue ASN[257]; chain=A>, 8], [<Residue GLU[258];  chain=A>, 9], [<Residue LYS[259]; chain=A>, 4], [<Residue TYR[260]; chain=A>, 12], [<Residue CYS[261]; chain=A>, 6], [<Residue ILE[262]; chain=A>, 8] , [<Residue PHE[263]; chain=A>, 11], [<Residue ALA[264]; chain=A>, 5], [<Residue ASN[265]; chain=A>, 8], [<Residue PHE[266]; chain=A>, 11], [<Residue  SER[267]; chain=A>, 6], [<Residue VAL[268]; chain=A>, 7], [<Residue THR[269]; chain=A>, 7], [<Residue GLY[270]; chain=A>, 4], [<Residue GLY[271]; ch ain=A>, 4], [<Residue LYS[272]; chain=A>, 4], [<Residue TRP[273]; chain=A>, 14], [<Residue ILE[274]; chain=A>, 8], [<Residue VAL[275]; chain=A>, 7],  [<Residue SER[276]; chain=A>, 6], [<Residue GLY[277]; chain=A>, 4], [<Residue SER[278]; chain=A>, 6], [<Residue GLU[279]; chain=A>, 9], [<Residue ASP [280]; chain=A>, 8], [<Residue ASN[281]; chain=A>, 8], [<Residue LEU[282]; chain=A>, 8], [<Residue VAL[283]; chain=A>, 7], [<Residue TYR[284]; chain= A>, 12], [<Residue ILE[285]; chain=A>, 8], [<Residue TRP[286]; chain=A>, 14], [<Residue ASN[287]; chain=A>, 8], [<Residue LEU[288]; chain=A>, 8], [<R esidue GLN[289]; chain=A>, 4], [<Residue THR[290]; chain=A>, 7], [<Residue LYS[291]; chain=A>, 9], [<Residue GLU[292]; chain=A>, 9], [<Residue ILE[29 3]; chain=A>, 8], [<Residue VAL[294]; chain=A>, 7], [<Residue GLN[295]; chain=A>, 9], [<Residue LYS[296]; chain=A>, 9], [<Residue LEU[297]; chain=A>,  8], [<Residue GLN[298]; chain=A>, 4], [<Residue GLY[299]; chain=A>, 4], [<Residue HIS[300]; chain=A>, 10], [<Residue THR[301]; chain=A>, 7], [<Resid ue ASP[302]; chain=A>, 8], [<Residue VAL[303]; chain=A>, 7], [<Residue VAL[304]; chain=A>, 7], [<Residue ILE[305]; chain=A>, 8], [<Residue SER[306];  chain=A>, 6], [<Residue THR[307]; chain=A>, 7], [<Residue ALA[308]; chain=A>, 5], [<Residue CYS[309]; chain=A>, 6], [<Residue HIS[310]; chain=A>, 10] , [<Residue PRO[311]; chain=A>, 7], [<Residue THR[312]; chain=A>, 7], [<Residue GLU[313]; chain=A>, 9], [<Residue ASN[314]; chain=A>, 8], [<Residue I LE[315]; chain=A>, 8], [<Residue ILE[316]; chain=A>, 8], [<Residue ALA[317]; chain=A>, 5], [<Residue SER[318]; chain=A>, 6], [<Residue ALA[319]; chai n=A>, 5], [<Residue ALA[320]; chain=A>, 5], [<Residue LEU[321]; chain=A>, 8], [<Residue GLU[322]; chain=A>, 4], [<Residue ASN[323]; chain=A>, 8], [<R esidue ASP[324]; chain=A>, 8], [<Residue LYS[325]; chain=A>, 9], [<Residue THR[326]; chain=A>, 7], [<Residue ILE[327]; chain=A>, 8], [<Residue LYS[32 8]; chain=A>, 9], [<Residue LEU[329]; chain=A>, 8], [<Residue TRP[330]; chain=A>, 14], [<Residue LYS[331]; chain=A>, 4], [<Residue SER[332]; chain=A> , 6], [<Residue ASP[333]; chain=A>, 8], [<Residue CYS[334]; chain=A>, 6], [<Residue LYS[32]; chain=A; segid=A--->, 5], [<Residue ASN[34]; chain=A; se gid=A--->, 4], [<Residue LYS[38]; chain=A; segid=A--->, 5], [<Residue LYS[46]; chain=A; segid=A--->, 5], [<Residue LYS[78]; chain=A; segid=A--->, 5],  [<Residue LYS[87]; chain=A; segid=A--->, 5], [<Residue LYS[120]; chain=A; segid=A--->, 5], [<Residue LYS[159]; chain=A; segid=A--->, 5], [<Residue L YS[162]; chain=A; segid=A--->, 5], [<Residue ARG[181]; chain=A; segid=A--->, 7], [<Residue GLN[204]; chain=A; segid=A--->, 5], [<Residue LYS[207]; ch ain=A; segid=A--->, 5], [<Residue ASP[211]; chain=A; segid=A--->, 4], [<Residue ASP[212]; chain=A; segid=A--->, 4], [<Residue ASP[213]; chain=A; segi d=A--->, 4], [<Residue ASN[214]; chain=A; segid=A--->, 4], [<Residue LYS[227]; chain=A; segid=A--->, 5], [<Residue LYS[239]; chain=A; segid=A--->, 5] , [<Residue SER[244]; chain=A; segid=A--->, 2], [<Residue LYS[247]; chain=A; segid=A--->, 5], [<Residue LYS[259]; chain=A; segid=A--->, 5], [<Residue  LYS[272]; chain=A; segid=A--->, 5], [<Residue GLN[289]; chain=A; segid=A--->, 5], [<Residue GLN[298]; chain=A; segid=A--->, 5], [<Residue GLU[322];  chain=A; segid=A--->, 5], [<Residue LYS[331]; chain=A; segid=A--->, 5]
INFO:WaterKit receptor preparation:Removed all hydrogen atoms
INFO:WaterKit receptor preparation:Removed all water molecules
INFO:WaterKit receptor preparation:Removed all non-standard Amber residues: LIG, NMA
Traceback (most recent call last):
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 494, in _find_gaps
    n_atom = [atom for atom in next_residue.atoms if atom.name == 'N'][0]
             ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 997, in <module>
    main()
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 981, in main
    pr.prepare(pdb_filename, lib_files, frcmod_files, clean)
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 777, in prepare
    gaplist = _find_gaps(pdbfixer.parm, RESPROT, self._fill_gaps)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 496, in _find_gaps
    gaprecord = (9999.999, c_atom.residue.name, i + 1, n_atom.residue.name, i + 2)
                                                       ^^^^^^
UnboundLocalError: cannot access local variable 'n_atom' where it is not associated with a value

It is not clear to me why waterkit/ambertools failed to properly assign n_atom? Thanks for the help with this.

How to calculate the dG, dH, dS and density all the hydration sites?

Dear developers,

How to calculate the dG, dH, dS and density all the hydration sites?
With the script on Waterkit homepage, only dG is calculated.
https://github.com/forlilab/waterkit#identification-of-hydration-sites

Best
Zhenting

ad_map attribute

Hi, the program crashes on line 95 in run_waterkit.py because ad_map has no add_map attribute and later on in utils.py no info attribute. (I wondered if it meant to use sw_map since this does have the attributes, but then there are issues with electrostatics missing for tip3p water model?!)

cutoff of desolvation free energy for a favorable hydration site?

Dear developers,

Inside the Waterkit paper, a cutoff of -5 kcal/mol is used to select favorable Waterkit hydration site for comparison with GIST result.
Is -5 kcal/mol a recommended cutoff to define a favorable hydration site, any data to support this?

What is the range of the desolvation free energy(DFE)? For the paper and my tests, DFE is always below 0, is it correct?

Best
Zhenting

run_waterkit.py failed

We tried to use waterkit to predict the distribution of water molecules near the binding site of protein-small molecule complex, but if the small molecule contains Cl atoms, when run the step run_waterkit.py, an error occurred.

run_waterkit.py -i com_prepared_amber.pdbqt -c 26.95 6.16 4.15 -s 23 23 23 -n 10000 -j 16 -o traj

PDBQT parsing error: Atom type CL is not a valid AutoDock type (atom types are case-sensitive).

ATOM 2610 Cl1 BNZ A 163 29.032 7.968 2.789 1.00 1.00 -0.097 CL

in our input amber_pdbqt file:
ATOM 2611 Cl1 BNZ A 163 29.032 7.968 2.789 1.00 1.00 -0.097 cl

com_prepared_amber.zip

use in docking or virtual screening

Hi,
Thanks for sharing this.
This is not an issue, but more like requesting for more help.
Is there a quick tutorial or example where it is illustrated how this is used in docking/virtual screening using either AutoDock Vina/GPU?
Thanks for your help.

Working directly with parm7 and coordinate files

Hello.

Thanks for this awesome work.

I have a protein-ligand system parameterized with openff. I can save it as amber format parm7 or prmtop files. But, it's not very straightforward to save frcmod and lib files https://github.com/openforcefield/openff-toolkit/issues/304

Do you think there could be a way to work directly with parm7/prmtop files?

I would be really grateful for any suggestions.

Best,
Amin.

Installation issues

Dear developer

cd autogrid
autoreconf -i
mkdir x86_64Linux2 # for x86_64architecture
cd x86_64Linux2
../configure
make

After the above commands, compilation of autogrid stops with error
make: *** No rule to make target '../../autodock/paramdat2h.csh', needed by 'default_parameters.h'. Stop.

Could you help on this?

Zhenting

wk_minimize_trajectory.py fails: kcal not accepted unit by openmm

The minimize step for running waterkit fails because openmm will not read in units of energy in kcal/mol:

python /media/bak11/binaries/git/waterkit/scripts/wk_minimize_trajectory.py -p protein_system.prmtop -t protein_system.nc -o protein_min.nc
Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.
Traceback (most recent call last):
  File "/media/bak11/binaries/git/waterkit/scripts/wk_minimize_trajectory.py", line 145, in <module>
    main()
  File "/media/bak11/binaries/git/waterkit/scripts/wk_minimize_trajectory.py", line 141, in main
    m.minimize_trajectory(prmtop_filename, traj_filename, output_filename)
  File "/media/bak11/binaries/git/waterkit/scripts/wk_minimize_trajectory.py", line 110, in minimize_trajectory
    simulation.minimizeEnergy(maxIterations=self._n_steps, tolerance=tolerance)
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/openmm/app/simulation.py", line 143, in minimizeEnergy
    mm.LocalEnergyMinimizer.minimize(self.context, tolerance, maxIterations, reporter)
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/openmm/openmm.py", line 4422, in minimize
    tolerance = tolerance.value_in_unit(unit.kilojoules_per_mole/unit.nanometer)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/openmm/unit/quantity.py", line 626, in value_in_unit
    val = self.in_units_of(unit)
          ^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/openmm/unit/quantity.py", line 662, in in_units_of
    raise TypeError('Unit "%s" is not compatible with Unit "%s".' % (self.unit, other_unit))
TypeError: Unit "kilocalorie/mole" is not compatible with Unit "kilojoule/(nanometer*mole)".

Looks like the culprit is in the import calls at the beginning of wk_miniimize_trajectory.py (simtk is deprecated in openmm 7/8 btw):

     17 from simtk.unit import Quantity, picoseconds, kilocalories_per_mole, angstroms, kelvin
     18 from simtk.openmm import CustomExternalForce, LangevinIntegrator, Platform
     19 from simtk.openmm.app import AmberPrmtopFile, HBonds, CutoffNonPeriodic, Simulation

and here is the code snippet from wk_minimize_trajectory.py that is causing it to fail:

     63     def minimize_trajectory(self, prmtop_filename, traj_filename, output_filename):
     64         nonbondedMethod = CutoffNonPeriodic
     65         nonbondedCutoff = 9 * angstroms
     66         rigidWater = True
     67         constraints = HBonds
     68         dt = 0.002 * picoseconds
     69         temperature = 300 * kelvin
     70         friction = 1.0 / picoseconds
     71         tolerance = 1 * kilocalories_per_mole
     72         K = self._restraint * kilocalories_per_mole / angstroms**2

Should be a simple enough fix -- import kilojoules per mol instead of kcal (http://docs.openmm.org/latest/userguide/theory/01_introduction.html) and change the minimize_trajectory function for tolerance and K to use kJ/mol instead of kcal/mol. However, are there any other steps in the waterkit workflow that would be affected that utilize kcal/mol? wk_minimize_trajectory.py is the only script that explicitly references kcal/mol.

EDIT: on second thought, it must be something in how tolerance is being defined as a variable, since I don't get an error in python when I import unit from simtk.openmm. Is it something similar to this issue? https://github.com/openmm/openmm/issues/4069

Receptor preparation falied

I have run constant pH simulations with Gromacs and Now I want to trace the water molecules at different pH in my trajectory that processed in PDB format with last 1000 frames.

When I run receptor preparation command, I get the following error

**raise PDBError(f'Coordinate mismatch in model {self._current_model_number}')

parmed.exceptions.PDBError: Coordinate mismatch in model 1**

It seems that issue should be with Parmed
and I have also added chain identifiers but still not working

Is Waterkit suitable for ligand bound system?

Hi,

In the paper about Waterkit, at the "Structure preparation" section, it reads
"Crystallographic water molecules and Amber nonstandard residues were removed from the receptor, while ions were kept."

• Does that mean that ligand is also removed before Waterkit calculation?
• Why ligand is removed?
• Does Waterkit parameters work well with ligand (organic compound)?

Best
Zhenting

python 3.12 breaks waterkit

Just an FYI to waterkit users -- I just installed waterkit using the installation instructions, but my conda environment pulled python 3.12.1. 3.12 dropped support for the imp module and so until that issue is resolved within your codebase, you need to restrict the python version to be <= 3.11: https://stackoverflow.com/questions/77274572/multiqc-modulenotfounderror-no-module-named-imp

I was getting the following error when trying to run wk_prepare_receptor.py:

python /media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py -i MDH-258_molrep7_refmac3-coot-0.pdb -o MDH-258-xtal_prepared --pdb --amber_pdbqt
Traceback (most recent call last):
  File "/media/bak11/binaries/git/waterkit/scripts/wk_prepare_receptor.py", line 20, in <module>
    from waterkit import utils
  File "/media/bak11/binaries/git/waterkit/waterkit/__init__.py", line 8, in <module>
    from .autogrid import AutoGrid
  File "/media/bak11/binaries/git/waterkit/waterkit/autogrid.py", line 16, in <module>
    from .molecule import Molecule
  File "/media/bak11/binaries/git/waterkit/waterkit/molecule.py", line 9, in <module>
    import imp
ModuleNotFoundError: No module named 'imp'

Specifying python=3.11 when creating the environment fixed the issue.

documentation on output files from waterkit

After getting waterkit to successfully complete on a sample protein, the following files are produced:

gist-Esw-dens.dx
gist-Eww-dens.dx
gist-dTStrans-dens.dx
gist-dTSorient-dens.dx
gist-dTSsix-dens.dx
gist-neighbor-norm.dx
gist-dipole-dens.dx
gist-order-norm.dx
gist-dipolex-dens.dx
gist-dipoley-dens.dx
gist-dipolez-dens.dx
gist-gO.dx
gist-output.dat
gist-gH.dx
hydration.py
hydration_sites_dG_smoothed.pdb
gist-dG-dens_smoothed.dx

I know how to load and visualize the pdb/dx files, but I was wondering if there is any documentation describing each of the dx files, or if there is none, could it be provided?

wk_prepare_receptor.py failed

python ../scripts/wk_prepare_receptor.py -i 7u9y_pro.pdb -o 7u9y_pro_prep --pdb --amber_pdbqt
Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.
INFO:WaterKit receptor preparation:Removed all hydrogen atoms
INFO:WaterKit receptor preparation:Removed all water molecules
Traceback (most recent call last):
File "/tools/solvation/waterkit-master/test/../scripts/wk_prepare_receptor.py", line 854, in
main()
File "/tools/solvation/waterkit-master/test/../scripts/wk_prepare_receptor.py", line 838, in main
pr.prepare(pdb_filename, lib_files, frcmod_files, clean)
File "/tools/solvation/waterkit-master/test/../scripts/wk_prepare_receptor.py", line 655, in prepare
ns_names = _find_non_standard_resnames(pdbfixer.parm, AMBER_SUPPORTED_RESNAMES + nonstandard_resnames)
TypeError: can only concatenate tuple (not "list") to tuple

7u9y_pro.zip

numpy bool call throws exception for identification of hydration sites

I went to use the provided python script at the end of the example workflow and it threw an exception. Code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
from gridData import Grid
from waterkit.analysis import HydrationSites
from waterkit.analysis import blur_map

gO = Grid("gist-gO.dx")
esw = Grid('gist-Esw-dens.dx')
eww = Grid('gist-Eww-dens.dx')
tst = Grid('gist-dTStrans-dens.dx')
tso = Grid('gist-dTSorient-dens.dx')
dg = (esw + 2 * eww) - (tst + tso)

# Identification of hydration site positions using gO
hs = HydrationSites(gridsize=0.5, water_radius=1.4, min_water_distance=2.5, min_density=2.0)
hydration_sites = hs.find(gO) # can pass "gist-gO.dx" directly also

# Get Gaussian smoothed energy for hydration sites only
dg_energy = hs.hydration_sites_energy(dg, water_radius=1.4)
hs.export_to_pdb("hydration_sites_dG_smoothed.pdb", hydration_sites, dg_energy)

# ... or get the whole Gaussian smoothed map
map_smooth = blur_map(dg, radius=1.4)
map_smooth.export("gist-dG-dens_smoothed.dx")

Error when trying to run:

python hydration.py 
Traceback (most recent call last):
  File "/media/bak11/modulus/MDH/waterkit/mdh-258-xtal/hydration.py", line 17, in <module>
    hydration_sites = hs.find(gO) # can pass "gist-gO.dx" directly also
                      ^^^^^^^^^^^
  File "/media/bak11/binaries/git/waterkit/waterkit/analysis/hydration_sites.py", line 117, in find
    hydration_sites, self._isocontour, self._cluster_ids = _hydration_sites(xyz, values, self._min_density,
                                                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/git/waterkit/waterkit/analysis/hydration_sites.py", line 26, in _hydration_sites
    mask = np.ones(tmp_values.shape, dtype=np.bool)
                                           ^^^^^^^
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/numpy/__init__.py", line 324, in __getattr__
    raise AttributeError(__former_attrs__[attr])
AttributeError: module 'numpy' has no attribute 'bool'.
`np.bool` was a deprecated alias for the builtin `bool`. To avoid this error in existing code, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
    https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations. Did you mean: 'bool_'?

I changed np.bool to bool in line 26 of hydration_sites.py and it completed successfully.