Calculate the ionic dipole moment using ZVAL from pseudopotential

    site: PeriodicSite
    structure: Structure
    zval: Charge value for ion (ZVAL for VASP pseudopotential)

    Returns polarization in electron Angstroms.
975950
def calc_ionic(site, structure, zval):
    
    norms = structure.lattice.lengths_and_angles[0]
    return np.multiply(norms, -site.frac_coords * zval)
743590
def delta13c_santrock(r45sam, r46sam, d13cstd, r45std, r46std,
                      ks='Santrock', d18ostd=23.5):
    
    k = delta13c_constants()[ks]

    # function for calculating 17R from 18R
    def c17(r):
        return k['K'] * r ** k['A']
    rcpdb, rosmow = k['S13'], k['S18']

    # known delta values for the ref peak
    r13std = (d13cstd / 1000. + 1) * rcpdb
    r18std = (d18ostd / 1000. + 1) * rosmow

    # determine the correction factors
    c45 = r13std + 2 * c17(r18std)
    c46 = c17(r18std) ** 2 + 2 * r13std * c17(r18std) + 2 * r18std

    # correct the voltage ratios to ion ratios
    r45 = (r45sam / r45std) * c45
    r46 = (r46sam / r46std) * c46

    def rf(r18):
        return -3 * c17(r18) ** 2 + 2 * r45 * c17(r18) + 2 * r18 - r46
    # r18 = scipy.optimize.root(rf, r18std).x[0]  # use with scipy 0.11.0
    r18 = fsolve(rf, r18std)[0]
    r13 = r45 - 2 * c17(r18)
    return 1000 * (r13 / rcpdb - 1)
1126488
private double calculateCurrentMass() {
        double mass = 0;
        for (int i = 0; i < isotopes.length; i++) {
            mass += currentCounts[i] * isotopes[i].getExactMass();
        }

        return mass;
    }
898091
private static double getRadiusForNucl(NucleotideImpl nuc, Atom atom) {

		if (atom.getElement().equals(Element.H)) return Element.H.getVDWRadius();
		if (atom.getElement().equals(Element.D)) return Element.D.getVDWRadius();

		if (atom.getElement()==Element.C) return NUC_CARBON_VDW;

		if (atom.getElement()==Element.N) return NUC_NITROGEN_VDW;

		if (atom.getElement()==Element.P) return PHOSPHOROUS_VDW;

		if (atom.getElement()==Element.O) return OXIGEN_VDW;

		logger.info("Unexpected atom "+atom.getName()+" for nucleotide "+nuc.getPDBName()+", assigning its standard vdw radius");
		return atom.getElement().getVDWRadius();
	}
201189
def _nac(self, q_direction):
        
        num_atom = self._pcell.get_number_of_atoms()
        nac_q = np.zeros((num_atom, num_atom, 3, 3), dtype='double')
        if (np.abs(q_direction) < 1e-5).all():
            return nac_q

        rec_lat = np.linalg.inv(self._pcell.get_cell())
        nac_factor = self._dynmat.get_nac_factor()
        Z = self._dynmat.get_born_effective_charges()
        e = self._dynmat.get_dielectric_constant()
        q = np.dot(rec_lat, q_direction)

        B = self._B(e, q)
        for i in range(num_atom):
            A_i = self._A(q, Z, i)
            for j in range(num_atom):
                A_j = self._A(q, Z, j)
                nac_q[i, j] = np.outer(A_i, A_j) / B

        num_satom = self._scell.get_number_of_atoms()
        N = num_satom // num_atom

        return nac_q * nac_factor / N
1310321
def rotational_invariance(force_constants,
                          supercell,
                          primitive,
                          symprec=1e-5):
    

    print("Check rotational invariance ...")

    fc = force_constants
    p2s = primitive.get_primitive_to_supercell_map()

    smallest_vectors, multiplicity = primitive.get_smallest_vectors()

    abc = "xyz"

    for pi, p in enumerate(p2s):
        for i in range(3):
            mat = np.zeros((3, 3), dtype='double')
            for s in range(supercell.get_number_of_atoms()):
                m = multiplicity[s, pi]
                vecs = smallest_vectors[s, pi, :m]
                v = np.dot(vecs.sum(axis=0) / m, primitive.get_cell())
                for j in range(3):
                    for k in range(3):
                        mat[j, k] += (fc[p, s, i, j] * v[k] -
                                      fc[p, s, i, k] * v[j])

            print("Atom %d %s" % (p + 1, abc[i]))
            for vec in mat:
                print("%10.5f %10.5f %10.5f" % tuple(vec))
