Intro to Molecular Modeling- Sketching, MM, MD Annealing and QM


1) The instructor will provide you with student accounts and passwords. These accounts will remain active for two weeks in case you would like to try additional things on your own.


2) Log onto the SGI workstations with your account and password


3) Materials for the training exercises will be available in the directory "MGMcourse".  To access this directory, type the following command at the command prompt:

·  cd MGMcourse

and hit [enter]


4) For these exercises, you will be using the SYBYL suite of software.  To run SYBYL, type:

·  sybyl

and hit [enter]



In this exercise you will be asked to draw and save the following molecule:






1) To begin, locate the {Build/Edit} button on the top bar of the SYBYL interface and click it.


2) On the pull-down menu, you will see a button for {Sketch Molecule...} -- click that too.  Note:  in the future, I will use the following shorthand notation:

*  {Build/Edit} >> {Sketch Molecule...}

to simplify the description of steps like 1) and 2)


3) SYBYL can handle as many as 10 different molecules in memory at any given time.  It does this by storing the different molecules in bins called "Molecular areas" numbered M1 through M10.  SYBYL is asking you which Molecular area you would like to use for your first molecule.  M1 is highlighted by default.  Simply click {OK} to select M1.


4) You now have the three widgets used for building molecules.  On the middle widget, click the {Draw} button.  You are now ready to being drawing.

Here are some hints for drawing:

*  In "Draw" mode, every time you click on the black screen, you will draw an atom

*  If the atom is enclosed in a small yellow square, it means that the next atom you draw will form a bond to the atom with the yellow square

*  If you do not want your next atom to be bound to the atom with the yellow square, click inside the yellow square and it will go away

*  If you want your next atom to be bound to an atom that does not have a yellow square around it, click on the atom you want to bind to and that will give it a yellow square

*  If you want to form a bond to two already existing atoms, simply click the two atoms in succession

*  If you want to form a multiple bond between two atoms, click on the two atoms in succession again


5) Begin inputting atoms where you would like to see them in your molecule.


*  When you first input your atoms, they will all be Carbon by default.  Leave them that way for now -- it will be easier to change them from Carbon to other elements later than it will be while you're sketching.

*  You should not bother drawing in any atoms that will correspond to Hydrogens when sketching.  Hydrogens can be added automatically later.


6) Sketching molecules can be tricky and many people end up with a few bonds that they did not intend to sketch.  If you would like to delete such spurious bonds:

·  Go to the middle drawing widget and click {REMOVE_BOND}

·  Click on the two atoms whose bond you wanted to delete


*  As soon as you click {REMOVE_BOND}, you are no longer in "Draw" mode -- you are in "Delete Bond" mode and will keep on deleting bonds with every click until you get back to "Draw" mode.

*  To get back to "Draw" mode, go to the middle drawing widget and click {Draw}


7) If you have any extra atoms that you don't want:

·  Go to the middle drawing widget and click {REMOVE_ATOM}

·  click on the atoms that you would like to delete

·  Don't forget to go back to "Draw" mode when you are done deleting!  Otherwise you will continue to delete atoms


8) Once you have placed atoms corresponding to all of the heavy atoms in flavopiridol, and have the right bonds between all of them, you are  ready to give them the right atomic type.

·  Flavopiridol has 3 Oxygen atoms.  In order to change atoms from Carbon to Oxygen, click on the {O} in the bottom widget and then click on each of the atoms you would like to change.

·  Do the same for the Chlorine and Nitrogen atoms.


9) You are now ready to add H's.  To do this, simply go to the middle widget and click on {ADDH}.


*  If you have made any mistakes in specifying your bonds (e.g., drawing a single bond where you wanted a double), you will get the wrong number of Hydrogens.  Check your structure carefully.


10) There are numerous other molecular drawing options in the sketch tool.  These are described in detail in the online documentation that comes with our Tripos license.  To access this online documentation, minimize your SYBYL window, click on the Netscape icon, and find the Bookmark entitled {Sketcher Menu Items}. 


