group_baryon_operator_w.h

Go to the documentation of this file.

// -*- C++ -*-
/*! \file
 *  \brief Construct group baryon operators
 */
 
#ifndef __group_baryon_operator_w_h__
#define __group_baryon_operator_w_h__
 
#include "handle.h"
#include "meas/hadron/baryon_operator.h"
#include "meas/smear/quark_smearing.h"
#include "io/xml_group_reader.h"
#include "io/qprop_io.h"
#include "io/cfgtype_io.h"
#include "meas/hadron/baryon_operator_aggregate_w.h"
#include "meas/hadron/baryon_operator_factory_w.h"
#include "meas/sources/source_smearing_aggregate.h"
#include "meas/sources/source_smearing_factory.h"
#include "meas/sinks/sink_smearing_aggregate.h"
#include "meas/sinks/sink_smearing_factory.h"
#include "meas/smear/quark_source_sink.h"
#include "meas/sources/dilutezN_source_const.h"
#include "meas/sources/zN_src.h"
 
namespace Chroma
{
//#define MAKE_SOURCE_OPERATORS
#define MAKE_SINK_OPERATORS
#define PARALLEL_RUN
//#define REDUCETOTIMEDILUTION 1
 
  //! Name and registration
  /*! @ingroup hadron */
  namespace GroupBaryonOperatorEnv
  {
    bool registerAll();
    //! Group baryon operator parameters
    /*! @ingroup hadron */
        
          struct BaryonOperator_t
          {
            //! Serialize generalized operator object
            multi1d<Complex> serialize();
            //! Baryon operator
            struct BaryonOperatorInsertion_t
            {
              //! Possible operator index
              struct BaryonOperatorIndex_t
              {
                //! Baryon operator element
                struct BaryonOperatorElement_t
                {
                  multi2d<DComplex> elem; 
                };
                multi1d<BaryonOperatorElement_t> ind;
              };
              multi3d<BaryonOperatorIndex_t> op; 
            };
        Seed seed_l;            /*!< Id of left quark */
        Seed seed_m;            /*!< Id of middle quark */
        Seed seed_r;            /*!< Id of right quark */
 
      multi1d< multi1d<int> > quarkOrderings; // [ index ][ LMR ]=0,1,2
        int mom2_max;
        int j_decay;
                        int version;
            multi1d<BaryonOperatorInsertion_t> orderings; 
          };
                
                int DilSwap( int ord, int i, int j, int k, int which );
                int DilSwapInv( int ord, int i, int j, int k, int which );
 
    struct Params
    {
      Params();
      Params( XMLReader& in, const std::string& path );
      void writeXML( XMLWriter& in, const std::string& path ) const;
    
                        PropSourceSmear_t source_smearing;
        PropSinkSmear_t sink_smearing;
 
      struct param_t
                        {
                                std::string baryon_operator; // NUCLEON, DELTA
                        } 
                        param;
                        
                        std::string name; // ?
      int j_decay;  
                        
      std::string InputFileName; // Input generated by "Gen_Input_For_All2All_Baryons2.pl"
      multi1d<std::string> Names; // Names of operators ... "G1g_L3_TDT_15"
      multi1d<int> nrow; // lattice size Nx Ny Nz Nt
      int Nmomenta, mom2_max; 
      int Noperators; // total number of baryon operators to be constructed
      int NQQQs;      // total number of QQQ's
      int NsrcOrderings, NsnkOrderings; 
      multi1d< multi1d<int> > SrcOrderings; // [ index ][ LMR ]=0,1,2
      multi1d< multi1d<int> > SnkOrderings; // [ index ][ LMR ]=0,1,2
                                                
      // the hybrid-list (dilution) sizes
      multi1d< multi1d<int> > NH; // NH[ ord ][ 0,1,2 ] 
                                
                        struct Qprop_t // Prop_t
                        {
                          struct Solutions_t // Operator_t
                          {
                                        multi1d<std::string> soln_file_names; // soln_files
                                        Seed seed;
                          };
                          multi1d<Solutions_t> solns; // quark index (0,1,2)  op
        std::string op_file;
                        } 
                        qprop; // prop
                        
                //! Structure holding a source and its solutions
                struct QuarkSourceSolutions_t
                {
                  //! Structure holding solutions
                  struct QuarkSolution_t
                  {
                    LatticeFermion source;
                    LatticeFermion soln;
                    PropSourceConst_t source_header;
                    ChromaProp_t prop_header;
                  };
                  int j_decay;
                  Seed seed;
                  multi1d<QuarkSolution_t> dilutions;
                };
                        
                        struct gaugestuff_t
                        {
                                Cfg_t cfg;
                                multi1d<LatticeColorMatrix> u;
                                GroupXML_t link_smearing;
                        } 
                        gaugestuff;
                        
