Research:Amoeba par tut
Quantum Calculations
From initial structure to TINKER simulation
(uses Gaussian03, TINKER, and formatting scripts written by P. Ren and C. King)
Copy these files into your ~/bin directory from ~kris/bin:
avgmpoles.pl, mkMP2spcom.pl, indexyz.pl, indexmpole.pl, poledit.x, MP2template.com, template.gdmain
You need to set up your .cshrc file so you can run these programs:
Gaussian g03, TINKER analyze.x, formchk
Type “source ~/.cshrc” to update.
Initial Structure: Make molecule in your favorite program. Save file in Cartesian coordinate format. For this tutorial, the molecule's arbitrary name is “molecule.*”, but you can use any name wherever “molecule” appears. Try to get the molecule in a roughly accurate conformation to reduce computing time (e.g. “Minimize” function in ChemDraw3D).
MP2 Structure Optimization: Open a new text file (or use MP2template.com). Use this format:
-------------------------------------------------------------------------------------------------------
%Mem=900MB
%Nosave
%Nproc=2
%Chk=moleculeMP2.chk
#MP2/6-31G* OPT MaxDisk=960MW
molecule
0 1
C1 0.000 0.000 0.000
H2 1.000 1.000 1.000
-----------------------------------------------------------------------------------------------------------
The newlines are part of the syntax. Make sure you leave one blank line at the bottom of the file also, or the program will not run! Then just copy-paste the Cartesian coordinates of your molecule from step 1 into the bottom section. Save this file as moleculeMP2.com . This will be your input for quantum structural optimization. Next, login to a node that has no processes running (you will need BOTH processors), navigate to the directory where you saved the moleculeMP2.com file, and type “g03 moleculeMP2.com &” from the command line. Wait approximately 12-24 hours.
MP2 Single Point Calculation: Next, you will use the Cartesian coordinate output from the MP2 optimization as the input for MP2 single point calculation. The g03 program should have generated moleculeMP2.chk and moleculeMP2.log. Open moleculeMP2.log and look at the bottom of the file to make sure it has finished. It should have some inspirational quote and the job CPU time and some other stuff. If it is done, you can generate a new .com file using a script. Type “mkMP2spcom.pl moleculeMP2.log”. This will make a new .com file called moleculeMP2sp.com. Like step 2, navigate to a free node and type “g03 moleculeMP2sp.com &”. Wait ~24hrs.
(In the absence of a script, you can just make the new .com file yourself by replacing the basis set line with “#MP2/6-311++G(2d,2p) Sp Density=MP2 MaxDisk=960MW” and copy/pasting the new Cartesian coordinates out of the end of moleculeMP2.log and into the new .com file and rewriting the atom type and number in the correct syntax.)
4. File Conversions: When it has finished running, the g03 program should have generated moleculeMP2sp.chk and moleculeMP2sp.log. Type “formchk moleculeMP2sp.chk”. This will generate moleculeMP2sp.fchk. Now open a new text file (or use template.gdmain) and make this:
--------------------------------------------------------------------------------------------------------
Title "MoleculeMP2sp gdmain"
File moleculeMP2sp.fchk density MP2
Angstrom
AU
Multipoles
Limit 2
Limit 2 H
Radius H 0.35
Punch moleculeMP2sp.punch
Start
Finish
--------------------------------------------------------------------------------------------------------
You should overwrite “molecule” with your own molecule name in this file so it points to your .fchk file and generates the filenames you want. Save this .gdmain file. Now type “/opt/gdma-2.2/bin/gdma < moleculeMP2sp.gdmain > moleculeMP2sp.gdmaout &”. Wait about 30 minutes. If the new .gdmaout file ends with the CPU runtime, then it is done. Now type “poledit.x moleculeMP2sp.gdmaout” then type “1” at the options prompt. Press enter a couple times. This should have generated moleculeMP2sp.key and moleculeMP2sp.xyz. The moleculeMP2sp.key is not a key file, but contains TINKER format atomic multipoles for use in creating a new parameter file. To prevent confusion, rename this file moleculeMP2sp.mpole.
(see Appendix 1: How does amoeba.prm formatting work?)
Parameter File Editing Part 1a: Open up amoeba.prm for editing. First, you need to create a new atom type for each atom in your molecule. You can just append the new types to the end of the “Atom Type Definitions” section. To make things easy on yourself later, make the atom types correspond sequentially to the index numbers of the associated .xyz file. (For example, atom 1 becomes type 300, atom 2 becomes type 301, and so on.) Assign each atom a preexisting class number. Multiple types can have the same class number. Just look at the preexisting atom type definitions and use an atom class number of the same species that is in a similar local environment. Alternatively, if you have some atoms in local environments that are not well defined by any Amoeba classes, you can create new classes for them. Also make sure to redefine the binding partner number for each atom. Since your molecule probably has atoms with class numbers that are not usually connected in amoeba.prm, you will need to create new definitions for bond, angle, and torsion parameters at the very least. You may also need to define new opbend parameters for ring systems. You also need to redefine the ideal bond lengths and bond angles, since the Gaussian quantum calculations gave you an .xyz file with optimized coordinates. Use your favorite visualization program (e.g. PyMol) and measure all bond lengths and bond angles. Rewrite these values into the corresponding sections of amoeba.prm, averaging multiple values for the same class. If you have created any new classes, you will also need to make new vdw parameters for each new class. To get all these new parameters, you can just copy/paste from other definitions in amoeba.prm and then rewrite the atom class connectivity, assuming there are some preexisting definitions that are relatively similar to the situation you are defining. If amoeba.prm doesn't have any entries in a given section that are close enough to one of your atomic environments, you can take paramters from mm3.prm, which has a more extensive list of atomic environments. However, if you use mm3.prm parameters, you need to convert the k constant for bond length, angle, and torsion by multiplying by a factor of 71.94. Next, you need to modify the “Polarization” section. Unlike the other parameters you have modified, polarization is defined in terms of atom type, not class. Thus you need to make a new entry for each of the atoms in your molecule. Like before, you can just copy/paste polarization constants from elsewhere in amoeba.prm or mm3.prm. You also need to define connectivity of each atoms. If your molecule is less than about 15 atoms, just list the type numbers of every atom that each atom is connected to (for now, each new type number corresponds to a single atom, even if you have multiple similar atoms.) However, if you have a larger molecule, you need to split up the molecule into polarization groups, each with less than 15 atoms. In amoeba.prm, this is done by omitting the connections (just in the “Polarization” section) between atoms and all of their binding partners that are in different polarization groups. Thus, if atom type “100” is connected to atom type “101”, but “101” is in a different polarization group, you would leave out that connectivity in both atoms' polarization definitions.
Parameter File Editing Part 1b: Now you are ready to define new multipoles for each new atom type. These were generated by poledit.x and named moleculeMP2sp.mpole at the end of step 4. Since the atom type numbers in moleculeMP2sp.mpole are just the index numbers of the .xyz file (starting at 1), you need to change them to your new atom type numbers. You can do it manually, or use indexmpole.pl. Let us say that the atom of index number “1” in molecule.xyz corresponds to atom type number “287” in the modified amoeba.prm. To use the script, type “indexmpole.pl moleculeMP2sp.mpole 287 > moleculeMP2sp-index.mpole”. This will copy moleculeMP2sp.mpole into moleculeMP2sp-index.mpole, but will offset the type numbers to the correct numbers for use in your new amoeba.prm. For this to work right, the new atom type numbers must be sequential (the script simply adds a user-defined number to each atom type/index number). Now you can just open up moleculeMP2sp-index.mpole and copy/paste directly into amoeba.prm at the end of the preexisting “Atomic Multipole Parameters” section.
Analyze Check: To make sure you did everything right, you will run analyze.x using your new parameter file. Before doing this, you need to change the atom type definitions in moleculeMP2sp.xyz to your own type numbers. Because Perl is awesome, you can use a script to do this. Assuming (for example) that your new type numbers start at type “287”, type “indexyz.pl moleculeMP2sp.xyz 287 > moleculeMP2sp-index.xyz”. Then create moleculeMP2sp-index.key that points to your modified amoeba.prm. Your *.key file for the analysis needs to have this term somewhere in the body: “polarizeterm NONE”. This will turn off polarization so you can accurately compare the output with your *.log file from before. Now you can run TINKER's analyze.x on moleculeMP2sp-index.xyz using the new amoeba.prm. Type “analyze.x moleculeMP2sp-index.xyz -k moleculeMP2sp-index.key EM”. Look at the “Total Electric Charge” on the output. If you have a neutral molecule, this should be very near to 0. Open up moleculeMP2sp.log and scroll down to the bottom of the file. Check and see if the dipoles and dipole moment magnitude generated by g03 in the .log file are the same as those generated by TINKER in the analyze.x output. If they are the same or very nearly the same, you can proceed to step 8.
Remove Induced Dipole from Total Dipole: Since everything is working correctly so far, it is time to change the dipole component of the multipole parameters such that only the permanent dipoles contribute to the magnitudes of the dipole directions. So far, the dipole moments of your atoms, which were calculated by Gaussian, have consisted of a sum of induced and permanent dipole factors (Gaussian does not know the difference). This has been OK so far because you turned off polarization in moleculeMP2sp-index.key (in step 7). Since you want TINKER to be able to calculate new induced dipole moments from intermolecular interactions, you need to remove this extra component from the amoeba.prm dipole moments before doing simulations with polarization turned back on. To do this, you will use poledit.x again. First, open up moleculeMP2sp-index.key and turn polarization back on by deleting or commenting out the line you wrote earlier that says “polarizeterm NONE”. Now create a file that looks like this:
2 1mer_amideMP2sp-indx.xyz blankline
Save the file as “molecule.poleditin”. This will be the input file for poledit.x (it will automatically select option 2 in poledit.x). Make sure that your *.key file which points to amoeba.prm has the same name as your *.xyz file (e.g. “moleculeMP2sp-index.key”) or the program will not work correctly! Now type “poledit.x < molecule.poleditin > molecule.poleditout”. The contents of molecule.poleditout are the output of poledit.x (multipoles with polarization terms added in).
9. Atom Type Grouping: Now you will reduce the number of different atom types by grouping some of the similar atoms in your molecule together as single atom types and averaging their multipole moments. You might want to make a copy of molecule.poleditout before you edit it. Look at your molecule and identify which atoms look like they are in essentially identical environments (they must at least be in the same class!). Open up molecule.poleditout which was generated by poledit.x in step 8. Look at the multipole moments for the atoms you suspect are similar. If the numbers are fairly close (within 0.1), then you can probably label them as the same type number. However, the atoms defining the X and Z axes must also be the same. Let's say that you want to group type 1 and type 2 into a single type. If their multipoles are going to be described by a single entry in the amoeba.prm, then the axes of both atoms must be defined by the same atom types. That is, if atom type 1 axes are defined by atom types (5, 6), and atom type 2 axes are defined by atom types (7,8), then types 5 and 7 must also become the same type, and types 6 and 8 must become the same type. If it is not possible to make the atoms defining their axes the same type (they are connected to different kinds of atoms), then you cannot group them together. For each group of atoms that you wish to assign to the same type, just overwrite the atom numbers in the multipole moments file molecule.poleditout with the same type number. To keep things simpler, just use the lowest type number in the group for all the atoms. It's OK if you end up having gaps in the sequence of numbers.
10. Atom Type Averaging: Now that you have eliminated some atom types by grouping, you need to delete all the instances of those unused atom types in amoeba.prm. This includes changing the connectivity numbers in the “Polarizability Parameters” section to atom types that have not been deleted. Using your edited molecule.poleditout file with the same number of multipole entries but a lesser number of different atom types, you can average the multipole moments of identical atom types using a script. For example, say that you want the lowest atom type number (corresponding to type “1” in molecule.poleditout) to be type 287 in your .prm file. If your molecule is, say, 2 atoms big, then the range would be 287 288. Since our molecule is 31 atoms big, it is 287 317. Type “avgmpoles.pl molecule.poleditout 287 317 > AVGmolecule.poleditout ”. This will output a file with multipoles averaged over identical type numbers and indexed to the correct range. Most likely, there will be gaps in the sequence of numbers, but that doesn't matter. You can now copy/paste these into amoeba.prm. Double check to make sure that none of the multipole axis types are numbers that you deleted from the amoeba.prm file. If there are, change them to the atom type number that represents that group which was not deleted.
Appendix 1: How does amoeba.prm formatting work?
In amoeba.prm, there are atom types and atom classes. Atom types are more specific (there are more types than classes). Multipoles and polarization parameters are defined for each atom type. All other paramteters are defined for entire classes. Below, you can find the format for each section.
Atom Type Definitions:
Type# Class# Element "Description" Atomic# Mass ConnectedAtoms#
Van der Waals Parameters:
Class# EnergyWellDepth ?
Bond Stretching Parameters:
Class#1 Class#2 ForceConstant(k) IdealLength
Angle Bending Parameters:
Class#1 Class#2 Class#3 (k) IdealAngle
Torsional Paramters:
Class#1 Class#2 Class#3 Class#4 (k) CisOrTrans Term#1 (then two more similar terms)
Atomic Multipole Parameters:
Type# XaxisType# ZaxisType# charge
dipoleX Y Z
multipoleXX
YX YZ
ZX ZY ZZ
Dipole Polarizability Parameters:
Type# Polarizability ConnectedAtomType#'s(in same polarization group)
</BODY> </HTML>