11) For now, however, let's exit the molecular sketcher by clicking the {Exit} button on the middle widget.


12) You now have your molecule alone on the black screen.  You can look it over carefully by manipulating it in the screen using your mouse.      Specifically:

·  If you hold down the right mouse button and move your mouse, you will rotate the molecule

·  If you hold down the middle mouse button and move your mouse, you will translate molecule

·  If you hold down both the middle and right buttons you will focus in, or pan away from the molecule


13) If you notice any serious chemical errors (incorrect atoms or bonds) in your molecule, go back to the sketcher as follows:

·  {Build/Edit}  >>  {Sketch Molecule}  >>  select M1  >>  {OK}    and make the relevant corrections.


14) In step 12 you probably noticed that the molecule does not have the proper three dimensional structure.  To correct this, you are going to need to do a geometry optimization as part of the next tutorial.



In this exercise you will extract the taxol molecule from a crystal structure for it docked in the tubulin protein.  You will then minimize the structure and try to predict how an Na+ ion binds to the molecule.


1) Read in 1JFF.pdb and extract taxol:

·  File  >>  Read...  >>  1JFF.pdb  >>  {OK}  >>  Center the Molecule = No

·  Build/Edit  >>  Extract...  >>  Substructures...  >>  scroll down to TA1601 and select it  >>  {OK}  >>  {OK}  >>  select M2 as the molecular area  >>  {OK}


2) Delete 1JFF:

·  Build/Edit  >>  Zap (Delete) Molecule  >>  select M1  >>  {OK}


3) Color taxol by atom types:

·  View  >>  Color  >>  By Atom Type

·  View  >>  Label  >>  Atom Id...  >>  All  >>  {OK}


4) Look at atoms 7-12 (should be a phenyl) and atoms 44,N43 (should be part of an

amino group) and observe the erroneous bond assignments.  In order to properly

represent the molecules, we are going to need to change these bonds by changing

the specified atom types:

·  Build/Edit  >>  Modify  >>  Atom...  >>  Type  >>  {OK}

·  Click on atoms 7-12,N43,44  >>  {OK}

·  As the program cycles through the different atoms selected, specify 7-12 to be type "", atom N43 to be "N.3", and 44 to be "C.2"


5) Add Hydrogens:

·  Build/Edit  >>  Add  >>  Hydrogens


6) Save starting structure:

·  File  >>  Save As...

·  Specify "File:" as "taxol"  >>  leave "Format:" as "MOL2"  >>  {Save}


7) Optimize the structure:

·  Compute  >>  Minimize...

·  change value for "Max Iterations:" from 100 to 2000

·  For "Energy Setup:", click {Modify...}

·  For "Charges:" change the value from "None" to "Gasteiger-Huckel"

·  Change the "NB Cutoff:" from 8.0 to 15.0  >>  {OK}

·  On the "Minimize" widget, change the "Job Name:" value to "taxol_opt"

·  {OK}


8) Save the structure as taxol_opt.mol2:

·  File  >>  Save As...

·  Specify "File:" as "taxol_opt"  >>  leave "Format:" as "MOL2"  >>  {Save}


9) Add a Na to the system:

·  Build/Edit  >>  Sketch Molecule...  >>  choose "M1:"  >>  {OK}

·  Click on the screen (near to, or inside, the molecule) to add an atom

·  On the elements widget click on "Na" and then click on the atom you just added -- this turns it into a Na atom

·  Click "End Select" on the top widget


10) Change the Na atom into a Na+ ion:

·  Build/Edit  >>  Modify  >>  Atom...  >>  CHARGE  >>  {OK}

·  Click on the Na you just added  >>  {OK}

·  Change the charge from "0.000" to "1.000"  >>  {OK}


11) Reoptimize the structure:

·  Compute  >>  Minimize...

·  Leave all settings the same except "Job Name:" (e.g., try "taxol+na1_opt")

·  {OK}

·  Record the energy value from the blue console window (the line that reads "Optimization completed, Total Energy : .....")


12) Move the Na+ to another location and reoptimize:

·  Build/Edit  >>  Sketch Molecule...  >>  choose "M1:"  >>  {OK}

·  Click on the "Move" button on the middle widget

·  Click on the Na atom on the screen

·  Click on a spot somewhere else in, or around, the molecule

·  Click "End Select" on the top widget

·  Compute  >>  Minimize...

·  Leave all settings the same except "Job Name:" ("taxol+na2_opt")  >>  {OK}

·  Record the total energy


13) Compare your two optimized structures:

·  File  >>  Read...  >>  choose "taxol+na1_opt"  >>  choose  "m2:"  >>  {OK}


*  You should find that your two "minimized" structures are not actually the same structure and have different "minimum" energies.  What does this tell you about the molecular mechanics "minimization" algorithm?  If you wish, you likely have time to minimize a third taxol+Na structure to compare with.



The lesson in the last exercise is that in complex molecular systems, there are often many "local minimum" structures -- metastable conformations in which a molecular mechanics algorithm can get stuck, instead of finding the absolute minimum structure.  The reason that molecular mechanics gets stuck in this way is that the algorithm accounts for only potential energy (not kinetic), and thus can't climb out of any local potential wells in order to find other lower energy configurations.

Molecular dynamics (MD) is a mechanism for adding kinetic energy to a system, by giving atoms velocities that obey statistical thermodynamic averages for a given temperature.  Although MD is mainly intended as a means for studying temperature- and time-dependent effects, it is also quite useful for sampling conformational space in search of lower energy minima.  In this exercise, you will use MD to refine the MM results obtained in the previous exercise.


1) Clear your workspace:

·  Build/Edit  >>  Zap (Delete) Molecule

·  if you have only one molecule on screen, it will disappear

·  If you had multiple molecules to begin with, a "Molecular Expression" widget will appear: click {All}  >>  {OK}


2) Choose and load your starting structure:

·  Refer to your notes from the previous exercise and identify the lowest energy minimized structure that you had been able to generate

·  File  >>  Read...  >>  (choose the mol2 file corresponding to your lowest energy structure)  >>  {OK}  >>  Center the Molecule = No


3) Perform a MD run:

·  Compute  >>  Dynamics  >>  Setup Dynamics...

·  click the {Modify...} button on the top right of the "Dynamics" widget

·  Change "Length:" to 5000

·  Change "Ensemble:" to "NTP"

·  change "Pressure:" to 1 atm

·  change "Initial Velocity:" to "Boltzmann"

·  On the "Dynamics" widget, change "Job Name:" to "taxol+Na_dyn"  >>  {OK}


*  As you run the MD job, you'll see that it samples lots of relatively high energy conformations.  This is natural.  Comparing MD to MM is like comparing how a molecule behaves at room temperature with how it behaves at absolute zero!  However, most instantaneous conformations are higher in energy, MD tends to statistically favor low energy conformations, so the average structure over a longer time intervals tends to have low energy attributes.


4) Average the MD structure:

·  Analyze  >>  Dynamics  >>  Average Structure...

·  Click on your MD history file ("taxol+Na_dyn.his")  >>  {OK}

·  Click {OK} for each of the two analysis prompts

·  The resulting structure on your screen will be the MD-average structure


5) Reoptimize the structure:

·  Compute  >>  Minimize...

·  leave all settings the same except "Job Name:" ("taxol+na_dyn_opt")

·  {OK}

·  Record the energy value from the blue console window


*  How does the final minimized energy compare with those you obtained in the previous exercise?  Most of the time, you should get a lower energy than you got from just straight MM minimization.



This exercise uses quantum chemistry methods to optimize two different molecules (propenol and propanal) and to determine the transition state structure between the two.


1) Bring up SYBYL and sketch the following molecule (propenol):


2) Minimize the molecule with all of the default settings, except for the maximum number of iterations (use 200 instead of 100).  Take note of the final energy and save the final structure as "propenol.mol2".

3) Generate a Gaussian input deck as follows:

