TINKER-OPENMM:Tinker-tut

From TINKER-OPENMM
Jump to: navigation, search

Download & Install

Download from github: https://github.com/jayponder/tinker Tips about compiling Tinker.

http://biomol.bme.utexas.edu/tinker-openmm/index.php/Poltype:Poltype#Quik_notes_about_compiling_tinker_7_or_8.C2.A0with_openMP:.C2.A0 

Input files

You can find examples of input files in tinker distribution under tinker/bench, tinker/example or tinker/test. Two files are required to run TINKER calculations: *.xyz and *.key For example, open butane.xyz and butane.key in tinker/example/ (available in all official distributions) to see what’s inside.

*.xyz file

  • The first number on the 1st line is how many atoms total.
  • There may be a second line that specify the box dimensions if the system is periodic (newer format since tinker 6)
  • The first column is the atomic index
  • The second column is the atomic symbol
  • The 3 -5th columns are x,y,z coordinates in Angstrom
  • The 6th is the “atom type” defined in the *.key file. This is the index tinker uses to assign parameters from the key/parameter file.
  • The 7th – last columns are lists of atoms that are connected to the current atom

*.key file

The key file may have all the actual parameters or a link to the actual parameters file specified in the first line. The parameters specify the bond, angle, torsion, vdW and electrostatic interactions between atoms based on the “atom type”. If you see an error related to OMP, please set the OPENMP-THREADS in the key file to a number less than the # of CPU cores on your computer. This sets how many CPU cores are used in the parallel execution.

For example, in protein.key below, borrowed from tinker/bench/bench7.key, the first line specific the actual parameters are contain in the amoebapro13.prm. More examples can be found in tinker/bench, tinker/example, or tinker/test.

 

parameters           $TINKERDIR/params/amoebapro13.prm
# 
# the amoebapro13.prm is the AMOEBA protein force field in tinker/params/. If you have a ligand you can add the parameters below. More info below
#
# verbose                                  # printing info for every step, for debugging mostly
#
randomseed            123456789
integrator            respa                #this multi-time step integrator allows TINKER to use 2 fs time step
HEAVY-HYDROGEN                             This will increase H atom mass automatically so 3 or 3.5 fs time step can be used. Use at your own risk
neighbor-list                              # this below requires your box is twice the cutoff plus 2-3 Ang.
                                           If your box is too small for vdw cutoff but OK for Ewald, you can use "mpole-list" here.
openmp-threads    16                       # how many core you want to use on the node.
#
#  Define the Periodic Box and Cutoffs
#
a-axis                62.23
vdw-cutoff            12.0
vdw-correction
#
#  Set Parameters for Ewald Summation
#
ewald
ewald-cutoff          7.0                  # we can use such small cutoff because dipole/quadrupole die off faster than point charge.
pme-grid  64 64 64                         #the grid size should be slight bigger than box size  62.23 (1.2x is the best). 
                                           Must be even with factors of only 2, 3 and 5. see source/kewald.f for a list.           
fft-package           FFTW
#
#  Set Parameters for Induced Dipole Convergence
#
polar-eps             0.0001                # the induced dipole convergence threshold
polar-predict
#
# the bond parameters between atom classes 1 and 4 are redefined and 
# will overwrite those already in the amoebapro13.prm. 
#this is only for illustration purpose
bond          1    4          200.00     1.1

How to generate input files of your own

proteins, nucleic acids, common organics

In tinker/params, you can find pre-existing parameters for certain molecular systems, both AMOEBA, amber, charmm, opls, mmff, and mm2/3.

AMOEBA library

For AMOEBA, please use amoeba09.prm for common small molecules, amoebapro13.prm for proteins. Nucleic acid parameters coming soon (end of 2017). More information here: http://biomol.bme.utexas.edu/tinker-openmm/index.php/TINKER-OPENMM:Amoeba

If you have a protein, you use "pdbxyz" to convert it to tinker xyz file. It will ask you for a key file (1bty.key below) or you can create a key file with "parameters $TINKERDIR/params/amoebapro13.prm" in it:

