Mod:Hunt Research Group/dipole moment
Jump to navigation
Jump to search
code to calculate the dipole moment
- What the program does?
- 1. Calculate the dipole moment using CHELPG charges and check against the gaussian dipole moment.
- 2. Calculate the dipole moment using other charge schemes.
- 3. Calculate the dipole moment relative to the centre of positive nuclear charge (CoP) vs. the centre of mass (CoM).
- The dipole moment (μ) wrt the CoP and the μ wrt the CoM is the same for neutral molecules but differs for charged molecules.
- Only the lowest non-vanishing multipole moment remains unchanged with the choice of the origin.
- preparing the input file
- the first line of the input file must specify in the following order separated by whitespaces: the number of atoms, the total charge(with a decimal point), the gaussian calculated dipole moment and its component (ux,uy,uz), then leave a blank line
- next lines must contain the atomic number and the corresponding XYZ coordinates
- leave a blank line
- the next line must specify the type of charge as, NBO, CHELPG, or AIM
- the following lines must include the atomic number and the corresponding charge
- leave a blank line after the first set of charges
- do steps 4-6 for the next type of charges
- sample input file
- for a positively charged molecule (emim+): Media:emim.log
- for water molecule is below
3 0.0 2.1993 0.0000 0.0000 -2.1993 8 -0.000000 0.000000 0.114720 1 0.000000 0.754031 -0.458881 1 -0.000000 -0.754031 -0.458881 NBO 8 -0.95658 1 0.47829 1 0.47829 CHELPG 8 -0.798274 1 0.399137 1 0.399137 AIM 8 -1.18105730290000 1 0.59052792621000 1 0.59052707781000
running the python script
- if the names of your script and input file are dipole_moment.py and water.inp, respectively, and if you want to save the output as a .txt file, run the script by typing, python dipole_moment.py water.inp txt
#!/usr/bin/env python3
#
#dipole_moment.py
import sys, re, os
dir=os.getcwd()
if len(sys.argv) < 3:
print('to run the script please type: python dipole_moment.py inp_name (with extension) output_file_extension')
sys.exit()
else:
inp=str(sys.argv[1])
split_inp=inp.split('.')
base=split_inp[0]
out_ext=str(sys.argv[2])
out=base+'.'+out_ext
s='directory is '
print('{0:}{1:}'.format(s,dir))
c=open(out,"w")
c.write('{0:}{1:} \n'.format(s,dir))
s='input file is '
print('{0:}{1:}'.format(s,inp))
c.write('{0:}{1:} \n'.format(s,inp))
s='output file is '
print('{0:}{1:} \n'.format(s,out))
c.write('{0:}{1:}\n \n'.format(s,out))
f=open(inp,'r')
line=f.readline()
a1,b1,c1,d1,e1,f1=line.rstrip().split()
Natoms=int(a1)
tQ=float(b1)
count=0
dipole=[]
mux=[]
muy=[]
muz=[]
dipole.append(float(c1))
mux.append(float(d1))
muy.append(float(e1))
muz.append(float(f1))
count=count+1
s='total charge is '
print('{0:}{1:}'.format(s,tQ))
c.write('{0:}{1:} \n'.format(s,tQ))
s='dipole moment, u from gaussian is '
print('{0:}{1:}'.format(s,dipole[0]))
c.write('{0:}{1:} \n'.format(s,dipole[0]))
s='ux from gaussian is '
print('{0:}{1:}'.format(s,mux[0]))
c.write('{0:}{1:} \n'.format(s,mux[0]))
s='uy from gaussian is '
print('{0:}{1:}'.format(s,muy[0]))
c.write('{0:}{1:} \n'.format(s,muy[0]))
s='uz from gaussian is '
print('{0:}{1:}\n'.format(s,muz[0]))
c.write('{0:}{1:}\n \n'.format(s,muz[0]))
atom_num=[]
x=[]
y=[]
z=[]
line=f.readline()
for i in range(0,Natoms):
line=f.readline()
a1,b1,c1,d1=line.rstrip().split()
atom_num.append(int(a1))
x.append(float(b1))
y.append(float(c1))
z.append(float(d1))
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}
# Helper function to detect charge types in the input file
def detect_charge_types(inp):
charge_types = {"CHELPG": False, "NBO": False, "AIM": False}
with open(inp, 'r') as f:
for line in f:
if "CHELPG" in line:
charge_types["CHELPG"] = True
elif "NBO" in line:
charge_types["NBO"] = True
elif "AIM" in line:
charge_types["AIM"] = True
return charge_types
def nboread(inp):
global charge
global Natoms
charge=[]
go=False
f=open(inp,'r')
line=f.readline()
while line:
go=True
line=f.readline()
if "NBO" in line:
s='NBO charges found'
print('{0:}'.format(s))
c.write('{0:} \n'.format(s))
for i in range(0, Natoms):
line=f.readline()
a1,b1=line.rstrip().split()
charge.append(float(b1))
if line=='':go=False
return charge
def chelpgread(inp):
global charge
charge=[]
go=False
f=open(inp,'r')
line=f.readline()
while line:
go=True
line=f.readline()
if "CHELPG" in line:
s='CHELPG charges found'
print('{0:}'.format(s))
c.write('{0:} \n'.format(s))
for i in range(0, Natoms):
line=f.readline()
a1,b1=line.rstrip().split()
charge.append(float(b1))
if line=='':go=False
return charge
def aimread(inp):
global charge
charge=[]
go=False
f=open(inp,'r')
line=f.readline()
while line:
go=True
line=f.readline()
if "AIM" in line :
s='AIM charges found'
print('{0:}'.format(s))
c.write('{0:} \n'.format(s))
for i in range(0, Natoms):
line=f.readline()
a1,b1=line.rstrip().split()
charge.append(float(b1))
if line=='':go=False
return charge
def main(charge):
global DMp
global mup
CoM=[0,0,0]
mass=0
CoP=[0,0,0]
pcharge=0
n=0
mu=[0,0,0]
mup=[0,0,0]
mum=[0,0,0]
CoC=[0,0,0]
COCC=0
uCOCC=0
DM=0
DMp=0
DMm=0
total_charge=0
charge_arm=0
length=0
ca_vector=[0,0,0]
#calculates for the center of pcharge
while n < Natoms:
CoM[0]=CoM[0]+x[n]*masses[atom_num[n]]
CoM[1]=CoM[1]+y[n]*masses[atom_num[n]]
CoM[2]=CoM[2]+z[n]*masses[atom_num[n]]
mass=mass+masses[atom_num[n]]
CoP[0]=CoP[0]+x[n]*atom_num[n]
CoP[1]=CoP[1]+y[n]*atom_num[n]
CoP[2]=CoP[2]+z[n]*atom_num[n]
pcharge=pcharge+atom_num[n]
n=n+1
CoM[0] = CoM[0] / mass
CoM[1] = CoM[1] / mass
CoM[2] = CoM[2] / mass
s='center of mass: '
print('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f}'.format(s,CoM))
c.write('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f} \n'.format(s,CoM))
CoP[0] = CoP[0] / pcharge
CoP[1] = CoP[1] / pcharge
CoP[2] = CoP[2] / pcharge
s='center of pcharge: '
print('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f}'.format(s,CoP))
c.write('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f} \n'.format(s,CoP))
n=0
while n < Natoms:
mup[0]=mup[0]+(charge[n]*(x[n]-CoP[0]))
mup[1]=mup[1]+(charge[n]*(y[n]-CoP[1]))
mup[2]=mup[2]+(charge[n]*(z[n]-CoP[2]))
mum[0]=mum[0]+(charge[n]*(x[n]-CoM[0]))
mum[1]=mum[1]+(charge[n]*(y[n]-CoM[1]))
mum[2]=mum[2]+(charge[n]*(z[n]-CoM[2]))
total_charge=total_charge+charge[n]
n=n+1
DMp = DMp + (((mup[0])**2+(mup[1])**2+(mup[2])**2)**0.5)
DMm = DMm + (((mum[0])**2+(mum[1])**2+(mum[2])**2)**0.5)
s='total charge: '
c.write('{0:}{1:5.2f}\n'.format(s,total_charge))
print('{0:} {1:5.2f}'.format(s,total_charge))
length=((CoP[0]-CoM[0])**2+(CoP[1]-CoM[1])**2+(CoP[2]-CoM[2])**2)**0.5
charge_arm=charge_arm + (total_charge*length)
charge_arm=charge_arm/0.208194
ca_vector[0]=charge_arm*(CoP[0]-CoM[0])/length
ca_vector[1]=charge_arm*(CoP[1]-CoM[1])/length
ca_vector[2]=charge_arm*(CoP[2]-CoM[2])/length
DMp = DMp/0.208194
mup[0]=mup[0]/0.208194
mup[1]=mup[1]/0.208194
mup[2]=mup[2]/0.208194
s='dipole vector wrt CoP (D): '
print('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f}'.format(s,mup))
c.write('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f} \n'.format(s,mup))
s='dipole moment wrt CoP (D): '
print('{0:}{1:>5.2f}'.format(s,DMp))
c.write('{0:}{1:>5.2f} \n'.format(s,DMp))
DMm = DMm/0.208194
mum[0]=mum[0]/0.208194
mum[1]=mum[1]/0.208194
mum[2]=mum[2]/0.208194
s='dipole vector wrt CoM (D): '
print('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f}'.format(s,mum))
c.write('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f} \n'.format(s,mum))
s='dipole moment wrt CoM (D): '
print('{0:}{1:>5.2f}'.format(s,DMm))
c.write('{0:}{1:>5.2f} \n'.format(s,DMm))
s='charge arm vector (D): '
print('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f}'.format(s,ca_vector))
c.write('{0:}{1[0]:>5.2f}, {1[1]:5.2f}, {1[2]:5.2f} \n'.format(s,ca_vector))
s='charge arm (D): '
print('{0:}{1:>5.2f}\n'.format(s,charge_arm))
c.write('{0:}{1:>5.2f} \n\n'.format(s,charge_arm))
# Detect available charge types
charge_types = detect_charge_types(inp)
# Calculate dipole moments only for the detected charge types
methods = []
methods.append("Gaussian")
if charge_types["CHELPG"]:
chelpg_charges = chelpgread(inp)
methods.append("CHELPG")
main(chelpg_charges)
dipole.append(float(DMp))
mux.append(float(mup[0]))
muy.append(float(mup[1]))
muz.append(float(mup[2]))
count=count+1
if charge_types["NBO"]:
nbo_charges = nboread(inp)
methods.append("NBO")
main(nbo_charges)
dipole.append(float(DMp))
mux.append(float(mup[0]))
muy.append(float(mup[1]))
muz.append(float(mup[2]))
count=count+1
if charge_types["AIM"]:
aim_charges = aimread(inp)
methods.append("AIM")
main(aim_charges)
dipole.append(float(DMp))
mux.append(float(mup[0]))
muy.append(float(mup[1]))
muz.append(float(mup[2]))
count=count+1
m=0
s='\nsummary of dipole moments wrt to CoP:'
print('{0:}'.format(s))
c.write('{0:}'.format(s))
s='\nMethod, Dipole Moment u, ux,uy,uz\n'
print('Method, Dipole Moment u, ux,uy,uz')
c.write('{0:}'.format(s))
while m<count:
print('{0:},{1: >05.2f} ,{2: >05.2f} ,{3: >05.2f}, {4: >05.2f}'.format(methods[m],dipole[m],mux[m],muy[m],muz[m]))
c.write('{0:},{1: >05.2f} ,{2: >05.2f} ,{3: >05.2f}, {4: >05.2f} \n'.format(methods[m],dipole[m],mux[m],muy[m],muz[m]))
m=m+1
quit()
sample output of the script above
directory is /nfs/home/ebardoir/bin/test input file is water.inp output file is water.txt total charge is 0.0 dipole moment, u from gaussian is 2.1993 ux from gaussian is 0.0 uy from gaussian is 0.0 uz from gaussian is -2.1993 CHELPG charges found center of mass: 0.00, 0.00, 0.05 center of pcharge: 0.00, 0.00, -0.00 total charge: 0.00 dipole vector wrt CoP (D): 0.00, 0.00, -2.20 dipole moment wrt CoP (D): 2.20 dipole vector wrt CoM (D): 0.00, 0.00, -2.20 dipole moment wrt CoM (D): 2.20 charge arm vector (D): 0.00, 0.00, -0.00 charge arm (D): 0.00 NBO charges found center of mass: 0.00, 0.00, 0.05 center of pcharge: 0.00, 0.00, -0.00 total charge: 0.00 dipole vector wrt CoP (D): 0.00, 0.00, -2.64 dipole moment wrt CoP (D): 2.64 dipole vector wrt CoM (D): 0.00, 0.00, -2.64 dipole moment wrt CoM (D): 2.64 charge arm vector (D): 0.00, 0.00, -0.00 charge arm (D): 0.00 AIM charges found center of mass: 0.00, 0.00, 0.05 center of pcharge: 0.00, 0.00, -0.00 total charge: -0.00 dipole vector wrt CoP (D): 0.00, 0.00, -3.25 dipole moment wrt CoP (D): 3.25 dipole vector wrt CoM (D): 0.00, 0.00, -3.25 dipole moment wrt CoM (D): 3.25 charge arm vector (D): -0.00, -0.00, 0.00 charge arm (D): -0.00 summary of dipole moments wrt to CoP: Method, Dipole Moment u, ux,uy,uz Gaussian, 2.20 , 0.00 , 0.00, -2.20 CHELPG, 2.20 , 0.00 , 0.00, -2.20 NBO, 2.64 , 0.00 , 0.00, -2.64 AIM, 3.25 , 0.00 , 0.00, -3.25
