alg_pot2.C

Go to the documentation of this file.

#include<config.h>
CPS_START_NAMESPACE
//--------------------------------------------------------------------
//  CVS keywords
//
//  $Author: chulwoo $
//  $Date: 2006/05/30 20:32:27 $
//  $Header: /space/cvs/cps/cps++/src/alg/alg_pot2/alg_pot2.C,v 1.2 2006/05/30 20:32:27 chulwoo Exp $
//  $Id: alg_pot2.C,v 1.2 2006/05/30 20:32:27 chulwoo Exp $
//  $Name: v5_0_8 $
//  $Locker:  $
//  $RCSfile: alg_pot2.C,v $
//  $Revision: 1.2 $
//  $Source: /space/cvs/cps/cps++/src/alg/alg_pot2/alg_pot2.C,v $
//  $State: Exp $
//
//--------------------------------------------------------------------
// last modification: manke, Feb 1, 2001
// 1. changes Vst --> Vts (for consistency)
// 2. explicitly transfer max_XYZT in pot_arg
// 3. print out potentials V(x,y,z,t) separately for x,y,z
// 4. changed V?? (potential) --> W?? (wilson loop)
//------------------------------------------------------------------
//
// alg_pot.C
//
// AlgPot is derived from Alg and it measures the potential
// into different propagation directions
//
//------------------------------------------------------------------

CPS_END_NAMESPACE
#include <util/qcdio.h>
#include <alg/alg_pot2.h>
#include <util/lattice.h>
#include <util/gjp.h>
#include <util/smalloc.h>
#include <util/vector.h>           // Tr() and ReTr()
#include <util/verbose.h> 
#include <util/error.h>
#include <comms/glb.h>
CPS_START_NAMESPACE

//------------------------------------------------------------------
// Constructor 
//------------------------------------------------------------------
AlgPot2::AlgPot2(Lattice& latt, 
             CommonArg *c_arg,
             PotArg *arg) : 
             Alg(latt, c_arg) 
{
  cname = "AlgPot2";
  char *fname = "AlgPot2(L&,CommonArg*,PotArg*)";
  VRB.Func(cname,fname);

  // Initialize the argument pointer
  //----------------------------------------------------------------
  if(arg == 0)
    ERR.Pointer(cname,fname, "arg");
  alg_pot_arg = arg;

  // Calculate normalization factor for each processor
  //----------------------------------------------------------------

  int total_sites    = GJP.VolNodeSites() * GJP.Xnodes() *
                       GJP.Ynodes() * GJP.Znodes() * GJP.Tnodes();
  int int_norm       = GJP.Colors()*total_sites;

  norm_fac           = 1.0 / ( Float(int_norm) );
  xiB2               = GJP.XiBare()*GJP.XiBare();
}


//------------------------------------------------------------------
// Destructor
//------------------------------------------------------------------
AlgPot2::~AlgPot2() {
  char *fname = "~AlgPot2()";
  VRB.Func(cname,fname);
}