pdbxyz 1bty.pdb -k ./1bty.key

AMOEBA for a new ligand

Note that pdbxyz recognize proteins, water (res name HOH and some ions). The ligand (benzamidine above) is stripped. For that you need to derive your own parameters. For AMOEBA this can be done using POLTYPE: https://biomol.bme.utexas.edu/wiki/index.php/Research_ex:Poltype

Input: ligand.sdf or ligand.pdb (for example, you can truncate the ligand "BEN" out of the 1bty.pdb) 
Output: will produce the xyz and corresponding key files (ttt.xyz and ttt.key):

If you will merge the ligand xyz file with another molecule, make sure you set the suitable range for atom types (an option for pOLTYPE) so that they won't overlap. See the "Check the results section" on POLTYPE website before using them. Use "analyze.x" to make sure the xyz and key files work correctly: "analyze.x ttt.xyz -k ttt.key ep"

Combining two xyz files

Tinker "xyzedit.x" program (option 20) can be used to combine two xyz files into one (matching key file for each xyz is required). You can merge the key or parameter files by appending one to the other, but please make sure the atom types are not overlapping between the two.

For the above 1BTY example, poltype will translate/rotate the ligand in standard orientation. You may need a script to copy the coordinates from ligand.pdb into ttt.xyz (hwich you get out of the POLTYPE run).

Then you can use xyzedit.x to combine 1bty.xyz and lig.xyz into one xyz file with the two molecules orient/position as in PDB. Merge the key file by appending ttt.key (except the first line which contains header already in amoebapro13.prm of 1bty.key) to 1bty.key; again avoid overlapping atom types between the ligand and protein.

Soaking solute in solvent

To soak a molecule/complex in a box water, you can use this utility: http://www.ime.unicamp.br/~martinez/packmol/home.shtml

Alternatively, tinker "xyzedit.x" program (option 20) can be used for this as well. You just need a molecule.xyz and waterbox.xyz along with matching key files (actual names do not matter). One can make a water box of different size using this program too (starting from one water or a cluster of water and use option 19).

How to run Tinker and tinker-openmm

Command line

Tinker programs can be run interactively, which is the best way to learn what are the required inputs. Tinker programs can also run in background with all parameters specified, for the purpose of automation:

analyze.x ttt.xyz -k ttt.key ep
dynamic.x bench7 100000 3.0 6.0 2 298.0 N > ben7.log &
# the bench7.xyz above can be found in tinker official distribution inside tinker/bench/
#2-fs time step for MD here is ok because of the "integrator respa" in the key file

Or for tinker-openmm:

 #export CUDA_VISIBLE_DEVICES=0
 #stable 6:
 nohup dynamic_omm.x bench7 100000 3.0 6.0 2 298.0 N > ben7.log &

Additional notes for "tinker-openmm"

*It s recommended to use the respa inetgrator and 2-fs time step
*"heavy-hydrogen" in key file allows a 3-fs time step
*Default and only thermostat is Anderson
*Only MC barostat is available for now. We are adding virial to openmm.

How to specify different ensembles for MD simulation

For TINKER-openmm

Not many options available. See above.

Available thermostat and barostat in TINKER (2/4/2016 Ponder)

Thermostats:

METHOD                     KEYWORD
Bussi-Parrinello           thermostat bussi
Berendsen                  thermostat berendsen
Andersen Stochastic        thermostat andersen
Nose-Hoover                thermostat nose-hoover

There are only two barostats available via the “barostat” keyword:

METHOD                     KEYWORD
Berendsen                  barostat berendsen
Monte-Carlo                barostat montecarlo

The above thermostat and barostats are available for the Verlet, Beeman and RESPA integrators, and can be used in any combination with those integrators. These integrators are available via:

METHOD                     KEYWORD
Verlet                     integrator verlet
Beeman                     integrator beeman
RESPA                      integrator respa

Note that the defaults are Bussi for thermostat, Berendsen for barostat, and Beeman for integrator.

