Research ex:Tinker OLD

Old version of Tinker tut page with some defunct stuff.

How to specify different ensembles, choose thermostat etc in TINKER for MD simulations 

You need a structure (e.g. protein.xyz) and matching key file (protein.key or tinker.key) to run MD. A typical key file for MD simulations looks like this (borrowed from tinker/bench/bench6.key in any distribution)

parameters            ../params/amoebapro13.prm
# the amoebapro13.prm is the built in protein force field in tinker/params. If you can a ligand you can add the parameters below in the key file instead of built-in prm file.
verbose                                # printing info for every step, for debugging mostly
randomseed            123456789
neighbor-list                       # this 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 quadrupoles die off faster.
pme-grid              64 64 64   #this 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.00001   # the induced dipole convergence threshold
polar-predict

Available thermostat and barostat (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:

integrator respa
thermostat bussi
Barostat MONTECARLO

or *this is questionable; berendsen doesn't give correct ensemble fluctuation

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. "Integrator Bussi” is available and behaves similar to “integrator nose-hoover” in the sense that it will trigger Bussi T and P control, but it is not as well tested.

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.

OSRW in TINKER

OSRW in tinker

TINKER binding free energy calculations

 Background reading

Alchemical (Sara)

Example Files for CCP Free Energy Calculations located on

http://biomol.bme.utexas.edu/~syc384/example_CCP_BAR_abs_solvation_free_energy/

Tutorial on how to calculate absolute solvation free energy for c9 ligand

SET UP TINKER/AMOEBA simulations

PDB to xyz

• Make sure the PDB files are standard, the columns are in the right place etc.
• Check the residue names. For example, LYP is not recognized. Change it to LYS. CYN is not, change to CYS. CYD is cya in disulfide bonds.
• Use HOH for water residue name.
• Missing H can be automatically added (not thoroughly tested).

TINKER 4 to 6 format changes

README First

Certain functional forms have changed from TINKER v4.2 to v5.x, so does the parameter formats in AMOEBA. A script is available to inter-convert between the two: donwload here

header in amoeba.prm and amoebapro.prm

The only changes seem related to opbend

v4.2

opbendunit              0.02191418 (1./radian**2 is the default; radian=57.29577951308232088d0)

v5

opbendtype              ALLINGER  
 !! tinker6 defaults to W-D-C; so you must have this keyword when using amoebapro.prm or poltype parameters
#opbendunit                 
!! Use default. POLTYPE ootput uses header file that include  "opbendunit    0.02191418".  As a result the force constant is about 71.94 time larger in tinker v6 prm files (which uses default for opbendunit) than that of POLTYPE. To combine POLTYPE with amoebapro etc, please multiply your opbend constant by 71.94 (and no "opbendunit              0.02191418" in the header)
opbend-cubic          -0.014    !! these coefficients were not in v4.2 where angle-cubic etc  but not opebend-cubic etc were used
opbend-quartic       0.000056 
opbend-pentic         -0.0000007
opbend-sextic          0.000000022

stretch-bend

Parameter file that works with tinker 4.3: strbnd 33 18.70 11.50 0.00. For a middle atom of type "33", there is one force constant for both stetch-bend coupling terms. One of the three numbers is used depending on the no of H atoms attached to 33.

Parameter file that works with tinker v5: strbnd 34 33 35 -4.50 38.00. The 33 is the middle atom here and 34 and 35 are bonded to 33; the two force constants are for the two stretch-bend coupling (34-33 and 35-33) respectively.

In tinker v4.2 estrbnd.f, estrbnd1.f (force & virial)

e = stbnunit * force * dt * dr
force...
virial

In tinker v5.0

e = stbnunit * (force1*dr1+force2*dr2) * dt
force ...
virial...

In the current sander amoeba_valence.f

Frc_K = degrees_to_radians*force_constant(it)
energy = energy + Frc_K*arg1(ind)*arg2(ind) 
...

In pmemd (amber v9?)

energy = energy + Frc_K * arg1(n) * arg2(n)
...

opbend

The main change is to make sure all 4 atom types involved in one opbend match the definition in the parameter file (see below). Also note the change in force constant units as per opbendunit in the header above; to convert force constants in the old units to the new units, multiply by 71.94.

tinker v4: opbend 78 30 0.200.

tinker v5: opbend 30 78 0 0 0.200

In the new format, the second of the four atoms classes is the central atom. The first atom class is the one that is out-of-plane. The third and fourth classes are the other two classes of atoms attached to the central atom ( "0" is a wildcard value that matches any atom class). So the old opbend parameter that was "B A" is now given as "A B 0 0".

Tinker source: v4:

e = opbunit * force * dt2
  &                * (1.0d0+cang*dt+qang*dt2+pang*dt3+sang*dt4)

v5:

 e = opbunit * force * dt2
  &                * (1.0d0+copb*dt+qopb*dt2+popb*dt3+sopb*dt4)

cange, qang, or copb, qopb... are the coefficients in the header.

Also in amoebapro4.prm we have the header “opbendunit 0.02191418”, so the actual force constant is “smaller” than those in amobeabiov4.prm which doesn’t have the opbendunit in the header file. It then uses the default: 1/radian^2, where radian =57.29.

Note that to use amber/PMEMD, the “opbendunit 0.02191418” is required in the tinker prm file since amber is hardwired to this.

opbend_force_in_prm= actual_opbend_force_used* (radian**2/opbunit)

When your opbendunit is 0.02191418 instead of 1/57.3^2 (by not defining the opbendunit in the header), your force constant in prm file should be reduced by a factor of 71.9.

In amoebapro13.prm, however, there is no "opbendunit 0.02191418" keyword, so the force constant is increased by a factor of 71.9.

When mixing prm files with and without this "opbendunit 0.02191418" keyword, please watch out the force constants.

Also it is imperative to use "opbendtype ALLINGER" keyword in prm file when using the latest tinker v6 since the default type is set to something else "opbtyp = 'W-D-C'"

torsion-torsion

protein backbone continues to use the original periodic 2-d spline to match AMOEBA with QM on the alanine dipeptide conformational energy. for sugar in nucleic acids, there is now a non-periodic 2-d spline used.

parameter definition/functional form/subroutine?

polarize

damping factor (0.39 except for few divalent ions) now explicitly included in the polarizability definition

polarize   337 0.696 0.390  333

new local frame

We added a new "z-bisector" frame definition in TINKER (since v6). It is not yet on PMEMD. This is mainly used for molecules such as N in H2-N-CH3 to avoid the use of dummy lone-pair solely for frame definition. The the N-C bond defines Z, and the bisector between the two NH bond defines the X.

What needs to be done:

• read the parameter correctly
• rotate the local multipoles to the global frame (energy/force are all computed using global multipoles so no change there)
• the hard part (at least for me) is to compute the force due to torque. In tinker, the torque is onverted back to force (I will post the block of code but is not useful here) but Tom uses the infinitesimal rotation to compute these forces directly. If I understand correctly we need derivatives of the rotation matrix (which we have of course) here. It would be much easier if Tom can help here.

The 4th  atom in the z-then-x definition is  used to indicate its chiral and flip the signs of some multipoles 
(see chkpole.f in TINKER) if necessary. It should not matter unless you actually have a right handed residue. 
You can just change the sign of atomic multipoles if you do have a right handed residue and 
you don’t even need the 4th frame defining atom type.

Also we add some new frame definitions in TINKER since tinker v6 such as the Z-Bisect 
which is not yet implemented in pmemd. You need to be careful about these even though they are not common. 
Here is how you tell which frame is used based on the sign and number of atom types (in kmpole.f of tinker):

Multipole_type   kz_type    kx_type   ky_type_if_needed
            if (kz.ne.0 .and. kx.eq.0)  axt = 'Z-Only'
            if (kz.lt.0 .or. kx.lt.0)  axt = 'Bisector'
            if (kx.lt.0 .and. ky.lt.0)  axt = 'Z-Bisect'
            if (max(kz,kx,ky) .lt. 0)  axt = '3-Fold'

Note you can always convert atomic multipoles between frame definition using TINKER poledit. 

Last, the best way to make sure is to compare TINKER and pemmd energy and make sure they match.

amoeba_par

This is the program that reads tinker analyze output and create the amber prmtop etc. Once I know what are the new variables and format that will go into sander/pmemd, I can change this accordingly.

additional mod that might be useful

• soft-core for buffered-14-7
• a script that automates the hydration free energy calculation of a small molecule using amobea and FEP; analyze the results using Bennett Acceptance ratio.
• a script for protein-ligand binding energy calculation is under way. The FEP calculation is similar to existing one in amber only the electrostatic part is not pairwise so we have to scale the parameters directly.

converting tinker to amber (get prmtop and inpcrd) trouble shooting

• The default max number of atoms in each polarization group is 20. If you have a big ligand with more than 20 or 22 atoms in a single polarization group, you should change the prmgrp.i file and analyze.f file and then recompile it.

In polgrp.i file, increase the parameters, i.e.

     integer maxp11,maxp12
   integer maxp13,maxp14
   parameter (maxp11=50)
   parameter (maxp12=100)
   parameter (maxp13=100)
   parameter (maxp14=100)

In the analyze.f file

 690             format (/,' Dipole Polarizability Parameters :',
  &                    //,10x,'Atom Number',9x,'Alpha',8x,
  &                       'Polarization Group',/)
            end if
            write (iout,700)  i,ia,polarity(i),
  &                           (ip11(j,ia),j=1,np11(ia))
700          format (i6,3x,i6,10x,f10.4,5x,22i6)
         end if
      end do
   end if

Change the line" format (i6,3x,i6,10x,f10.4,5x,20i6)" to "format (i6,3x,i6,10x,f10.4,5x,30i6)"

AMOEBA with Implicit solvent

PMPB

• PMPB
  *GK
  *Apolar (cavity+dispersion)
  *Simple SA approach

* Previous work using GK:
** Ramachandran map from all-atom+GK http://biomol.bme.utexas.edu/wiki/index.php/Research:All_atom_dipeptide_sol_implicit
** Anything else from biomol wiki?

In tinker5 from Mike 2007 (~schnied/Software), Bondi offset is used to scale or the atomic radii rsolv(i); in ksolv.f: rsolv(i) = rsolv(i) * offset

In tinker5.1,rscale replaces BONDI (default radtype) offset; in ksolv.f:

rscale = 1.0d0
if (radtyp .eq. 'VDW') rscale = 1.0d0
if (radtyp .eq. 'MACROMODEL') rscale = 1.0d0
if (radtyp .eq. 'AMOEBA') rscale = 1.0d0
if (radtyp .eq. 'BONDI') rscale = 1.03d0 #this is default 
if (radtyp .eq. 'TOMASI') rscale = 1.24d0
do i = 1, n
rsolv(i) = rsolv(i) * rscale
end do

Mike's comment 2010:
Yes - the apolar term is in TINKER5 and is done. It is computed as a
>> sum of cavitation and dispersion terms. The cavitation term is
>> proportional to volume for small molecules (very rigorous) and is
>> switched to surface area for protein sized solutes. The dispersion
>> term is based on integrating the Weeks-Chandler-Andersen integral
>> over a water continuum. There are 2 versions of the dispersion term
>> in TINKER5, one based on a GB like pairwise descreening integral
>> (subroutine ewca) and one based on numerical integration via Onion
>> shell (subroutine ewcax). For small molecules the two dispersion
>> methods will be nearly identical, but for proteins they may be
>> different due to interstitial spaces in the protein. For binding calculations, this sort of error will largely cancel.
>>
>> - Mike

>The apolar terms can use the vdW radii (ie. rsolv(i) = rad(class(i))).
>
>For Cavitation:
>There is the volume based cavitation coefficient, which is like a
>"solvent pressure" and has units of kcal/mole/A^3 (solvprs = 0.0327d0).
>There is the limiting surface tension for large molecules in units of
>kcal/mole/A^2 (surften = 0.080d0). Note that this is intentionally
>lower than the macroscopic, experimental value of 0.103d0.
>
>For Dispersion:
>There is a parameter called "dispersion offset" and units of Angstroms
>( dispoff = 0.26d0). The continuum water dispersion integral starts at
>"rsolv
>+ dispoff".
>
>For Electrostatics:
>The radii for continuum electrostatics are usually a little
>non-physical to implicitly capture hydrogen bonding or other non-continuum effects.
>
>Basically what I think needs to be done is we need a large table of
>experimental solvation energies to compare with Cav + Disp + PMPB/GK
>and with explicit water AMOEBA FEP. The question is how many
>electrostatic radii parameters do we need to get a reasonable unsigned
>error. Some people have gotten mean errors of only 0.35 for fixed
>charge models (AGNP), but I have a hard time believing this is not
>over-fit. I think going for 0.5 kcal/mole mean unsigned error would be
>terrific; although this is would still may be a lower error than the FEP results.
>
>Any thoughts?
>
>Best,
>Mike

to use simple surface area term in GK

For binding affinity a simple non-polar model is a good idea; there are currently no hooks in TINKER to let one do this from the key file. It is pretty easy to modify the code to do it though – just add “solvtyp .eq. GK” to the first if statement that selects SA only instead of the cav + disp code (in the esolv*.f files).

c c total solvation energy for surface area only models c

     if (solvtyp .eq. 'GK' .or. solvtyp.eq.'ASP'
  &   .or. solvtyp.eq.'SASA') then
      call surface (es,aes,rsolv,asolv,probe)
      ecav = es
      edisp = 0.0d0

For and SA only model, the surface tension needs to smaller than for a disp + cav model so I added the following to the key file:

1. Macroscopic is 0.103 kcal/mole/A^2; most SA only models use something less than 0.01

surface-tension 0.007

Bio force fields and revision history

amoeba bio forcefields