Research: AMOEBA Ligand Par
Build your ligand
(1) Make your molecule with Chem3D Ultra or Arguslab (windows)
You can also use other program too.
(2) Minimize the structure. Once the molecule is built, minimize energy with any of the methods available in Chem3D. In this case MOPAC method was used. There are other ways to minimize energy: MM2, Gaussian, Mechanics. You can use any of them as long as it works.
(3) Save the structure as xyz file.
Save as Tinker MM3 Input (*.xyz), benz.xyz
Optimize structure with QM
(1) Convert .xyz to .com
Just copy the atoms and their coordinates into the .com file with following format:
%chk=benz.chk %mem=50MW %nproc=2 # opt hf/6-31g* opt energy 1 1 C -0.37610769 0.74520000 0.98303077 H -0.43910769 1.74320000 0.99203077 C -0.72410769 0.05120000 2.12303077 H -1.03010769 0.54320000 2.93803077 C -0.64810769 -1.30880000 2.13203077 H -0.89810769 -1.82380000 2.95203077 C -0.22510769 -1.96280000 1.00303077 H -0.17010769 -2.96080000 1.00603077 C 0.12589231 -1.26980000 -0.13596923 H 0.43389231 -1.77480000 -0.94296923 C 0.05989231 0.11420000 -0.18296923 C 0.40489231 0.82420000 -1.33896923 N 0.00089231 2.06920000 -1.47396923 H 0.24089231 2.58620000 -2.29496923 H -0.54610769 2.49820000 -0.75496923 N 1.12789231 0.22220000 -2.26396923 H 1.38389231 0.71320000 -3.09696923 H 1.41989231 -0.72480000 -2.13196923 <blank line>
keyword "opt" indicates it is running optimization.
You have to leave one blank line at the bottom of the file, and then save it as benz.com which will be the input file of Gaussian03
(2) Run Gaussian program to optimize the structure.
Do "g03 benz.com". benz.chk and benz.log will be generated. benz.log is the output file which
has the optimized structure and energy.
Single Point calculation
(1) Make the single point file (.com)
Use the script mkspcom2.pl, do "mkspcom2.pl benz.log"
This will extract the optimized coordinates and make it a single point file benzsp.com which
looks like this:
%Mem=500MB %Nosave %Chk=benhsp.chk %Nproc=2 #MP2/6-311++G(2d,2p) Sp Density=MP2 MaxDisk=960MW sp energy 1 1 C1 -0.219172 0.772015 1.045507 H2 -0.124354 1.841828 1.092625 C3 -0.556685 0.065459 2.183830 H4 -0.742000 0.589477 3.102333 C5 -0.648127 -1.318255 2.138044 H6 -0.916424 -1.863862 3.023554 C7 -0.395888 -2.003092 0.957999 H8 -0.479004 -3.072957 0.925387 C9 -0.043176 -1.306084 -0.181602 H10 0.125831 -1.839237 -1.099697 C11 0.041747 0.084714 -0.138928 C12 0.409652 0.833072 -1.353567 N13 -0.200543 1.960083 -1.635531 H14 0.060474 2.534099 -2.409732 H15 -0.999989 2.242440 -1.111118 N16 1.345776 0.368932 -2.147446 H17 1.581717 0.805673 -3.013818 H18 1.901226 -0.409704 -1.866286 <blank line>
We use higher energy level and basis set MP2/6-311++G(2d,2p) in order to get accurate energy result. Lower energy level has been used for optimization since the energy is sensitive to energy level while structure is not.
Keyword "sp" means this is a single point calculation.
Now run "g03 benzsp.com", again this will generate benzsp.chk benzsp.com
Note the gaussian input structure is rotated into standard orientation unless "NoSymm" keyword is specified. SO when you compare the molecular dipoles from gdma using TINKER with Gaussian log file, remember to use the coordinates (xyz file produced by poledit is correct) in standard orientation.
Multipoles calculation.
(1) Make fchk file from benzsp.chk
Do "formchk benzsp.chk". This will make benzsp.fchk file.
(2) Make a gdma input file.
Open up a new text file. Type following text in the file.
Title "benzamidine gdmain" File benzsp.fchk density MP2 Angstrom AU Multipoles Limit 2 Radius H 0.31 Punch benzsp.punch Start Finish
(3) Run GDMA to calculate multipoles of this ligand.
Do "/opt/gdma-2.2/bin/gdma < gdmain > benzsp.gdmaout". Wait for about 20 min. There are two output files benzsp.punch and benzsp.gdmaout, the latter one has the multipoles.
(4) Edit the multipoles from the GDMA output.
Do "poledit benzsp.gdmaout". Two options:
(1) Atomic Multipoles from GDMA Output File (2) Removal of Intramolecular Polarization
Choose the first one. Then it will ask you to "Enter Altered Local Frame Definition". If you don't need to change the local frame definition, Press Enter. If you have symmetric atoms that you will average the multipoles later, please make sure the frame definition is symmetric. This will generate benzsp.xyz and benzsp.key. The .key file has all the multipole values.
If you have a large molecule that is flexible, go to step 8.
Average the multipoles
(e.g over methyl H atoms, or C atoms of benzene carry very similiar multipoles but may be slightly different)
(1) Change atom index.
Remove the negative signs of the local-frame indices of the atoms. The negative means atom number (so you can run analyze without defining new types). Group the atoms manually, by making the indices of the same type of atoms identical (including reference indices).
(2) Run avgmpoles.pl to average multipoles
The script is in /home/oscar/bin/.
There are five arguments for this script. Multipule key file, index of first atom in key file, index of last atom in the key file, the starting index of this ligand that will appear in the TINKER parameter file. (For example, the last atom number of the preexisting parameter file is 273, then the first atom number of this new ligand is 274).
In the case of benzamidine, do "avgmpoles.pl benzsp.key_2 1 14 274 benzsp.key_3" This script also zero out dy, qxy, qyz components of multipoles for z-then-x. For bisector frame multipoles, please zero out dx also. ALl these components should be very small (0.0xxx), otherwise your bisector definition is improper, or it is a chiral center (like Calpha). We have special definition for chiral atom (check Calpha in amoeba.prm).
6. Create parameters of new ligand in TINKER parameter file amoeba.prm.
(1) Atom definition.
Number all the atoms of the new ligand sequentially at the end of amoeba.prm. Atomic types should also numbered following the last atom type above. Copy the atomic number and mass from existed parameters.
(2) vdw, bond, angle, out-of-plane and torsion parameters
Most parameters can be obtained from the existed AMOEBA force field (or MM3). Try to find the closest atom type, for benzamidine, it is benzene + amidine (from arginine). You may also fine-tune the valence parameters such as bond length and angle values by adjust them to match the QM-optimized structure that you already have. Unless you worry about vibrational spectroscopy, we leave the force constants unchanged. They can be improved however by matching the Vib Freq from QM.
Torsion may or may not be transferable. In our hand it has been OK as long as you distinguish single and double bonds and transfer from similar atom_types. If you really care about particular torsions, e.g. peptide backbone or sugar ring puckering, you can refine the trosion parameters by fitting to the QM energy profiles. This is done at the last step (i.e. after you get all the other parameters), then you can calculation the torsional profile using QM (restrained optimization at MP2/6-311++G** level, for example), and repeat the same calculation using TINKER (with the particular torsion parameter set to zero). Now fit the 3-term Fourier series to the energy difference E(QM-TINKER). Here is a gnuplot script we use to do this. This example fit two torsion (phi-psi) simultaneously, and can be easily modified for one torsion:
f(x,y)=a/2.0*(1+cos(x))+b/2.0*(1-cos(2*x))+c/2.0*(1+cos(3*x))+a1/2.0*(1+cos(y))+b1/2.0*(1-cos(2 *y))+c1/2.0*(1+cos(3*y))+g a=1;b=1;c=1;a1=1;b1=1;c1=1;g=1 f(x,y)=a/2.0*(1+cos(x))+b/2.0*(1-cos(2*x))+c/2.0*(1+cos(3*x))+a1/2.0*(1+cos(y))+b1/2.0*(1-cos(2 *y))+c1/2.0*(1+cos(3*y))+g fit f(x,y) 'lacphipsi2.csv' using 1:2:10:($7) via a,b,c,a1,b1,c1,g #splot f(x,y), 'lacphipsi2.csv' using 1:2:10 splot 'lacphipsi2.csv' using 1:2:10, f(x,y) pause -1
(3) Polarizability
Polarization parameters can be also tranfered from existed parameters. Normally, if you have a large molecule, you need to split up the molecule into polarization groups (each with less than 15 atoms but can be larger if needed).You also need to modify the connectivity of each atom. Remove the connection between polarization groups. But benzamdine only has 18 atoms, so we still consider it as a whole group. SEE step 8 for necessary modification of permanent multipoles if you have 2 or more groups in a molecule.
(4) Append the multipole parameters from step 5 to the amoeba.prm.
7. Analyze check.
(1) Assign the new atom types to each atom in xyz file.
(2) Run "analyze ligand.xyz em" to print out the potential energy and electrostatic properties. If any parameter like bond, angle is missing, it will complain parameters are not defined. Compare the electrostatic properties with QM calculation espesially total charge and dipole moment and x y z components.
8. Intramolecular polarization group.
If the the molecule is large and flexible, it is necessary to definite a few groups so that the permanent multipoles will polarize atoms in other groups but not atoms within the same groups. See JCC 2002 paper for detailed discussion. For example, a polypeptide backbone will consists of amide groups. Each groups is essentially a rigid fragment/functional group.
The group is defined through polarize keyword in tinker parameter file:
polarize atom_type_integer polarizability_real connected atom_types_integer
The connected atom_types are atoms connected to the current atoms AND belong to the same group. TINKER uses this information to build up groups.
Because in QM DMA, the multipoles would include both permanent and induced dipole moments (depending on how you define the group). Once the groups is defined in the tinker key or parameter files, you can use "poledit" option 2 to extract the permanent multipoles that go in the final tinker parameter file.
* Moficaton of step 5: don;t goup the same atoms (symmetric) yet. Assign uniuqe atom_type to each atom. Define poalrize and group. use poledit to extract (option2). * You can now use ~pren/avgmpoles.pl extracted_key first_atom_type last 1 new_key to average mutipoles of the same.
* Now continue with step 6 and 7.
9. If you want to fine tune the torsion see Step 6.(2). Otherwise you are done. Again please check the parameters to make sure it is right before jump ahead. Check the total dipole moments against QM (using the same structure of course). Check the polarization energy and group definition (analyuze ed and analyze ed). Superimpose the optimize dstructure against QM structure etc.
Note please use tinker 4.3 parameter style and executable (in ~pren/tinker4.3) if you want to use with AMBER.
We also have in-house modification of tinker 4.3 for handling Ca/Mg and Zn. Jay used different style which is not understood by AMBER yet. The main differences is that TINKER5 has explicit damping coefficient in "polarize" definition. We used polar-damp-ion in tinker4.3. We have modified amber pre-9 to be compatible with our approach.