                        struct dilution_t
                        {
                                int N; /*!< Z(N) */
                                multi1d<int> spatial_mask_size; /*!< Spatial size of periodic mask */
                                multi1d< multi1d<int> > spatial_mask; /*!< Sites included in site mask */
                                multi1d<int> color_mask; /*!< Color size of periodic mask */
                                multi1d<int> spin_mask; /*!< Spin size of periodic mask */
                                int j_decay;  // is 3
                                int t_source; // assuming time-dilution
                                Seed ran_seed;
                        };
                        
                        multi1d<dilution_t> dilution;
                        
                        struct NamedObject_t
        {
          QuarkSourceSolutions_t vecs;
                                Qprop_t qprop;
                                BaryonOperator_t Bops;
          std::string gauge_id;
        } 
                        named_obj;                                              
 
                        int NdilReduce;
                        
    }; // end struct Params
                
                
          //! Structure holding a source and its solutions
          struct QuarkSourceSolutions_t
          {
            //! Structure holding solutions
            struct QuarkSolution_t
            {
              LatticeFermion source;
              LatticeFermion soln;
              PropSourceConst_t source_header;
              ChromaProp_t prop_header;
            };
            int j_decay;
            Seed seed;
            multi1d<QuarkSolution_t> dilutions;
          };
 
 
    class GroupBaryonOp //: public BaryonOperator<LatticeFermion>
    {
      public:
                                GroupBaryonOp()
        {}
                                ~GroupBaryonOp()
                                {}
                                //! Baryon operator
                        BaryonOperator_t baryonoperator;
                                
                                std::string name;
        Params myparams;
    }; // end class GroupBaryonOp
 
 
    /*! @ingroup hadron
     *
     * Create a group theoretical construction baryon
     */
    class GroupBaryonQQQ : public BaryonOperator<LatticeFermion>
    {
      public:
                void init( const Params& p );
        //! Full constructor
        GroupBaryonQQQ( const Params& p );
        ~GroupBaryonQQQ(){}
                                // ?
                                std::string name;                       
 
        //! Compute the operator
        multi1d<LatticeComplex> operator() ( const LatticeFermion& quark1,
                                             const LatticeFermion& quark2,
                                             const LatticeFermion& quark3,
                                                                                                                                                                                 int* qindices,
                                             enum PlusMinus isign 
                                           ) const;
 
        LatticeComplex operator() ( const LatticeFermion& quark1,
                                    const LatticeFermion& quark2,
                                    const LatticeFermion& quark3,
                                                                                                                                                int* qindices
                                  ) const;
        //! Compute the operator
        multi1d<LatticeComplex> operator() ( const LatticeFermion& quark1,
                                             const LatticeFermion& quark2,
                                             const LatticeFermion& quark3,
                                             enum PlusMinus isign 
                                           ) const;
 
        LatticeComplex operator() ( const LatticeFermion& quark1,
                                    const LatticeFermion& quark2,
                                    const LatticeFermion& quark3
                                  ) const;
 
        struct QuarkTerm_t
        {
          int displacement;    /*!< Orig plus/minus 1-based directional displacements */
          int spin;            /*!< 0-based spin index */
          int disp_dir;        /*!< 0-based direction */
          int disp_len;        /*!< 0-based length */
          int disp_ind;        // disp index 0=0, 1=1, 2=2, 3=3, 4=-1, 5=-2, 6=-3
        };
        multi1d<QuarkTerm_t> quark; /*!< Displacement and spin for each quark */
        
        // Number of Baryon operators requiring this QQQ
        int NBaryonOps;  // this is the size of the arrays below
        multi1d<DComplex> coef;      // coef[ local_BOpIndex ]
        multi1d<int> whichBaryonOps; // whichBaryonOps[local_BopIndex] of size NBaryonOps
        multi1d<std::string> Names;  // and their names
        multi1d<GroupBaryonOp*> baryon; // pointer to the GroupBaryonOp array
 
        //private: 
        //! Hide partial constructor
        GroupBaryonQQQ();
        Params myparams;
                                SpinMatrix spin_rotate_mat;
 
        //protected:
        //Handle< QuarkSourceSink<LatticeFermion> > sourceSmearing;
        //Handle< QuarkSourceSink<LatticeFermion> > sinkSmearing;
 
