Chemistry 版 (精华区)

发信人: zjliu (秋天的萝卜), 信区: Chemistry
标  题: 配分函数计算器V1.0 BY boliwa
发信站: BBS 哈工大紫丁香站 (Thu Jul  8 16:34:52 2004)

发信站: 日月光华  

可以用来从gaussian输出文件计算配分函数以及反应速率常数

gaussian98 03都可以使用

欢迎报告BUG
我可以尽快修改


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
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 过渡态文件名, Temp, qt(1), qv(1), qr(1),
Epoint(1), qv2(1), TotQ(1), TotQ2(1)
printq 产物1文件名, Temp, qt(2), qv(2), qr(2), Epoint(2), qv2(2),
TotQ(2), TotQ2(2)
printq 产物2文件名, Temp, qt(3), qv(3), qr(3), Epoint(3), qv2(3),
TotQ(3), TotQ2(3)
'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 "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,
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 *
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 UCASE$(LEFT$(GetStr, 20)) = UCASE$(" Rotational constant")
THEN
        RotConst.x = VAL(MID$(GetStr, 36, 11))
        IF RotConst.x <> 0 THEN
        RotConst.x = Planck / (8 * Pi * Pi * RotConst.x * 1000000000)
        ELSE
        ELSE
        RotConst.x = -1
        END IF
        RotConst.y = VAL(MID$(GetStr, 48, 11))
        IF RotConst.y <> 0 THEN
        RotConst.y = Planck / (8 * Pi * Pi * RotConst.y * 1000000000)
        ELSE
        RotConst.y = -1
        END IF
        RotConst.z = VAL(MID$(GetStr, 60, 11))
        IF RotConst.z <> 0 THEN
        RotConst.z = Planck / (8 * Pi * Pi * RotConst.z * 1000000000)
        ELSE
        RotConst.z = -1
        END IF
           PRINT RotConst.x, RotConst.y, RotConst.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
            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
IF RotConst.x > 0 AND RotConst.y > 0 THEN
    qr = 8 * Pi * Pi * SQR(2 * Pi * Boltzmann * Temp) ^ 3 * SQR(RotConst.
x * RotConst.y * RotConst.z) / RotSymm / Planck ^ 3
ELSE
    qr = 8 * Pi * Pi * RotConst.x * Boltzmann * Temp / RotSymm / Planck ^
2
END IF
qv = 1
qv2 = 1
FOR I = 1 TO VibFree
    qv = qv * 1 / (1 - EXP(-Planck * VibFreq(I) * 100 * C / Boltzmann /
Temp))
FOR I = 1 TO VibFree
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.*]