~ubuntu-branches/ubuntu/trusty/cdk/trusty-proposed

« back to all changes in this revision

Viewing changes to src/org/openscience/cdk/io/Gaussian98Reader.java

  • Committer: Bazaar Package Importer
  • Author(s): Paul Cager
  • Date: 2008-04-09 21:17:53 UTC
  • Revision ID: james.westby@ubuntu.com-20080409211753-46lmjw5z8mx5pd8d
Tags: upstream-1.0.2
ImportĀ upstreamĀ versionĀ 1.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* $Revision: 8201 $ $Author: egonw $ $Date: 2007-04-16 10:40:19 +0200 (Mon, 16 Apr 2007) $
 
2
 *
 
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>
 
6
 *
 
7
 *  Contact: cdk-devel@lists.sourceforge.net
 
8
 *
 
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.
 
13
 *
 
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.
 
18
 *
 
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.
 
22
 */
 
23
package org.openscience.cdk.io;
 
24
 
 
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;
 
34
 
 
35
import javax.vecmath.Point3d;
 
36
 
 
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;
 
54
 
 
55
/**
 
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>).
 
58
 * <p/>
 
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
 
62
 * of coordinates.
 
63
 * <p/>
 
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
 
67
 * Gaussian98.
 
68
 *
 
69
 * @author Bradley A. Smith <yeldar@home.com>
 
70
 * @author Egon Willighagen
 
71
 * @author Christoph Steinbeck
 
72
 * @cdk.module io
 
73
 */
 
74
public class Gaussian98Reader extends DefaultChemObjectReader {
 
75
 
 
76
    private BufferedReader input;
 
77
    private LoggingTool logger;
 
78
    private int atomCount = 0;
 
79
    private String lastRoute = "";
 
80
 
 
81
    /**
 
82
     * Customizable setting
 
83
     */
 
84
    private BooleanIOSetting readOptimizedStructureOnly;
 
85
 
 
86
 
 
87
    /**
 
88
     * Constructor for the Gaussian98Reader object
 
89
     */
 
90
    public Gaussian98Reader() {
 
91
        this(new StringReader(""));
 
92
    }
 
93
 
 
94
    public Gaussian98Reader(InputStream input) {
 
95
        this(new InputStreamReader(input));
 
96
    }
 
97
 
 
98
    public IResourceFormat getFormat() {
 
99
        return Gaussian98Format.getInstance();
 
100
    }
 
101
 
 
102
    /**
 
103
     * Sets the reader attribute of the Gaussian98Reader object
 
104
     *
 
105
     * @param input The new reader value
 
106
     * @throws CDKException Description of the Exception
 
107
     */
 
108
    public void setReader(Reader input) throws CDKException {
 
109
        if (input instanceof BufferedReader) {
 
110
            this.input = (BufferedReader) input;
 
111
        } else {
 
112
            this.input = new BufferedReader(input);
 
113
        }
 
114
    }
 
115
 
 
116
    public void setReader(InputStream input) throws CDKException {
 
117
        setReader(new InputStreamReader(input));
 
118
    }
 
119
 
 
120
 
 
121
    /**
 
122
     * Create an Gaussian98 output reader.
 
123
     *
 
124
     * @param input source of Gaussian98 data
 
125
     */
 
126
    public Gaussian98Reader(Reader input) {
 
127
        logger = new LoggingTool(this);
 
128
        if (input instanceof BufferedReader) {
 
129
            this.input = (BufferedReader) input;
 
130
        } else {
 
131
            this.input = new BufferedReader(input);
 
132
        }
 
133
        initIOSettings();
 
134
    }
 
135
 
 
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;
 
140
        }
 
141
        return false;
 
142
    }
 
