I have the following code that gets some input files and calculates some values. The issue is, the code is working for Mus_df values 1,2,3,4,6,12 but not 5,7,8,9,10,11. The error indicates that CCvecs is None-type. Can you help me?
import numpy as np
import matplotlib.pyplot as plt
import math
import itertools
from math import sqrt, pi, sin, cos, asin, acos, exp
import pandas as pd
cmap1 = plt.get_cmap('tab10')
ang2au = 1.8897259886
eV2au = 3.6749322175655e-2
debye2au = 0.393456
def readParams(NpAndMol_xyz_path, Mol_xyz_path, Np_Mu_path, Nmol, t=0.0):
# read NP+pyr and pyr XYZ into pd df
# read TDMs of Ag174 NP system, select electronic states? based on cutoff for f..???
NpAndMol_xyz_tmp = pd.read_csv(
NpAndMol_xyz_path, sep='\s+', header=None, skiprows=2)
Mol_xyz_df = pd.read_csv(Mol_xyz_path, sep='\s+', header=None, skiprows=2)
Mus_tmp_df = pd.read_csv(Np_Mu_path, sep='\s+', header=None,
skiprows=2, names=['no', 'E', 'f', 'x', 'y', 'z'])
# parse NpAndMol_xyz
NpAndMol_xyz_df = NpAndMol_xyz_tmp.drop(
NpAndMol_xyz_tmp[NpAndMol_xyz_tmp[0] == 'Ag'].index)
NP_xyz = NpAndMol_xyz_tmp[NpAndMol_xyz_tmp[0] == 'Ag'].to_numpy()
MolsOnNP_xyz = np.array_split(NpAndMol_xyz_df, Nmol)
# get only xyz of TDMs
# threshold t for oscillator strength f optional
Mus_tmp_dfRed = Mus_tmp_df.loc[(Mus_tmp_df['f'] >= t)]
Mus_df = Mus_tmp_dfRed.drop(Mus_tmp_dfRed.columns[[1, 2]], axis=1)
elStateEs = Mus_tmp_dfRed['E'].tolist()
return NP_xyz, MolsOnNP_xyz, Mol_xyz_df, Mus_df, elStateEs
def getVecs(MolsOnNP_xyz, Mol_xyz_df):
Nmol = len(MolsOnNP_xyz)
CMOnNPvecs = []
NNvecs = []
RotMats = []
CCvecs = []
# get CoM of pyr molecules, CM(NP)=(0,0,0)
Mol_xyzMat = Mol_xyz_df.drop(Mol_xyz_df.columns[[0]], axis=1).to_numpy()
NNOrig = Mol_xyzMat[1, :] - Mol_xyzMat[3, :]
#testvecOrig = Mol_xyzMat[2, :] - Mol_xyzMat[5, :]
for i in range(Nmol):
MolOnNP_xyz = MolsOnNP_xyz[i]
MolOnNP_xyzMat = MolsOnNP_xyz[i].drop(
MolOnNP_xyz.columns[[0]], axis=1).to_numpy()
# compute NN vector
NsXYZ = MolOnNP_xyz[MolOnNP_xyz[0] == 'N']
NsXYZ_mat = NsXYZ.drop(NsXYZ.columns[[0]], axis=1).to_numpy()
NNOnNPvec = NsXYZ_mat[0, :] - NsXYZ_mat[1, :]
NNOnNPvec_normal = NNOnNPvec * 1/np.linalg.norm(NNOnNPvec)
CMOnNPvecTmp = (MolOnNP_xyzMat[2, :] + MolOnNP_xyzMat[5, :])*0.5
# compute CC vec parallel to TDM vector
# CCvec = MolOnNP_xyzMat[3, :] - MolOnNP_xyzMat[1, :]
CCvec = findCCVec(MolOnNP_xyz,NNOnNPvec)
# print(np.linalg.norm(CCvec))
CCvecs.append(CCvec)
#print(len(CCvecs))
NNvecs.append(NNOnNPvec)
CMOnNPvecs.append(CMOnNPvecTmp)
return CCvecs, CMOnNPvecs, NNvecs
def alignAntiParallelVec(u, v):
if np.dot(u, v) < 0:
print('antiparallel')
return -v
elif np.dot(u, v) > 0:
print('parallel')
return v
elif np.dot(u, v) == 0:
print('vectors are orthogonal')
return v
def getCoupls(CCvecs, CMOnNPvecs, Mus_df, mydist=None, printOut=False, unit='au', multiplier=1.0):
#muMolBf = np.array([0.00000, 0.78457, 0.0000])
muMolBf = np.array([0.00000, -0.8915, 0.0000])
muMolBfAbs = np.linalg.norm(muMolBf)
musNP = Mus_df[['x', 'y', 'z']].to_numpy()
print('take mu vectors from input:')
print('\n###Cluster mu (modulus) in au###')
print([np.linalg.norm(v) for v in musNP])
print('\n###Molecule mu (modulus) in au###')
print(muMolBfAbs)
elStateInds = Mus_df['no'].to_list()
Nmol = len(CCvecs)
print(Nmol)
Nel = len(Mus_df.index)
gcoupls = np.zeros((Nmol, Nel))
uvec_out = []
for i in range(Nmol):
uvec = 1/np.linalg.norm(CMOnNPvecs[i]*ang2au) * CMOnNPvecs[i]*ang2au
print(CCvecs[i])
muMolVec = CCvecs[i]*ang2au
#muMolVec = alignAntiParallelVec(uvec, CCvecs[i])
muMolSf = muMolBfAbs/np.linalg.norm(muMolVec) * muMolVec
# uvec_out.append(muMolVec)
if mydist == None:
dist = np.linalg.norm(CMOnNPvecs[i])
dist_au = dist*ang2au
print('###Calculate dip-dip distance from coordinate file###')
print('dist = {} å = {} au'.format(dist, dist_au))
else:
print('###Use custom dist###')
dist = mydist
dist_au = dist*ang2au
for j in range(len(musNP)):
gcoupls[i, j] = 1/dist_au**3 * (np.dot(musNP[j], muMolSf) -
3*np.dot(musNP[j], uvec)*np.dot(muMolSf, uvec))
# g_ij matrix, i=mol index, j=el state index
if printOut:
print('\n\n###coupling matrix g_ij; i=mol, j=el state###')
for i in range(Nmol):
for j in range(Nel):
if unit == 'ev':
print('g_{}_{} = {:.5f}, ev'.format(
i+1, j+1, gcoupls[i, j]/eV2au*multiplier))
elif unit == 'au':
print('g_{}_{} = {:.5f}'.format(
i+1, j+1, gcoupls[i, j]*multiplier))
else:
print('unknown unit!')
print('\n###metal mus###')
for i in range(Nel):
# read and print maximum TDM component
muInd = np.argmax(np.abs(musNP[i]))
print("mu_1_{} = {:.4f}".format(i+1, musNP[i][muInd]))
return gcoupls
def findCCVec(molXYZ_df,NNvec):
CsXYZ = molXYZ_df[molXYZ_df[0] == 'C']
CsXYZ_Mat = CsXYZ.drop(CsXYZ.columns[[0]], axis=1).to_numpy()
Nmols = len(CsXYZ_Mat[:])
Inds = range(Nmols)
iters = itertools.combinations(Inds,2)
failCount = 0
for iter in iters:
CCtestVec = CsXYZ_Mat[iter[0]] - CsXYZ_Mat[iter[1]]
projCCNN = CCtestVec.dot(NNvec)
if np.abs(projCCNN) <= 1e-7:
return CCtestVec
else:
failCount += 1
if failCount == Nmols:
print("No suitable CC axis found")
return None
def printEs(elist, shift=0.0):
print("omega_0 = 0.0000, ev")
for i, e in enumerate(elist):
print("omega_{} = {:.4f}, ev".format(i+2, e+shift))
def printCouplOp(MolDofs, MolTags, NElStates, elDof):
for molTag, molDof in zip(MolTags, MolDofs):
for j in range(NElStates):
print("g_{}_{} |{} F |{} S1&{}".format(
molTag, j+1, molDof, elDof, j+2))
for j in range(int(NElStates/2)):
print("omega_1 |{} S{}&{}".format(elDof, 2*j+2, 2*j+2))
print("omega_2 |{} S{}&{}".format(elDof, 2*j+3, 2*j+3))
for j in range(int(NElStates/2)):
print("omega_1 |{} S{}&{}".format(elDof, 2*j+2, 2*j+2))
print("omega_2 |{} S{}&{}".format(elDof, 2*j+3, 2*j+3))
def visVecs(cms, vecs, l=0.1, angle=0, pts = None, plotPts = False):
ax = plt.figure().add_subplot(projection='3d')
ax.view_init(30, angle)
for i in range(len(cms)):
ax.scatter(cms[i,0].tolist(),cms[i,1].tolist(),cms[i,2].tolist(),color=cmap1(i),marker='*',s=50)
ax.quiver(cms[:,0],cms[:,1], cms[:,2], vecs[:,0], vecs[:,1], vecs[:,2], length=l, normalize=True, color='black')
if plotPts:
#print(pts[:,1])
#print(pts[:,2])
#print(pts[:,3])
#ax.scatter(pts[:,1],pts[:,2],pts[:,3])
ax.scatter(pts[:,1].tolist(),pts[:,2].tolist(),pts[:,3].tolist())
plt.show()
#visVecs(np.array(CMOnNPvecs), np.array(CCvecs), l=1, angle=40)
Np_xyz , MolsOnNP_xyz, Mol_xyz_df, Mus_df, elStateEs = readParams(
'Ag147-py12.xyz', 'Pyr_xyz.dat', 'Ag147_mu_sym_e1u_xy.dat',8 , t=0.000)
CCvecs, CMOnNPvecs, NNvecs = getVecs(MolsOnNP_xyz, Mol_xyz_df)
coupls = getCoupls(CCvecs, CMOnNPvecs, Mus_df, mydist = None , printOut = True, unit = 'ev', multiplier=1.0)
The Mus_df value is prompted in this line Np_xyz , MolsOnNP_xyz, Mol_xyz_df, Mus_df, elStateEs = readParams( 'Ag147-py12.xyz', 'Pyr_xyz.dat', 'Ag147_mu_sym_e1u_xy.dat',8 , t=0.000), which is set to be 8 in this case.
I tried debugging getVecs, getCoupls, findCCvec functions but had no luck. Can you please help me?