        //! Construct array of std::maps of displacements
        void displaceQuarks( multi1d< std::map<int, LatticeFermion> >& disp_quarks,
                             const multi1d<LatticeFermion>& q,
                                                                                                                 int* qindices
                           ) const;
        //! First smear then displace the quarks
        void smearDisplaceQuarks( multi1d< std::map<int, LatticeFermion> >& disp_quarks,
                                  const LatticeFermion& q1,
                                  const LatticeFermion& q2,
                                  const LatticeFermion& q3,
                                                                                                                                        int* qindices
                                ) const;
        //! Manipulate the quark fields
        void quarkManip( multi1d< std::map<int, LatticeFermion> >& disp_quarks,
                         const LatticeFermion& q1,
                         const LatticeFermion& q2,
                         const LatticeFermion& q3,
                                                                                                 int* qindices
                       ) const;
        //! The spin basis matrix to goto Dirac
        const SpinMatrix& rotateMat() const
        {
          return spin_rotate_mat;
        }
                                LatticeFermion operator() ( Params::dilution_t ) const;
 
    }; // end class GroupBaryonQQQ
                
    //
    // read in input parameters from an ASCII file
    //
    void ReadTextInput( Params&, multi1d<GroupBaryonOp>&, multi1d<GroupBaryonOp>&,
                        multi1d<GroupBaryonQQQ>&, multi1d<GroupBaryonQQQ>& );
                
        //! Reader
        /*! @ingroup hadron */
 
        void read( XMLReader& xml, const std::string& path, 
                           GroupBaryonOperatorEnv::Params::Qprop_t::Solutions_t& input );
        void read( XMLReader& xml, const std::string& path, 
                           GroupBaryonOperatorEnv::Params::Qprop_t& input );
        void read( XMLReader& xml, const std::string& path, 
                           GroupBaryonOperatorEnv::Params& input );
 
  }  // end namespace GroupBaryonOperatorEnv
 
  //! Writer
  /*! @ingroup hadron */
  void write( XMLWriter& xml, const std::string& path, 
                    const GroupBaryonOperatorEnv::Params& param );
 
}  // end namespace Chroma
 