143
 
 
144
    public IChemObject read(IChemObject object) throws CDKException {
 
145
        customizeJob();
 
146
 
 
147
        if (object instanceof IChemFile) {
 
148
            IChemFile file = (IChemFile) object;
 
149
            try {
 
150
                file = readChemFile(file);
 
151
            } catch (IOException exception) {
 
152
                throw new CDKException(
 
153
                        "Error while reading file: " + exception.toString(),
 
154
                        exception
 
155
                );
 
156
            }
 
157
            return file;
 
158
        } else {
 
159
            throw new CDKException("Reading of a " + object.getClass().getName() +
 
160
                    " is not supported.");
 
161
        }
 
162
    }
 
163
 
 
164
    public void close() throws IOException {
 
165
        input.close();
 
166
    }
 
167
 
 
168
 
 
169
    /**
 
170
     * Read the Gaussian98 output.
 
171
     *
 
172
     * @return a ChemFile with the coordinates, energies, and
 
173
     *         vibrations.
 
174
     * @throws IOException  if an I/O error occurs
 
175
     * @throws CDKException Description of the Exception
 
176
     */
 
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;
 
182
        String description;
 
183
        int modelCounter = 0;
 
184
 
 
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) {
 
188
 
 
189
                // Found a set of coordinates
 
190
                model = chemFile.getBuilder().newChemModel();
 
191
                readCoordinates(model);
 
192
                break;
 
193
            }
 
194
            line = input.readLine();
 
195
        }
 
196
        if (model != null) {
 
197
 
 
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
 
204
                    lastRoute = line;
 
205
                    modelCounter = 0;
 
206
 
 
207
                } else if (line.indexOf("Standard orientation:") >= 0) {
 
208
 
 
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);
 
213
                    } else {
 
214
                        logger.info("Skipping frame, because I was told to do");
 
215
                    }
 
216
                    fireFrameRead();
 
217
                    model = chemFile.getBuilder().newChemModel();
 
218
                    modelCounter++;
 
219
                    readCoordinates(model);
 
220
                } else if (line.indexOf("SCF Done:") >= 0) {
 
221
 
 
222
                    // Found an energy
 
223
                    model.setProperty(CDKConstants.REMARK, line.trim());
 
224
                } else if (line.indexOf("Harmonic frequencies") >= 0) {
 
225
 
 
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) {
 
231
 
 
232
                    // Found NMR data
 
233
                    readNMRData(model, line);
 
234
 
 
235
                } else if (line.indexOf("GINC") >= 0) {
 
236
 
 
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);
 
242
                } else {
 
243
                    //logger.debug("Skipping line: " + line);
 
244
                }
 
245
                line = input.readLine();
 
246
            }
 
247
 
 
248
            // Add last frame to file
 
249
            sequence.addChemModel(model);
 
250
            fireFrameRead();
 
251
        }
 
252
        chemFile.addChemSequence(sequence);
 
253
 
 
254
        return chemFile;
 
255
    }
 
256
 
 
257
 
 
258
    /**
 
259
     * Reads a set of coordinates into ChemFrame.
 
260
     *
 
261
     * @param model Description of the Parameter
 
262
     * @throws IOException  if an I/O error occurs
 
263
     * @throws CDKException Description of the Exception
 
264
     */
 
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)) {
 
275
                break;
 
276
            }
 
277
            int atomicNumber;
 
278
            StringReader sr = new StringReader(line);
 
279
            StreamTokenizer token = new StreamTokenizer(sr);
 
280
            token.nextToken();
 
281
 
 
282
            // ignore first token
 
283
            if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
 
284
                atomicNumber = (int) token.nval;
 
285
                if (atomicNumber == 0) {
 
286
 
 
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.
 
291
                    continue;
 
292
                }
 
293
            } else {
 
294
                throw new CDKException("Error while reading coordinates: expected integer.");
 
295
            }
 
296
            token.nextToken();
 
297
 
 
298
            // ignore third token
 
299
            double x;
 
300
            double y;
 
301
            double z;
 
302
            if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
 
303
                x = token.nval;
 
