Python scipy.constants.Boltzmann() Examples

def _ThCritical(self, rho, T, fase):
        """Younglove-Hanley thermal conductivity critical enhancement"""
        rhom = rho/self.M
        Tc = 150.86
        Pc = 4.9058e6
        rhocm = 13.41
        rhoc = rhocm*self.M

        # Eq 9
        T_ = T/Tc
        rho_ = rhom/rhocm
        DelT = (T-Tc)/Tc
        Delrho = (rhom-rhocm)/rhocm
        Xt = rho*fase.drhodP_T*Pc/rhoc**2

        # Eq 10
        tc = 1.02*Boltzmann*Pc/6.0795e-1/6/pi/fase.mu*(T_/rho_)**2 * \
            fase.dpdT_rho**2*Xt**0.46807*exp(-39.8*DelT**2-5.45*Delrho**4)

        return tc 

def _check_incar_ele_temp(testCase, idx, ele_temp):
    fp_path = os.path.join('iter.%06d' % idx, '02.fp')
    tasks = glob.glob(os.path.join(fp_path, 'task.*'))
    cwd = os.getcwd()
    for ii in tasks :
        os.chdir(ii)
        bname = os.path.basename(ii)
        sidx = int(bname.split('.')[1])
        tidx = int(bname.split('.')[2])
        with open('INCAR') as fp:
            incar = fp.read()
            incar0 = Incar.from_string(incar)
            # make_fake_md: the frames in a system shares the same ele_temp
            incar1 = Incar.from_string(vasp_incar_ele_temp_ref%(ele_temp[sidx][0] * pc.Boltzmann / pc.electron_volt))
            for ii in incar0.keys():
                # skip checking nbands...
                if ii == 'NBANDS':
                    continue
                testCase.assertAlmostEqual(incar0[ii], incar1[ii], msg = 'key %s differ' % (ii), places = 5)
        os.chdir(cwd) 

def _thermo0(self, rho, T, fase):
        """Custom dilute-gas limit form of thermal conductivity"""
        X1 = 0.95185202
        X2 = 1.0205422
        M = 28.013

        # Used the ideal isochoric heat capacity of paper because differ in
        # about a 20% of values using the ideal correlation in meos
        def cv(T):
            ni = [-0.837079888737e3, 0.379147114487e2, -0.601737844275,
                  0.350418363823e1, -0.874955653028e-5, 0.148968607239e-7,
                  -0.256370354277e-11, 0.100773735767e1, 0.335340610e4]
            sum1 = 0
            for i, n in enumerate(ni[:-2]):
                sum1 += n*T**(i-3)

            u = ni[8]/T
            eu1 = exp(u)-1
            sum2 = (ni[7]*u**2*(eu1+1))/eu1**2
            cv = 8.31434*(sum1+sum2-1)
            return cv

        muo = self._Visco0(T, self.visco2)
        F = Boltzmann*Avogadro*muo/M                                     # Eq 9
        ltr = 2.5*F*(1.5-X1)                                             # Eq 7
        lint = F*X2*(cv(T)/Boltzmann/Avogadro+X1)                        # Eq 8

        return (ltr+lint)*1e-3 

def _ThCondCritical(self, rho, T, fase):
        # Custom Critical enhancement

        # The paper use a diferent rhoc value to the EoS
        rhoc = 235

        t = abs(T-405.4)/405.4
        dPT = 1e5*(2.18-0.12/exp(17.8*t))
        nb = 1e-5*(2.6+1.6*t)

        DL = 1.2*Boltzmann*T**2/6/pi/nb/(1.34e-10/t**0.63*(1+t**0.5)) * \
            dPT**2 * 0.423e-8/t**1.24*(1+t**0.5/0.7)

        # Add correction for entire range of temperature, Eq 10
        DL *= exp(-36*t**2)

        X = 0.61*rhoc+16.5*log(t)
        if rho > 0.6*rhoc:
            # Eq 11
            DL *= X**2/(X**2+(rho-0.96*rhoc)**2)
        else:
            # Eq 14
            DL = X**2/(X**2+(0.6*rhoc-0.96*rhoc)**2)
            DL *= rho**2/(0.6*rhoc)**2

        return DL 

