Removed some dependencies on libraries/utilities.h

libraries/lsu3shell/lsu3shell_rme.cpp

@@ -454,7 +454,7 @@ namespace lsu3shell
               double residual = std::fabs(rme2-rme1);
               ++entries_compared;
               max_residual = std::max(max_residual,residual);
-              total_sqr_residual += sqr(residual);
+              total_sqr_residual += mcutils::sqr(residual);
               bool entries_agree = (residual <= tolerance);
               // std::cout<<std::endl<<fmt::format("residual {}  tolerance {}  bool {}", residual,tolerance,entries_agree)<<std::endl<<std::endl;
               success &= entries_agree;

libraries/moshinsky/moshinsky_xform.cpp

@@ -8,10 +8,10 @@
 ****************************************************************/
 #include <cmath>
 #include "fmt/format.h"
-
+#include "am/halfint.h"
 #include "sp3rlib/u3coef.h"
 #include "moshinsky/moshinsky_xform.h"
-
+#include "mcutils/gsl.h"
 extern double zero_threshold;
 
 namespace u3shell
@@ -22,11 +22,11 @@ namespace u3shell
     int Kmax=std::max(std::max(int(J+M),int(J-M)),int(J-Mp));
 
     for(int K=0; K<=Kmax; K++)
-      moshinsky_coef+=parity(K)*Choose(int(J+M),int(J-Mp-K))*Choose(int(J-M),K);
+      moshinsky_coef+=ParitySign(K)*mcutils::Choose(int(J+M),int(J-Mp-K))*mcutils::Choose(int(J-M),K);
     moshinsky_coef
-      *=parity(int(J-Mp))
-      *sqrt(Factorial(int(J+Mp))*Factorial(int(J-Mp))
-           /(pow(2.,TwiceValue(J))*Factorial(int(J+M))*Factorial(int(J-M)))
+      *=ParitySign(int(J-Mp))
+      *sqrt(mcutils::Factorial(int(J+Mp))*mcutils::Factorial(int(J-Mp))
+           /(pow(2.,TwiceValue(J))*mcutils::Factorial(int(J+M))*mcutils::Factorial(int(J-M)))
           );
     return moshinsky_coef;
   }
@@ -116,7 +116,7 @@ namespace u3shell
           // int exchange_phase=phase_even?1:-1;
           // overall factor for the bra
           // the factor of 1/sqrt(1+delta) comes from the normalization for particles in the same shell
-          double coef=norm?(1./std::sqrt(1.+KroneckerDelta(eta1p,eta2p))):1;
+          double coef=norm?(1./std::sqrt(1.+mcutils::KroneckerDelta(eta1p,eta2p))):1;
           bra_moshinky_12(bra_state_index,0)=coef*MoshinskyCoefficient(etap, eta_cm, eta1p, eta2p,xp);
         } 
       // iterate over ket subspace
@@ -129,7 +129,7 @@ namespace u3shell
 
           // overall factor for the ket
           // the factor of 1/sqrt(1+delta) comes from the normalization for particles in the same shell
-          double coef=norm?(1./std::sqrt(1.+KroneckerDelta(eta1,eta2))):1;
+          double coef=norm?(1./std::sqrt(1.+mcutils::KroneckerDelta(eta1,eta2))):1;
           ket_moshinky_12(0,ket_state_index)=coef*MoshinskyCoefficient(eta, eta_cm,eta1,eta2,x);
         }  
       ////////////////////////////////////////////////////////////////////////////////////////////////////////////////

libraries/moshinsky/moshinsky_xform.h

@@ -19,6 +19,7 @@
                     TransformRelativeUnitTensorToTwobodyTensor.
                  Renamed MoshinskyTransformUnitTensor to 
                     MoshinskyTransformTensor.
+  3/19/21 (aem): Removed dependence on utilities/utilities.h in favor of mcutils
 ****************************************************************/
 
 #ifndef MOSHINSKY_H_

libraries/moshinsky/relative_cm_xform.cpp

@@ -85,7 +85,7 @@ namespace u3shell
                     u3shell::RelativeCMBraketNLST braket_rel_cm(L0,S0,T0,bra_rel_cm,ket_rel_cm);
                     rel_cm_lst_map[braket_rel_cm]
                       =am::Unitary9J(L0,0,L0,Lr,Lcm,L,Lrp,Lcm,Lp)
