How to use Mopac with PCMODEL

Format of this document:
Things you click on or type are in bold. Don't forget the spaces.
Things you substitute are in italics.
For example, if the notes list filename.dat, then you must type in your file name with a ".dat" extension.

To optimize the geometry of a structure:

1. Draw and minimize your structure in PCMODEL. It works best if you draw a small part of your structure, add hydrogens and PI atoms, and minimize; then add to this fragment and minimize again; repeat until the molecule is finished. The atom positions, the bonding pattern, and the geometery (bond angles, dihedrals, stereochemistry) must be close to what you want. Mopac will do a final geometery optimization for you. The atom types do NOT have to be correct. See also Sherry Chemler's atom ordering trick for transition state calculations.
2. WRITEF your structure in mopac format. The file name must be filename.dat. Enter the keywords PM3 CHARGE=## MMOK EF HESS=1 GNORM=0.01 PRECISE NOINTER T=###### where CHARGE=## is the total charge of the molecule (if charge=0, then you don't need this) and T=#### is the maximum number of seconds (in CPU time, not our time) allowed for the calculation (T=86400 is 24 hours). Enter a title for your structure.
3. Open a UNIX shell. Then type:
   cd pcmodel
   ls to see your directory contents.
   jot filename.dat
4. Edit this file so that the atom types are correct (e.g., if you constructed a molecule with Co, and you want this atom to be Si, change Co to Si). See also the geometry optmization example dat file.
5. Save your file and exit jot.
6. Type /remote/server/mopac93/mopac93.exe filename &. Do NOT add the ".dat" extension to the file name.
7. Type top. Make sure that your mopac calculation is using mondo %CPU. If it isn't, try step 6 again. If your calculation is running, you can log out.
8. Wait a while, then log back in.
9. Type cd pcmodel
10. Type more filename.out. If the end of this file lists that the calculation has successfully completed, you may proceed. Otherwise, log out and try again later.
11. Type cp filename.out filename_opt.dat
12. Type jot filename_opt.dat and edit this file so that only the last "PM3..." line, the next two lines, and the atom list following these lines are in this file. Delete everything else, including any lines after the atom list. Save the file and exit. See also the geometry optmization example_opt.dat file.
13. Type pcmod
14. Read the mopac file filename_opt.dat into pcmodel.

To calculate a transition state:

1. Optimize the geometry of the reactant and product following the instructions above. You should then have a file with the reactant coordinates and another file with the product coordinates, both named filename_opt.dat.
2. Open a UNIX shell. To combine the reactant and product into one file, type:
   cd pcmodel
   babel -iai reactant_opt.dat -oac reactant_optc.dat
   babel -iai product_opt.dat -oac product_optc.dat
   cp reactant_optc.dat reactant_product.dat
   cat product_optc.dat >> reactant_product.dat
   
3. Edit reactant_product.dat so that it has ONLY two sets of coordinates at the bottom (with a blank line in-between) and these three lines at the top:
   PM3 CHARGE=## SADDLE XYZ BAR=0.05 T=#####
   your title here
   leave this line blank
   start coordinates here
   where CHARGE=## is the total charge of the molecule (if charge=0, then you don't need this) and T=#### is the maximum number of seconds (in CPU time, not our time) allowed for the calculation (T=270000 is 72 hours). See also the transition state example dat file.
4. Organize the list of atoms in the reactant and the product so that the order is the same. See the Transition State example or Sherry Chemler's atom ordering trick for transition state calculations.

   Sherry Chemler's atom ordering trick for transition state calculations

   Sherry Chemler uses the following instructions to create atom lists for the reactant and product that have the same order:
   1. Draw your reactant in PC-MODEL. Minimize. Save the reactant as a mopac file.
   2. Draw your product in PC-MODEL. Minimize. Add and delete bonds so that the product's bonds are the SAME configuration as the reactant's bonds, but DON'T change the coordinates of the atoms. Then save the product as a mopac file.
   This method works because mopac does NOT care about the bonding arrangement when it does a calculation---it only cares about the distances between atoms to decided where the bonds are. However, PCMODEL uses the bonding arrangement to determine the order of the atom list. If your reactant and product have the same bonding arrangement, PCMODEL will create a reactant and product atom list with the same order.

5. Type /remote/server/mopac93/mopac93.exe reactant_product &. Do NOT add the ".dat" extension to the file name.
6. Type top. Make sure that your mopac calculation is using mondo %CPU. If it isn't, try step 6 again. If your calculation is running, you can log out.
7. Wait a long while, then log back in.
8. Type cd pcmodel
9. Type tail reactant_product.out. If the end of this file lists that the calculation has successfully completed, you may proceed. Otherwise, log out and try again later.
10. Type cp reactant_product.arc reactant_ts.dat
11. Type jot reactant_ts.dat and edit this file so that only the (second) "PM3..." line, the next two lines, and the atom list immediately following these lines are in this file. Delete everything else, including any lines after the first atom list. Save the file and exit. See also the transition state example arc file.
12. Type babel -iac reactant_ts.dat -oai reactant_tsi.dat
13. Type cp reactant_product.arc product_ts.dat
14. Type jot product_ts.dat and edit this file so that only the last "PM3..." line, the next two lines, and the atom list immediately following these lines are in this file. Delete everything else, including any lines after this last atom list. Save the file and exit. See also the transition state example arc file.
15. Type babel -iac product_tsi.dat -oai product_tsi.dat
16. Type pcmod
17. Read the mopac files reactant_ts.dat and product_tsi.dat into pcmodel.
18. To make a "video" of the reaction pathway from the reactant_product.out file:
    1. Make a subdirectory with coordinate files of each step of the movie.
    2. If the coordinate files are "internal" mopac coordinates, type:
       babel -iai file1.dat -op file1.pdb
       If the coordinate files are "cartesian" xyz coordinates, type:
       babel -iac file1.dat -op file1.pdb
       
    3. Start InsightII. Read pdb2arc_dir/reactant_product_1.pdb file into the prgram: Click on molecule/get, select PDB format, and click on your file name.
    4. Select Object/Source_file and source /remote/server/mopac93/PDB_to_Arc.bcl. You also need car2arc.awk in the same directory asn PDB_to_Arc.bcl, but this is set up for you here at IU MolViz Facility.
    5. Select Custom/PDB_to_Arc.bcl (see step B3). The PDB files in pdb2arc_dir will be converted to 1 'mopac.arc' file. pdb2arc_dir will not be deleted, but the files temp*.car and temp*.mdf in pdb2arc_dir will be deleted. You may want to rename 'mopac.arc' in a UNIX shell after it is created.
    6. Select Analysis (under the Biosym icon).
    7. Select Trajectory/Get and read in the archive file. You need to have a 'target object' on the screen before you read in an archive file.
    8. Select Trajectory/Animate to watch your "video". you can also perform other analyses, Such as graphing bond distances as a function of reaction "time".