X_u_template[0:2] are teff, logg, vturb in km/s
    X_u_template[:,3] -> onward, put atomic number 
    atomic_number is 6 for C, 7 for N
496469
def get_model_spec_ting(atomic_number):
    
    DATA_DIR = "/Users/annaho/Data/LAMOST/Mass_And_Age"
    temp = np.load("%s/X_u_template_KGh_res=1800.npz" %DATA_DIR)
    X_u_template = temp["X_u_template"]
    wl = temp["wavelength"]
    grad_spec = X_u_template[:,atomic_number]
    return wl, grad_spec
1270774
def perform_es_corr(self, lattice, q, step=1e-4):
        
        logger.info("Running Freysoldt 2011 PC calculation (should be " "equivalent to sxdefectalign)")
        logger.debug("defect lattice constants are (in angstroms)" + str(lattice.abc))

        [a1, a2, a3] = ang_to_bohr * np.array(lattice.get_cartesian_coords(1))
        logging.debug("In atomic units, lat consts are (in bohr):" + str([a1, a2, a3]))
        vol = np.dot(a1, np.cross(a2, a3))  # vol in bohr^3

        def e_iso(encut):
            gcut = eV_to_k(encut)  # gcut is in units of 1/A
            return scipy.integrate.quad(lambda g: self.q_model.rho_rec(g * g)**2, step, gcut)[0] * (q**2) / np.pi

        def e_per(encut):
            eper = 0
            for g2 in generate_reciprocal_vectors_squared(a1, a2, a3, encut):
                eper += (self.q_model.rho_rec(g2)**2) / g2
            eper *= (q**2) * 2 * round(np.pi, 6) / vol
            eper += (q**2) * 4 * round(np.pi, 6) \
                * self.q_model.rho_rec_limit0 / vol
            return eper

        eiso = converge(e_iso, 5, self.madetol, self.energy_cutoff)
        logger.debug("Eisolated : %f", round(eiso, 5))

        eper = converge(e_per, 5, self.madetol, self.energy_cutoff)

        logger.info("Eperiodic : %f hartree", round(eper, 5))
        logger.info("difference (periodic-iso) is %f hartree", round(eper - eiso, 6))
        logger.info("difference in (eV) is %f", round((eper - eiso) * hart_to_ev, 4))

        es_corr = round((eiso - eper) / self.dielectric * hart_to_ev, 6)
        logger.info("Defect Correction without alignment %f (eV): ", es_corr)
        return es_corr
1283720
def Planck_frequency(nu, T):
    '''
    '''
    return 2*hPlanck*nu**3/c_light**2/(exp(hPlanck*nu/kBoltzmann/T)-1)
1218145
def as_phononwebsite(self):
        
        d = {}

        #define the lattice
        d["lattice"] = self.structure.lattice._matrix.tolist()
      
        #define atoms
        atom_pos_car = []
        atom_pos_red = []
        atom_types = []
        for site in self.structure.sites:
            atom_pos_car.append(site.coords.tolist())
            atom_pos_red.append(site.frac_coords.tolist())
            atom_types.append(site.species_string)

        #default for now
        d["repetitions"] = get_reasonable_repetitions(len(atom_pos_car))
    
        d["natoms"] = len(atom_pos_car)
        d["atom_pos_car"] = atom_pos_car
        d["atom_pos_red"] = atom_pos_red
        d["atom_types"] = atom_types
        d["atom_numbers"] = self.structure.atomic_numbers
        d["formula"] = self.structure.formula
        d["name"] = self.structure.formula

        #get qpoints
        qpoints = []
        for q in self.qpoints:
            qpoints.append(list(q.frac_coords))
        d["qpoints"] = qpoints

        # get labels
        hsq_dict = collections.OrderedDict()
        for nq,q in enumerate(self.qpoints):
            if q.label is not None:
                hsq_dict[nq] = q.label

        #get distances
        dist = 0
        nqstart = 0
        distances = [dist]
        line_breaks = []
        for nq in range(1,len(qpoints)):
            q1 = np.array(qpoints[nq])
            q2 = np.array(qpoints[nq-1])
            #detect jumps
            if ((nq in hsq_dict) and (nq-1 in hsq_dict)):
                if (hsq_dict[nq] != hsq_dict[nq-1]):
                    hsq_dict[nq-1] += "|"+hsq_dict[nq]
                del hsq_dict[nq]
                line_breaks.append((nqstart,nq))
                nqstart = nq
            else:
                dist += np.linalg.norm(q1-q2)
            distances.append(dist)
        line_breaks.append((nqstart,len(qpoints)))
        d["distances"] = distances
        d["line_breaks"] = line_breaks
        d["highsym_qpts"] = list(hsq_dict.items())

        #eigenvalues
        thz2cm1 = 33.35641
        bands = self.bands.copy()*thz2cm1
        d["eigenvalues"] = bands.T.tolist()

        #eigenvectors
        eigenvectors = self.eigendisplacements.copy()
        eigenvectors /= np.linalg.norm(eigenvectors[0,0])
        eigenvectors = eigenvectors.swapaxes(0,1)
        eigenvectors = np.array([eigenvectors.real, eigenvectors.imag])
        eigenvectors = np.rollaxis(eigenvectors,0,5)
        d["vectors"] = eigenvectors.tolist()
 
        return d
