Merge pull request #1 from zzh535897/zzh535897-patch-1

Zzh535897 patch 1

README.md

@@ -1,11 +1,34 @@
-### QPC-TDSE-CPC
+# QPC-TDSE-CPC
 The CPC library version of QPC-TDSE[[1]](https://doi.org/10.1016/j.cpc.2023.108787). 
 
 Stable version. 
 
 Please note this repo is for documentation. Usual developement will NOT be made to this repo. Bug fixing and new features are tested elsewhere, and will be gathered to this repo if a new stable version is formulated.
 
-### References
+
+## Issues
+2023.9.11<br>
+(1) The missing scaling factor sqrt(2E) in the analysis script plot_part.m has been corrected. 
+
+(2) The format of output "proj_val" for the real-time projection onto bound states has been corrected. Now it follows the format in the user guide.
+
+## Important
+2023.9.11<br>
+(1) The variable "lmd1" in an output file actually differs a phase factor to the P_{lm}(E) as
+
+    LMD(k,m,l)= (-i)^l \exp(i\Delta_{l}(k)) P_{lm}(E)
+
+Here P_{lm}(E) is defined in eq.(58). This feature will not affect the use of eq.(65), but could cause confusion when calculating physical quantities such as the time-delay.
+
+(2) The gaussian envelope in the original release is really "gaussian", which means it is non-zero everywhere. The simulation starts and ends at 2FWHM far from its central peak. If the PCS method is applied for this case, some unphysical results would occur since the instantaneous A(t_f) is non-zero. 
+
+A modified gaussian envelope with its tail cut should be customized instead. One may refer to include/structure/field.h for its format.
+
+(3) The projection onto field-free eigenstates could differ a sign between two runs using different box size parameters in the paper version code, since the LAPACK diagonalization routines randomly choose the sign of the eigenvectors. This may cause some confusion.
+
+In the latest version we always fix the sign of eigen states to be positive near the origin.
+
+## References
 [1] Zhao-Han Zhang, Yang Li, Yi-Jia Mao, Feng He,
 *QPC-TDSE: A parallel TDSE solver for atoms and small molecules in strong lasers*
 [Comput. Phys. Comm., **290** 108787 (2023)]

analysis/plot_part.m

@@ -3,18 +3,18 @@
     figure;hold on;
     for im=1:nm
         for il=1:nl
-            plot(pr.^2/2,squeeze(abs(lmd(im,il,:))).^2,'LineWidth',1.1);
+            plot(pr.^2/2,pr.*squeeze(abs(lmd(im,il,:))).^2,'LineWidth',1.1);
         end
     end
     xlabel('$E$(a.u.)'        ,'Interpreter','latex');
     ylabel('$P_{lm}(E)$(a.u.)','Interpreter','latex');
     title('Partial-Wave PMD');
 
     figure;hold on;
-    pes=sum(sum(abs(lmd).^2,1),2);
+    pes=pr.*sum(sum(abs(lmd).^2,1),2);
     plot(pr.^2/2,pes,'k-','LineWidth',1.1);
 
     xlabel('$E$(a.u.)'   ,'Interpreter','latex');
     ylabel('$P(E)$(a.u.)','Interpreter','latex');
     title('PES');
-end
+end

binsrc/sphTDSE.cpp

@@ -438,8 +438,11 @@ namespace proj
 			for(size_t im=0;im<dims_t::m_dims;++im)
 			for(size_t il=0;il<dims_t::l_dims;++il)
 			{
-				for(size_t ii=0;ii<nnz[im*dims_t::l_dims+il];++ii)
-				val.push_back(tmp[im*dims_t::l_dims*dims_t::n_dims+il*dims_t::n_dims+ii]);
+				auto beg = tmp.begin()+im*dims_t::l_dims*dims_t::n_dims+il*dims_t::n_dims;
+                                auto end = beg + nnz[im*dims_t::l_dims+il];
+                                val.insert(val.end(),beg,end);
+				//for(size_t ii=0;ii<nnz[im*dims_t::l_dims+il];++ii)
+				//#val.push_back(tmp[im*dims_t::l_dims*dims_t::n_dims+il*dims_t::n_dims+ii]);
 			}
 		}
 	}//end of calc_proj
@@ -449,7 +452,14 @@ namespace proj
 	{
 		if constexpr(calculate_proj==1)
 		{
-			proj->projection(coef,val,nnz);//calculate projection onto eigen states
+			proj->projection(coef,tmp,nnz);//calculate projection onto eigen states
+			for(size_t im=0;im<dims_t::m_dims;++im)
+                        for(size_t il=0;il<dims_t::l_dims;++il)
+                        {
+                                auto beg = tmp.begin()+im*dims_t::l_dims*dims_t::n_dims+il*dims_t::n_dims;
+                                auto end = beg + nnz[im*dims_t::l_dims+il];
+                                val.insert(val.end(),beg,end);
+                        }
 			file.save((double*)val.data(),val.size()*2ul,"/proj_val");
 			file.save((size_t*)nnz.data(),nnz.size()    ,"/proj_nnz");
 		}

include/spherical/spectrum/population_1e.hpp

@@ -92,7 +92,7 @@ struct	eigenstate_sph final//to store all required eigenstates and eigenvalues (
 				for(size_t i=0;i<dims_t::n_dims;++i)
 				{	
 					if(crit(l,temp_val[i],temp_vec.data()+i*dims_t::n_dims))
-					crit_result.push_back(i);//requirement: crit(m,l,eig,vec)
+					crit_result.push_back(i);//requirement: crit(l,eig,vec)
 				}
 				//dump all required states
 				auto& vec_this	=	vec(l);
@@ -103,6 +103,7 @@ struct	eigenstate_sph final//to store all required eigenstates and eigenvalues (
 				{
 					size_t s=crit_result[i];
 					val_this[i]=temp_val[s];
+					double sig = temp_vec[s*dims_t::n_dims]>0.? 1.:-1.;
 					//to store the left hand side basis, multiply by S
 					intrinsic::symb_mul_vecd<dims_t::n_dims,dims_t::n_elem,0>
 					(
@@ -221,41 +222,3 @@ struct	eigenstate_bic final
 	}//end of initialize
 
 };//end of eigenstate_bic
-
-
-/*template<class dims_t>
-struct	[[deprecated]] statistics
-{
-	static inline auto	generate_axis_m	()
-	{//get the m value on (im,il)-grid
-		auto	axis	=	std::vector<long>(dims_t::m_dims*dims_t::l_dims);
-		for(size_t im=0;im<dims_t::m_dims;++im)
-		for(size_t il=0;il<dims_t::l_dims;++il)
-		{
-			axis[im*dims_t::l_dims+il]=dims_t::in_m(im);
-		}
-		return	axis;//use NRVO or move
-	}//end of generate_axis_m
-	static inline auto	generate_axis_l	()
-	{//get the l value on (im,il)-grid
-		auto	axis	=	std::vector<long>(dims_t::m_dims*dims_t::l_dims);
-		for(size_t im=0;im<dims_t::m_dims;++im)
-                for(size_t il=0;il<dims_t::l_dims;++il)
-                {
-			axis[im*dims_t::l_dims+il]=dims_t::in_l(im,il);
-                }
-                return  axis;//use NRVO or move
-	}//end of generate_axis_l
-	static inline auto	generate_dist_ml(const coefficient_view<dims_t>& coef,const operator_sc<dims_t>& oper)
-	{//evaluate the polulation |P(m,l)|^2 on (im,il)-grid
-		auto	dist	=	std::vector<double>(dims_t::m_dims*dims_t::l_dims);
-		#pragma omp parallel for collapse(2)
-		for(size_t im=0;im<dims_t::m_dims;++im)
-                for(size_t il=0;il<dims_t::l_dims;++il)
-		{
-			dist[im*dims_t::l_dims+il]=norm(oper.observe_rsub(coef(im,il),coef(im,il)));
-		}
-		return 	dist;
-	}//end of generate_dist_ml
-};//end of statistics*/
-