304
            } else {
 
305
                throw new IOException("Error reading x coordinate");
 
306
            }
 
307
            if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
 
308
                y = token.nval;
 
309
            } else {
 
310
                throw new IOException("Error reading y coordinate");
 
311
            }
 
312
            if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
 
313
                z = token.nval;
 
314
            } else {
 
315
                throw new IOException("Error reading z coordinate");
 
316
            }
 
317
            String symbol = "Du";
 
318
            try {
 
319
                symbol = IsotopeFactory.getInstance(model.getBuilder()).getElementSymbol(atomicNumber);
 
320
            } catch (Exception exception) {
 
321
                throw new CDKException("Could not determine element symbol!", exception);
 
322
            }
 
323
            IAtom atom = model.getBuilder().newAtom(symbol);
 
324
            atom.setPoint3d(new Point3d(x, y, z));
 
325
            molecule.addAtom(atom);
 
326
        }
 
327
        /*
 
328
           *  this is the place where we store the atomcount to
 
329
           *  be used as a counter in the nmr reading
 
330
           */
 
331
        atomCount = molecule.getAtomCount();
 
332
        moleculeSet.addMolecule(molecule);
 
333
        model.setMoleculeSet(moleculeSet);
 
334
    }
 
335
 
 
336
 
 
337
    /**
 
338
     * Reads partial atomic charges and add the to the given ChemModel.
 
339
     *
 
340
     * @param model Description of the Parameter
 
341
     * @throws CDKException Description of the Exception
 
342
     * @throws IOException  Description of the Exception
 
343
     */
 
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");
 
355
                break;
 
356
            }
 
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;
 
361
 
 
362
                tokenizer.nextToken();
 
363
                // ignore the symbol
 
364
 
 
365
                double charge;
 
366
                if (tokenizer.nextToken() == StreamTokenizer.TT_NUMBER) {
 
367
                    charge = tokenizer.nval;
 
368
                    logger.debug("Found charge for atom " + atomCounter +
 
369
                            ": " + charge);
 
370
                } else {
 
371
                    throw new CDKException("Error while reading charge: expected double.");
 
372
                }
 
373
                IAtom atom = molecule.getAtom(atomCounter - 1);
 
374
                atom.setCharge(charge);
 
375
            }
 
376
        }
 
377
    }
 
378
 
 
379
    /**
 
380
     *  Reads a set of vibrations into ChemFrame.
 
381
     *
 
382
     *@param  model            Description of the Parameter
 
383
     *@exception IOException  if an I/O error occurs
 
384
     */
 
385
//      private void readFrequencies(IChemModel model) throws IOException
 
386
//      {
 
387
    /*
 
388
          *  FIXME: this is yet to be ported
 
389
          *  String line;
 
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);
 
402
          *  }
 
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);
 
413
          *  token.nextToken();
 
414
          *  / ignore first token
 
415
          *  token.nextToken();
 
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) {
 
420
          *  v[0] = token.nval;
 
421
          *  } else {
 
422
          *  throw new IOException("Error reading frequency");
 
423
          *  }
 
424
          *  if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
 
425
          *  v[1] = token.nval;
 
426
          *  } else {
 
427
          *  throw new IOException("Error reading frequency");
 
428
          *  }
 
429
          *  if (token.nextToken() == StreamTokenizer.TT_NUMBER) {
 
430
          *  v[2] = token.nval;
 
431
          *  } else {
 
432
          *  throw new IOException("Error reading frequency");
 
433
          *  }
 
434
          *  ((Vibration) currentVibs.elementAt(j)).addAtomVector(v);
 
435
          *  }
 
436
          *  }
 
437
          *  for (int i = 0; i < currentVibs.size(); ++i) {
 
438
          *  frame.addVibration((Vibration) currentVibs.elementAt(i));
 
439
          *  }
 
440
          *  line = input.readLine();
 
441
          *  line = input.readLine();
 
442
          *  line = input.readLine();
 
443
          *  }
 
444
          */
 
