FindQRC
Script to run the program
#!/bin/csh
if (-e $1) then
java FindQRC $*
else
echo This class reads a minimisation file produced via FreqQRC and
echo writes out a QRC in filename.qrc
echo
echo If given two arguments, it is assumed that the second is an .fdt
echo file from FreqQRC, and this should be used for the first structure
echo and energy in the sequence
echo
endif
Java code
import java.io.*;
class FindQRC {
/**
*
* Copyright J M Goodman, University of Cambridge, 2003
*
*
* This class reads a Jaguar output file or a GAMESS output file
* extracts the geometries and energies and prints out the
* energies and the distance between successive steps of the
* minimisation in mass weighted coordinates (Bohr (amu)0.5)
*
* If a second argument is given, it will be assumed that this is
* a .fdt file created by FreqQRC. The geometry and energy in this
* file will be used as the starting point, instead of the first
* structure in the minimisation procedure.
*
* These figures can be pasted into Excel to form a quick
* alternative to a full IRC calculation. It will not follow
* the reaction path so precisely but it will identify starting
* materials and products more rapidly.
*
* The output is recorded in a file called INPUTNAME.qrc
*
* @author J M Goodman
*/
public static void main(String[] args) {
/**
* The names of the input file and the output file
*/
String mname, outName;
if (args[0].lastIndexOf(".") < 0) {
mname = args[0]+".out";
outName = args[0]+".qrc";
} else {
mname = args[0];
outName = args[0].substring(0,args[0].lastIndexOf("."))+".qrc";
}
boolean fdt = false;
String fdtFile="fdt";
if (args.length > 1) {
// Second argument will be a .fdt file
// so must read initial energy and geometry from there
fdt = true;
fdtFile = args[1];
System.out.println("Found fdt file: "+fdtFile);
}
String fileType = "Jaguar";
String writeString;
String currLine="XXXX";
/**
* The number of atoms in the structure
*/
int niatms = 0;
BufferedReader buffRead;
BufferedWriter buffWrite;
/**
* The list of atom types in the molecule
*/
String[] atomtypes;
/**
* The list of atomic weights
*/
double[] atomweight;
/**
* The current coordinates of the atoms
*/
double[] x;
double[] y;
double[] z;
double energy = 0.0;
boolean lasttime = false;
/**
* ox, oy and oz record the original coordinates of the molecule
*/
double[] ox;
double[] oy;
double[] oz;
double oenergy = 0.0;
/**
* lx, ly and lz record the coordinates of the last structure
* before each step
*/
double[] lx;
double[] ly;
double[] lz;
double lenergy = 0.0;
double osos = 0.0;
double lsos = 0.0;
double totlsos = 0.0;
double quatosos = 0.0;
double quatlsos = 0.0;
double totquatlsos = 0.0;
System.out.println("(c) J M Goodman, 2003");
System.out.println("Cambridge University");
System.out.println("All rights reserved");
System.out.println("");
System.out.println("File name: "+mname);
try {
buffRead = new BufferedReader(new FileReader(new File(mname)));
for (int i = 0; i < 21 ; i++) {
currLine = buffRead.readLine();
while (currLine.indexOf("setenv") >= 0) currLine = buffRead.readLine();
if (currLine.indexOf("Jaguar") > -1) fileType = "Jaguar";
if (currLine.indexOf("GAMESS") > -1) fileType = "GAMESS";
}
// First read all through the input files, determining whether it is a Jaguar
// files or a GAMESS file, and counting the atoms
niatms = 0;
if (fileType.equalsIgnoreCase("Jaguar")) {
while (currLine.indexOf("Input geometry") < 0) currLine = buffRead.readLine();
currLine = buffRead.readLine();
currLine = buffRead.readLine();
while (currLine.length() > 10) {
currLine = buffRead.readLine();
niatms++;
}
niatms--;
} else if (fileType.equalsIgnoreCase("GAMESS")) {
while (currLine.indexOf("COORDINATES (BOHR)") < 0) currLine = buffRead.readLine();
currLine = buffRead.readLine();
while (currLine.length() > 10) {
currLine = buffRead.readLine();
niatms++;
}
niatms--;
while (currLine.indexOf("LOCATED *****") < 0) currLine = buffRead.readLine();
}
buffRead.close();
} catch(Exception e) {;}
// Now we know how many atoms there are, set up arrays
atomtypes = new String[niatms];
atomweight = new double[niatms];
x = new double[niatms];
y = new double[niatms];
z = new double[niatms];
ox = new double[niatms];
oy = new double[niatms];
oz = new double[niatms];
lx = new double[niatms];
ly = new double[niatms];
lz = new double[niatms];
// Read through file again, picking up energies and coordinates
try {
if (fdt) {
// System.out.println("fdt 1: ");
buffRead = new BufferedReader(new FileReader(new File(fdtFile)));
currLine = buffRead.readLine();
// System.out.println(currLine);
// System.out.println("01234567890123456789012345678901234567890");
oenergy = Double.valueOf(currLine.substring(29).trim()).doubleValue();
for (int i = 0; i < niatms; i++) {
currLine = buffRead.readLine();
// System.out.println(currLine);
if (fileType.equalsIgnoreCase("Jaguar")) {
atomtypes[i] = currLine.substring(0,9);
atomweight[i] = atomtype2atomweight(atomtypes[i]);
ox[i] = Double.valueOf(currLine.substring(9,19).trim()).doubleValue();
oy[i] = Double.valueOf(currLine.substring(20,32).trim()).doubleValue();
oz[i] = Double.valueOf(currLine.substring(33).trim()).doubleValue();
} else {
atomtypes[i] = currLine.substring(0,9);
atomweight[i] = atomtype2atomweight(atomtypes[i]);
ox[i] = Double.valueOf(currLine.substring(17,28).trim()).doubleValue();
oy[i] = Double.valueOf(currLine.substring(29,40).trim()).doubleValue();
oz[i] = Double.valueOf(currLine.substring(41).trim()).doubleValue();
}
}
buffRead.close();
}
buffRead = new BufferedReader(new FileReader(new File(mname)));
buffWrite = new BufferedWriter(new FileWriter(new File(outName)));
buffWrite.write("Pseudo-IRC analysis");
buffWrite.newLine();
buffWrite.newLine();
currLine = "XXXX";
if (fileType.equalsIgnoreCase("Jaguar")) {
while (currLine.indexOf("Input geometry") < 0) currLine = buffRead.readLine();
currLine = buffRead.readLine();
currLine = buffRead.readLine();
for (int i = 0; i < niatms; i++) {
currLine = buffRead.readLine();
atomtypes[i] = currLine.substring(0,9);
atomweight[i] = atomtype2atomweight(atomtypes[i]);
lx[i] = Double.valueOf(currLine.substring(9,25).trim()).doubleValue();
ly[i] = Double.valueOf(currLine.substring(25,43).trim()).doubleValue();
lz[i] = Double.valueOf(currLine.substring(43).trim()).doubleValue();
if (!fdt) {
ox[i] = lx[i];
oy[i] = ly[i];
oz[i] = lz[i];
}
}
} else if (fileType.equalsIgnoreCase("GAMESS")) {
while (currLine.indexOf("COORDINATES (BOHR)") < 0) currLine = buffRead.readLine();
currLine = buffRead.readLine();
// System.out.println(currLine);
for (int i = 0; i < niatms; i++) {
currLine = buffRead.readLine();
atomtypes[i] = currLine.substring(0,16);
atomweight[i] = (double)atomWeightInt[(int)Double.valueOf(atomtypes[i].substring(6).trim()).doubleValue()];
lx[i] = 0.52918*Double.valueOf(currLine.substring(17,33).trim()).doubleValue();
ly[i] = 0.52918*Double.valueOf(currLine.substring(34,53).trim()).doubleValue();
lz[i] = 0.52918*Double.valueOf(currLine.substring(54).trim()).doubleValue();
if (!fdt) {
ox[i] = lx[i];
oy[i] = ly[i];
oz[i] = lz[i];
}
}
}
if (fileType.equalsIgnoreCase("Jaguar")) {
// System.out.println("Looking for energy");
while (currLine.indexOf("Total energy") < 0) currLine = buffRead.readLine();
// System.out.println("found:"+currLine.substring(38,56)+"::");
lenergy = Double.valueOf(currLine.substring(38,56).trim()).doubleValue();
if (!fdt) oenergy = lenergy;
// System.out.println("Found energy"+oenergy);
} else if (fileType.equalsIgnoreCase("GAMESS")) {
while (currLine.indexOf("AFTER") < 0 || currLine.indexOf("ENERGY") < 0 ) currLine = buffRead.readLine();
lenergy = Double.valueOf(currLine.substring(currLine.indexOf("ENERGY")+11,currLine.indexOf("AFTER")).trim()).doubleValue
();
if (!fdt) oenergy = lenergy;
// System.out.println("Found energy: "+oenergy);
}
writeString = "Energy in hartrees; reaction coordinate in bohr (square root of amu)";
System.out.println(writeString);
buffWrite.write(writeString);
buffWrite.newLine();
writeString = "Energy Mass-weighted Non-mass-weighted";
System.out.println(writeString);
buffWrite.write(writeString);
buffWrite.newLine();
writeString = " origin last structure origin last structure";
System.out.println(writeString);
buffWrite.write(writeString);
buffWrite.newLine();
writeString = doubFormat(oenergy,12,6)+" "+doubFormat(0.0,12,6)+" "+doubFormat(0.0,12,6);
System.out.println(writeString);
buffWrite.write(writeString);
buffWrite.newLine();
setCentroid(ox, oy, oz, atomweight);
setCentroid(lx, ly, lz, atomweight);
if (fdt) {
quatosos = qtrfitPrep(lx, ly, lz, ox, oy, oz, niatms, atomweight);
// System.out.println("fdt 3: "+quatosos);
totquatlsos = quatosos;
energy = lenergy;
writeString = doubFormat(energy,12,6)+" "+doubFormat(quatosos/0.52918,12,6)+" "+doubFormat(totquatlsos/0.52918,12,
6);
System.out.println(writeString);
buffWrite.write(writeString);
buffWrite.newLine();
buffWrite.flush();
}
while (!lasttime) {
if (fileType.equalsIgnoreCase("Jaguar")) {
while (currLine.indexOf("new geometry:") < 0 && currLine.indexOf("final geometry:") < 0) currLine = buffRead.readLine(
);
if (currLine.indexOf("final geometry:") >= 0) lasttime = true;
// System.out.println("Found geometry "+lasttime);
currLine = buffRead.readLine();
currLine = buffRead.readLine();
for (int i = 0; i < niatms; i++) {
currLine = buffRead.readLine();
x[i] = Double.valueOf(currLine.substring(9,25).trim()).doubleValue();
y[i] = Double.valueOf(currLine.substring(25,43).trim()).doubleValue();
z[i] = Double.valueOf(currLine.substring(43).trim()).doubleValue();
}
while (currLine.indexOf("Total energy") < 0) currLine = buffRead.readLine();
energy = Double.valueOf(currLine.substring(38,56).trim()).doubleValue();
} else if (fileType.equalsIgnoreCase("GAMESS")) {
// System.out.println("Found geometry "+lasttime);
while (currLine.indexOf("COORDINATES") < 0 && currLine.indexOf("TERMINATED") < 0 || currLine.indexOf("NUCLEAR") >= 0)
currLine = buffRead.readLine();
if (currLine.indexOf("TERMINATED") > -1) {
lasttime = true;
} else {
currLine = buffRead.readLine();
currLine = buffRead.readLine();
// System.out.println("currLine: "+currLine);
for (int i = 0; i < niatms; i++) {
currLine = buffRead.readLine();
// System.out.println("Loop currLine: "+currLine);
x[i] = Double.valueOf(currLine.substring(17,32).trim()).doubleValue();
y[i] = Double.valueOf(currLine.substring(33,47).trim()).doubleValue();
z[i] = Double.valueOf(currLine.substring(48).trim()).doubleValue();
}
while (currLine.indexOf("AFTER") < 0 || currLine.indexOf("ENERGY") < 0 ) currLine = buffRead.readLine();
energy = Double.valueOf(currLine.substring(currLine.indexOf("ENERGY")+11,currLine.indexOf("AFTER")).trim()).doubleVa
lue();
// System.out.println("Found energy new "+energy);
}
}
// Now work out RMS distance for o and l
if (energy < lenergy) {
lenergy = energy;
osos = 0.0;
lsos = 0.0;
// should do rotation here
setCentroid(x, y, z, atomweight);
quatosos = qtrfitPrep(x, y, z, ox, oy, oz, niatms, atomweight);
quatlsos = qtrfitPrep(x, y, z, lx, ly, lz, niatms, atomweight);
totquatlsos += quatlsos;
for (int i = 0; i < niatms; i++) {
osos += (x[i]-ox[i])*(x[i]-ox[i]);
osos += (y[i]-oy[i])*(y[i]-oy[i]);
osos += (z[i]-oz[i])*(z[i]-oz[i]);
lsos += (x[i]-lx[i])*(x[i]-lx[i]);
lsos += (y[i]-ly[i])*(y[i]-ly[i]);
lsos += (z[i]-lz[i])*(z[i]-lz[i]);
}
osos = Math.sqrt(osos);
totlsos += Math.sqrt(lsos);
writeString = doubFormat(energy,12,6)+" "+doubFormat(quatosos/0.52918,12,6)+" "+doubFormat(totquatlsos/0.52918,12,
6)+" "+doubFormat(osos/0.52918,12,6)+" "+doubFormat(totlsos/0.52918,12,6);
System.out.println(writeString);
buffWrite.write(writeString);
buffWrite.newLine();
buffWrite.flush();
for (int i = 0; i < niatms; i++) {
lx[i] = x[i];
ly[i] = y[i];
lz[i] = z[i];
}
}
}
buffRead.close();
buffWrite.close();
} catch(Exception e) {;}
}
/**
* Create formatted strings for output
*/
static String doubFormat(double number, int length, int decPlace) {
String temp = " "+Double.toString(number)+"0000000000";
int j = temp.indexOf('.');
temp = temp.substring(0,j+decPlace+1);
int i = temp.length();
return temp.substring(i-length);
}
/**
* Convert atomtype to atomic weight
*/
static double atomtype2atomweight(String atomname) {
int atomnum;
String parsename, process;
process = atomname.trim();
if (process.length() > 1) {
if (Character.isLetter(process.charAt(1))) parsename = process.substring(0,2);
else parsename = process.substring(0,1);
} else parsename = process;
// System.out.println("Atomname: "+atomname+": :"+process+": :"+parsename);
atomnum = sym2num(parsename);
if (atomnum > atomWeightInt.length) System.out.println("Atomic weight too high");
return atomWeightInt[atomnum];
}
/**
* Convert atomic symbol to atomic number
*/
static int sym2num(String symbol) {
int atomnum = -1;
for (int i=0;i 0.0) {
dma = d[j] - d[i];
if (Math.abs(dma)+Math.abs(b) < Math.abs(dma)) {
t = b / dma;
} else {
q = 0.5 * dma / b;
t = fortranSign(1.0 / (Math.abs(q) + Math.sqrt(1.0 + q * q)), q);
}
c = 1.0 / Math.sqrt(t * t + 1.0);
s = t * c;
a[i][j] = 0.0;
for (k=0; k < i-1; k++) {
atemp = c * a[k][i] - s * a[k][j];
a[k][j] = s * a[k][i] + c * a[k][j];
a[k][i] = atemp;
}
for (k=i+1; k < j; k++) {
atemp = c * a[i][k] - s * a[k][j];
a[k][j] = s * a[i][k] + c * a[k][j];
a[i][k] = atemp;
}
for (k=j; k < n; k++) {
atemp = c * a[i][k] - s * a[j][k];
a[j][k] = s * a[i][k] + c * a[j][k];
a[i][k] = atemp;
}
for (k=0; k < n; k++) {
vtemp = c * v[k][i] - s * v[k][j];
v[k][j] = s * v[k][i] + c * v[k][j];
v[k][i] = vtemp;
}
dtemp = c * c * d[i] + s * s * d[j] - 2.0 * c * s * b;
d[j] = s * s * d[i] + c * c * d[j] + 2.0 * c * s * b;
d[i] = dtemp;
}
}
}
}
nrot = l;
for (j = 0; j < n-1; j++) {
k = j;
dtemp = d[k];
for (i = j; i < n; i++) {
if (d[i] < dtemp) {
k = i;
dtemp = d[k];
}
}
if (k > j) {
d[k] = d[j];
d[j] = dtemp;
for (i = 0; i < n; i++) {
dtemp = v[i][k];
v[i][k] = v[i][j];
v[i][j] = dtemp;
}
}
}
// System.out.println("Jacobi passes: "+nrot);
}
/**
* Q2MAT: Generate a left rotation matrix
* from a normalized quaternion;
* Input: q - normalized quaternion;
* Output: u - the rotation matrix;
*/
public static void q2mat (double[] q, double[][] u) {
u[0][0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3];
u[1][0] = 2.0 *(q[1] * q[2] - q[0] * q[3]);
u[2][0] = 2.0 *(q[1] * q[3] + q[0] * q[2]);
u[0][1] = 2.0 *(q[2] * q[1] + q[0] * q[3]);
u[1][1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3];
u[2][1] = 2.0 *(q[2] * q[3] - q[0] * q[1]);
u[0][2] = 2.0 *(q[3] * q[1] - q[0] * q[2]);
u[1][2] = 2.0 *(q[3] * q[2] + q[0] * q[1]);
u[2][2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3];
}
/**
* Java equivalent of Fortran SIGN function
*/
public static double fortranSign(double a, double b) {
if (b >= 0) return a;
else return -a;
}
}
© 2003 J M Goodman, Cambridge
|