·  {Compute}  >>  {Interfaces}  >>  {Gaussian...}

·  Change "Geometry Optimization:" from "None" to "Full"

·  Type something in the "Comment:" box (e.g., "Propenol Geometry optimization")

·  Enter a "Job Name:"  (e.g., propenol)  >>  {OK}

4) The console window should say that it has written the input file "propenol.dat".  You aren’t going to run the job yet, so on the "Option" widget that pops up, click {End}.

5) Sketch the following molecule (propanal):



*  You may find it easier just to rerun the sketcher on the propenol you sketched in step 1.

6) Minimize the propanal and then generate a Gaussian input deck in the same way as Step 3, except for substituting "propanal" for "propenol" for the comment and job name.  Again, note the final energy and save the structure as "propanal.mol2".

7) Shrink your SYBYL window down to icon size (i.e., click on the grey dot on the top right of the window), open up a new UNIX shell (Toolchest  >>  Desktop  >>  Open UNIX Shell) and cd to MGMcourse.

8) Rename propenol.dat to "" and propanal.dat to "".

9) ssh to and cd to the directory ./SGI/MGMcourse:

10) Submit the propenol job by typing:

·  g98s propenol

11) The screen will give you a readout that looks something like:
    Current processor usage is as follows:
    ku1           up      2:06,     2 users,  load 0.02, 0.40, 0.60
    ku2           up      2:03,     0 users,  load 0.00, 0.02, 0.23
    ku3           up      2:03,     0 users,  load 1.96, 1.80, 0.04
    ku4           up      2:03,     0 users,  load 1.95, 1,81, 0.04
    ku5           up      2:03,     0 users,  load 0.00, 0.00, 0.02
    ku6           up      2:03,     0 users,  load 0.00, 0.00, 0.02
    ku7           up      2:03,     0 users,  load 0.00, 0.00, 0.02
    ku8           up      2:03,     0 users,  load 0.86, 0.34, 0.17
    Please specify a low-load processor for your job:

Look at the three last columns which are the processor "load" columns (averaged over the last minute, last 5 minutes and the last hour respectively) choose a processor with a small value for load over the past minute (e.g., in the case above ku2, ku5, ku6 and ku7 have all been inactive), type in this processor and hit the "Enter" key.

12) In a similar fashion, submit the propanal job.

13) The jobs will take several minutes to complete, at which point the output files "propenol.log" and "propanal.log" will appear on your directory.  You can run the command

·  ls *.log

every once in a while to spot them as they arrive.

14) Once the output files are there you can view their contents with the command:

·  cat propanal.log

(and likewise for propenol).  First of all, check to see that the job ran correctly: the very last line of a successful run should read "Normal termination of Gaussian 98."  Next peruse and take note of some of the other output information you get from your run, including multipole moments, charges, spatial extent (related to molecular volume), orbital information, and the final optimized atomic coordinates.

15) The total energy is your calculation is somewhat buried in your output and can be hard to find.  To see this more easily, use the command:

·  grep "SCF Done:" propanal.log

and you will get back a bunch of numbers which reflect the total energies of the molecule at different steps during the geometry optimization.  The last value in this list will be your final total energy.  These energies are expressed in hartrees, and can be converted to kcal/mol (1 hartree = 627.5095 kcal/mol).

16) Analyze the output from both propenol.log and propanal.log, and compare the relative energies.  Note that the total energies will bear no resemblance to those you determined in steps 2 & 6, however the relative energies can be compared.  From your SCF results, compute the heat of reaction for propenol reacting to form propanal and convert this energy from hartrees to kcal/mol.  Compare this with the heat of reaction you would predict from the simple molecular mechanics calculations (steps 2,6).  How do they differ?  Why do you think they differ?

