Calculation of Birefringence in Liquid Crystals
Tutorial based on winning solution for InnoCentive
Challenge
"Calculation of the Intrinsic (Maximum) Birefringence of Molecules"
Updated (February 2015): Now the polarizability components for different types of bonds (Table 2 below) are contained in a separate file (bond_polarizabilities.txt), which must be present in the working directory. The file has the following format:
#Comment lines (the head lines only)
A1, A2, TRIPOS bond order, (alpha_isotropic, alpha_parallel, alpha_perpendicular)
END
For each bond type A1 and A2 are the atomic numbers in increasing order; bond order according to TRIPOS mol2 format (available values are 1 = single; 2 = double; 3 = triple;
am = amide; ar = aromatic; du = dummy; un = unknown; nc = not connected); alpha (polarizability) components are in Å^{3}. The user may now add his own types of bonds, or edit the existing ones. For example, the record for double C=O bond is looking like this:
6,8,2,(2.0,1.333,1.0)
This file now contains polarizability data for CF bonds (thanks to Bruno Cramer, State University of Campinas, Brasil), based on Miller KJ, Additivity Methods in Molecular Polarizability, JACS,1990, 112, p. 85338542
In brief the challenge was formulated as following: predict experimental birefringence (Dn, with 10% error) for pure liquid crystal substances, at the wavelength of 550 nm, using only their molecular structure as input.
It is not possible to apply some form of LorentzLorenz
equation to a single molecule to obtain birefringence value (Dn)
for a simple reason that the polarizability (a) is a molecular
property, but refraction coefficients (n) and Dn are the properties of the condensed phase. And this is
explicitly expressed in LL equation by the term N (number density, a number of
molecules per unit volume).
Therefore, in order to apply LL equation certain assumptions about the phase need to be taken. E.g. in a crystal we know that all the molecules arranged in a certain orientation in relation to each other, and fixed. So, LL equation can produce adequate n for certain crystals.
In the proposed solution the assumption was made that the liquid crystals form nematic phase, meaning that they generally keep their orientation along the director (x axis), but their positions are not fixed and they are able to rotate about their internal x axis:
The first figure is from ref. [1]; the second one is the representation of the nematic phase used in the challenge solution.
Brief Summary
Total 270 variations of five selected methods were tested for the calculation of birefringences (Dn). None of them were capable of direct reproduction of experimental values with satisfactory accuracy. However, many of these methods produced results, which excellently correlate with the experiment (r>0.95). Three best linear correlation functions were selected as a solution for the challenge. One method uses only properly optimized molecular structure (mol2 file); other two require additional calculation of molecular polarizabilities (a) using time dependent HartreeFock (TDHF) method as realized in Firefly (former PC GAMESS) and GAMESS programs. Application of these correlation functions allows for calculation of Dn with error less than 10% for nematic or smectic liquid crystals for diverse molecular structures.
The selected methods are nicknamed “Vuks”, “Direct”, and “Order”. The correlation charts between calculated and experimental Dn are shown in Fig. 1, and the test results are in Table 1.
Figure 1. Correlation graphs for the three methods deemed suitable for the calculation. Circles — data points of five substances used for regression analysis. Crosses — data points for substances used to test the calculation method.
Table 1. The results of calculations for the three selected methods.
Substance ID 
Structure 
Calculated Dn (Error, %) 
Experimental Dn 

“Vuks” 
“Direct” 
“Order” 

(45) 

0.225 
0.101 
0.059 
0.045[1] 
(121) 

0.132 
0.137 
0.128 
0.121[1] 
(125) 

0.305 
0.131 
0.110 
0.125[1] 
(175) 

0.161 
0.178 
0.176 
0.175[1] 
(194) 

0.196 
0.203 
0.201 
0.194[1] 
(197) 

0.185 
0.177 
0.179 
0.197[1] 
(215) 