-                      *parity(Lr+Lrp+L+Lp)*rme_rel;
+                      *ParitySign(Lr+Lrp+L+Lp)*rme_rel;
                   }
       }
   }
@@ -124,7 +124,7 @@ namespace u3shell
                   {
                     u3shell::RelativeCMUnitTensorLabelsU3ST tensor_cm(x0,S0,T0,rho0,bra,ket);
                     relative_cm_u3st_map[tensor_cm]
-                      =//parity(x.lambda()+x.mu()+xp.lambda()+xp.mu()+Np+N)
+                      =//ParitySign(x.lambda()+x.mu()+xp.lambda()+xp.mu()+Np+N)
                         u3::U(x0,xr,xp,x_cm,xrp,1,1,x,1,rho0);
                   }
                 }

libraries/sp3rlib/u3.h

@@ -28,7 +28,8 @@
 #include "am/halfint.h"
 #include "am/am.h"
 #include "sp3rlib/multiplicity_tagged.h"
-#include "utilities/utilities.h"
+#include "mcutils/arithmetic.h"
+// #include "utilities/utilities.h"
 
 namespace u3 
 {
@@ -156,7 +157,7 @@ namespace u3
   inline double Casimir2( const u3::SU3& x)
   //Second order Casimir 
   {
-    return 2./3*(sqr(x.lambda())+x.lambda()*x.mu()+sqr(x.mu())+3*x.lambda()+3*x.mu());
+    return 2./3*(mcutils::sqr(x.lambda())+x.lambda()*x.mu()+mcutils::sqr(x.mu())+3*x.lambda()+3*x.mu());
   } 
 
   inline double Casimir3(const u3::SU3& x)

libraries/sp3rlib/vcs.h

@@ -29,7 +29,8 @@
 
 #include <eigen3/Eigen/Eigen>
 #include <unordered_map>
-
+#include "basis/operator.h"
+#include "mcutils/arithmetic.h"
 #include "sp3rlib/sp3r.h"  
 
 namespace vcs
@@ -56,9 +57,9 @@ namespace vcs
   //   (double) : Omega factor
   {
     double value=0;
-    value += (2*sqr(double(omega.f1()))-sqr(double(n.f1()))+8*(double(omega.f1())-double(n.f1()))-2*(2*double(omega.f1())-double(n.f1())));
-    value += (2*sqr(double(omega.f2()))-sqr(double(n.f2()))+8*(double(omega.f2())-double(n.f2()))-4*(2*double(omega.f2())-double(n.f2())));
-    value += (2*sqr(double(omega.f3()))-sqr(double(n.f3()))+8*(double(omega.f3())-double(n.f3()))-6*(2*double(omega.f3())-double(n.f3())));
+    value += (2*mcutils::sqr(double(omega.f1()))-mcutils::sqr(double(n.f1()))+8*(double(omega.f1())-double(n.f1()))-2*(2*double(omega.f1())-double(n.f1())));
+    value += (2*mcutils::sqr(double(omega.f2()))-mcutils::sqr(double(n.f2()))+8*(double(omega.f2())-double(n.f2()))-4*(2*double(omega.f2())-double(n.f2())));
+    value += (2*mcutils::sqr(double(omega.f3()))-mcutils::sqr(double(n.f3()))+8*(double(omega.f3())-double(n.f3()))-6*(2*double(omega.f3())-double(n.f3())));
     return value/4.;
 
   }

libraries/spncci/recurrence.cpp

@@ -18,6 +18,7 @@
 #include "sp3rlib/u3coef.h"
 #include "spncci/spncci_common.h"
 #include "spncci/transform_basis.h"
+#include "utilities/utilities.h"
 
 extern double zero_threshold;
 

libraries/u3shell/interaction_truncation.cpp

@@ -241,7 +241,7 @@ namespace u3shell
       // Read in the interaction from file
       basis::RelativeSpaceLSJT relative_space_lsjt;
       std::array<basis::RelativeSectorsLSJT,3> isospin_component_sectors_lsjt;
-      std::array<basis::MatrixVector,3> isospin_component_blocks_lsjt;
+      std::array<basis::OperatorBlocks<double>,3> isospin_component_blocks_lsjt;
 
       basis::RelativeOperatorParametersLSJT operator_labels;
 

libraries/u3shell/relative_branching.cpp

@@ -49,7 +49,7 @@ namespace u3shell
               u3shell::RelativeStateLabelsNLST ket_nlst(N,L,S,T);
               u3shell::RelativeStateSectorNLST sector(L0,S0,T0,bra_nlst,ket_nlst);
               // double coef=u3::WCached(w_cache,u3::SU3(N,0), 1,L, x0, kappa0, L0, u3::SU3(Np,0),1,Lp,1);
-              rmes_nlst[sector]+=parity(n+np)*u3::WCached(w_cache,u3::SU3(N,0), 1,L, x0, kappa0, L0, u3::SU3(Np,0),1,Lp,1)*rme;
+              rmes_nlst[sector]+=ParitySign(n+np)*u3::WCached(w_cache,u3::SU3(N,0), 1,L, x0, kappa0, L0, u3::SU3(Np,0),1,Lp,1)*rme;
               // std::cout<<fmt::format("{} {} {} {}  {} {} {} {}  {} {} {} {} {}  {}  {} ",Np,Lp,Sp,Tp,N,L,S,T, x0, kappa0,L0,S0,T0, )
             }
       }
