</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

``````

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)
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 [ ]:

``````