0.213 
0.215 
0.213 
0.215[1] 
Toolbox required
1. General purpose molecular modeling program, which allows to build the molecule, read and write mol2 files, and align the molecule along x and y coordinate axes. E.g. Chimera or Accelrys Visualizer.
2.Open Babel to prepare the input files and extract optimized structures.
3. Ab initio computational chemistry program: Firefly (former PC GAMESS) or GAMESS. The solution program is designed to use the output files from these packages only.
4. The alignment utility, mol2align.exe.
5. The birefringence calculator, dn_corr.exe.
Steps to calculate birefringence
1. Create the structure of interest in your molecular modeling package and check that it has “reasonable” 3D structure (no exaggerated angles or bonds). If it doesn’t look right, perform a “cheap” structure optimization (e.g. by using molecular mechanics or semiempirical PM3 method available in most modeling software).
As an example I will use structure (194) from Table 1:
Distorted (for some reason) structure:
“Reasonable” structure (if the substance has biphenyl fragment like the one below, make sure that the dihedral angle between the rings is nonzero, otherwise the optimization may end up in incorrect planar configuration):
By convention, the resulted file is saved as 194.mol2 (substance ID is used as the file name).
The example mol2 files for the compounds in Table 1 could be downloaded here: substances.7z.
2. Prepare the input for fine structure optimization using RHF method, 6311+G(d,p) basis set [2] in Firefly. This is done by Babel with the following command:
babel.exe imol2 194.mol2 ogamin 194.inp xf header_opt.inp
Where “header_opt.inp” is the header file for the RHF method above, and is could be obtained via the link.
3. Run the structure optimization with Firefly using the input file (194.inp) produced in step 2. For the sample compounds in Table 1 the time of calculation is about 6–8 hr. on an Intel Core i7 processor with 4 cores.
4. Extract the optimized structure into mol2 file from Firefly output using Babel:
babel.exe igamout Firefly_output_file omol2 194.mol2
5. Very important! Open the optimized mol2 file from step 4 in your modeling program, and align the longest direction of the molecule with the x axis, and “widest” direction with the y axis, so that the molecule will be oriented approximately along the principal axes of the tensor of molecular polarizability, like on the following illustrations taken from [1]:
Save mol2 file, overwriting the old one.
Preferred method is to use the supplied mol2align program using the command
mol2align.exe xxx.mol2
producing xxx_aligned.mol2 file, which needs to be checked e.g. in Chimera package. The properly aligned molecule will have overall molecular plane lying in xyplane, the longest molecular axis aligned with xaxis, and molecular center of mass at the origin.
The program is not tested thoroughly, so on rare occasions, the alignment may fail. In this case one needs to manually rotate the original molecule at random and repeat the procedure. If there will be too many failures with mol2align, please report it to me.
6. After step 5 you already can do the Dn calculation by “Vuks” method using the following command:
dn_corr.exe 194.mol2
The result will be written to “194.mol2.out” file (more on the result format below).
7. To use “Direct” and “Order” methods it is necessary to calculate molecular polarizabilities at 550 nm using TDHF method in Firefly. So, in this step you will prepare the input file for Firefly using optimized and aligned mol2 file from step 5:
babel.exe imol2 194.mol2 ogamin 194.inp xf header_tdhf.inp
Where “header_tdhf.inp” is the header for polarizability calculation at 0 and 550 nm (you may get it via the link).
8. Run Firefly with the input file produced in step 7. Depending on the structure, full calculation may take 2–4 hrs on Intel Core i7, 4 cores run in parallel. The following strings from the Firefly output must be copypasted into 194.out for further calculations:
*****************************************************************
CALCULATION FOR A FREQUENCY OF 2.25425 EV ( 0.08284 A.U.)
WAVELENGTH OF 550.00 NM
*****************************************************************
......................................................................
ALPHA AT 0.08284 A.U.
......................................................................
..... INITIATING DIIS PROCEDURE .....
MAXIMUM RESPONSE VECTOR ERROR= 1.866282E+00
MAXIMUM RESPONSE VECTOR ERROR= 1.233805E+00
MAXIMUM RESPONSE VECTOR ERROR= 1.673674E01
ALPHA XX = 357.96100098340820 UA CONV. TO 4.70E06 IN 17 ITER.
ALPHA YX = 14.22455326214941
ALPHA ZX = 7.30887739164954
..... INITIATING DIIS PROCEDURE .....
MAXIMUM RESPONSE VECTOR ERROR= 3.708951E01
MAXIMUM RESPONSE VECTOR ERROR= 7.088355E06
ALPHA YY = 181.54755787106600 UA CONV. TO 6.47E06 IN 16 ITER.
ALPHA XY = 14.22453535016288
ALPHA ZY = 6.52190907992419
..... INITIATING DIIS PROCEDURE .....
MAXIMUM RESPONSE VECTOR ERROR= 7.319860E02
MAXIMUM RESPONSE VECTOR ERROR= 4.290535E02
ALPHA ZZ = 106.71731105211180 UA CONV. TO 9.86E06 IN 14 ITER.
ALPHA XZ = 7.30887464741015
ALPHA YZ = 6.52190977584003
ISOTROPIC AVERAGE ALPHA = 215.40862330219530 A.U.
The example “.out” files for substances from Table 1 are in the substances.7z file.
If the seeker wishes to use other quantum chemical programs (Gaussian, MOPAC etc.), the results of polarizability calculation must be extracted and put in the “.out” file manually according to the following format (blanks are important!):
ALPHA AT 0.08284 A.U.
ALPHA XX = 357.96100098340820
ALPHA YX = 14.22455326214941
ALPHA ZX = 7.30887739164954
ALPHA YY = 181.54755787106600
ALPHA XY = 14.22453535016288
ALPHA ZY = 6.52190907992419
ALPHA ZZ = 106.71731105211180
ALPHA XZ = 7.30887464741015
ALPHA YZ = 6.52190977584003
9. Calculate birefringence using all three methods using mol2 file from step 5 and .out file from step 8:
dn_corr.exe 194.mol2 194.out
The calculation produces short output file “194.mol2.out”, containing the following information for three methods (comments are in red):
Solution for InnoCentive challenge #9932811 by chern