@@ -62,7 +62,7 @@ namespace u3shell
     const basis::RelativeSpaceLSJT& relative_space_lsjt,
     const std::unordered_map<u3shell::RelativeStateSectorNLST,double,boost::hash<u3shell::RelativeStateSectorNLST>>& rmes_nlst,
     std::array<basis::RelativeSectorsLSJT,3>& isospin_component_sectors_lsjt,
-    std::array<basis::MatrixVector,3>& isospin_component_blocks_lsjt
+    std::array<basis::OperatorBlocks<double>,3>& isospin_component_blocks_lsjt
     )
 
   {  // setting up J sectors 
@@ -139,7 +139,7 @@ namespace u3shell
     const u3shell::RelativeRMEsU3ST& interaction_u3st,
     const basis::RelativeSpaceLSJT& relative_space_lsjt,
     std::array<basis::RelativeSectorsLSJT,3>& isospin_component_sectors_lsjt,
-    std::array<basis::MatrixVector,3>& isospin_component_blocks_lsjt
+    std::array<basis::OperatorBlocks<double>,3>& isospin_component_blocks_lsjt
     )
   {  
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

libraries/u3shell/relative_branching.h

@@ -33,15 +33,15 @@ namespace u3shell
     const basis::RelativeSpaceLSJT& relative_space_lsjt,
     const std::unordered_map<u3shell::RelativeStateSectorNLST,double,boost::hash<u3shell::RelativeStateSectorNLST>>& rmes_nlst,
     std::array<basis::RelativeSectorsLSJT,3>& isospin_component_sectors_lsjt,
-    std::array<basis::MatrixVector,3>& isospin_component_blocks_lsjt
+    std::array<basis::OperatorBlocks<double>,3>& isospin_component_blocks_lsjt
     );
    // Branching relative matrix elements from NLST to NLSJT.  J-branched rmes stored in block container
 
 	void BranchRelativeRMEs(const basis::OperatorLabelsJT& operator_labels,int Nmax, int Jmax, 
     const u3shell::RelativeRMEsU3ST& interaction_u3st,
     const basis::RelativeSpaceLSJT& relative_space_lsjt,
     std::array<basis::RelativeSectorsLSJT,3>& isospin_component_sectors_lsjt,
-    std::array<basis::MatrixVector,3>& isospin_component_blocks_lsjt
+    std::array<basis::OperatorBlocks<double>,3>& isospin_component_blocks_lsjt
     );
 	// Branch relative rmes from U3ST to NLSJT
 

libraries/u3shell/two_body_operator.cpp

@@ -80,8 +80,8 @@ namespace u3shell {
                  );
         int biquad_grade = eta1+eta2+ConjugationGrade(x)+int(S)+int(T);
         double biquad_norm_factor = 1/(
-                             sqrt(1+KroneckerDelta(eta1p,eta2p))
-                             *sqrt(1+KroneckerDelta(eta1,eta2))
+                             sqrt(1+mcutils::KroneckerDelta(eta1p,eta2p))
+                             *sqrt(1+mcutils::KroneckerDelta(eta1,eta2))
                              );
 
         // iterate over outer multiplicity on cross projector

