1 /* $RCSfile$ 2 * $Author: hansonr $ 3 * $Date: 2006-09-12 00:46:22 -0500 (Tue, 12 Sep 2006) $ 4 * $Revision: 5501 $ 5 * 6 * Copyright (C) 2004-2005 The Jmol Development Team 7 * 8 * Contact: jmol-developers@lists.sf.net 9 * 10 * This library is free software; you can redistribute it and/or 11 * modify it under the terms of the GNU Lesser General Public 12 * License as published by the Free Software Foundation; either 13 * version 2.1 of the License, or (at your option) any later version. 14 * 15 * This library is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 * Lesser General Public License for more details. 19 * 20 * You should have received a copy of the GNU Lesser General Public 21 * License along with this library; if not, write to the Free Software 22 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 23 */ 24 25 package org.jmol.adapter.readers.quantum; 26 27 import org.jmol.adapter.smarter.Atom; 28 import org.jmol.adapter.smarter.Bond; 29 30 import javajs.util.AU; 31 import javajs.util.Lst; 32 import javajs.util.PT; 33 import javajs.util.V3; 34 35 import java.util.Hashtable; 36 import java.util.Map; 37 38 import org.jmol.api.JmolAdapter; 39 40 import org.jmol.util.Escape; 41 import org.jmol.util.Logger; 42 43 /** 44 * Reader for Gaussian fchk files 45 * for vibrational modes, add Freq=(SaveNormalModes,Raman,VibRot) 46 * 47 * 48 * also allows appended freq data 49 * 50 * 51 * @author hansonr Bob Hanson hansonr@stolaf.edu 52 * 53 **/ 54 public class GaussianFchkReader extends GaussianReader { 55 56 private Map<String, Object> fileData; 57 private int atomCount; 58 59 @Override initializeReader()60 protected void initializeReader() throws Exception { 61 super.initializeReader(); 62 energyUnits = ""; 63 fileData = new Hashtable<String, Object>(); 64 fileData.put("title", rd().trim()); 65 calculationType = PT.rep(rd(), " ", " "); 66 asc.newAtomSet(); 67 asc.setCurrentModelInfo("fileData", fileData); 68 readAllData(); 69 readAtoms(); 70 readBonds(); 71 readDipoleMoment(); 72 readPartialCharges(); 73 readBasis(); 74 readMolecularObitals(); 75 checkForFreq(); 76 continuing = false; 77 } 78 checkForFreq()79 private void checkForFreq() throws Exception { 80 Integer n = (Integer) fileData.get("Vib-NDim"); 81 if (n == null) { 82 // NumAtom 9 83 // NumFreq 21 84 // 1 2 3 85 // A2 B1 A2 86 // Frequencies -- 613.8891 622.9996 722.2497 87 // Red. masses -- 3.1195 3.8445 1.3156 88 // Frc consts -- 0.6927 0.8791 0.4043 89 readFrequencies("NumFreq", false); // if freq file appended 90 return; 91 } 92 try { 93 int nModes = n.intValue(); 94 float[] vibE2 = (float[]) fileData.get("Vib-E2"); 95 float[] modes = (float[]) fileData.get("Vib-Modes"); 96 float[] frequencies = fillFloat(vibE2, 0, nModes); 97 float[] red_masses = fillFloat(vibE2, nModes, nModes); 98 float[] frc_consts = fillFloat(vibE2, nModes * 2, nModes); 99 float[] intensities = fillFloat(vibE2, nModes * 3, nModes); 100 int ac = asc.getLastAtomSetAtomCount(); 101 boolean[] ignore = new boolean[nModes]; 102 int fpt = 0; 103 for (int i = 0; i < nModes; ++i) { 104 ignore[i] = !doGetVibration(++vibrationNumber); 105 if (ignore[i]) 106 continue; 107 int iAtom0 = asc.ac; 108 asc.cloneAtomSetWithBonds(true); 109 // set the properties 110 String name = asc.setAtomSetFrequency(vibrationNumber, "Calculation " + calculationNumber, null, "" + frequencies[i], null); 111 appendLoadNote("model " + asc.atomSetCount + ": " + name); 112 namedSets.set(asc.iSet); 113 asc.setAtomSetModelProperty("ReducedMass", 114 red_masses[i]+" AMU"); 115 asc.setAtomSetModelProperty("ForceConstant", 116 frc_consts[i]+" mDyne/A"); 117 asc.setAtomSetModelProperty("IRIntensity", 118 intensities[i]+" KM/Mole"); 119 for (int iAtom = 0; iAtom < ac; iAtom++) { 120 asc.addVibrationVectorWithSymmetry(iAtom0 121 + iAtom, modes[fpt++], modes[fpt++], modes[fpt++], false); 122 } 123 } 124 } catch (Exception e) { 125 Logger.error("Could not read Vib-E2 section: " + e.getMessage()); 126 } 127 128 } 129 fillFloat(float[] f0, int i, int n)130 private float[] fillFloat(float[] f0, int i, int n) { 131 float[] f = new float[n]; 132 for (int i1 = 0, ilast = i + n; i < ilast; i++, i1++) 133 f[i1] = f0[i]; 134 return f; 135 } 136 readAllData()137 private void readAllData() throws Exception { 138 while ((line == null ? rd() : line) != null) { 139 if (line.length() < 40) { 140 if (line.indexOf("NumAtom") == 0) { 141 // freq file appended 142 // NumAtom 9 143 // NumFreq 21 144 return; 145 } 146 continue; 147 } 148 String name = PT.rep(line.substring(0, 40).trim(), " ", ""); 149 char type = line.charAt(43); 150 boolean isArray = (line.indexOf("N=") >= 0); 151 String v = line.substring(50).trim(); 152 Logger.info(name + " = " + v + " " + isArray); 153 Object o = null; 154 if (isArray) { 155 switch (type) { 156 case 'I': 157 case 'R': 158 o = fillFloatArray(null, 0, new float[parseIntStr(v)]); 159 line = null; 160 break; 161 default: // C H L 162 v = rd().trim(); 163 while (rd() != null && line.indexOf(" N= ") < 0) 164 v += " " + line.trim(); 165 o = v; 166 break; 167 } 168 } else { 169 switch (type) { 170 case 'I': 171 o = Integer.valueOf(parseIntStr(v)); 172 break; 173 case 'R': 174 o = Double.valueOf(Double.parseDouble(v)); 175 break; 176 case 'C': 177 case 'L': 178 o = v; 179 break; 180 } 181 line = null; 182 } 183 if (o != null) 184 fileData.put(name, o); 185 } 186 } 187 188 @Override readAtoms()189 protected void readAtoms() throws Exception { 190 float[] atomNumbers = (float[]) fileData.get("Atomicnumbers"); 191 float[] data = (float[]) fileData.get("Currentcartesiancoordinates"); 192 String e = "" + fileData.get("TotalEnergy"); 193 asc.setAtomSetEnergy(e, parseFloatStr(e)); 194 atomCount = atomNumbers.length; 195 float f = ANGSTROMS_PER_BOHR; 196 for(int i = 0, pt = 0; i < atomCount; i++) { 197 Atom atom = asc.addNewAtom(); 198 atom.elementNumber = (short) atomNumbers[i]; 199 if (atom.elementNumber < 0) 200 atom.elementNumber = 0; // dummy atoms have -1 -> 0 201 setAtomCoordXYZ(atom, data[pt++] * f, data[pt++] * f, data[pt++] * f); 202 } 203 } 204 205 /* 206 MxBond I 3 207 NBond I N= 11 208 3 3 2 3 3 3 209 1 1 1 1 1 210 IBond I N= 33 211 2 3 7 1 4 8 212 1 6 0 2 5 9 213 4 6 10 3 5 11 214 1 0 0 2 0 0 215 4 0 0 5 0 0 216 6 0 0 217 RBond R N= 33 218 1.50000000E+00 1.50000000E+00 1.00000000E+00 1.50000000E+00 1.50000000E+00 219 1.00000000E+00 1.50000000E+00 1.50000000E+00 0.00000000E+00 1.50000000E+00 220 1.50000000E+00 1.00000000E+00 1.50000000E+00 1.50000000E+00 1.00000000E+00 221 1.50000000E+00 1.50000000E+00 1.00000000E+00 1.00000000E+00 0.00000000E+00 222 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 223 0.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 224 1.00000000E+00 0.00000000E+00 0.00000000E+00 225 */ 226 readBonds()227 protected void readBonds() { 228 try { 229 float[] nBond = (float[]) fileData.get("NBond"); 230 float[] iBond = (float[]) fileData.get("IBond"); 231 if (nBond.length == 0) 232 return; 233 float[] rBond = (float[]) fileData.get("RBond"); 234 // MxBond record is not critical here 235 int mxBond = rBond.length / nBond.length; 236 for (int ia = 0, pt = 0; ia < atomCount; ia++) 237 for (int j = 0; j < mxBond; j++, pt++) { 238 int ib = (int) iBond[pt] - 1; 239 if (ib <= ia) 240 continue; 241 float order = rBond[pt]; 242 int iorder = (order == 1.5f ? JmolAdapter.ORDER_AROMATIC 243 : (int) order); 244 asc.addBond(new Bond(ia, ib, iorder)); 245 } 246 addJmolScript("connect 1.1 {_H} {*} "); 247 } catch (Exception e) { 248 Logger.info("GaussianFchkReader -- bonding ignored"); 249 } 250 } 251 252 @Override readDipoleMoment()253 protected void readDipoleMoment() throws Exception { 254 float[] data = (float[]) fileData.get("DipoleMoment"); 255 if (data == null) 256 return; 257 V3 dipole = V3.new3(data[0], data[1], data[2]); 258 Logger.info("Molecular dipole for model " + asc.atomSetCount 259 + " = " + dipole); 260 asc.setCurrentModelInfo("dipole", dipole); 261 } 262 263 @Override readPartialCharges()264 protected void readPartialCharges() throws Exception { 265 float[] data = (float[]) fileData.get("Mulliken Charges"); 266 if (data == null) 267 return; 268 Atom[] atoms = asc.atoms; 269 for (int i = 0; i < atomCount; ++i) { 270 float c = data[i]; 271 atoms[i].partialCharge = c; 272 if (Math.abs(c) > 0.8f) 273 atoms[i].formalCharge = Math.round(c); 274 } 275 Logger.info("Mulliken charges found for Model " + asc.atomSetCount); 276 } 277 278 279 /* SAMPLE BASIS OUTPUT */ 280 /* 281 * see also http://www.gaussian.com/g_ur/k_gen.htm -- thank you, Rick Spinney 282 283 Standard basis: VSTO-3G (5D, 7F) 284 AO basis set: 285 Atom O1 Shell 1 SP 3 bf 1 - 4 0.000000000000 0.000000000000 0.216790088607 286 0.5033151319D+01 -0.9996722919D-01 0.1559162750D+00 287 0.1169596125D+01 0.3995128261D+00 0.6076837186D+00 288 0.3803889600D+00 0.7001154689D+00 0.3919573931D+00 289 Atom H2 Shell 2 S 3 bf 5 - 5 0.000000000000 1.424913022638 -0.867160354429 290 0.3425250914D+01 0.1543289673D+00 291 0.6239137298D+00 0.5353281423D+00 292 0.1688554040D+00 0.4446345422D+00 293 Atom H3 Shell 3 S 3 bf 6 - 6 0.000000000000 -1.424913022638 -0.867160354429 294 0.3425250914D+01 0.1543289673D+00 295 0.6239137298D+00 0.5353281423D+00 296 0.1688554040D+00 0.4446345422D+00 297 There are 3 symmetry adapted basis functions of A1 symmetry. 298 There are 0 symmetry adapted basis functions of A2 symmetry. 299 There are 1 symmetry adapted basis functions of B1 symmetry. 300 There are 2 symmetry adapted basis functions of B2 symmetry. 301 302 or 303 304 AO basis set in the form of general basis input (Overlap normalization): 305 1 0 306 S 3 1.00 0.000000000000 307 0.1508000000D 01 -0.1754110411D 00 308 0.5129000000D 00 -0.4465363900D 00 309 0.1362000000D 00 0.1295841966D 01 310 S 3 1.00 0.000000000000 311 0.2565000000D 01 -0.1043105923D 01 312 0.1508000000D 01 0.1331478902D 01 313 0.5129000000D 00 0.5613064585D 00 314 S 1 1.00 0.000000000000 315 0.4170000000D-01 0.1000000000D 01 316 P 3 1.00 0.000000000000 317 0.4859000000D 01 -0.9457549473D-01 318 0.1219000000D 01 0.7434797586D 00 319 0.4413000000D 00 0.3668143796D 00 320 P 2 1.00 0.000000000000 321 0.5725000000D 00 -0.8808640317D-01 322 0.8300000000D-01 0.1028397037D 01 323 P 1 1.00 0.000000000000 324 0.2500000000D-01 0.1000000000D 01 325 D 3 1.00 0.000000000000 326 0.4195000000D 01 0.4857290090D-01 327 0.1377000000D 01 0.5105223094D 00 328 0.4828000000D 00 0.5730028106D 00 329 D 1 1.00 0.000000000000 330 0.1501000000D 00 0.1000000000D 01 331 **** 332 2 0 333 ... 334 */ 335 336 //S,X,Y,Z,XX,YY,ZZ,XY,XZ,YZ,XXX,YYY,ZZZ,XYY,XXY,XXZ,XZZ,YZZ,YYZ,XYZ 337 private static String[] AO_TYPES = {"F7", "D5", "L", "S", "P", "D", "F", "G", "H"}; 338 @Override readBasis()339 protected void readBasis() throws Exception { 340 float[] types = (float[]) fileData.get("Shelltypes"); 341 gaussianCount = 0; 342 shellCount = 0; 343 if (types == null) 344 return; 345 shellCount = types.length; 346 shells = new Lst<int[]>(); 347 float[] pps = (float[]) fileData.get("Numberofprimitivespershell"); 348 float[] atomMap = (float[]) fileData.get("Shelltoatommap"); 349 float[] exps = (float[]) fileData.get("Primitiveexponents"); 350 float[] coefs = (float[]) fileData.get("Contractioncoefficients"); 351 float[] spcoefs = (float[]) fileData.get("P(S=P)Contractioncoefficients"); 352 gaussians = AU.newFloat2(exps.length); 353 for (int i = 0; i < shellCount; i++) { 354 String oType = AO_TYPES[(int) types[i] + 3]; 355 int nGaussians = (int) pps[i]; 356 int iatom = (int) atomMap[i]; 357 int[] slater = new int[4]; 358 slater[0] = iatom; // 1-based 359 if (oType.equals("F7") || oType.equals("D5")) 360 slater[1] = BasisFunctionReader.getQuantumShellTagIDSpherical(oType.substring(0, 1)); 361 else 362 slater[1] = BasisFunctionReader.getQuantumShellTagID(oType); 363 slater[2] = gaussianCount + 1; 364 slater[3] = nGaussians; 365 if (debugging) 366 Logger.debug("Slater " + shells.size() + " " + Escape.eAI(slater)); 367 shells.addLast(slater); 368 for (int j = 0; j < nGaussians; j++) { 369 float[] g = gaussians[gaussianCount] = new float[3]; 370 g[0] = exps[gaussianCount]; 371 g[1] = coefs[gaussianCount]; 372 if (spcoefs != null) 373 g[2] = spcoefs[gaussianCount]; 374 gaussianCount++; 375 } 376 } 377 Logger.info(shellCount + " slater shells read"); 378 Logger.info(gaussianCount + " gaussian primitives read"); 379 } 380 readMolecularObitals()381 protected void readMolecularObitals() throws Exception { 382 if (shells == null) 383 return; 384 int nElec = ((Integer) fileData.get("Numberofelectrons")).intValue(); 385 int nAlpha = ((Integer) fileData.get("Numberofalphaelectrons")).intValue(); 386 int nBeta = ((Integer) fileData.get("Numberofbetaelectrons")).intValue(); 387 //int mult = ((Integer) fileData.get("Multiplicity")).intValue(); 388 float[] aenergies = (float[]) fileData.get("AlphaOrbitalEnergies"); 389 float[] benergies = (float[]) fileData.get("BetaOrbitalEnergies"); 390 float[] acoefs = (float[]) fileData.get("AlphaMOcoefficients"); 391 float[] bcoefs = (float[]) fileData.get("BetaMOcoefficients"); 392 if (acoefs == null) 393 return; 394 int occ = (bcoefs == null ? 2 : 1); 395 int n = (bcoefs == null ? nElec : nAlpha); 396 getOrbitals(aenergies, acoefs, occ, n); 397 if (bcoefs != null) 398 getOrbitals(benergies, bcoefs, occ, nBeta); 399 setMOData(false); 400 } 401 getOrbitals(float[] e, float[] c, int occ, int nElec)402 private void getOrbitals(float[] e, float[] c, int occ, int nElec) { 403 int nOrb = e.length; 404 int nCoef = c.length; 405 nCoef /= nOrb; 406 alphaBeta = (occ == 2 ? "" : alphaBeta.equals("alpha") ? "beta" : "alpha"); 407 int pt = 0; 408 int n = 0; 409 for (int i = 0; i < nOrb; i++) { 410 float[] coefs = new float[nCoef]; 411 for (int j = 0; j < nCoef; j++) 412 coefs[j] = c[pt++]; 413 Map<String, Object> mo = new Hashtable<String, Object>(); 414 mo.put("coefficients", coefs); 415 mo.put("occupancy", Float.valueOf(occ)); 416 n += occ; 417 if (n >= nElec) 418 occ = 0; 419 mo.put("energy", Float.valueOf(e[i])); 420 mo.put("type", alphaBeta); 421 setMO(mo); 422 } 423 } 424 } 425