17) Next you're going to find the transition state for the reaction of propenol forming propanal.  To do this, you will need to create a file that contains both the propenol and propanal starting geometries, looking like:
    # HF/3-21G OPT=(QST2)

    Geometry Optimization of Propenol

    0 1
       6      -8.0096951       0.6920720      -0.0557060
       6      -6.9229841       1.7384430      -0.0069330
       6      -5.6248922       1.4333850       0.1515210
       8      -4.6687498       2.3735509       0.1935450
       1      -8.7301207       0.8777330       0.7552720
       1      -7.6163831      -0.3299420       0.0554500
       1      -8.5367546       0.7593480      -1.0194790
       1      -7.2453079       2.7761090      -0.1082960
       1      -5.2990088       0.3985550       0.2527640
       1      -4.9636469       3.2736850       0.1049300

    Geometry Optimization of Propanol

    0 1
       6      -9.8537302       1.3308721      -0.5413280
       6      -8.9657078       2.4280190       0.0967940
       6      -7.5428600       2.2972870      -0.3834240
       8      -7.0096159       3.1974399      -1.0120389
       1     -10.8927994       1.4302990      -0.1895590
       1      -9.4835052       0.3305460      -0.2669700
       1      -9.8417864       1.4244860      -1.6386230
       1      -9.3647881       3.4220450      -0.1645020
       1      -6.9730101       1.3921860      -0.1740790
       1      -8.9807167       2.3281460       1.1940880

To generate this file carry out the following commands:

·  cp

·  cat >>

Note that the ">>" in the above command causes the contents of the file "" to be written to the end of the file "" instead of to the screen.

18) Once you have done this, you will need to edit the input deck (i.e., vi and make the following changes:

·  Line 1:  change "%chk=propenol.chk" to "%chk=transition.chk"

·  Line 2:  change "# HF/3-21G OPT" to "# HF/3-21G OPT=(QST2)"

·  Line 18:  delete entirely (i.e., use the "dd" command)

·  Line 19:  delete entirely

·  Line 20:  delete entirely

Check your file "" to ensure it looks basically the same as the template shown in step 19.  Minor variations in atomic coordinates or comment lines are fine.

19) Once your input deck is finished, save it and quit from vi.

20) Submit the transition state calculation as a Gaussian as follows:  "g98s transition".

21) Analyze the resulting output in the same manner as you did in step 16.  Compute the energy difference between your transition state structure and the final propenol reaction.  This is the activation barrier for the propenol ® propanal reaction.

22) Bring up your SYBYL window again, zap your molecular areas, and read in the following files:

·  Read propenol.mol2 into "M1:"

·  Reread propenol.mol2 into "M2:"

·  Read propanal.mol2 into "M3:"

23) Import your Gaussian output files for analysis as follows:

·  {Analyze}  >>  {Interfaces}  >>  {Gaussian...}

·  For "Gaussian output file:" type in "propenol.log"

·  For "Molecule:" select "m1:"

·  Click the little button next to "Charges".  If nothing happens, click it again, and the quantum chemical charges for each atom in the propenol structure should appear.

·  Click the button next to "Energies" and the text window should give you details such as the total SCF energy, orbital energies, and the dipole moment.

·  Click the button for "Geometry" and after a few seconds, the molecule should be redrawn with the geometry corresponding to the quantum simulation.

Note:  don't click the {OK} button until you're ready to leave the Gaussian analysis widget -- this will close it prematurely.

24) Do the same for your transition state structure by importing "transition.log" into molecule "m2:", and then importing “propanal.log" into "m3:"

25) Compare the relative charges on the atoms of the different structures.  If you have trouble remembering which structure corresponds to which molecule, you may wish to color each of the three structures a different color with the:

·  {View}  >>  {Color}  >>  {Atoms...}

command (remember to specify which molecule you want to color and specify "All" atoms for that molecule).

26) Compare the dipole moments for the three structures.  Why do you think the propanal has the largest dipole moment?  (Hint: in this case it related to the relative positions of different atoms and not due to changes in atomic charges).

David Johnson
Access to more than 40 computational chemistry software programs and databases
High-performance computational tools accelerate drug discovery and minimize costs
Analyze complex, multidimensional data sets to rapidly generate biological insights
KU Today