1 package org.jmol.adapter.readers.xtal; 2 3 import java.util.Hashtable; 4 import java.util.Map; 5 6 7 import org.jmol.adapter.smarter.AtomSetCollectionReader; 8 import org.jmol.adapter.smarter.Atom; 9 import org.jmol.api.SymmetryInterface; 10 11 import javajs.util.PT; 12 import javajs.util.V3; 13 14 /** 15 * Problems identified (Bob Hanson) -- 16 * 17 * -- Coordinates for the asymmetric unit are conventional. 18 * Default right now is to read conventional cell, not primitive celll 19 * 20 * -- Frequency data number of atoms does not correspond to initial atom count. 21 * It looks like there is a missing report of symmetry-generated atoms. 22 * 23 * see https://projects.ivec.org/gulp/ 24 * 25 * @author Pieremanuele Canepa, Room 104, FM Group School of Physical Sciences, 26 * Ingram Building, University of Kent, Canterbury, Kent, CT2 7NH United 27 * Kingdom, pc229@kent.ac.uk 28 * 29 * 30 * 31 * @version 1.0 32 */ 33 34 public class GulpReader extends AtomSetCollectionReader { 35 36 private boolean isSlab; 37 private boolean isPolymer; 38 private boolean isPrimitive; 39 private String sep = "-------"; 40 private boolean coordinatesArePrimitive; 41 private Map<String, Float> atomCharges; 42 43 @Override initializeReader()44 protected void initializeReader() throws Exception { 45 isPrimitive = !checkFilterKey("CONV"); 46 coordinatesArePrimitive = true; 47 setFractionalCoordinates(readDimensionality()); 48 } 49 50 @Override finalizeSubclassReader()51 protected void finalizeSubclassReader() { 52 if (atomCharges == null) 53 return; 54 Atom[] atoms = asc.atoms; 55 Float f; 56 for (int i = asc.ac; --i >= 0;) 57 if ((f = atomCharges.get(atoms[i].atomName)) != null 58 || (f = atomCharges.get(atoms[i].getElementSymbol())) != null) 59 atoms[i].partialCharge = f.floatValue(); 60 } 61 62 private boolean bTest; 63 64 @Override checkLine()65 protected boolean checkLine() throws Exception { 66 if (line.contains("Space group ")) { 67 readSpaceGroup(); 68 return true; 69 } 70 71 if (isSlab ? line.contains("Surface cell parameters") 72 : isPolymer ? line.contains("Polymer cell parameter") 73 : (bTest = line.contains("Cartesian lattice vectors")) 74 || line.contains("Cell parameters (Angstroms/Degrees)") 75 || line.contains("Primitive cell parameters")) { 76 readCellParameters(bTest); 77 return true; 78 } 79 if (line.contains("Monopole - monopole (total)")) { 80 readEnergy(); 81 return true; 82 } 83 84 if (line.contains("Fractional coordinates of asymmetric unit :") 85 || (bTest = line.contains("Final asymmetric unit coordinates")) 86 || (bTest = line.contains("Final fractional coordinates ")) 87 || line 88 .contains("Mixed fractional/Cartesian coordinates") 89 || line.contains("Cartesian coordinates of cluster ") 90 || line.contains("Final cartesian coordinates of atoms :") 91 && isMolecular) { 92 if (doGetModel(++modelNumber, null)) 93 readAtomicPos(!bTest); 94 return true; 95 } 96 97 if (line.contains("Species output for all configurations")) { 98 readPartialCharges(); 99 return true; 100 } 101 102 // past this point, we must already have coordinates defined 103 if (!doProcessLines) 104 return true; 105 if (line.contains("Final cell parameters and derivatives")) { 106 // this line comes AFTER atom positions 107 readFinalCell(); 108 return true; 109 } 110 111 //if (line.contains(" Phonon Calculation : ")) { 112 // readFrequency(); 113 114 return true; 115 } 116 readDimensionality()117 private boolean readDimensionality() throws Exception { 118 discardLinesUntilContains("Dimensionality"); 119 String[] tokens = getTokens(); 120 switch (parseIntStr(tokens[2])) { 121 case 0: 122 isMolecular = true; 123 isPrimitive = false; 124 return false; 125 case 1: 126 isPolymer = true; 127 isPrimitive = false; 128 break; 129 case 2: 130 isSlab = true; 131 isPrimitive = false; 132 break; 133 } 134 return true; 135 } 136 readSpaceGroup()137 private void readSpaceGroup() throws Exception { 138 sgName = line.substring(line.indexOf(":") + 1).trim(); 139 } 140 141 private float a, b, c, alpha, beta, gamma; 142 private float[] primitiveData; 143 final private static String[] tags = { "a", "b", "c", "alpha", "beta", "gamma" }; 144 parameterIndex(String key)145 private static int parameterIndex(String key) { 146 for (int i = tags.length; --i >= 0;) 147 if (tags[i].equals(key)) 148 return i; 149 return -1; 150 } 151 setParameter(String key, float value)152 private void setParameter(String key, float value) { 153 switch (parameterIndex(key)) { 154 case 0: 155 a = value; 156 break; 157 case 1: 158 b = value; 159 break; 160 case 2: 161 c = value; 162 break; 163 case 3: 164 alpha = value; 165 break; 166 case 4: 167 beta = value; 168 break; 169 case 5: 170 gamma = value; 171 break; 172 } 173 } 174 newAtomSet(boolean doSetUnitCell)175 private void newAtomSet(boolean doSetUnitCell) { 176 asc.newAtomSet(); 177 if (doSetUnitCell) { 178 setModelParameters(coordinatesArePrimitive); 179 if (totEnergy != null) 180 setEnergy(); 181 } 182 } 183 setModelParameters(boolean isPrimitive)184 private void setModelParameters(boolean isPrimitive) { 185 if (sgName != null) 186 setSpaceGroupName(isPrimitive ? "P1" : sgName); 187 if (isPrimitive && primitiveData != null) { 188 addExplicitLatticeVector(0, primitiveData, 0); 189 addExplicitLatticeVector(1, primitiveData, 3); 190 addExplicitLatticeVector(2, primitiveData, 6); 191 } else if (a != 0) { 192 if (isSlab) { 193 c = -1; 194 beta = gamma = 90; 195 } else if (isPolymer) { 196 b = c = -1; 197 alpha = beta = gamma = 90; 198 } 199 setUnitCell(a, b, c, alpha, beta, gamma); 200 } 201 } 202 203 /* 204 205 Cartesian lattice vectors (Angstroms) : 206 207 10.944693 0.000000 0.000000 208 3.123705 4.493221 0.000000 209 3.123705 1.632784 4.186054 210 211 Cell parameters (Angstroms/Degrees): 212 213 a = 10.9447 alpha = 55.1928 214 b = 5.4723 beta = 55.1928 215 c = 5.4723 gamma = 55.1928 216 217 Primitive cell parameters : Full cell parameters : 218 219 a = 5.1295 alpha = 55.2915 a = 4.7602 alpha = 90.0000 220 b = 5.1295 beta = 55.2915 b = 4.7602 beta = 90.0000 221 c = 5.1295 gamma = 55.2915 c = 12.9933 gamma = 120.0000 222 223 */ 224 readCellParameters(boolean isLatticeVectors)225 private void readCellParameters(boolean isLatticeVectors) throws Exception { 226 if (isLatticeVectors) { 227 rd(); 228 primitiveData = fillFloatArray(null, 0, new float[9]); 229 a = 0; 230 return; 231 } 232 int i0 = (line.indexOf("Full cell") < 0 ? 0 : 4); 233 coordinatesArePrimitive = (i0 == 0); 234 //if (!coordinatesArePrimitive) 235 //isPrimitive = false; 236 rd(); 237 while (rd() != null && line.contains("=")) { 238 String[] tokens = PT.getTokens(line.replace('=', ' ')); 239 for (int i = i0; i < i0 + 4; i += 2) 240 if (tokens.length > i + 1) 241 setParameter(tokens[i], parseFloatStr(tokens[i + 1])); 242 } 243 } 244 245 /* 246 Final cell parameters and derivatives : 247 248 -------------------------------------------------------------------------------- 249 a 5.153230 Angstrom dE/de1(xx) 0.000090 eV/strain 250 b 5.153230 Angstrom dE/de2(yy) 0.000000 eV/strain 251 c 5.153230 Angstrom dE/de3(zz) 0.000078 eV/strain 252 alpha 55.766721 Degrees dE/de4(yz) 0.000000 eV/strain 253 beta 55.766721 Degrees dE/de5(xz) 0.000000 eV/strain 254 gamma 55.766721 Degrees dE/de6(xy) 0.000000 eV/strain 255 -------------------------------------------------------------------------------- 256 257 Non-primitive lattice parameters : 258 259 a = 4.820055 b = 4.820055 c = 13.011659 260 alpha= 90.000000 beta= 90.000000 gamma= 120.000000 261 262 */ 263 readFinalCell()264 private void readFinalCell() throws Exception { 265 // problem here is if we have primitive data, we have to change 266 // that data by only changing vector lengths 267 discardLinesUntilContains(sep); 268 String tokens[]; 269 while (rd() != null && (tokens = getTokens()).length >= 2) 270 setParameter(tokens[0], parseFloatStr(tokens[1])); 271 if (primitiveData != null) { 272 scalePrimitiveData(0, a); 273 scalePrimitiveData(3, b); 274 scalePrimitiveData(6, c); 275 if (!coordinatesArePrimitive) 276 // we have a conventional cell -- get a, b, and c for it now 277 while (rd() != null && line.indexOf("Final") < 0) 278 if (line.indexOf("Non-primitive lattice parameters") > 0) { 279 rd(); 280 for (int i = 0; i < 2; i++) { 281 tokens = PT.getTokens(rd().replace('=', ' ')); 282 setParameter(tokens[0], parseFloatStr(tokens[1])); 283 setParameter(tokens[2], parseFloatStr(tokens[3])); 284 setParameter(tokens[4], parseFloatStr(tokens[5])); 285 } 286 break; 287 } 288 } 289 setModelParameters(coordinatesArePrimitive); 290 applySymmetryAndSetTrajectory(); 291 if (totEnergy != null) 292 setEnergy(); 293 } 294 scalePrimitiveData(int i, float value)295 private void scalePrimitiveData(int i, float value) { 296 V3 v = V3.new3(primitiveData[i], primitiveData[i + 1], primitiveData[i + 2]); 297 v.normalize(); 298 v.scale(value); 299 primitiveData[i++] = v.x; 300 primitiveData[i++] = v.y; 301 primitiveData[i++] = v.z; 302 } 303 304 @Override applySymmetryAndSetTrajectory()305 public void applySymmetryAndSetTrajectory() throws Exception { 306 if (coordinatesArePrimitive && iHaveUnitCell && doCheckUnitCell && primitiveData != null && !isPrimitive) { 307 setModelParameters(false); 308 SymmetryInterface symFull = symmetry; 309 setModelParameters(true); 310 // Full cell -- must convert primitive to conventional 311 312 Atom[] atoms = asc.atoms; 313 int i0 = asc.getLastAtomSetAtomIndex(); 314 int i1 = asc.ac; 315 for (int i = i0; i < i1; i++) { 316 Atom atom = atoms[i]; 317 symmetry.toCartesian(atom, true); 318 symFull.toFractional(atom, true); 319 if (fixJavaFloat) 320 PT.fixPtFloats(atom, PT.FRACTIONAL_PRECISION); 321 } 322 setModelParameters(false); 323 } 324 applySymTrajASCR(); 325 } 326 /* Fractional coordinates of asymmetric unit : 327 328 -------------------------------------------------------------------------------- 329 No. Atomic x y z Charge Occupancy 330 Label (Frac) (Frac) (Frac) (e) (Frac) 331 -------------------------------------------------------------------------------- 332 1 Ba c 0.572390 0.144780 0.144780 0.1690 1.000000 333 2 Ba c 0.677610 * 0.355219 * 0.355219 * 0.1690 1.000000 334 3 Fe c 0.822391 * 0.644781 * 0.644781 * 1.9710 1.000000 335 4 Fe c 0.927610 * 0.855219 * 0.855219 * 1.9710 1.000000 336 5 Fe c 0.072390 * 0.144780 * 0.144780 * 1.9710 1.000000 337 6 Fe c 0.177610 * 0.355219 * 0.355219 * 1.9710 1.000000 338 7 Fe c 0.322391 * 0.644781 * 0.644781 * 1.9710 1.000000 339 8 Fe c 0.427610 * 0.855219 * 0.855219 * 1.9710 1.000000 340 9 O c 0.778081 * 0.943838 * 0.250000 * 0.5130 1.000000 341 */ 342 343 /* 344 345 Cartesian coordinates of cluster : 346 347 -------------------------------------------------------------------------------- 348 No. Atomic x y z Charge Occupancy 349 Label (Angs) (Angs) (Angs) (e) (Frac) 350 -------------------------------------------------------------------------------- 351 1 C1 c -2.3420 * 1.0960 * -0.0010 * -0.1819 1.000000 352 2 C1 c -1.0490 * 1.5890 * 0.0010 * -0.2684 1.000000 353 3 C2 c 0.0730 0.7330 -0.0010 0.0661 1.000000 354 4 C3 c -0.1850 * -0.6550 * -0.0030 * 0.1492 1.000000 355 5 C1 c -1.4790 * -1.1720 * -0.0050 * -0.1988 1.000000 356 6 C4 c -2.5770 * -0.3000 * -0.0040 * 0.2336 1.000000 357 7 H1 c -3.1860 * 1.7820 * -0.0070 * 0.1768 1.000000 358 8 H1 c -0.8980 * 2.6650 * 0.0030 * 0.1892 1.000000 359 9 C5 c 1.4540 * 1.1720 * 0.0010 * 0.1307 1.000000 360 10 H1 c -1.6110 * -2.2500 * -0.0110 * 0.1768 1.000000 361 11 C7 c 2.1730 * -1.1940 * 0.0010 * 0.7946 1.000000 362 12 C6 c 2.4410 * 0.2330 * 0.0020 * -0.4195 1.000000 363 13 H1 c 3.4910 * 0.5060 * 0.0030 * 0.1718 1.000000 364 14 O1 c 0.8280 * -1.5660 * -0.0020 * -0.3694 1.000000 365 15 O2 c 2.9950 * -2.0810 * 0.0040 * -0.6874 1.000000 366 16 C8 c 1.7830 * 2.6420 * 0.0020 * -0.2362 1.000000 367 17 H2 c 1.3640 * 3.1400 * -0.8810 * 0.0960 1.000000 368 18 H2 c 1.3630 * 3.1390 * 0.8850 * 0.0960 1.000000 369 19 H2 c 2.8640 * 2.7990 * 0.0020 * 0.0960 1.000000 370 20 N1 c -3.8700 * -0.7920 * -0.0560 * -0.8553 1.000000 371 21 H3 c -4.0140 * -1.7570 * 0.2000 * 0.4202 1.000000 372 22 H3 c -4.6160 * -0.1740 * 0.2260 * 0.4202 1.000000 373 -------------------------------------------------------------------------------- 374 375 */ 376 readAtomicPos(boolean finalizeSymmetry)377 private void readAtomicPos(boolean finalizeSymmetry) throws Exception { 378 newAtomSet(finalizeSymmetry); 379 discardLinesUntilContains(sep); 380 discardLinesUntilContains(sep); 381 while (rd() != null) { 382 if (line.indexOf(sep) >= 0 && rd().indexOf("Region") < 0) 383 break; 384 if (line.indexOf("Region") >= 0) { 385 rd(); 386 continue; 387 } 388 line = line.replace('*', ' '); 389 String[] tokens = getTokens(); 390 if (tokens[2].equals("c")) 391 addAtomXYZSymName(tokens, 3, null, tokens[1]); 392 } 393 if (finalizeSymmetry) 394 applySymmetryAndSetTrajectory(); 395 } 396 397 /* 398 399 Species output for all configurations : 400 401 -------------------------------------------------------------------------------- 402 Species Type Atomic Atomic Charge Radii (Angs) Library 403 Number Mass (e) Cova Ionic VDW Symbol 404 -------------------------------------------------------------------------------- 405 Fe Core 26 55.85 1.9710 1.340 0.000 2.000 406 Fe Shell 26 55.85 1.0290 1.340 0.000 2.000 407 Ba Core 56 137.33 0.1690 1.340 0.000 2.000 408 Ba Shell 56 137.33 1.8310 1.340 0.000 2.000 409 O Core 8 16.00 0.5130 0.730 0.000 1.360 410 O Shell 8 16.00 -2.5130 0.730 0.000 1.360 411 -------------------------------------------------------------------------------- 412 */ readPartialCharges()413 private void readPartialCharges() throws Exception { 414 atomCharges = new Hashtable<String, Float>(); 415 discardLinesUntilContains(sep); 416 discardLinesUntilContains(sep); 417 String[] tokens; 418 while ((tokens = PT.getTokens(rd())).length > 5) { 419 String species = tokens[0]; 420 Float charge = atomCharges.get(species); 421 float f = (charge == null ? 0 : charge.floatValue()); 422 atomCharges.put(species, Float.valueOf((f + parseFloatStr(tokens[4])))); 423 } 424 } 425 426 427 /* 428 -------------------------------------------------------------------------------- 429 Total lattice energy = -386.17106576 eV 430 -------------------------------------------------------------------------------- 431 Total lattice energy = -37259.6047 kJ/(mole unit cells) 432 Total lattice energy = -8905.2592 kcal/(mole unit cells) 433 -------------------------------------------------------------------------------- 434 435 436 **** Optimisation achieved **** 437 438 439 Final energy = -557.53367977 eV 440 Final Gnorm = 0.00035566 441 442 */ 443 444 private Double totEnergy; 445 private String energyUnits; 446 readEnergy()447 private void readEnergy() throws Exception { 448 // question: Why read monopole-monopole energy as "totEnergy"? 449 // note that in some cases this is in Kcal/mol 450 if (line.indexOf("=") < 0) 451 discardLinesUntilContains("="); 452 String[] tokens = PT.getTokens(line.substring(line.indexOf("="))); 453 totEnergy = Double.valueOf(Double.parseDouble(tokens[1])); 454 energyUnits = tokens[2]; 455 discardLinesUntilContains(sep); 456 } 457 setEnergy()458 private void setEnergy() { 459 asc.setAtomSetEnergy("" + totEnergy, totEnergy.floatValue()); 460 asc.setInfo("Energy", totEnergy); 461 asc.setAtomSetName("E = " + totEnergy + " " + energyUnits); 462 totEnergy = null; 463 } 464 465 /* 466 private void setAtomSetInfo(){ 467 468 }*/ 469 470 /* 471 Frequency 0.0000 0.0000 0.0000 133.9963 134.1490 227.2719 472 IR Intensity 0.0000 0.0000 0.0000 0.0030 0.0030 0.0000 473 in X 0.0000 0.0000 0.0000 0.0008 0.0023 0.0000 474 in Y 0.0000 0.0000 0.0000 0.0023 0.0008 0.0000 475 in Z 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 476 Raman Intsty 0.0000 0.0000 0.0000 0.0000 0.0000 0.3276 477 478 1 x 0.390326 -0.056775 -0.015611 -0.039949 0.069193 -0.305899 479 1 y 0.054908 0.388754 -0.040954 0.091620 0.052898 0.000000 480 1 z -0.021265 -0.038324 -0.392302 -0.013870 -0.008008 0.000000 481 2 x 0.390326 -0.056775 -0.015611 0.099320 -0.011213 0.152950 482 2 y 0.054908 0.388754 -0.040954 0.011214 -0.086372 -0.264916 483 2 z -0.021265 -0.038324 -0.392302 0.013870 -0.008007 0.000000 484 3 x 0.390326 -0.056775 -0.015611 -0.039948 -0.091621 0.152949 485 3 y 0.054908 0.388754 -0.040954 -0.069194 0.052897 0.264916 486 3 z -0.021265 -0.038324 -0.392302 0.000000 0.016016 0.000000 487 4 x 0.294586 -0.042849 -0.011782 -0.120330 0.243402 -0.030861 488 4 y 0.041440 0.293399 -0.030909 0.086109 0.358453 0.294237 489 4 z -0.016049 -0.028924 -0.296077 0.125378 0.314751 0.179865 490 5 x 0.294586 -0.042849 -0.011782 0.381440 0.203588 -0.239386 491 5 y 0.041440 0.293399 -0.030909 0.046295 -0.143317 -0.173845 492 5 z -0.016049 -0.028924 -0.296077 -0.335271 -0.048796 0.179865 493 6 x 0.294586 -0.042849 -0.011782 0.096075 -0.211051 0.270247 494 6 y 0.041440 0.293399 -0.030909 -0.368344 0.142048 -0.120392 495 6 z -0.016049 -0.028924 -0.296077 0.209894 -0.265956 0.179865 496 7 x 0.294586 -0.042849 -0.011782 -0.150630 0.225910 -0.030862 497 7 y 0.041440 0.293399 -0.030909 0.353487 -0.104649 -0.294238 498 7 z -0.016049 -0.028924 -0.296077 0.335275 -0.048791 -0.179865 499 8 x 0.294586 -0.042849 -0.011782 -0.367031 -0.228549 -0.239387 500 8 y 0.041440 0.293399 -0.030909 -0.100971 0.111751 0.173846 501 8 z -0.016049 -0.028924 -0.296077 -0.209891 -0.265961 -0.179865 502 9 x 0.294586 -0.042849 -0.011782 0.134742 -0.188728 0.270249 503 9 y 0.041440 0.293399 -0.030909 -0.061150 -0.390021 0.120392 504 9 z -0.016049 -0.028924 -0.296077 -0.125383 0.314752 -0.179865 505 */ 506 /* 507 private void readFrequency() throws Exception { 508 discardLines(1); 509 readLine(); 510 511 //for the time being it reads only phonons at the gamma point k (0, 0, 0) 512 if (line.contains("gamma")) { 513 discardLines(9); 514 if (line.contains("Frequency")) 515 ; 516 } else { 517 return; 518 } 519 } 520 */ 521 } 522