1 /* $RCSfile$ 2 * $Author: hansonr $ 3 * $Date: 2006-05-13 19:17:06 -0500 (Sat, 13 May 2006) $ 4 * $Revision: 5114 $ 5 * 6 * Copyright (C) 2003-2005 Miguel, Jmol Development, www.jmol.org 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 package org.jmol.quantum; 25 26 import java.io.BufferedReader; 27 import java.util.Hashtable; 28 import java.util.Map; 29 import java.util.Map.Entry; 30 31 import javajs.util.BS; 32 import javajs.util.Lst; 33 import javajs.util.Measure; 34 import javajs.util.PT; 35 import javajs.util.SB; 36 import javajs.util.V3; 37 38 import org.jmol.api.JmolNMRInterface; 39 import org.jmol.modelset.Atom; 40 import org.jmol.modelset.Bond; 41 import org.jmol.modelset.MeasurementData; 42 import org.jmol.modelset.Model; 43 import org.jmol.quantum.NMRNoeMatrix.NOEParams; 44 import org.jmol.util.Edge; 45 import org.jmol.util.Escape; 46 import org.jmol.util.Logger; 47 import org.jmol.util.Tensor; 48 import org.jmol.viewer.FileManager; 49 import org.jmol.viewer.Viewer; 50 51 /* 52 * 53 * Bob Hanson hansonr@stolaf.edu 7/4/2013 54 * 55 * Added NOE and solution-phase J calcs 11/2019 56 * 57 */ 58 59 /* 60 * NOTE -- THIS CLASS IS INSTANTIATED USING Interface.getOptionInterface 61 * NOT DIRECTLY -- FOR MODULARIZATION. NEVER USE THE CONSTRUCTOR DIRECTLY! 62 * 63 */ 64 65 public class NMRCalculation implements JmolNMRInterface { 66 67 private static final int MAGNETOGYRIC_RATIO = 1; 68 private static final int QUADRUPOLE_MOMENT = 2; 69 private static final double e_charge = 1.60217646e-19; //C 70 private static final double h_planck = 6.62606957e-34; //J*s 71 private static final double h_bar_planck = h_planck / (2 * Math.PI); 72 private static final double DIPOLAR_FACTOR = h_bar_planck * 1E37; // 1e37 = (1/1e-10)^3 * 1e7 * 1e7 / 1e-7 73 private static final double J_FACTOR = h_bar_planck / (2 * Math.PI) * 1E33; 74 private static final double Q_FACTOR = e_charge * (9.71736e-7) / h_planck; 75 76 public static final int MODE_CALC_INVALID = 0; 77 public final static int MODE_CALC_2JHH = 1; 78 public final static int MODE_CALC_3JHH = 2; 79 public final static int MODE_CALC_JHH = 3; 80 public final static int MODE_CALC_3JCH = 4; 81 public final static int MODE_CALC_J = 7; 82 public final static int MODE_CALC_NOE = 8; 83 public static final int MODE_CALC_ALL = 0xF; 84 85 86 87 private Viewer vwr; 88 NMRCalculation()89 public NMRCalculation() { 90 } 91 92 @Override setViewer(Viewer vwr)93 public JmolNMRInterface setViewer(Viewer vwr) { 94 this.vwr = vwr; 95 getData(); 96 return this; 97 } 98 99 @Override getQuadrupolarConstant(Tensor efg)100 public float getQuadrupolarConstant(Tensor efg) { 101 if (efg == null) 102 return 0; 103 Atom a = vwr.ms.at[efg.atomIndex1]; 104 return (float) (getIsotopeData(a, QUADRUPOLE_MOMENT) * efg.eigenValues[2] * Q_FACTOR); 105 } 106 107 /** 108 * Returns a list of tensors that are of the specified type and have both 109 * atomIndex1 and atomIndex2 in bsA. If there is just one atom specified, then 110 * the list is "all tensors involving this atom". 111 * 112 * We have to use atom sites, because interaction tensors are not duplicated. 113 * 114 * @param type 115 * @param bsA 116 * @return list of Tensors 117 */ 118 @SuppressWarnings("unchecked") getInteractionTensorList(String type, BS bsA)119 private Lst<Tensor> getInteractionTensorList(String type, BS bsA) { 120 if (type != null) 121 type = type.toLowerCase(); 122 BS bsModels = vwr.ms.getModelBS(bsA, false); 123 BS bs1 = getAtomSiteBS(bsA); 124 int iAtom = (bs1.cardinality() == 1 ? bs1.nextSetBit(0) : -1); 125 Lst<Tensor> list = new Lst<Tensor>(); 126 for (int i = bsModels.nextSetBit(0); i >= 0; i = bsModels.nextSetBit(i + 1)) { 127 Lst<Tensor> tensors = (Lst<Tensor>) vwr.ms.getInfo(i, 128 "interactionTensors"); 129 if (tensors == null) 130 continue; 131 int n = tensors.size(); 132 for (int j = 0; j < n; j++) { 133 Tensor t = tensors.get(j); 134 if (type == null || t.type.equals(type) && t.isSelected(bs1, iAtom)) 135 list.addLast(t); 136 } 137 } 138 return list; 139 } 140 141 /** 142 * Interaction tensors are not repeated for every possible combination. They 143 * are just for the base atom set. These are identified as a.atomIndex == 144 * models[b.modelIndex].firstAtomIndex + b.atomSite - 1 145 * 146 * @param bsA 147 * @return new bs in terms of atom sites 148 */ getAtomSiteBS(BS bsA)149 private BS getAtomSiteBS(BS bsA) { 150 if (bsA == null) 151 return null; 152 BS bs = new BS(); 153 Atom[] atoms = vwr.ms.at; 154 Model[] models = vwr.ms.am; 155 156 for (int i = bsA.nextSetBit(0); i >= 0; i = bsA.nextSetBit(i + 1)) { 157 if (!bsA.get(i)) 158 continue; 159 Atom a = atoms[i]; 160 bs.set(models[a.mi].firstAtomIndex - 1 + a.atomSite); 161 } 162 return bs; 163 } 164 165 // private int getOtherAtom(Tensor t, int iAtom) { 166 // return (t.atomIndex1 == iAtom ? t.atomIndex2 : t.atomIndex1); 167 // } 168 169 @Override getUniqueTensorSet(BS bsAtoms)170 public BS getUniqueTensorSet(BS bsAtoms) { 171 BS bs = new BS(); 172 Atom[] atoms = vwr.ms.at; 173 for (int i = vwr.ms.mc; --i >= 0;) { 174 BS bsModelAtoms = vwr.getModelUndeletedAtomsBitSet(i); 175 bsModelAtoms.and(bsAtoms); 176 // exclude any models without symmetry 177 if (vwr.ms.getUnitCell(i) == null) 178 continue; 179 // exclude any symmetry- 180 for (int j = bsModelAtoms.nextSetBit(0); j >= 0; j = bsModelAtoms 181 .nextSetBit(j + 1)) 182 if (atoms[j].atomSite != atoms[j].i + 1) 183 bsModelAtoms.clear(j); 184 bs.or(bsModelAtoms); 185 // march through all the atoms in the model... 186 for (int j = bsModelAtoms.nextSetBit(0); j >= 0; j = bsModelAtoms 187 .nextSetBit(j + 1)) { 188 Object[] ta = atoms[j].getTensors(); 189 if (ta == null) 190 continue; 191 // go through all this atom's tensors... 192 for (int jj = ta.length; --jj >= 0;) { 193 Tensor t = (Tensor) ta[jj]; 194 if (t == null) 195 continue; 196 // for each tensor in A, go through all atoms after the first-selected one... 197 for (int k = bsModelAtoms.nextSetBit(j + 1); k >= 0; k = bsModelAtoms 198 .nextSetBit(k + 1)) { 199 Object[] tb = atoms[k].getTensors(); 200 if (tb == null) 201 continue; 202 // for each tensor in B, go through all this atom's tensors... 203 for (int kk = tb.length; --kk >= 0;) { 204 // if equivalent, reject it. 205 if (t.isEquiv((Tensor) tb[kk])) { 206 bsModelAtoms.clear(k); 207 bs.clear(k); 208 break; 209 } 210 } 211 } 212 } 213 } 214 } 215 return bs; 216 } 217 getJCouplingHz(Atom a1, Atom a2, String type, Tensor isc)218 public float getJCouplingHz(Atom a1, Atom a2, String type, Tensor isc) { 219 return getIsoOrAnisoHz(true, a1, a2, type, isc); 220 } 221 222 @Override getIsoOrAnisoHz(boolean isIso, Atom a1, Atom a2, String units, Tensor isc)223 public float getIsoOrAnisoHz(boolean isIso, Atom a1, Atom a2, String units, 224 Tensor isc) { 225 if (isc == null) { 226 String type = getISCtype(a1, units); 227 if (type == null || a1.mi != a2.mi) { 228 if (!units.equals("hz")) 229 return 0; 230 double[] data = calc2or3JorNOE(vwr, new Atom[]{a1, null, null, a2}, null, MODE_CALC_JHH); // H-H only here; no NOE, no JCH 231 return (data == null ? Float.NaN : (float) data[1]); 232 } 233 BS bs = new BS(); 234 bs.set(a1.i); 235 bs.set(a2.i); 236 Lst<Tensor> list = getInteractionTensorList(type, bs); 237 if (list.size() == 0) 238 return Float.NaN; 239 isc = list.get(0); 240 } else { 241 a1 = vwr.ms.at[isc.atomIndex1]; 242 a2 = vwr.ms.at[isc.atomIndex2]; 243 } 244 return (float) (getIsotopeData(a1, MAGNETOGYRIC_RATIO) 245 * getIsotopeData(a2, MAGNETOGYRIC_RATIO) 246 * (isIso ? isc.isotropy() : isc.anisotropy()) * J_FACTOR); 247 } 248 249 @SuppressWarnings("unchecked") getISCtype(Atom a1, String type)250 private String getISCtype(Atom a1, String type) { 251 Lst<Tensor> tensors = (Lst<Tensor>) vwr.ms.getInfo(a1.mi, 252 "interactionTensors"); 253 if (tensors == null) 254 return null; 255 type = (type == null ? "" : type.toLowerCase()); 256 int pt = -1; 257 if ((pt = type.indexOf("_hz")) >= 0 || (pt = type.indexOf("_khz")) >= 0 258 || (pt = type.indexOf("hz")) >= 0 || (pt = type.indexOf("khz")) >= 0) 259 type = type.substring(0, pt); 260 if (type.length() == 0) 261 type = "isc"; 262 return type; 263 } 264 265 @Override getDipolarConstantHz(Atom a1, Atom a2)266 public float getDipolarConstantHz(Atom a1, Atom a2) { 267 if (Logger.debugging) 268 Logger.debug(a1 + " g=" + getIsotopeData(a1, MAGNETOGYRIC_RATIO) + "; " 269 + a2 + " g=" + getIsotopeData(a2, MAGNETOGYRIC_RATIO)); 270 float v = (float) (-getIsotopeData(a1, MAGNETOGYRIC_RATIO) 271 * getIsotopeData(a2, MAGNETOGYRIC_RATIO) / Math.pow(a1.distance(a2), 3) * DIPOLAR_FACTOR); 272 return (v == 0 || a1 == a2 ? Float.NaN : v); 273 } 274 275 @Override getDipolarCouplingHz(Atom a1, Atom a2, V3 vField)276 public float getDipolarCouplingHz(Atom a1, Atom a2, V3 vField) { 277 V3 v12 = V3.newVsub(a2, a1); 278 double r = v12.length(); 279 double costheta = v12.dot(vField) / r / vField.length(); 280 return (float) (getDipolarConstantHz(a1, a2) * (3 * costheta - 1) / 2); 281 } 282 283 /** 284 * isotopeData keyed by nnnSym, for example: 1H, 19F, etc.; and also by 285 * element name itself: H, F, etc., for default 286 */ 287 private Map<String, double[]> isotopeData; 288 289 /** 290 * NOTE! Do not change this txt file! Instead, edit 291 * trunk/Jmol/_documents/nmr_data.xls and then clip its contents to 292 * org/jmol/quantum/nmr_data.txt. 293 * 294 */ 295 private final static String resource = "nmr_data.txt"; 296 297 /** 298 * Get magnetogyricRatio (gamma/10^7 rad s^-1 T^-1) and quadrupoleMoment 299 * (Q/10^-2 fm^2) for a given isotope or for the default isotope of an 300 * element. 301 * 302 * @param a 303 * 304 * @param iType 305 * @return g or Q 306 */ getIsotopeData(Atom a, int iType)307 private double getIsotopeData(Atom a, int iType) { 308 int iso = a.getIsotopeNumber(); 309 String sym = a.getElementSymbolIso(false); 310 double[] d = isotopeData.get(iso == 0 ? sym : "" + iso + sym); 311 return (d == null ? 0 : d[iType]); 312 } 313 314 /** 315 * Creates the data set necessary for doing NMR calculations. Values are 316 * retrievable using getProperty "nmrInfo" "Xx"; each entry is 317 * float[+/-isotopeNumber, g, Q], where [0] < 0 for the default value. 318 * 319 */ 320 @SuppressWarnings("resource") 321 // Resource leak: 'br' is not closed at this location -- Nonsense -- there's a finally{} block getData()322 private void getData() { 323 BufferedReader br = null; 324 try { 325 boolean debugging = Logger.debugging; 326 br = FileManager.getBufferedReaderForResource(vwr, this, 327 "org/jmol/quantum/", resource); 328 // #extracted by Simone Sturniolo from ROBIN K. HARRIS, EDWIN D. BECKER, SONIA M. CABRAL DE MENEZES, ROBIN GOODFELLOW, AND PIERRE GRANGER, Pure Appl. Chem., Vol. 73, No. 11, pp. 1795–1818, 2001. NMR NOMENCLATURE. NUCLEAR SPIN PROPERTIES AND CONVENTIONS FOR CHEMICAL SHIFTS (IUPAC Recommendations 2001) 329 // #element atomNo isotopeDef isotope1 G1 Q1 isotope2 G2 Q2 isotope3 G3 Q3 330 // H 1 1 1 26.7522128 0 2 4.10662791 0.00286 3 28.5349779 0 331 isotopeData = new Hashtable<String, double[]>(); 332 String line; 333 while ((line = br.readLine()) != null) { 334 if (debugging) 335 Logger.info(line); 336 if (line.indexOf("#") >= 0) 337 continue; 338 String[] tokens = PT.getTokens(line); 339 String name = tokens[0]; 340 String defaultIso = tokens[2] + name; 341 if (debugging) 342 Logger.info(name + " default isotope " + defaultIso); 343 for (int i = 3; i < tokens.length; i += 3) { 344 int n = Integer.parseInt(tokens[i]); 345 String isoname = n + name; 346 double[] dataGQ = new double[] { n, 347 Double.parseDouble(tokens[i + 1]), 348 Double.parseDouble(tokens[i + 2]) }; 349 if (debugging) 350 Logger.info(isoname + " " + Escape.eAD(dataGQ)); 351 isotopeData.put(isoname, dataGQ); 352 } 353 double[] defdata = isotopeData.get(defaultIso); 354 if (defdata == null) { 355 Logger.error("Cannot find default NMR data in nmr_data.txt for " 356 + defaultIso); 357 throw new NullPointerException(); 358 } 359 defdata[0] = -defdata[0]; 360 isotopeData.put(name, defdata); 361 } 362 } catch (Exception e) { 363 Logger.error("Exception " + e.toString() + " reading " + resource); 364 } finally { 365 try { 366 br.close(); 367 } catch (Exception ee) { 368 // ignore 369 } 370 } 371 } 372 373 @Override getInfo(String what)374 public Object getInfo(String what) { 375 if (what.equals("all")) { 376 Map<String, Object> map = new Hashtable<String, Object>(); 377 map.put("isotopes", isotopeData); 378 map.put("shiftRefsPPM", shiftRefsPPM); 379 return map; 380 } 381 if (PT.isDigit(what.charAt(0))) 382 return isotopeData.get(what); 383 Lst<Object> info = new Lst<Object>(); 384 for (Entry<String, double[]> e : isotopeData.entrySet()) { 385 String key = e.getKey(); 386 if (PT.isDigit(key.charAt(0)) && key.endsWith(what)) 387 info.addLast(e.getValue()); 388 } 389 return info; 390 } 391 392 @Override getChemicalShift(Atom atom)393 public float getChemicalShift(Atom atom) { 394 float v = getMagneticShielding(atom); 395 if (Float.isNaN(v)) 396 return v; 397 Float ref = shiftRefsPPM.get(atom.getElementSymbol()); 398 return (ref == null ? 0 : ref.floatValue()) - v; 399 } 400 401 @Override getMagneticShielding(Atom atom)402 public float getMagneticShielding(Atom atom) { 403 Tensor t = vwr.ms.getAtomTensor(atom.i, "ms"); 404 return (t == null ? Float.NaN : t.isotropy()); 405 } 406 407 private Map<String, Float> shiftRefsPPM = new Hashtable<String, Float>(); 408 409 @Override getState(SB sb)410 public boolean getState(SB sb) { 411 if (shiftRefsPPM.isEmpty()) 412 return false; 413 for (Entry<String, Float> nuc : shiftRefsPPM.entrySet()) 414 sb.append(" set shift_").append(nuc.getKey()).append(" ") 415 .appendO(nuc.getValue()).append("\n"); 416 return true; 417 } 418 419 @Override setChemicalShiftReference(String element, float value)420 public boolean setChemicalShiftReference(String element, float value) { 421 if (element == null) { 422 shiftRefsPPM.clear(); 423 return false; 424 } 425 element = element.substring(0, 1).toUpperCase() + element.substring(1); 426 shiftRefsPPM.put(element, Float.valueOf(value)); 427 return true; 428 } 429 430 @Override getTensorInfo(String tensorType, String infoType, BS bs)431 public Lst<Object> getTensorInfo(String tensorType, String infoType, BS bs) { 432 if ("".equals(tensorType)) 433 tensorType = null; 434 infoType = (infoType == null ? ";all." : ";" + infoType + "."); 435 Lst<Object> data = new Lst<Object>(); 436 Lst<Object> list1; 437 if (";dc.".equals(infoType)) { 438 // tensorType is irrelevant for dipolar coupling constant 439 Atom[] atoms = vwr.ms.at; 440 for (int i = bs.nextSetBit(0); i >= 0; i = bs.nextSetBit(i + 1)) 441 for (int j = bs.nextSetBit(i + 1); j >= 0; j = bs.nextSetBit(j + 1)) { 442 list1 = new Lst<Object>(); 443 list1.addLast(Integer.valueOf(atoms[i].i)); 444 list1.addLast(Integer.valueOf(atoms[j].i)); 445 list1 446 .addLast(Float.valueOf(getDipolarConstantHz(atoms[i], atoms[j]))); 447 data.addLast(list1); 448 } 449 return data; 450 } 451 if (tensorType == null || tensorType.startsWith("isc")) { 452 boolean isJ = infoType.equals(";j."); 453 boolean isEta = infoType.equals(";eta."); 454 Lst<Tensor> list = getInteractionTensorList(tensorType, bs); 455 int n = (list == null ? 0 : list.size()); 456 for (int i = 0; i < n; i++) { 457 Tensor t = list.get(i); 458 list1 = new Lst<Object>(); 459 list1.addLast(Integer.valueOf(t.atomIndex1)); 460 list1.addLast(Integer.valueOf(t.atomIndex2)); 461 list1.addLast(isEta || isJ ? Float.valueOf(getIsoOrAnisoHz(isJ, null, 462 null, null, t)) : t.getInfo(infoType)); 463 data.addLast(list1); 464 } 465 if (tensorType != null) 466 return data; 467 } 468 boolean isChi = tensorType != null && tensorType.startsWith("efg") 469 && infoType.equals(";chi."); 470 boolean isFloat = (isChi || Tensor.isFloatInfo(infoType)); 471 for (int i = bs.nextSetBit(0); i >= 0; i = bs.nextSetBit(i + 1)) { 472 if (tensorType == null) { 473 Object[] a = vwr.ms.getAtomTensorList(i); 474 if (a != null) 475 for (int j = 0; j < a.length; j++) 476 data.addLast(((Tensor) a[j]).getInfo(infoType)); 477 } else { 478 Tensor t = vwr.ms.getAtomTensor(i, tensorType); 479 data.addLast(t == null ? (isFloat ? Float.valueOf(0) : "") 480 : isChi ? Float.valueOf(getQuadrupolarConstant(t)) : t 481 .getInfo(infoType)); 482 } 483 } 484 return data; 485 } 486 487 @Override getMinDistances(MeasurementData md)488 public Map<String, Integer> getMinDistances(MeasurementData md) { 489 BS bsPoints1 = (BS) md.points.get(0); 490 int n1 = bsPoints1.cardinality(); 491 if (n1 == 0 || !(md.points.get(1) instanceof BS)) 492 return null; 493 BS bsPoints2 = (BS) md.points.get(1); 494 int n2 = bsPoints2.cardinality(); 495 if (n1 < 2 && n2 < 2) 496 return null; 497 Map<String, Integer> htMin = new Hashtable<String, Integer>(); 498 Atom[] atoms = vwr.ms.at; 499 for (int i = bsPoints1.nextSetBit(0); i >= 0; i = bsPoints1 500 .nextSetBit(i + 1)) { 501 Atom a1 = atoms[i]; 502 String name = a1.getAtomName(); 503 for (int j = bsPoints2.nextSetBit(0); j >= 0; j = bsPoints2 504 .nextSetBit(j + 1)) { 505 Atom a2 = atoms[j]; 506 int d = (int) (a2.distanceSquared(a1) * 100); 507 if (d == 0) 508 continue; 509 String name1 = a2.getAtomName(); 510 String key = (name.compareTo(name1) < 0 ? name + name1 : name1 + name); 511 Integer min = htMin.get(key); 512 if (min == null) { 513 min = Integer.valueOf(d); 514 htMin.put(key, min); 515 continue; 516 } 517 if (d < min.intValue()) 518 htMin.put(key, Integer.valueOf(d)); 519 } 520 } 521 return htMin; 522 } 523 524 //// non-quantum calcs 525 526 /** 527 * Calculate an H-X-X-H or C-X-X-H coupling constant. 528 * 529 * Use Altoona equation (Tetrahedron 36, 2783-2792) 530 * 531 * If there are fewer than three substituents on each central atom, or if 532 * either central atom is not carbon, defaults to general Karplus equation. 533 * 534 * 535 */ 536 537 static Hashtable<String, Double> deltaElectro = new Hashtable<String, Double>(); 538 539 static { 540 // Values taken from www.spectroscopynow.com/Spy/tools/proton-proton.html website. Need to find original source 541 // 2.20 is eneg of H 542 double enegH = 2.20; 543 deltaElectro.put("C", new Double(2.60 - enegH)); 544 deltaElectro.put("O", new Double(3.50 - enegH)); 545 deltaElectro.put("N", new Double(3.05 - enegH)); 546 deltaElectro.put("F", new Double(3.90 - enegH)); 547 deltaElectro.put("Cl", new Double(3.15 - enegH)); 548 deltaElectro.put("Br", new Double(2.95 - enegH)); 549 deltaElectro.put("I", new Double(2.65 - enegH)); 550 deltaElectro.put("S", new Double(2.58 - enegH));// Pauling 551 deltaElectro.put("Si", new Double(1.90 - enegH));// Pauling 552 } 553 554 static double[][] pAltona = new double[5][8]; 555 static { 556 for (int nNonH = 0; nNonH < 5; nNonH++) { 557 double[] p = pAltona[nNonH]; 558 switch (nNonH) { 559 case 0: 560 case 1: 561 case 2: 562 // Janocchio - original 563 // p[1] = 13.89; 564 // p[2] = -0.98; 565 // p[3] = 0; 566 // p[4] = 1.02; 567 // p[5] = -3.40; 568 // p[6] = 14.9; 569 // p[7] = 0.24; 570 // but 571 // Haasnoot and Altona 1981 https://www.researchgate.net/publication/230009488_The_relationship_between_proton-proton_NMR_coupling_constants_and_substituent_electronegativities_II-conformational_analysis_of_the_sugar_ring_in_nucleosides_and_nucleotides_in_solution_using_a_genera 572 // April 1981Magnetic Resonance in Chemistry 15(1):43 - 52 573 // see also http://www.colby.edu/chemistry/NMR/scripts/altona/altonaref.html 574 // and https://www.spectroscopynow.com/userfiles/sepspec/file/specNOW/HTML%20files/proton-proton2.htm 575 p[1] = 13.7; 576 p[2] = -0.73; 577 p[3] = 0; 578 p[4] = 0.56; 579 p[5] = -2.47; 580 p[6] = 16.9; 581 p[7] = 0.14; 582 break; 583 case 3: 584 p[1] = 13.22; 585 p[2] = -0.99; 586 p[3] = 0; 587 p[4] = 0.87; 588 p[5] = -2.46; 589 p[6] = 19.9; 590 p[7] = 0; 591 break; 592 case 4: 593 p[1] = 13.24; 594 p[2] = -0.91; 595 p[3] = 0; 596 p[4] = 0.53; 597 p[5] = -2.41; 598 p[6] = 15.5; 599 p[7] = 0.19; 600 break; 601 } 602 p[6] = p[6] * Math.PI / 180.0; 603 604 } 605 } 606 calcJKarplus(double theta)607 public static double calcJKarplus(double theta) { 608 // Simple Karplus equation for 3JHH, ignoring differences in C-substituents 609 double j0 = 8.5; 610 double j180 = 9.5; 611 double jconst = 0.28; 612 double cos = Math.cos(theta); 613 double jab = 0; 614 if (theta >= -90.0 && theta < 90.0) { 615 jab = j0 * cos * cos - jconst; 616 } else { 617 jab = j180 * cos * cos - jconst; 618 } 619 620 return jab; 621 } 622 getInitialJValue(int nNonH, double theta)623 private static double getInitialJValue(int nNonH, double theta) { 624 double[] p = pAltona[nNonH]; 625 double cos = Math.cos(theta); 626 return p[1] * cos * cos + p[2] * cos + p[3]; 627 } 628 getIncrementalJValue(int nNonH, String element, V3 sA_cA, V3 v21, V3 v23, double theta, int f)629 private static double getIncrementalJValue(int nNonH, String element, 630 V3 sA_cA, V3 v21, V3 v23, 631 double theta, int f) { 632 if (nNonH < 0 || nNonH > 5) 633 return 0; 634 Double de = deltaElectro.get(element); 635 if (de == null) 636 return 0; 637 double e = de.doubleValue(); 638 // f here accounts for still using v23 instead of v32 for cB. 639 int sign = getSubSign(sA_cA, v21, v23, f); 640 double[] p = pAltona[nNonH]; 641 double cos = Math.cos(sign * theta + p[6] * Math.abs(e)); 642 return e * (p[4] + p[5] * cos * cos); 643 644 // but see also http://www.colby.edu/chemistry/NMR/scripts/altona/altonaref.html 645 646 // note that Dave Evans did not add p7. 647 648 } 649 650 /** 651 * Look for sign of (v23 x v21).dot.(sA_cA). 652 * 653 * But note that for the second 654 * carbon, we must reverse this. 655 * 656 * @param sA_cA C to sub 657 * @param v21 C to H 658 * @param v23 C to other C 659 * @param f 1 for carbon A; -1 for carbon B 660 * @return f or -f (+1 or -1) 661 */ getSubSign(V3 sA_cA, V3 v21, V3 v23, int f)662 private static int getSubSign(V3 sA_cA, V3 v21, V3 v23, int f) { 663 V3 cross = new V3(); 664 cross.cross(v23, v21); 665 return (cross.dot(sA_cA) > 0 ? f : -f); 666 } 667 668 /** 669 * 670 * 671 * @param subElements 672 * int[2][3] with element names 673 * @param subVectors 674 * V3[2][4] with vectors TO these substituents from their respective 675 * centers 676 * @param v21 677 * vector from cA to hA 678 * @param v34 679 * vector from cB to hB 680 * @param v23 681 * vector from cA to cB 682 * @param theta 683 * dihedral angle hA-cA-cB-hB 684 * @param is23Double 685 * @return estimated coupling constant 686 */ calc3JHHOnly(String[][] subElements, V3[][] subVectors, V3 v21, V3 v34, V3 v23, double theta, boolean is23Double)687 private static double calc3JHHOnly(String[][] subElements, V3[][] subVectors, 688 V3 v21, V3 v34, V3 v23, double theta, 689 boolean is23Double) { 690 691 // Substituents on atoms A and B 692 // Check number of substituents 693 694 int nNonH = 0; 695 696 // Count substituents of carbons A and B 697 for (int i = 0, n = (is23Double ? 2 : 3); i < n; i++) { 698 if (!subElements[0][i].equals("H")) { 699 nNonH++; 700 } 701 if (!subElements[1][i].equals("H")) { 702 nNonH++; 703 } 704 } 705 706 double jvalue = getInitialJValue(nNonH, theta); 707 708 for (int i = 0, n = (is23Double ? 2 : 3); i < n; i++) { 709 String element = subElements[0][i]; 710 if (!element.equals("H")) { 711 jvalue += getIncrementalJValue(nNonH, element, subVectors[0][i], v21, 712 v23, theta, 1); 713 } 714 element = subElements[1][i]; 715 if (!element.equals("H")) { 716 jvalue += getIncrementalJValue(nNonH, element, subVectors[1][i], v34, 717 v23, theta, -1); 718 } 719 } 720 if (is23Double) { 721 if (Math.abs(theta) < Math.PI/2) 722 jvalue *= 0.75; // BH my simple approximation; just a guess to get us in the ballpark 723 else 724 jvalue *= 1.33; // just a guess 725 } 726 return jvalue; 727 } 728 729 // static { 730 // for (int n = 0; n < 360; n+= 10) 731 // System.out.println(n + "\t" + getInitialJValue(2, n*Math.PI/180)); 732 // System.out.println("NMRCalc???"); 733 // } 734 final public static String JCH3_NONE = "none"; 735 final public static String JCH3_WASYLISHEN_SCHAEFER = "was"; 736 final public static String JCH3_TVAROSKA_TARAVEL = "tva"; 737 final public static String JCH3_AYDIN_GUETHER = "ayd"; 738 739 /** 740 * 741 * @param CHequation 742 * 743 * 'was' Simple equation for 3JCH, from Wasylishen and Schaefer Can J 744 * Chem (1973) 51 961 used in Kozerski et al. J Chem Soc Perkin 2, 745 * (1997) 1811 746 * 747 * 'tva' Tvaroska and Taravel Adv. Carbohydrate Chem. Biochem. (1995) 748 * 51, 15-61 749 * 750 * 'ayd' Aydin and Guether Mag. Res. Chem. (1990) 28, 448-457 751 * 752 * @param theta 753 * dihedral 754 * @param is23Double 755 * 756 * @return 3JCH prediction 757 */ calc3JCH(String CHequation, double theta, boolean is23Double)758 public static double calc3JCH(String CHequation, double theta, 759 boolean is23Double) { 760 761 if (CHequation == null) 762 CHequation = "was"; 763 if (!PT.isOneOf(CHequation, ";was;tva;ayd;")) 764 throw new IllegalArgumentException( 765 "Equation must be one of was, tva, or ayd"); 766 767 if (CHequation.equals("was")) { 768 769 final double A = 3.56; 770 final double C = 4.26; 771 772 double j = A * Math.cos(2 * theta) - Math.cos(theta) + C; 773 return j; 774 775 } else if (CHequation.equals("tva")) { 776 double j = 4.5 - 0.87 * Math.cos(theta) + Math.cos(2.0 * theta); 777 return j; 778 } else if (CHequation.equals("ayd")) { 779 double j = 5.8 * Math.pow(Math.cos(theta), 2) - 1.6 * Math.cos(theta) 780 + 0.28 * Math.sin(2.0 * theta) - 0.02 * Math.sin(theta) + 0.52; 781 return j; 782 } else { 783 return 0.0; 784 } 785 } 786 calcNOE(Viewer viewer, Atom atom1, Atom atom2)787 public static double[] calcNOE(Viewer viewer, Atom atom1, Atom atom2) { 788 return calc2or3JorNOE(viewer, new Atom[] {atom1, null, null, atom2}, null, MODE_CALC_J); 789 } 790 791 792 /** 793 * Calculate a 2-bond (geminal) or 3-bond (vicinal) coupling constant or an 794 * NOE; 795 * @param viewer 796 * 797 * @param atoms 798 * required Atom[4]; can be just two atoms, then in atom[0] and atom[4] 799 * @param CHEquation 800 * 'none' or 'was' or 'tva' or 'ayd' 801 * @param mode 802 * 803 * @return [theta, jvalue, atom2.i, atom3.i] for 3JHH; [theta, jvalue, 804 * center.i] for 2JHH; [distance, noe] for NOE 805 */ calc2or3JorNOE(Viewer viewer, Atom[] atoms, String CHEquation, int mode)806 public static double[] calc2or3JorNOE(Viewer viewer, Atom[] atoms, String CHEquation, 807 int mode) { 808 809 if (CHEquation == null || CHEquation.equals("none")) 810 mode &= ~MODE_CALC_3JCH; 811 String[] elements = new String[4]; 812 mode = getCalcType(atoms, elements, mode); 813 switch (mode) { 814 default: 815 case MODE_CALC_INVALID: 816 return null; 817 case MODE_CALC_NOE: 818 return calcNOEImpl(viewer, atoms[0], atoms[3]); 819 case MODE_CALC_2JHH: 820 return calc2JHH(atoms[0], atoms[1], atoms[3]); 821 case MODE_CALC_3JCH: 822 case MODE_CALC_3JHH: 823 break; 824 } 825 String[][] subElements = new String[2][3]; 826 V3[][] subVectors = new V3[2][3]; 827 828 V3 v23 = V3.newVsub(atoms[2], atoms[1]); 829 V3 v21 = V3.newVsub(atoms[0], atoms[1]); 830 V3 v34 = V3.newVsub(atoms[3], atoms[2]); 831 832 Lst<Atom> subs = new Lst<Atom>(); 833 834 Bond[] bonds = atoms[1].bonds; 835 boolean is23Double = false; 836 for (int pt = 0, i = Math.min(bonds.length, 4); --i >= 0;) { 837 Atom sub = bonds[i].getOtherAtom(atoms[1]); 838 if (sub == atoms[2]) { 839 is23Double = (bonds[i].order == Edge.BOND_COVALENT_DOUBLE); 840 continue; 841 } 842 subElements[0][pt] = sub.getElementSymbol(); 843 subVectors[0][pt] = V3.newVsub(sub, atoms[1]); 844 pt++; 845 } 846 subs.clear(); 847 bonds = atoms[2].bonds; 848 for (int pt = 0, i = Math.min(bonds.length, 4); --i >= 0;) { 849 Atom sub = bonds[i].getOtherAtom(atoms[2]); 850 if (sub == atoms[1]) 851 continue; 852 subElements[1][pt] = sub.getElementSymbol(); 853 subVectors[1][pt] = V3.newVsub(sub, atoms[2]); 854 pt++; 855 } 856 857 double theta = Measure.computeTorsion(atoms[0], atoms[1], atoms[2], 858 atoms[3], false); 859 double jvalue = Double.NaN; 860 if (is23Double || subElements[0][2] != null && subElements[1][2] != null) { 861 switch (mode) { 862 case MODE_CALC_3JHH: 863 jvalue = calc3JHHOnly(subElements, subVectors, v21, v34, v23, theta, 864 is23Double); 865 break; 866 case MODE_CALC_3JCH: 867 if (is23Double) 868 return null; 869 jvalue = calc3JCH(CHEquation, theta, is23Double); 870 break; 871 } 872 } else { 873 jvalue = calcJKarplus(theta); 874 } 875 return new double[] { theta, jvalue, atoms[1].i, atoms[2].i }; 876 } 877 getCalcType(Atom[] atoms, String[] elementsToFill, int mode)878 public static int getCalcType(Atom[] atoms, String[] elementsToFill, int mode) { 879 Atom atom1 = atoms[0]; 880 Atom atom4 = atoms[3]; 881 Bond[] bonds1 = atom1.bonds; 882 Bond[] bonds4 = atom4.bonds; 883 if (bonds1 == null || bonds4 == null || atom1.isCovalentlyBonded(atom4)) 884 return MODE_CALC_INVALID; 885 boolean allowNOE = ((mode & MODE_CALC_NOE) == MODE_CALC_NOE); 886 boolean allow3JHH = ((mode & MODE_CALC_3JHH) == MODE_CALC_3JHH); 887 boolean allow2JHH = ((mode & MODE_CALC_2JHH) == MODE_CALC_2JHH); 888 boolean allow3JCH = ((mode & MODE_CALC_3JCH) == MODE_CALC_3JCH); 889 boolean isGeminal = false; 890 Atom atom2 = atoms[1]; 891 Atom atom3 = atoms[2]; 892 if (atom2 == null) { 893 for (int i = 0; i < bonds1.length; i++) { 894 atom2 = bonds1[i].getOtherAtom(atom1); 895 if (atom2.isCovalentlyBonded(atom4)) { 896 isGeminal = true; 897 break; 898 } 899 for (int j = 0; j < bonds4.length; j++) { 900 atom3 = bonds4[j].getOtherAtom(atom4); 901 if (atom2.isCovalentlyBonded(atom3)) 902 break; 903 atom3 = null; 904 } 905 } 906 atoms[1] = atom2; 907 atoms[2] = atom3; 908 } else if (atom2.isCovalentlyBonded(atom4)) { 909 isGeminal = true; 910 } 911 String e1 = atom4.getElementSymbol(); 912 String e2 = (atom2 == null ? null : atom2.getElementSymbol()); 913 String e3 = (atom3 == null ? null : atom3.getElementSymbol()); 914 String e4 = atom1.getElementSymbol(); 915 916 boolean isHH = e1.equals("H") && e4.equals("H"); 917 if (isGeminal) { 918 mode = (allow2JHH && isHH && e2.equals("C") ? MODE_CALC_2JHH 919 : MODE_CALC_INVALID); 920 } else if (atom3 == null) { 921 mode = (allowNOE && isHH ? MODE_CALC_NOE : MODE_CALC_INVALID); 922 } else if (allow3JHH && isHH) { 923 mode = MODE_CALC_3JHH; 924 } else if (allow3JCH 925 && e2.equals("C") 926 && e3.equals("C") 927 && (e1.equals("H") && e4.equals("C") || e1.equals("C") 928 && e4.equals("H"))) { 929 mode = MODE_CALC_3JCH; 930 } else { 931 mode = MODE_CALC_INVALID; 932 } 933 if (mode != MODE_CALC_INVALID && elementsToFill != null) { 934 elementsToFill[0] = e1; 935 elementsToFill[1] = e2; 936 elementsToFill[2] = e3; 937 elementsToFill[3] = e4; 938 } 939 return mode; 940 } 941 calc2JHH(Atom h1, Atom c, Atom h2)942 private static double[] calc2JHH(Atom h1, Atom c, Atom h2) { 943 // I don't actually have an algorithm for this. Just getting a rough idea here. 944 double val = Double.NaN; 945 switch (c.getCovalentBondCount()) { 946 case 3: 947 // alkene 948 val = 1.5; 949 break; 950 case 4: 951 val = 12.0; 952 break; 953 default: 954 return null; 955 } 956 float angle = Measure.computeAngle(h1, c, h2, new V3(), new V3(), false); 957 return new double[] { angle, val, c.i }; 958 } 959 calcNOEImpl(Viewer viewer, Atom atom1, Atom atom2)960 private static double[] calcNOEImpl(Viewer viewer, Atom atom1, Atom atom2) { 961 try { 962 NMRNoeMatrix noeMatrix = (NMRNoeMatrix) viewer.ms.getInfo(atom1.mi, 963 "noeMatrix"); 964 double dist = 0, noe = Double.NaN; 965 if (noeMatrix == null) { 966 noeMatrix = NMRNoeMatrix.createMatrix(viewer, 967 viewer.getModelUndeletedAtomsBitSet(atom1.mi), 968 (String[]) viewer.ms.getInfo(atom1.mi, "noeLabels"), 969 (NOEParams) viewer.ms.getInfo(atom1.mi, "noeParams")); 970 } 971 dist = noeMatrix.getJmolDistance(atom1.i, atom2.i); 972 noe = noeMatrix.getJmolNoe(atom1.i, atom2.i); 973 return (Double.isNaN(noe) ? null : new double[] { dist, noe }); 974 } catch (Exception e) { 975 e.printStackTrace(); 976 return null; 977 } 978 } 979 980 @Override getNOEorJHH(Atom[] atoms, int mode)981 public double[] getNOEorJHH(Atom[] atoms, int mode) { 982 return calc2or3JorNOE(vwr, atoms, null, mode); 983 } 984 985 } 986