libraries/u3shell/upcoupling.cpp

@@ -192,7 +192,7 @@ namespace u3shell
                         double u3_coef=u3::WCached(w_cache,ket.x(), 1,L, x0, kappa0, L0, bra.x(),1,Lp,1)
                           *u3::dim(x0)*am::dim(Lp)/(u3::dim(bra.x())*am::dim(L0));
                         std::tuple<u3shell::RelativeUnitTensorLabelsU3ST,int,int> key(operator_labels,kappa0,L0);
-                        rme_map[key]+=u3_coef*rme_nlst*parity(n+np);
+                        rme_map[key]+=u3_coef*rme_nlst*ParitySign(n+np);
                         /*if (fabs(rme_map[key]) <= zero_threshold) {
                           rme_map[key] = 0;
                         }*/

programs/lgi/get_spncci_seed_blocks.cpp

@@ -35,8 +35,6 @@ files in directory seeds containing:
 
   SPDX-License-Identifier: MIT
 
-
-
 12/29/17 (aem): Created. 
 ****************************************************************/
 #include <fstream>  
@@ -149,7 +147,7 @@ int main(int argc, char **argv)
   ////////////////////////////////////////////////////////////////
   // std::cout<<"Generate tensors "<<std::endl;
   // Get list of unit tensor labels between lgi's 
