CArm code
- make a file called CArm.py and copy the following into it
- this code needs some editing to make it general!!
- it reads data from cube files, nbo analysis, mulliken analysis ...
#!/usr/bin/python
import sys, re
def main():
global atoms
global charge
#print(atoms)
#print(charge)
print("=====================main=====================")
CoM=[0,0,0]
mass=0
for i in atoms:
CoM[0]=CoM[0]+i[2]*masses[i[0]]
CoM[1]=CoM[1]+i[3]*masses[i[0]]
CoM[2]=CoM[2]+i[4]*masses[i[0]]
mass=mass+masses[i[0]]
#print(mass)
CoM[0] = CoM[0] / mass
CoM[1] = CoM[1] / mass
CoM[2] = CoM[2] / mass
print "center of mass:", CoM
CoC=[0,0,0]
COCC=0
uCOCC=0
for i in charge:
#print(i)
CoC[0]=CoC[0]+i[0]*i[3]
CoC[1]=CoC[1]+i[1]*i[3]
CoC[2]=CoC[2]+i[2]*i[3]
COCC=COCC+i[3]
uCOCC=uCOCC+abs(i[3])
#print(COCC)
print "charge:", COCC
CoC[0] = (CoC[0] / COCC) - CoM[0]
CoC[1] = (CoC[1] / COCC) - CoM[0]
CoC[2] = (CoC[2] / COCC) - CoM[0]
#print(CoC)
norm=(CoC[0]**2+CoC[1]**2+CoC[2]**2)**.5*0.529177
print "length:", norm
CoC[0] = CoC[0] / norm
CoC[1] = CoC[1] / norm
CoC[2] = CoC[2] / norm
print "direction:", CoC
def cuberead():
print("===================cuberead===================")
global atoms
global charge
cube=open("OTf-_FREQ.cube",'r')
line=cube.readline()
#print(line)
line=cube.readline()
#print(line)
line=cube.readline()
items=line.split()
NAtoms=int(items[0])
origin=[0,0,0]
origin[0]=float(items[1])
origin[1]=float(items[2])
origin[2]=float(items[3])
print(origin)
line=cube.readline()
items=line.split()
Xstep=int(items[0])
Xvec=[0,0,0]
Xvec[0]=float(items[1])
Xvec[1]=float(items[2])
Xvec[2]=float(items[3])
print(Xvec)
line=cube.readline()
items=line.split()
Ystep=int(items[0])
Yvec=[0,0,0]
Yvec[0]=float(items[1])
Yvec[1]=float(items[2])
Yvec[2]=float(items[3])
print(Yvec)
line=cube.readline()
items=line.split()
Zstep=int(items[0])
Zvec=[0,0,0]
Zvec[0]=float(items[1])
Zvec[1]=float(items[2])
Zvec[2]=float(items[3])
print(Zvec)
atoms=[]
for i in range(0,NAtoms):
line=cube.readline()
items=line.split()
atoms.append([int(items[0]),float(items[1]),float(items[2]),float(items[3]),float(items[4])])
charge=[[0,0,0,0]]
CSum=0
XfaceA=0
XfaceB=0
YfaceA=0
YfaceB=0
ZfaceA=0
ZfaceB=0
for x in range(0,Xstep):
#print(x)
for y in range(0,Ystep):
z=0
while True:
line=cube.readline()
items=line.split()
for j in items:
#if x==0 :
# XfaceA=XfaceA + float(j)
#if x==Xstep-1 :
# XfaceB=XfaceB + float(j)
#if y==0 :
# YfaceA=YfaceA + float(j)
#if y==Ystep-1 :
# YfaceB=YfaceB + float(j)
#if z==0 :
# ZfaceA=ZfaceA + float(j)
#if z==Zstep-1 :
# ZfaceB=ZfaceB + float(j)
#charge.append([x*Xvec[0]+origin[0]+y*Yvec[0]+z*Zvec[0],x*Xvec[1]+y*Yvec[1]+z*Zvec[1]+origin[1],x*Xvec[2]+y*Yvec[2]+z*Zvec[2]+origin[2],float(j)])
charge=[[charge[0][0]+(x*Xvec[0]+y*Yvec[0]+z*Zvec[0]+origin[0])*float(j),charge[0][1]+(x*Xvec[1]+y*Yvec[1]+z*Zvec[1]+origin[1])*float(j),charge[0][2]+(x*Xvec[2]+y*Yvec[2]+z*Zvec[2]+origin[2])*float(j),float(j)+charge[0][3]]]
z=z+1
#CSum=CSum-float(j)*Xvec[0]*Yvec[1]*Zvec[2]
CSum=CSum-float(j)
if z == Zstep :
break
#print(x+1,'of', Xstep)
charge[0]=[charge[0][0]/abs(CSum),charge[0][1]/abs(CSum),charge[0][2]/abs(CSum),CSum*Xvec[0]*Yvec[1]*Zvec[2]]
#print(XfaceA)
#print(XfaceB)
#print(YfaceA)
#print(YfaceB)
#print(ZfaceA)
#print(ZfaceB)
#print(CSum)
#print(charge)
for i in atoms:
charge.append([i[2],i[3],i[4],i[1]])
#print(charge)
#print([i[2],i[3],i[4],float(items[2])])
#print(CSum/Xvec[0]/Yvec[1]/Zvec[2])
return()
def mullread():
print("===================mullread===================")
global atoms
global charge
mull=open("OTf-_FREQ.log",'r')
while True:
line=mull.readline()
if 'Standard orientation' in line:
#print("hurra")
atoms=[]
for i in range(0,5):
line=mull.readline()
#print(line)
while True:
if '--' in line:
break
items=line.split()
atoms.append([int(items[1]),float(items[1]),float(items[3])/0.529177,float(items[4])/0.529177,float(items[5])/0.529177])
line=mull.readline()
#print(atoms)
if "Mulliken atomic charges" in line:
charge=[]
for i in range(0,2):
line=mull.readline()
#print(line)
for i in atoms:
items=line.split()
charge.append([i[2],i[3],i[4],float(items[2])])
line=mull.readline()
#print(charge)
if line=='':
break
return()
def helpread():
print("===================helpread===================")
global atoms
global charge
atoms=[]
charge=[]
mull=open("OTf-_FREQ.log",'r')
while True:
line=mull.readline()
if 'Standard orientation' in line:
#print("hurra")
atoms=[]
for i in range(0,5):
line=mull.readline()
#print(line)
while True:
if '--' in line:
break
items=line.split()
atoms.append([int(items[1]),float(items[1]),float(items[3])/0.529177,float(items[4])/0.529177,float(items[5])/0.529177])
line=mull.readline()
#print(atoms)
if "Charges from ESP fit" in line:
charge=[]
for i in range(0,3):
line=mull.readline()
#print(line)
for i in atoms:
items=line.split()
charge.append([i[2],i[3],i[4],float(items[2])])
line=mull.readline()
#print(charge)
if line=='':
break
return()
def naturead():
print("===================naturead===================")
global atoms
global charge
mull=open("OTf-_FREQ.log",'r')
while True:
line=mull.readline()
if 'Standard orientation' in line:
#print("hurra")
atoms=[]
for i in range(0,5):
line=mull.readline()
#print(line)
while True:
if '--' in line:
break
items=line.split()
atoms.append([int(items[1]),float(items[1]),float(items[3])/0.529177,float(items[4])/0.529177,float(items[5])/0.529177])
line=mull.readline()
#print(atoms)
if "Summary of Natural Population Analysis" in line:
charge=[]
for i in range(0,6):
line=mull.readline()
#print(line)
for i in atoms:
items=line.split()
charge.append([i[2],i[3],i[4],float(items[2])])
#print([i[2],i[3],i[4],float(items[2])])
line=mull.readline()
#print(charge)
if line=='':
break
return()
masses = {7:14.00674, 16:32.066, 9:18.9984032, 1:1.00794, 6:12.0107, 14:28.0855, 17:35.4527, 8:15.9994}
cuberead()
main()
#print(atoms)
#print(charge)
mullread()
main()
#print(atoms)
#print(charge)
naturead()
main()
#print(atoms)
#print(charge)
helpread()
main()
#print(atoms)
#print(charge)
quit()