979350
def _parse_radic_header(header_group, dipole_mode="all"):
    
    header = ''.join(header_group[1])

    nr_of_electrodes = _get_nr_of_electrodes(header.split('\n'))
    groups = itertools.groupby(
        header.split('\n'),
        lambda line: (
            line.startswith('[End Readings]') or
            line.startswith('[Begin Readings]')
        )
    )

    for key, group in groups:
        # there is only one group we are interested in: the first one
        if key:
            tmp = np.array(
                [
                    np.fromstring(
                        x, sep=' ', dtype=int
                    )
                    for x in next(groups)[1]
                ],
            )

            # curpots = tmp.T.copy()
            # curpots.resize(nr_of_electrodes + 2, tmp.shape[0])
            break

    # now we have the raw readings, pad them with zeros
    readings = []
    for raw_reading in tmp:
        voltage_rus = raw_reading[3:]
        # note that we now have the last electrode in here, too (so one more
        # than we have RUs. This electrode is always 0
        normed_reading = np.zeros(nr_of_electrodes)
        for i in voltage_rus[np.where(voltage_rus > 0)]:
            normed_reading[i - 1] = i
        readings.append(np.hstack((raw_reading[0:3], normed_reading)))

    # now generate the available configurations
    # add calibration
    reading_configs = {}
    for nr, reading in enumerate(readings):
        # nr = reading[0]
        A = reading[1]
        B = reading[2]

        voltages = []

        # loop over the RU-settings
        # we always measure to the next zero electrode
        for i in range(3, nr_of_electrodes + 3):
            for j in range(i + 1, nr_of_electrodes + 3):
                if reading[j] == 0:
                    M = i - 2
                    N = j - 2
                    all_active = ((reading[i] == 0) and (reading[j] == 0))
                    # print('a', 'b', 'm, n', A, B, M, N, all_active)
                    voltages.append((A, B, M, N, all_active))
                    break

        reading_configs[nr + 1] = np.array(voltages)
    return reading_configs
1374150
def get_zmatrix(self):
        
        output = []
        outputvar = []
        for i, site in enumerate(self._mol):
            if i == 0:
                output.append("{}".format(site.specie))
            elif i == 1:
                nn = self._find_nn_pos_before_site(i)
                bondlength = self._mol.get_distance(i, nn[0])
                output.append("{} {} B{}".format(self._mol[i].specie,
                                                 nn[0] + 1, i))
                outputvar.append("B{}={:.6f}".format(i, bondlength))
            elif i == 2:
                nn = self._find_nn_pos_before_site(i)
                bondlength = self._mol.get_distance(i, nn[0])
                angle = self._mol.get_angle(i, nn[0], nn[1])
                output.append("{} {} B{} {} A{}".format(self._mol[i].specie,
                                                        nn[0] + 1, i,
                                                        nn[1] + 1, i))
                outputvar.append("B{}={:.6f}".format(i, bondlength))
                outputvar.append("A{}={:.6f}".format(i, angle))
            else:
                nn = self._find_nn_pos_before_site(i)
                bondlength = self._mol.get_distance(i, nn[0])
                angle = self._mol.get_angle(i, nn[0], nn[1])
                dih = self._mol.get_dihedral(i, nn[0], nn[1], nn[2])
                output.append("{} {} B{} {} A{} {} D{}"
                              .format(self._mol[i].specie, nn[0] + 1, i,
                                      nn[1] + 1, i, nn[2] + 1, i))
                outputvar.append("B{}={:.6f}".format(i, bondlength))
                outputvar.append("A{}={:.6f}".format(i, angle))
                outputvar.append("D{}={:.6f}".format(i, dih))
        return "\n".join(output) + "\n\n" + "\n".join(outputvar)
462192
def calcMzFromMass(mass, charge):
    
    mz = (mass + (maspy.constants.atomicMassProton * charge)) / charge
    return mz
224893
def show_cationpi(self):
        """"""
        for i, p in enumerate(self.plcomplex.pication):
            cmd.pseudoatom('ps-picat-1-%i' % i, pos=p.ring_center)
            cmd.pseudoatom('ps-picat-2-%i' % i, pos=p.charge_center)
            if p.protcharged:
                cmd.pseudoatom('Chargecenter-P', pos=p.charge_center)
                cmd.pseudoatom('Centroids-L', pos=p.ring_center)
                pilig_ids = '+'.join(map(str, p.ring_atoms))
                cmd.select('PiCatRing-L', 'PiCatRing-L or (id %s & %s)' % (pilig_ids, self.ligname))
                for a in p.charge_atoms:
                    cmd.select('PosCharge-P', 'PosCharge-P or (id %i & %s)' % (a, self.protname))
            else:
                cmd.pseudoatom('Chargecenter-L', pos=p.charge_center)
                cmd.pseudoatom('Centroids-P', pos=p.ring_center)
                pires_ids = '+'.join(map(str, p.ring_atoms))
                cmd.select('PiCatRing-P', 'PiCatRing-P or (id %s & %s)' % (pires_ids, self.protname))
                for a in p.charge_atoms:
                    cmd.select('PosCharge-L', 'PosCharge-L or (id %i & %s)' % (a, self.ligname))
            cmd.distance('PiCation', 'ps-picat-1-%i' % i, 'ps-picat-2-%i' % i)
        if self.object_exists('PiCation'):
            cmd.set('dash_color', 'orange', 'PiCation')
            cmd.set('dash_gap', 0.3, 'PiCation')
            cmd.set('dash_length', 0.6, 'PiCation')
