Intra-molecular rotation

Polyatomic molecules often composed of two entities that are connected through a single bond. These two entities can rotate around the single bond, creating various conformations. Due to quantum mechanical interactions between these two entities, the conformations would have different energies.

This tutorial will show how to scan the dihedral angle in C₂H₆

In order to simplify the dihedral, we need to represent the geometry in Z-Matrix format. Each line in the Z-Matrix start with an atom. Followed by a connecting atom and the bond, then a third atom and angle, and a fourth atom and dihedral, etc.

Z-Matrix input

Below is an example for Z-Matrix The first line would only contain an atom A (since not atom to connect at this point yet) The second line would contain another atom B and the distance_AB. No angle yet, since there are only two atoms upto this point The third line would contain another atom C and the distance_AC and angle_CAB No diheral yet. The fourth line would contains another atom D and distance_AD, angle_DAC, and dihedral_DABC The fifth line would contain another atom E, distance, angle, dihedral

So the totoal degree of freedom is 3N - 6, which is the total degree of freedom 3N minus the translation degree of freedom 3 and rotational degree of freedom 3. This just corresponds to the vibrational degree of freedom of the molecule.

A
B      1      dist_21
C      1      dist_31         2             angle_213
D      1      dist_41         2             angle_214               3                dihedral_3214
....

Psi4 calculation

Below is the input to Psi4.

molecule c2h6 {
  C 
  C 1 R1
  H 1 R2 2 A1
  H 1 R2 2 A1 3 A2
  H 1 R2 2 A1 4 A2
  H 2 R2 1 A1 5 A
  H 2 R2 1 A1 6 A2
  H 2 R2 1 A1 7 A2
}

Avals=range(0,360,2)

table=Table(rows=["A"], cols=["E(SCF)","E(SCS)","E(DFMP2)"])

set basis cc-pvdz

c2h6.R1 = 1.54
c2h6.R2 = 1.0
c2h6.A1 = 109.5
c2h6.A2 = 120.0

for A in Avals:
    c2h6.A = A
    energy('mp2')
    escf = get_variable('SCF TOTAL ENERGY')
    edfmp2 = get_variable('MP2 TOTAL ENERGY')
    escsmp2 = get_variable('SCS-MP2 TOTAL ENERGY')
    table[A] = [escf, escsmp2, edfmp2]

print(table)
relative=table.copy()
relative.absolute_to_relative()
print(relative)

Results

If we plot the escf energy as a function of the dihedral A, we obtain the following PES (plot generated using gnuplot).