def _thermo0(self, rho, T, fase):
        # λ1 in Eq 3 is always 0

        GT = [-2.903423528e5, 4.680624952e5, -1.8954783215e5, -4.8262235392e3,
              2.243409372e4, -6.6206354818e3, 8.9937717078e2, -6.0559143718e1,
              1.6370306422]
        lo = 0
        for i in range(-3, 6):
            lo += GT[i+3]*T**(i/3.)

        l2, lc = 0, 0
        if rho:
            tita = (rho-221)/221
            k = [-1.304503323e1, 1.8214616599e1, -9.903022496e3, 7.420521631e2,
                 -3.0083271933e-1, 9.6456068829e1, 1.350256962e4]
            l2 = exp(k[0]+k[3]/T) * (
                exp(rho.gcc**0.1*(k[1]+k[2]/T**1.5) +
                    tita*rho.gcc**0.5*(k[4]+k[5]/T+k[6]/T**2))-1)

            # Critical enhancement
            deltarho = (rho-221)/221
            deltaT = (T-282.34)/282.34

            xt = rho**2*fase.kappa*5.039/221**2
            B = abs(deltarho)/abs(deltaT)**1.19                         # Eq 11
            Gamma = xt*abs(deltaT)**1.19                                # Eq 12
            xi = 0.69/(B**2*5.039/Gamma/Boltzmann/282.34)**0.5          # Eq 14

            # Eq 19
            F = exp(-18.66*deltaT**2) * exp(-4.25*deltarho**4)

            # Eq 18
            c = (self.M/rho.gcc/Avogadro/Boltzmann/T)**0.5
            d = Boltzmann*T**2/6/pi/fase.mu.muPas/xi
            lc = c*d*fase.dpdT_rho**2*fase.kappa**0.5*F

        return unidades.ThermalConductivity(lo+l2+lc, "mWmK") 

def Tension_Miqueu(T, Tc, Vc, M, w):
    r"""Calculates surface tension of a liquid using the Miqueu et al.
    correlation

    .. math::
        \sigma = kT_c\left(\frac{N_A}{V_c}\right)^{2/3}
        (4.35+4.14\omega)t^{1.26}(1+0.19t^{0.5}-0.487t)

    Parameters
    ----------
    T : float
        Temperature, [K]
    Tc : float
        Critical temperature, [K]
    Vc : float
        Critical volume, [m^3/kg]
    M : float
        Molecular weight, [g/mol]
    w : float
        Acentric factor, [-]

    Returns
    -------
    sigma : float
        Liquid surface tension, [N/m]
    """
    t = 1 - T/Tc

    # Eq 13
    sigma = Boltzmann * Tc * (Avogadro/Vc/1000/M)**(2/3) * (4.35+4.14*w) * \
        t**1.26 * (1+0.19*t**0.5-0.25*t)
    return unidades.Tension(sigma, "mNm") 

def make_vasp_incar_ele_temp(jdata, filename, ele_temp, nbands_esti = None):
    with open(filename) as fp:
        incar = fp.read()
    incar = incar_upper(Incar.from_string(incar))
    incar['ISMEAR'] = -1
    incar['SIGMA'] = ele_temp * pc.Boltzmann / pc.electron_volt
    incar.write_file('INCAR')
    if nbands_esti is not None:
        nbands = nbands_esti.predict('.')
        with open(filename) as fp:
            incar = Incar.from_string(fp.read())
        incar['NBANDS'] = nbands
        incar.write_file('INCAR') 

def _get_res(self, res_dir):
        res = {}
        sys = dpdata.System(os.path.join(res_dir, 'POSCAR'))
        res['natoms'] = (sys['atom_numbs'])
        res['vol'] = np.linalg.det(sys['cells'][0])
        res['nvalence'] = (self._get_potcar_nvalence(os.path.join(res_dir, 'POTCAR')))
        res['ele_temp'] = self._get_incar_ele_temp(os.path.join(res_dir, 'INCAR')) * pc.electron_volt / pc.Boltzmann
        res['nbands'] = self._get_incar_nbands(os.path.join(res_dir, 'INCAR'))
        return res