666667
def get_band_structure(self):
        """"""

        kpath = HighSymmKpath(self.data.structure)
        kpt_line = [Kpoint(k, self.data.structure.lattice) for k
                    in
                    kpath.get_kpoints(coords_are_cartesian=False)[
                        0]]
        kpoints = np.array(
            [kp.frac_coords for kp in kpt_line])

        labels_dict = {l: k for k, l in zip(
            *kpath.get_kpoints(coords_are_cartesian=False)) if l}

        lattvec = self.data.get_lattvec()
        egrid, vgrid = fite.getBands(kpoints, self.equivalences, lattvec, self.coeffs)

        bands_dict = {Spin.up: (egrid / units.eV)}

        sbs = BandStructureSymmLine(kpoints, bands_dict,
                                    self.data.structure.lattice.reciprocal_lattice,
                                    self.efermi / units.eV,
                                    labels_dict=labels_dict)
        return sbs
42526
def _parse_forces(line, lines):
    """"""
    units = line.split()[4].rstrip(":")
    next(lines)
    newline = next(lines)
    total = []; non_local = []; ionic = []; local = []; core_correction = []
    hubbard = []; scf = []; types = []
    while (not "The non-local contrib." in newline) and len(newline.split()) > 0:
        if "=" in newline:
            total.append([float(x) for x in newline.partition("=")[2].split()])
            types.append(int(newline.split()[3]))
        newline = next(lines)

    if len(newline.split()) > 0:
        while not "The ionic contribution" in newline:
            if "=" in newline:
                non_local.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while not "The local contribution" in newline:
            if "=" in newline:
                ionic.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while not "The core correction contribution" in newline:
            if "=" in newline:
                local.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while not "The Hubbard contrib." in newline:
            if "=" in newline:
                core_correction.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while not "The SCF correction term" in newline:
            if "=" in newline:
                hubbard.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while len(newline.split()) > 0:
            if "=" in newline:
                scf.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
    newline = next(lines)
    total_force = float(newline.split()[3])
    total_scf = float(newline.split()[8])
    return {
        "force units": units,
        "total force": total_force,
        "total SCF correction": total_scf,
        "forces": total,
        "non-local contribution to forces": non_local,
        "ionic contribution to forces": ionic,
        "local contribution to forces": local,
        "core corrections to forces": core_correction,
        "Hubbard contribution to forces": hubbard,
        "SCF correction term to forces": scf,
        "Atomic species index for forces": types
    }
763259
def cal_pst(self, v):
        
        p_h = self._hugoniot_p(v)
        p_th_h = self._hugoniot_pth(v)
        p_st = p_h - p_th_h
        return p_st
424402
def term_symbols(self):
        
        L_symbols = 'SPDFGHIKLMNOQRTUVWXYZ'

        L, v_e = self.valence

        # for one electron in subshell L
        ml = list(range(-L, L + 1))
        ms = [1 / 2, -1 / 2]
        # all possible configurations of ml,ms for one e in subshell L
        ml_ms = list(product(ml, ms))

        # Number of possible configurations for r electrons in subshell L.
        n = (2 * L + 1) * 2
        # the combination of n_e electrons configurations
        # C^{n}_{n_e}
        e_config_combs = list(combinations(range(n), v_e))

        # Total ML = sum(ml1, ml2), Total MS = sum(ms1, ms2)
        TL = [sum([ml_ms[comb[e]][0] for e in range(v_e)])
              for comb in e_config_combs]
        TS = [sum([ml_ms[comb[e]][1] for e in range(v_e)])
              for comb in e_config_combs]
        comb_counter = Counter([r for r in zip(TL, TS)])

        term_symbols = []
        while sum(comb_counter.values()) > 0:
            # Start from the lowest freq combination,
            # which corresponds to largest abs(L) and smallest abs(S)
            L, S = min(comb_counter)

            J = list(np.arange(abs(L - S), abs(L) + abs(S) + 1))
            term_symbols.append([str(int(2 * (abs(S)) + 1))
                                 + L_symbols[abs(L)]
                                 + str(j) for j in J])
            # Without J
            # term_symbols.append(str(int(2 * (abs(S)) + 1)) \
            #                     + L_symbols[abs(L)])

            # Delete all configurations included in this term
            for ML in range(-L, L - 1, -1):
                for MS in np.arange(S, -S + 1, 1):
                    if (ML, MS) in comb_counter:

                        comb_counter[(ML, MS)] -= 1
                        if comb_counter[(ML, MS)] == 0:
                            del comb_counter[(ML, MS)]
        return term_symbols