-  // TODOL change to restrict N0 after testing
+  // TODO: Can we restrict N0?
   bool restrict_positive_N0=false;
   std::vector<u3shell::RelativeUnitTensorLabelsU3ST> lgi_unit_tensor_labels;
   u3shell::GenerateRelativeUnitTensorLabelsU3ST(

programs/operators/generate_relative_u3st_operators.cpp

@@ -10,6 +10,7 @@
   2/21/17 (aem,mac): Update input parsing.  Add parsing checks.
   5/14/19 (aem): Updated   basis::OperatorLabelsJT to 
     basis::RelativeOperatorParametersLSJT for reading in operators
+  3/19/21 (aem) : Removed dependency on utilities/utilities.h
 ****************************************************************/
 #include <iostream>
 #include <fstream>
@@ -107,7 +108,7 @@ namespace u3shell
     int kappa0=1; 
     int L0=0;
 
-    double intrinsic_factor=2.*(1.+KroneckerDelta(moshinsky_convention,false))/A;
+    double intrinsic_factor=2.*(1.+mcutils::KroneckerDelta(moshinsky_convention,false))/A;
 
     for(int N=0; N<=Nmax; N++)
       for(int S=0; S<=1; ++S)
@@ -280,7 +281,7 @@ namespace u3shell
   void Tintr(int Nmax,u3shell::RelativeRMEsU3ST& Tintr, int A, double hbar_omega, double coef=1.0, bool moshinsky_convention=false)
   {
     // Something weird with coef factor.  Won't work for moshinsky_convention=true
-    coef*=hbar_omega/(4*(KroneckerDelta(moshinsky_convention,false)));
+    coef*=hbar_omega/(4*(mcutils::KroneckerDelta(moshinsky_convention,false)));
     k2intr(Nmax,Tintr, A, coef, moshinsky_convention);
   }
 
@@ -291,7 +292,7 @@ namespace u3shell
       // Read in the interaction from file
       basis::RelativeSpaceLSJT relative_space_lsjt(Nmax, Jmax);
       std::array<basis::RelativeSectorsLSJT,3> isospin_component_sectors_lsjt;
-      std::array<basis::MatrixVector,3> isospin_component_matrices_lsjt;
+      std::array<basis::OperatorBlocks<double>,3> isospin_component_matrices_lsjt;
       basis::RelativeOperatorParametersLSJT operator_labels;
       basis::ReadRelativeOperatorLSJT(
         interaction_filename,relative_space_lsjt,operator_labels,

programs/rotor/triaxial_rotor.cpp

@@ -35,7 +35,7 @@ namespace rotor
     A.resize(3);
     for(int k=1; k<=3; ++k)
       {
-        A[k-1]=alpha/sqr(sin(gamma-2*PI/3*k));
+        A[k-1]=alpha/mcutils::sqr(sin(gamma-2*PI/3*k));
       }
   }
 
@@ -71,18 +71,18 @@ namespace rotor
       GetD(x,lambda,D);
 
       //First coefficient      
-      coefs[0]=lambda[1]*lambda[2]/(2*sqr(lambda[0])+lambda[1]*lambda[2])*A[0];
-      coefs[0]=coefs[0]+lambda[0]*lambda[2]/(2*sqr(lambda[1])+lambda[0]*lambda[2])*A[1];
-      coefs[0]=coefs[0]+lambda[0]*lambda[1]/(2*sqr(lambda[2])+lambda[0]*lambda[1])*A[2];
+      coefs[0]=lambda[1]*lambda[2]/(2*mcutils::sqr(lambda[0])+lambda[1]*lambda[2])*A[0];
+      coefs[0]=coefs[0]+lambda[0]*lambda[2]/(2*mcutils::sqr(lambda[1])+lambda[0]*lambda[2])*A[1];
+      coefs[0]=coefs[0]+lambda[0]*lambda[1]/(2*mcutils::sqr(lambda[2])+lambda[0]*lambda[1])*A[2];
 
       //Second coefficient 
-      coefs[1]=lambda[0]/(2*sqr(lambda[0])+lambda[1]*lambda[2])*A[0];
-      coefs[1]=coefs[1]+lambda[1]/(2*sqr(lambda[1])+lambda[0]*lambda[2])*A[1];
-      coefs[1]=coefs[1]+lambda[2]/(2*sqr(lambda[2])+lambda[0]*lambda[1])*A[2];
+      coefs[1]=lambda[0]/(2*mcutils::sqr(lambda[0])+lambda[1]*lambda[2])*A[0];
+      coefs[1]=coefs[1]+lambda[1]/(2*mcutils::sqr(lambda[1])+lambda[0]*lambda[2])*A[1];
+      coefs[1]=coefs[1]+lambda[2]/(2*mcutils::sqr(lambda[2])+lambda[0]*lambda[1])*A[2];
 
-      coefs[2]=1/(2*sqr(lambda[0])+lambda[1]*lambda[2])*A[0];
-      coefs[2]=coefs[2]+1/(2*sqr(lambda[1])+lambda[0]*lambda[2])*A[1];
-      coefs[2]=coefs[2]+1/(2*sqr(lambda[2])+lambda[0]*lambda[1])*A[2];
+      coefs[2]=1/(2*mcutils::sqr(lambda[0])+lambda[1]*lambda[2])*A[0];
+      coefs[2]=coefs[2]+1/(2*mcutils::sqr(lambda[1])+lambda[0]*lambda[2])*A[1];
+      coefs[2]=coefs[2]+1/(2*mcutils::sqr(lambda[2])+lambda[0]*lambda[1])*A[2];
 
 
     }
@@ -132,7 +132,7 @@ namespace rotor
                 // std::cout<<"    "
                 // <<am::Unitary6J(L,m,Lbar,m,L,0)<<"  "
                 // <<ParitySign(L+Lbar+m)*Hat(Lbar)/Hat(L)/Hat(m)<<std::endl;
-                // // <<sqr(am::Unitary6J(L,1,L,2,Lbar,m))<<"  "
+                // // <<mcutils::sqr(am::Unitary6J(L,1,L,2,Lbar,m))<<"  "
                 // std::cout<<"SU3 coefs "
                 // <<L<<"  "<<Lbar<<"  "<<m<<std::endl
                 // <<kappa<<" "<<kappa_bar<<"  "<<kappap<<std::endl
@@ -142,7 +142,7 @@ namespace rotor
 
                 rme+=ParitySign(m)*Hat(m)
                       *am::Unitary6J(L,m,Lbar,m,L,0)*L*(L+1)
-                      *sqr(am::Unitary6J(L,1,L,2,Lbar,m))
+                      *mcutils::sqr(am::Unitary6J(L,1,L,2,Lbar,m))
                       *u3::W(x,kappa_bar,Lbar,u3::SU3(1,1),1,2,x,kappap,L,1)
                       *u3::W(x,kappa,L,u3::SU3(1,1),1,2,x,kappa_bar,Lbar,1);
                 }