About Issues Tutorials Documentation </span>

Example 5: Calculating torsional barriers with relaxation

This workflow calculates the enthalpic barrier of a small alkane.

Author: Aaron Virshup, Autodesk Research
Created on: September 23, 2016
Tags: reaction path, constrained minimization, torsion, enthalpic



In [ ]:

    
import moldesign as mdt
from moldesign import units as u

%matplotlib notebook
from matplotlib.pyplot import *
try: import seaborn  # optional, makes graphs look better
except ImportError: pass

u.default.energy = u.kcalpermol  # use kcal/mol for energy

I. Create and minimize the molecule
II. Select the torsional bond
III. Rigid rotation scan
IV. Relaxed rotation scan
V. Plot the potential energy surfaces
VI. Investigate conformational changes

I. Create and minimize the molecule



In [ ]:

    
mol = mdt.from_smiles('CCCC')
mol.draw()



In [ ]:

    
mol.set_energy_model(mdt.models.GaffSmallMolecule)
mol.energy_model.configure()



In [ ]:

    
minimization = mol.minimize(nsteps=40)
minimization.draw()

II. Select the torsional bond

Next, we use the BondSelector to pick the bond that we'll rotate around.



In [ ]:

    
bs = mdt.widgets.BondSelector(mol)
bs



In [ ]:

    
twist = mdt.DihedralMonitor(bs.selected_bonds[0])
twist

III. Rigid rotation scan

First, we'll perform a simple energy scan, simply by rotating around the bond and calculating the energy at each point.

This gives us only an upper bound on the enthalpic rotation barrier. This is because we keep the molecule rigid, except for the single rotating bond.



In [ ]:

    
angles = np.arange(-150, 210, 5) * u.degree



In [ ]:

    
rigid = mdt.Trajectory(mol)
for angle in angles:
    twist.value = angle
    mol.calculate()
    rigid.new_frame(annotation='angle: %s, energy: %s' % (twist.value.to(u.degrees), mol.potential_energy))



In [ ]:

    
rigid.draw()



In [ ]:

    
figure()
plot(angles, rigid.potential_energy)
xlabel(u'dihedral / º'); ylabel('energy / kcal/mol')
xticks(np.arange(-120,211,30))

IV. Relaxed rotation scan

Next, we'll get the right barrier (up to the accuracy of the energy model).

Here, we'll rotate around the bond, but then perform a constrained minimization at each rotation point. This will allow all other degrees of freedom to relax, thus finding lower energies at each point along the path.

Note: In order to break any spurious symmetries, this loop also adds a little bit of random noise to each structure before performing the minimization.



In [ ]:

    
constraint = twist.constrain()



In [ ]:

    
relaxed = mdt.Trajectory(mol)
for angle in angles:
    print(angle,':')
    
    #add random noise to break symmetry
    mol.positions += np.random.random(mol.positions.shape) * 0.01*u.angstrom
    mol.positions -= mol.center_of_mass
    
    twist.value = angle
    constraint.value = angle
    
    t = mol.minimize(nsteps=100)
    relaxed.new_frame(annotation='angle: %s, energy: %s' % (twist.value.to(u.degrees), mol.potential_energy))



In [ ]:

    
relaxed.draw()

V. Plot the potential energy surfaces

If you plotted butane's rotation around its central bond, you'll see three stable points: two at about ±60º (the gauche conformations), and one at 180º (the anti conformation).

You will likely see a large differences in the energetics of the relaxed and rigid scans; depending on the exact starting conformation, the rigid scan can overestimate the rotation barrier by as much as 5 kcal/mol!



In [ ]:

    
figure()
plot(angles, rigid.potential_energy, label='rigid')
plot(angles, relaxed.potential_energy, label='relaxed')
plot(angles, abs(rigid.potential_energy - relaxed.potential_energy), label='error')
xlabel(u'dihedral / º'); ylabel('energy / kcal/mol'); legend()
xticks(np.arange(-120,211,30))

VI. Investigate conformational changes

This cell illustrates a simple interactive "app" - select the bonds you're interested in, then click the "show_dihedral" button to show their relaxed angles as a function of the central twist during the relaxed scan.



In [ ]:

    
from ipywidgets import interact_manual

bs = mdt.widgets.BondSelector(mol)
def show_dihedral():
    figure()
    for bond in bs.selected_bonds:
        dihemon = mdt.DihedralMonitor(bond)
        plot(angles, dihemon(relaxed).to(u.degrees), label=str(bond))
    legend(); xlabel(u'central twist / º'); ylabel(u'bond twist / º')
interact_manual(show_dihedral)
bs



In [ ]: