Chemistry 版 (精华区)

发信人: zjliu (秋天的萝卜), 信区: Chemistry
标  题: 配分函数计算器V0.1 BY boliwa&auntyellow
发信站: BBS 哈工大紫丁香站 (Wed Jun 30 10:16:16 2004)

发信站: 日月光华  

CONST Pi = 3.141592653589793#
CONST C = 299792458#
CONST Boltzmann = 1.380662D-23
CONST mole = 6.022045D+23
CONST Planck = 6.62616D-34
DIM Temp  AS DOUBLE, qt(1 TO 256) AS DOUBLE, qv(1 TO 256) AS DOUBLE, qr(1
 TO 256) AS DOUBLE, Epoint(1 TO 256) AS DOUBLE, qv2(1 TO 256) AS DOUBLE,
TotQ(1 TO 256) AS DOUBLE, TotQ2(1 TO 256) AS DOUBLE
DIM TotQ AS DOUBLE
DIM k AS DOUBLE
DIM E0 AS DOUBLE
TYPE Coordinate
    x AS DOUBLE
    y AS DOUBLE
    z AS DOUBLE
END TYPE
DECLARE SUB printq (FileName AS STRING, Temp AS DOUBLE, qt AS DOUBLE, qv
AS DOUBLE, qr AS DOUBLE, Epoint AS DOUBLE, qv2 AS DOUBLE, TotQ AS DOUBLE,
 TotQ2 AS DOUBLE)
Temp = 298.15#
printq "f:/OUT23.OUT", Temp, qt(1), qv(1), qr(1), Epoint(1), qv2(1), TotQ
(1), TotQ2(1)
printq "f:/C4H8.OUT ", Temp, qt(2), qv(2), qr(2), Epoint(2), qv2(2), TotQ
(2), TotQ2(2)
printq "f:/C2H4.OUT", Temp, qt(3), qv(3), qr(3), Epoint(3), qv2(3), TotQ(
3), TotQ2(3)
PRINT
FOR i = 1 TO 3
    PRINT "qt="; qt(i),
    PRINT "qv="; qv(i),
    PRINT "qr="; qr(i),
    PRINT "ep="; Epoint(i),
    PRINT "qv2="; qv2(i),
    PRINT "totq="; TotQ(i),
    PRINT "totq2="; TotQ2(i)
NEXT
E0 = (Epoint(1) - (Epoint(2) + Epoint(3))) * 27.2 * 96485
PRINT "E0="; E0
k = Boltzmann * Temp / Planck * TotQ2(1) / TotQ(2) / TotQ(3) * EXP(-E0 /
8.31 / Temp)
PRINT "k="; k

SUB printq (FileName AS STRING, Temp AS DOUBLE, qt AS DOUBLE, qv AS DOUBLE
, qr AS DOUBLE, Epoint AS DOUBLE, qv2 AS DOUBLE, TotQ AS DOUBLE, TotQ2 AS
 DOUBLE)
CONST Pi = 3.141592653589793#
CONST C = 299792458#
CONST Boltzmann = 1.3806503D-23
CONST u = 1.66053873D-27
CONST Planck = 6.62616D-34
CONST Pressure = 101325
DIM VibFreq(1 TO 256) AS DOUBLE
DIM MolMass AS DOUBLE, RotConst AS Coordinate, RotSymm AS INTEGER, VibFree
 AS INTEGER
DIM i AS INTEGER, FileSpec AS INTEGER, GetStr AS STRING
FileSpec = FREEFILE
OPEN FileName FOR INPUT AS #FileSpec
i = 0
DO WHILE NOT EOF(FileSpec)
    LINE INPUT #FileSpec, GetStr
    IF LEFT$(GetStr, 16) = " Molecular mass:" THEN
        MolMass = VAL(MID$(GetStr, 18, 11))
        qt = Boltzmann * Temp / Pressure * SQR(2 * Pi * MolMass * u * Boltzman
n
 * Temp) ^ 3 / Planck ^ 3
    ELSEIF LEFT$(GetStr, 43) = " Sum of electronic and zero-point Energies="
 THEN
        Epoint = VAL(MID$(GetStr, 50, 20))
    ELSEIF LEFT$(GetStr, 28) = " Rotational constants (GHZ):" THEN
        RotConst.x = VAL(MID$(GetStr, 36, 11))
        RotConst.x = Planck / 8 / Pi / Pi / RotConst.x * .000000001#
        RotConst.y = VAL(MID$(GetStr, 48, 11))
        RotConst.y = Planck / 8 / Pi / Pi / RotConst.y * .000000001#
        RotConst.z = VAL(MID$(GetStr, 60, 11))
        RotConst.z = Planck / 8 / Pi / Pi / RotConst.z * .000000001#
    ELSEIF LEFT$(GetStr, 27) = " ROTATIONAL SYMMETRY NUMBER" THEN
        RotSymm = VAL(MID$(GetStr, 29, 2))
        qr = 8 * Pi * Pi * SQR(2 * Pi * Boltzmann * Temp) ^ 3 * SQR(RotConst.x

* RotConst.y * RotConst.z) / RotSymm / Planck ^ 3
    ELSEIF LEFT$(GetStr, 16) = " Frequencies -- " THEN
        i = i + 1
        VibFreq(i) = VAL(MID$(GetStr, 17, 10))
        IF LEN(GetStr) > 26 THEN
            i = i + 1
            VibFreq(i) = VAL(MID$(GetStr, 40, 10))
            IF LEN(GetStr) > 49 THEN
                i = i + 1
                VibFreq(i) = VAL(MID$(GetStr, 63, 10))
            END IF
        END IF
    END IF

LOOP
VibFree = i

CLOSE #FileSpec
qv = 1
qv2 = 1
FOR i = 1 TO VibFree
    qv = qv * 1 / (1 - EXP(-Planck * VibFreq(i) * 100 * C / Boltzmann / Temp
))
NEXT
FOR i = 1 TO VibFree
    IF VibFreq(i) > 0 THEN
    qv2 = qv2 * 1 / (1 - EXP(-Planck * VibFreq(i) * 100 * C / Boltzmann / Temp

))
    END IF
NEXT
TotQ2 = qt * qv2 * qr

END SUB


--
╔═══════════════════╗
║★★★★★友谊第一  比赛第二★★★★★║
╚═══════════════════╝


※ 来源:·哈工大紫丁香 http://bbs.hit.edu.cn·[FROM: 202.118.229.*]