From 0867c8f21f91d55011a84f38d81606ee704faadd Mon Sep 17 00:00:00 2001
From: Yuichi Motoyama <y-motoyama@issp.u-tokyo.ac.jp>
Date: Fri, 22 Mar 2024 17:15:27 +0900
Subject: [PATCH] remove old DLA tools

---
 CMakeLists.txt                    |   25 +-
 src/dla/CMakeLists.txt            |    5 -
 src/dla/generators/CMakeLists.txt |   26 -
 src/dla/generators/boson_B.h      |  127 ---
 src/dla/generators/canonical.h    |   78 --
 src/dla/generators/cfgene.cc      |  140 ----
 src/dla/generators/dla_alg.cc     | 1216 -----------------------------
 src/dla/generators/dla_alg.h      |  356 ---------
 src/dla/generators/exact_B.cc     |  216 -----
 src/dla/generators/exact_H.cc     |  193 -----
 src/dla/generators/lattgene_C.cc  |  211 -----
 src/dla/generators/lattgene_T.cc  |  218 ------
 src/dla/generators/matrix.h       |  618 ---------------
 src/dla/generators/sfgene.cc      |  194 -----
 src/dla/generators/spin_H.h       |  113 ---
 15 files changed, 9 insertions(+), 3727 deletions(-)
 delete mode 100644 src/dla/generators/CMakeLists.txt
 delete mode 100644 src/dla/generators/boson_B.h
 delete mode 100644 src/dla/generators/canonical.h
 delete mode 100644 src/dla/generators/cfgene.cc
 delete mode 100644 src/dla/generators/dla_alg.cc
 delete mode 100644 src/dla/generators/dla_alg.h
 delete mode 100644 src/dla/generators/exact_B.cc
 delete mode 100644 src/dla/generators/exact_H.cc
 delete mode 100644 src/dla/generators/lattgene_C.cc
 delete mode 100644 src/dla/generators/lattgene_T.cc
 delete mode 100644 src/dla/generators/matrix.h
 delete mode 100644 src/dla/generators/sfgene.cc
 delete mode 100644 src/dla/generators/spin_H.h

diff --git a/CMakeLists.txt b/CMakeLists.txt
index aa8cd011..3133a7b6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -22,9 +22,6 @@ option(Document "Build HTML document" OFF)
 
 option(USE_SYSTEM_BOOST "use Boost installed in system" OFF)
 
-option(BUILD_NEW_GENERATORS "build new file-generators" ON)
-option(BUILD_OLD_GENERATORS "build old file-generators" OFF)
-
 enable_language(C CXX)
 
 set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
@@ -69,16 +66,14 @@ include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src/third-party/plog)
 set(python_version_required 3.6)
 
 if(NOT PYTHON_EXECUTABLE)
