# This is a calculating program for second-order nonlinear
# optical parameters based on Gaussian 09 output file
# Dongdong Qi, NLO Calculator, Version 0.2
# University of Science and Technology Beijing (USTB)
# Beijing 100083, China
# E-mail: qdd-china@foxmail.com
import math
freq = {}
xxx = {}
xxy = {}
xyx = {}
xyy = {}
xxz = {}
xzx = {}
xyz = {}
xzy = {}
xzz = {}
yxx = {}
yyx = {}
yxy = {}
yyy = {}
yxz = {}
yzx = {}
yyz = {}
yzy = {}
yzz = {}
zxx = {}
zxy = {}
zyx = {}
zyy = {}
zxz = {}
zzx = {}
zyz = {}
zzy = {}
zzz = {}
x={}
y={}
z={}
tot={}
R1={}
R2={}
R3={}
R4={}
R5={}
R6={}
R7={}
R8={}
R9={}
R10={}
R11={}
beta2zzz={}
beta2xzz={}
betaHRS={}
DR={}
rho={}
betaJ1={}
betaJ3={}
phiJ1={}
phiJ3={}
vec={}
Bers={}
IHV={}
IHVjx={}
IHVjy={}

def Calc(file,line):
    # define global var
    global freq, xxx, xxy, xyx, xyy, xxz, xzx, xyz, xzy, xzz, yxx, yyx, yxy, \
           yyy, yxz, yzx, yyz, yzy, yzz, zxx, zxy, zyx, zyy, zxz, zzx, zyz, zzy, zzz
    global i, x, y, z, tot, R1, R2, R3, R4, R5, R6, R7, R8, R9, R10, R11, beta2zzz, \
           beta2xzz, betaHRS, DR, rho, betaJ1, betaJ3, phiJ3, phiJ1, ux, uy, uz, vec, \
           Bers, IHV, IHVjx, IHVjy
    # read log file
    freq[i] = float(line[49:60])
    file.readline()
    GL1 = file.readline() #xxx Line in Gausian output file
    GL2 = file.readline() #xxy/xyx Line in Gausian output file
    GL3 = file.readline() #xyy Line in Gausian output file
    GL4 = file.readline() #xxz/xzx Line in Gausian output file
    GL5 = file.readline() #xyz/xzy Line in Gausian output file
    GL6 = file.readline() #xzz Line in Gausian output file
    GL7 = file.readline() #yxx Line in Gausian output file
    GL8 = file.readline() #yyx/yxy Line in Gausian output file
    GL9 = file.readline() #yyy Line in Gausian output file
    GL10 = file.readline() #yxz/yzx Line in Gausian output file
    GL11 = file.readline() #yyz/yzy Line in Gausian output file
    GL12 = file.readline() #yzz Line in Gausian output file
    GL13 = file.readline() #zxx Line in Gausian output file
    GL14 = file.readline() #zxy/zyx Line in Gausian output file
    GL15 = file.readline() #zyy Line in Gausian output file
    GL16 = file.readline() #zxz/zzx Line in Gausian output file
    GL17 = file.readline() #zyz/zzy Line in Gausian output file
    GL18 = file.readline() #zzz Line in Gausian output file
    xxx[i] = -float(GL1[8:18])*math.pow(10,int(GL1[19:22]))
    xxy[i] = -float(GL2[8:18])*math.pow(10,int(GL2[19:22]))
    xyy[i] = -float(GL3[8:18])*math.pow(10,int(GL3[19:22]))
    xxz[i] = -float(GL4[8:18])*math.pow(10,int(GL4[19:22]))
    xyz[i] = -float(GL5[8:18])*math.pow(10,int(GL5[19:22]))
    xzz[i] = -float(GL6[8:18])*math.pow(10,int(GL6[19:22]))
    yxx[i] = -float(GL7[8:18])*math.pow(10,int(GL7[19:22]))
    yyx[i] = -float(GL8[8:18])*math.pow(10,int(GL8[19:22]))
    yyy[i] = -float(GL9[8:18])*math.pow(10,int(GL9[19:22]))
    yxz[i] = -float(GL10[8:18])*math.pow(10,int(GL10[19:22]))
    yyz[i] = -float(GL11[8:18])*math.pow(10,int(GL11[19:22]))
    yzz[i] = -float(GL12[8:18])*math.pow(10,int(GL12[19:22]))
    zxx[i] = -float(GL13[8:18])*math.pow(10,int(GL13[19:22]))
    zxy[i] = -float(GL14[8:18])*math.pow(10,int(GL14[19:22]))
    zyy[i] = -float(GL15[8:18])*math.pow(10,int(GL15[19:22]))
    zxz[i] = -float(GL16[8:18])*math.pow(10,int(GL16[19:22]))
    zyz[i] = -float(GL17[8:18])*math.pow(10,int(GL17[19:22]))
    zzz[i] = -float(GL18[8:18])*math.pow(10,int(GL18[19:22]))
    xyx[i] = -float(GL2[8:18])*math.pow(10,int(GL2[19:22]))
    xzx[i] = -float(GL4[8:18])*math.pow(10,int(GL4[19:22]))
    xzy[i] = -float(GL5[8:18])*math.pow(10,int(GL5[19:22]))
    yxy[i] = -float(GL8[8:18])*math.pow(10,int(GL8[19:22]))
    yzx[i] = -float(GL10[8:18])*math.pow(10,int(GL10[19:22]))
    yzy[i] = -float(GL11[8:18])*math.pow(10,int(GL11[19:22]))
    zyx[i] = -float(GL14[8:18])*math.pow(10,int(GL14[19:22]))
    zzx[i] = -float(GL16[8:18])*math.pow(10,int(GL16[19:22]))
    zzy[i] = -float(GL17[8:18])*math.pow(10,int(GL17[19:22]))
    # print Gaussian beta data
    print "\n============= Freq", i, "============================================="
    if i == 0 : print 'Freq:', freq[i], 'a.u.\nbeta( lambda = inf , Static)'
    else : print 'Freq:', freq[i], 'a.u. (', int(round(freq[i]*219474.6307)), \
         'cm-1 )\nbeta( lambda =', int(round(10000000.0/(freq[i]*219474.6307))), \
         'nm , Dynamic)'
    print "------------ Input -----------------------------------------------"
    print "beta(xxx) =", int(round(xxx[i])), "; beta(xxy) =", int(round(xxy[i])), \
          "; beta(xxz) =", int(round(xxz[i]))
    print "beta(xyx) =", int(round(xyx[i])), "; beta(xyy) =", int(round(xyy[i])), \
          "; beta(xyz) =", int(round(xyz[i]))
    print "beta(xzx) =", int(round(xzx[i])), "; beta(xzy) =", int(round(xzy[i])), \
          "; beta(xzz) =", int(round(xzz[i]))
    print "beta(yxx) =", int(round(yxx[i])), "; beta(yxy) =", int(round(yxy[i])), \
          "; beta(yxz) =", int(round(yxz[i]))
    print "beta(yyx) =", int(round(yyx[i])), "; beta(yyy) =", int(round(yyy[i])), \
          "; beta(yyz) =", int(round(yyz[i]))
    print "beta(yzx) =", int(round(yzx[i])), "; beta(yzy) =", int(round(yzy[i])), \
          "; beta(yzz) =", int(round(yzz[i]))
    print "beta(zxx) =", int(round(zxx[i])), "; beta(zxy) =", int(round(zxy[i])), \
          "; beta(zxz) =", int(round(zxz[i]))
    print "beta(zyx) =", int(round(zyx[i])), "; beta(zyy) =", int(round(zyy[i])), \
          "; beta(zyz) =", int(round(zyz[i]))
    print "beta(zzx) =", int(round(zzx[i])), "; beta(zzy) =", int(round(zzy[i])), \
          "; beta(zzz) =", int(round(zzz[i]))
    print "------------ Output ----------------------------------------------"
    # Start calculation
    # Step 1: beta(total)
    x[i] = ((xyy[i] + xzz[i] + yxy[i] + zxz[i] + yyx[i] + zzx[i]) / 3.0) +  xxx[i]
    y[i] = ((yxx[i] + yzz[i] + xyx[i] + zyz[i] + xxy[i] + zzy[i]) / 3.0) +  yyy[i]
    z[i] = ((zxx[i] + zyy[i] + xzx[i] + yzy[i] + xxz[i] + yyz[i]) / 3.0) +  zzz[i]
    print 'beta(x) =', int(round(x[i])), '; beta(y) =', int(round(y[i])), \
          '; beta(z) =', int(round(z[i]))
    tot[i] = math.sqrt(math.pow(x[i],2) + math.pow(y[i],2) + math.pow(z[i],2))
    print 'beta(total) =', int(round(tot[i]))
    # step 2: and
    R1[i] = (math.pow(xxx[i],2) + math.pow(yyy[i],2) + math.pow(zzz[i],2))
    R2[i] = (math.pow(xxy[i],2) + math.pow(xxz[i],2) + math.pow(yyx[i],2) + \
             math.pow(yyz[i],2) + math.pow(zzx[i],2) + math.pow(zzy[i],2))
    R3[i] = (xxx[i]*xyy[i] + xxx[i]*xzz[i] + yyy[i]*yxx[i] + yyy[i]*yzz[i] + \
             zzz[i]*zxx[i] + zzz[i]*zyy[i])
    R4[i] = (xyy[i]*yyx[i] + xzz[i]*zzx[i] + yxx[i]*xxy[i] + yzz[i]*zzy[i] + \
             zxx[i]*xxz[i] + zyy[i]*yyz[i])
    R5[i] = (xxx[i]*yyx[i] + xxx[i]*zzx[i] + yyy[i]*xxy[i] + yyy[i]*zzy[i] + \
             zzz[i]*xxz[i] + zzz[i]*yyz[i])
    R6[i] = (math.pow(yxx[i],2) + math.pow(zxx[i],2) + math.pow(xyy[i],2) + \
             math.pow(zyy[i],2) + math.pow(xzz[i],2) + math.pow(yzz[i],2))
    R7[i] = (xxy[i]*yzz[i] + xxz[i]*zyy[i] + yyx[i]*xzz[i] + yyz[i]*zxx[i] + \
             zzx[i]*xyy[i] + zzy[i]*yxx[i])
    R8[i] = (xyy[i]*xzz[i] + yxx[i]*yzz[i] + zxx[i]*zyy[i]) * 2.0
    R9[i] = (xxy[i]*zzy[i] + xxz[i]*yyz[i] + yyx[i]*zzx[i]) * 2.0
    R10[i] = (math.pow(xyz[i],2) + math.pow(xzy[i],2) + math.pow(yxz[i],2) + \
              math.pow(yzx[i],2) + math.pow(zyx[i],2) + math.pow(zxy[i],2))
    R11[i] = (xyz[i]*yxz[i] + xzy[i]*zxy[i] + yzx[i]*zyx[i]) * 2.0
    beta2zzz[i] = (R1[i]*15.0 + R2[i]*12.0 + R3[i]*6.0 + R4[i]*12.0 + \
                   R5[i]*12.0 + R6[i]*3.0 + R7[i]*4.0 + R8[i] + R9[i]*4.0 + \
                   R10[i]*2.0 + R11[i]*4.0)/105.0
    beta2xzz[i] = (R1[i]*3.0 + R2[i]*8.0 + R3[i]*4.0 - R4[i]*6.0 - \
                   R5[i]*6.0 + R6[i]*9.0 - R7[i]*2.0 + R8[i]*3.0 - \
                   R9[i]*2.0 + R10[i]*6.0 - R11[i]*2.0)/105.0
    betaHRS[i] = math.sqrt(beta2zzz[i] + beta2xzz[i])
    DR[i] = beta2zzz[i]/beta2xzz[i]
    if DR[i] < 1.5 :
        print "Be careful! DR is lower than 1.5."
    print " =", int(round(beta2zzz[i])), "; =", \
          int(round(beta2xzz[i]))
    print "beta(HRS)(-2w;w,w) =", int(round(betaHRS[i])), "; DR =", DR[i]
    if DR[i] > 1.5 :
        rho[i] = math.sqrt(7.0*(9.0-DR[i])/(12.0*DR[i]-18.0))
        print "rho =", rho[i]
        betaJ1[i] = betaHRS[i]/math.sqrt((2.0/9.0)+(2.0*rho[i]*rho[i]/21.0))
        betaJ3[i] = rho[i]*betaJ1[i]
        phiJ3[i] = rho[i]/(1.0+rho[i])
        phiJ1[i] = 1.0-phiJ3[i]
        print "|beta(J=1)| =", betaJ1[i], "; |beta(J=3)| =", betaJ3[i]
        print "phi[beta(J=1)] =", phiJ1[i], "; phi[beta(J=3)] =", phiJ3[i]
    else :
        print "\n      ____ Attention ___________________________________\n"
        print "       If DR is lower than 1.5, the calculations of"
        print "       rho, |beta(J=1)|, |beta(J=3)|, phi[beta(J=1)],"
        print "       and phi[beta(J=3)] are unreliable. In addition,"
        print "       please check whether the half of the incident"
        print "       light freqency (lambda/2) is near to the"
        print "       absorption band, which is the direct reason of"
        print "       the resonance for SFG. For details, please"
        print "       check any NLO textbook."
        print "      __________________________________________________\n"
    # step 3: IHV list
    Bers[i] = 7.0 * beta2xzz[i] - beta2zzz[i]
    angle = -180
    while (angle < 181):
        IHV[i,angle] = (math.pow(math.cos(float(angle)*(2.0*math.pi)/360.0),4))*beta2xzz[i] + \
                       (math.pow(math.sin(float(angle)*(2.0*math.pi)/360.0),4))*beta2zzz[i] + \
                       (math.pow(math.cos(float(angle)*(2.0*math.pi)/360.0),2))*\
                       (math.pow(math.sin(float(angle)*(2.0*math.pi)/360.0),2))*Bers[i]
        IHVjx[i,angle] = IHV[i,angle] * math.cos(float(angle)*(2.0*math.pi)/360.0)
        IHVjy[i,angle] = IHV[i,angle] * math.sin(float(angle)*(2.0*math.pi)/360.0)
        angle = angle + 1
    print "If you need the I(2w,PsiV) values, please check them in I2wPsiV.txt."
    i=i+1
   