//
// Perl script to generate the input file
//
/*
#! /usr/bin/perl
#------------------------------------------------------------------------
# The output file of this script ... the name is kind of important 
# because the xml input file is not ready to go yet and hence the name 
# will be used to identify the record of the inputs
$OutputFileName = "Chroma_Input.txt"; 
#------------------------------------------------------------------------
# Some other files that will be generated ...
$SortedListFileName = "SortedChannelList.txt";
       $ListOfNames = "NamesOfOperators.txt";
#------------------------------------------------------------------------
$ProjCoeffDir = "$ENV{HOME}/src2/work2/projection_coefficients/Nucleons";
#------------------------------------------------------------------------
# Input file for operator indices (note that the range op - doesn't work)
# an example is given below ...
$OperatorIndicesFile = "from_point2all.in";
#------------------------------------------------------------------------
# Lattice Size
@Nsize = ( 12, 12, 12, 48 ); 
$Nt = $Nsize[ 3 ];
#------------------------------------------------------------------------
# Hybrid list size
for(my $q=1; $q<=3; ++$q) 
{
        $NTdil[$q] = $Nt; $NCdil[$q] = 1; $NSdil[$q] = 1; $NXdil[$q] = 1;
        #$NTdil[$q] = $Nt; $NCdil[$q] = 3; $NSdil[$q] = 4; $NXdil[$q] = 1;
}
for(my $q=1; $q<=3; ++$q){ $NH[$q]=$NTdil[$q]*$NCdil[$q]*$NSdil[$q]*$NXdil[$q]; }
#
@Ndil = ( $NH[1], $NH[2], $NH[3] );
#------------------------------------------------------------------------
# Different displacement lengths
@DispLengths = ( 3 );
# all the channels
@Channels = ( 'G1g', 'G1u', 'G2g', 'G2u', 'Hg', 'Hu' );
#------------------------------------------------------------------------
# Solution file names  arg: stub, time-dil, colour-dil, spin-dil, space-dil, ".lime"
$qstub = "wl_6p1_12_48_m0p305_nu0p902"; $qext = ".lime"; 
@qpropFileNames = ();
for(my $q=1; $q<=3; ++$q)
{
        @temp = &make_qprop_names( "${qstub}${q}",$NTdil[$q],$NCdil[$q],$NSdil[$q],$NXdil[$q],$qext );
        push( @qpropFileNames, @temp );
}
#------------------------------------------------------------------------
# Nucleon correlation function is re-written as,
#
#    (N+)         (I)   ijk               _ijk     _ijk      _ijk               
#   C        =   c     B      (t) {     2 B  ___ - B  ___  - B  ___             
#    IJ[012]      mnt   I<mnt>             J<mnt>   J<mtn>    J<nmt>            
#                                         _ijk     _ijk      _ijk               
#                                   + ( 2 B  ___ - B  ___  - B  ___  )          
#                                          J<tnm>   J<tmn>    J<ntm>            
#                                       _(J)                                                                                                                       
#                                 }(t0) c___                                                                                                                       
#                                        mnt                                                                                                                       
#------------------------------------------------------------------------
# Input File Example: (from_point2all.in)
#
#    <G1>
#        SS  0  1  2
#        SD  0 11 17 20 22
#        DDI 0  4  5  9 12
#        DDL 1  4 10 15 21
#        TDT 3  5  9 11 25
#        <G2>
#        SS   
#        SD   0  1  2  6  7
#        DDI  1  4  5  6  7
#        DDL 32 37 41 52 61
#        TDT  1 33 45 51 61
#        <H>
#        SS   0
#        SD   7  9 10  13  31
#        DDI  8 15 17  23  31
#        DDL 47 54 84 113 124
#        TDT 35 71 86  95 104
#
#------------------------------------------------------------------------
#                                       Nothing from here needs to be touched
#------------------------------------------------------------------------
#
# Permutations of noise indices
#
$NumSrcPerm = 6; # number of perm of 0,1,2 for SOURCE
        @srcOrder0 = ( 0, 1, 2 ); $SrcPerm{0} = "srcOrder0";
        @srcOrder1 = ( 2, 1, 0 ); $SrcPerm{1} = "srcOrder1";
        @srcOrder2 = ( 0, 2, 1 ); $SrcPerm{2} = "srcOrder2"; 
        @srcOrder3 = ( 1, 0, 2 ); $SrcPerm{3} = "srcOrder3";
        @srcOrder4 = ( 1, 2, 0 ); $SrcPerm{4} = "srcOrder4";
        @srcOrder5 = ( 2, 0, 1 ); $SrcPerm{5} = "srcOrder5";
$NumSnkPerm = 1; # number of perm of 0,1,2 for SINK
        @snkOrder0 = ( 0, 1, 2 ); $SnkPerm{0} = "snkOrder0";
#------------------------------------------------------------------------
@ListOfIrreps  = ( 'G1', 'G2', 'H' ); # in the input file
@OperatorTypes = ( 'SS', 'SD', 'DDI', 'DDL', 'TDT' );
$LongForm{'SS'}  = "Single_Site";
$LongForm{'SD'}  = "Singly_Displaced";
$LongForm{'DDI'} = "Doubly_Displaced_I";
$LongForm{'DDL'} = "Doubly_Displaced_L";
$LongForm{'TDT'} = "Triply_Displaced_T";
$QNS{'G1g'} = 'G1'; $QNS{'G1u'} = 'G1';
$QNS{'G2g'} = 'G2'; $QNS{'G2u'} = 'G2';
 $QNS{'Hg'} = 'H';   $QNS{'Hu'} = 'H';
#------------------------------------------------------------------------
#
# First get all of the operator indices
#
if ( 1 ) {
  $Noperators = 0;
  foreach $irrep ( @ListOfIrreps ) {
    $name = $irrep.'ops'; # ex. G1ops
    &GetOperatorList( $OperatorIndicesFile, ListOfIrreps, $irrep, $name );
    foreach $type ( @OperatorTypes ) {
      $name2 = $name.'List_'.$type; # ex. G1opsList_SS
      @$name2 = &GetNums( $$name{"$type"} ); # ex. $G1ops{"SS"}
      foreach $channel ( @Channels ) {
                                if ( ($QNS{$channel} eq $irrep) && (@$name2 ne {}) ) {
                                        if ( $type eq 'SS' ) {
                                                $Noperators += (scalar(@$name2));
                                        } else {
                                                $Noperators += (scalar(@$name2)) * (scalar(@DispLengths));
        } } } } }
}
#
# Generate the hash table
#
if ( 1 ) {
  open(HASH,">Hashes");
  foreach my $QN ( @Channels ) 
        {
        $QNS = $QNS{"$QN"};
          foreach my $OperatorType ( @OperatorTypes ) 
                {
      $InputFileName = "${ProjCoeffDir}/${QN}/$LongForm{\"$OperatorType\"}";
      $ListName = "${QNS}opsList_${OperatorType}";
      foreach my $OperatorIndx ( @${ListName} ) 
                        {
              $Nterms = 0;
        ( $Nterms, %list ) = &ReadProjCoeffFile( $InputFileName, $OperatorIndx );
              foreach my $DispLength ( @DispLengths ) 
                                {
                #
                # SS L1
                #
                if ( (${OperatorType}eq'SS')and(${DispLength}==${DispLengths[0]}) ) 
                                        {
                    for (my $i=0; $i<$Nterms; ++$i) 
                                                        {
                @IJK = &GetNums( $list{$i} );
                      $line = '';
                # DO NOT use the coefficients here
                for (my $j=0; $j<($#IJK-1); ++$j) { $line = $line . "$IJK[$j] "; }
                      $line = $line . "0";
                      #$line = $line . ${DispLength};
                      $hash = &HashCode($line);
                      printf(HASH "${hash}\n");
                    }
                }
                #
                # SD L1 L2 : DDI L1 L2 : DDL L1 L2 : TDT L1 L2
                #
                if ( ${OperatorType} ne 'SS' ) 
                                        {
                    for (my $i=0; $i<$Nterms; ++$i) 
                                                        {
                @IJK = &GetNums( $list{$i} );
                      $line = '';
                # DO NOT use the coefficients here
                for (my $j=0; $j<($#IJK-1); ++$j) { $line = $line . "$IJK[$j] "; }
                      $line = $line . ${DispLength};
                      $hash = &HashCode($line);
                      printf(HASH "${hash}\n");
                    }
                }
              } # end foreach $DispLength ( @DispLengths )
            } # end foreach $OperatorIndx ( @${ListName} )
          } # end foreach $OperatorType ( @OperatorTypes )
        } # end foreach $QN ( @Channels )
  close(HASH);
 
        system( "sort Hashes | uniq > Hashes2" );
  open(FILE,"<Hashes2");
        @UniqHashes = ();
        while (<FILE>) { chop($_); push( @UniqHashes, $_ ); }
  close(FILE);
        system("rm Hashes Hashes2");
}
#
# Main Loop over Irreps
#
foreach $QN ( @Channels ) 
{ 
        &MakeOpTable($QN); 
}
#
# Sorted List (according to channel name)
#   "Name re im # s1 s2 s3 d1 d2 d3 dL"
#
if ( 1 ) {
        open(FILE,">$SortedListFileName");
  open(FILENAMES,">$ListOfNames");
  foreach my $QN ( @Channels ) 
        {
    $QNS = $QNS{"$QN"};
    foreach my $OperatorType ( @OperatorTypes ) 
                {
      $ListName = "${QNS}opsList_${OperatorType}";
      $InputFileName = "${ProjCoeffDir}/${QN}/$LongForm{\"$OperatorType\"}";
      foreach my $OperatorIndx ( @${ListName} ) 
                        {
        ( $Nterms, %list ) = &ReadProjCoeffFile( $InputFileName, $OperatorIndx );
        for (my $i=0; $i<$Nterms; ++$i) 
                                { 
          @IJK = &GetNums( $list{$i} );
          foreach my $DispLength ( @DispLengths ) 
                                        {
            $line = '';
            for (my $j=0; $j<($#IJK-1); ++$j) { $line=$line."$IJK[$j] "; }
            if(($OperatorType eq 'SS')&&(${DispLength}==${DispLengths[0]})) 
                                                { 
                                                        $line=$line." 0"; 
            } else { 
                                                        $line=$line." ${DispLength}"; 
                                                }
            $hash = &HashCode($line);
            for (my $r=0; $r <= $#{$IrrepsInHash{$hash}}; $r=$r+3) 
                                                {
              $newline = "${$IrrepsInHash{$hash}}[$r] ${$IrrepsInHash{$hash}}[$r+1] ${$IrrepsInHash{$hash}}[$r+2] ${line}";
                                                        printf(FILE "$newline\n");
                                                        printf(FILENAMES "${$IrrepsInHash{$hash}}[$r]\n");
  } } } } } }
  close(FILE);
  close(FILENAME);
  system("sort ${SortedListFileName} | uniq > mytemp");
  system("mv mytemp ${SortedListFileName}");
  system("sort ${ListOfNames} | uniq > myDistinctOperators");
        open(FILE,"<myDistinctOperators");
        open(OUTFILE,">DistinctOperators");
  $GlobalNumbering = 0;
        while (<FILE>) {
                chop($_);
                printf(OUTFILE "%-3i %s\n",$GlobalNumbering,$_);
                $OperatorNames{"$_"} = $GlobalNumbering;
                $GlobalNumbering++;
        }
        close(FILE);
        close(OUTFILE);
        system("wc -l DistinctOperators > mytemp2");
        system("cat mytemp2 DistinctOperators > ${ListOfNames}");
        system("rm -f mytemp2 myDistinctOperators DistinctOperators");
}
#
# Output Results to File
#   "# s1 s2 s3 d1 d2 d3 dL NtermsNeedingThisComb"
#   "index" corresponding to "Name" (check NamesOfOperators.txt)
#   "re im "
#
if ( 1 ) {
  open(FILE,">${OutputFileName}");
        printf(FILE "${Nsize[0]} ${Nsize[1]} ${Nsize[2]} ${Nsize[3]}\n");
        printf(FILE "$NumSrcPerm\n");
        for (my $i=0; $i < $NumSrcPerm; ++$i)
        {
                printf(FILE "${$SrcPerm{$i}}[0] ${$SrcPerm{$i}}[1] ${$SrcPerm{$i}}[2]\n");
        }
        # added so that XML and TXT input files are more similar
        printf(FILE "$NumSnkPerm\n");
        for (my $i=0; $i < $NumSnkPerm; ++$i)
        {
                printf(FILE "${$SnkPerm{$i}}[0] ${$SnkPerm{$i}}[1] ${$SnkPerm{$i}}[2]\n");
        }
        printf(FILE "${Ndil[0]} ${Ndil[1]} ${Ndil[2]}\n");
        printf(FILE "${Noperators}\n");
        # moved here so that XML and TXT input files are more similar
        open(INFILE,"<${ListOfNames}");
                $line = <INFILE>;
                while ( <INFILE> ) 
                {
                        printf(FILE "$_");
                }
                close(INFILE);
        #
        $NElemOps = scalar(@UniqHashes);
        printf(FILE "${NElemOps}\n");
        foreach $hash ( @UniqHashes ) 
        {
                $line = $InvHashCode{$hash}; # all 7 indices
                @indices = &GetNums( $line );
                printf(FILE "${hash} ");
                PrintList2FILE( FILE, @indices );
                my $N = ( scalar(@{$IrrepsInHash{$hash}}) ) / 3;
                printf(FILE "${N}\n"); # number of ops which need this qqq
                for (my $r=0; $r <= $#{$IrrepsInHash{$hash}}; $r=$r+3) 
                {
                 #printf(FILE                "${$IrrepsInHash{$hash}}[$r]\n");
                        printf(FILE "$OperatorNames{${$IrrepsInHash{$hash}}[$r]}\n");
                        printf(FILE "${$IrrepsInHash{$hash}}[$r+1]");
                        printf(FILE " ${$IrrepsInHash{$hash}}[$r+2]\n");
                }
        }
        for (my $i=0; $i <= $#{qpropFileNames}; ++$i) 
        {
                printf(FILE "$qpropFileNames[$i]\n");
        }
        close(FILE);
}
 
#===================================================================================
 
# name s1 s2 s3 d1 d2 d3  h
#   0   1  2  3  4  5  6  7  (use with "sort list_by_hash_value keys %hasharray")
sub list_by_hash_value { ($shash{$a}[7]) <=> ($shash{$b}[7]); };
sub PrintList{ foreach $i (@_) { printf("%2s ",$i); } print "\n";};
sub PrintList2FILEN
{ local *HANDLE = shift;
        foreach$i(@_){printf(HANDLE "%2s ",$i);}printf(HANDLE "\n");
}
sub PrintList2FILE 
{ local *HANDLE = shift;
        foreach$i(@_){printf(HANDLE "%2s ",$i);}
}
sub HashCode
{ my ( $line ) = @_;
  @in = &GetNums(@_); # s1 s2 s3 d1 d2 d3 dL
        $hash_code = ( ( $in[0] - 1 ) << 17 ) # spins from 1
             + ( ( $in[1] - 1 ) << 15 )
             + ( ( $in[2] - 1 ) << 13 )
             + ( ( $in[3] + 3 ) << 10 ) # displ dirs from -3
             + ( ( $in[4] + 3 ) << 7 )
             + ( ( $in[5] + 3 ) << 4 )
             + (   $in[6] );            # disp length index from 1
        $InvHashCode{$hash_code} = "$in[0] $in[1] $in[2] $in[3] $in[4] $in[5] $in[6]";
        $GetHashCode{$line} = $hash_code;
        return($hash_code);
}
sub MakeOpTable
{ # Coefficient File: ${ProjCoeffDir}/${QN}/$LongForm{\"$OperatorType\"}
        my($irrep,$type,$name,$name2,$IndexName,$ListName,$InputFileName);
        my($Nterms,$line,$hash);
        $QN = $_[0];
        $QNS = $QNS{"$QN"};
        print "$QNS $QN\n";
  foreach my $DispLength ( @DispLengths ) {
    foreach my $OperatorType ( @OperatorTypes ) {
      $InputFileName = "${ProjCoeffDir}/${QN}/$LongForm{\"$OperatorType\"}";
      $ListName = "${QNS}opsList_${OperatorType}";
                        if ( ${OperatorType} eq 'SS' ) {
                                if ( ${DispLength} == ${DispLengths[0]} ) {
                            foreach my $OperatorIndx ( @${ListName} ) {
            $Nterms = 0;
                              ( $Nterms, %list ) = &ReadProjCoeffFile($InputFileName, $OperatorIndx);
            $IndexName = "${OperatorType}_${OperatorIndx}";
            for (my $i=0; $i<$Nterms; ++$i) { 
              @IJK = &GetNums( $list{$i} );
              $line = '';
              for (my $j=0; $j<($#IJK-1); ++$j) { 
                $line = $line . "$IJK[$j] ";
              }
                                $line = $line . "0";
                                #$line = $line . ${DispLength};
              $hash = &HashCode( $line );
                                                        push( @{$IrrepsInHash{$hash}}, "${QN}_L0_${IndexName}" );
                                                        #push( @{$IrrepsInHash{$hash}}, "${QN}_L${DispLength}_${IndexName}" );
                                                        push( @{$IrrepsInHash{$hash}}, $IJK[($#IJK-1)] );
                                                        push( @{$IrrepsInHash{$hash}}, $IJK[($#IJK)] );
            } # end of going thru terms
                            } # end foreach OperatorIndx
                                } # end if (${DispLength} == ${DispLengths[0]})
                        } # end if SS 
                        if ( ${OperatorType} ne 'SS' ) {
                          foreach my $OperatorIndx ( @${ListName} ) {
          $Nterms = 0;
                            ( $Nterms, %list ) = &ReadProjCoeffFile($InputFileName, $OperatorIndx);
          $IndexName = "${OperatorType}_${OperatorIndx}";
          for (my $i=0; $i<$Nterms; ++$i) { 
            @IJK = &GetNums( $list{$i} );
            $line = '';
            for (my $j=0; $j<($#IJK-1); ++$j) { 
              $line = $line . "$IJK[$j] ";
            }
            $line = $line . "$DispLength";
            $hash = &HashCode( $line );
                                                push( @{$IrrepsInHash{$hash}}, "${QN}_L${DispLength}_${IndexName}" );
                                                push( @{$IrrepsInHash{$hash}}, $IJK[($#IJK-1)] );
                                                push( @{$IrrepsInHash{$hash}}, $IJK[($#IJK)] );
          } # end of $i<Nterms loop
                          } # end of : foreach OperatorIndx
                        } # end of : ne SS
    } # end OperatorType loop
  } # end DispLength loop
};
sub ReadProjCoeffFile 
{ my ( $FileName, $TargetIndx ) = @_;
        my($NumOps,$OpIndx,$Nterms_,$NTerms,$line,$k,$len);
  open( PROJCOEF, "< $FileName" );
  $NumOps = <PROJCOEF>; chop($NumOps); $NumOps*=1;
  $OpIndx = -1; $Nterms_ = 0; $NTerms = 0; %mylist = ();
  while( <PROJCOEF> )
  { chomp;
    $len = split;
    if( $len == 1 ) {
      ( $Nterms_ ) = split; $OpIndx++;
    } else {
      if( $OpIndx == $TargetIndx ) {
        $line = $_;
        my @entry = &GetNums($line);
                                $mylist{"0"}="$entry[0] $entry[1] $entry[2] $entry[3] $entry[4] $entry[5] $entry[6] $entry[7]";
        for($k=1; $k < $Nterms_; ++$k) {
          $line = <PROJCOEF>; chomp;
          @entry = &GetNums($line);
                                        $mylist{"$k"}="$entry[0] $entry[1] $entry[2] $entry[3] $entry[4] $entry[5] $entry[6] $entry[7]";
        }
                                $NTerms = $Nterms_;
                        }
                }
  } close( PROJCOEF );
        return( $NTerms, %mylist );
};
sub GetOperatorList
{ my ( $filename, $list, $channel, $outputref ) = @_;
        my ( $irrep, $line, $oplist );
        open(FILE,"<$filename");
        foreach $irrep ( @{$list} ) {
        $line = <FILE>; # irrep
    $line = <FILE>; chop($line);
                if ( ($irrep eq "$channel") and ( $line ) ) {
                        ${$outputref}{'SS'} = $line; }
    $line = <FILE>; 
                if ( ($irrep eq "$channel") and ( $line ) ) {
                        ${$outputref}{'SD'} = $line; }
    $line = <FILE>; 
                if ( ($irrep eq "$channel") and ( $line ) ) {
        ${$outputref}{'DDI'} = $line; }
    $line = <FILE>; 
                if ( ($irrep eq "$channel") and ( $line ) ) {
        ${$outputref}{'DDL'} = $line; }
    $line = <FILE>; 
                if ( ($irrep eq "$channel") and ( $line ) ) {
        ${$outputref}{'TDT'} = $line; }
        } close(FILE);
};
sub RmTrail 
{ my($i,$j,$len,$c,$Out);
        $i=0; $len=length($_[0]); $j=$len;
        while($i<$len){ $c=substr($_[0],$j-1,1);
        if($c =~ /[.0-9]/){$Out=substr($_[0],0,$len-$i);$i=$len;}$j--;$i++;}
        $_[0]=$Out;
};
sub GetNums # ASCII separated by space, comma, semi-colon, colon
{ my ( $line ) = @_; # can simplify this routine for ints only ...
  my($k,$char,$i,$j,$num,$prev);
        $k=0;$i=0;$j=0;$num='';$char='a';@ReturnThis=();# $j is counter
        while(($i<=length($line))and($char ne "\n")and($char ne "")){
                $char=substr($line,$i,1);if(($char eq ',')or($char eq ';')or
                ($char eq ':')){if(length($num)>0){&RmTrail($num);
                push(@ReturnThis,$num);}$num='';$j++;}else{if((($char eq '+')
                or($char eq '-'))and(length($num)==0)){$num=$char;$k++;}else
                {if($k>0){if($char=~/[0-9,.,E,e,D,d,\-,+]/){
                $prev=substr($num,$k-1,1);if(not(($char eq '0')and(
                ($prev eq '+')or($prev eq '-')))){if(($char eq 'D')or($prev eq 'd')
                ){$num=$num.'E';$k++;}else{$num=$num.$char;$k++;}}}else{if(
                length($num)>0){&RmTrail($num);push(@ReturnThis,$num);$num='';
                $j++;}}}else{if($char=~/[0-9,.,\-,+]/){$num=$num.$char;$k++;}}}}
        $i++;}
        return(@ReturnThis);
};
sub sign
{ my ( $line ) = @_;
  my $sign_ = 1;
        @spin = &GetNums(@_); # s1 s2 s3
        $count = 0;
  for (my $i=0; $i < 3; ++$i) {
                if ( $spin[ $i ] >= 3 ) { $count++; }
        }
        $odd = ( $count%2 );
        if ( $odd ) {   $sign_ = -1; }
        return( $sign_ );
}
sub make_qprop_names 
{ my ( $stub, $TD, $CD, $SD, $XD, $ext ) = @_;
        # TfSeCfGf
        if( ($CD >= 2)&&($SD >= 2)&&($XD >= 2) )
        { my $i = 0;
                for(my $t=0; $t<$TD; ++$t) {
                for(my $c=0; $c<$CD; ++$c) {
                for(my $s=0; $s<$SD; ++$s) {
                for(my $x=0; $x<$XD; ++$x) {
                        $FileNames[ $i ] = $stub . "_t${t}_c${c}_s${s}_x${x}" . $ext;
                        $i++; }}}}}
        # TfSxCfGf
        if( ($CD >= 2)&&($SD >= 2)&&($XD == 1) )
        { my $i = 0;
                for(my $t=0; $t<$TD; ++$t) {
                for(my $c=0; $c<$CD; ++$c) {
                for(my $s=0; $s<$SD; ++$s) {
                        $FileNames[ $i ] = $stub . "_t${t}_c${c}_s${s}" . $ext;
                        $i++; }}}}
        # TfSeCfGx
        if( ($CD >= 2)&&($SD == 1)&&($XD >= 2) )
        { my $i = 0;
                for(my $t=0; $t<$TD; ++$t) {
                for(my $c=0; $c<$CD; ++$c) {
                for(my $x=0; $x<$XD; ++$x) {
                        $FileNames[ $i ] = $stub . "_t${t}_c${c}_x${x}" . $ext;
                        $i++; }}}}
        # TfSeCxGf
        if( ($CD == 1)&&($SD >= 2)&&($XD >= 2) )
        { my $i = 0;
                for(my $t=0; $t<$TD; ++$t) {
                for(my $s=0; $s<$SD; ++$s) {
                for(my $x=0; $x<$XD; ++$x) {
                        $FileNames[ $i ] = $stub . "_t${t}_s${s}_x${x}" . $ext;
                        $i++; }}}}
        # TfSeCxGx
        if( ($CD == 1)&&($SD == 1)&&($XD >= 2) )
        { my $i = 0;
                for(my $t=0; $t<$TD; ++$t) {
                for(my $x=0; $x<$XD; ++$x) {
                        $FileNames[ $i ] = $stub . "_t${t}_x${x}" . $ext;
                        $i++; }}}
        # TfSxCfGx
        if( ($CD >= 2)&&($SD == 1)&&($XD == 1) )
        { my $i = 0;
                for(my $t=0; $t<$TD; ++$t) {
                for(my $c=0; $c<$CD; ++$c) {
                        $FileNames[ $i ] = $stub . "_t${t}_c${c}" . $ext;
                        $i++; }}}
        # TfSxCxGf
        if( ($CD == 1)&&($SD >= 2)&&($XD == 1) )
        { my $i = 0;
                for(my $t=0; $t<$TD; ++$t) {
                for(my $s=0; $s<$SD; ++$s) {
                        $FileNames[ $i ] = $stub . "_t${t}_s${s}" . $ext;
                        $i++; }}}
        # TfSxCxGx
        if( ($CD == 1)&&($SD == 1)&&($XD == 1) )
        { my $i = 0;
                for(my $t=0; $t<$TD; ++$t) {
                        $FileNames[ $i ] = $stub . "_t${t}" . $ext;
                        $i++; }}
        return( @FileNames );
}
 
1;
# end
*/
#endif