107280
def as_dict(self):
        
        d = {"vasp_version": self.vasp_version,
             "has_vasp_completed": True,
             "nsites": len(self.final_structure)}
        comp = self.final_structure.composition
        d["unit_cell_formula"] = comp.as_dict()
        d["reduced_cell_formula"] = Composition(comp.reduced_formula).as_dict()
        d["pretty_formula"] = comp.reduced_formula
        symbols = [s.split()[1] for s in self.potcar_symbols]
        symbols = [re.split(r"_", s)[0] for s in symbols]
        d["is_hubbard"] = self.is_hubbard
        d["hubbards"] = self.hubbards

        unique_symbols = sorted(list(set(self.atomic_symbols)))
        d["elements"] = unique_symbols
        d["nelements"] = len(unique_symbols)

        d["run_type"] = self.run_type

        vin = {"incar": {k: v for k, v in self.incar.items()},
               "crystal": self.final_structure.as_dict(),
               "kpoints": self.kpoints.as_dict()}
        actual_kpts = [{"abc": list(self.actual_kpoints[i]),
                        "weight": self.actual_kpoints_weights[i]}
                       for i in range(len(self.actual_kpoints))]
        vin["kpoints"]["actual_points"] = actual_kpts
        vin["potcar"] = [s.split(" ")[1] for s in self.potcar_symbols]
        vin["potcar_spec"] = self.potcar_spec
        vin["potcar_type"] = [s.split(" ")[0] for s in self.potcar_symbols]
        vin["parameters"] = {k: v for k, v in self.parameters.items()}
        vin["lattice_rec"] = self.lattice_rec.as_dict()
        d["input"] = vin

        vout = {"crystal": self.final_structure.as_dict(),
                "efermi": self.efermi}

        if self.eigenvalues:
            eigen = defaultdict(dict)
            for spin, values in self.eigenvalues.items():
                for i, v in enumerate(values):
                    eigen[i][str(spin)] = v
            vout["eigenvalues"] = eigen
            (gap, cbm, vbm, is_direct) = self.eigenvalue_band_properties
            vout.update(dict(bandgap=gap, cbm=cbm, vbm=vbm,
                             is_gap_direct=is_direct))

            if self.projected_eigenvalues:
                peigen = []
                for i in range(len(eigen)):
                    peigen.append({})
                for spin, v in self.projected_eigenvalues.items():
                    for kpoint_index, vv in enumerate(v):
                        if str(spin) not in peigen[kpoint_index]:
                            peigen[kpoint_index][str(spin)] = vv
                vout['projected_eigenvalues'] = peigen

        d['output'] = vout
        return jsanitize(d, strict=True)
1324510
def show_cationpi(self):
        """"""
        grp = self.getPseudoBondGroup("Cation-Pi-%i" % self.tid, associateWith=[self.model])
        grp.lineWidth = 3
        grp.lineType = self.chimera.Dash
        for i, cat in enumerate(self.plcomplex.pication):

            m = self.model
            r = m.newResidue("pseudoatoms", " ", 1, " ")
            chargecenter = m.newAtom("CHARGE", self.chimera.Element("CHARGE"))
            x, y, z = cat.charge_center
            chargecenter.setCoord(self.chimera.Coord(x, y, z))
            r.addAtom(chargecenter)

            centroid = m.newAtom("CENTROID", self.chimera.Element("CENTROID"))
            x, y, z = cat.ring_center
            centroid.setCoord(self.chimera.Coord(x, y, z))
            r.addAtom(centroid)

            b = grp.newPseudoBond(centroid, chargecenter)
            b.color = self.colorbyname('orange')

            if cat.protcharged:
                self.bs_res_ids += cat.charge_atoms
            else:
                self.bs_res_ids += cat.ring_atoms
644943
def final_energy(self):
        
        try:
            final_istep = self.ionic_steps[-1]
            if final_istep["e_wo_entrp"] != final_istep[
                'electronic_steps'][-1]["e_0_energy"]:
                warnings.warn("Final e_wo_entrp differs from the final "
                              "electronic step. VASP may have included some "
                              "corrections, e.g., vdw. Vasprun will return "
                              "the final e_wo_entrp, i.e., including "
                              "corrections in such instances.")
                return final_istep["e_wo_entrp"]
            return final_istep['electronic_steps'][-1]["e_0_energy"]
        except (IndexError, KeyError):
            warnings.warn("Calculation does not have a total energy. "
                          "Possibly a GW or similar kind of run. A value of "
                          "infinity is returned.")
            return float('inf')