Everyday Scientist

comparing gaussian structures

May 23, 2006 at 12:53 am | sam | everyday science, software

Now, there’s probably already a program out there that will compare the optimized structures of quantum-chemistry calculations (e.g. Gaussian). But I couldn’t find one, and I wanted control over how to compare the structures (for instance, I wanted to be able to ignore the errors in the dihedrals for light atoms). So I wrote a MatLab M-file that outputs the absolute differences in the bond lengths, angles, and dihedrals between two structures. It took me a while to get what I wanted out of Molden and make MatLab spit out the right set of numbers, but I finally got it to do what I wanted (I think).

It’s not perfect, I’m pretty happy with it. It just imports the variables from a Molden z-matrix and takes the absolute difference in them, sorting my type. Then I just throw those outputs into SigmaPlot or Excel (ew!) or whatever and do the error analysis I want. For instance, I remove the dihedrals that don’t really matter (some light atoms and some artificially inflated dihedral errors) and just look at the mean and standard deviation for different method/basis set combinations.

I dunno, not essential. But I find it pretty handy. If you want to check it out or try it, I’ll post it below. Let me know if you find errors or if you find ways to make it better.

%Sam Lord: sjlord (at) stanford.edu
%5/17/06
%
%This program compares the z-matrices for multiple QC calculations using
%different levels of theory or different basis sets. It calculates
%differences in the bond lengths, angles, and dihedrals directly from a
%z-matrix from Molden. Be sure to use “molden -A ‘structure'” when looking
%at a file; this prevents Molden from rewriting the z-matrix. Go into the
%Z-matrix Editor in Molden, click “US” and “Gaussian,” then save the
%z-matrix with a good name (e.g. ‘~_in’ or ‘~_out’ for before and after
%”movie,” respectively).
%
%Because I used the crystal structure as the starting z-matrix, I can
%compare the starting structure to the optimized structure (before and
%after pressing “movie” in Molden). This ensures that all the angles and
%dihedrals are defined the same (I was having trouble comparing the
%original crystal structure z-matrix to the calculation .log files’
%z-matrices, because things changed).
%
%Molden defines a dihedral of -190 as +170. In other words, I needed to
%correct these far negative angles to their corresponding positive angles
%(or the difference is 360 too big). So I take any angle over 180 and
%subtract it from 360.

clear all

%Input the initial z-matrix (e.g. ‘~_in’) from Molden (“first”) in for
%reference. Use the z-matrix that Molden outputs (e.g. ‘~_out’) after you
%run “movie” in calculation. These read in the variables at the end of the
%z-matrices:
reference=dlmread(‘file_in’,”,47,1);
calculation=dlmread(‘file_out’,”,47,1);

difference_uncorr=abs(reference-calculation);

%This loop corrects for angles larger than 180:
i=1
while i< (length(difference_uncorr)+1) if difference_uncorr(i)>180
differenceT(i)=360-difference_uncorr(i)
elseif difference_uncorr(i)<=180 differenceT(i)=difference_uncorr(i) end i=i+1 end %Need the transform of the result of the above loop: difference=differenceT' %This separates the three types of values by when in the column they %appear: diff_bondlengths=[difference(1); difference(2); difference(4); difference(7:3:end)] diff_angles=[difference(3); difference(5); difference(8:3:end)] diff_dihedrals=difference(6:3:end) %Now copy the three columns of absolute differences and paste into %SigmaPlot: