Quantum solver

.gitignore

@@ -0,0 +1 @@
+*.dat

README.md

@@ -3,4 +3,6 @@
 This is used to test the performance of TAB on Tully's model.
 
 To run TAB dynamics, simply run 
-'python3 name-main.py'
+'python3 namd-main.py'
+
+Initial conditions, calcH, and diffH must be set manually.

analysis.py

@@ -0,0 +1,48 @@
+import pandas as pd
+import matplotlib.pyplot as plt
+import glob
+import os
+# Specify the paths to the directories containing the ene*.dat files
+data_directory1 = '/home/adurden/Ari/TAB/tab-model/data6'
+data_directory2 = '/home/adurden/Ari/TAB/tab-model/data4/old'
+# Define the column names
+
+def process_data2(directory):
+    columns2 = ["t", "S0 pop", "S1 pop", "S2 pop", "poptot"]
+    file_paths = glob.glob(os.path.join(directory, 'dpop*.dat'))
+    dfs = []        
+    count=0
+    for file_path in file_paths:
+        df = pd.read_csv(file_path, delim_whitespace=True, names=columns2, skiprows=1)
+        dfs.append(df)
+        count+=1
+    print(count)
+    combined_df = pd.concat(dfs)
+    average_df = combined_df.groupby('t').mean().reset_index()
+    return average_df, len(file_paths)
+avg_df1, file_count1 = process_data2(data_directory1)
+avg_df2, file_count2 = process_data2(data_directory2)
+#avg_df3, file_count3 = process_data2(data_directory3)
+# Plot the data
+
+columns_pop = ["t", "S0 pop", "S1 pop", "S2 pop", "poptot"]
+pop_df = pd.read_csv('/home/adurden/Ari/TAB/tab-model/data4/pop.dat', delim_whitespace=True, names=columns_pop, index_col=False)
+pop_df['t'] = pop_df['t'] / 100
+
+
+plt.figure(figsize=(12, 6))
+plt.plot(avg_df1['t'], avg_df1['S0 pop'], label='new')
+plt.plot(avg_df2['t'], avg_df2['S0 pop'], label='old')
+#plt.plot(avg_df3['t'], avg_df3['S0 pop'], label='rescales all directions')
+plt.scatter(pop_df['t'], pop_df['S0 pop'], label='exact', color='black', marker='o', s=7)  # Plot as dots with smaller size
+plt.xlabel('Time (fs)', fontsize=24)
+plt.ylabel('S0 Population', fontsize=24)
+
+#plot x axis up to 300 fs
+plt.xlim(0, 300)
+
+plt.title('S0 pop vs Time', fontsize=24)
+plt.legend(fontsize=18)
+plt.grid(True)
+plt.savefig('momentum_rescale.pdf', format='pdf')
+#plt.show()

calcH.py

@@ -1,18 +1,17 @@
-def buildH(dimH, x1, x2, w1, w2, c, delta):
+def buildH(dimH, x, w1, w2, c, delta):
 #Build the Hamiltonian for at position (x1,x2)
 
 	import numpy as np
 
 	#initialize Hamiltonian matrix
 	H = np.zeros((dimH,dimH))
-
-	H[0][0] = -1.0*w1*x1   #diabatic potential
+	H[0][0] = -1.0*w1*x[0]   #diabatic potential
 
 	i = 1
 
 	while i < dimH:
-		H[i][i] = w2*x1 - (i-1)*delta  #diabatic potential
-		H[i][0] = c*x2		#diabatic coupling
+		H[i][i] = w2*x[0] - (i-1)*delta  #diabatic potential
+		H[i][0] = c*x[1]		#diabatic coupling
 		H[0][i] = H[i][0]
 		i = i + 1
 	pass

cgauss.py

@@ -1,4 +1,4 @@
-def gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol,npthresh,pehrptol,odotrho,tolodotrho,nta,dtw,zpop,dgscale): # Determines which coherent sub-block 
+def gcollapse(dimH,ndof,deltatn,aforce,poparray,dcp,nzthresh,errortol,npthresh,pehrptol,odotrho,tolodotrho,nta,dtw,zpop,dgscale): # Determines which coherent sub-block 
 	"""of the electronic density matrix the WF collapses into"""
 
 	# Standard library imports =====================================
@@ -26,7 +26,10 @@ def gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol,
 	while i < dimH:
 		j = i + 1
 		while j < dimH:
-			invtau[i][j] = ((aforce1[i]-aforce1[j])**(2.0)/(8.0*dcp1)+(aforce2[i]-aforce2[j])**(2.0)/(8.0*dcp2))**(0.50)
+			invtausum = 0.0
+			for k in range(ndof):
+				invtausum+= (aforce[k, i]-aforce[k, j])**2/(8.0*dcp[k])
+			invtau[i][j] = (invtausum)**(0.50)
 			invtau[j][i] = invtau[i][j]
 			j = j + 1
 		i = i + 1
@@ -370,7 +373,7 @@ def gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol,
 	# constructing npop from the density matrix
 
 	if (track == 0):
-		return poparray
+		return poparray, track
 	pass
 
 	k = 0
@@ -390,5 +393,5 @@ def gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol,
 #	print A[track]
 
 
-	return npop
+	return npop, track