-  if(BUILD_NEW_GENERATORS OR Testing OR Document)
-    if(${CMAKE_VERSION} VERSION_LESS 3.12)
-      find_package(PythonInterp ${python_version_required} REQUIRED)
-      set(python_version_mm "${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}")
-    else(${CMAKE_VERSION} VERSION_LESS 3.12)
-      find_package(Python3 ${python_version_required} COMPONENTS Interpreter REQUIRED)
-      set(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
-      set(python_version_mm "${Python3_VERSION_MAJOR}.${Python3_VERSION_MINOR}")
-    endif(${CMAKE_VERSION} VERSION_LESS 3.12)
-  endif()
+  if(${CMAKE_VERSION} VERSION_LESS 3.12)
+    find_package(PythonInterp ${python_version_required} REQUIRED)
+    set(python_version_mm "${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}")
+  else(${CMAKE_VERSION} VERSION_LESS 3.12)
+    find_package(Python3 ${python_version_required} COMPONENTS Interpreter REQUIRED)
+    set(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
+    set(python_version_mm "${Python3_VERSION_MAJOR}.${Python3_VERSION_MINOR}")
+  endif(${CMAKE_VERSION} VERSION_LESS 3.12)
 else()
   # check python version
   execute_process(COMMAND "${PYTHON_EXECUTABLE}" "-V" OUTPUT_VARIABLE result ERROR_VARIABLE result)
@@ -105,9 +100,7 @@ if(MPI_FOUND)
   add_subdirectory(src/pmwa)
 endif(MPI_FOUND)
 
-if(BUILD_NEW_GENERATORS)
-  add_subdirectory(tool)
-endif()
+add_subdirectory(tool)
 
 if (Testing)
   enable_testing()
diff --git a/src/dla/CMakeLists.txt b/src/dla/CMakeLists.txt
index b19bb003..dd2928dd 100644
--- a/src/dla/CMakeLists.txt
+++ b/src/dla/CMakeLists.txt
@@ -21,8 +21,3 @@ if(MPI_FOUND)
 endif(MPI_FOUND)
 
 set(dla_deps dla)
-
-if(BUILD_OLD_GENERATORS)
-add_subdirectory(generators)
-set(dla_deps ${dla_deps} dla_alg_old lattgene_C lattgene_T hamgen_H hamgen_B cfgene sfgene)
-endif(BUILD_OLD_GENERATORS)
diff --git a/src/dla/generators/CMakeLists.txt b/src/dla/generators/CMakeLists.txt
deleted file mode 100644
index dbece9b0..00000000
--- a/src/dla/generators/CMakeLists.txt
+++ /dev/null
@@ -1,26 +0,0 @@
-# include guard
-if(${CMAKE_PROJECT_NAME} STREQUAL "Project")
-  message(FATAL_ERROR "cmake should be executed not for 'src' subdirectory, but for the top directory of DLA.")
-endif(${CMAKE_PROJECT_NAME} STREQUAL "Project")
-
-find_package(LAPACK REQUIRED)
-
-include_directories(.)
-
-set(dla_alg_old_SOURCES dla_alg.cc)
-set(cfgene_SOURCES cfgene.cc)
-set(sfgene_SOURCES sfgene.cc)
-set(lattgene_C_SOURCES lattgene_C.cc)
-set(lattgene_T_SOURCES lattgene_T.cc)
-set(hamgen_H_SOURCES exact_H.cc)
-set(hamgen_B_SOURCES exact_B.cc)
-
-set(exelist dla_alg_old cfgene sfgene lattgene_C lattgene_T hamgen_H hamgen_B)
-
-foreach(exe IN LISTS exelist)
-  add_executable(${exe} ${${exe}_SOURCES})
-  install(TARGETS ${exe} RUNTIME DESTINATION bin)
-endforeach(exe)
-
-target_link_libraries(hamgen_H ${LAPACK_LIBRARIES})
-target_link_libraries(hamgen_B ${LAPACK_LIBRARIES})
diff --git a/src/dla/generators/boson_B.h b/src/dla/generators/boson_B.h
deleted file mode 100644
index cef44b87..00000000
--- a/src/dla/generators/boson_B.h
+++ /dev/null
@@ -1,127 +0,0 @@
-
-//============================================================================
-//    Spin Matrices
-//============================================================================
-
-class BosonOperator {
- public:
-  int K;
-  int D;
-  cmatrix I;
-  cmatrix UP;
-  cmatrix DN;
-  cmatrix X;
-  cmatrix Y;
-  cmatrix Z;
-  BosonOperator(int);
-  ~BosonOperator(){};
-};
-
-//----------------------------------------------------------------------------
-
-BosonOperator::BosonOperator(int K0) {
-  K = K0;  // K = 2 * S
-  D = K0 + 1;
-  I.resize(D, D);
-  UP.resize(D, D);
-  DN.resize(D, D);
-  X.resize(D, D);
-  Y.resize(D, D);
-  Z.resize(D, D);
-  I.zero();
-  UP.zero();
-  DN.zero();
-  X.zero();
-  Y.zero();
-  Z.zero();
-  for (int i = 0; i < D - 1; i++) {
-    UP.re(i + 1, i) = sqrt((double)(i + 1));  // sqrt((double)((i+1)*(K-i)));
-  }
-  DN = t(UP);
-  //  X = 0.5 * (UP + DN);
-  cmatrix temp = UP + DN;
-  X = 0.5 * temp;
-  Y = ((-0.5) * IUNIT) * (UP - DN);
-  for (int i = 0; i < D; i++) {
-    I.re(i, i) = 1.0;
-    Z.re(i, i) = (double)i;
-  }
-}
-
-//----------------------------------------------------------------------------
-
-class BosonOperatorSet {
- public:
-  int DS;     // the dimension of the 1-spin Hilbert space
-  int N;      // "N" in SU(N)
-  int K;      // the number of bosons on each site
-  int NSITE;  // the number of spins
-  int DIM;    // the dimension of the whole Hilbert space
-
-  cmatrix* UP;
-  cmatrix* DN;
-  cmatrix* X;
-  cmatrix* Y;
-  cmatrix* Z;
-  cmatrix I;
-
-  BosonOperatorSet(int K0, int NSITE0) {
-    printf("BosonSet> start.\n");
-
-    K = K0;
-    NSITE = NSITE0;
-
-    BosonOperator S(K);
-
-    DS = S.D;
-    UP = new cmatrix[NSITE];
-    DN = new cmatrix[NSITE];
-    X = new cmatrix[NSITE];
-    Y = new cmatrix[NSITE];
-    Z = new cmatrix[NSITE];
-
-    printf("  definig spins ...\n");
-    DIM = 1;
-
-    for (int i = 0; i < NSITE; i++) {
-      DIM *= (DS);
-      cmatrix up;
-      cmatrix dn;
-      cmatrix x;
-      cmatrix y;
-      cmatrix z;
-      up.unity();
-      dn.unity();
-      x.unity();
-      y.unity();
-      z.unity();
-      for (int j = 0; j < NSITE; j++) {
-        if (i == j) {
-          up = (S.UP) ^ up;
-          dn = (S.DN) ^ dn;
-          x = (S.X) ^ x;
-          y = (S.Y) ^ y;
-          z = (S.Z) ^ z;
-        } else {
-          up = (S.I) ^ up;
-          dn = (S.I) ^ dn;
-          x = (S.I) ^ x;
-          y = (S.I) ^ y;
-          z = (S.I) ^ z;
-        }
-
-        UP[i] = up;
-        DN[i] = dn;
-        X[i] = x;
-        Y[i] = y;
-        Z[i] = z;
-      }
-      printf("  ... done.\n");
-
-      I.resize(DIM);
-      I.identity();
-
-      printf("SpinSet> end.\n");
-    };
-  };
-};
diff --git a/src/dla/generators/canonical.h b/src/dla/generators/canonical.h
deleted file mode 100644
index 4f8e623c..00000000
--- a/src/dla/generators/canonical.h
+++ /dev/null
@@ -1,78 +0,0 @@
-
-//============================================================================
-//    Computation of Canonical Averages
-//============================================================================
-
-dgematrix DensityMatrix(double B, vector<double>& V, dgematrix& U) {
-  int DIM = V.size();
-  dgematrix UT(DIM, DIM);
-  dgematrix W(DIM, DIM);
-  UT = t(U);
-  W.zero();
-  double emin = V[0];
-  for (int i = 0; i < DIM; i++)
-    if (emin > V[i]) emin = V[i];
-  for (int i = 0; i < DIM; i++) W(i, i) = exp(-B * (V[i] - emin));
-  return U * W * UT;
-}
-
-//----------------------------------------------------------------------------
-
-double chi(double B, double E1, double E0) {
-  if (fabs(E1 - E0) < 1.0e-12) {
-    return B * exp(-B * E0);
-  } else {
-    return -(exp(-B * E1) - exp(-B * E0)) / (E1 - E0);
-  }
-}
-
-//----------------------------------------------------------------------------
-
-//
-// ans = (\beta)^{-1} \int_0^{\beta} dt ( <Q(t)Q(0)> - <Q(t)Q(0)> )
-//     --> <Q^2>-<Q>^2 (in the classical or the high-T limit)
-//
-double Susceptibility(double B, vector<double>& V, dgematrix& U, dgematrix& Q) {
-  int DIM = V.size();
-  double* E = new double[DIM];
-  double emin = V[0];
-  for (int i = 0; i < DIM; i++)
-    if (V[i] < emin) emin = V[i];
-  for (int i = 0; i < DIM; i++) E[i] = V[i] - emin;
-
-  dgematrix UT(DIM, DIM);
-  dgematrix W(DIM, DIM);
-  UT = t(U);
-  W = UT * Q * U;
-
-  double Z0 = 0.0;
-  double Z1 = 0.0;
-  double ZX = 0.0;
-  for (int i = 0; i < DIM; i++) {
-    Z0 += exp(-B * E[i]);
-    Z1 += exp(-B * E[i]) * W(i, i);
-    for (int j = 0; j < DIM; j++) {
-      ZX += chi(B, E[j], E[i]) * W(j, i) * W(i, j);
-    }
-  }
-
-  double ave = Z1 / Z0;
-  double ans = ZX / Z0 / B - ave * ave;
-  delete[] E;
-  return ans;
-}
-
-//----------------------------------------------------------------------------
-
-double CanonicalAverage(dgematrix& R, dgematrix& Q) {
-  double z0 = 0.0;
-  double z1 = 0.0;
-  int n = R.n;
-  dgematrix W(n, n);
-  W = R * Q;
-  for (int i = 0; i < n; i++) {
-    z0 += R(i, i);
-    z1 += W(i, i);
-  }
-  return z1 / z0;
-}
diff --git a/src/dla/generators/cfgene.cc b/src/dla/generators/cfgene.cc
deleted file mode 100644
index 85e7f6a2..00000000
--- a/src/dla/generators/cfgene.cc
+++ /dev/null
@@ -1,140 +0,0 @@
-/*---------------------------------------------
-
-   Generating cf.xml for a hypercubic lattice
-
-----------------------------------------------*/
-
-#include <cmath>
-#include <cstdlib>
-#include <cstring>
-#include <fstream>
-#include <iostream>
-#include <string>
-
-using namespace std;
-
-void ShowUsage(std::string const& exename) {
-  cout << "usage:\n";
-  cout << "    " << exename << " [-o filename] D L_1 ... L_D Ntau\n";
-  cout << "arguments:\n";
-  cout << "    D           ... dimension of lattice\n";
-  cout << "    L_1 ... L_D ... the liner size of the lattice. \n";
-  cout << "                    (must be even number. )\n";
-  cout << "    Ntau        ... the number of discretized imaginary time\n";
-  cout << "options:\n";
-  cout << "    -o filename ... output file (default: cf.xml)";
-}
-
-void WriteXML(int D, int L[], int Ntau, std::string const& filename) {
-  ofstream fout(filename.c_str());
-  fout.precision(15);
-  int N = 1;  // number of sites.
-  for (int i = 0; i < D; i++) {
-    N *= L[i];
-  }
-
-  int NumberOfElements = N * N;
-  int NumberOfKinds = N;
-
-  fout << "<CorrelationFunction>" << endl << endl;
-  fout << "<Comment>" << endl;
-  fout << "  " << D << "-dimension hypercubic lattice" << endl;
-  fout << "</Comment>" << endl << endl;
-
-  fout << "<Ntau>             " << Ntau << " </Ntau>" << endl;
-  fout << "<NumberOfElements> " << NumberOfElements << " </NumberOfElements>"
-       << endl;
-  fout << "<NumberOfKinds>    " << NumberOfKinds << " </NumberOfKinds>" << endl;
-  fout << endl;
-
-  fout << "<!-- <CF> [kind] [isite] [jsite] </CF> -->" << endl << endl;
-
-  int NB = 0;  // 3 * N ;   // number of bonds
-  int* x = new int[D];
-  int* dx = new int[D];
-
-  int kind = 0;
-
-  for (int di = 0; di < N; di++) {
-    int dk = di;
-    for (int q = 0; q < D; q++) {
-      dx[q] = dk % L[q];
-      dk /= L[q];
-    }
-
-    for (int i = 0; i < N; i++) {
-      int k = i;
-      for (int q = 0; q < D; q++) {
-        x[q] = k % L[q];
-        k /= L[q];
-      }
-      for (int q = 0; q < D; q++) {
-        x[q] = (x[q] + dx[q]) % L[q];
-      }
-      int j = 0;
-      for (int q = D - 1; q >= 0; q--) {
-        j *= L[q];
-        j += x[q];
-      }
-
-      fout << "<CF> " << kind << " " << i << " " << j << " </CF>" << endl;
-    }
-    kind++;
-  }
-  cout << "NumberOfKinds = " << kind << endl;
-  fout << endl;
-  fout << "</CorrelationFunction>" << endl;
-
-  delete[] x;
-}
-
-int main(int argc, char** argv) {
-  std::string exename(argv[0]);
-  std::string filename("cf.xml");
-  if (argc < 3) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  if (std::strcmp(argv[1], "-o") == 0) {
-    filename = argv[2];
-    argc -= 2;
-    argv += 2;
-  }
-  int NARG = 3;
-  if (argc < NARG + 1) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  int iarg = 1;
-  const int D = atoi(argv[iarg]);
-  iarg++;
-  if (argc != D + 3) {
-    ShowUsage(exename);
-    exit(0);
-  }
-
-  int* L = new int[D];
-
-  for (int i = 0; i < D; i++) {
-    L[i] = atoi(argv[iarg]);
-    iarg++;
-  }
-  int Ntau = atoi(argv[iarg]);
-
-  int EvenOrOdd = 0;
-  cout.precision(15);
-  cout << "D     = " << D << endl;
-  for (int i = 0; i < D; i++) {
-    cout << "L" << i << "    = " << L[i] << endl;
-    EvenOrOdd += L[i] % 2;
-  }
-
-  if (EvenOrOdd) {
-    cout << "Warnig: L should be an even number." << endl;
-  }
-
-  WriteXML(D, L, Ntau, filename);
-  cout << "... done." << endl;
-  delete[] L;
-  return 0;
-}
diff --git a/src/dla/generators/dla_alg.cc b/src/dla/generators/dla_alg.cc
deleted file mode 100644
index 198797d8..00000000
--- a/src/dla/generators/dla_alg.cc
+++ /dev/null
@@ -1,1216 +0,0 @@
-#include "dla_alg.h"
-
-bool isdiagonal(int* x, int NBODY) {
-  for (int i = 0; i < NBODY; ++i) {
-    if (x[2 * i] != x[2 * i + 1]) {
-      return false;
-    }
-  }
-  return true;
-}
-
-int main(int argc, char** argv) {
-  HFILE = "hamiltonian.xml";
-  AFILE = "algorithm.xml";
-  if (argc > 1) {
-    HFILE = argv[1];
-  }
-  if (argc > 2) {
-    AFILE = argv[2];
-  }
-  printf("HFILE= %s\n", HFILE);
-  printf("AFILE= %s\n", AFILE);
-
-  XML::Block X(HFILE, "Hamiltonian");
-  FALG = fopen(AFILE, "w");
-  // FALG = stdout;
-
-  XML::Block& XGEN = X["General"];
-  GENERAL G(XGEN);
-  Site = new SITE[NSTYPE];
-  Source = new SOURCE[NSTYPE];
-  Interaction = new INTERACTION[NITYPE];
-  Vertex = new VERTEX[NVTYPE];
-
-  for (int i = 0; i < X.NumberOfBlocks(); i++) {
-    XML::Block& B = X[i];
-    const std::string& name = B.getName();
-    if (name == "Site") {
-      int id = B["STYPE"].getInteger();
-      Site[id].load(B);
-      //      Site[id].dump();
-    } else if (name == "Source") {
-      int id = B["TTYPE"].getInteger();
-      Source[id].load(B);
-      //      Source[id].dump();
-      //      Vertex[Source[id].VTYPE].dump();
-    } else if (name == "Interaction") {
-      int id = B["ITYPE"].getInteger();
-      Interaction[id].load(B);
-      //      Interaction[id].dump();
-      //      Vertex[Interaction[id].VTYPE].dump();
-    }
-  }
-
-  // Scattering probs
-
-  for (int i = 0; i < NVTYPE; i++) {
-    Vertex[i].Grouping();
-  }
-  // Worm vertex has the shared EBASE
-  G.WeightDiagonal = -INF;
-
-  for (int i = 0; i < NSTYPE; i++) {
-    VERTEX& V = Site[i].V();
-    double eb = V.ComputeEBASE();
-    if (G.WeightDiagonal < 0.5 * eb) G.WeightDiagonal = 0.5 * eb;
-  }
-
-  for (int i = 0; i < NSTYPE; i++) {
-    VERTEX& V = Site[i].V();
-    V.EBASE = 2.0 * G.WeightDiagonal;
-  }
-
-  // Interaction vertex has the own EBASE
-
-  for (int i = 0; i < NITYPE; i++) {
-    VERTEX& V = Interaction[i].V();
-    V.EBASE = V.ComputeEBASE();
-  }
-
-  for (int i = 0; i < NVTYPE; i++) Vertex[i].ComputeScatteringProbability();
-
-  for (int i = 0; i < NSTYPE; i++) {
-    Site[i].SetInitialHeadTypeProbability();
-  }
-
-  for (int i = 0; i < NITYPE; i++) {
-    Interaction[i].SetVertexDensity();
-  }
-
-  fprintf(FALG, "<Algorithm>\n");
-  G.write();
-  for (int i = 0; i < NSTYPE; i++) Site[i].write();
-  for (int i = 0; i < NITYPE; i++) Interaction[i].write();
-  for (int i = 0; i < NVTYPE; i++) Vertex[i].write();
-  fprintf(FALG, "\n</Algorithm>\n");
-
-  delete[] Site;
-  delete[] Source;
-  delete[] Interaction;
-  delete[] Vertex;
-}
-
-//######################################################################
-
-GENERAL::GENERAL(XML::Block& X) {
-  NHTYPE = 0;
-  NXMAX = 0;
-
-  comment = X["Comment"].getJoinedString();
-  NSTYPE = X["NSTYPE"].getInteger();
-  NITYPE = X["NITYPE"].getInteger();
-  NXMAX = X["NXMAX"].getInteger();
-  NVTYPE = NSTYPE + NITYPE;
-}
-
-//======================================================================
-
-void GENERAL::write() {
-  fprintf(FALG, "\n");
-  fprintf(FALG, "  <Comment>\n");
-  fprintf(FALG, "    %s\n", comment.c_str());
-  fprintf(FALG, "  </Comment>\n");
-  fprintf(FALG, "\n");
-  fprintf(FALG, "  <General>\n");
-  fprintf(FALG, "    <NSTYPE> %2d </NSTYPE>\n", NSTYPE);
-  fprintf(FALG, "    <NITYPE> %2d </NITYPE>\n", NITYPE);
-  fprintf(FALG, "    <NVTYPE> %2d </NVTYPE>\n", NVTYPE);
-  fprintf(FALG, "    <NXMAX>  %2d </NXMAX>\n", NXMAX);
-  fprintf(FALG, "    <WDIAG>  %24.16lf </WDIAG>\n", WeightDiagonal);
-  fprintf(FALG, "  </General>\n");
-}
-
-//######################################################################
-
-void SITE::load(XML::Block& X) {
-  ID = X["STYPE"].getInteger();
-  TTYPE = X["TTYPE"].getInteger();
-  NX = X["NX"].getInteger();
-  _T = &Source[TTYPE];
-}
-
-//======================================================================
-
-void SITE::SetInitialHeadTypeProbability() {
-  VTYPE = Source[TTYPE].VTYPE;
-  // WormCreationNewState( xi, ch)
-  // "xi". .. the initial local state
-  // "ch" ... the scattering channel id
-  WormCreationNewState.init(2, NXMAX, 2 * NXMAX);
-  WormCreationNewState.set_all(STATE::UNDEF);
-  WormCreationDirection.init(2, NXMAX, 2 * NXMAX);
-  WormCreationDirection.set_all(DIR::UNDEF);
-  WormCreationProbability.init(2, NXMAX, 2 * NXMAX);
-  WormCreationProbability.set_all(0.0);
-  NumberOfChannels = new int[NX];
-  VERTEX& V = Vertex[VTYPE];
-  for (int i = 0; i < V.NICG; i++) {
-    InitialConfigurationGroup& icg = V.ICG(i);
-    for (int j = 0; j < icg.NIC; j++) {
-      InitialConfiguration& ic = icg.IC[j];
-      if (ic.INC == DIR::UNDEF) {
-        int x0 = ic.State[0];
-        int NCH = ic.NCH;
-        NumberOfChannels[x0] = NCH;
-        for (int c = 0; c < NCH; c++) {
-          int st = ic.FinalState[c];
-          int out = ic.FinalDirection[c];
-          double p = ic.ScatteringProbability[c];
-          WormCreationNewState(x0, c) = st;
-          WormCreationDirection(x0, c) = out;
-          WormCreationProbability(x0, c) = p;
-        }
-      }
-    }
-  }
-}
-
-//======================================================================
-
-void SITE::write() {
-  fprintf(FALG, "\n");
-  fprintf(FALG, "  <Site>\n");
-  fprintf(FALG, "    <STYPE> %d </STYPE>\n", ID);
-  fprintf(FALG, "    <NumberOfStates> %d </NumberOfStates>\n", NX);
-  fprintf(FALG, "    <VertexTypeOfSource> %d </VertexTypeOfSource>\n", VTYPE);
-  // printf("    <EBASE> %24.16f </EBASE>\n", V().EBASE );
-  // fprintf(FALG,"    <EBASE> %24.16f </EBASE>\n", V().EBASE );
-  for (int x0 = 0; x0 < NX; x0++) {
-    fprintf(FALG, "    <InitialConfiguration>\n");
-    fprintf(FALG, "      <State> %d </State>\n", x0);
-    int NCH = NumberOfChannels[x0];
-    fprintf(FALG, "      <NumberOfChannels> %d </NumberOfChannels>\n", NCH);
-    for (int c = 0; c < NCH; c++) {
-      int diri = WormCreationDirection(x0, c);
-      int xnew = WormCreationNewState(x0, c);
-      double p = WormCreationProbability(x0, c);
-      if (p == 0.0) continue;
-      fprintf(FALG, "      <Channel> %4d %4d %24.16lf </Channel>\n", diri, xnew,
-              p);
-    }
-    fprintf(FALG, "    </InitialConfiguration>\n");
-  }
-  fprintf(FALG, "  </Site>\n");
-}
-
-//======================================================================
-
-void SITE::dump() {
-  printf("\n");
-  printf("SITE[%d]> tt=%2d, NX=%2d\n", ID, TTYPE, NX);
-  if (NumberOfChannels == 0) {
-    printf(" ... Channels have not been defined yet.\n");
-  } else {
-    for (int x = 0; x < NX; x++) {
-      printf("  x=%2d : ", x);
-      int NCH = NumberOfChannels[x];
-      printf("  NCH= %d\n", NCH);
-      for (int c = 0; c < NCH; c++) {
-        printf("   %2d %2d %8.3lf\n", WormCreationDirection(x, c),
-               WormCreationNewState(x, c), WormCreationProbability(x, c));
-      }
-    }
-  }
-}
-
-//######################################################################
-
-void SOURCE::load(XML::Block& X) {
-  ID = X["TTYPE"].getInteger();
-  VTYPE = ID;
-  VERTEX& V = Vertex[VTYPE];
-  _V = &V;
-  V.ID = VTYPE;
-  V.NBODY = 1;
-  V.CATEGORY = VCAT::WORM;
-  V.load(X);
-}
-
-//======================================================================
-
-void SOURCE::dump() { printf("SOURCE[%d]> vt=%2d\n", ID, VTYPE); }
-
-//######################################################################
-
-void INTERACTION::load(XML::Block& X) {
-  ID = X["ITYPE"].getInteger();
-  VTYPE = ID + NSTYPE;
-  _V = &(Vertex[VTYPE]);
-  NBODY = X["NBODY"].getInteger();
-
-  int NLEG = 2 * NBODY;
-
-  Sign.init("Weight", NLEG, NXMAX, ARRAY::EOL);
-  Sign.set_all(0.0);
-
-  int* x = new int[NLEG];
-  for (int i = 0; i < X.NumberOfBlocks(); i++) {
-    XML::Block& B = X[i];
-    const std::string& name = B.getName();
-    if (name == "Weight") {
-      for (int i = 0; i < NLEG; i++) x[i] = B.getInteger(i);
-      double w = B.getDouble(NLEG);
-      if (w < 0.0 && !isdiagonal(x, NBODY)) {
-        Sign(x) = -1.0;
-      }
-    }
-  }
-  delete[] x;
-
-  V().ID = VTYPE;
-  V().NBODY = NBODY;
-  V().CATEGORY = VCAT::INT;
-  V().load(X);
-}
-
-//======================================================================
-
-void INTERACTION::SetVertexDensity() {
-  VERTEX& V = Vertex[VTYPE];
-  IndexSystem& I = V.Weight.index_system();
-  int* x = new int[I.dimension()];
-  VertexDensity.init(NBODY, NXMAX, ARRAY::EOL);
-  IndexSystem& J = VertexDensity.index_system();
-  int* y = new int[J.dimension()];  // edit sakakura
-  VertexDensity.set_all(0.0);
-  for (int g = 0; g < V.NICG; g++) {
-    InitialConfigurationGroup& G = V.ICG(g);
-    for (int c = 0; c < G.NIC; c++) {
-      InitialConfiguration& C = G.IC[c];
-      if (!C.isKink()) {
-        int sti = C.STI;
-        double w = C.vertex_weight();
-        I.coord(sti, x);
-        for (int i = 0; i < NBODY; i++) y[i] = x[2 * i];
-        VertexDensity(y) = w;
-      }
-    }
-  }
-  delete[] x;
-  delete[] y;
-}
-
-//======================================================================
-
-void INTERACTION::write() {
-  fprintf(FALG, "\n");
-  fprintf(FALG, "  <Interaction>\n");
-  fprintf(FALG, "    <ITYPE> %d </ITYPE>\n", ID);
-  fprintf(FALG, "    <VTYPE> %d </VTYPE>\n", VTYPE);
-  fprintf(FALG, "    <NBODY> %d </NBODY>\n", NBODY);
-  fprintf(FALG, "    <EBASE> %24.16lf </EBASE>\n", V().EBASE);
-  if (VertexDensity.isDefined()) {
-    IndexSystem& I = VertexDensity.index_system();
-    int* x = new int[NBODY];
-    for (int i = 0; i < I.size(); i++) {
-      double d = VertexDensity[i];
-      if (d > 0.0) {
-        I.coord(i, x);
-        fprintf(FALG, "    <VertexDensity> ");
-        for (int j = 0; j < NBODY; j++) fprintf(FALG, " %2d", x[j]);
-        fprintf(FALG, " %24.16lf </VertexDensity>\n", d);
-      }
-    }
-    delete[] x;
-  } else {
-    //  fprintf(FALG,"  VertexDensity is not defined.\n");
-  }
-
-  int NLEG = 2 * NBODY;
-  Array<double>& Weight = V().Weight;
-  IndexSystem& I = Weight.index_system();
-  int* x = new int[2 * NBODY];
-  for (int i = 0; i < I.size(); ++i) {
-    if (Weight[i] != 0.0) {
-      I.coord(i, x);
-      if (!isdiagonal(x, NBODY)) {
-        double sgn = Weight[i] > 0.0 ? 1.0 : -1.0;
-        fprintf(FALG, "    <Sign> ");
-        for (int j = 0; j < 2 * NBODY; ++j) {
-          fprintf(FALG, " %2d", x[j]);
-        }
-        fprintf(FALG, " %24.16lf </Sign>\n", sgn);
-      }
-    }
-  }
-  delete[] x;
-
-  fprintf(FALG, "  </Interaction>\n");
-}
-
-//======================================================================
-
-void INTERACTION::dump() {
-  printf("INTERACTION[%d] ", ID);
-  printf(", vt=%2d", VTYPE);
-  printf(", nb=%2d", NBODY);
-  printf("\n");
-  if (VertexDensity.isDefined()) {
-    IndexSystem& I = VertexDensity.index_system();
-    int* x = new int[NBODY];
-    for (int i = 0; i < I.size(); i++) {
-      if (VertexDensity[i] > 0.0) {
-        I.coord(i, x);
-        printf(" x= (");
-        for (int j = 0; j < I.dimension(); j++) printf(" %1d", x[j]);
-        printf(")");
-        printf(" density= %8.3f\n", VertexDensity[i]);
-        printf("\n");
-      }
-    }
-    delete[] x;
-  } else {
-    printf("  VertexDensity is not defined.\n");
-  }
-}
-
-//######################################################################
-
-void VERTEX::load(XML::Block& X) {
-  NLEG = 2 * NBODY;
-  Weight.init("Weight", NLEG, NXMAX, ARRAY::EOL);
-  setINDX(Weight.index_system());
-  Weight.set_all(0.0);
-
-  NST = Weight.size();
-  SiteTypeOfLeg.init("STYPE", 1, NLEG);
-  SiteTypeOfLeg.set_all(STYPE::UNDEF);
-  for (int i = 0; i < NBODY; i++) {
-    int st = X["STYPE"].getInteger(i);
-    SiteTypeOfLeg[2 * i] = st;
-    SiteTypeOfLeg[2 * i + 1] = st;
-  }
-  int* x = new int[NLEG];
-  for (int i = 0; i < X.NumberOfBlocks(); i++) {
-    XML::Block& B = X[i];
-    const std::string& name = B.getName();
-    if (name == "Weight") {
-      for (int i = 0; i < NLEG; i++) x[i] = B.getInteger(i);
-      double w = B.getDouble(NLEG);
-      if (w < 0.0 && !isdiagonal(x, NBODY)) {
-        w *= -1.0;
-      }
-      Weight(x) = w;
-    }
-  }
-  EBASE = 0.0;
-  delete[] x;
-}
-
-//======================================================================
-
-int VERTEX::NumberOfValidInitialConfigurations() {
-  int c = 0;
-  for (int i = 0; i < NICG; i++) {
-    int n = ICG(i).NumberOfValidInitialConfigurations();
-    c += n;
-  }
-  return c;
-}
-
-//======================================================================
-
-void VERTEX::Grouping() {
-  // printf("VERTEX::Grouping> Start ID=%d\n", ID);
-  // printf("VERTEX::Grouping> Start.\n");
-  // printf("VERTEX::Grouping> Pass 1\n");
-
-  int NICmax = 2 * NBODY * NXMAX + 1;
-  checked.init(3, NST, NLEG, NXMAX);
-  checked.set_all(0);
-  int* x = new int[NLEG];
-  //  int x[NLEG];
-  NICG = 0;
-  if (isWorm()) {
-    for (int sti = 0; sti < NST; sti++) {
-      if (testForbidden(sti)) continue;
-      if (testKink(sti)) continue;
-      int inc = DIR::UNDEF;
-      int xinc = STATE::UNDEF;
-      InitialConfigurationGroup& icg = *(new InitialConfigurationGroup);
-      icg.init(*this, NBODY, sti, inc, xinc);
-      add_ICG(icg);
-      // printf("VERTEX::Grouping> New initial configuration group. ID= %d\n",
-      // NICG);
-      NICG++;
-    }
-  } else {
-    for (int sti = 0; sti < NST; sti++) {
-      INDX().coord(sti, x);
-      if (testForbidden(sti)) continue;
-      for (int inc = 0; inc < NLEG; inc++) {
-        int xinc_old = x[inc];
-        int st = SiteTypeOfLeg[inc];
-        SOURCE& W = Source[Site[st].TTYPE];
-        for (int xinc = 0; xinc < NXMAX; xinc++) {
-          // W.V().Weight.index_system().dump();
-          double ww = W.Weight(xinc_old, xinc);
-          if (ww == 0.0) continue;
-          if (checked(sti, inc, xinc) == 1) continue;
-
-          // printf("VERTEX::Grouping> New initial configuration group. ID=
-          // %d\n", NICG);
-          InitialConfigurationGroup& icg = *(new InitialConfigurationGroup);
-          icg.init(*this, NBODY, sti, inc, xinc);
-          add_ICG(icg);
-          NICG++;
-        }
-      }
-    }
-    // printf(" NICG= %d, icg.size()= %d\n", NICG, icg.size());
-  }
-  delete[] x;
-}
-
-//======================================================================
-
-double VERTEX::ComputeEBASE() {
-  //
-  // Calculate EBASE as the minimum additional term to
-  // make probability of offdiagonal bounce of negative Hamiltonian (-H)
-  // be zero.
-  //
-  // - For worm creation/annihilation, avoid direct transition
-  //   between offdiagonal states.
-  //   In other words, let EBASE be so large that
-  //   the probability of a worm annihilation is 1.
-  // - For interaction vertex, minimize diagonal weights
-  //   as long as no offdiagonal bounce occurs.
-  //
-
-  for (int icg = 0; icg < NICG; icg++) ICG(icg).ResetWeight();
-
-  double eb;
-  if (CATEGORY == VCAT::TERM) {
-    eb = 1.0;
-  } else if (CATEGORY == VCAT::WORM) {
-    eb = -INF;
-    for (int icg = 0; icg < NICG; icg++) {
-      double e = ICG(icg).sum_of_all_weights();
-      if (eb < e) eb = e;
-    }
-  } else if (CATEGORY == VCAT::INT) {
-    eb = -INF;
-    for (int icg = 0; icg < NICG; icg++) {
-      double e = ICG(icg).ebase();
-      if (eb < e) eb = e;
-    }
-  } else {
-    printf("VERTEX::ComputeEBASE> Error.\n");
-    exit(0);
-  }
-  return eb;
-}
-
-//======================================================================
-
-void VERTEX::ComputeScatteringProbability() {
-  // printf("\nVERTEX::ComputeScatteringProbability> Start.\n");
-  NICV = 0;
-  for (int icg = 0; icg < NICG; icg++) {
-    ICG(icg).ResetWeight();
-    ICG(icg).numbering(NICV);
-    ICG(icg).ScatteringProbability();
-  }
-  // printf("\nVERTEX::ComputeScatteringProbability> End.\n");
-}
-
-//======================================================================
-
-void VERTEX::dump() {
-  printf("VERTEX[%d] cat= %d", ID, CATEGORY);
-  printf(" nb= %d", NBODY);
-  printf(" nl= %d", NLEG);
-  printf(" ncig= %d", NICG);
-  printf(" eb= %8.3f", EBASE);
-  printf(" st= (");
-  for (int i = 0; i < NLEG; i++) printf(" %1d", SiteTypeOfLeg(i));
-  printf(")\n");
-  int* x = new int[NLEG];  // edit sakakura
-                           //  int x[NLEG];
-  for (int i = 0; i < NST; i++) {
-    if ((!testKink(i)) || Weight[i] != 0.0) {
-      if (testKink(i)) {
-        printf(" kink:yes ");
-      } else {
-        printf(" kink:no  ");
-      }
-      if (Weight[i] == 0.0) {
-        printf(" w=0:yes ");
-      } else {
-        printf(" w=0:no  ");
-      }
-      INDX().coord(i, x);
-      printf(" (");
-      for (int j = 0; j < NLEG; j++) printf(" %1d", x[j]);
-      printf(") --> w= %8.3f\n", Weight[i]);
-    }
-  }
-  for (int i = 0; i < NICG; i++) {
-    printf("\n");
-    ICG(i).dump();
-  }
-  delete[] x;
-}
-
-//======================================================================
-
-void VERTEX::write() {
-  fprintf(FALG, "\n");
-  fprintf(FALG, "  <Vertex>\n");
-  fprintf(FALG, "    <VTYPE> %d </VTYPE>\n", ID);
-  fprintf(FALG, "    <VCATEGORY> %d </VCATEGORY>\n", CATEGORY);
-  fprintf(FALG, "    <NBODY> %d </NBODY>\n", NBODY);
-  int NICW = NumberOfValidInitialConfigurations();
-  fprintf(FALG,
-          "    <NumberOfInitialConfigurations> %d "
-          "</NumberOfInitialConfigurations>\n",
-          NICW);
-  for (int i = 0; i < NICG; i++) ICG(i).write();
-  fprintf(FALG, "\n  </Vertex>\n");
-}
-
-//======================================================================
-
-bool VERTEX::testKink(int ist) {
-  // printf("VERTEX::testKink> D = %d\n", D);
-  int* x = new int[NLEG];
-  INDX().coord(ist, x);
-  bool ans = false;
-  for (int i = 0; i < NLEG / 2; i++) {
-    if (x[2 * i] != x[2 * i + 1]) ans = true;
-  }
-  delete[] x;
-  return ans;
-}
-
-//======================================================================
-
-bool VERTEX::testForbidden(int ist) {
-  if (!testKink(ist)) return false;
-  if (Weight[ist] > 0.0) return false;
-  return true;
-}
-
-//######################################################################
-
-void InitialConfigurationGroup::init(VERTEX& v, int nbody, int sti, int inc,
-                                     int xinc) {
-  // handling a group of initial configurations
-  // that scatters into each other
-
-  // printf("InitialConfigurationGroup::init> Start. icg=%d\n", ID);
-
-  /*
-printf("nbody= %d, nicmax= %d, sti= %d, inc= %d, xinc= %d\n",
-nbody, NICmax, sti, inc, xinc);
-  */
-
-  _V = &v;
-  NBODY = nbody;
-  NLEG = 2 * nbody;
-
-  // int x[NLEG];
-  // int y[NLEG];
-  int* x = new int[NLEG];  // edit sakakura
-  int* y = new int[NLEG];
-  IndexSystem& INDX = V().INDX();
-  INDX.coord(sti, x);
-  if (inc != DIR::UNDEF) x[inc] = xinc;
-  // x[l] is the local state on the l-th leg
-  // just after the entrance of the incoming
-  // worm head (before its outgoing).
-
-  /*
-if (MON) {
-printf("InitialConfigurationGroup::init> ");
-printf("Grouping the configs related to ... \n");
-printf("  ...  nb= %d, sti= %d (", nbody, sti);
-for (int i=0; i<NLEG; i++) printf(" %1d", x[i]);
-printf("), hti= %d, inc= %d\n", hti, inc);
-}
-  */
-
-  int sitei = inc / 2;
-  int diri = inc % 2;
-
-  NIC = 0;
-  if (inc == DIR::UNDEF) NIC = 1;  // generate worm
-  for (int out = 0; out < NLEG; out++) {
-    for (int xout = 0; xout < NXMAX; xout++) {
-      // xout is the new local state on the outgoing leg
-      for (int d = 0; d < NLEG; d++) y[d] = x[d];
-      y[out] = xout;
-      int stf = INDX(y);
-      if (stf == STATE::UNDEF) continue;
-      if (V().testForbidden(stf)) continue;
-      NIC++;
-    }
-  }
-
-  U.init("U", 1, NIC);
-  U.set_all(0.0);
-  //  W.init("W",2,NIC,NIC);
-  //  W.set_all(0.0);
-  //  P.init("P",2,NIC,NIC);
-  //  P.set_all(0.0);
-  IC.init(1, NIC);
-
-  int ic = 0;
-  if (inc == DIR::UNDEF) {  // generate worm
-    IC[ic].init(V(), sti, inc, xinc, NIC);
-    ic++;
-  }
-  for (int out = 0; out < NLEG; out++) {
-    for (int xout = 0; xout < NXMAX; xout++) {
-      // xout is the new local state on the outgoing leg
-      for (int d = 0; d < NLEG; d++) y[d] = x[d];
-      y[out] = xout;
-      /*
-printf(" (");
-for (int i=0; i<INDX.dimension(); i++) printf(" %d", y[i]);
-printf(")\n");
-      */
-      int stf = INDX(y);
-      if (stf == STATE::UNDEF) continue;
-      if (V().testForbidden(stf)) continue;
-      IC[ic].init(V(), stf, out, x[out], NIC);
-      V().checked(stf, out, x[out]) = 1;
-      ic++;
-    }
-  }
-  delete[] x;
-  delete[] y;
-  // if (MON) printf("InitialConfigurationGroup::init> End.\n");
-}
-
-//======================================================================
-
-int InitialConfigurationGroup::NumberOfValidInitialConfigurations() {
-  int c = 0;
-  for (int i = 0; i < NIC; i++) {
-    if (IC[i].isValid()) c++;
-  }
-  return c;
-}
-
-//======================================================================
-
-void InitialConfigurationGroup::ResetWeight() {
-  for (int i = 0; i < NIC; i++) {
-    U[i] = IC[i].weight();
-    /*
-    for (int j=0; j<NIC; j++) {
-      W(i,j) = 0.0;
-      P(i,j) = 0.0;
-    }
-    */
-  }
-}
-
-//======================================================================
-
-int InitialConfigurationGroup::number_of_diagonal_states() {
-  int n = 0;
-  for (int i = 0; i < NIC; i++) {
-    if (!IC[i].isKink()) n++;
-  }
-  return n;
-}
-
-//======================================================================
-
-int InitialConfigurationGroup::number_of_offdiagonal_states() {
-  return NIC - number_of_diagonal_states();
-}
-
-//======================================================================
-
-double InitialConfigurationGroup::minimum_diagonal_weight() {
-  double dmin = INF;
-  for (int i = 0; i < NIC; i++) {
-    if (!IC[i].isKink()) {
-      if (U[i] < dmin) dmin = U[i];
-    }
-  }
-  return dmin;
-}
-
-//======================================================================
-
-double InitialConfigurationGroup::maximum_diagonal_weight() {
-  double dmax = -INF;
-  for (int i = 0; i < NIC; i++) {
-    if (!IC[i].isKink()) {
-      if (U[i] > dmax) dmax = U[i];
-    }
-  }
-  return dmax;
-}
-
-//======================================================================
-
-double InitialConfigurationGroup::maximum_offdiagonal_weight() {
-  double omax = -INF;
-  for (int i = 0; i < NIC; i++) {
-    if (IC[i].isKink()) {
-      if (U[i] > omax) omax = U[i];
-    }
-  }
-  return omax;
-}
-
-//======================================================================
-
-double InitialConfigurationGroup::sum_of_all_weights() {
-  double sum = 0.0;
-  for (int i = 0; i < NIC; i++) sum += U[i];
-  return sum;
-}
-
-//======================================================================
-
-double max(double& d0, double& d1, double& d2) {
-  double d = d0;
-  if (d1 > d) d = d1;
-  if (d2 > d) d = d2;
-  return d;
-}
-
-//======================================================================
-
-double InitialConfigurationGroup::ebase() {
-  // if ( MON ) printf("InitialConfigurationGroup::ebase> Start.\n");
-
-  double eb, eb0, eb1, eb2;
-
-  int nd = number_of_diagonal_states();
-  int no = NIC - nd;
-  if (nd == 0) {
-    eb = -INF;
-  } else {
-    double dmax = maximum_diagonal_weight();
-    double dmin = minimum_diagonal_weight();
-    double omax = maximum_offdiagonal_weight();
-    double sum = sum_of_all_weights();
-    // eb0 ... minimum EBASE to avoid negative diagonal elements
-    // eb1 ... minimum EBASE to avoid offdiagonal bounce
-    // eb2 ... minimum EBASE to avoid direct transition between offdiagonal
-    // states
-    eb0 = -dmin;
-    eb1 = (2.0 * omax - sum) / (double)nd;
-    //  eb2 = -INF;
-    //  eb2 = sum - 2.0 * dmax;
-    if (eb0 >= eb1) {
-      eb = eb0;
-    } else {
-      eb = eb1;
-    }
-
-    // Above we calculate EBASE from Hamilonian weights producted by worm
-    // weights
-    double ww = 0.0;
-    for (int i = 0; i < NIC; i++) {
-      if (!IC[i].isKink()) {
-        double wt = IC[i].worm_weight();
-        if (fabs(ww) > EPS && fabs(wt) > EPS && wt != ww) {
-          printf("InitialConfigurationGroup::ebase> Error.\n");
-          printf("  The worm weights of non-kinks are not equal.\n");
-          dump();
-          exit(0);
-        }
-        ww = IC[i].worm_weight();
-      }
-    }
-    eb /= ww;
-  }
-
-  // if ( MON ) printf("InitialConfigurationGroup::ebase> End.\n");
-
-  return eb;
-}
-
-//======================================================================
-
-void InitialConfigurationGroup::add_to_diagonal_element(double E) {
-  for (int i = 0; i < NIC; i++) {
-    if (!IC[i].isKink()) U[i] += E;
-  }
-}
-
-//======================================================================
-
-void InitialConfigurationGroup::numbering(int& ic) {
-  for (int i = 0; i < NIC; i++) {
-    if (U[i] > -EPS && U[i] < EPS) {
-      IC[i].id() = ICONF::UNDEF;
-    } else {
-      IC[i].id() = ic;
-      ic++;
-    }
-  }
-}
-
-//======================================================================
-
-void InitialConfigurationGroup::ScatteringProbability() {
-  // printf("\nInitialConfigurationGroup::ScatteringProbability> Start. ICG=
-  // %d\n", ID); if (MON) dump();
-
-  // NIC : # of initial state
-  // U   : weights of initial states
-  // UP  : weights of initial states (sorted)
-  Array<int> PERM("PERM");
-  PERM.init(1, NIC);
-  bubble_sort(NIC, U, PERM);
-
-  Array<double> UP("UP");
-  UP.init(1, NIC);
-  for (int i = 0; i < NIC; i++) {
-    UP[i] = U[PERM[i]];
-  }
-
-  Array<double> P("P");
-  Array<double> W("W");
-  Array<double> WP("WP");
-  P.init(2, NIC, NIC);
-  W.init(2, NIC, NIC);
-  WP.init(2, NIC, NIC);
-
-  SolveWeightEquation(NIC, UP, WP);
-
-  /*
-  for (int i=0; i<NIC; i++) {
-    for (int j=0; j<NIC; j++) W(PERM[i],PERM[j]) = WP(i,j);
-  }
-*/
-
-  IndexSystem& INDX = V().INDX();
-  for (int i = 0; i < NIC; i++) {
-    int I = PERM[i];
-    if (-EPS < U[I] && U[I] < EPS) continue;
-    int nc = 0;
-    for (int j = 0; j < NIC; j++) {
-      int J = PERM[j];
-      if (-EPS < U[J] && U[J] < EPS) continue;
-      W(I, J) = WP(i, j);
-      P(I, J) = W(I, J) / U[I];
-      if (P(I, J) > EPS) {
-        double p = P(I, J);
-        int stf = IC[J].STI;
-        int out = IC[J].INC;
-        int xout = STATE::UNDEF;
-        if (out != DIR::UNDEF) {
-          xout = INDX.coord(stf, out);
-        }
-        IC[I].FinalDirection[nc] = out;
-        IC[I].FinalState[nc] = xout;
-        IC[I].ScatteringProbability[nc] = p;
-        nc++;
-      }
-    }
-    IC[I].NCH = nc;
-  }
-
-  // if (MON) printf("  ==>\n");
-}
-
-//======================================================================
-
-void InitialConfigurationGroup::dump() {
-  for (int i = 0; i < NIC; i++) {
-    IC[i].dump();
-  }
-  for (int i = 0; i < NIC; i++) {
-    printf("(%d)", i);
-    printf(" ID= %2d", IC[i].id());
-    printf(" : U= %4.1f", U[i]);
-    /*
-    printf(" , W= ");
-    for (int j=0; j<NIC; j++) {
-      printf(" %4.1f", W(i,j));
-    }
-    printf(" , P= ");
-    for (int j=0; j<NIC; j++) {
-      printf(" %6.3f", P(i,j));
-    }
-    printf("\n");
-    */
-  }
-}
-
-//======================================================================
-
-void InitialConfigurationGroup::write() {
-  for (int i = 0; i < NIC; i++) {
-    IC[i].write();
-  }
-}
-
-//######################################################################
-
-void InitialConfiguration::init(VERTEX& v, int sti, int inc, int xinc,
-                                int NICmax) {
-  // define an initial configuration for which
-  // the initial vertex state is "sti"
-  // the incoming worm head is coming on the "inc"-th leg
-  // the new state on the leg after the passage of the worm head is "xinc"
-
-  setV(v);
-  STI = sti;
-  INC = inc;
-  XINC = xinc;
-  NLEG = V().NLEG;
-  NBODY = 2 * NLEG;
-  State = new int[NLEG];
-  ScatteringProbability = new double[NICmax];
-  FinalState = new int[NICmax];
-  FinalDirection = new int[NICmax];
-  V().INDX().coord(STI, State);
-  NCH = 0;
-}
-
-//======================================================================
-
-bool InitialConfiguration::isValid() {
-  if (INC == DIR::UNDEF) return false;
-  if (fabs(weight()) < EPS) return false;
-  return true;
-}
-
-//======================================================================
-
-bool InitialConfiguration::isKink() { return V().testKink(STI); }
-
-//======================================================================
-
-double InitialConfiguration::weight() {
-  double w = vertex_weight() * worm_weight();
-  return w;
-}
-
-//======================================================================
-
-double InitialConfiguration::vertex_weight() {
-  double ws = V().Weight[STI];
-  if (!isKink()) ws += V().EBASE;
-  return ws;
-}
-
-//======================================================================
-
-double InitialConfiguration::worm_weight() {
-  double ww;
-  if (INC == DIR::UNDEF) {
-    ww = 1.0;
-  } else {
-    int stype = V().SiteTypeOfLeg[INC];
-    int ttype = Site[stype].TTYPE;
-    SOURCE& T = Source[ttype];
-    ww = T.Weight(XINC, State[INC]);
-  }
-  return ww;
-}
-
-//======================================================================
-
-void InitialConfiguration::dump() {
-  double wv, ww;
-  if (INC == DIR::UNDEF) {
-    wv = 1.0;
-    ww = 1.0;
-  } else {
-    wv = vertex_weight();
-    ww = worm_weight();
-  }
-  printf(" IC(%3d) : nleg= %2d, sti= %4d (", id(), NLEG, STI);
-  for (int i = 0; i < NLEG; i++) {
-    printf(" %1d", State[i]);
-  }
-  printf("), inc= %2d, xinc= %2d ", INC, XINC);
-  printf("--> wv= %8.3f, ww= %8.3f, w= %8.3f\n", wv, ww, wv * ww);
-  if (NCH != 0) {
-    for (int i = 0; i < NCH; i++) {
-      printf("   Channel[%d] : out= %d, xout= %d, prob= %8.3f\n", i,
-             FinalDirection[i], FinalState[i], ScatteringProbability[i]);
-    }
-  }
-}
-
-//======================================================================
-
-void InitialConfiguration::write() {
-  if (!isValid()) return;
-  fprintf(FALG, "\n");
-  fprintf(FALG, "    <InitialConfiguration>\n");
-  fprintf(FALG, "      <State> ");
-  for (int i = 0; i < NLEG; i++) fprintf(FALG, " %d", State[i]);
-  fprintf(FALG, " </State>\n");
-  fprintf(FALG, "      <IncomingDirection> %d </IncomingDirection>\n", INC);
-  fprintf(FALG, "      <NewState> %d </NewState>\n", XINC);
-  fprintf(FALG, "      <NumberOfChannels> %d </NumberOfChannels>\n", NCH);
-  for (int ch = 0; ch < NCH; ch++) {
-    int out = FinalDirection[ch];
-    int xout = FinalState[ch];
-    double p = ScatteringProbability[ch];
-    fprintf(FALG, "      <Channel> %4d %4d %24.16lf </Channel>\n", out, xout,
-            p);
-  }
-  fprintf(FALG, "    </InitialConfiguration>\n");
-}
-
-//######################################################################
-
-void QUANTITY::load(XML::Block& X) {
-  ID = X["QTYPE"].getInteger();
-  NAME = X["Name"].getString();
-  Value.init(2, NSTYPE, NXMAX);
-  Value.set_all(0.0);
-  isDefined.init(2, NSTYPE, NXMAX);
-  isDefined.set_all(false);
-  for (int i = 0; i < X.NumberOfBlocks(); i++) {
-    XML::Block& B = X[i];
-    if (B.getName() == "Value") {
-      int st = B.getInteger(0);
-      int x = B.getInteger(1);
-      double v = B.getDouble(2);
-      isDefined(st, x) = true;
-      Value(st, x) = v;
-    }
-  }
-}
-
-//======================================================================
-
-void QUANTITY::write() {
-  fprintf(FALG, "\n");
-  fprintf(FALG, "  <Quantity>\n");
-  fprintf(FALG, "    <QTYPE> %d </QTYPE>\n", getID());
-  fprintf(FALG, "    <Name> %s </Name>\n", getName().c_str());
-  for (int st = 0; st < NSTYPE; st++) {
-    for (int x = 0; x < NXMAX; x++) {
-      if (isDefined(st, x)) {
-        fprintf(FALG, "    <Value> %d %d %24.16lf </Value>\n", st, x,
-                Value(st, x));
-      }
-    }
-  }
-  fprintf(FALG, "  </Quantity>\n");
-}
-
-//######################################################################
-
-//######################################################################
-
-void SolveWeightEquation(int N, Array<double>& V, Array<double>& W) {
-  /*
-  for (int i=0; i<N; i++) {
-    printf("V[%d] = %8.3f\n", i, V[i]);
-  }
-*/
-
-  W.set_all(0.0);
-
-  int p, q;
-  double V_first;
-  double V_second;
-  double V_third;
-  double V_first_new;
-  double V_second_new;
-  int N_first;   // the number of the largest elements
-  int N_second;  // the number of the second largest
-
-  // First three (unique) weights and indices
-  while (V[1] > EPS) {
-    V_first = V[0];
-    for (p = 0; p < N; p++)
-      if (V[p] != V_first) break;
-    N_first = p;
-    if (p == N) {
-      V_second = 0.0;
-      V_third = 0.0;
-      N_second = 0;
-    } else {
-      V_second = V[p];
-      for (q = p; q < N; q++)
-        if (V[q] != V_second) break;
-      if (q == N) {
-        V_third = 0.0;
-      } else {
-        V_third = V[q];
-      }
-      N_second = q - p;
-    }
-
-    // Calculation w_{ij} from V_i
-
-    double dw1;  // decrement of the first largest element
-    double dw2;  // decrement of the second largest element
-    if (N_first == 1) {
-      // When the maximum weight state is unique
-      // introduce a transition between the max state and the second
-      // and reduce weights of these states
-      double x = V_first - V_second;
-      double y = (double)(N_second - 1) * (V_second - V_third);
-      if (x < y) {
-        dw1 = (V_first - V_second) / (1.0 - 1.0 / (double)(N_second));
-        dw2 = dw1 / (double)N_second;
-        V_second_new = V_second - dw2;
-        V_first_new = V_second_new;
-      } else {
-        dw2 = V_second - V_third;
-        dw1 = dw2 * (double)N_second;
-        V_second_new = V_third;
-        V_first_new = V_first - dw1;
-      }
-      V[0] = V_first_new;
-      for (int i = 1; i < 1 + N_second; i++) {
-        W(0, i) += dw2;
-        W(i, 0) += dw2;
-        V[i] = V_second_new;
-      }
-    } else {
-      // When the maximum weight state is degenerated
-      // introduce a transition between these states
-      // and reduce weights of these states to the weight of the second largest.
-      dw1 = (V_first - V_second) / (double)(N_first - 1);
-      for (int i = 0; i < N_first; i++) {
-        V[i] = V_second;
-        for (int j = 0; j < N_first; j++) {
-          if (i == j) continue;
-          W(i, j) += dw1;
-        }
-      }
-    }
-  }
-  if (V[0] > 0.0) {
-    W(0, 0) += V[0];
-    V[0] = 0.0;
-  }
-}
-
-//======================================================================
-
-void bubble_sort(int N, Array<double>& V, Array<int>& I) {
-  for (int i = 0; i < N; i++) I[i] = i;
-  for (int i = 0; i < N - 1; i++) {
-    for (int j = i + 1; j < N; j++) {
-      if (V[I[i]] < V[I[j]]) {
-        int ii = I[i];
-        I[i] = I[j];
-        I[j] = ii;
-      }
-    }
-  }
-}
diff --git a/src/dla/generators/dla_alg.h b/src/dla/generators/dla_alg.h
deleted file mode 100644
index db27916f..00000000
--- a/src/dla/generators/dla_alg.h
+++ /dev/null
@@ -1,356 +0,0 @@
-//
-//     Algorithm File Generator for dla.cc
-//
-//                   Copy Right 2005, Naoki Kawashima
-//
-
-#include <cmath>
-#include <vector>
-
-#include "array.h"
-#include "io.h"
-#include "name.h"
-#include "xml.h"
-
-class InitialConfiguration;
-class InitialConfigurationGroup;
-class GENERAL;
-class SITE;
-class SOURCE;
-class INTERACTION;
-class VERTEX;
-class QUANTITY;
-
-//######################################################################
-
-void SolveWeightEquation(int N, Array<double>& V, Array<double>& W);
-bool testKink(int ist);
-bool testForbidden(int ist);
-void bubble_sort(int N, Array<double>& V, Array<int>& I);
-
-//######################################################################
-
-class GENERAL {
- public:
-  std::string comment;
-  // "WeightDiagonal"
-  // the artificial weight attached to the diagonal state
-  // in solving the weight equation for the worm creation/annihilation
-  // ( the same as 1/(N_{box} \eta^2) )
-  double WeightDiagonal;
-  GENERAL(XML::Block& X);
-  void write();
-};
-
-//######################################################################
-
-class SITE {
- public:
-  int ID;
-  int TTYPE;  // the SOURCE type of the worm tail
-  int VTYPE;  // the VERTEX type of the worm tail
-  int NX;
-  int* NumberOfChannels;
-  SOURCE* _T;
-  SOURCE& T() { return *_T; };
-  VERTEX& V();
-  Array<int> WormCreationNewState;
-  Array<int> WormCreationDirection;
-  Array<double> WormCreationProbability;
-
-  SITE() {
-    _T = 0;
-    ID = STYPE::UNDEF;
-    TTYPE = TTYPE::UNDEF;
-    VTYPE = VTYPE::UNDEF;
-    NX = 0;
-    NumberOfChannels = 0;
-  };
-
-  ~SITE() {
-    if (NumberOfChannels != 0) delete[] NumberOfChannels;
-    WormCreationNewState.reset();
-    WormCreationDirection.reset();
-    WormCreationProbability.reset();
-  };
-
-  void load(XML::Block& X);
-  void SetInitialHeadTypeProbability();
-  void write();
-  void dump();
-};
-
-//######################################################################
-
-class VERTEX {
- private:
-  IndexSystem* _INDX;
-
- public:
-  int ID;
-  std::vector<InitialConfigurationGroup*> icg;
-  InitialConfigurationGroup& ICG(int i) { return *(icg[i]); };
-  void add_ICG(InitialConfigurationGroup& x) { icg.push_back(&x); };
-
-  int NBODY;  // the number of sites interacting
-  int NLEG;
-  int NICG;      // number of initial configuration group
-  int NICV;      // total number of initial configurations
-  int NST;       // the number of initial states (not including worm type or
-                 // direction)
-  int CATEGORY;  // =0 (TERM), =1 (WORM), =2 (INTERACTION)
-  double EBASE;
-  Array<double> Weight;
-  Array<int> SiteTypeOfLeg;
-  Array<int> checked;
-
-  void setINDX(IndexSystem& I) { _INDX = &I; };
-  IndexSystem& INDX() { return *_INDX; };
-
-  VERTEX() {
-    ID = VTYPE::UNDEF;
-    NBODY = 0;
-    NLEG = 0;
-    NICG = 0;
-    NICV = 0;
-    NST = 0;
-    CATEGORY = VCAT::UNDEF;
-    EBASE = 0.0;
-  };
-
-  ~VERTEX(){};
-
-  bool isTerm() { return CATEGORY == VCAT::TERM; };
-  bool isWorm() { return CATEGORY == VCAT::WORM; };
-  bool isInteraction() { return CATEGORY == VCAT::INT; };
-  void Grouping();
-  double ComputeEBASE();
-  void ComputeScatteringProbability();
-  void dump();
-  void load(XML::Block& X);
-  void write();
-  int NumberOfValidInitialConfigurations();
-  bool testKink(int ist);
-  bool testForbidden(int ist);
-};
-
-//######################################################################
-
-class SOURCE {
- public:
-  int ID;
-  int STYPE;
-  int VTYPE;
-  VERTEX* _V;
-  SOURCE() {
-    ID = TTYPE::UNDEF;
-    STYPE = STYPE::UNDEF;
-    VTYPE = VTYPE::UNDEF;
-    _V = 0;
-  };
-  VERTEX& V() {
-    if (_V == 0) {
-      printf("SOURCE::V> Error. _V is not defined.\n");
-      exit(0);
-    }
-    return *_V;
-  };
-  double Weight(int x0, int x1) { return V().Weight(x0, x1); };
-  void load(XML::Block& X);
-  void dump();
-};
-
-//######################################################################
-
-class INTERACTION {
- public:
-  int ID;
-  int VTYPE;
-  int NBODY;
-  VERTEX* _V;
-  VERTEX& V() { return *_V; };
-  Array<double> VertexDensity;
-  Array<double> Sign;
-
-  INTERACTION() {
-    ID = ITYPE::UNDEF;
-    VTYPE = VTYPE::UNDEF;
-    NBODY = 0;
-  };
-  ~INTERACTION() { VertexDensity.reset(); };
-  int STYPE(int i);
-  void SetVertexDensity();
-  void load(XML::Block& X);
-  void write();
-  void dump();
-};
-
-//######################################################################
-
-class InitialConfiguration {
- private:
-  static int LastID;
-  int ID;
-
- public:
-  friend class InitialConfigurationGroup;
-
-  int NBODY;
-  int NLEG;
-  int NCH;
-  int STI;
-  int XINC;
-  int INC;
-
-  int* State;
-  int* FinalState;
-  int* FinalDirection;
-  double* ScatteringProbability;
-  VERTEX* _V;
-
-  InitialConfiguration() {
-    ID = LastID;
-    LastID++;
-    State = 0;
-    FinalState = 0;
-    FinalDirection = 0;
-    ScatteringProbability = 0;
-  };
-
-  ~InitialConfiguration() {
-    if (State != 0) delete[] State;
-    if (ScatteringProbability != 0) delete[] ScatteringProbability;
-    if (FinalState != 0) delete[] FinalState;
-    if (FinalDirection != 0) delete[] FinalDirection;
-  };
-
-  bool isValid();
-  void setV(VERTEX& V0) { _V = &V0; };
-  VERTEX& V() { return *_V; };
-  void init(VERTEX& V, int sti, int inc, int xinc, int NICmax);
-  bool isKink();
-  int& id() { return ID; }
-  double weight();
-  double vertex_weight();
-  double worm_weight();
-  void dump();
-  void write();
-};
-
-int InitialConfiguration::LastID = 0;
-
-//######################################################################
-
-class InitialConfigurationGroup {
- private:
-  static int LastID;
-  int ID;
-
- public:
-  int NIC;
-  int NDIAG;
-  Array<double> U;
-  //  Array<double> W;
-  //  Array<double> P;
-  VERTEX* _V;
-  Array<InitialConfiguration> IC;
-  int NBODY;
-  int NLEG;
-
-  InitialConfigurationGroup() {
-    ID = LastID;
-    LastID++;
-    _V = 0;
-  };
-
-  ~InitialConfigurationGroup() { IC.reset(); };
-
-  void setV(VERTEX& V0) { _V = &V0; };
-
-  VERTEX& V() {
-    if (_V == 0) {
-      printf("VERTEX::V> Error. _V has not been defined.\n");
-      exit(0);
-    }
-    return *_V;
-  };
-
-  int NumberOfValidInitialConfigurations();
-
-  //  void init_w( int nbody , int sti );
-  void init(VERTEX& v, int nbody, int sti, int inc, int xinc);
-  int number_of_diagonal_states();
-  int number_of_offdiagonal_states();
-  double minimum_diagonal_weight();
-  double maximum_diagonal_weight();
-  double maximum_offdiagonal_weight();
-  double sum_of_all_weights();
-  double ebase();  // the maximum base energy to make
-                   // all diagonal weights positive
-  void add_to_diagonal_element(double E);
-  void numbering(int& ICNF);
-  void ScatteringProbability();
-  void ResetWeight();
-
-  void dump();
-  void write();
-};
-
-int InitialConfigurationGroup::LastID = 0;
-
-//######################################################################
-
-class QUANTITY {
- private:
-  int ID;
-  std::string NAME;
-
- public:
-  Array<double> Value;
-  Array<bool> isDefined;
-  int getID() { return ID; };
-  void setID(int i) { ID = i; };
-  const std::string& getName() { return NAME; };
-  void setName(const std::string& s) { NAME = s; };
-  void load(XML::Block& X);
-  //  void dump();
-  void write();
-};
-
-//######################################################################
-
-int INTERACTION::STYPE(int i) { return _V->SiteTypeOfLeg(2 * i); }
-
-//######################################################################
-
-VERTEX& SITE::V() {
-  if (_T == 0) {
-    printf("SITE::V> Error. _T is not defined.\n");
-    exit(0);
-  }
-  return T().V();
-}
-
-//######################################################################
-
-const char* HFILE;  // file name of the hamiltonian and the worm
-const char* AFILE;  // file name of algorithm data file
-FileReader H;
-FILE* FALG;
-
-SITE* Site;
-SOURCE* Source;
-INTERACTION* Interaction;
-VERTEX* Vertex;
-QUANTITY* Quantity;
-
-bool isKink(int st);
-
-int NVTYPE;
-int NSTYPE;  // the number of site types
-int NITYPE;  // the number of interaction types
-int NHTYPE;  // the number of worm head types
-int NQTYPE;  // the number of observable types
-int NXMAX;
-
-bool MON = false;  // true for monitoring
diff --git a/src/dla/generators/exact_B.cc b/src/dla/generators/exact_B.cc
deleted file mode 100644
index 4156f4d8..00000000
--- a/src/dla/generators/exact_B.cc
+++ /dev/null
@@ -1,216 +0,0 @@
-#include <cstdio>
-#include <cstring>
-#include <iostream>
-#include <string>
-using namespace std;
-#include "boson_B.h"
-#include "canonical.h"
-#include "matrix.h"
-
-//----------------------------------------------------------------------
-
-// Bose-Hubbard Model
-// H= - J\sum_{<i,j>} b_i b_j + V\sum_{<i,j>} n_i n_j
-//    + u/2z \sum_{<i,j>}( n_i (n_i + 1) + n_j (n_j + 1) )
-//     - mu/z \sum_{<i,j>}( n_i + n_j ).
-// F=mu/z,
-// U=u/z
-
-class BoseHubbardModel {
- public:
-  int M;  // Number of bosons ( [M]-representaion )
-  int NSITE;
-  double J;  // the bilinear coupling for XY
-  double V;  // the bilinear coupling for Z
-  double U;  // the on-site coupling for Z
-  double F;  // the field
-  int DIM;   // dimension of the Hilbert space
-  dgematrix H;
-  dgematrix MXU;
-  dgematrix MZU;
-  dgematrix MXS;
-  dgematrix MZS;
-  dgematrix I;
-
-  BoseHubbardModel(int M0, int NSITE0, double J0, double V0, double U0,
-                   double F0) {
-    M = M0;
-    NSITE = NSITE0;
-    J = J0;
-    V = V0;
-    U = U0;
-    F = F0;
-    // double Vh=V*0.5;
-
-    BosonOperatorSet S(M, NSITE);
-
-    DIM = S.DIM;
-    cmatrix h(DIM);
-    h.zero();
-    printf("  defining the hamiltonian...\n");
-    for (int k = 0; k < NSITE; k++) {
-      int l = (k + 1) % NSITE;
-      printf("    k= %d\n", k);
-      // h += (-J) * ( S.X[k] * S.X[l] + S.Y[k] * S.Y[l] );
-      h += (-J) * (S.UP[k] * S.DN[l]);
-      h += (+V) * (S.Z[k] * S.Z[l]);
-    }
-    for (int k = 0; k < NSITE; k++) {
-      h += (+U) * (S.Z[k] * S.Z[k]);
-      h += (-F - U) * (S.Z[k]);
-    }
-    printf("  ... done.\n");
-
-    printf("  defining the magnetizations...\n");
-    cmatrix mxu(DIM);
-    cmatrix mzu(DIM);
-    cmatrix mxs(DIM);
-    cmatrix mzs(DIM);
-    mxu.zero();
-    mzu.zero();
-    mxs.zero();
-    mzs.zero();
-    for (int k = 0; k < NSITE; k++) {
-      printf("   k= %d\n", k);
-      double sgn = 1.0;
-      if (k % 2 == 1) sgn = -1.0;
-      mxu += (1.0 * S.X[k]);
-      mzu += (1.0 * S.Z[k]);
-      mxs += (sgn * S.X[k]);
-      mzs += (sgn * S.Z[k]);
-    }
-    printf("  ... done.\n");
-
-    H = h.re;
-    MXU = mxu.re;
-    MZU = mzu.re;
-    MXS = mxs.re;
-    MZS = mzs.re;
-    I = S.I.re;
-  };
-};
-
-//============================================================================
-
-void Average(int DIM, dgematrix& R, dgematrix& Q, double& Ave, double& Var) {
-  dgematrix W(DIM, DIM);
-  W = Q * Q;
-  Ave = CanonicalAverage(R, Q);
-  Var = CanonicalAverage(R, W);
-  Var = Var - Ave * Ave;
-}
-
-//============================================================================
-
-void WriteXML(int M, dgematrix& Q, dgematrix& H, std::string const& filename) {
-  FILE* FOUT = fopen(filename.c_str(), "w");
-  int D = M + 1;
-  int DD = D * D;
-  fprintf(FOUT, "<Hamiltonian>\n");
-  fprintf(FOUT, "  <General>\n");
-  fprintf(FOUT,
-          "    <Comment> Extended Bose-Hubbard model with NMAX=%d </Comment>\n",
-          M);
-  fprintf(FOUT, "    <NSTYPE> 1 </NSTYPE>\n");
-  fprintf(FOUT, "    <NITYPE> 1 </NITYPE>\n");
-  fprintf(FOUT, "    <NXMAX>  %d </NXMAX>\n", D);
-  fprintf(FOUT, "  </General>\n");
-  fprintf(FOUT, "\n");
-  fprintf(FOUT, "  <Site>\n");
-  fprintf(FOUT, "    <STYPE> 0 </STYPE>\n");
-  fprintf(FOUT, "    <TTYPE> 0 </TTYPE>\n");
-  fprintf(FOUT, "    <NX>   %d </NX>\n", D);
-  fprintf(FOUT, "  </Site>\n");
-  fprintf(FOUT, "\n");
-  fprintf(FOUT, "  <Source>\n");
-  fprintf(FOUT, "    <TTYPE> 0 </TTYPE>\n");
-  fprintf(FOUT, "    <STYPE> 0 </STYPE>\n");
-  for (int i = 0; i < D; i++) {
-    for (int j = 0; j < D; j++) {
-      double x = Q(i, j);
-      if (abs(x) > 1.0e-8) {
-        fprintf(FOUT, "    <Weight> %d %d %24.16f </Weight>\n", i, j, x);
-      }
-    }
-  }
-  fprintf(FOUT, "  </Source>\n");
-  fprintf(FOUT, "\n");
-  fprintf(FOUT, "  <Interaction>\n");
-  fprintf(FOUT, "    <ITYPE> 0 </ITYPE>\n");
-  fprintf(FOUT, "    <NBODY> 2 </NBODY>\n");
-  fprintf(FOUT, "    <STYPE> 0 0 </STYPE>\n");
-  for (int i = 0; i < DD; i++) {
-    int i0 = i % D;
-    int i1 = i / D;
-    for (int j = 0; j < DD; j++) {
-      int j0 = j % D;
-      int j1 = j / D;
-      double x = -H(i, j);
-      if (abs(x) > 1.0e-8) {
-        // if (i != j) x = abs(x);
-        fprintf(FOUT, "    <Weight> %d %d %d %d %24.16f </Weight>\n", i0, j0,
-                i1, j1, x);
-      }
-    }
-  }
-  fprintf(FOUT, "  </Interaction>\n");
-  fprintf(FOUT, "</Hamiltonian>\n");
-}
-
-void ShowUsage(std::string const& exename) {
-  printf("usage:\n");
-  printf("    %s [-o filename] M J V U F\n", exename.c_str());
-  printf("arguments:\n");
-  printf("    M ... the maxmum number of bosons on each site \n");
-  printf("    J ... the hopping constant (positive)\n");
-  printf(
-      "    V ... the nearest neighbor interaction (positive for "
-      "replusive)\n");
-  printf("    U ... the on-site interaction (positive for replusive)\n");
-  printf("          ( = u/z if the field u is the field per site.)\n");
-  printf("    F ... the chemical potential in the pair Hamiltonian\n");
-  printf("          ( = H/z if the field H is the field per site\n");
-  printf("            and H is shared equally by all pairs,\n");
-  printf("            e.g., F = H/4 for a square lattice. )\n");
-  printf("option:\n");
-  printf("    -o filename ... output file (default: hamiltonian.xml)\n");
-}
-
-//============================================================================
-//    Main
-//============================================================================
-
-int main(int argc, char** argv) {
-  string exename = argv[0];
-  string filename("hamiltonian.xml");
-  if (argc < 3) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  if (std::strcmp(argv[1], "-o") == 0) {
-    filename = argv[2];
-    argc -= 2;
-    argv += 2;
-  }
-
-  int NARG = 5;
-  if (argc != NARG + 1) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  int M = atoi(argv[1]);
-  double J = (double)atof(argv[2]);
-  double V = (double)atof(argv[3]);
-  double U = (double)atof(argv[4]);
-  double F = (double)atof(argv[5]);
-  printf(" M     = %4d\n", M);
-  printf(" J     = %8.3f\n", J);
-  printf(" V     = %8.3f\n", V);
-  printf(" U     = %8.3f\n", U);
-  printf(" F     = %8.3f\n", F);
-  BosonOperator S(M);
-  BoseHubbardModel MDL(M, 2, J, 0.5 * V, 0.5 * U, F);
-  WriteXML(M, S.X.re, MDL.H, filename);
-
-  return 0;
-}
diff --git a/src/dla/generators/exact_H.cc b/src/dla/generators/exact_H.cc
deleted file mode 100644
index b2feed9b..00000000
--- a/src/dla/generators/exact_H.cc
+++ /dev/null
@@ -1,193 +0,0 @@
-#include <cstdio>
-#include <cstring>
-#include <iostream>
-#include <string>
-using namespace std;
-#include "canonical.h"
-#include "matrix.h"
-#include "spin_H.h"
-
-//----------------------------------------------------------------------
-
-// Heisenberg Model
-
-class HeisenbergModel {
- public:
-  int M;  // Number of bosons ( [M]-representaion )
-  int NSITE;
-  double J;  // the bilinear coupling
-  double F;  // the field
-  int DIM;   // dimension of the Hilbert space
-  dgematrix H;
-  dgematrix MXU;
-  dgematrix MZU;
-  dgematrix MXS;
-  dgematrix MZS;
-  dgematrix I;
-
-  HeisenbergModel(int M0, int NSITE0, double J0, double F0) {
-    M = M0;
-    NSITE = NSITE0;
-    J = J0;
-    F = F0;
-
-    HeisenbergSpinSet S(M, NSITE);
-
-    DIM = S.DIM;
-    cmatrix h(DIM);
-    h.zero();
-    printf("  defining the hamiltonian...\n");
-    for (int k = 0; k < NSITE; k++) {
-      int l = (k + 1) % NSITE;
-      printf("    k= %d\n", k);
-      h += (-J) * (S.X[k] * S.X[l] + S.Y[k] * S.Y[l] + S.Z[k] * S.Z[l]);
-    }
-    for (int k = 0; k < NSITE; k++) {
-      h += (-F) * (S.Z[k]);
-    }
-    printf("  ... done.\n");
-
-    printf("  defining the magnetizations...\n");
-    cmatrix mxu(DIM);
-    cmatrix mzu(DIM);
-    cmatrix mxs(DIM);
-    cmatrix mzs(DIM);
-    mxu.zero();
-    mzu.zero();
-    mxs.zero();
-    mzs.zero();
-    for (int k = 0; k < NSITE; k++) {
-      printf("   k= %d\n", k);
-      double sgn = 1.0;
-      if (k % 2 == 1) sgn = -1.0;
-      mxu += (1.0 * S.X[k]);
-      mzu += (1.0 * S.Z[k]);
-      mxs += (sgn * S.X[k]);
-      mzs += (sgn * S.Z[k]);
-    }
-    printf("  ... done.\n");
-
-    H = h.re;
-    MXU = mxu.re;
-    MZU = mzu.re;
-    MXS = mxs.re;
-    MZS = mzs.re;
-    I = S.I.re;
-  };
-};
-
-//============================================================================
-
-void Average(int DIM, dgematrix& R, dgematrix& Q, double& Ave, double& Var) {
-  dgematrix W(DIM, DIM);
-  W = Q * Q;
-  Ave = CanonicalAverage(R, Q);
-  Var = CanonicalAverage(R, W);
-  Var = Var - Ave * Ave;
-}
-
-//============================================================================
-
-void WriteXML(int M, dgematrix& Q, dgematrix& H, std::string const& filename) {
-  FILE* FOUT = fopen(filename.c_str(), "w");
-  int D = M + 1;
-  int DD = D * D;
-  fprintf(FOUT, "<Hamiltonian>\n");
-  fprintf(FOUT, "  <General>\n");
-  fprintf(FOUT, "    <Comment> SU(2) Heisenberg model with S=%d/2 </Comment>\n",
-          M);
-  fprintf(FOUT, "    <NSTYPE> 1 </NSTYPE>\n");
-  fprintf(FOUT, "    <NITYPE> 1 </NITYPE>\n");
-  fprintf(FOUT, "    <NXMAX>  %d </NXMAX>\n", D);
-  fprintf(FOUT, "  </General>\n");
-  fprintf(FOUT, "\n");
-  fprintf(FOUT, "  <Site>\n");
-  fprintf(FOUT, "    <STYPE> 0 </STYPE>\n");
-  fprintf(FOUT, "    <TTYPE> 0 </TTYPE>\n");
-  fprintf(FOUT, "    <NX>   %d </NX>\n", D);
-  fprintf(FOUT, "  </Site>\n");
-  fprintf(FOUT, "\n");
-  fprintf(FOUT, "  <Source>\n");
-  fprintf(FOUT, "    <TTYPE> 0 </TTYPE>\n");
-  fprintf(FOUT, "    <STYPE> 0 </STYPE>\n");
-  for (int i = 0; i < D; i++) {
-    for (int j = 0; j < D; j++) {
-      double x = Q(i, j);
-      if (abs(x) > 1.0e-8) {
-        fprintf(FOUT, "    <Weight> %d %d %24.16f </Weight>\n", i, j, x);
-      }
-    }
-  }
-  fprintf(FOUT, "  </Source>\n");
-  fprintf(FOUT, "\n");
-  fprintf(FOUT, "  <Interaction>\n");
-  fprintf(FOUT, "    <ITYPE> 0 </ITYPE>\n");
-  fprintf(FOUT, "    <NBODY> 2 </NBODY>\n");
-  fprintf(FOUT, "    <STYPE> 0 0 </STYPE>\n");
-  for (int i = 0; i < DD; i++) {
-    int i0 = i % D;
-    int i1 = i / D;
-    for (int j = 0; j < DD; j++) {
-      int j0 = j % D;
-      int j1 = j / D;
-      double x = -H(i, j);
-      if (abs(x) > 1.0e-8) {
-        // if (i != j) x = abs(x);
-        fprintf(FOUT, "    <Weight> %d %d %d %d %24.16f </Weight>\n", i0, j0,
-                i1, j1, x);
-      }
-    }
-  }
-  fprintf(FOUT, "  </Interaction>\n");
-  fprintf(FOUT, "</Hamiltonian>\n");
-}
-
-void ShowUsage(std::string const& exename) {
-  printf("usage:\n");
-  printf("    %s [-o filename] M J F\n", exename.c_str());
-  printf("arguments:\n");
-  printf("    M ... the number of bosons on each site \n");
-  printf("            ( M= 1, 2, 3, ... for S= 1/2, 1, 3/2, ... ) \n");
-  printf("    J ... the coupling constant (positive for ferromagnets)\n");
-  printf("    F ... the magnetic field in the pair Hamiltonian\n");
-  printf("          ( = H/z if the field H is the field per site\n");
-  printf("            and H is shared equally by all pairs,\n");
-  printf("            e.g., F = H/4 for a square lattice. )\n");
-  printf("option:\n");
-  printf("    -o filename ... output file (default: hamiltonian.xml)\n");
-}
-
-//============================================================================
-//    Main
-//============================================================================
-
-int main(int argc, char** argv) {
-  string exename = argv[0];
-  string filename("hamiltonian.xml");
-
-  if (argc < 3) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  if (std::strcmp(argv[1], "-o") == 0) {
-    filename = argv[2];
-    argc -= 2;
-    argv += 2;
-  }
-
-  int NARG = 3;
-  if (argc != NARG + 1) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  int M = atoi(argv[1]);
-  double J = (double)atof(argv[2]);
-  double F = (double)atof(argv[3]);
-  printf(" M     = %4d\n", M);
-  printf(" J     = %8.3f\n", J);
-  printf(" F     = %8.3f\n", F);
-  HeisenbergSpin S(M);
-  HeisenbergModel MDL(M, 2, 0.5 * J, F);
-  WriteXML(M, S.X.re, MDL.H, filename);
-  return 0;
-}
diff --git a/src/dla/generators/lattgene_C.cc b/src/dla/generators/lattgene_C.cc
deleted file mode 100644
index 82372c71..00000000
--- a/src/dla/generators/lattgene_C.cc
+++ /dev/null
@@ -1,211 +0,0 @@
-/*---------------------------------------------
-
-   Generating lattice.xml for a square lattice.
-
-----------------------------------------------*/
-
-#include <cmath>
-#include <cstdlib>
-#include <cstring>
-#include <fstream>
-#include <iostream>
-#include <string>
-
-using namespace std;
-
-//--------------------------------------------------------------
-void ShowUsage(std::string const& exename) {
-  cout << "usage:\n";
-  cout << "    " << exename << " [-o output] D L\n";
-  cout << "    " << exename << " [-o output] D L_1 ... L_D\n";
-  cout << "arguments:\n";
-  cout << "    D ... dimension.\n";
-  cout << "    L ... the liner size of the lattice. \n";
-  cout << "          ( L, L1, ... must be even number. )\n";
-  cout << "option:\n";
-  cout << "    -o output ... output file (default: lattice.xml)\n";
-}
-//-------------------------------------------------------------
-void WriteXML(int D, int L[], std::string const& filename) {
-  ofstream fout(filename.c_str());
-  fout.precision(15);
-  int N = 1;  // number of sites.
-  for (int i = 0; i < D; i++) {
-    N *= L[i];
-  }
-
-  int NumberOfInteractions = N * D;
-  int NumberOfSiteTypes = 1;
-  int NumberOfInteractionTypes = 1;
-  int NumberOfEdgeInteractions = 0;
-
-  // Hypercubic
-  int* Lext = new int[D];
-
-  for (int di = 0; di < D; di++) {
-    Lext[di] = 1;
-    for (int dj = 0; dj < D; dj++) {
-      if (di != dj) Lext[di] *= L[dj];
-    }
-  }
-
-  for (int di = 0; di < D; di++) NumberOfEdgeInteractions += Lext[di];
-
-  int BD = D;  // Hypercubic
-
-  fout << "<LATTICE>" << endl << endl;
-  fout << "<Comment>" << endl;
-  fout << "  " << D << "-dimension square lattice" << endl;
-  fout << "</Comment>" << endl << endl;
-
-  fout << "<Dimension> " << D << " </Dimension>" << endl;
-  fout << "<BondDimension> " << BD << " </BondDimension>" << endl;
-  fout << "<LinearSize> ";
-  for (int i = 0; i < D; i++) {
-    fout << L[i] << " ";
-  }
-  fout << "</LinearSize>" << endl;
-  fout << "<NumberOfSites> " << N << " </NumberOfSites>" << endl;
-  fout << "<NumberOfInteractions> " << NumberOfInteractions
-       << " </NumberOfInteractions>" << endl;
-  fout << "<NumberOfSiteTypes> " << NumberOfSiteTypes << " </NumberOfSiteTypes>"
-       << endl;
-  fout << "<NumberOfInteractionTypes> " << NumberOfInteractionTypes
-       << " </NumberOfInteractionTypes>" << endl;
-  fout << "<NumberOfEdgeInteractions> " << NumberOfEdgeInteractions
-       << " </NumberOfEdgeInteractions>" << endl;
-
-  fout << endl;
-
-  fout << "<!-- <S> [id] [stype] [mtype] </S> -->" << endl;
-
-  int stype = 0;
-  int mtype = 0;
-
-  for (int id = 0; id < N; id++) {
-    int p = id;
-    int Nt = N;
-    mtype = 0;
-    for (int q = D - 1; q >= 0; q--) {
-      Nt /= L[q];
-      mtype += p / Nt;
-      p = p % Nt;
-    }
-
-    mtype = mtype % 2;
-    fout << "<S> " << id << " " << stype << " " << mtype << " </S>" << endl;
-  }
-
-  fout << endl;
-  fout << "<!-- <I> [id] [itype] [nbody] [s0] [s1] ... [eid] [edim] </I> -->"
-       << endl
-       << endl;
-
-  int NB = D * N;  // number of bonds
-  int* x = new int[D];
-  int itype = 0;
-  int nbody = 2;
-  NB = 0;
-  int etype;
-  int eid = 0;
-
-  for (int i = 0; i < N; i++) {
-    for (int p = 0; p < D; p++) {
-      int k = i;
-      for (int q = 0; q < D; q++) {
-        x[q] = k % L[q];
-        k /= L[q];
-      }
-
-      if (x[p] == L[p] - 1) {
-        etype = eid;
-        eid++;
-      } else
-        etype = -1;
-
-      x[p] = (x[p] + 1) % L[p];
-      int j = 0;
-      for (int q = D - 1; q >= 0; q--) {
-        j *= L[q];
-        j += x[q];
-      }
-
-      fout << "<I> " << NB << " " << itype << " " << nbody << " " << i << " "
-           << j << " " << etype << " " << p << " </I>" << endl;
-
-      NB++;
-    }
-  }
-
-  fout << endl;
-  fout << "<!-- <V> [bond dim] [unit vectorx]...  </V> -->" << endl << endl;
-  for (int p = 0; p < BD; p++) {
-    fout << "<V> " << p << " ";
-    for (int q = 0; q < D; q++) {
-      if (p == q)
-        fout << 1 << " ";
-      else
-        fout << 0 << " ";
-    }
-    fout << "</V>" << endl;
-  }
-  fout << "</LATTICE>" << endl;
-
-  delete[] x;
-  delete[] Lext;
-}
-//--------------------------------------------------------------
-
-int main(int argc, char** argv) {
-  std::string exename(argv[0]);
-  if (argc < 3) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  std::string filename("lattice.xml");
-  if (std::strcmp(argv[1], "-o") == 0) {
-    filename = argv[2];
-    argc -= 2;
-    argv = argv + 2;
-  }
-  int NARG = 3;
-  if (argc < NARG) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  const int D = atoi(argv[1]);
-  //  int        L[D] ;
-  int* L = new int[D];  // edit sakakura
-
-  if (argc == NARG) {
-    int lx = atoi(argv[2]);
-    for (int i = 0; i < D; i++) {
-      L[i] = lx;
-    }
-  } else if (argc == D + NARG - 1) {
-    for (int i = 0; i < D; i++) {
-      L[i] = atoi(argv[2 + i]);
-    }
-  } else {
-    cout << "error: D != number of L[]." << endl;
-    ShowUsage(exename);
-    exit(0);
-  }
-
-  int EvenOrOdd = 0;
-  cout.precision(15);
-  cout << "D     = " << D << endl;
-  for (int i = 0; i < D; i++) {
-    cout << "L" << i << "    = " << L[i] << endl;
-    EvenOrOdd += L[i] % 2;
-  }
-
-  if (EvenOrOdd) {
-    cout << "Warnig: L should be an even number." << endl;
-  }
-
-  WriteXML(D, L, filename);
-  cout << "... done." << endl;
-  delete[] L;
-  return 0;
-}
diff --git a/src/dla/generators/lattgene_T.cc b/src/dla/generators/lattgene_T.cc
deleted file mode 100644
index 18fe7215..00000000
--- a/src/dla/generators/lattgene_T.cc
+++ /dev/null
@@ -1,218 +0,0 @@
-/*---------------------------------------------
-
-   Generating lattice.xml for a triangular lattice.
-
-----------------------------------------------*/
-
-#include <cmath>
-#include <cstdlib>
-#include <cstring>
-#include <fstream>
-#include <iostream>
-#include <string>
-#include <vector>
-
-using namespace std;
-
-void ShowUsage(std::string const& exename) {
-  cout << "usage:\n";
-  cout << "    " << exename << " [-o output] Lx [Ly]\n";
-  cout << "arguments:\n";
-  cout << "    Lx ... the liner size of the lattice in x direction. \n";
-  cout << "    Ly ... the liner size of the lattice in y direction. \n";
-  cout << "           if omitted, Ly will be the same as Lx\n";
-  cout << "option:\n";
-  cout << "    -o output ... output file (default: lattice.xml)\n";
-}
-void WriteXML(std::vector<int> const& L, std::string const& filename) {
-  const int D = L.size();
-  ofstream fout(filename.c_str());
-  fout.precision(15);
-  int N = 1;  // number of sites.
-  for (int i = 0; i < D; i++) {
-    N *= L[i];
-  }
-
-  int BD = 3;  // Triangular
-
-  int NumberOfInteractions = N * BD;
-  int NumberOfSiteTypes = 1;
-  int NumberOfInteractionTypes = 1;
-  int NumberOfEdgeInteractions = 0;
-
-  std::vector<int> Lext(D);
-
-  for (int di = 0; di < D; di++) {
-    Lext[di] = 1;
-    for (int dj = 0; dj < D; dj++) {
-      if (di != dj) Lext[di] *= 2 * L[dj];
-    }
-  }
-
-  for (int di = 0; di < D; di++) {
-    NumberOfEdgeInteractions += Lext[di];
-  }
-
-  NumberOfEdgeInteractions--;  // double count of the corner site.
-
-  fout << "<LATTICE>" << endl << endl;
-  fout << "<Comment>" << endl;
-  fout << "  " << D << "-dimension triangular lattice" << endl;
-  fout << "</Comment>" << endl << endl;
-
-  fout << "<Dimension> " << D << " </Dimension>" << endl;
-  fout << "<BondDimension> " << BD << " </BondDimension>" << endl;
-  fout << "<LinearSize> ";
-  for (int i = 0; i < D; i++) {
-    fout << L[i] << " ";
-  }
-  fout << "</LinearSize>" << endl;
-  fout << "<NumberOfSites> " << N << " </NumberOfSites>" << endl;
-  fout << "<NumberOfInteractions> " << NumberOfInteractions
-       << " </NumberOfInteractions>" << endl;
-  fout << "<NumberOfSiteTypes> " << NumberOfSiteTypes << " </NumberOfSiteTypes>"
-       << endl;
-  fout << "<NumberOfInteractionTypes> " << NumberOfInteractionTypes
-       << " </NumberOfInteractionTypes>" << endl;
-  fout << "<NumberOfEdgeInteractions> " << NumberOfEdgeInteractions
-       << " </NumberOfEdgeInteractions>" << endl;
-
-  fout << endl;
-
-  fout << "<!-- <S> [id] [stype] [mtype] </S> -->" << endl;
-
-  std::vector<int> x(D);
-  for (int id = 0; id < N; id++) {
-    x[0] = id % L[0];
-    x[1] = id / L[0];
-    int stype = 0;
-    int mtype = (x[0] + x[1]) % 3;
-
-    fout << "<S> " << id << " " << stype << " " << mtype << " </S>" << endl;
-  }
-
-  fout << endl;
-  fout << "<!-- <I> [id] [itype] [nbody] [s0] [s1] ... [eid] [edim] </I> -->"
-       << endl
-       << endl;
-
-  int NB = 0;  // number of bonds
-  int itype = 0;
-  int nbody = 2;
-  int etype;
-  int eid = 0;
-  int p;
-  int j;
-
-  for (int i = 0; i < N; i++) {
-    int k = i;
-    x[0] = k % L[0];
-    x[1] = k / L[0];
-
-    etype = -1;
-    p = -1;
-    j = (x[0] + 1) % L[0] + x[1] * L[0];
-    if (x[0] == L[0] - 1) {
-      etype = eid;
-      eid++;
-      p = 0;
-    }
-    fout << "<I> " << NB << " " << itype << " " << nbody << " " << i << " " << j
-         << " " << etype << " " << p << " </I>" << endl;
-    NB++;
-
-    etype = -1;
-    p = -1;
-    j = x[0] + ((x[1] + 1) % L[1]) * L[0];
-    if (x[1] == L[1] - 1) {
-      etype = eid;
-      eid++;
-      p = 1;
-    }
-    fout << "<I> " << NB << " " << itype << " " << nbody << " " << i << " " << j
-         << " " << etype << " " << p << " </I>" << endl;
-    NB++;
-
-    etype = -1;
-    p = -1;
-    j = (x[0] + 1) % L[0] + ((x[1] + 1) % L[1]) * L[0];
-    if (x[0] == L[0] - 1 && x[1] == L[1] - 1) {
-      etype = eid;
-      eid++;
-      p = 2;
-    } else if (x[0] == L[0] - 1) {
-      etype = eid;
-      eid++;
-      p = 0;
-    } else if (x[1] == L[1] - 1) {
-      etype = eid;
-      eid++;
-      p = 1;
-    }
-
-    fout << "<I> " << NB << " " << itype << " " << nbody << " " << i << " " << j
-         << " " << etype << " " << p << " </I>" << endl;
-    NB++;
-  }
-
-  fout << endl;
-  fout << "<!-- <V> [bond dim] [unit vectorx]...  </V> -->" << endl << endl;
-
-  fout << "<V> " << 0 << " " << 1 << " " << 0 << " </V>" << endl;
-  fout << "<V> " << 1 << " " << 0 << " " << 1 << " </V>" << endl;
-  fout << "<V> " << 2 << " " << 1 << " " << 1 << " </V>" << endl;
-
-  fout << endl;
-  fout << "</LATTICE>" << endl;
-}
-//--------------------------------------------------------------
-
-int main(int argc, char** argv) {
-  std::string exename(argv[0]);
-  std::string filename("lattice.xml");
-  if (argc < 3) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  if (std::strcmp(argv[1], "-o") == 0) {
-    filename = argv[2];
-    argc -= 2;
-    argv = argv + 2;
-  }
-  int NARG = 2;
-  if (argc < NARG) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  const int D = 2;
-  std::vector<int> L(D);
-
-  if (argc == NARG) {
-    int lx = atoi(argv[1]);
-    for (int i = 0; i < D; i++) {
-      L[i] = lx;
-    }
-  } else if (argc == D + NARG - 1) {
-    for (int i = 0; i < D; i++) {
-      L[i] = atoi(argv[1 + i]);
-    }
-  } else {
-    cout << "error: D != number of L[]." << endl;
-    ShowUsage(exename);
-    exit(0);
-  }
-
-  int EvenOrOdd = 0;
-  cout.precision(15);
-  cout << "D     = " << D << endl;
-  for (int i = 0; i < D; i++) {
-    cout << "L" << i << "    = " << L[i] << endl;
-    EvenOrOdd += L[i] % 2;
-  }
-
-  //  if( EvenOrOdd ) { cout<<"Warnig: L should be an even number."<<endl;}
-
-  WriteXML(L, filename);
-  cout << "... done." << endl;
-  return 0;
-}
diff --git a/src/dla/generators/matrix.h b/src/dla/generators/matrix.h
deleted file mode 100644
index 4e3fd960..00000000
--- a/src/dla/generators/matrix.h
+++ /dev/null
@@ -1,618 +0,0 @@
-
-//============================================================================
-//    Exact Calculation of Finite Size Spin Systems
-//============================================================================
-
-#include <stdint.h>
-
-#include <algorithm>
-#include <cfloat>
-#include <cmath>
-#include <complex>
-#include <cstdlib>
-#include <fstream>
-#include <map>
-#include <string>
-#include <vector>  //only for ?geev, ?gegv, etc.
-
-//#include <ctime>
-
-//=============================================================================
-#ifdef __INTEL_COMPILER
-#define DSQSS_INT MKL_INT
-#define MKL_Complex16 std::complex<double>
-#else  //__INTEL_COMPILER
-#define DSQSS_INT int
-#endif  //__INTEL_COMPILER
-
-//=============================================================================
-
-#ifdef __INTEL_COMPILER
-#include <mkl.h>
-#define dsyev_ DSYEV
-#define dgemm_ DGEMM
-#else   //__INTEL_COMPILER
-extern "C" {
-
-void dsyev_(const char* jobz, const char* uplo, const DSQSS_INT* N, double* a,
-            const DSQSS_INT* lda, double* w, double* work,
-            const DSQSS_INT* lwork, DSQSS_INT* info);
-
-void dgemm_(const char* transa, const char* transb, const DSQSS_INT* M,
-            const DSQSS_INT* N, const DSQSS_INT* k, const double* alpha,
-            const double* a, const DSQSS_INT* lda, const double* b,
-            const DSQSS_INT* ldb, const double* beta, double* c,
-            const DSQSS_INT* ldc);
-}
-#endif  //__INTEL_COMPILER
-
-//============================================================================
-class dgematrix {
- private:
-  vector<double> index;
-
- public:
-  int n, m;
-
-  ////////////////
-  dgematrix(int _n, int _m) {
-    resize(_m, _n);
-    for (int i = 0; i < m * n; i++) {
-      index[i] = 0.0;
-    }
-  };
-
-  dgematrix() {
-    n = 0;
-    m = 0;
-  };
-  //  inline dgematrix(const _dgematrix&);
-  ~dgematrix() { index.clear(); };
-  //////////////////
-
-  void resize(int _m, int _n) {
-    n = _n;
-    m = _m;
-    index.resize(_m * _n);
-  };
-
-  inline void clear() { index.clear(); };
-
-  inline dgematrix& zero() {
-    for (long i = 0; i < m * n; i++) {
-      index[i] = 0.0;
-    }
-    return *this;
-  }
-
-  inline dgematrix& identity() {
-    for (int i = 0; i < m * n; i++) {
-      index[i] = 0.0;
-    }
-    for (int i = 0; i < m; i++) {
-      operator()(i, i) = 1.0;
-    }
-    return *this;
-  }
-
-  inline dgematrix& operator*=(const dgematrix&);
-  inline double& operator()(const int& i, const int& j) {
-    return index[i + j * m];
-  };
-  inline double operator()(const int& i, const int& j) const {
-    return index[i + j * m];
-  };
-  inline dgematrix& operator=(const dgematrix& A) {
-    resize(A.m, A.n);
-
-    for (int i = 0; i < A.m; i++) {
-      for (int j = 0; j < A.n; j++) {
-        index[i + j * m] = A(i, j);
-      }
-    }
-
-    return *this;
-  };
-  inline dgematrix& operator+=(const dgematrix& A) {
-    for (int i = 0; i < A.m; i++) {
-      for (int j = 0; j < A.n; j++) {
-        index[i + j * m] += A(i, j);
-      }
-    }
-
-    return *this;
-  };
-  inline dgematrix& operator-=(const dgematrix& A) {
-    for (int i = 0; i < A.m; i++) {
-      for (int j = 0; j < A.n; j++) {
-        index[i + j * m] -= A(i, j);
-      }
-    }
-
-    return *this;
-  };
-  ////////////
-};
-
-// transposed matrix
-inline dgematrix t(const dgematrix& A) {
-  dgematrix C(A.n, A.m);
-
-  for (int i = 0; i < C.m; i++) {
-    for (int j = 0; j < C.n; j++) {
-      C(i, j) = A(j, i);
-    }
-  }
-
-  return C;
-};
-
-//////////////
-
-dgematrix operator+(const dgematrix& A, const dgematrix& B) {
-  dgematrix C(A.m, A.n);
-
-  for (int i = 0; i < A.m; i++) {
-    for (int j = 0; j < A.n; j++) {
-      C(i, j) = A(i, j) + B(i, j);
-    }
-  }
-
-  return C;
-};
-
-dgematrix operator-(const dgematrix& A, const dgematrix& B) {
-  dgematrix C(A.m, A.n);
-
-  for (int i = 0; i < A.m; i++) {
-    for (int j = 0; j < A.n; j++) {
-      C(i, j) = A(i, j) - B(i, j);
-    }
-  }
-
-  return C;
-};
-dgematrix operator*(const double& a, const dgematrix& A) {
-  dgematrix C(A.m, A.n);
-
-  for (int i = 0; i < A.m; i++) {
-    for (int j = 0; j < A.n; j++) {
-      C(i, j) = a * A(i, j);
-    }
-  }
-
-  return C;
-};
-
-dgematrix operator*(const dgematrix& dA, const dgematrix& dB) {
-  DSQSS_INT M = dA.m;
-  DSQSS_INT N = dB.n;
-  DSQSS_INT K = dA.n;
-  double *A, *B, *C;
-
-  A = new double[dA.m * dA.n];
-  B = new double[dB.m * dB.n];
-  C = new double[dA.m * dB.n];
-
-  dgematrix dC(dA.m, dB.n);
-
-  double alpha = 1.0, beta = 0.0;
-
-  char transa = 'n';
-  char transb = 'n';
-
-  for (int i = 0; i < dA.m; i++) {
-    for (int j = 0; j < dA.n; j++) {
-      A[i + j * dA.m] = dA(i, j);
-    }
-  }
-
-  for (int i = 0; i < dB.m; i++) {
-    for (int j = 0; j < dB.n; j++) {
-      B[i + j * dB.m] = dB(i, j);
-    }
-  }
-
-  DSQSS_INT lda = dA.m;
-  DSQSS_INT ldb = dA.n;
-
-  //  cout<<"M, N, K ="<< M <<" "<<N<<" "<<K<<" "<<endl;
-
-  dgemm_(&transa, &transb, &M, &N, &K, &alpha, A, &lda, B, &ldb, &beta, C, &M);
-
-  for (int i = 0; i < dA.m; i++) {
-    for (int j = 0; j < dB.n; j++) {
-      dC(i, j) = C[i + j * dA.m];
-    }
-  }
-
-  delete[] A;
-  delete[] B;
-  delete[] C;
-
-  return dC;
-};
-
-//
-class dsymatrix {
- private:
-  double* index;
-
- public:
-  int n;
-
-  inline double& operator()(const int& i, const int& j) {
-    return index[i + j * n];
-  };
-  inline double operator()(const int& i, const int& j) const {
-    return index[i + j * n];
-  };
-
-  DSQSS_INT dsyev(vector<double>&);
-
-  dsymatrix(int N) {
-    n = N;
-    index = new double[N * N];
-  };
-
-  dsymatrix() { n = 0; };
-
-  ~dsymatrix() { delete[] index; };
-};
-
-DSQSS_INT dsymatrix::dsyev(vector<double>& E) {
-  ////////////////////
-  DSQSS_INT info;
-  char jobz = 'V';
-  char uplo = 'I';
-  DSQSS_INT w_size = -1;
-  DSQSS_INT size = n;
-  DSQSS_INT size2 = size * size;
-  double* work;
-  work = new double[size * 20];
-  double* vr;
-  vr = new double[size2];
-  double* wr;
-  wr = new double[size];
-  //  wr //eigen value
-  //  vr //eigen vector
-
-  for (int i = 0; i < (int)size2; i++) vr[i] = index[i];
-
-  dsyev_(&jobz, &uplo, &size, vr, &size, wr, work, &w_size, &info);
-
-  for (int i = 0; i < (int)size; i++) E[i] = wr[i];
-  for (int i = 0; i < (int)size2; i++) index[i] = vr[i];
-  ///////////////
-
-  delete[] work;
-  delete[] vr;
-  delete[] wr;
-
-  return info;
-};
-
-//============================================================================
-
-complex<double> IUNIT(0.0, 1.0);
-
-//============================================================================
-//    Display
-//============================================================================
-
-void dump(const vector<double>& V) {
-  printf("\n");
-  for (int i = 0; i < V.size(); i++) {
-    printf(" %8.3f", V[i]);
-  }
-  printf("\n");
-}
-
-//----------------------------------------------------------------------------
-
-void dump(char* s, const vector<double>& V) {
-  printf("\n");
-  printf("%s\n", s);
-  dump(V);
-}
-
-//----------------------------------------------------------------------------
-
-void dump(const dgematrix& A, int Mmax = 10) {
-  int M = A.m;
-  if (M > Mmax) M = Mmax;
-  int N = A.n;
-  if (N > Mmax) N = Mmax;
-  printf("\n");
-  for (int i = 0; i < M; i++) {
-    for (int j = 0; j < N; j++) {
-      printf(" %8.3f", A(i, j));
-      //    printf(" %2d", (int)(A(i,j)+0.1));
-    }
-    printf("\n");
-  }
-}
-
-//----------------------------------------------------------------------------
-
-void dump01(const dgematrix& A, int Mmax = 64) {
-  int M = A.m;
-  if (M > Mmax) M = Mmax;
-  int N = A.n;
-  if (N > Mmax) N = Mmax;
-  printf("\n");
-  for (int i = 0; i < M; i++) {
-    for (int j = 0; j < N; j++) {
-      char x = '.';
-      if (abs(A(i, j)) > 1.0e-8) x = 'X';
-      printf("%1c", x);
-    }
-    printf("\n");
-  }
-}
-
-//----------------------------------------------------------------------------
-
-void dump(char* s, const dgematrix& A) {
-  printf("\n");
-  printf("%s\n", s);
-  dump(A);
-}
-
-//============================================================================
-//    Diagonalization
-//============================================================================
-
-void diagonalize(dsymatrix& A, vector<double>& E, dgematrix& U) {
-  dsymatrix A0 = A;
-
-  A0.dsyev(E);
-
-  for (int i = 0; i < A.n; i++) {
-    for (int j = 0; j < A.n; j++) {
-      U(i, j) = A0(i, j);
-    }
-  }
-
-  return;
-}
-
-//----------------------------------------------------------------------------
-
-void diagonalize(dgematrix& A, vector<double>& E, dgematrix& U) {
-  dsymatrix A0(A.n);
-  for (int i = 0; i < A.n; i++) {
-    for (int j = 0; j <= i; j++) {
-      A0(i, j) = A(i, j);
-    }
-  }
-  diagonalize(A0, E, U);
-}
-
-//============================================================================
-//    Tensor Product
-//============================================================================
-
-dgematrix operator^(const dgematrix& A, const dgematrix& B) {
-  int _m = A.m * B.m;
-  int _n = A.n * B.n;
-  dgematrix C(_m, _n);
-
-  for (int i0 = 0; i0 < A.m; i0++) {
-    for (int i1 = 0; i1 < B.m; i1++) {
-      int i = i0 + A.m * i1;
-
-      for (int j0 = 0; j0 < A.n; j0++) {
-        for (int j1 = 0; j1 < B.n; j1++) {
-          int j = j0 + A.n * j1;
-          C(i, j) = A(i0, j0) * B(i1, j1);
-        }
-      }
-    }
-  }
-  return C;
-}
-
-//============================================================================
-//    Complex Matrix
-//============================================================================
-
-class cmatrix {
- public:
-  long m;
-  long n;
-
-  dgematrix re;
-  dgematrix im;
-
-  void resize(const long m0, const long n0) {
-    m = m0;
-    n = n0;
-    re.resize(m, n);
-    im.resize(m, n);
-  };
-
-  void resize(const long n0) { resize(n0, n0); }
-
-  cmatrix(){};
-  cmatrix(const long n0) { resize(n0, n0); };
-  cmatrix(const long m0, const long n0) { resize(m0, n0); };
-  cmatrix(const cmatrix& X) { resize(X.m, X.n); };
-
-  cmatrix& operator+=(const cmatrix& A) {
-    if (m != A.m) {
-      printf("cmatrix::operator+= >> ERROR;");
-      exit(0);
-    }
-    if (n != A.n) {
-      printf("cmatrix::operator+= >> ERROR;");
-      exit(0);
-    }
-    re += A.re;
-    im += A.im;
-    return *this;
-  };
-
-  cmatrix& operator-=(const cmatrix& A) {
-    if (m != A.m) {
-      printf("cmatrix::operator+= >> ERROR;");
-      exit(0);
-    }
-    if (n != A.n) {
-      printf("cmatrix::operator+= >> ERROR;");
-      exit(0);
-    }
-    re -= A.re;
-    im -= A.im;
-    return *this;
-  };
-
-  cmatrix& operator=(const cmatrix& A) {
-    m = A.m;
-    n = A.n;
-    re = A.re;
-    im = A.im;
-    return *this;
-  };
-
-  void clear() {
-    re.clear();
-    im.clear();
-  };
-
-  void zero();
-
-  void unity() {
-    resize(1, 1);
-    re(0, 0) = 1.0;
-    im(0, 0) = 0.0;
-  };
-
-  void identity() {
-    re.identity();
-    im.identity();
-  }
-
-  void dump(int Mmax);
-  void dump(char*);
-  void dump01(int Mmax);
-};
-
-//----------------------------------------------------------------------------
-
-void cmatrix::zero() {
-  for (int i = 0; i < m; i++) {
-    for (int j = 0; j < n; j++) {
-      re(i, j) = 0.0;
-      im(i, j) = 0.0;
-    }
-  }
-}
-
-//----------------------------------------------------------------------------
-
-void cmatrix::dump(int Mmax = 10) {
-  int M = m;
-  if (M > Mmax) M = Mmax;
-  int N = n;
-  if (N > Mmax) N = Mmax;
-  printf("\n");
-  for (int i = 0; i < M; i++) {
-    for (int j = 0; j < N; j++) {
-      //    if ( j!= 0) printf(" |");
-      printf(" %6.3f", re(i, j));
-      //    printf(" %4.1f %4.1f", re(i,j), im(i,j));
-    }
-    printf("\n");
-  }
-}
-
-//----------------------------------------------------------------------------
-
-void cmatrix::dump01(int Mmax = 64) {
-  int M = m;
-  if (M > Mmax) M = Mmax;
-  int N = n;
-  if (N > Mmax) N = Mmax;
-  printf("\n");
-  for (int i = 0; i < M; i++) {
-    for (int j = 0; j < N; j++) {
-      char x = '.';
-      if (re(i, j) * re(i, j) + im(i, j) * im(i, j) > 1.0e-16) x = 'X';
-      printf("%1c", x);
-    }
-    printf("\n");
-  }
-}
-
-//----------------------------------------------------------------------------
-
-void cmatrix::dump(char* s) {
-  printf("%s\n", s);
-  dump();
-}
-
-//----------------------------------------------------------------------------
-
-cmatrix operator+(const cmatrix& A, const cmatrix& B) {
-  cmatrix C(A.m, B.n);
-  C.re = A.re + B.re;
-  C.im = A.im + B.im;
-  return C;
-}
-
-//----------------------------------------------------------------------------
-
-cmatrix operator-(const cmatrix& A, const cmatrix& B) {
-  cmatrix C(A.m, B.n);
-  C.re = A.re - B.re;
-  C.im = A.im - B.im;
-  return C;
-}
-
-//----------------------------------------------------------------------------
-
-cmatrix operator*(const cmatrix& A, const cmatrix& B) {
-  cmatrix C(A.m, B.n);
-  C.re = A.re * B.re - A.im * B.im;
-  C.im = A.re * B.im + A.im * B.re;
-  return C;
-}
-
-//----------------------------------------------------------------------------
-
-cmatrix operator*(const double a, const cmatrix& A) {
-  cmatrix C(A.m, A.n);
-  C.re = a * A.re;
-  C.im = a * A.im;
-  return C;
-}
-
-//----------------------------------------------------------------------------
-
-cmatrix operator*(const complex<double> c, const cmatrix& A) {
-  cmatrix C(A.m, A.n);
-  C.re = c.real() * A.re - c.imag() * A.im;
-  C.im = c.real() * A.im + c.imag() * A.re;
-  return C;
-}
-
-//----------------------------------------------------------------------------
-
-cmatrix t(const cmatrix& A) {
-  cmatrix C(A.n, A.m);
-  C.re = t(A.re);
-  C.im = t(A.im);
-  return C;
-}
-
-//----------------------------------------------------------------------------
-
-cmatrix operator^(cmatrix& A, cmatrix& B) {
-  int m = A.m * B.m;
-  int n = A.n * B.n;
-  cmatrix C(m, n);
-  C.re = ((A.re) ^ (B.re)) - ((A.im) ^ (B.im));
-  C.im = ((A.re) ^ (B.im)) + ((A.im) ^ (B.re));
-  return C;
-}
diff --git a/src/dla/generators/sfgene.cc b/src/dla/generators/sfgene.cc
deleted file mode 100644
index 9657925f..00000000
--- a/src/dla/generators/sfgene.cc
+++ /dev/null
@@ -1,194 +0,0 @@
-/*---------------------------------------------
-
-   Generating cf.xml for a hypercubic lattice
-
-----------------------------------------------*/
-
-#include <cmath>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-#include <fstream>
-#include <iostream>
-#include <sstream>
-#include <stdexcept>
-#include <string>
-#include <vector>
-
-#include "../util.hpp"
-
-using namespace std;
-
-std::vector<int> index2coord(int i, std::vector<int> L) {
-  if (i < 0) {
-    std::stringstream ss;
-    ss << "invalid index i=";
-    ss << i;
-    util::ERROR(ss.str().c_str());
-  }
-  const int D = L.size();
-  int N = 1;
-  std::vector<int> r(D);
-  for (int d = 0; d < D; ++d) {
-    r[d] = i % L[d];
-    i /= L[d];
-    N *= L[d];
-  }
-  if (i >= N) {
-    std::stringstream ss;
-    ss << "ERROR: invalid index i=";
-    ss << i;
-    util::ERROR(ss.str().c_str());
-  }
-  return r;
-}
-
-void ShowUsage(std::string const& exename) {
-  cout << "usage:\n";
-  cout << "    " << exename
-       << " [-o filename] D L_1 ... L_D Ntau Ntau_cutoff KTYPE\n";
-  cout << "arguments:\n";
-  cout << "    D             ... dimension of lattice\n";
-  cout << "    L_d           ... the liner size of the lattice\n";
-  cout << "                      (must be even)\n";
-  cout << "    Ntau          ... number of discretized imaginary time\n";
-  cout << "    Ntau_cutoff   ... maximum distance in imaginary time between "
-          "two spacetime points\n";
-  cout << "    KTYPE         ... type of wavenumbers\n";
-  cout << "                      0: kx=pi n/L, n=0, 2, ..., L" << endl;
-  cout << "                      1: k/pi = (0,0), (1,0), (0,1), (1,1) [example "
-          "in D=2 case]"
-       << endl;
-  cout << "options:\n";
-  cout << "    -o filename   ... output file (default: sf.xml)";
-}
-
-void WriteXML(std::vector<int> const& L, int Ntau, int CutoffOfNtau, int KTYPE,
-              std::string const& filename) {
-  const int D = L.size();
-  ofstream fout(filename.c_str());
-  fout.precision(15);
-  int N = 1;  // number of sites.
-  for (int i = 0; i < D; i++) {
-    N *= L[i];
-  }
-  int KMAX;
-
-  std::vector<int> Q(D, 2);
-
-  if (KTYPE == 0) {
-    KMAX = L[0] / 2 + 1;
-  } else if (KTYPE == 1) {
-    KMAX = 1;
-    for (int d = 0; d < D; ++d) {
-      KMAX *= 2;
-    }
-  } else {
-    cout << " KTYPE should 0 or 1";
-    exit(1);
-  }
-
-  int NumberOfElements = N * KMAX;
-
-  fout << "<StructureFactor>" << endl << endl;
-  fout << "<Comment>" << endl;
-  fout << "  " << D << "-dimension hypercubic lattice" << endl;
-  fout << "</Comment>" << endl << endl;
-
-  fout << "<Ntau>                   " << Ntau << " </Ntau>" << endl;
-  fout << "<NumberOfElements>       " << NumberOfElements
-       << " </NumberOfElements>" << endl;
-  fout << "<CutoffOfNtau>           " << CutoffOfNtau << " </CutoffOfNtau>"
-       << endl;
-  fout << "<NumberOfInverseLattice> " << KMAX << " </NumberOfInverseLattice>"
-       << endl;
-
-  fout << endl;
-  fout << "<!-- <SF> [phase(cos)] [phase(sin)] [isite] [kindx] </SF> -->"
-       << endl
-       << endl;
-
-  int NB = 0;  // 3 * N ;   // number of bonds
-
-  for (int i = 0; i < N; i++) {
-    std::vector<int> r = index2coord(i, L);
-    for (int q = 0; q < KMAX; q++) {
-      double phase = 0;
-      if (KTYPE == 0) {
-        phase = (r[0] * 2 * q * M_PI) / L[0];  // r_x * q_x
-      } else if (KTYPE == 1) {
-        std::vector<int> k = index2coord(q, Q);
-        for (int d = 0; d < D; ++d) {
-          phase += r[d] * k[d];
-        }
-        phase *= M_PI;
-      }
-      double COSrk = cos(phase);
-      double SINrk = sin(phase);
-      if (fabs(COSrk) < 1.0e-12) COSrk = 0.0;
-      if (fabs(SINrk) < 1.0e-12) SINrk = 0.0;
-
-      fout << "<SF> " << COSrk << " " << SINrk << " " << i << " " << q
-           << " </SF>" << endl;
-    }
-  }
-
-  fout << endl;
-  fout << "</StructureFactor>" << endl;
-}
-
-int main(int argc, char** argv) {
-  std::string exename(argv[0]);
-  std::string filename("sf.xml");
-  if (argc < 3) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  if (std::strcmp(argv[1], "-o") == 0) {
-    filename = argv[2];
-    argc -= 2;
-    argv += 2;
-  }
-  int NARG = 3;
-  if (argc < NARG + 1) {
-    ShowUsage(exename);
-    exit(0);
-  }
-  int iarg = 1;
-  const int D = atoi(argv[iarg]);
-  iarg++;
-  if (argc != D + 5) {
-    ShowUsage(exename);
-    exit(0);
-  }
-
-  std::vector<int> L(D);
-
-  for (int i = 0; i < D; i++) {
-    L[i] = atoi(argv[iarg]);
-    iarg++;
-  }
-  int Ntau = atoi(argv[iarg]);
-  iarg++;
-
-  int Ntcut = atoi(argv[iarg]);
-  iarg++;
-
-  int KTYPE = atoi(argv[iarg]);
-
-  int EvenOrOdd = 0;
-  cout.precision(15);
-  cout << "D     = " << D << endl;
-  for (int i = 0; i < D; i++) {
-    cout << "L_" << i + 1 << "    = " << L[i] << endl;
-    EvenOrOdd += L[i] % 2;
-  }
-
-  if (EvenOrOdd) {
-    cout << "Warnig: L should be an even number." << endl;
-  }
-
-  WriteXML(L, Ntau, Ntcut, KTYPE, filename);
-  cout << "... done." << endl;
-  return 0;
-}
diff --git a/src/dla/generators/spin_H.h b/src/dla/generators/spin_H.h
deleted file mode 100644
index 38dbb18e..00000000
--- a/src/dla/generators/spin_H.h
+++ /dev/null
@@ -1,113 +0,0 @@
-
-//============================================================================
-//    Spin Matrices
-//============================================================================
-
-class HeisenbergSpin {
- public:
-  int K;
-  int D;
-  cmatrix I;
-  cmatrix UP;
-  cmatrix DN;
-  cmatrix X;
-  cmatrix Y;
-  cmatrix Z;
-  HeisenbergSpin(int);
-  ~HeisenbergSpin(){};
-};
-
-//----------------------------------------------------------------------------
-
-HeisenbergSpin::HeisenbergSpin(int K0) {
-  K = K0;  // K = 2 * S
-  D = K0 + 1;
-  I.resize(D, D);
-  UP.resize(D, D);
-  DN.resize(D, D);
-  X.resize(D, D);
-  Y.resize(D, D);
-  Z.resize(D, D);
-  I.zero();
-  UP.zero();
-  DN.zero();
-  X.zero();
-  Y.zero();
-  Z.zero();
-  for (int i = 0; i < D - 1; i++) {
-    UP.re(i + 1, i) = sqrt((double)((i + 1) * (K - i)));
-  }
-  DN = t(UP);
-  //  X = 0.5 * (UP + DN);
-  cmatrix temp = UP + DN;
-  X = 0.5 * temp;
-  Y = ((-0.5) * IUNIT) * (UP - DN);
-  for (int i = 0; i < D; i++) {
-    I.re(i, i) = 1.0;
-    Z.re(i, i) = -0.5 * (double)K + (double)i;
-  }
-}
-
-//----------------------------------------------------------------------------
-
-class HeisenbergSpinSet {
- public:
-  int DS;     // the dimension of the 1-spin Hilbert space
-  int N;      // "N" in SU(N)
-  int K;      // the number of bosons on each site
-  int NSITE;  // the number of spins
-  int DIM;    // the dimension of the whole Hilbert space
-
-  cmatrix* X;
-  cmatrix* Y;
-  cmatrix* Z;
-  cmatrix I;
-
-  HeisenbergSpinSet(int K0, int NSITE0) {
-    printf("SpinSet> start.\n");
-
-    K = K0;
-    NSITE = NSITE0;
-
-    HeisenbergSpin S(K);
-
-    DS = S.D;
-    X = new cmatrix[NSITE];
-    Y = new cmatrix[NSITE];
-    Z = new cmatrix[NSITE];
-
-    printf("  definig spins ...\n");
-    DIM = 1;
-
-    for (int i = 0; i < NSITE; i++) {
-      DIM *= (DS);
-      cmatrix x;
-      cmatrix y;
-      cmatrix z;
-      x.unity();
-      y.unity();
-      z.unity();
-      for (int j = 0; j < NSITE; j++) {
-        if (i == j) {
-          x = (S.X) ^ x;
-          y = (S.Y) ^ y;
-          z = (S.Z) ^ z;
-        } else {
-          x = (S.I) ^ x;
-          y = (S.I) ^ y;
-          z = (S.I) ^ z;
-        }
-
-        X[i] = x;
-        Y[i] = y;
-        Z[i] = z;
-      }
-      printf("  ... done.\n");
-
-      I.resize(DIM);
-      I.identity();
-
-      printf("SpinSet> end.\n");
-    };
-  };
-};