445
//      }
 
446
 
 
447
 
 
448
    /**
 
449
     * Reads NMR nuclear shieldings.
 
450
     */
 
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
 
455
                return;
 
456
        } // otherwise insert in the first AC
 
457
        
 
458
        IAtomContainer ac = (IAtomContainer)containers.get(0);
 
459
        // Determine label for properties
 
460
        String label;
 
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)";
 
465
        } else {
 
466
            label = "Magnetic shielding (Isotropic)";
 
467
        }
 
468
        int atomIndex = 0;
 
469
        for (int i = 0; i < atomCount; ++i) {
 
470
            try {
 
471
                String line = input.readLine().trim();
 
472
                while (line.indexOf("Isotropic") < 0) {
 
473
                    if (line == null) {
 
474
                        return;
 
475
                    }
 
476
                    line = input.readLine().trim();
 
477
                }
 
478
                StringTokenizer st1 = new StringTokenizer(line);
 
479
 
 
480
                // Find Isotropic label
 
481
                while (st1.hasMoreTokens()) {
 
482
                    if (st1.nextToken().equals("Isotropic")) {
 
483
                        break;
 
484
                    }
 
485
                }
 
486
 
 
487
                // Find Isotropic value
 
488
                while (st1.hasMoreTokens()) {
 
489
                    if (st1.nextToken().equals("=")) break;
 
490
                }
 
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));
 
494
                ++atomIndex;
 
495
            } catch (Exception exc) {
 
496
                logger.debug("failed to read line from gaussian98 file where I expected one.");
 
497
            }
 
498
        }
 
499
    }
 
500
 
 
501
 
 
502
    /**
 
503
     * Select the theory and basis set from the first archive line.
 
504
     *
 
505
     * @param line Description of the Parameter
 
506
     * @return Description of the Return Value
 
507
     */
 
508
    private String parseLevelOfTheory(String line) {
 
509
        StringBuffer summary = new StringBuffer();
 
510
        summary.append(line);
 
511
        try {
 
512
 
 
513
            do {
 
514
                line = input.readLine().trim();
 
515
                summary.append(line);
 
516
            } while (!(line.indexOf("@") >= 0));
 
517
        }
 
518
        catch (Exception exc) {
 
519
            logger.debug("syntax problem while parsing summary of g98 section: ");
 
520
            logger.debug(exc);
 
521
        }
 
522
        logger.debug("parseLoT(): " + summary.toString());
 
523
        StringTokenizer st1 = new StringTokenizer(summary.toString(), "\\");
 
524
 
 
525
        // Must contain at least 6 tokens
 
526
        if (st1.countTokens() < 6) {
 
527
            return null;
 
528
        }
 
529
 
 
530
        // Skip first four tokens
 
531
        for (int i = 0; i < 4; ++i) {
 
532
            st1.nextToken();
 
533
        }
 
534
 
 
535
        return st1.nextToken() + "/" + st1.nextToken();
 
536
    }
 
537
 
 
538
 
 
539
    private void initIOSettings() {
 
540
        readOptimizedStructureOnly = new BooleanIOSetting("ReadOptimizedStructureOnly", IOSetting.LOW,
 
541
                "Should I only read the optimized structure from a geometry optimization?",
 
542
                "false");
 
543
    }
 
544
 
 
545
    private void customizeJob() {
 
546
        fireIOSettingQuestion(readOptimizedStructureOnly);
 
547
    }
 
548
 
 
549
 
 
550
    /**
 
551
     *  Gets the iOSettings attribute of the Gaussian98Reader object
 
552
     *
 
553
     *@return The iOSettings value
 
554
     */
 
555
    public IOSetting[] getIOSettings() {
 
556
        IOSetting[] settings = new IOSetting[1];
 
557
        settings[0] = readOptimizedStructureOnly;
 
558
                return settings;
 
559
        }
 
560
 
 
561
}
 
562