1263605
def draw_multiplet(ax, fine_state, p, hmin, w, fside='right',
                   label_separation=1, label_fontsize=15, fsize=10,
                   deltanu_fontsize=6, proportional=False, text='',
                   text_pos='top', magnetic_lines=False, **kwds):
    """
    # We determine the vertical positions, calculated from p[1] up.
    hyperfine_states = make_list_of_states([fine_state], 'hyperfine')

    h_list = [ei.nu - hyperfine_states[0].nu for ei in hyperfine_states]
    h_list = [i/h_list[-1] for i in h_list]
    h_min = min([h_list[i+1]-h_list[i] for i in range(len(h_list)-1)])
    h_list = [hmin*i/h_min + p[1] for i in h_list]

    if proportional:
        h_list = [p[1]+i*hmin for i in range(len(hyperfine_states))]

    omegaij = [(hyperfine_states[i+1].nu-hyperfine_states[i].nu)/1e6
               for i in range(len(hyperfine_states)-1)]

    for i in range(len(h_list)):
        label = '$\mathrm{F}='+str(hyperfine_states[i].f)+'$'
        if magnetic_lines:
            maxf = max([eee.f for eee in hyperfine_states])
            f = hyperfine_states[i].f
            nm = 2*maxf+1
            for mf in range(-f, f+1):
                draw_state(ax, [p[0]+mf*w/nm, h_list[i]], "", w/nm*0.5,
                           alignment=fside, fontsize=fsize)

            if fside == 'right':
                ax.text(p[0]+w+label_separation, h_list[i], label,
                        fontsize=fsize, horizontalalignment="right",
                        verticalalignment="center")
            elif fside == 'left':
                ax.text(p[0]-w-label_separation, h_list[i], label,
                        fontsize=fsize, horizontalalignment="left",
                        verticalalignment="center")
        else:
            draw_state(ax, [p[0], h_list[i]], label, w,
                       alignment=fside, fontsize=fsize)

    for i in range(len(h_list)-1):
        hmid = (h_list[i+1]+h_list[i])/2.0-0.5
        nu = str(omegaij[i])[:5]
        if fside == 'left':
            ax.text(p[0]-w/2.0, hmid, r'$'+nu+' \ \mathrm{MHz}$',
                    fontsize=deltanu_fontsize,
                    horizontalalignment=fside, verticalalignment='bottom')
        else:
            ax.text(p[0]+w/2.0, hmid, r'$'+nu+' \ \mathrm{MHz}$',
                    fontsize=deltanu_fontsize,
                    horizontalalignment=fside, verticalalignment='bottom')

    a = label_separation

    if text != '':
        if text_pos == 'top':
            labelx = p[0]
            labely = h_list[-1]+a

            ax.text(labelx, labely, '$'+text+'$',
                    verticalalignment='bottom',
                    horizontalalignment='center', fontsize=label_fontsize)
        elif text_pos == 'right':
            labelx = p[0]+w/2+2.0*a
            if fside == 'right': labelx = labelx+a*5.0
            labely = (h_list[-1]+h_list[0])/2.0

            ax.text(labelx, labely, '$'+text+'$',
                    verticalalignment='center',
                    horizontalalignment='left', fontsize=label_fontsize)

        elif text_pos == 'left':
            labelx = p[0]-w/2-2.0*a
            if fside == 'left': labelx = labelx-a*5.0
            labely = (h_list[-1]+h_list[0])/2.0

            ax.text(labelx, labely, '$'+text+'$',
                    verticalalignment='center',
                    horizontalalignment='right', fontsize=label_fontsize)

    return [[p[0], i] for i in h_list]
821503
void buildPT_T() {
    int numTags = tagIndex.size();
    m_TT = new double[numTags][numTags];
    m_T = new double[numTags];
    double[] tmp = new double[numTags];
    for (IntTaggedWord word : words) {
      double tot = 0.0;
      for (int t = 0; t < numTags; t++) {
        IntTaggedWord iTW = new IntTaggedWord(word.word, t);
        tmp[t] = seenCounter.getCount(iTW);
        tot += tmp[t];
      }
      if (tot < 10) {
        continue;
      }
      for (int t = 0; t < numTags; t++) {
        for (int t2 = 0; t2 < numTags; t2++) {
          if (tmp[t2] > 0.0) {
            double c = tmp[t] / tot;
            m_T[t] += c;
            m_TT[t2][t] += c;
          }
        }
      }
    }
  }
1031608
def proc_atom(self):
        """"""
        atom_data = self.proc_line_coordinate(self.current_line)
        (at_type, at_ser, at_name, alt_loc, res_name, chain_id, res_seq,
         i_code, x, y, z, occupancy, temp_factor, element, charge) = atom_data
        # currently active state
        a_state = self.pdb_parse_tree['data'][self.state]
        res_id = (res_seq, i_code)
        if chain_id not in a_state:
            a_state[chain_id] = (set(), OrderedDict())
        if res_id not in a_state[chain_id][1]:
            a_state[chain_id][1][res_id] = (set(), OrderedDict())
        if at_type == 'ATOM':
            if res_name in standard_amino_acids.values():
                poly = 'P'
            else:
                poly = 'N'
        else:
            poly = 'H'
        a_state[chain_id][0].add((chain_id, at_type, poly))
        a_state[chain_id][1][res_id][0].add(
            (at_type, res_seq, res_name, i_code))
        if at_ser not in a_state[chain_id][1][res_id][1]:
            a_state[chain_id][1][res_id][1][at_ser] = [atom_data]
        else:
            a_state[chain_id][1][res_id][1][at_ser].append(atom_data)
        return
429316
protected void newAtomData() {
        atomCounter = 0;
        elsym = new ArrayList<String>();
        elid = new ArrayList<String>();
        eltitles = new ArrayList<String>();
        formalCharges = new ArrayList<String>();
        partialCharges = new ArrayList<String>();
        isotope = new ArrayList<String>();
        atomicNumbers = new ArrayList<String>();
        exactMasses = new ArrayList<String>();
        x3 = new ArrayList<String>();
        y3 = new ArrayList<String>();
        z3 = new ArrayList<String>();
        x2 = new ArrayList<String>();
        y2 = new ArrayList<String>();
        xfract = new ArrayList<String>();
        yfract = new ArrayList<String>();
        zfract = new ArrayList<String>();
        hCounts = new ArrayList<String>();
        atomParities = new ArrayList<String>();
        parityARef1 = new ArrayList<String>();
        parityARef2 = new ArrayList<String>();
        parityARef3 = new ArrayList<String>();
        parityARef4 = new ArrayList<String>();
        atomAromaticities = new ArrayList<String>();
        atomDictRefs = new ArrayList<String>();
        spinMultiplicities = new ArrayList<String>();
        occupancies = new ArrayList<String>();
        atomCustomProperty = new HashMap<Integer, List<String>>();
    }
491778
def write(fn, a, **kwargs):
    
    ext = fn[fn.rfind('.'):].split('@')
    if ext[0] == '.out' or ext[0] == '.dat':
        return write_atoms(fn, a)
    elif ext[0] == '.lammps':
        return write_lammps_data(fn, a, velocities=True, **kwargs)
    elif ext[0] == '.nc':
        return NetCDFTrajectory(fn, 'w').write(a)
    else:
        return ase.io.write(fn, a, **kwargs)
1138643
def give_zip_element_z_and_names(element_name):
    ''' '''

    #import numpy as np

    global z_bismuth
    z_bismuth = 83
    global z_for_elem
    z_for_elem = []
    global index_stable
    index_stable = []

    i_for_stable = 1
    i_for_unstable = 0
    for i in range(z_bismuth):
        z_for_elem.append(int(i+1))
        # the only elements below bismuth with no stable isotopes are Tc and Pm
        if i+1 == 43 or i+1 == 61:
            index_stable.append(i_for_unstable)
        else:
            index_stable.append(i_for_stable)

        dummy_index = np.zeros(len(element_name))
        for i in range(len(element_name)):
            if element_name[i] == 'Neutron':
                dummy_index[i] = 0
            elif element_name[i] == 'H':
                dummy_index[i] = 1
            elif element_name[i] == 'He':
                dummy_index[i] = 2
            elif element_name[i] == 'Li':
                dummy_index[i] = 3
            elif element_name[i] == 'Be':
                dummy_index[i] = 4
            elif element_name[i] == 'B':
                dummy_index[i] = 5
            elif element_name[i] == 'C':
                dummy_index[i] = 6
            elif element_name[i] == 'N':
                dummy_index[i] = 7
            elif element_name[i] == 'O':
                dummy_index[i] = 8
            elif element_name[i] == 'F':
                dummy_index[i] = 9
            elif element_name[i] == 'Ne':
                dummy_index[i] = 10
            elif element_name[i] == 'Na':
                dummy_index[i] = 11
            elif element_name[i] == 'Mg':
                dummy_index[i] = 12
            elif element_name[i] == 'Al':
                dummy_index[i] = 13
            elif element_name[i] == 'Si':
                dummy_index[i] = 14
            elif element_name[i] == 'P':
                dummy_index[i] = 15
            elif element_name[i] == 'S':
                dummy_index[i] = 16
            elif element_name[i] == 'Cl':
                dummy_index[i] = 17
            elif element_name[i] == 'Ar':
                dummy_index[i] = 18
            elif element_name[i] == 'K':
                dummy_index[i] = 19
            elif element_name[i] == 'Ca':
                dummy_index[i] = 20
            elif element_name[i] == 'Sc':
                dummy_index[i] = 21
            elif element_name[i] == 'Ti':
                dummy_index[i] = 22
            elif element_name[i] == 'V':
                dummy_index[i] = 23
            elif element_name[i] == 'Cr':
                dummy_index[i] = 24
            elif element_name[i] == 'Mn':
                dummy_index[i] = 25
            elif element_name[i] == 'Fe':
                dummy_index[i] = 26
            elif element_name[i] == 'Co':
                dummy_index[i] = 27
            elif element_name[i] == 'Ni':
                dummy_index[i] = 28
            elif element_name[i] == 'Cu':
                dummy_index[i] = 29
            elif element_name[i] == 'Zn':
                dummy_index[i] = 30
            elif element_name[i] == 'Ga':
                dummy_index[i] = 31
            elif element_name[i] == 'Ge':
                dummy_index[i] = 32
            elif element_name[i] == 'As':
                dummy_index[i] = 33
            elif element_name[i] == 'Se':
                dummy_index[i] = 34
            elif element_name[i] == 'Br':
                dummy_index[i] = 35
            elif element_name[i] == 'Kr':
                dummy_index[i] = 36
            elif element_name[i] == 'Rb':
                dummy_index[i] = 37
            elif element_name[i] == 'Sr':
                dummy_index[i] = 38
            elif element_name[i] == 'Y':
                dummy_index[i] = 39
            elif element_name[i] == 'Zr':
                dummy_index[i] = 40
            elif element_name[i] == 'Nb':
                dummy_index[i] = 41
            elif element_name[i] == 'Mo':
                dummy_index[i] = 42
            elif element_name[i] == 'Tc':
                dummy_index[i] = 43
            elif element_name[i] == 'Ru':
                dummy_index[i] = 44
            elif element_name[i] == 'Rh':
                dummy_index[i] = 45
            elif element_name[i] == 'Pd':
                dummy_index[i] = 46
            elif element_name[i] == 'Ag':
                dummy_index[i] = 47
            elif element_name[i] == 'Cd':
                dummy_index[i] = 48
            elif element_name[i] == 'In':
                dummy_index[i] = 49
            elif element_name[i] == 'Sn':
                dummy_index[i] = 50
            elif element_name[i] == 'Sb':
                dummy_index[i] = 51
            elif element_name[i] == 'Te':
                dummy_index[i] = 52
            elif element_name[i] == 'I':
                dummy_index[i] = 53
            elif element_name[i] == 'Xe':
                dummy_index[i] = 54
            elif element_name[i] == 'Cs':
                dummy_index[i] = 55
            elif element_name[i] == 'Ba':
                dummy_index[i] = 56
            elif element_name[i] == 'La':
                dummy_index[i] = 57
            elif element_name[i] == 'Ce':
                dummy_index[i] = 58
            elif element_name[i] == 'Pr':
                dummy_index[i] = 59
            elif element_name[i] == 'Nd':
                dummy_index[i] = 60
            elif element_name[i] == 'Pm':
                dummy_index[i] = 61
            elif element_name[i] == 'Sm':
                dummy_index[i] = 62
            elif element_name[i] == 'Eu':
                dummy_index[i] = 63
            elif element_name[i] == 'Gd':
                dummy_index[i] = 64
            elif element_name[i] == 'Tb':
                dummy_index[i] = 65
            elif element_name[i] == 'Dy':
                dummy_index[i] = 66
            elif element_name[i] == 'Ho':
                dummy_index[i] = 67
            elif element_name[i] == 'Er':
                dummy_index[i] = 68
            elif element_name[i] == 'Tm':
                dummy_index[i] = 69
            elif element_name[i] == 'Yb':
                dummy_index[i] = 70
            elif element_name[i] == 'Lu':
                dummy_index[i] = 71
            elif element_name[i] == 'Hf':
                dummy_index[i] = 72
            elif element_name[i] == 'Ta':
                dummy_index[i] = 73
            elif element_name[i] == 'W':
                dummy_index[i] = 74
            elif element_name[i] == 'Re':
                dummy_index[i] = 75
            elif element_name[i] == 'Os':
                dummy_index[i] = 76
            elif element_name[i] == 'Ir':
                dummy_index[i] = 77
            elif element_name[i] == 'Pt':
                dummy_index[i] = 78
            elif element_name[i] == 'Au':
                dummy_index[i] = 79
            elif element_name[i] == 'Hg':
                dummy_index[i] = 80
            elif element_name[i] == 'Tl':
                dummy_index[i] = 81
            elif element_name[i] == 'Pb':
                dummy_index[i] = 82
            elif element_name[i] == 'Bi':
                dummy_index[i] = 83
            elif element_name[i] == 'Po':
                dummy_index[i] = 84
            elif element_name[i] == 'At':
                dummy_index[i] = 85

        #if spe[0] == 'N   1':
        #       znum_int[0] = 0

        # here the index to connect name and atomic numbers.
        global index_z_for_elements
        index_z_for_elements = {}
        for a,b in zip(element_name,dummy_index):
            index_z_for_elements[a]=b
365301
def get_atomic_stoichiometry(references):
    
    n = len(references)
    stoichiometry = np.ndarray((n, n), float)
    stoichiometry[:] = 0.
    key_index = {}
    for i, (key, species) in enumerate(references):
        # in case species uses a suffix like _gas
        species = species.split('_')[0]
        key_index[key] = i

    for i, (key, species) in enumerate(references):
        composition = ase.symbols.string2symbols(species)
        for j, symbol in enumerate(composition):
            stoichiometry[i, key_index[symbol]] += 1
    istoichiometry = np.linalg.inv(stoichiometry)

    return istoichiometry.tolist()
1349739
function murmur3(data, seed) {
  const tail = data.length - (data.length % 4);
  const c1 = 0xcc9e2d51;
  const c2 = 0x1b873593;
  let h1 = seed;
  let k1;

  for (let i = 0; i < tail; i += 4) {
    k1 = (data[i + 3] << 24)
      | (data[i + 2] << 16)
      | (data[i + 1] << 8)
      | data[i];
    k1 = mul32(k1, c1);
    k1 = rotl32(k1, 15);
    k1 = mul32(k1, c2);
    h1 ^= k1;
    h1 = rotl32(h1, 13);
    h1 = sum32(mul32(h1, 5), 0xe6546b64);
  }

  k1 = 0;
  switch (data.length & 3) {
    case 3:
      k1 ^= data[tail + 2] << 16;
    case 2:
      k1 ^= data[tail + 1] << 8;
    case 1:
      k1 ^= data[tail + 0];
      k1 = mul32(k1, c1);
      k1 = rotl32(k1, 15);
      k1 = mul32(k1, c2);
      h1 ^= k1;
  }

  h1 ^= data.length;
  h1 ^= h1 >>> 16;
  h1 = mul32(h1, 0x85ebca6b);
  h1 ^= h1 >>> 13;
  h1 = mul32(h1, 0xc2b2ae35);
  h1 ^= h1 >>> 16;

  if (h1 < 0)
    h1 += 0x100000000;

  return h1;
}
60390
function () {
    var tmp_0 = this.get(2,2) * this.get(3,3);
    var tmp_1 = this.get(3,2) * this.get(2,3);
    var tmp_2 = this.get(1,2) * this.get(3,3);
    var tmp_3 = this.get(3,2) * this.get(1,3);
    var tmp_4 = this.get(1,2) * this.get(2,3);
    var tmp_5 = this.get(2,2) * this.get(1,3);
    var tmp_6 = this.get(0,2) * this.get(3,3);
    var tmp_7 = this.get(3,2) * this.get(0,3);
    var tmp_8 = this.get(0,2) * this.get(2,3);
    var tmp_9 = this.get(2,2) * this.get(0,3);
    var tmp_10 = this.get(0,2) * this.get(1,3);
    var tmp_11 = this.get(1,2) * this.get(0,3);
    var tmp_12 = this.get(2,0) * this.get(3,1);
    var tmp_13 = this.get(3,0) * this.get(2,1);
    var tmp_14 = this.get(1,0) * this.get(3,1);
    var tmp_15 = this.get(3,0) * this.get(1,1);
    var tmp_16 = this.get(1,0) * this.get(2,1);
    var tmp_17 = this.get(2,0) * this.get(1,1);
    var tmp_18 = this.get(0,0) * this.get(3,1);
    var tmp_19 = this.get(3,0) * this.get(0,1);
    var tmp_20 = this.get(0,0) * this.get(2,1);
    var tmp_21 = this.get(2,0) * this.get(0,1);
    var tmp_22 = this.get(0,0) * this.get(1,1);
    var tmp_23 = this.get(1,0) * this.get(0,1);

    var t0 = ((tmp_0 * this.get(1,1) + tmp_3 * this.get(2,1) + tmp_4 * this.get(3,1)) -
              (tmp_1 * this.get(1,1) + tmp_2 * this.get(2,1) + tmp_5 * this.get(3,1)));
    var t1 = ((tmp_1 * this.get(0,1) + tmp_6 * this.get(2,1) + tmp_9 * this.get(3,1)) -
              (tmp_0 * this.get(0,1) + tmp_7 * this.get(2,1) + tmp_8 * this.get(3,1)));
    var t2 = ((tmp_2 * this.get(0,1) + tmp_7 * this.get(1,1) + tmp_10 * this.get(3,1)) -
              (tmp_3 * this.get(0,1) + tmp_6 * this.get(1,1) + tmp_11 * this.get(3,1)));
    var t3 = ((tmp_5 * this.get(0,1) + tmp_8 * this.get(1,1) + tmp_11 * this.get(2,1)) -
              (tmp_4 * this.get(0,1) + tmp_9 * this.get(1,1) + tmp_10 * this.get(2,1)));

    var d = 1.0 / (this.get(0,0) * t0 + this.get(1,0) * t1 + this.get(2,0) * t2 + this.get(3,0) * t3);

    var out_00 = d * t0;
    var out_01 = d * t1;
    var out_02 = d * t2;
    var out_03 = d * t3;

    var out_10 = d * ((tmp_1 * this.get(1,0) + tmp_2 * this.get(2,0) + tmp_5 * this.get(3,0)) -
                      (tmp_0 * this.get(1,0) + tmp_3 * this.get(2,0) + tmp_4 * this.get(3,0)));
    var out_11 = d * ((tmp_0 * this.get(0,0) + tmp_7 * this.get(2,0) + tmp_8 * this.get(3,0)) -
                      (tmp_1 * this.get(0,0) + tmp_6 * this.get(2,0) + tmp_9 * this.get(3,0)));
    var out_12 = d * ((tmp_3 * this.get(0,0) + tmp_6 * this.get(1,0) + tmp_11 * this.get(3,0)) -
                      (tmp_2 * this.get(0,0) + tmp_7 * this.get(1,0) + tmp_10 * this.get(3,0)));
    var out_13 = d * ((tmp_4 * this.get(0,0) + tmp_9 * this.get(1,0) + tmp_10 * this.get(2,0)) -
                      (tmp_5 * this.get(0,0) + tmp_8 * this.get(1,0) + tmp_11 * this.get(2,0)));

    var out_20 = d * ((tmp_12 * this.get(1,3) + tmp_15 * this.get(2,3) + tmp_16 * this.get(3,3)) -
                      (tmp_13 * this.get(1,3) + tmp_14 * this.get(2,3) + tmp_17 * this.get(3,3)));
    var out_21 = d * ((tmp_13 * this.get(0,3) + tmp_18 * this.get(2,3) + tmp_21 * this.get(3,3)) -
                      (tmp_12 * this.get(0,3) + tmp_19 * this.get(2,3) + tmp_20 * this.get(3,3)));
    var out_22 = d * ((tmp_14 * this.get(0,3) + tmp_19 * this.get(1,3) + tmp_22 * this.get(3,3)) -
                      (tmp_15 * this.get(0,3) + tmp_18 * this.get(1,3) + tmp_23 * this.get(3,3)));
    var out_23 = d * ((tmp_17 * this.get(0,3) + tmp_20 * this.get(1,3) + tmp_23 * this.get(2,3)) -
                      (tmp_16 * this.get(0,3) + tmp_21 * this.get(1,3) + tmp_22 * this.get(2,3)));
    
    var out_30 = d * ((tmp_14 * this.get(2,2) + tmp_17 * this.get(3,2) + tmp_13 * this.get(1,2)) -
                      (tmp_16 * this.get(3,2) + tmp_12 * this.get(1,2) + tmp_15 * this.get(2,2)));
    var out_31 = d * ((tmp_20 * this.get(3,2) + tmp_12 * this.get(0,2) + tmp_19 * this.get(2,2)) -
                      (tmp_18 * this.get(2,2) + tmp_21 * this.get(3,2) + tmp_13 * this.get(0,2)));
    var out_32 = d * ((tmp_18 * this.get(1,2) + tmp_23 * this.get(3,2) + tmp_15 * this.get(0,2)) -
                      (tmp_22 * this.get(3,2) + tmp_14 * this.get(0,2) + tmp_19 * this.get(1,2)));
    var out_33 = d * ((tmp_22 * this.get(2,2) + tmp_16 * this.get(0,2) + tmp_21 * this.get(1,2)) -
                      (tmp_20 * this.get(1,2) + tmp_23 * this.get(2,2) + tmp_17 * this.get(0,2)));

    this.elements[0*4+0] = out_00;
    this.elements[0*4+1] = out_01;
    this.elements[0*4+2] = out_02;
    this.elements[0*4+3] = out_03;
    this.elements[1*4+0] = out_10;
    this.elements[1*4+1] = out_11;
    this.elements[1*4+2] = out_12;
    this.elements[1*4+3] = out_13;
    this.elements[2*4+0] = out_20;
    this.elements[2*4+1] = out_21;
    this.elements[2*4+2] = out_22;
    this.elements[2*4+3] = out_23;
    this.elements[3*4+0] = out_30;
    this.elements[3*4+1] = out_31;
    this.elements[3*4+2] = out_32;
    this.elements[3*4+3] = out_33;
    return this;
  }
194472
def mw(self):
        '''

        '''
        counter = collections.Counter(self.seq.lower())
        mw_a = counter['a'] * 313.2
        mw_t = counter['t'] * 304.2
        mw_g = counter['g'] * 289.2
        mw_c = counter['c'] * 329.2
        mw_u = counter['u'] * 306.2

        if self.material == 'dna':
            return mw_a + mw_t + mw_g + mw_c + 79.0
        else:
            return mw_a + mw_u + mw_g + mw_c + 159.0
360338
def energy_coefficients(m1, m2, s1z=0, s2z=0, phase_order=-1, spin_order=-1):
    
    implemented_phase_order = 7
    implemented_spin_order = 7
    if phase_order > implemented_phase_order:
        raise ValueError("pN coeffiecients of that order have not been implemented")
    elif phase_order == -1:
        phase_order = implemented_phase_order

    if spin_order > implemented_spin_order:
        raise ValueError("pN coeffiecients of that order have not been implemented")
    elif spin_order == -1:
        spin_order = implemented_spin_order

    qmdef1 = 1.0
    qmdef2 = 1.0

    M = m1 + m2
    dm = (m1-m2)/M
    m1M = m1 / M
    m2M = m2 / M

    s1z = s1z * m1M * m1M
    s2z = s2z * m2M * m2M

    _, eta = mass1_mass2_to_mchirp_eta(m1, m2)

    ecof = numpy.zeros(phase_order+1)
    # Orbital terms
    if phase_order >= 0:
        ecof[0] = 1.0
    if phase_order >= 1:
        ecof[1] = 0
    if phase_order >= 2:
        ecof[2] = -(1.0/12.0) * (9.0 + eta)
    if phase_order >= 3:
        ecof[3] = 0
    if phase_order >= 4:
        ecof[4] = (-81.0 + 57.0*eta - eta*eta) / 24.0
    if phase_order >= 5:
        ecof[5] = 0
    if phase_order >= 6:
        ecof[6] = - 675.0/64.0 + ( 34445.0/576.0    \
              - 205.0/96.0 * lal.PI * lal.PI ) * eta  \
              - (155.0/96.0) *eta * eta - 35.0/5184.0 * eta * eta
    # Spin terms

    ESO15s1 = 8.0/3.0 + 2.0*m2/m1
    ESO15s2 = 8.0/3.0 + 2.0*m1/m2

    ESS2 = 1.0 / eta
    EQM2s1 = qmdef1/2.0/m1M/m1M
    EQM2s1L = -qmdef1*3.0/2.0/m1M/m1M
    #EQM2s2 = qmdef2/2.0/m2M/m2M
    EQM2s2L = -qmdef2*3.0/2.0/m2M/m2M

    ESO25s1 = 11.0 - 61.0*eta/9.0 + (dm/m1M) * (-3.0 + 10.*eta/3.0)
    ESO25s2 = 11.0 - 61.0*eta/9.0 + (dm/m2M) * (3.0 - 10.*eta/3.0)

    ESO35s1 = 135.0/4.0 - 367.0*eta/4.0 + 29.0*eta*eta/12.0 + (dm/m1M) * (-27.0/4.0 + 39.0*eta - 5.0*eta*eta/4.0)
    ESO35s2 = 135.0/4.0 - 367.0*eta/4.0 + 29.0*eta*eta/12.0 - (dm/m2M) * (-27.0/4.0 + 39.0*eta - 5.0*eta*eta/4.0)

    if spin_order >=3:
        ecof[3] += ESO15s1 * s1z + ESO15s2 * s2z
    if spin_order >=4:
        ecof[4] += ESS2 * (s1z*s2z - 3.0*s1z*s2z)
        ecof[4] += EQM2s1*s1z*s1z + EQM2s1*s2z*s2z + EQM2s1L*s1z*s1z + EQM2s2L*s2z*s2z
    if spin_order >=5:
        ecof[5] = ESO25s1*s1z + ESO25s2*s2z
    if spin_order >=7:
        ecof[7] += ESO35s1*s1z + ESO35s2*s2z

    return ecof
1135095
function () {
        _gpfBinZ = 36969 * (_gpfBinZ & 65535) + (_gpfBinZ >> 16) & _gpfMax32;
        _gpfBinW = 18000 * (_gpfBinW & 65535) + (_gpfBinW >> 16) & _gpfMax32;
        return ((_gpfBinZ << 16) + _gpfBinW & _gpfMax32) + _gpfMax31;
    }
1270876
def lambda_tilde(mass1, mass2, lambda1, lambda2):
    
    m1, m2, lambda1, lambda2, input_is_array = ensurearray(
        mass1, mass2, lambda1, lambda2)
    lsum = lambda1 + lambda2
    ldiff, _ = ensurearray(lambda1 - lambda2)
    mask = m1 < m2
    ldiff[mask] = -ldiff[mask]
    eta = eta_from_mass1_mass2(m1, m2)
    p1 = (lsum) * (1 + 7. * eta - 31 * eta ** 2.0)
    p2 = (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta ** 2.0) * (ldiff)
    return formatreturn(8.0 / 13.0 * (p1 + p2), input_is_array)