# -*- coding: utf-8 -*- """ Created on Fri Oct 1 18:13:35 2021 @author: Thomas Schultz ----------------------------------------------------------------------- Python Code for Molecular Structure Analysis from Rotational Constants ----------------------------------------------------------------------- The code is based on the following literature sources: Kraitchman 1953 [ J. Kraitchman, Am. J. Phys., 1953, 21, 17–24.] Costain 1958 [C. C. Costain, The J. Chem. Phys., 1958, 29, 864-874.] Laurie 1962 [ V. W. Laurie and D. R. Herschbach, J. Chem. Phys., 1962, 37, 1687–1693.] Typke 1978 [ V. Typke, J. Mol. Spectr., 1978, 69, 173–178.] Watson 1999 [ J. K. Watson, A. Roytburg and W. Ulrich, J. Mol. Spectr., 1999, 196, 102–119.] W. Gordy and R. Cook [Microwave Molecular Spectra, John Wiley & Sons, New York, 1984] The STRFIT program [http://www.ifpan.edu.pl/~kisiel/struct/struct.htm] see also: [Z. Kisiel, J. Mol. Spectr., 2003, 218, 58–67] <<< The code is work in progress, expect some rough edges. >>> To Do: Implement flag to fix fit parameters (turn them into constants, e.g., rH:1.1f, ca:0.012f). Save/load .ini file from UI to remember last program state. (use configparser?) Optional Z-Matrix to define molecular geometry and fit internals. (add menu options to import zmat file implement Zmat2Cart converter, might be slow) Implement more plotting options for fit results / Monte-Carlo results. Abbreviations COM: center of mass MOI: moment of inertia PMI: principal moments of inertia PAS: principal axis system (corresponding to PMI) r0,re,rs,rm: effective, equilibrium, Kraitchman, mass-weighted geometry parameters Required dependencies (tested version numbers): numpy (1.21.2) (math and linear algebra functions) scipy (1.7.3) (scientific functions, natural constants) matplotlib (3.5.0) (plotting) pyqtgraph (0.11.0) (tools for graphical user interface; includes PyQt5 (5.9.2)) To run the code (if unfamiliar with Python): Install Miniconda [https://docs.conda.io/en/latest/miniconda.html]. Open an 'Anaconda prompt' and install dependencies: >>> conda install numpy scipy matplotlib pyqtgraph ipython Then run 'ipython' and execute the script by copy/pasting code from below. The code should be run sequentially (top to bottom), e.g., in an iPython environment. The sections are organized as follows: --- SECTION (1) Imports, natural constants and conversion functions --- SECTION (2) Functions for geometry fitting --- SECTION (3) Graphical User Interface for fit --- [SECTION (4) Isotopologue data] """ #------------------------------------------------------------------------# #--- SECTION (1): Imports, natural constants and conversion functions ---# #------------------------------------------------------------------------# # Natural constants and imports #--------------------------------------------------------------- import sys import scipy.constants as sc # Natural onstants from scientific python import numpy as np # Linear algebra an math functions from numerical python npa = np.array # Shortcut for conversion of list to numpy array nan = np.nan # Shortcut for not-a-number from scipy.optimize import leastsq # Wrapper for MINPACK Levenberg Marquard minimization algorithm from scipy.spatial.transform import Rotation as Rot # Coordinate rotation functions from scipy.stats import pearsonr, spearmanr # Pearson / Spearman correlation analysis rng = np.random.default_rng() # Numpy default (PCG) pseudo-random-number-generator for Monte-Carlo simulations # rng.normal(value,sigma) gives normal-distributed sample with mean=value and std=sigma. #from tqdm import tqdm # Progressbar; wrap around iterable: > for tqdm(iterable): ... import pyqtgraph as pg # Plotting (based on PyQt) -- used to visualize Monte-Carlo results import pyqtgraph.opengl as gl # Required for plotting of molecules from PyQt5 import QtGui, QtWidgets, QtCore # Graphical user interface from PyQt5.Qt import QApplication # More Qt imports from functools import partial # Call functions with partial arguments import matplotlib.pyplot as plt # Plotting of fit results #%plt.style.use('seaborn-whitegrid') # Alternate plot style # Use inline plots in Jupyter or non-blocking graphs in iPython (choose your evironment). try: env = get_ipython().__class__.__name__ except: env = False if env == 'TerminalInteractiveShell': # Ipython environment mgc = get_ipython().magic mgc(u'%matplotlib qt') # Enable interactive, non-blocking plots if env == 'ZMQInteractiveShell': # Jupyter environment mgc = get_ipython().magic mgc(u'%pylab inline') # Enable inline plots # Default directory _dir = r'C:\Users\UNIST\OneDrive\benzene\2.Structure_analysis' DEBUG = False # Activate DEBUG to store diagnostic output (use in ipython environment) if DEBUG: # (Type 'log' to move data to clipboard) class logFile(): def __init__(self): import io # use log.write(string) to log data; log to copy to clipboard and reset. self.logfile = io.StringIO("") def __repr__(self): self.clip() self.logfile = io.StringIO("") return 'Copied log content to clipboard and reset log.' def write(self,txt): self.logfile.write(txt+'\n') def clip(self): QtWidgets.QApplication.clipboard().setText(self.logfile.getvalue().strip()) def emptyLog(self): self.logfile.truncate(0) def close(self): self.logfile.close() log = logFile() # Atomic masses from NIST. #--------------------------------------------------------------- # Coursey, J.S., Schwab, D.J., Tsai, J.J., and Dragoset, R.A. (2015), # Atomic Weights and Isotopic Compositions (version 4.1). [Online] # Available at: http://physics.nist.gov/Comp [accessed Oct. 16 2020]. # National Institute of Standards and Technology, Gaithersburg, MD. # For compactness, I substitute single-letter codes: D=H(2), T=H(3), # X=C(13), P=O(17), Q=O(18). NISTmass = {'H':1.00782503223, 'D':2.01410177812, 'T':3.0160492779, 'C':12.0, 'X':13.00335483507, 'O':15.99491461957,'P':16.99913175650, 'Q':17.99915961286} # Utility functions for unit conversion, etc #--------------------------------------------------------------- # Conversion between rotational frequency units (nu in Hz, MHz, cm-1) # and corresponding moments-of-inertia (I in amu*Angstrom**2) (DEFAULT UNIT) # Moments of inertia are: I = h/(8*pi**2*c*nu) #cm2I = lambda nu: sc.h/(8*sc.pi**2 *nu*sc.c*100) /(1e-20 * sc.atomic_mass) # wavenumbers to amu*A**2 cm2I = lambda nu: 16.857629191640175/nu if nu else np.nan # wavenumbers to amu*A**2 I2cm = lambda I: 16.857629191640175/I if I else np.nan # amu*A**2 to wavenumbers #MHz2I = lambda nu: sc.h/(8*sc.pi**2 *nu*1e6) /(1e-20 * sc.atomic_mass) # MHz to amu*A**2 MHz2I = lambda nu: 505379.00914143625/nu if nu else np.nan # MHz to amu*A**2 I2MHz = lambda nu: 505379.00914143625/nu if nu else np.nan # amu*A**2 to MHz I2I = lambda nu: nu # Null conversion ;) convertToUnit = {'amu*AA':I2I,'MHz':I2MHz,'1/cm':I2cm} # Conversion from amu*AA to other unit convertFromUnit = {'amu*AA':I2I,'MHz':MHz2I,'1/cm':cm2I} # Conversion from other unit to amu*AA #rndSD = lambda val,dig=9: 0 if val==0 else round(val, -int(np.floor(np.log10(abs(val))))+(dig-1)) # Round to sig. digits #rndList = lambda lst,dig=9: npa([rndSD(val,dig) for val in lst]) # (... for number tables; default:9 digits) gSum = lambda l: np.sqrt(np.sum(npa(l)**2)) # sqrt(sum(squares(list))) for error propagation #------------------------------------------------------# #--- SECTION (2) Functions for geometry fitting ---# #------------------------------------------------------# # (2.1) Coordinate system transformations. # (2.2) Calculation of principal moments of inertia. # (2.3) Molecule-specific internal -> Cartesian coordinate transformations # (2.4) Function to print/plot fit results. # (2.5) Fit function. # (2.6) Residual functions for r0,rs, and rm geometry parameter fit. # The fits are based on the following steps: # 1 Calculate reference molecule PMI_ref and diagnose PAS_ref rotation. # 2-7: For each isotopomer: # (2o) If required, rotate coordinates to PAS_ref coordinates. # 3 Calculate isotopomer PMI_i and diagnose rotation. # (4o) If required, rotate PMI_i to PAS_ref. # 5 Apply rovibrational correction terms. # (6o) If required, rotate back to PMI_i # 7 Compare calculated and experimental PMI_i values. # (2.1) Rotate or translate molecular coordinate system or inertial moment tensor #-------------------------------------------------------------------- # Rovibrational correction terms should be applied along the same axes for all isotopologues. # We rotate all coordinates into the principal axis system (PAS) of a reference molecule. # This requires: A function rotR to rotate the isotopomer coordinates to the reference PAS. # A function rotI to rotate isotopomer PMIs to the reference PAS. def rotI(I, eVec): # Rotate inertial moment tensor (3x3) or array (3x1) to principal axes. """ Rotate inertial moment tensor according to WAT99 eq. (35). The rotation is based on Eigenvectors eVec that map a inertial moment tensor I onto its principal moments PMI (via PMI,eVec = np.linalg.eigh(I)). This rotates the principal moments PMI back onto the tensor I or a tensor I to PMI (the tensors and rotation matrix are symmetric). If eVec corresponds to a different ('parent')isotopologue, this rotates the moments into the axis system of the parent. /I/ can contain principal moments [Ia,Ib,Ic] or a (3x3) I-tensor [[Iaa,Iab,...]. Function returns an I-tensor. Use this to rotate the PMI of an isotopologue to the reference coordinate system. """ if len(npa(I).shape) == 1: I = np.diag(I) # Expand 1D to 2D, if required I_rot = np.round(np.linalg.inv(eVec) @ I @ eVec,12) # Rotate and round to maintain 0 values return I_rot def rotR(rc,eInv): # Rotate Cartesian coordinates to principal axes. """ Rotate coordinates to a new axis system (typically to the principal axis system of a reference molecule). The rotation is based on inverse Eigenvectors eInv. eInv = np.linalg.inv(eVec) for Eigenvectors eVec, typically the Eigenvectors that describe the principal axis coordinate system rotation for a reference species. rotR(r,np.linalg.inv(eVec)) will rotate Cartesian isotopologue coordinates into the reference principal axis system. Use this if the internal2cart coordinate system does not correspond to the reference isotopologue principal axis system. STRfit definition of eInv from e: eInv = npa([[e[1,1]*e[0,2]-e[1,2]*e[0,1], e[1,2]*e[0,0]-e[1,0]*e[0,2], e[1,0]*e[0,1]-e[1,1]*e[0,0]], [e[2,2]*e[0,1]-e[2,1]*e[0,2], e[2,0]*e[0,2]-e[2,2]*e[0,0], e[2,1]*e[0,0]-e[2,0]*e[0,1]], [e[2,1]*e[1,2]-e[1,1]*e[2,2], e[2,2]*e[1,0]-e[2,0]*e[1,2], e[2,0]*e[1,1]-e[2,1]*e[1,0]]]) WW = e[2,0]*eInv[0,0] + e[1,0]*eInv[1,0] + e[0,0]*eInv[2,0] eInv = eInv / WW """ # eInv = np.linalg.inv(eVec) r_principal = [(eInv @ r).tolist() for r in rc] # dot product of eInv and r for each atom return npa(r_principal) def rotVA(rotMatrix): # Give rotation axis and angle from rotation matrix """ Return Euler angles for rotation matrix. """ if 1 in [np.count_nonzero(v) for v in rotMatrix]: # Rotation is axially symmetric rotMatrix = rotMatrix * np.linalg.eigh(rotMatrix)[0] # Fix sign before analyzing the rotation angle #print(f'symmetric rotation in {rotMatrix}') # (order/sign of eigenvectors is arbitrary) r = Rot.from_matrix(rotMatrix) #print(f"r_as_Euler (deg): {r.as_euler('zyx', degrees=True)}") return np.round(r.as_euler('zyx', degrees=True),1) def COM(rc, masses): # Translate molecular coordinates to their center-of-mass coordinats. """ Return center-of-mass coordinates for molecule with atom masses m_i at Cartesian coordinates r. """ _COM = np.round(np.dot(masses, rc ) / np.sum(masses),16) # center of mass coordinates (x,y,z) (COM) return npa([r -_COM for r in rc]) # (2.2) Principal moments of inertia (PMI) calculation #------------------------------------------------------- # Calculation of principal inertial moments for a molecule based on Cartesian coordinates # and atom masses. This is molecule independent. def PMI(rc, m_i, printResults=False): """ Calculation of molecular principal moments of inertia from atomic masses m_i and coordinates rc. Parameter 'atoms' must be a string describing isotope composition. For a description of the MOI calculation, see Gordy and Cook, Chapter 13, or Kraitchman1953 [Am. J. Phys. 21, 17 (1953); doi: 10.1119/1.1933338]. I refer to equations in Kraitchman1953 as KRA(equation number). """ _COM = np.round(np.dot(m_i, rc ) / np.sum(m_i),16) # Center of mass (COM [x,y,z]) x,y,z = rc [:,0]-_COM[0], rc [:,1]-_COM[1], rc [:,2]-_COM[2] # COM coordinates; KRA(5) I_xx = np.sum(m_i*(y**2 + z**2)) # Diag. elements for the inertial tensor I (KRA 2a) I_yy = np.sum(m_i*(z**2 + x**2)) # I_zz = np.sum(m_i*(x**2 + y**2)) # I_xy = - np.sum(m_i*x*y) # Off-diag. elements for I (Kra 2b) I_xz = - np.sum(m_i*x*z) # I_yz = - np.sum(m_i*y*z) # I = np.array([[I_xx, I_xy, I_xz], # I tensor (KRA 2f.) [I_xy, I_yy, I_yz], # Solving the secular equation (SE) gives Eigenvalues Ix,Iy,Iz, [I_xz, I_yz, I_zz]]) # which correspond to the principal moments of inertia (PMI) #print(f'I: {I}') _PMI, e_vectors = np.linalg.eigh(I) # -> PMI and eigenvectors from solving the SE if np.sum(abs(I - np.diag(np.diag(I)))) > 1e-12: # Non-diagonal I tensor: rc are not principal axis coordinates sym = 0 # -> symmetry flag 0 else: # Diagonal I tensor: rc are principal axis coordinates if (I_yy-I_zz)<1e-12 and (I_xx-I_yy)<1e-12: # Axis mapping abc == xyz sym = 2 # No axis remapping required (symmetry flag 2) else: sym = 1 # -> axis remapping required (symmetry flag 1) if printResults: print(f'Center-of-mass [A]: COM = np.{repr(np.round(COM,7))}') print(f'Inertial moment [amu*A**2]: PMI = np.{repr(np.round(_PMI,5))}') print(f'Rot. Constants [MHz]): A,B,C = {I2MHz(_PMI[0]):10.6f}, {I2MHz(_PMI[1]):10.6f}, {I2MHz(_PMI[2]):10.6f}') return _PMI, e_vectors, sym, I # Return PMI, eigenvectors, symmetry flag, I-tensor # (2.3) Molecule-specific inerternal -> Cartesian coordinate transformations #-------------------------------------------------------------------- # This maps a minimal set of internal coordinates to Cartesian atomic coordintes, # accounting for molecular symmetry. This is different for every molecule. # The optional Laurie correction assumes that a large H/D isotope effect is localized and # can be described by one mass-weighted parameter dH. Typical values for the Laurie # correction parameter are close to dH = 10 sqrt(amu)*mA (corresponds to rH-rD ~= 3mA). # For benzene and H2O, I also tested a int2cart function that fits independent X-H and X-D # bond lengths. These functions use the 'Laurie parameter' to describe the X-D bond length instead. # (2.3a) Benzene, C6H6 (extra function for rC,rH,rD 3-parameter fit) def bz_internal2cart(r, m_i=None, Laurie=False): """ Benzene Cartesian atom coordinates from /r/ = [rH,rC,(dH)] internal coordinates. rH, rC are atom distances to center-of-mass. Optional parameter dH must be specified only when Laurie correction is /Laurie/=True. /m_i/ should contain a list of atom masses but is only required for the Laurie correction. """ if m_i == None: m_i = [] rC,rH = r[:2] # first parameters are rC, rH cos30 = np.cos(30/180*np.pi) # cos(30 deg) sin30 = 1/2 # sin(30 deg) if Laurie: # different bond length rH,rD,rT dH = r[-1] # optional last parameter is dH M = sum(m_i) # Isotopologue mass rHD = lambda atom: rH + dH* np.sqrt( M / (m_i[atom]*(M-m_i[atom])) ) # Laurie correction else: rHD = lambda atom: rH # No Laurie correction requested C0 = [ rC, 0, 0] # x,y,z coordinates for C1 C1 = [ rC*sin30, rC*cos30, 0] # x,y,z coordinates for C2 C2 = [ -rC*sin30, rC*cos30, 0] # ... C3 = [ -rC, 0, 0] C4 = [ -rC*sin30, -rC*cos30, 0] C5 = [ rC*sin30, -rC*cos30, 0] H6 = [ rHD(6), 0, 0] # x,y,z coordinates for H1 H7 = [ rHD(7)*sin30, rHD(7)*cos30, 0] # x,y,z coordinates for H2 H8 = [ -rHD(8)*sin30, rHD(8)*cos30, 0] # ... H9 = [ -rHD(9), 0, 0] H10 = [-rHD(10)*sin30, -rHD(10)*cos30, 0] H11 = [ rHD(11)*sin30, -rHD(11)*cos30, 0] rc = npa([C0,C1,C2,C3,C4,C5,H6,H7,H8,H9,H10,H11]) return rc # Return Cartesian coordinates rc def bz_internal2cart_rC_rH_rD_Laurie(r, m_i=None, Laurie=True): """ Benzene Cartesian atom coordinates from /r/ = [rH,rC,dH] internal coordinates. The Laurie dH parameter is used to fit a r(C-D) bond length independent from r(C-H). rH, rC, dH are atom distances of H,C,D atoms to the center-of-mass. Note that we must use the "Laurie Correction" version of the fit to ensure re-calculation of coordinates for each isotopologue, hence the dH parameter represents rD: We replaced the mass-dependent Laurie correction with separate H/D bond lengths. /m_i/ must list atom masses for all 12 atoms in clockwise direction and is required for the H/D correction, to determine if an atom is H or D. Note that this approach gives effective r0 bond lengths for rH, rD, whereas the Laurie correction gives approximate re bond lengths (plus the correction term dH to get to r0). """ if m_i == None: m_i = [] cos30 = np.cos(30/180*np.pi) # cos(30 deg) sin30 = 1/2 # sin(30 deg) rC,rH = r[:2] # first parameters are rC, rH rD = r[-1] # last parameter is rD rHD = lambda atom: rH if int(m_i[atom])==1 else rD # Distinguish rH and rD C0 = [ rC, 0, 0] # x,y,z coordinates for C1 C1 = [ rC*sin30, rC*cos30, 0] # x,y,z coordinates for C2 C2 = [ -rC*sin30, rC*cos30, 0] # ... C3 = [ -rC, 0, 0] C4 = [ -rC*sin30, -rC*cos30, 0] C5 = [ rC*sin30, -rC*cos30, 0] H6 = [ rHD(6), 0, 0] # x,y,z coordinates for H1 H7 = [ rHD(7)*sin30, rHD(7)*cos30, 0] # x,y,z coordinates for H2 H8 = [ -rHD(8)*sin30, rHD(8)*cos30, 0] # ... H9 = [ -rHD(9), 0, 0] H10 = [-rHD(10)*sin30, -rHD(10)*cos30, 0] H11 = [ rHD(11)*sin30, -rHD(11)*cos30, 0] rc = npa([C0,C1,C2,C3,C4,C5,H6,H7,H8,H9,H10,H11]) return rc # Return Cartesian coordinates rc def bz_internal2cart_rC_rH_r13C_Laurie(r, m_i=None, Laurie=True): """ Benzene Cartesian atom coordinates from /r/ = [rH,rC,dH] internal coordinates. The Laurie dH parameter is used to fit a r(C-D) bond length independent from r(C-H). rH, rC, dH are atom distances of H,C,D atoms to the center-of-mass. Note that we must use the "Laurie Correction" version of the fit to ensure re-calculation of coordinates for each isotopologue, hence the dH parameter represents rD: We replaced the mass-dependent Laurie correction with separate H/D bond lengths. /m_i/ must list atom masses for all 12 atoms in clockwise direction and is required for the H/D correction, to determine if an atom is H or D. Note that this approach gives effective r0 bond lengths for rH, rD, whereas the Laurie correction gives approximate re bond lengths (plus the correction term dH to get to r0). """ if m_i == None: m_i = [] cos30 = np.cos(30/180*np.pi) # cos(30 deg) sin30 = 1/2 # sin(30 deg) r12C,rH = r[:2] # first parameters are rC, rH rHD = lambda atom: rH r13C = r[-1] # last parameter is rD rC = lambda atom: r12C if int(m_i[atom])==12 else r13C # Distinguish rH and rD C0 = [ rC(0), 0, 0] # x,y,z coordinates for C1 C1 = [ rC(1)*sin30, rC(1)*cos30, 0] # x,y,z coordinates for C2 C2 = [ -rC(2)*sin30, rC(2)*cos30, 0] # ... C3 = [ -rC(3), 0, 0] C4 = [ -rC(4)*sin30, -rC(4)*cos30, 0] C5 = [ rC(5)*sin30, -rC(5)*cos30, 0] H6 = [ rHD(6), 0, 0] # x,y,z coordinates for H1 H7 = [ rHD(7)*sin30, rHD(7)*cos30, 0] # x,y,z coordinates for H2 H8 = [ -rHD(8)*sin30, rHD(8)*cos30, 0] # ... H9 = [ -rHD(9), 0, 0] H10 = [-rHD(10)*sin30, -rHD(10)*cos30, 0] H11 = [ rHD(11)*sin30, -rHD(11)*cos30, 0] rc = npa([C0,C1,C2,C3,C4,C5,H6,H7,H8,H9,H10,H11]) return rc # Return Cartesian coordinates rc # (2.3b) Water, H2O (symmetric and asymmetric coord. definitions for testing) def h2o_internal2cart_sym(r, m_i=None, Laurie=False): """ Water Cartesian atom coordinates from /r/ = [rOH,angle_HOH,(dH)] internal coordinates. Optional parameter dH must be specified only when Laurie correction is True. /m_i/ lists atom masses for the three atoms and is required for the Laurie correction. """ if m_i == None: m_i = [] rOH,aHOH = r[:2] # First 2 parameters are r(OH), angle(HOH) if Laurie: # different bond length rH,rD,rT dH = r[-1] # optional last parameter is dH M = sum(m_i) # Isotopologue mass rHD = lambda atom: rOH + dH* np.sqrt( M / (m_i[atom]*(M-m_i[atom])) ) # Laurie correction else: rHD = lambda atom: rOH # No Laurie correction requested sin_a = np.sin(aHOH*np.pi/360) # sin(alpha/2) for HOH bond angle (alpha/2 in radians) cos_a = np.cos(aHOH*np.pi/360) # cos(alpha/2) for HOH bond angle (alpha/2 in radians) H1 = [-rHD(0)*sin_a, rHD(0)*cos_a, 0] # x,y,z coordinates for hydrogen (atom 0) O2 = [ 0, 0, 0] # x,y,z coordinates for oxygen (atom 2) H3 = [ rHD(2)*sin_a, rHD(2)*cos_a, 0] # x,y,z coordinates for hydrogen (atom 1) #print(f'r0:{_r} + Laurie:{rHD(0)}') # Check bond angle + Laurie correction return np.array([H1,O2,H3]) # Return Cartesian coordinates rc def h2o_internal2cart_rHrD_sym(r, m_i=None, Laurie=True): """ Water Cartesian atom coordinates from /r/ = [rOH,angle_HOH,dH] internal coordinates. The Laurie parameter dH is used to fit a r(O-D) bond length independent from r(O-H). Note that we must use the "Laurie Correction" version of the fit to ensure re-calculation of coordinates for each isotopologue, hence the dH parameter represents r(O-D). /m_i/ lists atom masses for the three atoms and is required to distinguish H/D atoms. """ if m_i == None: m_i = [] rOH,aHOH = r[:2] # First 2 parameters are r(OH), angle(HOH) rOD = r[-1] if Laurie: # different bond length rH,rD,rT rHD = lambda atom: rOH if int(m_i[atom])==1 else rOD # Distinguish rOH and rOD else: rHD = lambda atom: rOH # No Laurie correction requested sin_a = np.sin(aHOH*np.pi/360) # sin(alpha/2) for HOH bond angle (alpha/2 in radians) cos_a = np.cos(aHOH*np.pi/360) # cos(alpha/2) for HOH bond angle (alpha/2 in radians) H1 = [-rHD(0)*sin_a, rHD(0)*cos_a, 0] # x,y,z coordinates for hydrogen (atom 0) O2 = [ 0, 0, 0] # x,y,z coordinates for oxygen (atom 2) H3 = [ rHD(2)*sin_a, rHD(2)*cos_a, 0] # x,y,z coordinates for hydrogen (atom 1) #print(f'r0:{_r} + Laurie:{rHD(0)}') # Check bond angle + Laurie correction return np.array([H1,O2,H3]) # Return Cartesian coordinates rc def h2o_internal2cart(r, m_i=None, Laurie=False): """ Water Cartesian atom coordinates from /r/ = [rOH,angle_HOH,(dH)] internal coordinates. Optional parameter dH must be specified only when Laurie correction is True. /atoms/ lists atom masses for the three atoms and is required for the Laurie correction. """ if m_i == None: m_i = [] rOH,aHOH = r[:2] # First 2 parameters are r(OH), angle(HOH) if Laurie: # different bond length rH,rD,rT dH = r[-1] # optional last parameter is dH M = sum(m_i) # Isotopologue mass rHD = lambda atom: rOH + dH* np.sqrt( M / (m_i[atom]*(M-m_i[atom])) ) # Laurie correction else: rHD = lambda atom: rOH # No Laurie correction requested _a = (aHOH-90)*np.pi/180 # Convenient conversion to degree angle cos_a = np.cos(aHOH*np.pi/360) # cos(alpha/2) for HOH bond angle (alpha/2 in radians) H1 = [ -rHD(0), 0, 0] # x,y,z coordinates for hydrogen (atom 0) O2 = [ 0, 0, 0] # x,y,z coordinates for oxygen (atom 2) H3 = [ rHD(2)*(np.sin(_a)), rHD(2)*np.cos(_a), 0] # x,y,z coordinates for hydrogen (atom 1) #print(f'r0:{_r} + Laurie:{rHD(0)}') # Check bond angle + Laurie correction return np.array([H1,O2,H3]) # Return Cartesian coordinates rc # (2.3c) Ozone, O3 (symmetric and asymmetric coord. definitions for testing) def O3_int2cart_sym(r, m_i=None, Laurie=False): """ Ozone Cartesian atom coordinates from /r/ = [rOO,angle_OOO] internal coordinates. Symmetric mapping to minimize coordinate rotations: O1 Symmetric structure with r(OO) and angle(OOO) / \ O2 O3 """ if m_i == None: m_i = [] _r,a = r[:2] _a = a*np.pi/180/2 # alpha/2 in radians O1 = [-_r*np.sin(_a), _r*np.cos(_a), 0] # x,y,z coordinates for O1 O2 = [ 0, 0, 0] # x,y,z coordinates for O2 O3 = [ _r*np.sin(_a), _r*np.cos(_a), 0] # ... return np.array([O1,O2,O3]) # Return Cartesian coordinates rc def O3_int2cart(r, m_i=None, Laurie=False): """ Ozone Cartesian atom coordinates from /r/ = [rOO,angle_OOO] internal coordinates. Asymmetric mapping as in the STRFIT example: ----> x O1--O2 Structure as defined in STRFIT, | \ with r(OO) and angle(OOO) V O3 y """ if m_i == None: m_i = [] _r,a = r[:2] _a = (a-90)*np.pi/180 # Conversion degrees to radians O1 = [ -_r, 0, 0] # x,y,z coordinates for O1 O2 = [ 0, 0, 0] # x,y,z coordinates for O2 O3 = [ _r*(np.sin(_a)), _r*np.cos(_a), 0] # ... return np.array([O1,O2,O3]) # Return Cartesian coordinates rc # (2.3d) HCN, linear molecule def HCN_int2cart(r, m_i=None, Laurie=False): """ HCN Cartesian atom coordinates from /r/ = [rHC,rCN] internal coordinates. """ if m_i == None: m_i = [] rHC,rCN = r[:2] if Laurie: # different bond length rH,rD,rT dH = r[-1] # optional last parameter is dH M = sum(m_i) # Isotopologue mass rHD = lambda atom: rHC + dH* np.sqrt( M / (m_i[atom]*(M-m_i[atom])) ) # Laurie correction else: rHD = lambda atom: rHC # No Laurie correction requested H1 = [ -rHD(0), 0, 0] # x,y,z coordinates for H C2 = [ 0, 0, 0] # x,y,z coordinates for OC N3 = [ rCN, 0, 0] # ... return np.array([H1,C2,N3]) # Return Cartesian coordinates rc # (2.3e) HNCO, nonlinear hydride def HNCO_int2cart(r, m_i=None, Laurie=False): """ HCN Cartesian atom coordinates from /r/ = [rHC,rCN] internal coordinates. """ if m_i == None: m_i = [] rHN,rNC,rCO,aHNC,aNCO = r[:5] if Laurie: # different bond length rH,rD,rT dH = r[-1] # optional last parameter is dH M = sum(m_i) # Isotopologue mass _rHN = lambda atom: rHN + dH* np.sqrt( M / (m_i[atom]*(M-m_i[atom])) ) # Laurie correction else: _rHN = lambda atom: rHN # No Laurie correction requested a1 = aHNC/180*np.pi a2 = aNCO/180*np.pi H1 = [ _rHN(0)*np.cos(a1), _rHN(0)*np.sin(a1), 0] # x,y,z coordinates for H N2 = [ 0, 0, 0] # x,y,z coordinates for OC C3 = [ rNC, 0, 0] # ... O4 = [rNC - rCO*np.cos(a2), -rCO*np.sin(a2), 0] # ... return np.array([H1,N2,C3,O4]) # Return Cartesian coordinates rc # (2.3f) furan # I give several int2cart functions to illustrate partial # fitting of parameters and the flexibility of fit parameters. # The C2h molecule is fully described by 8 parameters. def furan_int2cart(r, m_i=None, Laurie=False): """ Furan Cartesian atom coordinates from planar internal coordinates: H5 H6 ^ (y) /r/ = [rCO,rCC, C1 = C2 | rCH1,rCH2, / | | aCOC,aOCC, O0 | ------> (x) aOCH,aCCH] \ | C3 = C4 H7 H8 Notes: All coordinate assignments clockwise from O1. Angles are inner angles for lowest numbered atoms. Reduced coordinates based on C2h symmetry: 3 bond lengths + 1 angle define carbon frame or 2 bond length + 2 angles define carbon frame (used here). r = [1.36,1.35,1.1,1.1,106,110,120,120] p = plotMolecule(['OCCCCHHHH',furan_int2cart(r)]) """ if m_i == None: m_i = [] cos,sin,pi = np.cos, np.sin, np.pi rOC,rCC, rCH1,rCH2, *rest = r # Bond lengths from r aCOC, aOCC, aOCH,aCCH = [a/180*np.pi for a in rest[:4]] # Bond angles from r, degree -> radians a1 = aCOC/2 # Angle O1C2 versus x-axis a2 = a1 -pi +aOCC # Angle C2C3 versus x-axis a3 = a1 + pi-aOCH # Angle O1C2H6 versus x-axis a4 = a2 + pi-aCCH # Angle O1C2H6 versus x-axis if Laurie: # Laurie corrected bond length for rH,rD,rT dH = r[-1] # optional last parameter is dH M = sum(m_i) # Isotopologue masses _rCH5 = rCH1 + dH* np.sqrt( M / (m_i[5]*(M-m_i[atom])) ) # Laurie corrected rCH1 _rCH6 = rCH2 + dH* np.sqrt( M / (m_i[6]*(M-m_i[atom])) ) # Laurie corrected rCH2 _rCH7 = rCH1 + dH* np.sqrt( M / (m_i[7]*(M-m_i[atom])) ) # Laurie corrected rCH1 _rCH8 = rCH2 + dH* np.sqrt( M / (m_i[8]*(M-m_i[atom])) ) # Laurie corrected rCH2 else: # No Laurie correction requested _rCH5 = rCH1 _rCH6 = rCH2 _rCH7 = rCH1 _rCH8 = rCH2 O0 = npa([ 0, 0, 0]) # x,y,z coordinates for O1 C1 = O0 + npa([ rOC*cos(a1), rOC*sin(a1), 0]) # x,y,z coordinates for C2 C2 = C1 + npa([ rCC*cos(a2), rCC*sin(a2), 0]) # ... C3 = C1 * npa([1,-1,0]) # Remaining atoms from symmetry C4 = C2 * npa([1,-1,0]) # ... H5 = C1 + npa([ _rCH5*cos(a3), _rCH5*sin(a3), 0]) # ... H6 = C2 + npa([ _rCH6*cos(a4), _rCH6*sin(a4), 0]) # ... H7 = C3 + npa([ _rCH7*cos(a3), -_rCH7*sin(a3), 0]) # ... H8 = C4 + npa([ _rCH8*cos(a4), -_rCH8*sin(a4), 0]) # ... return np.array([O0,C1,C2,C3,C4,H5,H6,H7,H8]) # Return Cartesian coordinates rc def furan_int2cart_noH(r, m_i=None, Laurie=False): """ Wrapper for furan_int2cart with fixed C-H bond lengths and angles. E.g., set rCH1,rCH2, aOCH,aCCH to ab initio values. r = [1.36,1.35, 106,110] # fit parameters [rCO,rCC1,rCC2, aCOC,aOCC] plotMolecule(['OCCCCHHHH',furan_int2cart_noH(r)]) """ rCO,rCC1, aCOC,aOCC = r[:4] rCH1,rCH2, aOCH,aCCH = 1.077343,1.077343, 116.64,126.56 # Values from PBEh-3c calculation in ORCA return furan_int2cart([rCO,rCC1,rCH1,rCH2, aCOC,aOCC,aOCH,aCCH], m_i=m_i, Laurie=False) def furan_int2cart_scaledH(r, m_i=None, Laurie=False): """ Wrapper for furan_int2cart with scaled C-H bond lengths and fixed CH angles. Set rCH1,rCH2, aOCH,aCCH to ab initio values and scale r values with single scaling factor r = [1.36,1.36,1.43, 106,110, 1] # fit parameters [rCO,rCC1,rCC2, aCOC,aOCC, scaleH] plotMolecule(['OCCCCHHHH',furan_int2cart_scaledH(r)]) """ rCO,rCC1, aCOC,aOCC, scaledH = r[:5] rCH1,rCH2, aOCH,aCCH = 1.077343,1.077343, 116.64,126.56 # Values from PBEh-3c calculation in ORCA return furan_int2cart([rCO,rCC1,rCH1*scaledH,rCH2*scaledH, aCOC,aOCC,aOCH,aCCH], m_i=m_i, Laurie=False) def furan_int2cart_xy(r, m_i=None, Laurie=False): """ Furan Cartesian atom coordinates from minimal set of Cartesian coordinates: We can freely choose our 8 parameters as long as the parameters give a complete decription of the molecule. H1 H2 ^ (y) /r/ = [xC1,yc1, C1 = C2 | xC2,yC2, / | | xH1,yH1, O | ------> (x) xH2,yH2] \ | C1'= C2' H1' H2' Notes: Oxygen is at the origin [0,0,0]. Reduced coordinates based on C2h symmetry: We only require y,x coordinates for symmetry-unique atoms. r = [0.5,1, 2.,0.7, 0.4,2, 2.1,2] p = plotMolecule(['OCCCCHHHH',furan_int2cart_xy(r)]) """ if m_i == None: m_i = [] xC1,yc1, xC2,yC2, xH1,yH1, xH2,yH2 = r[:8] # Map onto explicit parameter names O1 = npa([ 0, 0, 0]) # x,y,z coordinates for O1 C2 = O1 + npa([ xC1, yC1, 0]) # x,y,z coordinates for C2 C3 = C2 + npa([ xC2, yC2, 0]) # ... H6 = C2 + npa([ xH1, yH1, 0]) # ... H7 = C3 + npa([ xH2, yH2, 0]) # ... C4 = C3 * npa([ 1, -1, 0]) # Remaining atoms from symmetry C5 = C2 * npa([ 1, -1, 0]) # ... H8 = H7 * npa([ 1, -1, 0]) # ... H9 = H6 * npa([ 1, -1, 0]) # ... return np.array([O1,C2,C3,C4,C5,H6,H7,H8,H9]) # Return Cartesian coordinates rc # (2.3g) Analysis of int2cart symmetry # For symmetric molecules, we should use a reduced number of rovibrational correction # parameters. The following function analyzes the molecular symmetry for a molecule # defined by an int2cart function (above), based on the molecular PMI. def analyzeSymmetry(isoList,isoMass, int2cart, fitParameters, nPars): """ Analyze molecular symmetry for all isotopologues through PMI. Returns isoList with appended symmetry and isotope mass information: isoList += symmetryString, [z,y,x] Euler rotation angles vs. parent, [a,b,c] <-> [x,y,z] axis mapping vs. parent, [atomMasses] """ retList=[] r = list(fitParameters.values()) # Fit parameters (starting with internal geometry pars.) # Analyze parent molecule (first in list) iso = isoList[nPars['parent']] # Reference molecule (analyze PMI rotation) m_ref = [isoMass[atom] for atom in iso[0]] # Atom masses for reference molecule rc0 = int2cart(r,m_ref,nPars['Laurie']) # Convert internal to Cartesian coordinates I0,eVec_ref,sym_ref,I = PMI(rc0,m_ref) # Evaluate reference PMI if not sym_ref==2: # Reference coordinate system is not principal axis system (PAS)? eInv_ref = np.dot(np.linalg.inv(eVec_ref),nPars['rotMat']) # -> Rotation matrix into reference PAS + user supplied angle rotEuler = rotVA(eInv_ref) # Euler zyx rotation vector/angle into reference PAS print(f'Reference isotopologue not symmetric, rotated by {rotEuler} degrees around axes [z,y,x]') else: eInv_ref = nPars['rotMat'] # -> Rotation matrix into reference PAS + user supplied angle rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS: a,b,c = x,y,z. for iso in isoList: # Now analyze all isotopologues isoSym = 'asymmetric' # Assumption of asymmetric top m = [isoMass[atom] for atom in iso[0]] # Create mass list for atoms in isotopologue [Ia,Ib,Ic],eVec,sym,I = PMI(rc,m) # Evaluate PMI #--- Analyze principal axes coordinate system for isotopologue rotAxes = [0,0,0] # Rotation between PMI and reference PAS coordinate system axisMap = [0,1,2] # Axis mapping between PMI and reference PAS coordinate system if sym == 0: # PMI coordinate system is rotated eInv = np.linalg.inv(eVec) # Rotation matrix into reference PAS rotAxes = rotVA(eInv) # --> [z,y,x] Euler rotation angles into reference PAS if sym == 1: # abc <-> xyz axis mapping differs between PMI and PAS axisMap = list(np.argsort( np.linalg.eig(I)[0] )) # Mapping of axes [a,b,c] to [x,y,z] #--- Analyze molecular symmetry for isotopologue if abs(Ia - Ib) < 1e-12: # Ia == Ib isoSym = 'oblate' # --> Oblate top if abs(Ia - Ic) < 1e-12: # Ia == Ib == Ic isoSym = 'spherical' # --> Spherical top if abs(Ia+Ib-Ic) < 1e-12: # Ia + Ib == Ic (planar) if isoSym == 'oblate': # also oblate? isoSym = 'oblate pl.' # --> Planar oblate top (Cnv with n>=3 ?) else: isoSym = 'planar' # Only planar molecule if Ia == 0: isoSym = 'linear' # --> Linear molecule (Ia==0, Ib==Ic) elif (Ib - Ic < 1e-12): # Ib == Ic isoSym = 'prolate' # --> Prolate top retList.append(iso[:4] + [isoSym,rotAxes,axisMap] + [m]) # Append symmetry, mass info to list return retList # (2.4) Function to print and plot fit results #-------------------------------------------------------------------- # Plotting the fit results helps to spot outliers in the experimental data. def printFitResults(fitList, fitParameters, Fit, title='Fit results'): """ Print fitted parameters and uncertainties and plot residuals. Beware: in weighted fits, outliers can greatly skew the Chi-squared result.""" #--- Unpack variables, etc. ---# pNames = list(fitParameters.keys()) # Unpack parameter names nPars = len (pNames) # Number of parameters (popt, pcov, infodict, errmsg, ier) = Fit # Unpack fit results residuals = infodict['fvec'] # Residuals for uncertainty calculation ... weights = [i for iso in fitList for i in iso[2] if i>0] scaledResiduals = residuals/weights # Remove weight scaling of residuals scaledChiSqr = (scaledResiduals**2).sum() # Chi squared for unweighted residuals chi_sq = (residuals**2).sum() # Chi squared (weigthed) DOF = len(residuals)-len(popt) # Fit degrees of freedom s_sq = chi_sq / DOF # ... reduced chi**2 try: rcov = pcov * s_sq # ... covariances (cf. scipy.optimize.curve_fit) pSigma = np.sqrt(np.diag(rcov)) # ... estimated 1-sigma parameter uncertainties except: print(f'failed to calculate covariances') rcov = np.full((nPars,nPars),np.NAN) pSigma = [np.NAN]*len(popt) # ... possibly insufficient data for covariances #--- Print fit results ---# pCostain = [np.nan if n[0]!='r' else np.sqrt((0.0015/r)**2 + u**2) for u,r,n in zip(pSigma,popt,pNames)] # Costain errors only for bond lengths resultString = f'Fit Function: {title}\n' resultString += 'Fitted parameters (1-sigma uncertainties) [Costain uncertainties]:\n' resultString += f' ------------------------------------\n' resultString += '\n'.join([f'{n:>10} = {p:>9.5f}({s:3.2e}) [{sC:.5f}]' for n,p,s,sC in zip(pNames,popt,pSigma,pCostain)]) # I shouold really analyze the parameter names and then decide on useful output. Here are some benzene-specific hacks. if 'r' in pNames[0] and 'r' in pNames[1]: # Quick and dirty: give r(atom1-atom2) resultString += f'\n ------------------------------------\n' #resultString += f' r12 = {popt[1]-popt[0]:>9.5f}({np.sqrt(pSigma[1]**2+pSigma[0]**2):.5f}) [{np.sqrt(pCostain[1]**2 + pCostain[0]**2):.5f}]\n' if 'rC' in pNames[0] and 'rH' in pNames[1] and 'dH' in pNames[-1]: # Quick and dirty: give Laurie-corrected CH, CD bond lengths in benzene rC,rH,dH = popt[0],popt[1],popt[-1] # Take first molecule in fitlist, calculate r0CH,r0CD, srC,srH,sdH = pSigma[0],pSigma[1],pSigma[-1] # for the corresponding singly-deuterated molecule. molMass = sum(fitList[0][-1]) r0CH = (rH + dH * np.sqrt( (molMass+1) / (molMass))) -rC # r_0(H) from rH and Laurie parameter dH sr0CH = (srH**2 + (sdH*np.sqrt( (molMass+1) / (molMass)))**2 + srC**2)**0.5 r0CD = rH + dH * np.sqrt( (molMass+1) / (2*(molMass-1))) -rC # r_0(D) from rm and Laurie parameter dH sr0CD = (srH**2 + (sdH*np.sqrt( (molMass+1) / (2*(molMass-1))))**2 + srC**2)**0.5 resultString += f' r0(H) = {r0CH:>9.5f}({sr0CH:.5f})\n' resultString += f' r0(D) = {r0CD:>9.5f}({sr0CD:.5f})\n' resultString += f'\n----------------------------------------------------------------------------\n' resultString += f'{" Chi-squared":>18} = {scaledChiSqr:9.6f} [Sum(Iobs-Icalc)**2]\n' #resultString += f'{" weighted Chi-sq":>18} = { chi_sq/np.mean(weights):9.6f} [Sum(Iobs-Icalc)**2]\n' resultString += f'{"Deviation of fit":>18} = {np.sqrt(scaledChiSqr):9.6f} uA^2 [Sqrt(scaledChiSqr/Ndegf), Ndegf= {DOF}]' resultString += f'\n----------------------------------------------------------------------------\n' resultString += f'\nParameter Correlation Coefficients\n{"":6}' resultString += ''.join([f'{pNames[i]:>8}' for i in range(len(fitParameters))]) for i in range(len(fitParameters)): # Print parameter correlation values resultString += f'\n{pNames[i]:>8}' for j in range(i+1): resultString += f'{rcov[i,j]/pSigma[i]/pSigma[j]:>8.3f}' #--- Plot residuals ---# plt.rcParams["figure.figsize"] = 12,4 #labels, sigmaI, weights = [],[],[] #for iso in fitList: # Iterate over all isotopologues # weight = npa(iso[2]) # Fit weights (prop. to 1/sigma uncertainty) # weights += list(npa(weight[weight>0])) # List of fit weights # sigmaI += list(npa(iso[1][1])[weight>0]) # 1-sigma experimental uncertainty (plot on residuals) # labels += [f'{iso[3]}\n{iso[0]}, {axis}' for axis in npa(['A','B','C'])[weight>0]] sigmaI = [s for iso in fitList for s,w in zip(iso[1][1],iso[2]) if w>0] if len(fitList)<15: labels = [f'{iso[3]}\n{iso[0]}, {axis}' for iso in fitList for axis in npa(['A','B','C'])[npa(iso[2])>0]] else: labels = [f'{iso[3]} {iso[0]}, {axis}' for iso in fitList for axis in npa(['A','B','C'])[npa(iso[2])>0]] Res = residuals / weights # Actual (unweighted) residuals for PMIs legend = f'{title}: ' + ' '.join([f'{n} = {p:.4f}({s:.4f})' for n,p,s in zip(pNames,popt,pSigma)][:2]) plt.errorbar(np.arange(len(Res)), Res, yerr=sigmaI, capsize=3,fmt='o',ecolor='C0',label=f"{legend}") plt.xticks(np.arange(len(labels)), labels, rotation=85) plt.ylabel('Fit residuals [amu*A**2]', fontsize = 11) plt.get_current_fig_manager().canvas.manager.set_window_title(title) plt.gca().yaxis.set_major_formatter(plt.FormatStrFormatter('%.0e')) # plt.title('\n'.join(fitResults), fontsize = 11, y=1.0, x=0.98, pad=-40, loc='right') plt.axhline(y = 0, color = 'r', linestyle = '-') plt.legend(loc='lower right') plt.tight_layout() plt.draw() return resultString # (2.5) Fit Function #------------------------------------------------------- # I implement a least-squares molecular structure anaysis methods as described in WAT99 # with the Levenberg Marquardt nonlinear fit method implemented via scipy.linalg.leastsq # (wrapper for MINPACK’s lmdif and lmder algorithms). # The residuals function determines the type of rovibrational corrections. def fit(fitList, residualFunc, int2cart=bz_internal2cart, fitParameters={}, nPars={'Laurie':0}): """ A scipy.optimize.leastsq fit using the residual function specified in Residuals. Parameters must hold a dictionary with all fit parameters (internals and roVib. correction terms). Parameter dH activates Laurie correction. Todo: Implement option to fit only selected parameters (fixedParameters as args.). """ if DEBUG: log.write(f'New Fit: {repr(residualFunc).split()[1]}{"L" if nPars["Laurie"]==1 else ""}, {repr(int2cart).split()[1]}') startPar = list(fitParameters.values()) # List of parameter start values Fit = leastsq(residualFunc, startPar, args=(fitList,nPars,int2cart), full_output=1) # Fit fitTitle = f'{repr(residualFunc).split()[1]}{"L" if nPars["Laurie"]==1 else ""}, {repr(int2cart).split()[1]}' resultString = printFitResults(fitList, fitParameters, Fit, title=fitTitle) return resultString,Fit # (2.6) Plotting of molecule and printing of molecular geometry #------------------------------------------------------- # To check int2cart conversion functions and guess parameters, # I implement a simplistic molecule viewer to plot a ball and stick model. def findBonds(atoms, maxDist=2): """ Returns list of atoms with distance < maxDist, i.e., chemical bonds. (atoms is numpy list of strings of atom label and cartesian coordinates, e.g. array([['C', '0.782463', '-1.250359', '-0.076314'],[...]]). """ bondList = [] atomLabels = atoms[:,0] coords = atoms[:,1:].astype(float) nAtoms = len(coords) distances = np.zeros([nAtoms, nAtoms]) for i in range(nAtoms): for j in range (i+1,nAtoms): distance = np.linalg.norm(coords[i] - coords[j]) if distance < maxDist: bondList.append(((i,j),distance)) # bondcenters.append( (coords[i] + coords[j]) /2) # Todo: write bondlengths on bond return bondList def bondList(r, maxDist=1.8): """ Returns list of atoms with distance < maxDist, i.e., chemical bonds based on Cartedian coordinates in r. """ bonds = [] nAtoms = len(r) for i in range(nAtoms): for j in range (i+1,nAtoms): distance = np.linalg.norm(r[i] - r[j]) if distance < maxDist: bonds.append(((i,j),distance)) # bondcenters.append( (coords[i] + coords[j]) /2) # Todo: write bondlengths on bond return bonds def tabulateStructure(atoms,coords, maxDist=2): """ Prints all angles between three bound atoms. Variable atoms contains atomic composition as string, e.g. 'CCHHHH', r contains Cart coordinates [[x1,y1,z,],[x2,...]...]. """ atoms = np.array([a for a in atoms]) # Convert string to numpy array coord = np.array(coords) # Atom coordinates nAtoms = atoms.size unitVector = lambda vector: vector / np.linalg.norm(vector) i,j,k = 0,0,0 returnString = f'{"-"*21}\nBond lengths r (Ang)\n{"-"*21}\n' print(f'nAtoms: {nAtoms}') for i in range(nAtoms): # List bond lengths ci = coord[i] for j in range (i+1, nAtoms): cj = coord[j] distance_ij = np.linalg.norm(ci - cj) if distance_ij < maxDist: returnString += f'r({atoms[i]}{i}-{atoms[j]}{j}) {distance_ij:.4f}\n' returnString += f'{"-"*28}\n\n{"-"*28}\nBond angles a (deg)\n{"-"*28}\n' for i in range(nAtoms): # List bond angles ci = coord[i] for j in range (i+1, nAtoms): cj = coord[j] distance_ij = np.linalg.norm(ci - cj) for k in range(j+1, nAtoms): ck = coord[k] distance_jk = np.linalg.norm(cj - ck) distance_ik = np.linalg.norm(ci - ck) if (distance_ij < maxDist and distance_ik < maxDist and j!=k): angle = np.arccos(np.clip(np.dot(unitVector(cj-ci), unitVector(ck-ci)), -1.0, 1.0)) returnString += 'a(%s%i-%s%i-%s%i) %.2f degrees\r\n' %(atoms[j],j, atoms[i],i, atoms[k],k, np.degrees(angle)) elif (distance_ij < maxDist and distance_jk < maxDist and j!=k): angle = np.arccos(np.clip(np.dot(unitVector(ci-cj), unitVector(ck-cj)), -1.0, 1.0)) returnString += 'a(%s%i-%s%i-%s%i) %.2f degrees\r\n' %(atoms[i],i, atoms[j],j, atoms[k],k, np.degrees(angle)) elif (distance_ik < maxDist and distance_jk < maxDist and j!=k): angle = np.arccos(np.clip(np.dot(unitVector(ci-ck), unitVector(cj-ck)), -1.0, 1.0)) returnString += 'a(%s%i-%s%i-%s%i) %.2f degrees\r\n'%(atoms[i],i, atoms[k],k, atoms[j],j, np.degrees(angle)) returnString += f'----------------------------\n' return returnString # (2.6) Plotting of molecules #------------------------------------------------------- # To check int2cart conversion functions and guess parameters, # I implement a simplistic molecule viewer to plot a ball and stick model. # class plotMolecule(QtGui.QApplication): # """ Simple plot of molecular structure based on pytgraph GLViewWidget. # Molecule must be specified as numpy array with strings [[atom,x,y,z]], # where atoms contains the atom label ('H','C',...) and x,y,z are Cartesian coordinates # or atoms contains [atomlist, coordinates] with float coordinates. # Example for use (requires 'import filetools as ft' to load ab initio results): # molecule = ft.Load_Orca(r'D:\TS\abinitio\substituted_benzenes\pyridine_dimer\H-Bound_Results\neutrals') # p = plotMolecule(molecule.atoms) # """ # def __init__(self, atoms): # QtGui.QApplication.__init__(self, sys.argv) # try: # Handle data from int2cart functions or from ab initio calculations # if len(atoms[0]) == 4: # Assume Cartesian coordinate list, e.g., [['C' '1.1' '0.2' '3.1'], [...]] # print('Cartesian coordinate list') # self.atoms = atoms # self.labels = atoms[:,0] # Atom symbols # self.coords = atoms[:,1:].astype(float) # Atom coordinates # elif len(atoms) == 2: # Assume atom string and separate coordinate list, e.g. ['OCCCCHHHH',[1.1, 0.0, 0.0], [...]] # print('molecular string, coordinate list') # lab, coo = atoms # print(f'labels: {lab}, coordinates{coo}') # self.labels = npa([a for a in lab]) # self.coords = npa(coo) # self.atoms = np.column_stack((npa([a for a in lab]), npa(coo).astype('|S8'))) # else: # Failed to recognize molecualr coordinates. # print('Atom coordinates not recognized') # return # except: # print('Exception: Atom coordinates not recognized') # return # # Some atom colors in (r,g,b,opacity). I tried to follow common plotting colors as used in Wikipedia. # # ToDo: adapt Jmol colors from https://sciencenotes.org/molecule-atom-colors-cpk-colors/ for complete periodic table? # atomColors = {'H': (0.92,0.92,0.92,1), 'B': (1.00,0.76,0.76,1), 'C': (0.25,0.25,0.25,1),'N': (0.25,0.25,0.87,1), # 'O': (0.95,0.00,0.00,1), 'Si':(0.57,0.66,0.66,1), 'P': (1.00,0.56,0.00,1), 'S': (1.00,1.00,0.00,1), # 'F': (0.86,0.95,0.31,1), 'Cl':(0.18,0.95,0.18,1), 'Br':(0.70,0.22,0.22,1), 'I': (0.64,0.00,0.64,1),} # # Some covalent Radii from [Dalton Trans., DOI: 10.1039/b801115j] # # ToDo: use covalent radii for bond recognition. # atomSizes = {'H': 0.31, 'B': 0.84, 'C': 0.76, 'N': 0.71, 'O': 0.66,'F': 0.57, 'Si':1.11, # 'P': 1.07, 'S': 1.05, 'Cl':1.02, 'Br':1.20, 'I': 1.39,} # self.w = gl.GLViewWidget() # self.w.setAttribute(QtCore.Qt.WA_DeleteOnClose) # self.w.show() # self.w.setWindowTitle('plotMolecule') # self.w.setCameraPosition(distance=15, azimuth=-90) # #cyl = gl.MeshData.cylinder(rows=10, cols=10, radius=[0.2, 0.2], length=1.0,) # Plot bonds as cylinders? # for bond in (findBonds(self.atoms)): # v1,v2 = self.coords[bond[0]], self.coords[bond[1]] # Vectors to bonded atoms # b = gl.GLLinePlotItem(pos=np.array([v1,v2]), color=(0.4,0.4,0.4,0.4), width=3) # self.w.addItem(b) # # need to debug the bond rotation ... # #b = gl.GLMeshItem(meshdata=cyl, smooth=True, shader='shaded', color=(1,0.0,0.0,0.5), glOptions='opaque') # #b.scale(1,1,np.linalg.norm(v2-v1)) # bond length # #b.rotate(np.arctan2( (v2[0]-v1[0]),(v2[2]-v1[2]) )/np.pi*180, 0,1,0) # Rotate around y in xz plane # #b.rotate(np.arctan2( (v2[1]-v1[1]),(v2[0]-v1[0]) )/np.pi*180 , 0,0,1) # Rotate around z in xy plane # #b.translate(*v1) # sphere = gl.MeshData.sphere(rows=10, cols=10, radius=0.7) # Plot atoms as spheres, scale covalent radii x0.7 # for atom,coord in zip(self.labels,self.coords): # a = gl.GLMeshItem(meshdata=sphere, smooth=True, shader='shaded', color=atomColors.get(atom,(0.5,0.5,0.5)), glOptions='opaque') # a.translate(*coord) # translate to atom position # a.scale(*[atomSizes[atom]]*3) # Scale proportional to atom size # self.w.addItem(a) def plotMolecule(molecule): """ Simple plot of molecular structure based on pytgraph GLViewWidget. Molecule must be specified as numpy array with strings [[atom,x,y,z]], where atoms contains the atom label ('H','C',...) and x,y,z are Cartesian coordinates or atoms contains [atomlist, coordinates] with float coordinates. Example for use (requires 'import filetools as ft' to load ab initio results): molecule = ft.Load_Orca(r'D:\TS\abinitio\substituted_benzenes\pyridine_dimer\H-Bound_Results\neutrals') p = plotMolecule(molecule.atoms) """ try: # Handle data from int2cart functions or from ab initio calculations if len(molecule[0]) == 4: # Assume Cartesian coordinate list, e.g., [['C' '1.1' '0.2' '3.1'], [...]] print('Cartesian coordinate list') atoms = molecule labels = atoms[:,0] # Atom symbols coords = atoms[:,1:].astype(float) # Atom coordinates elif len(molecule) == 2: # Assume atom string and separate coordinate list, e.g. ['OCCCCHHHH',[1.1, 0.0, 0.0], [...]] print('molecular string, coordinate list') lab, coo = molecule print(f'labels: {lab}, coordinates\n{coo}') labels = npa([a for a in lab]) coords = npa(coo) atoms = np.column_stack((npa([a for a in lab]), npa(coo).astype('|S8'))) else: # Failed to recognize molecualr coordinates. print('Atom coordinates not recognized') return except: print('Exception: Atom coordinates not recognized') return # Some atom colors in (r,g,b,opacity). I tried to follow common plotting colors as used in Wikipedia. # ToDo: adapt Jmol colors from https://sciencenotes.org/molecule-atom-colors-cpk-colors/ for complete periodic table? atomColors = {'H': (0.92,0.92,0.92,1), 'B': (1.00,0.76,0.76,1), 'C': (0.25,0.25,0.25,1),'N': (0.25,0.25,0.87,1), 'O': (0.95,0.00,0.00,1), 'Si':(0.57,0.66,0.66,1), 'P': (1.00,0.56,0.00,1), 'S': (1.00,1.00,0.00,1), 'F': (0.86,0.95,0.31,1), 'Cl':(0.18,0.95,0.18,1), 'Br':(0.70,0.22,0.22,1), 'I': (0.64,0.00,0.64,1),} # Some covalent Radii from [Dalton Trans., DOI: 10.1039/b801115j] # ToDo: use covalent radii for bond recognition. atomSizes = {'H': 0.31, 'B': 0.84, 'C': 0.76, 'N': 0.71, 'O': 0.66,'F': 0.57, 'Si':1.11, 'P': 1.07, 'S': 1.05, 'Cl':1.02, 'Br':1.20, 'I': 1.39,} w = gl.GLViewWidget() w.setAttribute(QtCore.Qt.WA_DeleteOnClose) w.show() w.setWindowTitle('plotMolecule') w.setCameraPosition(distance=15, azimuth=-90) #cyl = gl.MeshData.cylinder(rows=10, cols=10, radius=[0.2, 0.2], length=1.0,) # Plot bonds as cylinders? for bond in (findBonds(atoms)): v1,v2 = coords[bond[0][0]], coords[bond[0][1]] # Vectors to bonded atoms b = gl.GLLinePlotItem(pos=np.array([v1,v2]), color=(0.4,0.4,0.4,0.4), width=3) w.addItem(b) # need to debug the bond rotation ... #b = gl.GLMeshItem(meshdata=cyl, smooth=True, shader='shaded', color=(1,0.0,0.0,0.5), glOptions='opaque') #b.scale(1,1,np.linalg.norm(v2-v1)) # bond length #b.rotate(np.arctan2( (v2[0]-v1[0]),(v2[2]-v1[2]) )/np.pi*180, 0,1,0) # Rotate around y in xz plane #b.rotate(np.arctan2( (v2[1]-v1[1]),(v2[0]-v1[0]) )/np.pi*180 , 0,0,1) # Rotate around z in xy plane #b.translate(*v1) sphere = gl.MeshData.sphere(rows=10, cols=10, radius=0.7) # Plot atoms as spheres, scale covalent radii x0.7 for atom,coord in zip(labels,coords): a = gl.GLMeshItem(meshdata=sphere, smooth=True, shader='shaded', color=atomColors.get(atom,(0.5,0.5,0.5)), glOptions='opaque') a.translate(*coord) # translate to atom position a.scale(*[atomSizes[atom]]*3) # Scale proportional to atom size w.addItem(a) return w # (2.6) Residuals for Geometry fitting #------------------------------ # The residual function will be minimized in the fit and therefore determines the nature of the fit. # I define a number of different residual functions corresponding to common structure fit methods # described in the literature. A variable array /r/ holds the fit parameter list, beginning with the # internal structure parameters required for int2cart coordinate transformation, followed by # optional ci, di, and Laurie parameter dH (as last value): r = [r1,r2,...,ca,cb,...,(dH)]. # nPars specifies the number of parameters and the molecule symmetry to facilitate the use of reduced # parameter sets. # Symmetry-dependent unpacking of rovib correction parameters. def roVibParameters(par,nPars): """ Unpack rovibrational correction parameters based on molecular symmetry. return cij-tensor (3x3), di vector (3), Laurie correction. """ ni = nPars['int'] nc = nPars['ci'] nd = nPars['di'] sym = nPars['sym'] Laurie = nPars['Laurie'] di=[] if ('Linear' in sym or 'Spherical' in sym): # Only 1 parameters ci,di for rovib correction c = par[ni] cij = np.diag([c,c,c]) # cij tensor if nd >0: di = npa([par[ni+nc]]*3) # di vector present #print(f'linear') elif 'Planar Oblate' in sym: # Ia=Ib, ca=cb, planar: no cac,cbc ca = par[ni] cab = ca # assume ca=cb=cab if nc == 1: cc = ca*2 # Ic=2*Ia; assumption of cc=2*ca if nc > 1: cc = par[ni+1] # Separate diagonal element cc (ca,ca,cc) if nc == 3: cab = par[ni+2] # separate off-diagonal cij element cij = npa([[ca,cab,0],[cab,ca,0],[0,0,cc]]) # cij tensor if nd==1: di = npa([par[ni+nc],par[ni+nc],-par[ni+nc]]) # di vector (da,da,-da) if nd==2: di = npa([par[ni+nc],par[ni+nc],par[ni+nc+1]])# di vector (da,da,dc) #print(f'planar oblate, cij:{cij}, di:{di}') elif 'Oblate Top' in sym: # Ia=Ib ca,cc = par[ni:ni+2] # diagonal elements cii (ca,ca,cc) if nc == 5: cab,cac,cbc = par[ni+2:ni+5] # off-diagonal cij elements else: cab,cac,cbc = ca,(ca+cc)/2,(ca+cc)/2 # assume symmetry / average for cij cij = npa([[ca,cab,cac],[cab,ca,cbc],[cac,cbc,cc]]) # cij tensor if nd >0: di = npa([par[ni+nc],par[ni+nc],par[ni+nc+1]])# di vector (da,da,dc) #print(f'oblate, cij:{cij}') elif 'Prolate Top' in sym: # Ib=Ic ca,cb = par[ni:ni+2] # diagonal elements cii (ca,cb,cb) if nc == 5: cab,cac,cbc = par[ni+2:ni+5] # off-diagonal cij elements else: cab,cac,cbc = (ca+cb)/2,cb,cb # assume symmetry / average for cij cij = npa([[ca,cab,cac],[cab,cb,cbc],[cac,cbc,cb]]) # cij tensor if nd >0: di = npa([par[ni+nc],par[ni+nc+1],par[ni+nc+1]]) # di vector (da,db,db) print(f'prolate, cij:{cij}') elif 'Planar Molecule' in sym: # planar: no cac,cbc ca,cb,cc = par[ni:ni+3] # diagonal elements cii (ca,cb,cc) if nc == 4: cab = par[ni+3] # off-diagonal cij elements else: cab = 0 # Ignore off-diagonal element cij = npa([[ca,cab,0],[cab,cb,0],[0,0,cc]]) # cij tensor if nd >0: di = npa([par[ni+nc],par[ni+nc+1],par[ni+nc+2]]) # di vector (da,db,dc) #print(f'planar, cij:{cij}') else: # Non-symmetric molecule: full parameter set ca,cb,cc = par[ni:ni+3] # diagonal elements cii (ca,cb,cc) cij = np.diag([ca,cb,cc]) # diagonal cij tensor if nc == 6: # off-diagonal cij elements cab,cac,cbc = par[ni+3:ni+6] cij = npa([[ca,cab,cac],[cab,cb,cbc],[cac,cbc,cc]]) # cij tensor if nd >0: di = npa([par[ni+nc],par[ni+nc+1],par[ni+nc+2]]) # di vector (da,db,dc) #print(f'asymmetric, cij:{cij}') return cij,di,Laurie # (2.6a) Straightforward _r0_ fit #----------------------------------------------------------------------- # This is a direct fit of calculated inertial moments I0 to measured inertial moments Iexp, # assuming that geometry parameters are isotope-independent. Resulting bond parameters # correspond to 'effective' bond lengths. def Residuals_r0(r,isoList, nPars, int2cart=bz_internal2cart): """ Return residuals (I0 - PMI) for r_0 structure determination. (eq. 1 in Wat99) The structure is fitted to Iexp = I(r_0) with internal structure parameters r = [r1,r2,...,(dH)] as required by int2cart. """ residuals = [] ni,Laurie = nPars['int'],nPars['Laurie'] # Number of internal and Laurie parameters rc = int2cart(r,isoList[0][-1],Laurie) # Convert internal to Cartesian coordinates if DEBUG: log.write(f'Reference Coordinates: {rc[0]},{rc[1]},{rc[2]}') if DEBUG: log.write(f'Residuals_r0, Laurie: {Laurie}\n') for iso in isoList: # Iterate over isotopologues if Laurie: rc = int2cart(r,iso[-1],Laurie) # Convert internal to Cartesian coordinates with Laurie correction Iexp = iso[1][0] # Measured inertial moments weight = npa(iso[2]) # Fit weights I0,eVec,sym,I = PMI(rc,iso[-1]) # Calculate principal moments I0 from atom position and mass residuals.extend(( (Iexp - I0) *weight) [weight>0]) # Add weighted residuals to list if DEBUG: log.write(f'{iso[0]}, PMI:{I0}\n') return residuals # (0.0) Implementation of Typke fit (removed) # This is very simple to implement (see old versions of Benzene_isotopologue_geometry fit) -- # but requires a complete set of A,B,C rotational moments for the calculation of planar moments. # For symmetric molecules, missing rotational moments may be calculated (e.g., ignoring the inertial # defect for benzene we have A=B, C=2A for a planar oblate top). I didn't finish generalizing the # calculation of missing moments for all symetry cases and I suggest to just use the fitting of # Residuals_rs_differences, which corresponds to the same fit but implemented with principal moments # instead of planar moments. The two types of fit are mathematically identical. # def Residuals_rs_Typke(r,isoList, nPars, int2cart=bz_internal2cart): #-- Not Completed # """ Residual function for Typke type structure determination. # (Fitting of planar moment differences between a parent isotopologue # and all child isotopologues. The result is identical to that of a direct # Kraitchman analysis but allows to fit partial or overdetermined # data sets. # r: internal coordinates as required by int2cart functions. # """ # residuals = [] # ni,Laurie = nPars['int'],nPars['Laurie'] # Number of internal and Laurie parameters # if len(isoList[0]) < 8: # First run, calculate planar moments # print('Calculating planar moments') # for iso in isoList: # sym = iso[4] # Rot. symmetry # weight = npa(iso[2]) # Fit weights # Ia,Ib,Ic = npa(iso[1][0]) [weight>0] # Measured inertial moments # if not Ia: # Need to fill Ia # if sym in ['spherical','oblate','oblate pl.'] and Ib: Ia = Ib # elif sym in ['spherical'] and Ic: Ia = Ic # if sym in ['planar'] and Ib and Ic: Ia = Ic-Ib # if not Ib: # if sym in ['spherical','oblate','oblate pl.'] and Ia: Ib = Ia # elif sym in ['spherical'] and Ic: Ib = Ic # if sym in ['planar'] and Ia and Ic: Ib = Ic-Ia # if sym in ['linear','prolate']: Ib = Ic # if not Ic: # if sym in ['spherical','prolate','linear'] and Ib: Ic = Ib # elif sym in ['spherical'] and Ia: Ic = Ia # if sym in ['planar','oblate pl.'] and Ia and Ib: Ic = Ia+Ib # if Ia and Ib and Ic: # All principal axis moments available # Pa = (-Ia+Ib+Ic)/2 # --> calculate planar moments # Pb = ( Ia-Ib+Ic)/2 # Pc = ( Ia+Ib-Ic)/2 # iso.append([Pa,Pb,Pc]) # Append planar moments as col. 7 in isoList # iso = isoList[0] # Analyze parent molecule (first in list) # weight = npa(iso[2]) # Fit weights # P_exp_parent = iso[7] # Experimental planar moments # if np.any(npa(P_exp_parent)==0): print('Please select parent with nonzero Ia, Ib, Ic.'); return # atoms, m_parent = iso[0], npa(iso[6]) # rc = int2cart(r,m_parent,Laurie=0) # Convert parent internal to Cartesian coordinates # rcom = COM(rc,m_parent) # Convert parent internal to Cartesian coordinates # P_i_parent = m_parent* rcom.transpose()**2 # Atomic contributions to planar moment from center-of-mass coordinates and masses # P_calc_parent = np.sum(P_i_parent, axis=1) # Planar moments (x,y,z array), diagonal elements # residuals.extend( ((P_exp_parent - P_calc_parent)*weight)[weight>0] ) # Add residual to list # children = isoList[1:] # for iso in children: # weight = npa(iso[2]) # Fit weights # P_exp_child = iso[7] # Experimental planar moments # m_child = npa(iso[6]) # Child atom masses # dm_i = m_child - m_parent # Isotope mass differences child-parent # mu_i = sum(m_parent)*dm_i / (sum(m_parent)+dm_i) # Reduced mass differences (==0 when unsubstituted) # #dP_i = np.sum( [m*np.outer(r,r) for m,r in zip(mu,r_i)], axis=0) # (TYP 4); but we only need diag. elements # dP_i = mu_i * (rcom.transpose())**2 # Planar moment difference child-parent # P_calc_child = np.sum(P_i_parent + dP_i, axis=1) # Planar moments (x,y,z array), diagonal elements (TYP 4) # residuals.extend( ((P_exp_child - P_calc_child) *weight)[weight>0] ) # Add residual to list # return residuals # Return residuals # Untestested! # (2.6b) Nossberger-type fit of substitution structure rs # -- Beware when fitting more than one axis! Identical results to Residuals_rs. #----------------------------------------------------------------------- # A constant rovibrational coupling term ci is implied or fitted for each principal I tensor # element Iij. The original Nössberger approach fitted inertial moment differences between # isotopologue moments and a parent moment, which cancels the ci terms. Fit results then # depend on the choice of parent. This is equivalent to a Kraitchman analysis, but can # accounts for all available isotopologue inertial moments in a single fit instead of # pairwise isotopologue comparisons. Fitted parameters may be expected to agree with r0 # and re values within the Costain errors. The 'rs' fit explicitly fits ci terms and the # resulting fit does not depend on the choice of a parent species. The 'rsr' fit implements # principal axis rotation into the parent principal axis system before rovibrational # correction, so correction terms are applied along invariant intramolecular axes. def Residuals_rs_Nossberger(r,isoList, nPars, int2cart=bz_internal2cart): """ Return residuals for r_s structure determination (fitting of inertial moment differences) as proposed by Nössberger, Bauder and Günthardt [ChemPhys 1, 418 (1973)]. The structure is fitted to Iexp-pIexp = I(r)-pI(r) with p indicating the parent species values. Note that rovib. correction terms ci cancel in the difference I0(r)+ci - ( pI0(r)+ci ). r holds the fit parameter list, beginning with the internal structure parameters required for int2cart coordinate transformation, followed by ci parameters and with the optional Laurie parameter dH as last value. (r = [r1,r2,...,ca,cb,cc,(dH)]). """ residuals = [] Laurie = nPars['Laurie'] # Laurie correction is last parameter #--- Reference isotopologue --- iso = isoList[nPars['parent']] # Reference molecule (fit I-I_ref) rc = int2cart(r,iso[-1],Laurie) # Convert parent internal to Cartesian coordinates if DEBUG: log.write(f'Reference Coordinates: {rc[0]},{rc[1]},{rc[2]}') parent = iso[0] pIexp = iso[1][0] pWeight = npa(iso[2]) # Parent fit weights pPmi,*rest = PMI(rc,iso[-1]) # Evaluate parent / reference PMI if DEBUG: log.write(f'pPmi: {pPmi}, pIexp: {pIexp}\n') #--- Calculate residuals for all isotopologues --- for iso in isoList: # Iterate over isotopologues if Laurie: # Isotope-dependent bond lengths? rc = int2cart(r,iso[-1],Laurie) # -> Convert int2cart for each isotopologue with Laurie correction Iexp = iso[1][0] # Measured inertial moments, atom masses weight = npa(iso[2]) # Fit weights pmi,*rest = PMI(rc,iso[-1]) # Calculate principal moments I0(r) residuals.extend( (((Iexp - pmi) - (pIexp-pPmi)) *weight) [(weight>0) & (pWeight>0)]) # Add to residuals if DEBUG: log.write(f'{iso[0]}, Iexp:{Iexp}, PMI:{pmi}, pIexp:{pIexp}, pPMI:{pPmi}, residuals:{residuals[-3:]}\n') return residuals def Residuals_rs(r,isoList, nPars, int2cart=bz_internal2cart): """ Return residuals for r_s structure determination. Rovibration coupling is corected by a constant ci correction for each principal axis and the structure is fitted to Iexp = I0(r) + c (eqs. 20, 21 in Wat99). As opposed to the Nössberger fit, this fit does not depend on the choice of a parent species. """ residuals = [] cij,_,Laurie = roVibParameters(r,nPars) # Unpack rovibrational parameters #--- Reference isotopologue: check for MI axis rotation --- iso = isoList[nPars['parent']] # Reference molecule (analyze PMI rotation) rc = int2cart(r,iso[-1],Laurie) # Convert parent internal to Cartesian coordinates I0,eVec_ref,sym_ref,I = PMI(rc,iso[-1]) # Evaluate reference PMI if DEBUG: log.write(f'Reference Coordinates: {rc[0]},{rc[1]},{rc[2]}') if DEBUG: log.write(f'Parent symmetry: {sym_ref} (2:sym,aligned, 1:sym, 0:asym\n') if DEBUG: log.write(f'Residuals_rs, Laurie: {Laurie}, ci: {np.diag(cij)}n') #--- Calculate residuals for all isotopologues --- for iso in isoList: # Iterate over isotopologues if Laurie: # Isotope-dependent bond lengths? rc = int2cart(r,iso[-1],Laurie) # -> Convert int2cart for each isotopologue with Laurie correction Iexp = iso[1][0] # Measured inertial moments, atom masses weight = npa(iso[2]) # Fit weights Map = iso[6] # Mapping of a,b,c to reference x,y,z I0,eVec,sym,I = PMI(rc,iso[-1]) # Calculate principal moments I0(r) Is = I0 + np.diag(cij)[Map] # Inertial moment Im, rm1-corrected residuals.extend( ((Iexp - Is) *weight) [weight>0]) # Add to residuals if DEBUG: log.write(f'{iso[0]}, PMI:{I0}, Is:{Is}\n') return residuals def Residuals_rsr(r,isoList, nPars, int2cart=bz_internal2cart): """ Return residuals for r_sr structure determination. As compared to the rs fit (above), the rovibrational correction is applied to the internal structure parameters required for int2cart coordinate transformation, followed by ci parameters and with the optional Laurie parameter dH as last value. (r = [r1,r2,...,ca,cb,cc,(dH)]). """ residuals = [] cij,_,Laurie = roVibParameters(r,nPars) # Unpack rovibrational parameters #--- Reference isotopologue: check for MI axis rotation --- iso = isoList[nPars['parent']] # Reference molecule (analyze PMI rotation) rc0 = int2cart(r,iso[-1],Laurie) # Convert parent internal to Cartesian coordinates I0,eVec_ref,sym_ref,I = PMI(rc0,iso[-1]) # Evaluate reference PMI if 1: #not sym_ref==2: # Reference coordinate system is not principal axis system (PAS)? eInv_ref = np.dot(np.linalg.inv(eVec_ref),nPars['rotMat']) # -> Rotation matrix into reference PAS + user supplied angle else: eInv_ref = nPars['rotMat'] # -> Rotation matrix into reference PAS + user supplied angle rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS if DEBUG: log.write(f'Reference Coordinates: {rc[0]},{rc[1]},{rc[2]}') if DEBUG: log.write(f'Parent symmetry: {sym_ref} (2:sym,aligned, 1:sym, 0:asym\n') if DEBUG: log.write(f'Residuals_rsNr, Laurie: {Laurie}, cij: {cij[0]},{cij[1]},{cij[2]}\n') #--- Calculate residuals for all isotopologues --- for iso in isoList: # Iterate over isotopologues if Laurie: # Isotope-dependent bond lengths? rc0 = int2cart(r,iso[-1],Laurie) # -> Convert int2cart for each isotopologue with Laurie correction rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS Iexp = iso[1][0] # Measured inertial moments, atom masses weight = npa(iso[2]) # Fit weights Map = iso[6] # Mapping of a,b,c to reference x,y,z I0,eVec,sym,I = PMI(rc,iso[-1]) # Calculate principal moments I0(r) rotIs = I + cij # cij corrections to I tensor in reference coordinate system Is,e_vectors = np.linalg.eigh(rotIs) # Principal moments for corrected I tensor residuals.extend( ((Iexp - Is) *weight) [weight>0]) # Add weighted residuals to list if DEBUG: log.write(f'{iso[0]}, PMI:{I0}, Is:{Is}\n') return residuals # Watson's mass-dependent _rmIr_ fit #---------------------------------------------------------------- # The rmI fit scales the rovibrational coupling correction terms ci by sqrt(I). This # accounts for the expected first-order mass-dependence of the correction terms. # The mass dependence sqrt(I) = sum(sqrt(m)*r) approximates the expected dependence # For harmonic vibration, where the average ZPV displacement is proportional to sqrt(m): # I0 = sum(m*r0**2) = sum(m*(re+dr)**2)) = sum(m*re**2) + sum(m*2*re*dr) + sum(m*dr**2) # Ignoring the last term and substituting dr=k/sqrt(m), we find: I0 = Ie + sum(k/2*sqrt(Ie)). # We can fit the resulting I0 to I0_exp and obtain an approximation for Ie (and re). We call # the resulting values Im and rm to indicate the approximate nature of the solution: # I0 = Im + c*sqrt(Im) # The fitted ci correction factors represent the sum of the rovib. corrections between I0 and Im. # Note that the fit will optimize ci parameters not to match the actual harmonic rovib. effects, # but to account for all effects (harm., anharm, Coriolis). # The rmIr fit accounts for PMI axis rotation: If the isotopologue principal # moments are not parallel, then the rovibrational correction terms ci*sqrt(I) are # applied along different molecular axis orientations. This can be corrected by rotating the # isotopologue principal moments into a common coordinate system before applying corrections # to the I tensor. The I tensor is then rediagonalized (rotated to its pricipal axes) for the # fit. Because we correct the 3x3 I tensor instead of the 3 principal moments, we require # additional ci coefficients cab (or cab,cac,cbc) for planar (or non-planar) molecules. # This reproduces the STRFIT rmIr result (STRFIT always rotates to the parent axes). # When including Laurie corrections, we cannot perfectly reproduce STRFIT results (e.g., (10 muA, # .01 deg difference in fitted parameters for H2O), but we reproduce values in Watson1999 exactly. def Residuals_rmI(r,isoList, nPars, int2cart=bz_internal2cart): """ Residual function for r_m(I) structure determination. The rovibration coupling term is approximated by the mass-dependent correction terms [ci*sqrt(I)]. The structure is fitted to Iexp = PMI(r) + c*sqrt(PMI(r)) (eq. 22 in Wat99) with calculated rigid rotor principal moments I = PMI(r). There is one parameter ci for each axis (a,b,c). Parameters are r = rH,rC, ca,cb,cc. """ residuals = [] cij,_,Laurie = roVibParameters(r,nPars) # Unpack rovibrational parameters #--- Reference isotopologue: check for MI axis rotation --- iso = isoList[nPars['parent']] # Reference molecule (analyze PMI rotation) rc0 = int2cart(r,iso[-1],Laurie) # Convert parent internal to Cartesian coordinates pmi,eVec_ref,sym_ref,I = PMI(rc0,iso[-1]) # Evaluate reference PMI eInv_ref = np.dot(np.linalg.inv(eVec_ref),nPars['rotMat']) # -> Rotation matrix into reference PAS + user supplied angle rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS if DEBUG: log.write(f'Reference Coordinates: {rc0[0]},{rc0[1]},{rc0[2]}') if DEBUG: log.write(f'Parent symmetry: {sym_ref} (2:sym,aligned, 1:sym, 0:asym\n') if DEBUG: log.write(f'Residuals_rmI, Laurie: {Laurie}, cij: {cij[0]},{cij[1]},{cij[2]}\n') #--- Calculate residuals for all isotopologues --- for iso in isoList: # Iterate over isotopologues if Laurie: # Isotope-dependent bond lengths? rc0 = int2cart(r,iso[-1],Laurie) # -> Convert int2cart for each isotopologue with Laurie correction rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS Iexp = iso[1][0] # Measured inertial moments, atom masses weight = npa(iso[2]) # Fit weights Map = iso[6] # Mapping of a,b,c to reference x,y,z pmi,eVec,sym,I = PMI(rc,iso[-1]) # Calculate principal moments I0(r) cTerm = np.diag(cij)[Map]*np.sqrt(pmi) # First order (rmI) correction term (only diag. elements) I0 = pmi + cTerm # Inertial moment Im, rm1-corrected. residuals.extend( ((Iexp - I0) *weight) [weight>0]) # Add to residuals if DEBUG: log.write(f'{iso[0]}, I0:{I0}\n, PMI:{pmi}, cTerm:{cTerm}') return residuals def Residuals_rmIr(r,isoList, nPars, int2cart=bz_internal2cart): """ Residual function for r_m(I) structure determination. The rovibration coupling term is approximated by the mass-dependent correction terms [ci*sqrt(I)]. The structure is fitted to Iexp = PMI(r) + c*sqrt(PMI(r)) (eq. 22 in Wat99) with calculated rigid rotor principal moments I = PMI(r). There is one parameter ci for each axis (a,b,c). Parameters are r = rH,rC, ca,cb,cc. """ residuals = [] cij,_,Laurie = roVibParameters(r,nPars) # Unpack rovibrational parameters #--- Reference isotopologue: check for MI axis rotation --- iso = isoList[nPars['parent']] # Reference molecule (analyze PMI rotation) rc0 = int2cart(r,iso[-1],Laurie) # Convert parent internal to Cartesian coordinates pmi,eVec_ref,sym_ref,I = PMI(rc0,iso[-1]) # Evaluate reference PMI eInv_ref = np.dot(np.linalg.inv(eVec_ref),nPars['rotMat']) # -> Rotation matrix into reference PAS + user supplied angle rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS if DEBUG: log.write(f'Reference Coordinates: {rc[0]},{rc[1]},{rc[2]}') if DEBUG: log.write(f'Parent symmetry: {sym_ref} (2:sym,aligned, 1:sym, 0:asym\n') if DEBUG: log.write(f'Residuals_rmIr, Laurie: {Laurie}, cij: {cij[0]},{cij[1]},{cij[2]}\n') #--- Calculate residuals for all isotopologues --- for iso in isoList: # Iterate over isotopologues if Laurie: # Isotope-dependent bond lengths? rc0 = int2cart(r,iso[-1],Laurie) # -> Convert int2cart for each isotopologue with Laurie correction rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS Iexp = iso[1][0] # Measured inertial moments, atom masses weight = npa(iso[2]) # Fit weights Map = iso[6] # Mapping of a,b,c to reference x,y,z pmi,eVec,sym,I = PMI(rc,iso[-1]) # Calculate principal moments pmi(r) r_sqrtPMI = rotI(np.sqrt(pmi),eVec) # Rotate sqrt(pmi) to reference PAS rI0 = I + cij*r_sqrtPMI # rmI correction to I tensor in reference PAS I0,e_vectors = np.linalg.eigh(rI0) # Principal moments for corrected I tensor residuals.extend( ((Iexp - I0) *weight) [weight>0]) # Residuals if DEBUG: log.write(f'asym:{iso[0]}, I0:{I0}\n, PMI:{pmi}, cTensor:{(cij*r_sqrtPMI)[0]},{(cij*r_sqrtPMI)[1]},{(cij*r_sqrtPMI)[2]}') return residuals # Watson's mass-dependent _rmIIr_ fit #-------------------------------------------------------------------- # The rmII fit includes first and second order rovibrational correction terms based # on parameters ci and di. The mass dependence for the di-terms are given by Wat99 # as (prod(m)/sum(m))**1/(2N-2) and the fit minimizes the difference between # Iexp and I0 with I0 = Im + c*sqrt(Im) + d*(prod(m)/sum(m))**1/(2N-2). # For a description of PMI rotation ('r'), refer to the rmI # fit description, above. The rmIIr fits reproduce STRFIT rmIIr results (STRFIT always # rotates to parent axes before rovibbrational correction). Laurie correction leads # to a small difference (40 muA, .02 deg) between our and the STRFIT results, but gives # an exact reproduction of results published in Watson1999. def Residuals_rmII(r,isoList, nPars, int2cart=bz_internal2cart): """ Residual function for r_m(II) structure determination. The structure is fitted to (Watson1999, eq. 31): Iexp = Im + c*sqrt(Im_pq) + d*(m1*m2*...)/sum(m1,m2,...) Parameters c and d are constants for each axis (a,b,c); corrections are applied along the a,b,c principal axis of each isotopologue. Parameters are r = [r(OH),angle(HOH), ca,cb,cc,..., da,db,dc, dH]. """ residuals = [] cij,di,Laurie = roVibParameters(r,nPars) # Unpack rovibrational parameters #--- Reference isotopologue: check for MI axis rotation --- iso = isoList[nPars['parent']] # Reference molecule (analyze PMI rotation) rc0 = int2cart(r,iso[-1],Laurie) # Convert parent internal to Cartesian coordinates pmi,eVec_ref,sym_ref,I = PMI(rc0,iso[-1]) # Evaluate reference PMI eInv_ref = np.dot(np.linalg.inv(eVec_ref),nPars['rotMat']) # -> Rotation matrix into reference PAS + user supplied angle rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS if DEBUG: log.write(f'Reference Coordinates: {rc[0]},{rc[1]},{rc[2]}') if DEBUG: log.write(f'Parent symmetry: {sym_ref} (2:sym,aligned, 1:sym, 0:asym\n') if DEBUG: log.write(f'Residuals_rmII, Laurie: {Laurie}, cij: {cij[0]},{cij[1]},{cij[2]}, di:{di}\n') #--- Calculate residuals for all isotopologues --- for iso in isoList: # Iterate over isotopologues if Laurie: # Isotope-dependent bond lengths? rc0 = int2cart(r,iso[-1],Laurie) # -> Convert int2cart for each isotopologue rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS Iexp = iso[1][0] # Measured inertial moments weight = npa(iso[2]) # Fit weights m_i = iso[-1] # Atom masses Map = iso[6] # Mapping of a,b,c to reference x,y,z pmi,eVec,sym,I = PMI(rc,iso[-1]) # Calculate principal moments I0(r) dTerm = di*(np.prod(m_i)/np.sum(m_i))**(1/(2*len(m_i)-2)) # Second order (rmII) correction term # Watson1999 simplified the d-Term mass dependence. Here I checked how this approximation affected the fit results (small effect). # r_cbrtPMI = np.cbrt(pmi) # Cubed root (pmi) # dTerm = di* (np.prod(m_i)/np.sum(m_i)/r_cbrtPMI)**(1/(2*len(m_i)-4)) # d-Term using the proper Wilson's G-matrix mass dependence cTerm = np.diag(cij)[Map]*np.sqrt(pmi) # First order (rmI) correction term (only diag. elements) I0 = pmi + cTerm + dTerm[Map] # Inertial moment with rm1 and rm2 correction terms residuals.extend( ((Iexp - I0) *weight) [weight>0]) # Add to residuals if DEBUG: log.write(f'{iso[0]}, I0:{I0}\n, PMI:{pmi}, cTerm:{cTerm}, dTerm:{dTerm[Map]}') return residuals def Residuals_rmIIr(r,isoList, nPars, int2cart=bz_internal2cart): """ Residual function for r_m(II) structure determination. The structure is fitted to (Watson1999, eq. 31): Iexp = Im + c*sqrt(Im_pq) + d*(m1*m2*...)/sum(m1,m2,...) Parameters c and d are constants for each axis (a,b,c); corrections are applied along the parent principal axes. Parameters are r = [r(OH),angle(HOH), ca,cb,cc,..., da,db,dc, dH]. """ residuals = [] cij,di,Laurie = roVibParameters(r,nPars) # Unpack rovibrational parameters #--- Reference isotopologue: check for MI axis rotation --- iso = isoList[nPars['parent']] # Reference molecule (analyze PMI rotation) rc0 = int2cart(r,iso[-1],Laurie) # Convert parent internal to Cartesian coordinates pmi,eVec_ref,sym_ref,I = PMI(rc0,iso[-1]) # Evaluate reference PMI eInv_ref = np.dot(np.linalg.inv(eVec_ref),nPars['rotMat']) # -> Rotation matrix into reference PAS + user supplied angle rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS if DEBUG: log.write(f'Reference Coordinates: {rc[0]},{rc[1]},{rc[2]}') if DEBUG: log.write(f'Parent symmetry: {sym_ref} (2:sym,aligned, 1:sym, 0:asym\n') if DEBUG: log.write(f'Residuals_rmIIr, Laurie: {Laurie}, cij: {cij[0]},{cij[1]},{cij[2]}, di:{di}\n') #--- Calculate residuals for all isotopologues --- for iso in isoList: # Iterate over isotopologues if Laurie: # Isotope-dependent bond lengths? rc0 = int2cart(r,iso[-1],Laurie) # -> Convert int2cart for each isotopologue rc = rotR(rc0,eInv_ref) # -> Rotate rc into reference PAS Iexp = iso[1][0] # Measured inertial moments weight = npa(iso[2]) # Fit weights m_i = iso[-1] # Atom masses Map = iso[6] # Mapping of a,b,c to reference x,y,z pmi,eVec,sym,I = PMI(rc,m_i) # Calculate principal moments pmi(r) dTerm = di*(np.prod(m_i)/np.sum(m_i))**(1/(2*len(m_i)-2)) # second order (rmII) correction mass dependence # Watson1999 simplified the d-Term mass dependence. Here I checked if this approximation afected the fit results (it didn't). # r_cbrtPMI = np.cbrt(pmi) # Cubed root (pmi) # dTerm = di* (np.prod(m_i)/np.sum(m_i)/r_cbrtPMI)**(1/(2*len(m_i)-4)) # d-Term using the proper Wilson's G-matrix mass dependence r_sqrtPMI = rotI(np.sqrt(pmi),eVec) # Rotate sqrt(pmi) to reference PMI coordinate system rI0 = I + cij*r_sqrtPMI + np.diag(dTerm) # rmI,II corrections to I in reference PMI coordinate system I0,e_vectors = np.linalg.eigh(rI0) # Principal moments for corrected I tensor residuals.extend( ((Iexp - I0) *weight) [weight>0]) # Residuals if DEBUG: log.write(f'\nasym:{iso[0]}, I0:{I0}, PMI:{pmi}\n, cTensor:{(cij*r_sqrtPMI)[0]},{(cij*r_sqrtPMI)[1]},{(cij*r_sqrtPMI)[2]}, dTerm:{dTerm}') return residuals #------------------------------------------------------# #--- SECTION (3) Graphical User Interface for fit ---# #------------------------------------------------------# # Create user interface to handle isotopologue data/files, to select fit type # (r0,rs,rmI,rmII), fit options, etc. # Automatically activate / deavtivate fit parameters ci, di, dH as required. # Load/Save fit options, fit data, fit results. # ToDo: - Accept tabular data copy/paste from publication format # - Implement fixed parameters # (3.1) Utility functions #-------------------------------------------------------------------- # Convert between fitList and string representation of data # from calculation of fit-relevant parameters (fit weight, atom masses) def list2dataString(isoList, isoMasses={}, targetUnit='amu*AA', extended=False): """ Convert Isotopologue data (units of amu*AA): list --> string (using target units). """ dataString = f' Isotopologue, [ IA /{targetUnit:6}, IB /{targetUnit:6}, IC /{targetUnit:6}], [ uA /{targetUnit:6}, uB /{targetUnit:6}, uC /{targetUnit:6}], Fit[A,B,C], Author' if extended: dataString += ', Symmetry, AxisRot [z,y,x], abc<->xyz' dataString += '\n' + '-'*134 # Header line if extended: dataString += '-'*44 dataString += '\n' d = {'amu*AA':9,'MHz':6,'1/cm':11}[targetUnit] # Number of digits depending on unit for i,iso in enumerate(isoList): try: atoms,([A,B,C],[sA,sB,sC]),[wA,wB,wC],author,*rest = iso except: atoms,([A,B,C],[sA,sB,sC]),[wA,wB,wC],author = iso rest = [] if len(rest) > 2: sym,[rotZ,rotY,rotX],mapAxes = rest[:3] #print(f'sym:{sym}, rot:{rotZ,rotY,rotX}, mapAxes:{mapAxes}') else: sym,[rotZ,rotY,rotX] = 'unknown', [0.,0.,0.] if targetUnit != 'amu*AA': # Need to convert units sA,sB,sC = sA/A,sB/B,sC/C # To fractional errors A,B,C = [convertToUnit[targetUnit](val) for val in [A,B,C]] # To amu*A2 sA,sB,sC = [sA*A,sB*B,sC*C] # To errors in amu*A2 dataString += f'{i:2}: {atoms:12}, [{A:13.{d}f},{B:13.{d}f},{C:13.{d}f}], [{sA:12.{d}f},{sB:12.{d}f},{sC:12.{d}f}], ' fitA,fitB,fitC = [1 if w>0 else 0 for w in [wA,wB,wC]] # Fit flags dataString += f' [{fitA}, {fitB}, {fitC}], {author:>16}' if extended: dataString += f', {sym:>10}, [{rotZ:>5},{rotY:>4},{rotX:>4}], {mapAxes}' dataString += '\n' dataString += '-'*134 # Add masses at bottom of string if extended: dataString += '-'*44 dataString += '\nisotope masses ' + ', '.join([f'{key}:{value}' for key,value in isoMasses.items()]) +'\n' return dataString def dataString2list(dataString, weights=0): """ Convert Isotopologue data: string --> fit list (in units of amu*A**2) Fit weight==0,1,2,3 creates weights 1, I/uI, (I/uI)**2, sqrt(I/uI). """ if DEBUG: log.write('dataString2list') header,_,data = dataString.partition('---\n') # Header above, isotope data below line --- possibleUnits = convertToUnit.keys() units = [unit for unit in possibleUnits if unit in header] # Units specified in header line? if len(units) == 1: # Found unit string in header? unit = units[0] # --> Unit established if DEBUG: log.write(f'Unit established ({unit})') else: # No unit string found unit = 'amu*AA' # --> Assume default units print(u'Give one unit in header line {units}! (assuming amu*A\u00B2)') fitList = [] # Now create list of data isoMasses = {} for line in data.split('\n'): # Analyze each line of isotopologue data try: # Interpret line as isotope data line for char in ('[',']','(',')'): line = line.replace(char,'') # Remove brackets atoms,A,B,C,sA,sB,sC,fA,fB,fC,author,*rest = line.split(',') # Split isotopologue information fields A,B,C,sA,sB,sC = npa([A,B,C,sA,sB,sC]).astype(float) # Convert strings to floats try: # Handle possible symmetry data sym = npa(rest)[0] # Symmetry string rx,ry,rz = npa(rest)[1:4].astype(float) # Euler angles mapping principal axes to xyz mx,my,mz = npa(rest)[4:7].astype(int) # Mapping abc <-> xyz except: sym,[rx,ry,rz],[mx,my,mz] = '',[0,0,0],[0,0,0] # print('Exception: no symmetry information in Isotope data string.') if unit != 'amu*AA': # Need to convert units to amu*Angstrom**2 fsA,fsB,fsC = sA/A,sB/B,sC/C # Fractional errors A,B,C = [convertFromUnit[unit](val) for val in [A,B,C]] # ABC constants to amu*A2 sA,sB,sC = fsA*A,fsB*B,fsC*C # Errors in amu*A2 wABC = npa([int(fA), int(fB), int(fC)]) # Weights based on fit flags [0,1] if weights ==1: # Weights proposed by Eijck wABC = wABC/(npa([sA,sB,sC])**2 + npa([0.0005]*3)**2) # [J.Mol.Spec. 91,348(1982)] if weights > 1: try: wABC = wABC * npa([1/sA/A,1/sB/B,1/sC/C]) # Requested fit weight prop to 1/uncertainty except: # Catch divide by zero error if sA,sB,sC values are missing wABC = 0 print('Need non-zero rot. constant uncertainties to use weights.') if weights ==3: wABC= wABC**2 # Fit weights prop to 1/uncertainty**2 if weights ==4: wABC = np.sqrt(wABC) # Fit weights prop to 1/sqrt(uncertainty) wABC[~np.isfinite(wABC)] = 0 # Replace np.nan and np.inf wA,wB,wC = wABC if ':' in atoms: atoms = atoms.split(':')[1] fitList.append([atoms.strip(),([A,B,C],[sA,sB,sC]),[wA,wB,wC],author.strip(),sym.strip(),[rx,ry,rz]]) # Assemble fitList except: if 'mass' in line: # Read isotope masses from lines containing 'mass' if DEBUG: log.write(f"found masses in line: {line}") isoMasses = valueString2Dict(line,isoMasses) # Add masses to isoMasses dictionary elif '---' in line: pass # Separator lines contain '---' elif len(line)>10: print(f'Unrecognized data line: {line}') # Warn about unrecognized data line for iso in fitList: # Add missing masses from NISTmass for atom in iso[0].strip(): # Check is isoMasses contains all relevant masses if not atom in isoMasses.keys(): print(f"Taking mass for atom '{atom}' from NISTmass") try: isoMasses[atom] = NISTmass[atom] except: print(f"Can't find mass for atom '{atom}'") return fitList, isoMasses # Contract isotopologue list that conatins multiple entries for the same isotopologue # To a list with one entry and modified uncertainty for each isotopologue. (Shouldn't # affect fit results but will change Monte-Carlo parameter uncertainties.) def contractIsoList(isoList): """ Combine isotopologue list (amu*AA units) from different sources into single value, uncertainty pairs with the uncertainty reflecting the range of values. """ isoCount, isoConst, isoSigma, isoWeight = {},{},{},{} # Dicts of isotopologue species isoAuthor, isoRest = {},{} for i,iso in enumerate(isoList): # Create list of isotopologue species atoms,([A,B,C],[sA,sB,sC]),[wA,wB,wC],author,*rest = iso # Unpack current isotopologue information if not atoms in isoCount: # Add entry for new isotopologue isoCount[atoms] = 1 # ... isoConst[atoms] = [[A,B,C],] isoSigma[atoms] = [[sA,sB,sC],] isoWeight[atoms] = [[wA,wB,wC],] isoAuthor[atoms] = [author,] isoRest[atoms] = rest # Information on symmetry, atom masses else: # Add entry for existing isotopologue isoCount[atoms] += 1 # ... isoConst[atoms].append([A,B,C]) isoSigma[atoms].append([sA,sB,sC]) isoWeight[atoms].append([wA,wB,wC]) isoAuthor[atoms].append(author) fitList = [] # Collapsed isotopologue list for atoms in isoCount: # For each unique isotopologue if isoCount[atoms] == 1: # Single entry: Keep original values print(f'{atoms}: 1 entry') fitList.append([atoms,(isoConst[atoms][0],isoSigma[atoms][0]),isoWeight[atoms][0],isoAuthor[atoms][0],*isoRest[atoms]]) if isoCount[atoms] > 1: # Multiple entries: Calculate average, range, ... print(f'{atoms}: {isoCount[atoms]} entries') fitFlag = npa(isoWeight[atoms]) != 0 # Fit flags: don't use values where weight==0! ABC = npa(isoConst[atoms])*fitFlag # Rot constants for fit valsABC = np.count_nonzero(ABC,axis=0) # Number of constants for axis a,b,c weightABC = np.sum(isoWeight[atoms], axis=0) sigmaABC = [0,0,0] avgABC = [0,0,0] for i,(f,B) in enumerate(zip(fitFlag.T,ABC.T)): # Iterate over Axes a,b,c if np.count_nonzero(f) > 0: # Is there a fit constant for axis? avgABC[i] = np.average(B[f]) sigmaABC[i] = gSum(npa(isoSigma[atoms])[:,i][f]) # Propagated uncertainty if np.count_nonzero(f) >= 2: # Need to account for std of rot constants stdev = np.std(B[f]) sigmaABC[i] = gSum([stdev, sigmaABC[i]]) fitList.append([atoms,(avgABC,sigmaABC),list(weightABC),','.join(isoAuthor[atoms]),*isoRest[atoms]]) return fitList # Retrieve parameter values from a string (par1:val1 par2:val2 ...) def valueString2Dict(valueString, valueDict={}): """ Take sting containing space separated label value:pairs and add/replace values in valueDict. "A:1, B:2 C:3; other:4. Just text" --> {'A':1,'B':2,'C':3,'other':4} """ if len(valueString) == 0: return valueDict # Empty input terms = valueString.replace(',',' ').split() # User supplied value string (a:1, b:2, ...) for term in terms: try: label,value = term.split(':') # Get user supplied rm1, rm2 (c,d) correction terms while not value[-1].isdigit(): value = value[:-1] # Remove non-numerical chars if len(label) and len(value): valueDict[label] = float(value) # Add / uppdate value in valueDict except: pass # print(f'unpackValueString: Failed to recognize {term}') return valueDict # Utility function GREP: Return values after keywords. (For loading of data.) def grepLines(text, keywords): """ Return dictionary with {key:line} containing textlines after keywords. If a keyword is not found, line returns an empty string. Dictionary entry 'leftovers' contains the lines of string after the last keyword.""" returnDict = dict.fromkeys(keywords, '0') # Initialize return dictionary lastKey = 0 # Track location of last keyword in text for i,line in enumerate(text.splitlines()): # Parse text line by line for keyword in keywords: # Iterate over keywords if keyword in line: # Found keyword in line returnDict[keyword] = line.partition(keyword)[2].lstrip(' :') # Return string after 'keyword : ' lastKey = i break returnDict['lefttovers'] = '\n'.join(text.splitlines()[lastKey+1:]) return returnDict def grepFloats(text, keywords, skip=0): """ Return dictionary with {key:value} containing float values after keywords. If a keyword is not found, value np.NAN is returned. """ returnDict = dict.fromkeys(keywords, '0') # Initialize return dictionary for key in keywords: # Iterate over keywords try: value = text.split(key)[1][:100] # Look for float values in 100 characters after key pos = 0 # Start position to loook for float for s in range(skip+1): for i in range(pos,len(value)): # Parse string to cut leading non-numerical characters if value[i].isnumeric() or value[i] in ['-','.']: # Legal leading float characters #print(f'Key {key}: Leading char {i} is non-numeric') break for k in range(i,len(value)): # Cut trailing non-numerical characters if not (value[k].isnumeric() or value[k] in ['.','e','+','-']): # Legal float characters #print(f'Key {key}: Trailing char {k} is non-numeric') break returnDict[key] = float(value[i:k]) pos = k+1 except: returnDict[key] = np.NAN return returnDict # Calculate rotation matrix based on zyx rotation angles (deg) in 'rz:nn, ry:nn, rx:nn'. def makeRotMat(rotAxes): """ Return rotation matrix based on zyx rotation angles (deg) in 'rz:nn, ry:nn, rx:nn'. """ rot = valueString2Dict(rotAxes) rx = npa([[1,0,0],[0,1,0],[0,0,1]]) ry,rz = rx,rx if rot['x'] != 0: a = rot['x'] /180*np.pi rx = npa([[1,0,0],[0,np.cos(a),-np.sin(a)],[0,np.sin(a),np.cos(a)]]) if rot['y'] != 0: a = rot['y'] /180*np.pi ry = npa([[np.cos(a),np.sin(a),0],[0,1,0],[-np.sin(a),0,np.cos(a)]]) if rot['z'] != 0: a = rot['z'] /180*np.pi rz = npa([[np.cos(a),-np.sin(a),0],[np.sin(a),np.cos(a),0],[0,0,1],]) return np.dot(np.dot(rz,ry),rx), rot # Return rotation angles based on zyx rotation angles (deg) in 'rz:nn, ry:nn, rx:nn'. def rotAngles(): rot = valueString2Dict(GUI.customRot.text()) x,y,z = rot['x'],rot['y'],rot['z'] return x,y,z # Equivalent rotation matrices are obtained by: # R.as_matrix(R.from_euler('xyz', rotAngles(), degrees=True)) # makeRotMat(GUI.customRot.text()) # The second option is 3x faster. # Equivalent coordinate rotationis obtained by: # rot = R.from_euler('xyz', rotAngles(), degrees=True) # rot.apply(rc) # RM = makeRotMat(GUI.customRot.text()) # [(RM @ r).tolist() for r in rc] # The second option is 2x faster for rc containing 3 atom coordinates. # QTextEdit subclass with editingFinished signal, see: https://gist.github.com/hahastudio/4345418 class MyTextEdit(QtGui.QTextEdit): """ A TextEdit editor that sends editingFinished events when the text was changed and focus is lost. """ editingFinished = QtCore.pyqtSignal() receivedFocus = QtCore.pyqtSignal() def __init__(self, parent): super(MyTextEdit, self).__init__(parent) self._changed = False self.setTabChangesFocus( True ) self.textChanged.connect( self._handle_text_changed ) def focusInEvent(self, event): super(MyTextEdit, self).focusInEvent( event ) self.receivedFocus.emit() def focusOutEvent(self, event): if self._changed: self.editingFinished.emit() super(MyTextEdit, self).focusOutEvent( event ) def _handle_text_changed(self): self._changed = True def setTextChanged(self, state=True): self._changed = state def setHtml(self, html): QtGui.QTextEdit.setHtml(self, html) self._changed = False # (3.2) Function to print and plot Monte-Carlo analysis results #-------------------------------------------------------------------- # Plotting the MC results helps to spot correlations in the fit parameters. def plotMonteCarloResults(mcSim,pNames,name='Name'): """ Calculate and plot results for Monte Carlo Simulation. """ # Add normal distribution over histogram nSamples, nCoeff = mcSim.shape normalDist = lambda x,x0,sigma: (1 / (np.sqrt(2 * np.pi) * sigma)) *np.exp(-0.5 *((x-x0)/sigma)**2) parAVG = np.mean(mcSim,axis=0) parMED = np.median(mcSim,axis=0) parSTD = np.std(mcSim,axis=0) resultString = f"\nMonte-Carlo Simulation for Fit Function: {name}" resultString += f"\n------------------------------------------------\n" resultString += f"{'Parameter':10} {'Average':>9} {'Median':>9} (1-sigma)" resultString += f"\n------------------------------------------------\n" resultString += '\n'.join([f'{n:>10} = {a:>9.5f} {m:>9.5f} ({s:.5f})' for n,a,m,s in zip(pNames,parAVG,parMED,parSTD)]) resultString += f"\n------------------------------------------------\n" Pearson = np.corrcoef(mcSim.T) # Pearson's (linear) correlation coefficients resultString += f"\nPearson (Linear) Correlation Coefficient\n{'':6}" resultString += ''.join([f'{pName:>8}' for pName in pNames]) for i in range(len(pNames)-1): # Correlation coeff. (without chisqr) resultString += f'\n{pNames[i]:>8}' for j in range(i+1): resultString += f'{Pearson[i,j]:>8.3f}' resultString += f"\n-----------------------------------------------\n" Spearman = spearmanr(mcSim)[0] # Spearman's correlation coefficients resultString += f"\nSpearman Correlation Coefficient\n{'':6}" resultString += ''.join([f'{pNames[i]:>8}' for i in range(len(pNames))]) try: for i in range(len(pNames)): resultString += f'\n{pNames[i]:>8}' for j in range(i+1): resultString += f'{Spearman[i,j]:>8.3f}' except: pass #--- Plot results in default PyQtGraph correlation window. dtype,fields = [],[] # Prepare data for plotting for name in pNames: dtype.append((name,float)) if name[0] == 'r': # Assume distance parameter name starts with 'r' fields.append((name,{'units':u'\u212B'})) elif name[0] == 'a' or name == 'z-Rot': # Assume angle parameter name starts with 'a' fields.append((name,{'units':u'\u00B0 degrees'})) elif name == 'dH': # Assume Laurie parameter name is 'dH' fields.append((name,{'units':u'\u212B\u22C5amu0.5'})) else: fields.append((name, {'units':'arb. units'})) data = np.empty(nSamples, dtype=dtype) for i,name in enumerate(pNames): data[name] = mcSim.T[i] pg.mkQApp() # Create PyQtGraph Scatter Plot Widget in QT App spw = pg.ScatterPlotWidget() spw.setToolTip('Select parameter pairs to graph the parameter correlation.') spw.setFields(fields) spw.style['symbolPen']=pg.mkPen('b') spw.plot.showGrid(x=True, y=True, alpha=0.25) spw.setData(data) spw.setWindowTitle(name) font = QtGui.QFont('Arial', 11) xaxis = spw.plot.getAxis(['bottom'][0]) xaxis.setTickFont(font) xaxis.setHeight(int(11 * 1.7 + 22)) xaxis.setStyle(tickTextOffset=7) xaxis.label.setFont(font) yaxis = spw.plot.getAxis(['left'][0]) yaxis.setTickFont(font) yaxis.label.setFont(font) spw.show() return resultString, spw # (3.3) List fit functions, residual functions for the UI. #-------------------------------------------------------------------- # The following dict will populate the fit options in the fitUI. ResidualsList = {'Fit r0 geometry':Residuals_r0, #'Fit rs (Typke, Untested!)':Residuals_rs_Typke, 'Fit rs (Nossberger, I-I0)':Residuals_rs_Nossberger, 'Fit rs geometry':Residuals_rs,'Fit rsr geometry':Residuals_rsr, 'Fit rm(1) geometry':Residuals_rmI, 'Fit rm(1)r geometry':Residuals_rmIr, 'Fit rm(2) geometry':Residuals_rmII, 'Fit rm(2)r geometry':Residuals_rmIIr} # The following list will populate the interanel-> cartesian options in the fitUI. int2cartList = [bz_internal2cart,bz_internal2cart_rC_rH_rD_Laurie,h2o_internal2cart_sym, furan_int2cart,furan_int2cart_noH,furan_int2cart_scaledH,furan_int2cart_xy, h2o_internal2cart_rHrD_sym,h2o_internal2cart,O3_int2cart_sym,O3_int2cart, HCN_int2cart,HNCO_int2cart] # (3.4) Graphical User Interface. #-------------------------------------------------------------------- # I try to keep things simple, data is entered in string fields and then interpreted. # The UI uses widgets from PyQtgraph and PyQt. class fitUI(QtGui.QMainWindow): """ Simple GUI for fitting molecular structure from rovibrational data. I store some variables as parameter classes (self.xyz) to allow easy debugging :). """ def __init__(self): """ Create application, add widgets, and show GUI. """ QtGui.QMainWindow.__init__(self) self.setGeometry(80, 28, 1340, 500) self.setWindowTitle('Fit Molecular Structure from Rotational Constants') # QtGui.QToolTip.setFont(QtGui.QFont('Lucida Console', 10)) centralwidget = QtGui.QWidget() self.setCentralWidget(centralwidget) self.layout = QtGui.QGridLayout() centralwidget.setLayout(self.layout) self.addMenus() # This creates the GUI elements self.addWidgets() # This creates the GUI elements self.initConstants() # Initialize constants, data-lists self.show() self.dir = _dir def addMenus(self): #--- File menu ---# self.Menu = self.menuBar() fileMenu = self.Menu.addMenu('&File') menuitems = ["&Load Fit","&Save Fit","&Quit"] # Lazy way to add multiple menu items for item in menuitems: # Check FileMenu() for action implementations entry = fileMenu.addAction(item) entry.triggered.connect(partial(self.FileMenu,item=item)) def FileMenu(self, item): """ This implements functions for the FileMenu.""" if item == "&Load Fit": print("Loading ...") fname,_ = QtGui.QFileDialog.getOpenFileName(self, 'Load File', self.dir,'Fit Files(*.fit)') self.setWindowTitle(f"Fit Molecular Structure --- {fname.rsplit('/',1)[1]}") file = open(fname, 'r') dataRead = file.read() file.close() #i2c,r_internal,symmetry,fitType,Laurie,roVib,fitWeight,_,*data = dataRead.splitlines() keywords = 'int2cart,r_internal,symmetry,fitType,Laurie Correction,roVib parameters,Drop parameters,fitWeight'.split(',') vals = grepLines(dataRead, keywords) self.i2c.setCurrentIndex(self.i2c.findText(vals['int2cart'][4:])) self.r_internal.setText(vals['r_internal']) self.symmetry.setCurrentIndex(int(vals['symmetry'][0])) self.fitType.setCurrentIndex(self.fitType.findText(vals['fitType'][4:])) self.laurieBox.setChecked('True' in vals['Laurie Correction']) self.roVib.setText(vals['roVib parameters']) self.dropParams.setCurrentIndex(int(vals['Drop parameters'][0])) self.fitWeight.setCurrentIndex(int(vals['fitWeight'][0])) self.isoData.setText(vals['lefttovers']) self.updateData() if item == "&Save Fit": fname,_ = QtGui.QFileDialog.getSaveFileName(self, 'Save File', self.dir,'Fit Files(*.fit)') int2cart = f'int2cart: {self.i2c.currentIndex()} : {self.i2c.currentText()}\n' r_internal = f'r_internal:: {self.r_internal.text()}\n' symmetry = f'symmetry: {self.symmetry.currentIndex()} : {self.symmetry.currentText()}\n' fitType = f'fitType: {self.fitType.currentIndex()} : {self.fitType.currentText()}\n' Laurie = f'Laurie Correction: {self.laurieBox.isChecked()}\n' roVib = f'roVib parameters:: {self.roVib.text()}\n' dropParams = f'Drop parameters: {self.dropParams.currentIndex()} : {self.dropParams.currentText()}\n' fitWeight = f'fitWeight: {self.fitWeight.currentIndex()} : {self.fitWeight.currentText()}\n' data = f'\n{self.isoData.toPlainText()}' file = open(fname,'w') file.write(int2cart+ r_internal+ symmetry+ fitType+ Laurie+ roVib+ dropParams+ fitWeight+ data) file.close() if item == "&Quit": self.close() def addWidgets(self): """ Add GUI buttons, fields, etc.. """ #label = QtWidgets.QLabel('Fit Properties') #label.setAlignment(QtCore.Qt.AlignCenter | QtCore.Qt.AlignVCenter) #label.setFont(QtGui.QFont('Arial', 14)) #self.layout.addWidget(label, 1,1,1,1) self.fitType = QtGui.QComboBox() #--- Select Fit Type tt = 'This specifies the rovibrational corrections:\n'\ 'r0 denotes direct fit without corrections.\n'\ 'rs denotes Nossberger-type fit of substitution satructure.\n'\ 'rm(1,2) denotes first and second order mass-weighted corrections.\n'\ '..r denotes that corrections are applied in the invariant (rotated) molecular axis frame.' self.fitType.setToolTip(tt) for key,value in ResidualsList.items(): self.fitType.addItem(key,value) self.fitType.currentIndexChanged.connect(self.selectFitType) self.layout.addWidget(self.fitType, 1,1,1,1) self.laurieBox = QtGui.QCheckBox("Laurie Correction") #--- Select Laurie Correction self.laurieBox.setToolTip('This adds a localized correction for rH,rD bond lengths.') self.laurieBox.stateChanged.connect(self.selectFitType) self.layout.addWidget(self.laurieBox, 1,2,1,1) self.contractBox = QtGui.QCheckBox("Contract Datalist") #--- Select Laurie Correction self.contractBox.setToolTip('Calculate single average/uncertainty value for each isotope species before fitting (faster execution).') self.layout.addWidget(self.contractBox, 1,3,1,1) self.symmetry = QtGui.QComboBox() #--- Select molecular symmetry self.symmetry.setToolTip('Select the molecular symmetry to determine number of rovib. correction parameters') self.symmetry.addItem('Analyze Int2Cart symmetry') self.symmetry.addItem('Asymmetric Molecule') self.symmetry.addItem('Planar Molecule') self.symmetry.addItem('Prolate Top') self.symmetry.addItem('Oblate Top') self.symmetry.addItem('Planar Oblate Top') self.symmetry.addItem('Linear Molecule') self.symmetry.addItem('Spherical Top') self.symmetry.setItemData(0, QtGui.QFont("myFontFamily",italic=True), QtCore.Qt.FontRole) self.symmetry.currentIndexChanged.connect(self.selectFitType) self.layout.addWidget(self.symmetry, 1,4,1,1) self.i2c = QtGui.QComboBox() #--- Select int2cart function self.i2c.setToolTip('Select the mapping function for internal -> Cartesian coordinates') for func in int2cartList: self.i2c.addItem(repr(func).split()[1],func) # Names and pointers to int2cart functions self.layout.addWidget(self.i2c, 2,4,1,1) label = QtWidgets.QLabel('Internal Coordinates') #--- Enter Coordinate parameters label.setAlignment(QtCore.Qt.AlignRight | QtCore.Qt.AlignVCenter) self.layout.addWidget(label, 2,0,1,1) self.r_internal = QtGui.QLineEdit('r1:1, r2:1.5, r3:1.5, a1:120, a2:120 (geometry parameters for Internal -> Cartesian conversion)') self.r_internal.setToolTip('Internal structure parameters ri (in correct order for Internal -> Cartesian conversion)') self.r_internal.setDragEnabled(True) self.layout.addWidget(self.r_internal, 2,1,1,2) self.showMolecule= QtGui.QPushButton("Show Molecule") #--- Fit button self.showMolecule.setToolTip('Show the molecular structure based on int2Cart and Internal Coordinates.') self.showMolecule.clicked.connect(self._showMolecule) self.showMolecule.setStyleSheet("background-color: yellow") self.layout.addWidget(self.showMolecule,2,3,1,1) label = QtWidgets.QLabel('Rovibrational Correction Parameters') #--- Enter rovib. correction pars. label.setAlignment(QtCore.Qt.AlignRight | QtCore.Qt.AlignVCenter) self.layout.addWidget(label, 3,0,1,1) self.roVib = QtGui.QLineEdit('ca:value, cb:value, ... (as required for selected fit method)') self.roVib.setToolTip('rovibrational correction terms ci, di, as required for Fit Type and molecular symmetry.') self.roVib.setDragEnabled(True) #self.roVib.editingFinished.connect(self.updateRovibParameters) self.layout.addWidget(self.roVib, 3,1,1,3) self.dropParams = QtGui.QComboBox() #--- Reduce number of rovib. correction pars. self.dropParams.setToolTip('Reduce number of rovib. correction parameters.\nUse only diagonal parameters for rs,rm fits.\nUse full parameter set for rsr,rmr fits. (see Watson1999)') self.dropParams.addItem('Single c,d parameter') self.dropParams.addItem('Use only diagonal ci,di parameters') self.dropParams.addItem('Full parameter set (up to 6 cij, 3 di parameters)') self.dropParams.setCurrentIndex(2) self.dropParams.currentIndexChanged.connect(self.selectFitType) self.layout.addWidget(self.dropParams, 3,4,1,1) self.fitWeight = QtGui.QComboBox() #--- Select weighting of fitted data self.fitWeight.setToolTip('Choose how experimental data is weigthed for the fit.') self.fitWeight.addItem('Fit Weight 0,1') self.fitWeight.addItem('Fit Weight 1/(uI**2 + 0.0005uA**2') self.fitWeight.addItem('Fit Weight I/uI') self.fitWeight.addItem('Fit Weight (I/uI)**2') self.fitWeight.addItem('Fit Weight (I/uI)*0.5') self.fitWeight.setCurrentIndex(1) self.layout.addWidget(self.fitWeight,4,0,1,1) self.units = QtGui.QComboBox() #--- Select display unit for data self.units.setToolTip('Convert units of experimental data.\nCorrect unit string must appear in the data header (first line)!') for unit in convertToUnit.keys(): # See unit conversion functions, above self.units.addItem(unit) self.layout.addWidget(self.units,4,1,1,1) self.units.currentIndexChanged.connect(self.updateData) # Update data if units are changed self.customRot = QtGui.QLineEdit('z:0, y:0, x:0 (Reference frame rotation in deg.)') self.customRot.setToolTip('rovibrational correction terms are applied along the reference coordinate axes.\nuse z: zstart...zend for MC simulation with z angle distribution.') self.customRot.setDragEnabled(True) self.customRot.editingFinished.connect(self.analyzeSymmetry) self.layout.addWidget(self.customRot, 4,2,1,1) self.parent = QtGui.QComboBox() #--- Select parent molecule (defines PAS) self.parent.setToolTip("Select the parent molecule (defines principal axis system for '...r' fits).") self.parent.addItem('Enter isotopologue data') self.parent.currentIndexChanged.connect(self.updateParent) self.layout.addWidget(self.parent, 4,4,1,1) self.fitButton = QtGui.QPushButton("Run the Fit") #--- Fit button self.fitButton.setToolTip('Run the selected fit with durrent data.') self.fitButton.clicked.connect(self.Fit) self.fitButton.setStyleSheet("background-color: green") self.layout.addWidget(self.fitButton,1,5,1,1) self.MCButton = QtGui.QPushButton("Run MC sim.") self.MCButton.setToolTip('Run a Monte-Carlo simulation with selected settings.') self.MCButton.clicked.connect(self.McSim) self.MCButton.setStyleSheet("background-color: green") self.layout.addWidget(self.MCButton,2,5,1,1) self.pbar = QtWidgets.QProgressBar(self) # Progress bar for Monte-Carlo Simulation self.pbar.setFixedWidth(120) self.layout.addWidget(self.pbar, 3,5,1,1) self.pbar.setStyleSheet('text-align: center') self.extraInfo = QtGui.QCheckBox("Show extra info") #--- Toggle additional information in isoData window self.extraInfo.setToolTip('This adds a localized correction for rH,rD bond lengths.') self.extraInfo.setChecked(True) self.extraInfo.stateChanged.connect(self.analyzeSymmetry) self.layout.addWidget(self.extraInfo, 4,5,1,1) tt = 'Text describing isotopologue atomic composition, measured inertial moments, i.m. uncertainties, and fit flags.\n'\ 'Text above ---\\n is considered as header and must contain the unit for experimental data (MHz, amu*AA, 1/cm).\n'\ 'Experimental data lines contain comma-separated values for:'\ ' - atomic composition (single-letter code for each atoms, corresponding to the int2cart function; e.g.: HOH)\n'\ ' - inertial moments [Ia,Ib,Ic], and uncertainties [sIa,sIb,sIc]\n'\ ' - fit weights [wa,wb,wc] for axis a,b,c; set to zero to drop a value from the fit\n'\ ' - author string to denote the origin of the data\n'\ 'Mass information for each atom in the atomic composition must be given as atom:mass in a line containing the keyword "mass",\n'\ 'e.g., mass: H:1.00782503, D:2.0141022, O:15.9949146221\n' self.dataHeader = "Isotopologue, [ IA (MHz), IB (MHz), IC (MHz)], " self.dataHeader += "[uA (MHz), uB (MHz), uC (MHz)], Fit[A,B,C], Author\n\r" self.dataHeader += '-'*110 + '\n' self.isoData = MyTextEdit(self.dataHeader) self.isoData.setFont(QtGui.QFont('Lucida Console', 9)) self.isoData.setToolTip(tt) self.isoData.setLineWrapColumnOrWidth(1600) self.isoData.setLineWrapMode(QtGui.QTextEdit.FixedPixelWidth) self.isoData.editingFinished.connect(self.updateData) self.layout.addWidget(self.isoData, 5,0,-1,-1) self.cursor = self.isoData.textCursor() #self.dataSym = QtGui.QTextEdit("symmetry, rot axis, angle") #self.dataSym.setFont(QtGui.QFont('Lucida Console', 9)) #self.dataSym.setToolTip("""Symmetry of isotopologues and their rotation relative to the parent (first) i.""") #self.layout.addWidget(self.dataSym, 5,4,-1,1) def _showMolecule(self): int2cart = self.i2c.currentData() # User selected int2cart function fitParameters = valueString2Dict(self.r_internal.text(),{}) # Internal geometry parameter list rCart = int2cart(list(fitParameters.values())) # Cartesian atom coordinates parent = self.parent.currentIndex() # Parent molecule isoList,isoMass = dataString2list(self.isoData.toPlainText(), self.fitWeight.currentIndex()) # Create fitList from current data. atoms = isoList[parent][0] # Parent molecule string (e.g., 'CCCCOHHHH') self.plot = plotMolecule([moleculeString, rCart]) # Show molecule def initConstants(self): """ Initialize fitList and default rovibrational correction parameters. This is to add some persistence for user supplied parameters. """ self.rovibParams = {'ca':0.001, 'cb':0.001, 'cc':0.001, # Diagonal rm1 correction terms 'cab':0.001, 'cac':0.001, 'cbc':0.001, # off-diagonal rm1 correction terms 'da':0.002, 'db':0.002, 'dc':0.002, # rm2 correction terms 'dH':0.01} # Laurie correction self.nPars = {'int':0, 'ci':0, 'di':0, 'Laurie':0} # [Number of internal, rmI, rm2, Laurie parameters] self.contract = False def selectFitType(self): """ Update rovib correction parameter list for selected fit type. """ if 'Analyze' in self.symmetry.currentText(): # Analyze symmetry properties self.analyzeSymmetry(update=True) #pass #self.updateParent() # Update parent molecule selection list (implicit analyzeSymmetry) #self.analyzeSymmetry() fitType = self.fitType.currentText() # User selected fit type (r0, rs, rm1, rm2) terms = self.roVib.text().split(' ') # User supplied rovibrational correction terms self.nPars['ci'] = 0 # Reset number of rmI fit parameters (ci) self.nPars['di'] = 0 # Reset number of rm2 fit parameters (di) self.nPars['Laurie'] = 0 # Reset Laurie fit parameter self.nPars['sym'] = self.symmetry.currentText() # Set symmetry parameter if ('only diagonal' in self.dropParams.currentText() # Off-diagonal cij required? or not ('rsr' in fitType or 'rm(1)r' in fitType or 'rm(2)r' in fitType)): OD = 0 # Fit only diagonal correction cii else: OD = 1 # Fit with off-diagonal corection cij # Initialize rovibrational constants for symmetry cases if 'r0' in fitType or 'Typke' in fitType or 'Nossberger' in fitType: # r0 fit: no rovib correction terms self.roVib.setStyleSheet("color: grey;") # r0 fit: grey out rovib. parameters roVibText = 'No correction terms required for r0 fit.' if self.laurieBox.isChecked(): roVibText = '' else: # Need rovib. correction terms for rs, rm1, rm2 fits self.roVib.setStyleSheet("color: black;") # Rovib corrections required: reset color if ('Linear' in self.symmetry.currentText() or 'Spherical' in self.symmetry.currentText() or 'Single' in self.dropParams.currentText()): roVibText = f"cb:{self.rovibParams['cb']}" # Only 1 rovib. parameter self.nPars['ci'] = 1; # Update number of c parameters if 'rm(2)' in fitType: # rm2 correction d-Terms required roVibText += f", db:{self.rovibParams['db']}" self.nPars['di'] = 1; # Update number of d parameters elif 'Planar Oblate' in self.symmetry.currentText(): # Ia=Ib; planar: no cac,cbc; oblate: ca=cb (=cab?), roVibText = f"ca:{self.rovibParams['ca']}, cc:{self.rovibParams['cc']}" self.nPars['ci'] = 2; # Update number of c parameters if 'rm(2)' in fitType: # rm2 correction d-Terms required roVibText += f", da:{self.rovibParams['da']}, dc:{self.rovibParams['dc']}" self.nPars['di'] = 2; # Update number of d parameters elif 'Oblate Top' in self.symmetry.currentText(): # Ia=Ib roVibText = f"ca:{self.rovibParams['ca']}, cc:{self.rovibParams['cc']}" self.nPars['ci'] = 2; # Update number of c parameters if OD: # off-diagonal cij elements required roVibText += f", cab:{self.rovibParams['cab']}, cac:{self.rovibParams['cac']}, cbc:{self.rovibParams['cbc']}" self.nPars['ci'] = 5 # Update number of c parameters if 'rm(2)' in fitType: # rm2 correction d-Terms required roVibText += f", da:{self.rovibParams['da']}, dc:{self.rovibParams['dc']}" self.nPars['di'] = 2 # Update number of d parameters elif 'Prolate Top' in self.symmetry.currentText(): # Ib=Ic roVibText = f"ca:{self.rovibParams['ca']}, cb:{self.rovibParams['cb']}" self.nPars['ci'] = 2; # Update number of c parameters if OD: # rm2 correction d-Terms required roVibText += f", cab:{self.rovibParams['cab']}, cac:{self.rovibParams['cac']}, cbc:{self.rovibParams['cbc']}" self.nPars['ci'] = 5 # Update number of c parameters if 'rm(2)' in fitType: # rm2 correction d-Terms required roVibText += f", da:{self.rovibParams['da']}, db:{self.rovibParams['db']}" self.nPars['di'] = 2 # Update number of d parameters elif 'Planar Molecule' in self.symmetry.currentText(): # planar: no cac,cbc roVibText = f"ca:{self.rovibParams['ca']}, cb:{self.rovibParams['cb']}, cc:{self.rovibParams['cc']}" self.nPars['ci'] = 3; # Update number of c parameters if OD: # off-diagonal cij elements required roVibText += f", cab:{self.rovibParams['cab']}" self.nPars['ci'] = 4 # Update number of c parameters if 'rm(2)' in fitType: roVibText += f", da:{self.rovibParams['da']}, db:{self.rovibParams['db']}, dc:{self.rovibParams['dc']}" self.nPars['di'] = 3 # Update number of d parameters else: # Non-symmetric molecule: full parameter set roVibText = f"ca:{self.rovibParams['ca']}, cb:{self.rovibParams['cb']}, cc:{self.rovibParams['cc']}" self.nPars['ci'] = 3; # Update number of c parameters if OD: # off-diagonal cij elements required roVibText += f", cab:{self.rovibParams['cab']}, cac:{self.rovibParams['cac']}, cbc:{self.rovibParams['cbc']}" self.nPars['ci'] = 6 # Update number of c parameters if 'rm(2)' in fitType: roVibText += f", da:{self.rovibParams['da']}, db:{self.rovibParams['db']}, dc:{self.rovibParams['dc']}" self.nPars['di'] = 3 # Update number of d parameters if self.laurieBox.isChecked(): # Laurie correction selected self.roVib.setStyleSheet("color: black;") # Make sure text is not greyed out roVibText += f", dH:{self.rovibParams['dH']}" # Addy Laurie correction factor in correction terms self.nPars['Laurie'] = 1 # Update number of c parameters self.roVib.setText(roVibText) # Update parameters for selected fit type def Fit(self): fitInput = self.analyzeSymmetry(update=False) # fitList from isoList (masses!) self._FitStr,self._Fit = fit(*fitInput) # Perform fit self.showResults(self._FitStr,f"Results for {self.fitType.currentText().split()[1]}{'L' if self.laurieBox.isChecked() else ''}_fit_using_{self.i2c.currentText()}") def McSim(self): """ Run a Monte-Carlo simulation to test statistical bounds based on exp. uncertainties. """ samples,ok = QtGui.QInputDialog.getInt(self,"Number of samples", "Input number of samples for Monte Carlo simulation\nWarning: Large sampling is time-consuming.",100) # Number of Monte Carlo samples if not ok: return self.pbar.setValue(0) # Progress bar reset QtCore.QCoreApplication.processEvents() # Update screen # Fit settings and parameter start values. fitList, residualFunc, int2cart, fitParameters, nPars = self.analyzeSymmetry(update=False) # fitList from isoList (masses!) if self.contractBox.isChecked(): fitList = contractIsoList(fitList) # Reduce isoList to one entry per isoptopologue startValues = list(fitParameters.values()) # Initial fit parameter values Fit = leastsq(residualFunc, startValues, args=(fitList,nPars,int2cart), full_output=1) # Fit fitTitle = f'{repr(residualFunc).split()[1]}{"L" if nPars["Laurie"]==1 else ""}, {repr(int2cart).split()[1]}' printFitResults(fitList, fitParameters, Fit, title=fitTitle) optValues = Fit[0] # Optimized fit parameters as MC start values # Generate statistical sample of input values for MC simulation rot = False if '..' in self.customRot.text(): # Random reference frame rotation (z-axis) for rovib correction axes rot = True xyz = valueString2Dict(self.customRot.text()) x,y,z = xyz['x'],xyz['y'],xyz['z'], start,end = self.customRot.text().split('..') if 'z' in start[-5:]: # For now, allow z-axis variation zStart = float(start.split(':')[1]) # Minimum z-angle for Monte Carlo Simulation zEnd = float(end.split(',')[0]) # Max. z-angle for MCsim mcSim = [] # List to hold Monte-Carlo results for i in range(samples): # Run fits for MC sampling fitList_ = [] # List of isotopologues to fit for isotope in fitList: # Generate randomized input values within uncertainties for each isotopologue iso = [ isotope[0],([rng.normal(isotope[1][0][i],isotope[1][1][i]) for i in [0,1,2]], isotope[1][1])] iso.extend(isotope[2:]) fitList_.append(iso) if rot: # Rotation of reference frame for rovib. corrections z = np.random.uniform(zStart,zEnd) rotMat,rot = makeRotMat(f'z:{z}, y:{y}, x:{x}') nPars['rotMat'] = rotMat (params, pcov, infodict, errmsg, ok) = leastsq(residualFunc, optValues, args=(fitList_,nPars,int2cart), full_output=1) # Perform fit with current values if ok: # Successful fit if rot: params = np.append(params,z) # Add reference frame rotation parameter to list residuals = infodict['fvec'] # Residuals weights = [i for iso in fitList for i in iso[2] if i>0] scaledResiduals = residuals/weights # Remove weight scaling of residuals scaledChiSqr = (scaledResiduals**2).sum() # Chi squared for unweighted residuals params = np.append(params,scaledChiSqr) # Add unweighted residuals parameter to list mcSim.append(list(params)) # Store parameter list in mcSim self.pbar.setValue(int(i/(samples-1)*100)) # Progress bar update QtCore.QCoreApplication.processEvents() # Progress bar update self.mcSim = np.array(mcSim) # Convert results to numpy array title = f"Results for {self.fitType.currentText().split()[1]}{'L' if self.laurieBox.isChecked() else ''}_fit_using_{self.i2c.currentText()}" parNames = list(fitParameters.keys()) if rot: parNames = np.append(parNames,'z-Rot') parNames = np.append(parNames,'ChiSqr') mcResult,self.graph = plotMonteCarloResults(self.mcSim,parNames,name=title) self.showResults(mcResult,title) def showResults(self, resultString, title): print(resultString) self.resultWin = QtGui.QMessageBox() # ToDo: Replace with Qdialog for no-blocking popup (https://www.tutorialspoint.com/pyqt/pyqt_qdialog_class.htm) self.resultWin.setWindowTitle(f"Results for {self.r_internal.text().split()[1]}{'L' if self.laurieBox.isChecked() else ''} fit using {self.i2c.currentText()}") self.resultWin.setFont(QtGui.QFont('Lucida Console', 9)) # Add bond lengths and angles to results string. parent = self.parent.currentIndex() # Parent molecule isoList,isoMass = dataString2list(self.isoData.toPlainText(), self.fitWeight.currentIndex()) # Create fitList from current data. atoms = isoList[parent][0] # Parent molecule string (e.g., 'CCCCOHHHH') int2cart = self.i2c.currentData() # User selected int2cart function coord = int2cart(GUI._Fit[0]) resultString += '\n\n' + tabulateStructure(atoms, coord) self.resultWin.setDetailedText(resultString) self.resultWin.setStandardButtons(QtGui.QMessageBox.Ok | QtGui.QMessageBox.Save) for button in self.resultWin.buttons(): # Expand detail view of message box if self.resultWin.buttonRole(button) == QtGui.QMessageBox.ActionRole: button.click() break textWin = self.resultWin.findChildren(QtGui.QTextEdit)[0] textWin.setFixedSize(800,500) retval = self.resultWin.exec_() if retval == 2048: # Save to file name = self.dir + '\\' + title fname,_ = QtGui.QFileDialog.getSaveFileName(self, 'Save File', name,'Fit Files(*.fitresult)') with open(fname,'w') as file: file.write(resultString) def analyzeSymmetry(self, update=True): """ Analyze symmetry properties of molecules and assemble parameters required for a fit. If requested, update isotopologue text field. """ residualFunc = self.fitType.currentData() # Selected residuals function for fit fitWeight = self.fitWeight.currentIndex() # User selected weights for fit data int2cart = self.i2c.currentData() # User selected int2cart function fitParameters = valueString2Dict(self.r_internal.text(),{}) # Create internal geometry parameter list self.nPars['int'] = len(fitParameters) # Update number of internal parameters self.nPars['parent'] = self.parent.currentIndex() # Parent molecule (determining PAS axis system) rotMat,rot = makeRotMat(self.customRot.text()) # User selected PAS rotation (affecting rovib corrections) if rot['z']>0: print(f"custom {rot['z']} deg. z-axis rotation for rovib. correction.") self.nPars['rotMat'] = rotMat fitParameters = valueString2Dict(self.roVib.text(), fitParameters) # Add roVib parameters to list self.isoList,isoMass = dataString2list(self.isoData.toPlainText(), fitWeight) # Create fitList from current data. fitList = analyzeSymmetry(self.isoList,isoMass,int2cart,fitParameters,self.nPars) # Analysis of symmetry if update: # Update symmetry in symmetry Combobox self.isoData.editingFinished.disconnect(self.updateData) self.isoData.setText(list2dataString(fitList,isoMass,self.units.currentText(), self.extraInfo.isChecked())) # Convert back to string and show in isoData field self.isoData.editingFinished.connect(self.updateData) symList = npa(fitList, dtype=object)[:,4] # List of molecular symmetries self.symmetry.currentIndexChanged.disconnect() # Set correct symmetry option if 'linear' in symList: self.symmetry.setCurrentText('Linear Molecule') elif 'spherical' in symList: self.symmetry.setCurrentText('Spherical Top') elif 'oblate pl.' in symList: self.symmetry.setCurrentText('Planar Oblate Top') elif 'oblate' in symList: self.symmetry.setCurrentText('Oblate Top') elif 'prolate' in symList: self.symmetry.setCurrentText('Prolate Top') elif 'planar' in symList: self.symmetry.setCurrentText('Planar Molecule') else: self.symmetry.setCurrentText('Asymmetric Molecule') self.symmetry.currentIndexChanged.connect(self.selectFitType) return fitList, residualFunc, int2cart, fitParameters, self.nPars def updateData(self): """ Update data from isoData field. """ self.isoList,isoMass = dataString2list(self.isoData.toPlainText()) # Convert to list, fill information (e.g., NISTmass) #fitParameters = valueString2Dict(self.r_internal.text(),{}) # Create internal geometry parameter list #fitParameters = valueString2Dict(self.roVib.text(), fitParameters) # Add roVib parameters to list #self.nPars['parent'] = self.parent.currentIndex() # Parent molecule (determining PAS axis system) #self.nPars['rotMat'] = makeRotMat(self.customRot.text())[0] # User selected PAS rotation (affecting rovib corrections) #fitList = analyzeSymmetry(self.isoList,isoMass,self.i2c.currentData(),fitParameters,self.nPars) # Analysis of symmetry self.updateParent() self.analyzeSymmetry(update=True) def updateParent(self): try: self.parent.currentIndexChanged.disconnect() # Disconnect Combobox to allow changes currentParent = self.parent.currentIndex() # Remember parent species except: currentParent = 0 #if currentParent == -1: currentParent = 1 # Initial parent species is first in list self.parent.clear() # Reset combobox list for i,iso in enumerate(self.isoList): # Refill list with marked parent species if i == currentParent: self.parent.addItem(f'{i}: {iso[0]} (current parent molecule)') else: self.parent.addItem(f'{i}: {iso[0]}') self.parent.setCurrentIndex(currentParent) # Reset selection self.parent.currentIndexChanged.connect(self.updateParent) # Reconnect: update upon user selection self.analyzeSymmetry(update=True) def setData(self, isoList, massList={}): # Update isotopologue data from list """ Set isotopologue data from command line. """ self.isoData.setText(list2dataString(isoList,massList,self.units.currentText())) # Show data in isoData field self.updateData() def isotopeRM(self, isoNumber): # Return coordinates and atom masses for isotope (testing) """ Return Cartesians r and masses m for isotope iso. """ fitWeight = self.fitWeight.currentIndex() # User selected weights for fit data int2cart = self.i2c.currentData() # User selected int2cart function fitParameters = valueString2Dict(self.r_internal.text(),{}) # Create internal geometry parameter list isoList,isoMass = dataString2list(self.isoData.toPlainText(), fitWeight) # Create fitList from current data. iso = isoList[isoNumber] # Reference molecule (first in isoList) r = list(fitParameters.values()) # internal geometry parameters rc = int2cart(r) # Cartesian coordinates m = [isoMass[atom] for atom in iso[0]] # Atom masses return rc,m def UR(self): # Update residual functions (testing) """ Utility function to update residuals functions. """ self.fitType.clear() ResidualsList = {'Fit r0 geometry':Residuals_r0, 'Fit rs (Nossberger, I-I0)':Residuals_rs_Nossberger, 'Fit rs geometry':Residuals_rs,'Fit rsr geometry':Residuals_rsr, 'Fit rm(1) geometry':Residuals_rmI, 'Fit rm(1)r geometry':Residuals_rmIr, 'Fit rm(2) geometry':Residuals_rmII, 'Fit rm(2)r geometry':Residuals_rmIIr} for key,value in ResidualsList.items(): self.fitType.addItem(key,value) """ #----------------------------------------------# #---- Plotting of covariance distribution ----# #----------------------------------------------# # This code creates 2-D correlation plots for Monte-Carlo simulation results. # I didn't implement this code into a function because it requires quite specific # selection of plot parameters and correlations. I'll work this into the GUI at some point. #--- Run Monte-Carlo simulation and evaluate results ---# GUI = fitUI() # Execute graphical user interface #! --> load/enter data, set up desired fit (here: rm(2) and rm(2)L fitting of benzene data) #! --> run Monte-Carlo Simulation for plot. #--- Functions for visualizing covariance distributions ---# def density_estimation(xi, yi, X=None,Y=None): ''' Estimate density of points xi,yi on grid x,y (for plotting contour lines). ''' from scipy import stats if not X.any() or not X.any(): # No X,Y grid xmin,xmax = xi.min(), xi.max() # Default grid given by data range ymin,ymax = yi.min(),yi.max() # ... X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] # Create grid positions = np.vstack([X.ravel(), Y.ravel()]) values = np.vstack([xi, yi]) kernel = stats.gaussian_kde(values) # Calculate kernel density estimate Z = np.reshape(kernel(positions).T, X.shape) # Map kde on grid return X, Y, Z def cmFromColorList(colorList, tInitial=1, tFinal=1, N=256, gamma=1.): ''' Create colormap from colorList (e.g., ['red','green','blue'] or ['#0000ff','#ff0000']).''' from matplotlib.colors import ListedColormap,LinearSegmentedColormap cmap = LinearSegmentedColormap.from_list("cmFromColorList", colorList, N=N, gamma=gamma) # Create colormap from colorList if not (tInitial==1 and tFinal==1): tCmap = cmap(np.arange(cmap.N)) # Prepare colormap with transparency tCmap[:,-1] = np.linspace(tInitial, tFinal, cmap.N) # Append transparency values (alpha channel) cmap=ListedColormap(tCmap) return cmap #-- Get uncertainty ellipse from covariance matrix of fit. def covarEllipse(X,Y,sigma1,sigma2,mu1,mu2,angle): ''' Just brute force calculate the z-values (amplitude) for an elliptical cone with widths corresponding to sigma uncertainties on the plotting grid. ''' CE = [] # covariance ellipse cone a = -np.deg2rad(angle) # Angle to radians for x in X[:,0]: # For each point along x r = [] # New row for y in Y[0,:]: # For each point along y x0 = (x-mu1)*np.cos(a) + (y-mu2)*np.sin(a) # Calculate x pos. on unrotated grid y0 = (y-mu2)*np.cos(a) - (x-mu1)*np.sin(a) # Calculate y pos. on unrotated grid r.append( ((abs(x0)/sigma1/2)**2 + (abs(y0)/sigma2/2)**2)**0.5) # Value of sigma-cone at current position CE.append(r) return np.array(CE) # Calculate ellipse on grid, with widths sigma, center mu. def getCE(mu_x=0,mu_y=0): ''' Calculate covariance ellipse for rC, rH parameters from fitUI() fit. Optional displacement of center by mu_x, mu_y. ''' global GUI (popt, pcov, infodict, errmsg, ier) = GUI._Fit # Unpack fit results rCf,rHf = popt[:2] # Fitted rC, rH bond lengths residuals = infodict['fvec'] # Residuals for uncertainty calculation ... chi_sq = (residuals**2).sum() # Chi squared (weigthed) DOF = len(residuals)-len(popt) # Fit degrees of freedom s_sq = chi_sq / DOF # ... reduced chi**2 rcov = pcov * s_sq # ... covariances (cf. scipy.optimize.curve_fit) pSigma = np.sqrt(np.diag(rcov)) # ... estimated 1-sigma parameter uncertainties eVal, eVec = np.linalg.eigh(rcov[:2,:2]) # Variance (eigenvalues, vectors) for rCC, rCH sig1,sig2 = eVal[0]**0.5, eVal[1]**0.5 # Width of distribution sigma (along principal axis) #mu_x,mu_y = 0,0 # Center of normal distribution; allow offset through call angle = np.rad2deg(np.arctan2(*eVec[:,0][::-1])) # covariance angle (in degrees for sanity check) return covarEllipse(X,Y, sig1*1000,sig2*1000, mu_x,mu_y, angle) # Calculate covariance ellipse for rC, rH parameters from last GUI fit. #--- McSim distribution from Monte Carlo Simulation (Data from last MC-Sim. performed in GUI) ---# mcSim = npa(GUI.mcSim) # Get Monte-Carlo results from GUI program rC, rH, dH = mcSim[:,0], mcSim[:,1], mcSim[:,6] # Unpack rC, rH, dH values arC, arH = np.average(rC), np.average(rH) # Average results for rC, rH rmC, rmH = 1.39139, 2.47280 # Best estimate for rC,rH from rm2 fit (use as graph origin) rH0 = mcSim[:,1] + dH * np.sqrt( 79 / (1*78)) # r_0(H) from rm and Laurie parameter, in C6H5D rD0 = mcSim[:,1] + dH * np.sqrt( 79 / (2*77)) # r_0(D) from rm and Laurie parameter, in C6H5D rD0 = mcSim[:,1] + dH * np.sqrt( 84 / (2*82)) # r_0(D) from rm and Laurie parameter, in C6D6 # Comment error propagation for rH0,rD0 via simple geometric sum of uC + uH + udH*LaurieFactor deltaH = np.mean(rH0) - np.mean(rH) # Rovib correction between r0(H) - re(H) deltaD = np.mean(rD0) - np.mean(rH) # Rovib correction between r0(D) - re(H) np.mean(rH0)-np.mean(rD0) # Mean rCH - rCD bond length difference from MC sim. rC_,rH0_,rD0_ = (rC-rmC)*1000, (rH0-rmH)*1000, (rD0-rmH)*1000 # bond lengths relative to r_m (in mA) _x,_y = 0.2,1 # Border for graph (in mA) xmin,xmax = -0.25, rC_.max()+_x # Graph size ymin,ymax = -2.5, rH0_.max()+_y # ... X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] # Create plotting grid (x,y) Xh, Yh, Zh = density_estimation(rC_, rH0_,X,Y) # Density estimate for rH vs. rC, mapped on grid Xd, Yd, Zd = density_estimation(rC_, rD0_,X,Y) # Density estimate for rD vs. rC, mapped on grid Zhd = Zh/(Zh.max()*0.8) + Zd/(Zd.max()*0.8) # Joint normalized density map for H,D #--- Plot MC distribution ---# plt.rcParams["figure.figsize"] = (6,4) from matplotlib.ticker import MultipleLocator fig, ax = plt.subplots() ax.imshow(np.rot90(Zhd), cmap=cmFromColorList(['#000000','#B00505'],0,0.5), extent=[xmin, xmax, ymin, ymax], aspect='auto', interpolation='bilinear') # Show McSim density map plt.contour(Xh, Yh, Zh, zorder=1, cmap=cmFromColorList(['#0000a0','#c00808'],0,1)) # Show McSim H contour lines plt.contour(Xd, Yd, Zd, zorder=1, cmap=cmFromColorList(['#0000a0','#c00808'],0,1)) # Show McSim D contour lines #CH = plt.contour(X, Y, Zch, zorder=3, cmap=cmFromColorList(['#0000a0','#04047099'],0,0.6)) # Show fit H contour lines #ax.clabel(CH, inline=True, fontsize=10) ax.scatter((rC_), (rH0_), c='#04047099', s=12, edgecolor='none',label=f'rCH', zorder=2) # Show McSim H scatter ax.scatter((rC_), (rD0_), c='#04550499', s=12, edgecolor='none',label=f'rCD', zorder=2) # Show McSim D scatter fig.canvas.manager.set_window_title('Parameter correlation') ax.set(xlabel=f'r(C\u2013C) [mA] - {rmC:.4f} A', ylabel=f'r(C\u2013H,D) [mA] - {rmH-rmC:.4f} A') ax.legend() fig.tight_layout() plt.show() plt.savefig(GUI.dir+'/parameter_correlation.png', dpi=300) #--- Estimate fit parameter covariance and prepare covariance ellipses for plot ---# #! --> Make sure the fit data and fit type is properly set-up in the GUI (as for MC-Sim.) # We now run rm(2), rm(2)L fits in GUI to calculate the uncertanity ellipses GUI.laurieBox.setChecked(False) # rm(2) fit setup; do not fit Laurie parameter GUI.Fit() # Run the fit CErm2 = getCE() # Covariance ellipse for rm2 fit (center at origin) GUI.laurieBox.setChecked(True) # rm(2)L fit setup; fit the Laurie parameter GUI.Fit() # Run the fit CErm2L = getCE() # Covariance ellipse for rm2L fit #--- Plot MC distribution and covariance ellipses for rm2, rm2L fits ---# plt.rcParams["figure.figsize"] = (5,4) fig, ax = plt.subplots() ax.imshow(np.rot90(Zhd), cmap=cmFromColorList(['#000000','#B00505'],0,0.4), extent=[xmin, xmax, ymin, ymax], aspect='auto', interpolation='bilinear', zorder=5) # Show McSim density map CA = plt.contour(X, Y, Zh, linewidths=1, cmap=cmFromColorList(['#0000a000','#400408A0','#500507B0','#600808C0']), zorder=2) # Show McSim H contour lines CB = plt.contour(X, Y, Zd, linewidths=1, cmap=cmFromColorList(['#00a00000','#400704A0','#500605B0','#600808C0']), zorder=2) # Show McSim D contour lines ax.scatter((rC_), (rH0_), c='#04047060', s=12, edgecolor='none',label=f'r(C-H)', zorder=1) # Show McSim H scatter ax.scatter((rC_), (rD0_), c='#04550460', s=12, edgecolor='none',label=f'r(C-D)', zorder=1) # Show McSim D scatter Crm2L = plt.contour(X, Y, CErm2L, cmap=cmFromColorList(['#A00000','#A00000'],0.6,0.2),linewidths=1.6,levels=[1,2,3,4,5], zorder=0) # Show sigma-ellipse countour lines ax.clabel(Crm2L, inline=True, fontsize=10, fmt='%d-\u03c3') plt.contour(X, Y, CErm2L, cmap=cmFromColorList(['#A00000','#A00000'],0.6,0.2),linewidths=0.5,levels=[0.5,1.5,2.5,3.5,4.5], zorder=-2) # Show odd sigma-ellipse countour lines #Crm2 = plt.contour(X, Y, CErm2, cmap=cmFromColorList(['#333333','#ebebeb']),linewidths=1.,levels=[5,10,15,20], zorder=-2) # Show even sigma-ellipse countour lines Crm2 = plt.contour(X, Y, CErm2, cmap=cmFromColorList(['#333333','#ebebeb']),linewidths=1.,levels=np.arange(0,14,2), zorder=-2) # Show even sigma-ellipse countour lines ax.clabel(Crm2, inline=True, fontsize=10, fmt='%d-\u03c3') #plt.contour(X, Y, CErm2, cmap=cmFromColorList(['#222222','#f0f0f0']), linewidths=0.3, levels=[2.5,7.5,12.5,17.5,22.5], zorder=-3) # Show odd sigma-ellipse countour lines plt.contour(X, Y, CErm2, cmap=cmFromColorList(['#222222','#f0f0f0']), linewidths=0.3, levels=np.arange(0,18,1), zorder=-3) # Show odd sigma-ellipse countour lines #levels=[1,2,3,5,6,7,9,10,11,13,14,15,17,18,19,21,22,23] fig.canvas.manager.set_window_title('Parameter correlation') ax.set(xlabel=f'r(C\u2013C) [mA] - {rmC:.4f} A', ylabel=f'r(C\u2013H,D) [mA] - {rmH-rmC:.4f} A') leg = ax.legend(loc=2) # Legend for scatterplot data leg.legendHandles[0].set_color('#04047090') # Give less transparent legend markers for scatterplot datapoints leg.legendHandles[1].set_color('#04550490') ax.set_xlim(left=-0.19, right=1.48) # Plot origin inportant to have legible sigma-ellipse markers in the field-of-view ax.yaxis.set_minor_locator(MultipleLocator(0.5)) ax.xaxis.set_minor_locator(MultipleLocator(0.1)) ax.set_ylim(bottom=-1.9, top=16.4) fig.tight_layout() plt.show() """ #------------------------------------------------------# #--- SECTION (4) Isotopologue data ---# #------------------------------------------------------# # To allow the easy reproduction of literature data, I give some rotational constant data-sets below # and include a method to import the data into the fit program. # (4.0) Function to convert rot. frequencies to inertial moments and prepare a fit list compatible with fitUI. def expandFitList(isoList, fitAxis=[1,1,1], fitWeights=True): """ Create fitting list from isotopeList with one entry for each isotopologue (a,b,c axes). Concurrently convert from B [MHz] to I [amu*A**2] and calculate fit-weights. Expanded list: [ isotope, ([inertial moments, uncertainties]), [fit weights], author], [ isotope, ( [Ia, Ib, Ic], [sIa, sIb, sIa]), [wa, wb, wc], author], e.g.: ['XCCCCCHHHHHH', ([88.8,90.7,179.0], [2e-4,2e-4,0.1]), [316,309,0.],'In2021'], """ fitList = [] for row in isoList: # Cycle through rows in isoList name, [const,sigma], fitABC, author = row # Parse row # Calculate Ia, Ib, Ic and uncertainties from rotational constants A,B,C, sA,sB,sC Ib = MHz2I(const[1]) # Ib from rot.constant B Ia = MHz2I(const[0]) if const[0] else np.nan #Ib # Ia (=Ib for symmetric isotopologue) Ic = MHz2I(const[2]) if const[2] else np.nan #2*Ib # Ic (=2Ib for symmetric isotopologue) I = [Ia,Ib,Ic] sIb = sigma[1]/const[1]*Ib # Uncertainty sIb from sigma_B sIa = sigma[0]/const[0]*Ia if sigma[0] else np.nan#sIb # sIa (=sIb for symmetric isotopologue) sIc = sigma[2]/const[2]*Ic if sigma[2] else np.nan#2*sIb # sIc (=2*sIb for symmetric isotopologue) sI = [sIa,sIb,sIc] if not fitWeights: # Don't use fit weights? weights = [1 if w>0 else 0 for w in npa(fitABC)*fitAxis] # Replace weights with 1 else: weights = [i/si for i,si in zip(npa(I)*fitABC*fitAxis,sI)] # Calculate fit weights fitList.append([name, ([Ia,Ib,Ic],[sIa,sIb,sIc]), weights, author]) return fitList # (4.1) Data for benzene (incl. Literature) #-------------------------------------------------------------------- # We try to give a reasonably complete account of literaeture constants for benzene isotopologues. # -> For C6H6, an excessively large and diverse number of published rotational constants are available # (see Lee2019 for a summary). Here we only give values from the last 40 years and with error bounds # below 100 kHz. # -> Kunishige2014 values for CCCCCCDHHHHH, CCCCCCDDHHHH, and CCCCCCDHDHHH are taken from Oldani but were # re-fit with fixed (estimated) distortion constants and, in one case, including an additional transition. # Kunishige did not fit distortion constants but fixed those to averged values between C6H6 and C6D6 and # stated that "their [Oldani's] rotational constants are considered to be more accurate than ours". We # included bothe sets of values. # -> The Doi2004 B-constant for C6D6 is significantly lower than other literature values. We therefore # considered this values as an outliers and did not include it in the fit. # -> Pliva 1989, 1989a and 1990 analyzed different vibrational bands for C6D6 (1989a, 1990) and 13C6H6 (1989, 1990). # For C6D6, the later manuscript states "The ground state constants from our previous work (I) [i.e., 1989a] # were fixed in all calculations", but the manuscript nevertheless contains slightly different ground state # constants with larger error bars. Lacking further information, we give relavant constants from all manuscripts. # # We give a table with symmetric isotopologue definitions, i.e., the resulting Cartesian coordinates generated by # the corresponding int2cart function are symmetric with respect to the x,y axes. This allows us to apply rovibrational # correction terms, in rm(1), rm(2) fits, along identical molecular coordinate axes for all isotopologues without # the need for coordinate system rotation. Fits for asymmetric coordinate definitions should use rm(1)r, rm(2)r # fit functions to actively rotate the coordinate system before rovibrational corrections -- but this can lead to # artifacts because 1st order, 2nd order and Laurie corrections are not applied along identical axes. We found # that 'random' asymmetric coordinate definition gives similar results as the symmetric coordinate definition, # but the chosen symmetric definition is better, because it removes the need for additional off-diagonal 1st order # correction constants and removes all potential fit coordinate rotation artifacts. #[ Isotopologue, [ (rot._constants_A,B,C) , (uncertainty_A,B,C) ],[fit_A,B,C], Author ], cList_sym = [ # List of [isotopologe composition, [(rotational constants), (1-sigma uncertainties)], [fit-mask], author string] ['CCCCCCHHHHHH', [( nan, 5689.248415, nan), ( nan, 0.045000, nan)], [0, 1, 0], 'Lindenmayer1988'], ['CCCCCCHHHHHH', [( nan, 5689.266400, nan), ( nan, 0.006000, nan)], [0, 1, 0], 'Hollenstein1990'], ['CCCCCCHHHHHH', [( nan, 5689.241000, nan), ( nan, 0.013000, nan)], [0, 1, 0], 'Domenech1991'], ['CCCCCCHHHHHH', [( nan, 5689.278100, nan), ( nan, 0.001000, nan)], [0, 1, 0], 'Juntilla1991'], ['CCCCCCHHHHHH', [( nan, 5689.242000, nan), ( nan, 0.015000, nan)], [0, 1, 0], 'Okruss1999'], ['CCCCCCHHHHHH', [( nan, 5689.212440, nan), ( nan, 0.008994, nan)], [0, 1, 0], 'Doi2004'], ['CCCCCCHHHHHH', [( nan, 5689.267100, nan), ( nan, 0.005200, nan)], [0, 1, 0], 'Lee2019'], ['CCCCCCHHHHHH', [( nan, 5689.285500, nan), ( nan, 0.005400, nan)], [0, 1, 0], 'In2021'], ['XCCCCCHHHHHH', [( 5689.474000, 5568.473000, 2868.600000), ( 0.018000, 0.023000, 7.300000)], [1, 1, 0], 'In2021'], ['CCCCCCDHHHHH', [( 5689.144000, 5323.934000, 2749.674000), ( 0.006000, 0.006000, 0.006000)], [1, 1, 1], 'Oldani1984'], ['CCCCCCDHHHHH', [( 5689.143488, 5323.933318, 2749.675439), ( 0.005996, 0.005996, 0.005996)], [1, 1, 1], 'Kunishige2015'], ['CCCCCCHDDHHH', [( 5498.062000, 5164.242000, 2662.496000), ( 0.004000, 0.004000, 0.004000)], [1, 1, 1], 'Oldani1988'], ['CCCCCCHDDHHH', [( 5498.031792, 5164.212890, 2662.465813), ( 0.008994, 0.008994, 0.005996)], [1, 1, 1], 'Kunishige2015'], ['CCCCCCHHDHDH', [( 5502.669000, 5152.057000, 2660.358000), ( 0.007000, 0.005996, 0.005996)], [1, 1, 1], 'Oldani1988'], ['CCCCCCHHDHDH', [( 5502.666583, 5152.053308, 2660.328293), ( 0.008994, 0.008994, 0.008994)], [1, 1, 1], 'Kunishige2015'], ['CCCCCCHHDDDH', [( 5168.017256, 5151.933391, 2579.579194), ( 0.014990, 0.059958, 0.005996)], [1, 1, 1], 'Kunishige2015'], ['CCCCCCDDDDHH', [( 5163.715234, 4846.813621, 2499.792430), ( 0.005996, 0.005996, 0.002998)], [1, 1, 1], 'Kunishige2015'], ['CCCCCCDHDDDH', [( 5151.993349, 4850.312199, 2497.900739), ( 0.119917, 0.119917, 0.029979)], [1, 1, 1], 'Kunishige2015'], ['CCCCCCHDDDDD', [( 4998.169839, 4707.221259, 2423.971919), ( 0.149896, 0.149896, 0.005996)], [1, 1, 1], 'Kunishige2015'], ['CCCCCCDDDDDD', [( nan, 4707.125325, nan), ( nan, 0.005996, nan)], [0, 0, 0], 'Doi2004'], ['CCCCCCDDDDDD', [( nan, 4707.311796, nan), ( nan, 0.051864, nan)], [0, 1, 0], 'Pliva1989'], ['CCCCCCDDDDDD', [( nan, 4707.311196, nan), ( nan, 0.149896, nan)], [0, 1, 0], 'Pliva1990'], ['CCCCCCDDDDDD', [( nan, 4707.327685, nan), ( nan, 0.005996, nan)], [0, 1, 0], 'Snels2002'], ['CCCCCCDDDDDD', [( nan, 4707.317500, nan), ( nan, 0.003400, nan)], [0, 1, 0], 'In2021'], ['XCCCCCDDDDDD', [( 4707.541000, 4624.188400, 2332.000000), ( 0.036000, 0.031000, 16.000000)], [1, 1, 0], 'In2021'], ['XXXXXXHHHHHH', [( nan, 5337.915638, nan), ( nan, 0.020985, nan)], [0, 1, 0], 'Pliva1989a'], ['XXXXXXHHHHHH', [( nan, 5337.924632, nan), ( nan, 0.059958, nan)], [0, 1, 0], 'Pliva1990'], ['XXXXXXHHHHHH', [( nan, 5337.884000, nan), ( nan, 0.051000, nan)], [0, 1, 0], 'In2021'], ['XXXXXXDDDDDD', [( nan, 4464.371380, nan), ( nan, 0.023983, nan)], [0, 1, 0], 'Pliva1991'], ] benzeneData = expandFitList(cList_sym) # Convert cList data to correct format #GUI.setData(benzeneData) # Write data into GUI program interface #GUI.r_internal.setText('rC:1 rH:2') # Set reasonable geometry parameters # (4.2) Literature data for H2O #-------------------------------------------------------------------- # Test data to reproduce fit results in Watson1999 and STRFIT. # H2O isotopologue data; H=H(1), D=H(2), T=H(3), O=O(16), P=O(17), Q=O(18) # Kyro1981 [E. Kyrö, J. Mol. Spectr. 88, 167 (1981).] # Toth1993a [R.A. Toth, J.Mol.Spectr. 162, 20 (1993).] # Toth1993b [R.A. Toth, J.Mol.Spectr. 162, 41 (1993).] # Bellet1973 [J.B. Walter, J.Lafferty, G. Steenbeckeliers, J.Mol.Spectr. 47, 388 (1973).] # Lucia 1973 [F.C. DeLucia, P. Helminger, W. Gordy, H.W. Morgan, P.A. Staats, PRA 8, 2785 (1973).] # Helminger1974 [P. Helminger, F.C. DeLucia, W. Gordy, P.A. Staats, H.W. Morgan, PRA 10, 1072 (1974).] # Values may include rounding errors in the last digits. Some rot. constant values in the STRFIT H2O example # deviate significantly from values calculated here. (Mostly values that involve wavenumber to MHz # conversion, e.g., Kyro1981 values.) # #[Isotopologue, [ (rot._constants_A,B,C), (uncertainty_A,B,C)], [fit_A,B,C], Author ], h2oList0 = [ ['HOH', [(835839.4 , 435346.7 , 278138.5 ), ( 8.7, 4.2, 4.2)], [1,1,1], 'Kyro1981'], ['HPH', [(830279.4 , 435349.1 , 277506.7 ), ( 12.6, 6.3, 6.3)], [1,1,1], 'Kyro1981'], ['HQH', [(825365.7 , 435353.5 , 276948.8 ), ( 10.1, 5.0, 5.0)], [1,1,1], 'Kyro1981'], ['HOD', [(701931.327, 272904.297, 192063.387), (0.092, 0.008, 0.009)], [1,1,1], 'Toth1993a'], ['HOT', [(677849.040, 198197.489, 150462.412), (0.170, 0.128, 0.128)], [1,1,1], 'Helminger1974'], ['DOT', [(410174.145, 172101.952, 119127.850), (0.078, 0.045, 0.045)], [1,1,1], 'Helminger1974'], ['DOD', [(462278.878, 218037.946, 145258.324), (0.092, 0.008, 0.009)], [1,1,1], 'Toth1993b'], ['TOT', [(338810.923, 145665.417, 100259.415), (0.076, 0.044, 0.044)], [1,1,1], 'Lucia1973'], ['DPD', [(456766.9 , 218041.0 , 144701.5 ), ( 10, 10, 10)], [1,1,1], 'Bellet1973'], ['DQD', [(451891.823, 218043.800, 144203.299), (0.191, 0.017, 0.018)], [1,1,1], 'Toth1993b'], ['HQD', [(692830.401, 271580.359, 190701.798), (0.119, 0.105, 0.119)], [1,1,1], 'Toth1993a'], ['HPD', [(697089.019, 272208.086, 191341.460), (0.225, 0.195, 0.210)], [1,1,1], 'Toth1993a'],] # Including a duplicate value for DQD destroys agreement with WAT99 literature fit. # ['DQD', [(451891.9 , 218045.2 , 144201.7 ), ( 10, 10, 10)], [1,1,1], 'Bellet1973'],] # These are the values in the STRFIT input example h2oList = [ ['HOH', [(835839.100, 435347.353, 278139.826), ( 0, 0, 0)], [1,1,1], 'Kyro1981'], ['HPH', [(830283.720, 435350.739, 277511.307), ( 0, 0, 0)], [1,1,1], 'Kyro1981'], ['HQH', [(825367.428, 435353.585, 276950.565), ( 0, 0, 0)], [1,1,1], 'Kyro1981'], ['HOD', [(701931.502, 272912.599, 192055.245), ( 0, 0, 0)], [1,1,1], 'Toth1993a'], ['HOT', [(677849.040, 198197.489, 150462.412), ( 0, 0, 0)], [1,1,1], 'Helminger1974'], ['DOT', [(410174.145, 172101.952, 119127.850), ( 0, 0, 0)], [1,1,1], 'Helminger1974'], ['TOT', [(338810.923, 145665.417, 100259.415), ( 0, 0, 0)], [1,1,1], 'Lucia1973'], ['DOD', [(462278.854, 218038.233, 145258.022), ( 0, 0, 0)], [1,1,1], 'Toth1993b'], ['DPD', [(456766.850, 218041.011, 144701.487), ( 0, 0, 0)], [1,1,1], 'Bellet1973'], ['DQD', [(451891.921, 218045.181, 144201.681), ( 0, 0, 0)], [1,1,1], 'Toth1993b'], ['HQD', [(692872.910, 271598.610, 190714.170), ( 0, 0, 0)], [1,1,1], 'Toth1993a'], ] # Create Watson 1999 fit list from data table h2oFitListWA = expandFitList(h2oList0, fitWeights=False) # Watson1999 fitlist with constant fit weights. # Now use this data. Remember to select the int2cart function, molecular symmetry, parameters, ... #> GUI = fitUI() #> GUI.setData(h2oFitListWA,NISTmass) # Create STRFIT example fit list from data table h2oFitListST = expandFitList(h2oList, fitWeights=False) # STRFIT fitlist with STRFIT example data. STRmass = {'H': 1.0078250, 'D':2.0141018, 'T':3.0160493, 'O': 15.9949146, 'P':16.9991315, 'Q':17.9991604} #> GUI.setData(h2oFitListST,STRmass) # Write data into GUI program interface # (4.3) Literature data for HCN #-------------------------------------------------------------------- # Data to reproduce STRFIT example for HCN (also: Watson1999, Table 7) HCNList = """ Isotopic species 1: B = 44315.97570 HCN Isotopic species 2: B = 43170.14000 HXN Isotopic species 3: B = 43027.64400 HCM Isotopic species 4: B = 41863.95300 HXM Isotopic species 5: B = 36207.46270 DCN Isotopic species 6: B = 35587.61900 DXN Isotopic species 7: B = 35169.79100 DCM Isotopic species 8: B = 34531.27600 DXM """ STRmass = {'H': 1.0078252, 'D':2.0141022, 'C': 12.0000000, 'X':13.0033544, 'N': 14.0030744, 'M':15.0001077} HCNdata, pos = [],1 lines = HCNList.splitlines() M2I = lambda nu: 505379.00914143625/float(nu) # rot. constant string (MHz) to inertial moment in amu*A**2 while (1): # Conversion to fitList format line = lines[pos].split(':')[1].replace("="," ") axisB,B,isotopologue = line.split() assert axisB == 'B' HCNdata.append([isotopologue,([M2I(np.nan),M2I(B),M2I(np.nan)],[0,0,0]),[0,1,0],'STRFIT_values']) pos += 1 if pos >= len(lines): break #> GUI.setData(HCNdata,STRmass) # Write data into GUI program interface # (4.4) Literature data for Ozone, O3 #-------------------------------------------------------------------- # Data to reproduce STRFIT example for O3 (also: Watson1999, Table 8) O3List = """ Parent species 1: A =106536.23526 # OOO B = 13349.25468 C = 11834.36144 Isotopic species 2: A =102351.08900 # OPO B = 13351.08930 C = 11781.71900 Isotopic species 3: A =105490.95800 # POO B = 12951.27360 C = 11508.02890 Isotopic species 4: A =103511.75600 # POQ B = 12210.30620 C = 10897.39910 Isotopic species 5: A = 97601.72740 # PQO B = 12954.26740 C = 11408.73130 Isotopic species 6: A =100388.70700 # OPQ B = 12592.52830 C = 11162.88320 Isotopic species 7: A = 95622.15900 # PQQ B = 12213.19460 C = 10804.87590 Isotopic species 8: A = 98393.97900 # QPQ B = 11867.17230 C = 10566.24730 Isotopic species 9: A =104437.20400 # POP B = 12562.17840 C = 11187.63470 Isotopic species 10: A = 99326.43500 # PPQ B = 12211.83480 C = 10849.98860 Isotopic species 11: A = 96547.36680 # PQP B = 12565.47870 C = 11092.13060 Isotopic species 12: A =100251.75490 # PPP B = 12563.89603 C = 11138.66285 Isotopic species 13: A =104573.15785 # OOQ B = 12591.53290 C = 11212.50594 Isotopic species 14: A = 98646.40911 # OQO B = 13352.73300 C = 11731.76807 Isotopic species 15: A = 94689.33394 # QQQ B = 11868.63952 C = 10522.74526 Isotopic species 16: A =102579.22989 # QOQ B = 11865.48869 C = 10611.88255 Isotopic species 17: A = 96684.92942 # QQO B = 12593.34881 C = 11115.69277 """ STRmass = {'O':15.9949146221, 'P':16.99913150, 'Q':17.9991604 } # Masses given in STRFIT input O3data, pos = [],1 lines = O3List.splitlines() M2I = lambda nu: 505379.00914143625/float(nu) # rot. constant string (MHz) to inertial moment in amu*A**2 while (1): # Conversion to fitList format line = lines[pos].split(':')[1].replace("="," ") axisA,A,c,isotopologue = line.split() assert axisA == 'A' line = lines[pos+1][14:].replace("="," ") axisB,B, = line.split() assert axisB == 'B' line = lines[pos+2][14:].replace("="," ") axisC,C, = line.split() assert axisC == 'C' O3data.append([isotopologue,([M2I(A),M2I(B),M2I(C)],[0,0,0]),[1,1,1],'STRFIT_values']) pos += 3 if pos >= len(lines): break # Conversion to fitList format #> GUI.setData(O3data,STRmass) # Write data into GUI program interface # (4.5) Literature data for HNCO #-------------------------------------------------------------------- # Data to reproduce STRFIT example for HNCO # also: Watson1999, Table 9; but note that the values used by Watson are sligthly different # (probably different unit conversion) and rot. constants should be taken directly from the # literature for an exact reproduction of Watson values. HNCOList =""" Isotopic species 1: A =918504.40000 # HNCO B = 11071.00825 C = 10910.57553 Isotopic species 2: A =908968.20000 # HMCO B = 10737.83040 C = 10585.46530 Isotopic species 3: A =916293.80000 # HNXO B = 11071.47980 C = 10910.73070 Isotopic species 4: A =918416.51000 # HNCP B = 10470.89559 C = 10327.24228 Isotopic species 5: A =512472.30000 # DNCO B = 10313.71310 C = 10079.67613 """ STRmass = {'H': 1.00782503, 'D':2.0141022,'N': 14.00307401, 'M':15.0001077,"C": 12.0000000, 'X':13.0033548,'O': 15.9949146221, 'P':17.9991604} HNCOdata, pos = [],1 lines = HNCOList.splitlines() M2I = lambda nu: 505379.00914143625/float(nu) # rot. constant string (MHz) to inertial moment in amu*A**2 while (1): line = lines[pos].split(':')[1].replace("="," ") axisA,A,c,isotopologue = line.split() assert axisA == 'A' line = lines[pos+1][14:].replace("="," ") axisB,B, = line.split() assert axisB == 'B' line = lines[pos+2][14:].replace("="," ") axisC,C, = line.split() assert axisC == 'C' HNCOdata.append([isotopologue,([M2I(A),M2I(B),M2I(C)],[0,0,0]),[1,1,1],'STRFIT_values']) pos += 3 if pos >= len(lines): break # Conversion to fitList format #> GUI.setData(HNCOdata,STRmass) # Write data into GUI program interface #------------------------------------------------------# #--- SECTION (5) Execute from command line ---# #------------------------------------------------------# # Run the GUI program when this file is executed # I prepare the program for the rm(2)L fit of benzene (publication Table 1, row 6). # This fit is quite badly determined because we fit three rovibrational # corrections (rm1, rm2 and Laurie parameter) for only two distinct bonds. # In this case, fitted parameters are very sensitive to the nuderlying data # set and to rounding errors. The exact numerical values in Table 1, row 6 # can be reproduced with 'benzeneData' and with MHz units. if __name__ == "__main__": app = QtWidgets.QApplication(sys.argv) GUI = fitUI() GUI.setData(benzeneData) # Substitute with other data if desired. GUI.r_internal.setText('rC:1 rH:2') # Set reasonable geometry parameters for benzene GUI.units.setCurrentIndex(1) # Fit in MHz (reproduce rounding errors in underdetermined fits) GUI.fitType.setCurrentIndex(6) # Select rm(2) fit GUI.laurieBox.setChecked(True) # Select Laurie correction sys.exit(app.exec_()) if __name__ != "__main__": GUI = fitUI()