Chemistry 版 (精华区)

发信人: zjliu (秋天的萝卜), 信区: Chemistry
标  题: 配分函数计算器 V0.2 BY boliwa
发信站: BBS 哈工大紫丁香站 (Thu Jul  8 16:32:25 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
    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 "C:\WINDOWS\DESKTOP\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)
    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, RotTemp 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
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 *
Boltzmann * Temp) ^ 3 / Planck ^ 3
    ELSEIF LEFT$(GetStr, 43) = " Sum of electronic and zero-point
Energies=" THEN
        Epoint = VAL(MID$(GetStr, 50, 20))
    ELSEIF UCASE$(LEFT$(GetStr, 27)) = " ROTATIONAL SYMMETRY NUMBER" THEN
        RotSymm = VAL(MID$(GetStr, 29, 2))
    ELSEIF LEFT$(GetStr, 33) = " ROTATIONAL TEMPERATURES (KELVIN)" OR
LEFT$(GetStr, 32) = " Rotational temperature (Kelvin)" THEN
        RotTemp.x = VAL(MID$(GetStr, 34, 11))
        RotTemp.y = VAL(MID$(GetStr, 48, 11))
        IF RotTemp.y = 0 THEN
            RotTemp.y = RotTemp.x
        END IF
        RotTemp.z = VAL(MID$(GetStr, 60, 11))
        IF RotTemp.z = 0 THEN
            RotTemp.z = 1
            RotTemp.z = 1
        END IF
        qr = Temp / RotSymm / SQR(RotTemp.x * RotTemp.y * RotTemp.z)
    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
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
TotQ = qt * qv * qr
TotQ2 = qt * qv2 * qr

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


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