Research:AMOEBA Amber Free
README FIRST
Right now, only tinker v4 parameter files can be used by amber. Use these to convert from TINKER to amber:
/home/pren/tinker/bin/analyze
/opt/amber/bin/amoeba_parm
The prm tinker key has to be in tinker4.3 format (no 0.39 in polarize key)
Afterwords you can use amber9, pmemd9 exe with softcore/ion modifications. ) The executable may not always be in /opt/amber/bin, see oscar's document on AMBER mod for exact location.
To go from tinker to AMBER, you need:
mol.xyz mol.key mol.pdb
((pdb must match xyz)
1. Make a water box to soak the protein
You can build a big box with small boxes using the command: crystal.x. Make sure the box is big enough to contain the protein--side length should be at least 12 angstroms longer than protein's dimension.
If you want to speed up the simulation you can also chop off the corners of the cubic box to make a octohedral water box (see 10.)
2. If the box is a little too big, delete some water molecules in the end of the xyz file of waterbox proportionally.
Do not forget to change the atom number in the first line. Then run command dynamic to make the box shape back to normal.
Put this keyword "octahedron" in the tinker.key file or whatever key file you are using before dynamic or minimize in order to transform the box from cubic to octahedron.
3. Soak the protein in the water box.
Command: xyzedit protein.xyz
Option 11: to translate the center of mass to origin in order that the protein soaked in the box completely.
Option 17: input waterbox.xyz, soak the protein in the water box.
This step will generate a new xyz file of the whole system, we name it pro_wat.xyz
Make sure the a-zxis in the tinker key file will give the right density (1.00 g/cc). later when you run MD, if the box size change din dyn file, make sure to change the key accordingly.
4. Optional: Run minimize pro_wat.xyz, and then dynamic pro_wat.xyz.
REN: freeze the protein position (see manual for positional restrain or fix atoms). This is to reduce the disruption of protein structure from non-equilibrated water.
You can specify the atoms/residues you want to fix during the simulation in the mdin file by adding two keywords restraint_wt and restraintmask.
restraint_wt is the strength of the restraint. restraintmask is the atom selection. For example if you want the strength to be 10 and the protein you are to fix has 224 residues, include the following in the cntrl block of mdin file:
restraint_wt=10.0,restraintmask=':1-224'
How to relax a system using TINKER:
- For acceleration, turn mpole and polarization off at first by adding to the key file and them "minimize" to a convergence of 5.0:
neighbor-list
mpoleterm none
polarizeterm none
#polar-eps 0.00001 # uncomment this when uncomment polarization
#polar-predict # uncomment this when uncomment polarization
openmp-threads 8
printout 1 #so you can see the move
cutoff 9.0 #it is also faster not use Ewald for large box > 60 A
#but in case if you have convergence issue with polarization in minimization, you may have to use Ewald.
#Ewald (with this commented out, no Ewald will be used. Faster for initial refinement but needs uncomment for final relaxation)
#position-restrain proteins & ligands ("inactive” has trouble)
restrain-position 1 x y z 10.0 (restrain atom 1 at x y z location with a force constant=10 kcal/mol/A; one for each atom, write a script to restrain all proteins)
##Note you may be able to skip the protein-restrain step if heat up slowly (below).
....
- Minimize to 5.0 with protein/ligand inactive/restrained, with mpole but no polarization
- Turn both mpole and polarization (polar-eps=0.05) back on and minimize to 5.0
- remove restraints on P/L and minimize to 5.0
- Uncomment Ewald, polarization keywords, run minimization to 5.0
- Run NPT MD for 10000 steps to gradually heat up the system from 50K, to 110K, 150K, ...300K (eg. 5000 steps at each T)
- Run NPT at 300K for 500ps-1ns to equilibrate and calculate average volume before switching to NVT
5. Take the last structure of tinker and make it the xyz file. Generate PDB file that macthes the xyz
Do xyzpdb to this xyz file with the protein's seq file which is generated when doing pdbxyz to the protein.
It turns out atom order will be changed within residues in pdb file. There is one way to correct this. Do xyzpdb to the whole structure without seq file. This will give the right atom order. Now copy the atom column together with coordinate columns to substitute these parts in the former pdb file.
Also you can use the script ~yues/bin/pdb_crt to fix this error.
6. Now we have dyn (dyn is optional) pdb and xyz files. Ready for amber
Do analyze [ ] PC > [].analout
comment out the additional keywords you added for restraint and turning off polarization etc. Only analyze from TINKER4.3 work for now (/home/pren/tinker/bin
Note: when completely turning off electrostatic parameters of a group of atoms (or a molecule), the polarizability can be set to exact 0.0, but the permanent multipole needs a tiny charge has of 0.0000001 (all dipole and Q are eaxct 0), otherwise analyze PC won;t print out any multipoles/polarizability for theses atoms. Thus amber is confused becuase natom misatches the parameters and will bomb. Solution: Use a tiny charge with 6 zeros; This small charge will make tinker print out 0.00000 so they are exactly zero but it is printed in the analout.
When there is only one molecule and/or all atomic polarizabilities are zero; amoeba_parm will also bomb, complaining no polarizability found. Can use 0.0000001 for poalrizability Check the analout!
7. /opt/amber/bin/amoeba_parm -name [ ] -title "...." > [].amoeba_parm.out
you may need to set intel v9 path to use the amber_pre9 from Darden
This step will also generate input coordinate file(.inpcrd) and parameter file (.prmtop). Remember you need to have the key file with the same name as amoeba_parm.out, and also note that in the key file, the keyword "parameter" can't have capital P.
8. Write a mdin file and then run sander to do simulation
/opt/amber/bin/sander -i [].mdin -c [].inpcrd -p [].prmtop -o mdout &
9. Use AMBER to further relax the system before the production run.
First minimize the structure well
Then gradually heat up the protein from 0K to 300K in a series of MD (10 ps at 100K, another 10 ps at 200K...). to restrain pro/lig in PMEMD: Group input for restraint residue 2.0 RES 1 365 END END
Run NPT dynamics at 300K for at 100 ps or more to get the equilibrium density. Then ready for NVT using the average density. Please check your protein/ligand structure at every step by comparing against original PDB structure, to make sure nothing is screwed up.
Octahedron
If the water box is octahedron, there are some differences from cubic box. [Ren: this is because TINKER and amber define octahedron differently so the density drops when moving from TINKER to amber].
10. Add the keyword "octahedron" in the tinker key file. Soak the protein in an octahedron box and go through step 1-6 above.
11. Since the definition of octahedron size is different in tinker and amber we need to modify some files.
Change the last two lines in inpcrd file. (1) Angle: from 90 to 109.4712; 0.109471219000E+03 0.109471219000E+03 0.109471219000E+03
(2) Make the size big at the beginning and then shrink it otherwise SANDER bomb complaining too many subcells.
for example, the tinker size is 29 A, change it to 30 A in .inpcrd will let sander to run, but density is only 0.56 (see shrinking box below).
In the prmtop file there are two places need to be changed. In the last second row of FLAG POINTERS,
change the value 1 (the last third) to 2 meaning octahedron.
%FLAG POINTERS
%FORMAT(10I8)
1170 1 1 1 1 1 1 1 0 0
1170 390 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 2 0 0
0 0
%FLAG ATOM_NAME
ALso change the angle (search BOX) from 90 to 109 0.90E2 to 0.10947122E+03
12. Shrink the box before long simulation.
(1) Run NPT simulations(300 fs) to shrink the box (remember to freeze everything but water); stop when density = 1.0.
To run NPT: In the mdin file add this: ntb=2, ntp=1, taup=1.0 press0=500 to run NPT instead of NVT. 500 atom is used to squeeze the box. Alternatively, one may run a series of short NVT and manually reduce the box size in inpcrd: Every time overwrite inpcrd file with restrt file and decrease the size a little bit.
To FREEZE THE PROTEIN and ION and LIGANDS with positional constrain on protein/ion/ligand:'
from AMBER 9.0 Drug/DNA Complex Tutorial tutorial:
Initial restrained minimization on DNA & DRG, 10 cut &cntrl imin = 1, maxcyc = 1000, ncyc = 250, ntb = 1, ntr = 1, cut = 10 / Hold the DNA & DRG fixed 500.0 RES 1 25 END END
Hold the DNA and the DRG fixed 500.0 (This is the force in kcal/mol used to restrain the atom positions.) RES 1 25 (Tells AMBER to apply this force to residue #’s 1 to 25). Then let go the restrain and relax the system for another period
(2)Let go restrain and Follow instructions at the step 9 above to equilibrate the system.
Compare your current structure (protein/peptide/ligands) and original structure to make sure RMSD is small (ie. npt didn't ruin the structure)
(3)Check the structure with crystal frequently to see if the ligand is still in the pocket, me important ions are still in the same place and the RMSD of the protein is reasonable. You can use software as VMD to visualize the two structures at the same time and compute their RMSD. Normally backbone RMSD should be around 1.0 Angstrom.
(4)It is time for long NVT simulation.
In the mdin file, delete ntb=2, ntp=1, pres0=1.0 taup=50. BE CAREFUL AT EACH STEP. CHECK EVERYTHING THAT CAN INDUCE ERROR!!!
How to calculate free energy between step i and j:
four energy files are needed:
ii.e ij.e jj.e ji.e
The first letter indicates the trajectory and the 2nd is the parameter.
To generate the above energy files: For given trajectory, run
- get energy from current parameters
/opt/amber/bin/sander -i ele1.min -c ../ele1.inpcrd -x ../mdcrd -p ../ele1.prmtop -o ele1.e &
- get energy from perturbed parameters
/opt/amber/bin/sander -i ele1.min -c ../ele1.inpcrd -x ../mdcrd -p ../ele2/ele2.prmtop -o ele2.e &
ele1.min is the same as ele2.min. No need for restrain as they cancel.
Then run "BAR-amber ii.e ij.e ji.e jj.e start_frame end_frame guessed_init"
- energy file with first trajectory and first prm file
$file00=$ARGV[0];
- energy file with first trajectory and second prm file
$file01=$ARGV[1];
- energy file with second trajectory and first prm file
$file10=$ARGV[2];
- energy file with second trajectory and second prm file
$file11=$ARGV[3];
We also have a program to evaluate the statistical error of dG from BAR (/home/pren/bin/barerror_bigstep.oscar.amber). It also tells you if the step is too big (error too big).