# Main Programe
print "       __________________ NLO Calculator _________________\n"
print "                       v 0.2 ( June 21, 2012 )\n"
print "       Theory: [1] J. Chem. Phys. 2012, 136, 024506."
print "               [2] Phys. Chem. Chem. Phys. 2008, 10, 6223.\n"
print "Cite:   Dongdong Qi, NLO Calculator, Version 0.2"
print "        University of Science and Technology Beijing (USTB)"
print "        Beijing 100083, China\nE-mail: qdd-china@foxmail.com"
filename = raw_input('Please input the name of Gaussian log/out file: ')
file = open(filename, 'r')
line = None
i = 0
freqN = input('How many special frequencies in this log/out file? ')
while (line != "EOF") and (i < freqN+1) :
    line = file.readline()
    if line.startswith(' Property number 6 -- Beta(w,w,-2w) frequency'):
        Calc(file,line)
file.close()
f = open('I2wPsiV.txt','w')
f.write("       __________________ NLO Calculator _________________\n\n")
f.write("                       v 0.2 ( June 21, 2012 )\n")
f.write("       Theory: [1] J. Chem. Phys. 2012, 136, 024506.\n")
f.write("               [2] Phys. Chem. Chem. Phys. 2008, 10, 6223.\n")
f.write("       ___________________________________________________\n")
f.write("\n___________________________________________________________________\n")
f.write("\nCite:   Dongdong Qi, NLO Calculator, Version 0.2")
f.write("        University of Science and Technology Beijing (USTB)")
f.write("        Beijing 100083, China\nE-mail: qdd-china@foxmail.com")
f.write("___________________________________________________________________\n\n")
i = -1
while (i < freqN):
    i = i+1
    f.write("Freq ")
    f.write(repr(i))
    f.write(": ")
    f.write(repr(freq[i]))
    f.write(" a.u.\nlambda = ")
    f.write(repr(int(round(freq[i]*219474.6307))))
    f.write(" cm-1 \n")
    if (i!=0): f.write(repr(int(round(10000000.0/(freq[i]*219474.6307)))))
    else: f.write("inf")
    f.write(" nm\n\n")
    if (i!=0): f.write("Dynamic beta(HRS)(-2w;w,w): ")
    else: f.write("Static beta(HRS)(-2w;w,w): ")
    f.write(repr(int(round(betaHRS[i]))))
    f.write("\n\n")
    anglex = -180
    while (anglex < 181):
        f.write("Psi = ")
        f.write(repr(anglex))
        f.write(" ; [Cartesian Coord.] I(2w,PsiV) = ")
        f.write(repr(int(round(IHV[i,anglex]))))
        f.write(" ; [Polar Coord.] I(2w,PsiV,x) = ")
        f.write(repr(int(round(IHVjx[i,anglex]))))
        f.write(" and I(2w,PsiV,y) = ")
        f.write(repr(int(round(IHVjy[i,anglex]))))
        if (anglex == 0): f.write(" ===== Horizontal [I(2w,HV)] =====")
        elif (anglex == 90): f.write(" ===== Vertical [I(2w,VV)] =====")
        elif (anglex == 45): f.write(" ===== Right-circular =====")
        elif (anglex == -45): f.write(" ===== Left-circular =====")
        f.write("\n")
        anglex = anglex + 1
    f.write("________________________________________\n\n")
f.close()
raw_input('Press to exit.')

3樓: Originally posted by 諸葛神鷹 at 2012-11-18 22:19:20
好牛啊，那天要計算么，找你教教怎么使用，呵呵