24小时热门版块排行榜

返回列表

本帖产生 1 个 QC强帖，点击这里进行查看

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

yjcmwgk

禁虫 (文坛精英)

密度泛函·小卒

QC强帖: 38
应助: 165 (高中生)
贵宾: 16.999
金币: 553.8
散金: 39880
红花: 395
沙发: 2
帖子: 10406
在线: 1528.1小时
虫号: 448368
注册: 2007-11-01
专业: 理论和计算化学
管辖: 科研工具资源

[交流] 【小卒个人文集】用NLO-Calculator来计算二阶非线性光学响应系数张量及其各向异性已有23人参与

用NLO-Calculator来计算二阶非线性光学响应系数张量及其各向异性
密度泛函小卒

通过Gaussian 09的CPDFT或CPHF计算超极化率，进而再用NLO-Calculator来计算二阶非线性光学响应系数张量及其各向异性。这件事情已经完成了。Gaussian 09是商业程序。NLO Calculator (Version 0.2)是我自己写的程序，用的编程语言是Python 2.7。
可执行程序可以在美国化学会网站下载。网址http://pubs.acs.org/doi/suppl/10.1021/jp3079293，程序很小，zip压缩包的体积是4.9M，解压直接能用，不用安装。
论文是J. Phys. Chem. A, 2012, 116(41), pp 10249-10256。这个论文第一作者是我老婆，第二作者才是我。老婆指挥我干活儿，我只负责把老婆的要求翻译成程序而已。
程序说明书也同时附在可执行程序的文件夹里，解压就能看到。还有一个例子，也在程序文件夹里
特别鸣谢中山大学的Denis，南开大学的coolrainbow！(排名不分先后，都感谢)
源码公布如下：

CODE:

# 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.')