# Implicit Ligand Sampling run script
# ===================================

# Running ILS calculation:

# vmd -dispdev text -e <name_of_this_file>

# You will need VMD 1.8.7. or higher.

# Change the input parameters below to your liking.
# The filenames used in this script are relative to the directory
# for which it was generated but you can of course change then.

# If you have a CUDA enables Nvidia GPU VMD will use the GPU for
# the computation. Since the all the GPU resources will then be
# used for ILS your graphical display will freeze up, so don't be
# surprised. After finishing each frame the display will briefly
# be updated and freeze again.

# Comment this line out to prevent the use of the CUDA implementation:
#set env(VMDCUDAILS) 1

# You might want to do a quick test with 1 frame first to see if
# the syntax is correct and to determine the approximate runtime
# per frame.

package require ilstools

# Adjustable parameters:
# ----------------------

# First and last frames to process
set first 0
set last  5001

# Resolution of final downsampled map in Angstrom
set res    1.0

# Subsampling of each dimension during computation
# i.e. each gridpoint of the final map will actually
# be downsampled from subres^3 points.
set subres 3

# Control of the angular spacing of probe orientation vectors,
# i.e. the number of probe conformers generated.
#
#   1: use 1 orientation only
#   2: use 6 orientations (vertices of octahedron)
#   3: use 8 orientations (vertices of hexahedron)
#   4: use 12 orientations (faces of dodecahedron)
#   5: use 20 orientations (vertices of dodecahedron)
#   6: use 32 orientations (faces+vert. of dodecahedron)
#  >6: geodesic subdivisions of icosahedral faces
#      with frequency 1, 2, ...
#
#   For each orientation a number of rotamers will be
#   generated. The angular spacing of the rotations
#   around the orientation vectors is chosen to be about
#   the same as the angular spacing of the orientation
#   vector itself.
#   If the probe ha at least one symmetry axis then the
#   rotations around the orientation vectors are reduced
#   accordingly. If there is an infinite oder axis (linear
#   molecule) the rotation will be omitted.
#   In case there is an additional perpendicular C2 axis
#   the half of the orientations will be ignored so that
#   there are no antiparallel pairs.
#
#   Probes with tetrahedral symmetry:
#   Here conf denotes the number of rotamers for each of
#   the 8 orientations defined by the vertices of the
#   tetrahedron and its dual tetrahedron.
set orient   5

# Cutoff energy above which the occupancy is regarded zero
# For GPUs energies of more than 87 always correspond to
# floating point values of zero for the occupancy. Hence
# there is no point going higher than that.
set maxen  85

# Temperature of the MD simualtion
set T   300

# Nonbonded interaction cutoff
set cutoff 12.0

# The minmax box defining the free energy map
# (two opposite corners of the grid)
set minmax {{-8 0 -12} {36.0 40.0 27.0}}

# The DX file containing the free energy map
set dxfile apo1a6gils_out.dx

# -------------------------------------------------------

# Set up the probe

# WARNING: The probe parameters have only been verified to
#          reproduce experimental solvation free energies for
#          xenon, oxygen, nitric oxide and carbon monoxide probes.
#          Use parameters for other probes as a starting point
#          for optimization.
set pmol [mol new ilsprobe_oxygen.xyz]
set psel [atomselect $pmol all]
$psel set radius     1.7;  # VDW radius
$psel set occupancy -0.12; # VDW epsilon

# -------------------------------------------------------


# Load the molecule
set molid [mol new /home/newhoir/apo1a6g/apo.psf]
mol addfile "/home/newhoir/apo1a6g/apo.pdb" waitfor all
mol addfile "/home/newhoir/apo1a6g/apo1a6g-last5.dcd" waitfor all

# Selection used for alignment
set asel [atomselect $molid "protein and name CA" frame 0]

# Align all frames of the protein, but only the transformation
# part since a rotated system cannot be wrapped properly.
ILStools::shift_to_center $asel

# Now we have a system were the protein (or whatever you used
# as basis of the alignment) does not shift anymore.
# It still rotates but we have to live with that.
# We rewrap the solvent so that our PBC box is always centered
# at the COM of the protein.
package require pbctools
pbc wrap -molid $molid -orthorhombic -compound residue -all \
         -center com -centersel [$asel text]

# Set the radius and occupancy field for each atom to the
# VDW rmin and epsilon values from the force field
package require ilstools
ILStools::readcharmmparams ../par_all27_4impligsamp.inp
ILStools::assigncharmmparams $molid

# Run ILS
volmap ils $molid $minmax -cutoff $cutoff -pbc \
    -res $res -subres $subres -probesel $psel -orient $orient \
    -alignsel $asel  -maxenergy $maxen \
    -T $T -first $first -last $last \
    -o $dxfile

# Quit VMD when done with ILS calculation
quit