Then there are two special integrators, a stochastic one, and a Nose-Hoover that does NPT. The stochastic integrator uses a kind of Langevin temperature bath for thermostating, and does listen to the barostat keyword. The Nose-Hoover integrator uses a separate code branch and does only NPT with Nose-Hoover methods following Martyna-Tuckerman-Klein. You can get these via:

METHOD                      KEYWORD
Stochastic                   integrator stochastic (no need for other T control)
Nose-Hoover NPT             integrator nose-hoover (no need other keywords for T or P; starting structures need to reasonable)

Recommended keywords (add somewhere in the .key file) for NVT

integrator respa
thermostat bussi

Command:

dynamic xxx.xyz 100000 2.0 1.0 2 298  (100000 steps, 2.0 fs time step, dump structure every 1.0 ps, option 2 is NVT, 298 is the target T)

It is also possible to combine "integrator Beeman" or "thermostat Berendsen” or "thermostat Andersen". But RESPA allows large time steps (2.0 or 2.5 fs) than Beeman. Berendsen thermostat does not provide canonical ensemble fluctuation.


Recommended keywords for NPT

Integrator nose-hoover # this will trigger also nose-hoover P and T control

Note with Nose-Hoover you can not use a time step bigger than 1. 0.5fs is the most stable

Good alternative: Berendsen produces correct average but not "correct" ensemble fluctuation (see 2003 water paper)

integrator stochastic
thermostat bussi
barostat berendsen

Command:

dynamic xxx.xyz 100000 2.0 1.0 4 298 1.0 (option 4 is NPT, 298 is T and the last 1.0 is the target pressure 1atm)

Barostat Berendsen is available and useful for equilibrate system but again not recommended for production since the ensemble fluctuation is incorrect.

NVE

integrator respa

command:

dynamic xxx.xyz 100000 2.0 1.0 1

No need for thermostat or barostat of course. It is best to use smaller time step such as 1.0fs to conserve energy better. May even use smaller polar-eps (10^-6) than default (10^-5) in the key file.

Non periodic system, e.g. gas molecules not in a box

integrate stochastic

If there is no box (a-axis) in the key file, or no box dimensions in the den file, the system is non-periodic. This will set the stochastic temperature control along with the stochastic MD integration. For gas phase molecular cluster (very few atoms), the recommended time step is 0.1 fs.

Free energy calculations

BAR

Alchemical free energy calculations are available in TINKER and TINKER-OpenMM. One needs to specify the ligand or solute in the key file (example below). The lambda scaling schedule can be specified by user, automated by the "bar.x" in TINKER. See this reference for examples: J Comput Chem. 2017 Sep 5;38(23):2047-2055

For host-guest, one can also apply a bond restrain between host and guest. Two corrections may be needed for this:

*when lambda=1 for both ele and vdw. You can remove the effect of restraint from FEP or BAR; or you can avoid this correction if the restraint strength is gradually set to 0 when the lambda approaches 1.
*The correction for lambda=0 from RTln(C0*V), where C0=1/1660 A^3 and V=integrate{4*pi*r^2*exp[-k*(r-r0)^2/RT]dr}. k=15 kcal/mol in the example below. If you are using a harmonic restraint (example below) and the equilibrium r0 is not 0, numerical integration is necessary. A good reference is JACS v126, NO. 24, 2004.
ligand                  -1 24, 26    #negative means range: atoms 1 to 24 plus 26
vdw-lambda              1.00
ele-lambda              0.50
GROUP 1 -1 24
GROUP 2 -25 150
RESTRAIN-GROUPS  1  2  15.0  2.0 2.0
#lower bound & upper bound do not have to be the same

OSRW

Only implemented in tinker CPU. GPU version is under development.

Visualization

Force Field Explorer

By Mike Schneider and Jay Ponder, http://dasher.wustl.edu Can visualize the xyz and arc (MD trajectory) files; create and start TINKER calcualtions

VMD

Choose TINKER format when open a xyz file. Trajectory file (.arc) also works.

Pymol

Sometimes the xyz file can not be displayed correctly

Other resources

https://sites.google.com/site/biomolsimstw/workshop-materials/tutorials

Please email me if you have tutorials related to AMOEBA or Tinker you would like to share.