//------------------------------------------------------------------
// run()
//   check is used to set r2max
//------------------------------------------------------------------
void AlgPot2::run()
{
#if TARGET==cpsMPI
  using MPISCU::fprintf;
#endif
  
  char *fname = "run()";
  VRB.Func(cname,fname);
  
  int sweep    = alg_pot_arg->sweep;
  int prop_dir = alg_pot_arg->prop_dir;
  int check    = alg_pot_arg->check;

  // maximal loop extent in (X,Y,Z,T) direction 
  // running over all the biggest loops would often be very expensive

  int max_ext[4];
 
  int max_X = alg_pot_arg->max_X;
  int max_Y = alg_pot_arg->max_Y;
  int max_Z = alg_pot_arg->max_Z;
  int max_T = alg_pot_arg->max_T;

  max_ext[0] = max_X;
  max_ext[1] = max_Y;
  max_ext[2] = max_Z;
  max_ext[3] = max_T;

  int max_prop = max_ext[prop_dir];

  /*
  // to use pot_arg, r2max is defined here 
  int r2max = 1;
  for (int i = 0; i < 4; i++)
    if( i != prop_dir ) {
      if( r2max <= max_ext[i] ) r2max = max_ext[i];
    }
  r2max = r2max*r2max;
  */
  int r2max = check*check;

  // orthogonal directions to prop_dir 
  int dir_orth[3];
  int idir = 0;
  for (int i = 0; i < 4; i++)
    if( i != prop_dir ) {
      dir_orth[idir] = i;
      idir++;
    }

  int max_1 = max_ext[dir_orth[0]];
  int max_2 = max_ext[dir_orth[1]];
  int max_3 = max_ext[dir_orth[2]];

  int num_loop = 0;
  for (int ext1 = 0; ext1 <= max_1; ext1++){
    for (int ext2 = 0; ext2 <= max_2; ext2++){
      for (int ext3 = 0; ext3 <= max_3; ext3++){  

        int cmax = ext3; 
        if (cmax <= ext2) {cmax = ext2;}
        if (cmax <= ext1) {cmax = ext1;}
        int cmin = ext1;
        if (cmin >= ext2) {cmin = ext2;}
        if (cmin >= ext3) {cmin = ext3;}
          
        int cmid = ext1+ext2+ext3 - cmax - cmin;
          
        int r2 = ext1*ext1 + ext2*ext2 + ext3*ext3;
        if( r2 <= r2max || cmid == 0 ) num_loop+=max_prop;
      }
    }
  }
                              
  Float aniso_factor;
  Float pot_real;
  Float pot_imag;

  // Set the Lattice pointer
  //----------------------------------------------------------------
  Lattice& lat = AlgLattice();

  // static QQ-pair propagates in temporal direction
  // calculate potential for all possible distances
  // on-axis and off-axis
  
  // allocate space for potential and multiplicity
  // actually the following is bigger than necessary
  // as the loops below are only UP TO < MAX_R2
  // its safe for now though
  //  Complex *Wts = new Complex[max_r2];
  //  int *number  = new int[max_r2];
  Complex *Wloop = new Complex[num_loop];
  //  int (*qq_sep)[6];
  //  qq_sep = new int[num_loop][6];
  int *qq_sep0 = new int[num_loop];
  int *qq_sep1 = new int[num_loop];
  int *qq_sep2 = new int[num_loop];
  int *qq_sep3 = new int[num_loop];
  int *qq_sep4 = new int[num_loop];
  int *qq_sep5 = new int[num_loop];


  for (int i = 0; i < num_loop; i++){
    Wloop[i] = 0.;
    qq_sep0[i] = -1;  // x
    qq_sep1[i] = -1;  // y
    qq_sep2[i] = -1;  // z
    qq_sep3[i] = -1;  // t
    qq_sep4[i] =  0;  // prop 
    qq_sep5[i] =  0;  // r2
  }

  int iloop = 0;
  for (int ext_prop=1; ext_prop <= max_prop ; ext_prop++){

    // aniso_factor = xi0^{-2* (number of temporal links) }
    aniso_factor = power(xiB2,-1.*ext_prop);
    
    // vary the extent of static Q\bar Q source along the three
    // different spatial directions
    // there is a certain restrcition here in that we only
    // consider spatial paths like (x,x,x,x,y,y,z,z,z,z,z)
    // but not the Oh-related ones
    //  e.g. not: (x,y,z,z,x,x,y,z,z,x,z)
    
    for (int ext1 = 0; ext1 <= max_1; ext1++){
      for (int ext2 = 0; ext2 <= max_2; ext2++){
        for (int ext3 = 0; ext3 <= max_3; ext3++){  

          int cmax = ext3; 
          int dir_max = 2;
          if (cmax <= ext2) {cmax = ext2; dir_max=1;}
          if (cmax <= ext1) {cmax = ext1; dir_max=0;}
          int cmin = ext1;
          int dir_min = 0;
          if (cmin >= ext2) {cmin = ext2; dir_min=1;}
          if (cmin >= ext3) {cmin = ext3; dir_min=2;}
          
          int cmid = ext1+ext2+ext3 - cmax - cmin;
          int dir_mid = 0+1+2 - dir_max - dir_min;
          
          int r2 = ext1*ext1 + ext2*ext2 + ext3*ext3;
          //      if( r2 <= r2max ) {
          if( r2 <= r2max || cmid == 0 ) {
            
            int check1 = 
              dir_max*dir_max + dir_mid*dir_mid + dir_min*dir_min;
            if ( check1 != 5 ) {
              cout << " sum dir*dir should be 5 not " << check1 << endl;
              exit(-1);
            }

            int cmax2 = 2*cmax;
            int cmid2 = 2*cmid;
            int cmin2 = 2*cmin;
            int chi_mid = cmid2 - cmax;
            int chi_min = cmin2 - cmax;

            // squared distance between poitn and origin
            // int r2 = ext1*ext1 + ext2*ext2 + ext3*ext3;
            // number of links in Wilson Loop
            int nlinks = 2*(ext1+ext2+ext3+ext_prop);
            
            // allocate space for path
            int *path = new int[nlinks]; 
            int ipath = 0;

            // path between Q and Qbar
            for (int i = 1; i <= cmax; i++) {
              path[ipath] = dir_orth[dir_max];
              ipath++;
              if (chi_mid >= 0) {
                chi_mid -= cmax2;
                path[ipath] = dir_orth[dir_mid];
                ipath++;
              }
              if (chi_min >= 0) {
                chi_min -= cmax2;
                path[ipath] = dir_orth[dir_min];
                ipath++;
              }
              chi_mid += cmid2;
              chi_min += cmin2;
            }
            
            // path for propagation
            for (int i = 0; i < ext_prop; i++) {
              path[ipath] = prop_dir;
              ipath++;
            }
            
            // x = 0 , -x = 4 = x + 4
            // y = 1 , -y = 5 = y + 4
            // z = 2 , -z = 6 = z + 4
            // t = 3 , -t = 7 = t + 4
            // inverse path between Q and Qbar
            int qq_path = ext1 + ext2 + ext3;
            int inv_path = qq_path - 1;
            for (int i = 0; i < qq_path; i++) {
              path[ipath] = path[inv_path] + 4;  // inverse direction
              ipath++;
              inv_path--;
            }
            
            // inverse path for propagation
            for (int i = 0; i < ext_prop; i++) {
              path[ipath] = prop_dir + 4;
              ipath++;
            }
            
            // check for path length
            if (ipath != nlinks) {
              cout << "illegal path !!!" << endl;
              cout << "ipath  = " << ipath << endl;
              cout << "nlinks = " << nlinks << endl;
              exit(-1);
            }
            
            // calculated and accumulate the LOOP
            // over all possible origins on the lattice
            // printf("Looping over all sources \n");
            
            int x[4]={0,0,0,0};
            Matrix LOOP;
            //    LOOP.ZeroMatrix();
            Complex pot=0.;
            
            int isite = 0;
            for(x[0]=0;x[0]<GJP.XnodeSites();x[0]++)
              for(x[1]=0;x[1]<GJP.YnodeSites();x[1]++)
                for(x[2]=0;x[2]<GJP.ZnodeSites();x[2]++)
                  for(x[3]=0;x[3]<GJP.TnodeSites();x[3]++){
                    
                    LOOP.ZeroMatrix();
                    lat.PathOrdProdPlus(LOOP, x, path, nlinks);
                    pot += LOOP.Tr();
                    isite++;
                    
                  } // end loop over all local origins
            // accumulate result for each distance and its multiplicity.
            // Notice that some distances can not be
            // represented on the 3D spatial lattice ==
            // some integers can not be represented as
            // a sum of 3 squares, in this case number[r2] remains zero
            Wloop[iloop] = pot;
            int tmp_sep[4] = {-1,-1,-1,-1};
            tmp_sep[dir_orth[0]] = ext1;
            tmp_sep[dir_orth[1]] = ext2;
            tmp_sep[dir_orth[2]] = ext3;
            qq_sep0[iloop] = tmp_sep[0];
            qq_sep1[iloop] = tmp_sep[1];
            qq_sep2[iloop] = tmp_sep[2];
            qq_sep3[iloop] = tmp_sep[3];
            qq_sep4[iloop] = ext_prop;
            qq_sep5[iloop] = r2;
            iloop++;
            
            // deallocate space for path
            delete [] path;
            //    cout << "after delete path" << endl;
          } // end for if(r2<=r2max||cmid==0)
        } // end for ext3
      } // end for ext2
    } // end for ext1
    
  } // end for ext_prop=0; ext_prop <= max_T
  
  if (iloop != num_loop ) {
    cout << "something is wrong about loop counting" << endl;
    cout << "iloop = " << iloop << endl;
    cout << "num_loop = " << num_loop << endl;
    exit(-1);
  }
  
  // dump the averaged potential for all r2 at this timeslice
  char wfname[100];
  sprintf(wfname, "Wloop%04dsme%03d.dat", sweep/1000, sweep%1000);
  FILE *fp= Fopen(wfname, "w");
  //  FILE *fp= Fopen("Wloop_new.dat", "w");
  
  Fprintf(fp, "# max_X= %d max_Y= %d max_Z= %d max_T= %d \n",
          max_X, max_Y, max_Z, max_T);
  Fprintf(fp, "# sweep= %d  prop_dir= %d  num_loop= %d\n",
          sweep/1000, prop_dir, num_loop);
  Fprintf(fp, "# smear= %d \n", sweep%1000);
  Fprintf(fp, "# x  y  z  t prop   pot_real       pot_imag\n");
  
  for (int i = 0; i < num_loop; i++){
    Complex W =  Wloop[i];
    pot_real = real(W);
    pot_imag = imag(W);
    glb_sum(&pot_real) ;
    glb_sum(&pot_imag) ;
    
    // normalise wrt colour and total space-time
    // take into account the anisotropy factor for each temporal link
    pot_real *= aniso_factor*norm_fac;
    pot_imag *= aniso_factor*norm_fac;

    Fprintf(fp," %d  %d  %d  %d  %d    %e   %e  \n",
            //      qq_sep[i][0],qq_sep[i][1],qq_sep[i][2],qq_sep[i][3],
            //      qq_sep[i][4], pot_real, pot_imag);
            qq_sep0[i],qq_sep1[i],qq_sep2[i],qq_sep3[i],
            qq_sep4[i], pot_real, pot_imag);
    //    Fprintf(fp," %d  %d  %d  %d  %d  %d   %e   %e  \n",
    //      qq_sep[i][0],qq_sep[i][1],qq_sep[i][2],qq_sep[i][3],
    //      qq_sep[i][4], qq_sep[i][5], pot_real, pot_imag);
  }
  
  Fclose(fp);
  //  cout << "total_site = " << (1.0 / norm_fac)/3.0 << endl;

  delete[] Wloop;
  //  delete[] qq_sep;
  delete[] qq_sep0;
  delete[] qq_sep1;
  delete[] qq_sep2;
  delete[] qq_sep3;
  delete[] qq_sep4;
  delete[] qq_sep5;

}
CPS_END_NAMESPACE

---