Performing DFT Calculations on Polymers using Periodic Boundary Conditions

We will use DFT calculations to obtain properties of the polymer polyethylene.

Density (0.898 – 0.962 g/cm³)

Bulk Modulus for medium density polyethylene (0.14 GPa). This is for a polyethylene that has a density of 0.93 g/cm³.

Setting up the DFT calculation

Starting with a periodic assembly as that shown for polyethylene below you will want to calculate the energy of extended system using DMol3. You will need to copy the template input file and rename it so that the root has the same name as the car file. For example, to run the polyethylene_pbc1 job you will want to have the files:

polyethylene_pbc1.car

polyethylene_pbc1.input

The input file contains a number of options. The fact that the calculation is a PBC calculation is not specific in the input file, but rather in the PBC=ON statement in the car file. All other information regarding the type of calculation etc. is to be found in the input file. We will run with the following.

Calculation type: optimize

Note that the optimization does not change the PBC conditions, but permits the nuclei to rearrange themselves within the periodic box in order to find the minimum energy. A geometry optimization is an essential step if one is to determine whether the geometry imposed by the artificial procedure used to obtain the periodic box limits is in any way valid for true polyethylene.

Calculation type: frequency

A frequency calculation should only be submitted on a structure that has been fully optimized. The vibrational frequencies are calculated by a numerical procedure.

Calculation type: energy

If you have an optimized structure and you wish to calculate the absorption spectrum, molecular orbitals, charges, dipoles or other parameters you may choose an energy calculation. No further optimization will be performed, but a set of SCF coefficients (that minimize the energy according to the variational principle) will be determined and the appropriate selected properties will be output in the .outmol file.

The experimental structure should be run using a geometry optimization. The optimization will permit the atoms to relax in the structure and to find their optimal geometries. For the calculation of elastic constants there are two possibilities for the calculation type:

energy = using the optimized geometry the energy will be calculated formed lattice structures with no rearrangement of the atoms in the lattice. In other words if the lattice (and all the atoms in it) are expanded along the z-direction by 2% the energy of this static structure is calculated.

optimize = reoptimize the structure for each deformed lattice. This option permits the molecules to rearrange in the lattice. This is the more realistic option, however, as shown in the section on elastic constants we can learn about the structure by comparing these two calculation types.

Submitting the Job

Once you have modified the input file with the appropriate parameters, you will need to create a .job file so that you can launch the calculation in the queue system of the machine.

Examining the Output

The output is found in the .outmol file. You can examine the output using the vi editor. Alternatively if you wish to extract information from a group of .outmol files the grep UNIX command is very useful. For example, if you have a group of .outmol files that have been generated for calculated elastic constants and you wish to retrieve the energies from those files you will want to find the energy on the line that contains the phrase "binding energy". To obtain this you can use the following command:

Ø grep binding *.outmol

When you enter this command you will see output to the screen all of the occurrences of the work "binding" in all of the .outmol files in the directory. To write these to a file you simply use the redirect.

Ø grep binding *.outmol > pe.energy

Now the file pe.energy contains all of the energies that were output to the monitor in the first grep command. You can examine them by typing:

Ø vi pe.energy

You will notice that there may be several occurrences of the word "binding". The reason for this is that we optimized the geometry in each modified unit cell. Therefore, a binding energy was calculated for each geometry. We want only the final binding energy. You can use the "dd" command in vi to delete a row. Simply delete all of the rows except the final one for each file. In the end your pe.energy file will contain the binding energies of each of the distorted geometries. These can be plotted in Excel or another plotting program to permit you to examine the output. These plots can be used to determine the minimum energy and the energy U(0) of distortion U(D) along x, y, or z for calculation of the elastic constants.