1
/* $Revision: 8201 $ $Author: egonw $ $Date: 2007-04-16 10:40:19 +0200 (Mon, 16 Apr 2007) $
3
* Copyright (C) 2002-2003 Bradley A. Smith <yeldar@home.com>
4
* Copyright (C) 2003-2007 Egon Willighagen <egonw@users.sf.net>
5
* Copyright (C) 2003-2007 Christoph Steinbeck <steinbeck@users.sf.net>
7
* Contact: cdk-devel@lists.sourceforge.net
9
* This library is free software; you can redistribute it and/or
10
* modify it under the terms of the GNU Lesser General Public
11
* License as published by the Free Software Foundation; either
12
* version 2.1 of the License, or (at your option) any later version.
14
* This library is distributed in the hope that it will be useful,
15
* but WITHOUT ANY WARRANTY; without even the implied warranty of
16
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17
* Lesser General Public License for more details.
19
* You should have received a copy of the GNU Lesser General Public
20
* License along with this library; if not, write to the Free Software
21
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
23
package org.openscience.cdk.io;
25
import java.io.BufferedReader;
26
import java.io.IOException;
27
import java.io.InputStream;
28
import java.io.InputStreamReader;
29
import java.io.Reader;
30
import java.io.StreamTokenizer;
31
import java.io.StringReader;
32
import java.util.List;
33
import java.util.StringTokenizer;
35
import javax.vecmath.Point3d;
37
import org.openscience.cdk.CDKConstants;
38
import org.openscience.cdk.config.IsotopeFactory;
39
import org.openscience.cdk.exception.CDKException;
40
import org.openscience.cdk.interfaces.IAtom;
41
import org.openscience.cdk.interfaces.IAtomContainer;
42
import org.openscience.cdk.interfaces.IChemFile;
43
import org.openscience.cdk.interfaces.IChemModel;
44
import org.openscience.cdk.interfaces.IChemObject;
45
import org.openscience.cdk.interfaces.IChemSequence;
46
import org.openscience.cdk.interfaces.IMolecule;
47
import org.openscience.cdk.interfaces.IMoleculeSet;
48
import org.openscience.cdk.io.formats.Gaussian98Format;
49
import org.openscience.cdk.io.formats.IResourceFormat;
50
import org.openscience.cdk.io.setting.BooleanIOSetting;
51
import org.openscience.cdk.io.setting.IOSetting;
52
import org.openscience.cdk.tools.LoggingTool;
53
import org.openscience.cdk.tools.manipulator.ChemModelManipulator;
56
* A reader for Gaussian98 output. Gaussian 98 is a quantum chemistry program
57
* by Gaussian, Inc. (<a href="http://www.gaussian.com/">http://www.gaussian.com/</a>).
59
* <p>Molecular coordinates, energies, and normal coordinates of vibrations are
60
* read. Each set of coordinates is added to the ChemFile in the order they are
61
* found. Energies and vibrations are associated with the previously read set
64
* <p>This reader was developed from a small set of example output files, and
65
* therefore, is not guaranteed to properly read all Gaussian98 output. If you
66
* have problems, please contact the author of this code, not the developers of
69
* @author Bradley A. Smith <yeldar@home.com>
70
* @author Egon Willighagen
71
* @author Christoph Steinbeck
74
public class Gaussian98Reader extends DefaultChemObjectReader {
76
private BufferedReader input;
77
private LoggingTool logger;
78
private int atomCount = 0;
79
private String lastRoute = "";
82
* Customizable setting
84
private BooleanIOSetting readOptimizedStructureOnly;
88
* Constructor for the Gaussian98Reader object
90
public Gaussian98Reader() {
91
this(new StringReader(""));
94
public Gaussian98Reader(InputStream input) {
95
this(new InputStreamReader(input));
98
public IResourceFormat getFormat() {
99
return Gaussian98Format.getInstance();
103
* Sets the reader attribute of the Gaussian98Reader object
105
* @param input The new reader value
106
* @throws CDKException Description of the Exception
108
public void setReader(Reader input) throws CDKException {
109
if (input instanceof BufferedReader) {
110
this.input = (BufferedReader) input;
112
this.input = new BufferedReader(input);
116
public void setReader(InputStream input) throws CDKException {
117
setReader(new InputStreamReader(input));
122
* Create an Gaussian98 output reader.
124
* @param input source of Gaussian98 data
126
public Gaussian98Reader(Reader input) {
127
logger = new LoggingTool(this);
128
if (input instanceof BufferedReader) {
129
this.input = (BufferedReader) input;
131
this.input = new BufferedReader(input);
136
public boolean accepts(Class classObject) {
137
Class[] interfaces = classObject.getInterfaces();
138
for (int i = 0; i < interfaces.length; i++) {
139
if (IChemFile.class.equals(interfaces[i])) return true;
144
public IChemObject read(IChemObject object) throws CDKException {
147
if (object instanceof IChemFile) {
148
IChemFile file = (IChemFile) object;
150
file = readChemFile(file);
151
} catch (IOException exception) {
152
throw new CDKException(
153
"Error while reading file: " + exception.toString(),
159
throw new CDKException("Reading of a " + object.getClass().getName() +
160
" is not supported.");
164
public void close() throws IOException {
170
* Read the Gaussian98 output.
172
* @return a ChemFile with the coordinates, energies, and
174
* @throws IOException if an I/O error occurs
175
* @throws CDKException Description of the Exception
177
private IChemFile readChemFile(IChemFile chemFile) throws CDKException, IOException {
178
IChemSequence sequence = chemFile.getBuilder().newChemSequence();
179
IChemModel model = null;
180
String line = input.readLine();
181
String levelOfTheory;
183
int modelCounter = 0;
185
// Find first set of coordinates by skipping all before "Standard orientation"
186
while (input.ready() && (line != null)) {
187
if (line.indexOf("Standard orientation:") >= 0) {
189
// Found a set of coordinates
190
model = chemFile.getBuilder().newChemModel();
191
readCoordinates(model);
194
line = input.readLine();
198
// Read all other data
199
line = input.readLine().trim();
200
while (input.ready() && (line != null)) {
201
if (line.indexOf("#") == 0) {
202
// Found the route section
203
// Memorizing this for the description of the chemmodel
207
} else if (line.indexOf("Standard orientation:") >= 0) {
209
// Found a set of coordinates
210
// Add current frame to file and create a new one.
211
if (!readOptimizedStructureOnly.isSet()) {
212
sequence.addChemModel(model);
214
logger.info("Skipping frame, because I was told to do");
217
model = chemFile.getBuilder().newChemModel();
219
readCoordinates(model);
220
} else if (line.indexOf("SCF Done:") >= 0) {
223
model.setProperty(CDKConstants.REMARK, line.trim());
224
} else if (line.indexOf("Harmonic frequencies") >= 0) {
226
// Found a set of vibrations
227
// readFrequencies(frame);
228
} else if (line.indexOf("Total atomic charges") >= 0) {
229
readPartialCharges(model);
230
} else if (line.indexOf("Magnetic shielding") >= 0) {
233
readNMRData(model, line);
235
} else if (line.indexOf("GINC") >= 0) {
237
// Found calculation level of theory
238
levelOfTheory = parseLevelOfTheory(line);
239
logger.debug("Level of Theory for this model: " + levelOfTheory);
240
description = lastRoute + ", model no. " + modelCounter;
241
model.setProperty(CDKConstants.DESCRIPTION, description);
243
//logger.debug("Skipping line: " + line);
245
line = input.readLine();
248
// Add last frame to file
249
sequence.addChemModel(model);
252
chemFile.addChemSequence(sequence);
259
* Reads a set of coordinates into ChemFrame.
261
* @param model Description of the Parameter
262
* @throws IOException if an I/O error occurs
263
* @throws CDKException Description of the Exception
265
private void readCoordinates(IChemModel model) throws CDKException, IOException {
266
IMoleculeSet moleculeSet = model.getBuilder().newMoleculeSet();
267
IMolecule molecule = model.getBuilder().newMolecule();
268
String line = input.readLine();
269
line = input.readLine();
270
line = input.readLine();
271
line = input.readLine();
272
while (input.ready()) {
273
line = input.readLine();
274
if ((line == null) || (line.indexOf("-----") >= 0)) {
278
StringReader sr = new StringReader(line);
279
StreamTokenizer token = new StreamTokenizer(sr);
282
// ignore first token
283
if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
284
atomicNumber = (int) token.nval;
285
if (atomicNumber == 0) {
287
// Skip dummy atoms. Dummy atoms must be skipped
288
// if frequencies are to be read because Gaussian
289
// does not report dummy atoms in frequencies, and
290
// the number of atoms is used for reading frequencies.
294
throw new CDKException("Error while reading coordinates: expected integer.");
298
// ignore third token
302
if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
305
throw new IOException("Error reading x coordinate");
307
if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
310
throw new IOException("Error reading y coordinate");
312
if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
315
throw new IOException("Error reading z coordinate");
317
String symbol = "Du";
319
symbol = IsotopeFactory.getInstance(model.getBuilder()).getElementSymbol(atomicNumber);
320
} catch (Exception exception) {
321
throw new CDKException("Could not determine element symbol!", exception);
323
IAtom atom = model.getBuilder().newAtom(symbol);
324
atom.setPoint3d(new Point3d(x, y, z));
325
molecule.addAtom(atom);
328
* this is the place where we store the atomcount to
329
* be used as a counter in the nmr reading
331
atomCount = molecule.getAtomCount();
332
moleculeSet.addMolecule(molecule);
333
model.setMoleculeSet(moleculeSet);
338
* Reads partial atomic charges and add the to the given ChemModel.
340
* @param model Description of the Parameter
341
* @throws CDKException Description of the Exception
342
* @throws IOException Description of the Exception
344
private void readPartialCharges(IChemModel model) throws CDKException, IOException {
345
logger.info("Reading partial atomic charges");
346
IMoleculeSet moleculeSet = model.getMoleculeSet();
347
IMolecule molecule = moleculeSet.getMolecule(0);
348
String line = input.readLine();
349
// skip first line after "Total atomic charges"
350
while (input.ready()) {
351
line = input.readLine();
352
logger.debug("Read charge block line: " + line);
353
if ((line == null) || (line.indexOf("Sum of Mulliken charges") >= 0)) {
354
logger.debug("End of charge block found");
357
StringReader sr = new StringReader(line);
358
StreamTokenizer tokenizer = new StreamTokenizer(sr);
359
if (tokenizer.nextToken() == StreamTokenizer.TT_NUMBER) {
360
int atomCounter = (int) tokenizer.nval;
362
tokenizer.nextToken();
366
if (tokenizer.nextToken() == StreamTokenizer.TT_NUMBER) {
367
charge = tokenizer.nval;
368
logger.debug("Found charge for atom " + atomCounter +
371
throw new CDKException("Error while reading charge: expected double.");
373
IAtom atom = molecule.getAtom(atomCounter - 1);
374
atom.setCharge(charge);
380
* Reads a set of vibrations into ChemFrame.
382
*@param model Description of the Parameter
383
*@exception IOException if an I/O error occurs
385
// private void readFrequencies(IChemModel model) throws IOException
388
* FIXME: this is yet to be ported
390
* line = input.readLine();
391
* line = input.readLine();
392
* line = input.readLine();
393
* line = input.readLine();
394
* line = input.readLine();
395
* while ((line != null) && line.startsWith(" Frequencies --")) {
396
* Vector currentVibs = new Vector();
397
* StringReader vibValRead = new StringReader(line.substring(15));
398
* StreamTokenizer token = new StreamTokenizer(vibValRead);
399
* while (token.nextToken() != StreamTokenizer.TT_EOF) {
400
* Vibration vib = new Vibration(Double.toString(token.nval));
401
* currentVibs.addElement(vib);
403
* line = input.readLine();
404
* line = input.readLine();
405
* line = input.readLine();
406
* line = input.readLine();
407
* line = input.readLine();
408
* line = input.readLine();
409
* for (int i = 0; i < frame.getAtomCount(); ++i) {
410
* line = input.readLine();
411
* StringReader vectorRead = new StringReader(line);
412
* token = new StreamTokenizer(vectorRead);
414
* / ignore first token
416
* / ignore second token
417
* for (int j = 0; j < currentVibs.size(); ++j) {
418
* double[] v = new double[3];
419
* if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
422
* throw new IOException("Error reading frequency");
424
* if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
427
* throw new IOException("Error reading frequency");
429
* if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
432
* throw new IOException("Error reading frequency");
434
* ((Vibration) currentVibs.elementAt(j)).addAtomVector(v);
437
* for (int i = 0; i < currentVibs.size(); ++i) {
438
* frame.addVibration((Vibration) currentVibs.elementAt(i));
440
* line = input.readLine();
441
* line = input.readLine();
442
* line = input.readLine();
449
* Reads NMR nuclear shieldings.
451
private void readNMRData(IChemModel model, String labelLine) throws CDKException {
452
List containers = ChemModelManipulator.getAllAtomContainers(model);
453
if (containers.size() == 0) {
454
// nothing to store the results into
456
} // otherwise insert in the first AC
458
IAtomContainer ac = (IAtomContainer)containers.get(0);
459
// Determine label for properties
461
if (labelLine.indexOf("Diamagnetic") >= 0) {
462
label = "Diamagnetic Magnetic shielding (Isotropic)";
463
} else if (labelLine.indexOf("Paramagnetic") >= 0) {
464
label = "Paramagnetic Magnetic shielding (Isotropic)";
466
label = "Magnetic shielding (Isotropic)";
469
for (int i = 0; i < atomCount; ++i) {
471
String line = input.readLine().trim();
472
while (line.indexOf("Isotropic") < 0) {
476
line = input.readLine().trim();
478
StringTokenizer st1 = new StringTokenizer(line);
480
// Find Isotropic label
481
while (st1.hasMoreTokens()) {
482
if (st1.nextToken().equals("Isotropic")) {
487
// Find Isotropic value
488
while (st1.hasMoreTokens()) {
489
if (st1.nextToken().equals("=")) break;
491
double shielding = Double.valueOf(st1.nextToken()).doubleValue();
492
logger.info("Type of shielding: " + label);
493
ac.getAtom(atomIndex).setProperty(CDKConstants.ISOTROPIC_SHIELDING, new Double(shielding));
495
} catch (Exception exc) {
496
logger.debug("failed to read line from gaussian98 file where I expected one.");
503
* Select the theory and basis set from the first archive line.
505
* @param line Description of the Parameter
506
* @return Description of the Return Value
508
private String parseLevelOfTheory(String line) {
509
StringBuffer summary = new StringBuffer();
510
summary.append(line);
514
line = input.readLine().trim();
515
summary.append(line);
516
} while (!(line.indexOf("@") >= 0));
518
catch (Exception exc) {
519
logger.debug("syntax problem while parsing summary of g98 section: ");
522
logger.debug("parseLoT(): " + summary.toString());
523
StringTokenizer st1 = new StringTokenizer(summary.toString(), "\\");
525
// Must contain at least 6 tokens
526
if (st1.countTokens() < 6) {
530
// Skip first four tokens
531
for (int i = 0; i < 4; ++i) {
535
return st1.nextToken() + "/" + st1.nextToken();
539
private void initIOSettings() {
540
readOptimizedStructureOnly = new BooleanIOSetting("ReadOptimizedStructureOnly", IOSetting.LOW,
541
"Should I only read the optimized structure from a geometry optimization?",
545
private void customizeJob() {
546
fireIOSettingQuestion(readOptimizedStructureOnly);
551
* Gets the iOSettings attribute of the Gaussian98Reader object
553
*@return The iOSettings value
555
public IOSetting[] getIOSettings() {
556
IOSetting[] settings = new IOSetting[1];
557
settings[0] = readOptimizedStructureOnly;