Birefringence calculation for: 194.mol2
S parameter: 0.910 Used for “Order” method
D parameter: 0.000 Used for “Order” method
Number density (N): 2.414 See the technical details below
Contribution of ny to no: 0.900 See the technical details below
Contribution of nz to no: 0.100

Polarization tensor from structure (nm^3): Derived from mol2 file
0.0271 0.0000 0.0000
0.0000 0.0245 0.0000
0.0000 0.0000 0.0245
###
Vuks calculations
NX: 2.1708 x component of n
NY: 2.0845 y component of n
NZ: 2.0845 z component of n
N_AVERAGE: 2.1132 <n>
NO: 2.0845 n_{┴}
NE: 2.1708 n_{}
DN: 0.0863 Dn calculated
Predicted birefringence from mol2 structure: 0.196 Final predicted Dn based on the correlation, your main result

TDHF calculation at 0.08284 a.u. (550.0 nm)
Polarization tensor from TDHF (nm^3): Derived from TDHF calc.
0.0530 0.0021 0.0011
0.0021 0.0269 0.0010
0.0011 0.0010 0.0158
###
DIRECT calculations
NX: 2.1146
NY: 1.4564
NZ: 1.2534
N_AVERAGE: 1.6082
NO: 1.4361
NE: 2.1146
DN: 0.6784
DN_ERROR: 0.0 %
Predicted birefringence from DIRECT TDHF: 0.203
###
ORDER calculations
NX: 2.1146
NY: 1.4564
NZ: 1.2534
N_AVERAGE: 1.6082
NO: 1.0743
NE: 1.1104
DN: 0.0361
DN_ERROR: 0.0 %
Predicted birefringence from ORDER TDHF: 0.201
Technical details
In the course of algorithm development general approach was to:
1. Obtain molecular polarizabilities (a) for a given molecule;
2. Calculate parallel and perpendicular components of the refractive index (n_{ }and n_{┴} correspondingly).
Only the computation methods selected for the final solution will be described below.
Molecular polarizabilities(a)
Two procedures were tested for the calculation of the molecular polarizability tensor
:
1. The bond increment procedure described in [3, 4], with slight modifications.
Each chemical bond in the molecule is assumed to contribute certain increment to the total molecular polarizability. The increments, which are currently in the program database, are listed in Table 2. This increment depends on the bond’s orientation in relation to the main directions of the polarizability tensor (hence the importance of the step 5 of molecule alignment above), and the final increments were calculated as following:
For each bond i:
(1)
(2)
Where (x,y,z) are the angles formed by the bond i and axes x, y, z correspondingly, and a_{}, a_{┴ }are the bond values from Table 2.
The final molecular polarizability is the sum of all bonds’ increments:
(3)
This calculation requires only properly optimized and aligned molecular structure, but is limited only to the types of bonds listed in Table 2 (if necessary, it could be extended), and very sensitive to the orientation of the molecule and its conformation.
Table 2. Bond types and their polarizability increments.
Bond type 
a_{isotropic, Ǻ}^{3} 
a_{, Ǻ}^{3} 
a_{┴, Ǻ}^{3} 
Reference 
C–H 
0.82 
0.673 
0.6 
[4] 
C–O 
1.46 
0.6 
0.17 
[4] 
C=O 
2.0 
1.333 
1.0 
[3, 4] 
C–C 
0.97 
0.497 
0.25 
[3, 4] 
(ar.) 
2.25 
1.07 
0.48 
[4] 
C=C 
2.9 
1.68 
1.07 
[4] 
C≡C 
2.03 
3.54 
1.27 
[5] 
C–N 
1.38 
0.607 
0.22 
[4] 
C=N 
2.23 
1.429 
0.9 
[4] and extrapolations 
C≡N 
1.97 
3.1 
1.4 
[5] 
C–S 
3.42 
1.97 
1.25 
Extrapolations from [4] 
C=S 
7.57 
4.37 
2.77 
[4] 
N=N 
4.04 
2.343 
1.49 
[4] and extrapolations 
N–H 
0.75 
0.58 
0.84 
[5] 
C–Cl 
4.32 
3.67 
2.6 
[6] 
C–Br 
6.24 
5.04 
3.71 
[6] 
C–I 
9.57 
8.09 
5.77 
[6] 
The last fact could be the main reason why it didn’t work for conformationally labile compounds, like cyclohexane derivatives (45) and (125) from Table 1. For rigid molecules made of aromatic rings conjugated through C=C bonds and heteroatoms (O,S,N) this method works surprisingly well.
2. Calculation of molecular polarizabilities using timedependent HartreeFock method at 550 nm.
This method is implemented in high quality quantum chemical packages Firefly and GAMESS, and described in [7, 8]. It allows for calculation of at a given experimental wavelength (550 nm), and does not appear to be that sensitive to the molecule orientation. The input for this method is the same molecule as for the previous one (see calculation steps above).
The basis set was adopted from ref. [2], where it had been used to calculate nonlinear optical properties of azulenes. As long as it produced acceptable results, I did not investigate any further.
What I found so far most important in practical predictions of Dn is not the accuracy of polarizability, but the consistency of the procedures. By the latter I mean that all computation procedures must be performed in exactly the same manner. No matter how accurate are polarizabilities, the uncertainty of the phase structure will nullify that. Yet, as I mentioned before, none of the methods used are able to reproduce experimental Dn directly. However, many methods produce Dn in excellent correlation with the observed values. The equation below is of the utmost importance for the success of birefringence prediction.
Dn_{experimental }= aDn_{calculated} + b
Where a and b are the linear regression coefficients. For the time being these are just phenomenological coefficients resulting from the regression analysis. They cover any uncertainty in the calculations of polarizabilities, and lack of understanding of molecules interaction in the mesophase.
If one wishes to pursue fundamental investigation of the mesophase optics then, probably, robust DFT calculations will help to get rid of the first uncertainty (polarizabilities will be very accurate). If there’s only need to practically predict the birefringence, then there’s no reason to use more costly methods.
However, in any case, the regression coefficients above must be found separately for any method used. And, to be on the safe side, even for the same method used by different quantum chemistry package (e.g. Gaussian).
Refractive indices (n_{ }and n_{┴})
1. “Direct” calculation.
“Direct” method uses well known LorentzLorenz equation [1, 9, 10] to calculate the components of refractive indices:
(4)
Where n_{i} are refractive indices along x, y, or z axes, and N is the number density (number of absorbing molecules). As long as the main unit of polarizability in calculations is nm^{3}, N is measured as the number of absorbing molecules per nm^{3}.
Two ways of calculating N were tried: N = 1/molecular volume (nm^{3}), and N calculated from the molecular mass assuming that density of most liquid crystals is around 1 g/ml: N = 602/M_{r}. The last formula demonstrated better results , so it was used throughout to determine N.
Three calculated refractive indices were reduced to the main two n_{ }and n_{┴}_{ }by the following approach:
n_{ }= n_{x}_{ }(5)
_{ }
(6)
Where f_{n,y} is the fraction contributed by n_{y}, and 1 f_{n,y} is the fraction contributed by n_{z}. The values of f_{n,y} in the range from 0.1 to 0.9 were tested for all calculation methods, but they had little or no influence on the final correlation functions. In simple terms, the exact contribution of n_{y} or n_{z} to n_{┴} is not important for Dn.
2. “Vuks” calculation.
It uses a variation of LorentzLorenz equation first suggested by Vuks [3, 911], assuming that the internal field in the condensed phase is isotropic:
(7)
Here n_{i}^{2}+2 is replaced by an average of squares of all n components. Therefore, in order to find n_{i} from Vuks equation it is necessary to solve a system of three equations.
All other details for calculation of N, n_{ }and n_{┴ }are the same.
3. “Order” calculation.
The method is adopted from [1], pg. 116 onwards. A modified form of LorentzLorenz equation is used:
(8)
For this equation a_{} and a_{┴} are calculated using liquid crystals order parameters S and D:
(9)
(10)
The order parameters S and D reflect the anisotropy of liquid crystals, decreasing with the temperature, and finally zeroing out at the clearing point. For the present solution a consensus value of S = 0.91at 25ºC has been chosen [12], while D was ignored and left out to be zero [1].
The calculations for N were done as in “Direct” method.
Regression analysis
As it has been noted in the summary, none of the methods was capable of reproducing the experimental values of birefringence with acceptable accuracy. However, when the experimental Dn were plotted against calculated ones, many methods produced very clear linear correlation plots (Fig. 1) with regression coefficient r>0.95.
In the linear regression analysis, the data of only five compounds from Table 1 were included. The other compounds were used to test the final selected algorithms.
Below the final list of the three methods selected for the solution is given, together with corresponding linear regression equations.
1. “Vuks” method (adopted from [3])
a: calculated from bond increments;
n: calculated according to Vuks [3, 11].
Regression equation: Dn_{experimental }= 1.9392Dn_{calculated} + 0.0282; r = 0.983
Comments: Requires only properly optimized structure (mol2 file), works fine for rigid molecules (error <10%, see Table 1); doesn’t work for flexible molecules like cyclohexane derivatives (error>100%); sensitive to molecule’s orientation; suitable only for limited set of chemical bonds (Table 2).
2. “Direct” method
a: calculated from TDHF;
n: calculated according to LorentzLorenz equation.
Regression equation: Dn_{experimental }= 0.2249Dn_{calculated} + 0.0505; r = 0.999
Comments: Requires properly optimized structure, and the output from TDHF calculation; works well for all but very flexible molecules like dicyclohexane derivatives (typical error <10%, see Table 1).
3. “Order” method
a: calculated from TDHF;
n: calculated using order parameter S.
Regression equation: Dn_{experimental }= 5.8643Dn_{calculated} 0.0112; r = 0.994
Comments: Could be recommended as the best of all three methods. Requires properly optimized structure, and the output from TDHF calculation; worked for all tested compounds (error <10%, with one exception, when error is 12%, see Table 1). The large error for compound (45) of 31% is due to small experimental value of Dn (0.045), otherwise the absolute error in this case is comparable to other substances.
Supplement files
The following files are available to aid the calculations above:
dn_corr.exe Compiled birefringence calculator
mol2align.exe Compiled alignment utility
header_opt.inp Header for structure optimization with
Firefly or GAMESS (RHF 6311+G(d,p))
header_tdhf.inp Header for TDHF polarizabilities
calculation with Firefly or GAMESS
SUBSTANCE_ID\
id.mol2 Optimized and aligned mol2 structure
id.mol2.inp Input file for TDHF calculation
id.out Required excerpt from TDHF output file
id.mol2.out Output of “dn_corr id.mol2 id.out”
command
References
1. Demus, D., Physical properties of liquid crystals1999, Weinheim ; New York: WileyVCH. 503 pp.
5. Lippincott, E. and J. Stutman, Polarizabilities from dFunction Potentials. J. Phys. Chem., 1964. 68(10): p. 29262940.