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& 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 ( - ) -// --> -^2 (in the classical or the high-T limit) -// -double Susceptibility(double B, vector& 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 -#include -#include -#include -#include -#include - -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 << "" << endl << endl; - fout << "" << endl; - fout << " " << D << "-dimension hypercubic lattice" << endl; - fout << "" << endl << endl; - - fout << " " << Ntau << " " << endl; - fout << " " << NumberOfElements << " " - << endl; - fout << " " << NumberOfKinds << " " << endl; - fout << endl; - - fout << "" << 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 << " " << kind << " " << i << " " << j << " " << endl; - } - kind++; - } - cout << "NumberOfKinds = " << kind << endl; - fout << endl; - fout << "" << 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, "\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\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, " \n"); - fprintf(FALG, " %s\n", comment.c_str()); - fprintf(FALG, " \n"); - fprintf(FALG, "\n"); - fprintf(FALG, " \n"); - fprintf(FALG, " %2d \n", NSTYPE); - fprintf(FALG, " %2d \n", NITYPE); - fprintf(FALG, " %2d \n", NVTYPE); - fprintf(FALG, " %2d \n", NXMAX); - fprintf(FALG, " %24.16lf \n", WeightDiagonal); - fprintf(FALG, " \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, " \n"); - fprintf(FALG, " %d \n", ID); - fprintf(FALG, " %d \n", NX); - fprintf(FALG, " %d \n", VTYPE); - // printf(" %24.16f \n", V().EBASE ); - // fprintf(FALG," %24.16f \n", V().EBASE ); - for (int x0 = 0; x0 < NX; x0++) { - fprintf(FALG, " \n"); - fprintf(FALG, " %d \n", x0); - int NCH = NumberOfChannels[x0]; - fprintf(FALG, " %d \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, " %4d %4d %24.16lf \n", diri, xnew, - p); - } - fprintf(FALG, " \n"); - } - fprintf(FALG, " \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, " \n"); - fprintf(FALG, " %d \n", ID); - fprintf(FALG, " %d \n", VTYPE); - fprintf(FALG, " %d \n", NBODY); - fprintf(FALG, " %24.16lf \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, " "); - for (int j = 0; j < NBODY; j++) fprintf(FALG, " %2d", x[j]); - fprintf(FALG, " %24.16lf \n", d); - } - } - delete[] x; - } else { - // fprintf(FALG," VertexDensity is not defined.\n"); - } - - int NLEG = 2 * NBODY; - Array& 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, " "); - for (int j = 0; j < 2 * NBODY; ++j) { - fprintf(FALG, " %2d", x[j]); - } - fprintf(FALG, " %24.16lf \n", sgn); - } - } - } - delete[] x; - - fprintf(FALG, " \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, " \n"); - fprintf(FALG, " %d \n", ID); - fprintf(FALG, " %d \n", CATEGORY); - fprintf(FALG, " %d \n", NBODY); - int NICW = NumberOfValidInitialConfigurations(); - fprintf(FALG, - " %d " - "\n", - NICW); - for (int i = 0; i < NICG; i++) ICG(i).write(); - fprintf(FALG, "\n \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 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 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 PERM("PERM"); - PERM.init(1, NIC); - bubble_sort(NIC, U, PERM); - - Array UP("UP"); - UP.init(1, NIC); - for (int i = 0; i < NIC; i++) { - UP[i] = U[PERM[i]]; - } - - Array P("P"); - Array W("W"); - Array 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 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 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, " \n"); - fprintf(FALG, " "); - for (int i = 0; i < NLEG; i++) fprintf(FALG, " %d", State[i]); - fprintf(FALG, " \n"); - fprintf(FALG, " %d \n", INC); - fprintf(FALG, " %d \n", XINC); - fprintf(FALG, " %d \n", NCH); - for (int ch = 0; ch < NCH; ch++) { - int out = FinalDirection[ch]; - int xout = FinalState[ch]; - double p = ScatteringProbability[ch]; - fprintf(FALG, " %4d %4d %24.16lf \n", out, xout, - p); - } - fprintf(FALG, " \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, " \n"); - fprintf(FALG, " %d \n", getID()); - fprintf(FALG, " %s \n", getName().c_str()); - for (int st = 0; st < NSTYPE; st++) { - for (int x = 0; x < NXMAX; x++) { - if (isDefined(st, x)) { - fprintf(FALG, " %d %d %24.16lf \n", st, x, - Value(st, x)); - } - } - } - fprintf(FALG, " \n"); -} - -//###################################################################### - -//###################################################################### - -void SolveWeightEquation(int N, Array& V, Array& W) { - /* - for (int i=0; i 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& V, Array& 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 -#include - -#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& V, Array& W); -bool testKink(int ist); -bool testForbidden(int ist); -void bubble_sort(int N, Array& V, Array& 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 WormCreationNewState; - Array WormCreationDirection; - Array 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 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 Weight; - Array SiteTypeOfLeg; - Array 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 VertexDensity; - Array 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 U; - // Array W; - // Array P; - VERTEX* _V; - Array 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 Value; - Array 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 -#include -#include -#include -using namespace std; -#include "boson_B.h" -#include "canonical.h" -#include "matrix.h" - -//---------------------------------------------------------------------- - -// Bose-Hubbard Model -// H= - J\sum_{} b_i b_j + V\sum_{} n_i n_j -// + u/2z \sum_{}( n_i (n_i + 1) + n_j (n_j + 1) ) -// - mu/z \sum_{}( 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, "\n"); - fprintf(FOUT, " \n"); - fprintf(FOUT, - " Extended Bose-Hubbard model with NMAX=%d \n", - M); - fprintf(FOUT, " 1 \n"); - fprintf(FOUT, " 1 \n"); - fprintf(FOUT, " %d \n", D); - fprintf(FOUT, " \n"); - fprintf(FOUT, "\n"); - fprintf(FOUT, " \n"); - fprintf(FOUT, " 0 \n"); - fprintf(FOUT, " 0 \n"); - fprintf(FOUT, " %d \n", D); - fprintf(FOUT, " \n"); - fprintf(FOUT, "\n"); - fprintf(FOUT, " \n"); - fprintf(FOUT, " 0 \n"); - fprintf(FOUT, " 0 \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, " %d %d %24.16f \n", i, j, x); - } - } - } - fprintf(FOUT, " \n"); - fprintf(FOUT, "\n"); - fprintf(FOUT, " \n"); - fprintf(FOUT, " 0 \n"); - fprintf(FOUT, " 2 \n"); - fprintf(FOUT, " 0 0 \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, " %d %d %d %d %24.16f \n", i0, j0, - i1, j1, x); - } - } - } - fprintf(FOUT, " \n"); - fprintf(FOUT, "\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 -#include -#include -#include -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, "\n"); - fprintf(FOUT, " \n"); - fprintf(FOUT, " SU(2) Heisenberg model with S=%d/2 \n", - M); - fprintf(FOUT, " 1 \n"); - fprintf(FOUT, " 1 \n"); - fprintf(FOUT, " %d \n", D); - fprintf(FOUT, " \n"); - fprintf(FOUT, "\n"); - fprintf(FOUT, " \n"); - fprintf(FOUT, " 0 \n"); - fprintf(FOUT, " 0 \n"); - fprintf(FOUT, " %d \n", D); - fprintf(FOUT, " \n"); - fprintf(FOUT, "\n"); - fprintf(FOUT, " \n"); - fprintf(FOUT, " 0 \n"); - fprintf(FOUT, " 0 \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, " %d %d %24.16f \n", i, j, x); - } - } - } - fprintf(FOUT, " \n"); - fprintf(FOUT, "\n"); - fprintf(FOUT, " \n"); - fprintf(FOUT, " 0 \n"); - fprintf(FOUT, " 2 \n"); - fprintf(FOUT, " 0 0 \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, " %d %d %d %d %24.16f \n", i0, j0, - i1, j1, x); - } - } - } - fprintf(FOUT, " \n"); - fprintf(FOUT, "\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 -#include -#include -#include -#include -#include - -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 << "" << endl << endl; - fout << "" << endl; - fout << " " << D << "-dimension square lattice" << endl; - fout << "" << endl << endl; - - fout << " " << D << " " << endl; - fout << " " << BD << " " << endl; - fout << " "; - for (int i = 0; i < D; i++) { - fout << L[i] << " "; - } - fout << "" << endl; - fout << " " << N << " " << endl; - fout << " " << NumberOfInteractions - << " " << endl; - fout << " " << NumberOfSiteTypes << " " - << endl; - fout << " " << NumberOfInteractionTypes - << " " << endl; - fout << " " << NumberOfEdgeInteractions - << " " << endl; - - fout << endl; - - fout << "" << 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 << " " << id << " " << stype << " " << mtype << " " << endl; - } - - fout << endl; - fout << "" - << 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 << " " << NB << " " << itype << " " << nbody << " " << i << " " - << j << " " << etype << " " << p << " " << endl; - - NB++; - } - } - - fout << endl; - fout << "" << endl << endl; - for (int p = 0; p < BD; p++) { - fout << " " << p << " "; - for (int q = 0; q < D; q++) { - if (p == q) - fout << 1 << " "; - else - fout << 0 << " "; - } - fout << "" << endl; - } - fout << "" << 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 -#include -#include -#include -#include -#include -#include - -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 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 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 << "" << endl << endl; - fout << "" << endl; - fout << " " << D << "-dimension triangular lattice" << endl; - fout << "" << endl << endl; - - fout << " " << D << " " << endl; - fout << " " << BD << " " << endl; - fout << " "; - for (int i = 0; i < D; i++) { - fout << L[i] << " "; - } - fout << "" << endl; - fout << " " << N << " " << endl; - fout << " " << NumberOfInteractions - << " " << endl; - fout << " " << NumberOfSiteTypes << " " - << endl; - fout << " " << NumberOfInteractionTypes - << " " << endl; - fout << " " << NumberOfEdgeInteractions - << " " << endl; - - fout << endl; - - fout << "" << endl; - - std::vector 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 << " " << id << " " << stype << " " << mtype << " " << endl; - } - - fout << endl; - fout << "" - << 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 << " " << NB << " " << itype << " " << nbody << " " << i << " " << j - << " " << etype << " " << p << " " << 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 << " " << NB << " " << itype << " " << nbody << " " << i << " " << j - << " " << etype << " " << p << " " << 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 << " " << NB << " " << itype << " " << nbody << " " << i << " " << j - << " " << etype << " " << p << " " << endl; - NB++; - } - - fout << endl; - fout << "" << endl << endl; - - fout << " " << 0 << " " << 1 << " " << 0 << " " << endl; - fout << " " << 1 << " " << 0 << " " << 1 << " " << endl; - fout << " " << 2 << " " << 1 << " " << 1 << " " << endl; - - fout << endl; - fout << "" << 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 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."< - -#include -#include -#include -#include -#include -#include -#include -#include -#include //only for ?geev, ?gegv, etc. - -//#include - -//============================================================================= -#ifdef __INTEL_COMPILER -#define DSQSS_INT MKL_INT -#define MKL_Complex16 std::complex -#else //__INTEL_COMPILER -#define DSQSS_INT int -#endif //__INTEL_COMPILER - -//============================================================================= - -#ifdef __INTEL_COMPILER -#include -#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 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 <<" "<&); - - dsymatrix(int N) { - n = N; - index = new double[N * N]; - }; - - dsymatrix() { n = 0; }; - - ~dsymatrix() { delete[] index; }; -}; - -DSQSS_INT dsymatrix::dsyev(vector& 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 IUNIT(0.0, 1.0); - -//============================================================================ -// Display -//============================================================================ - -void dump(const vector& V) { - printf("\n"); - for (int i = 0; i < V.size(); i++) { - printf(" %8.3f", V[i]); - } - printf("\n"); -} - -//---------------------------------------------------------------------------- - -void dump(char* s, const vector& 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& 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& 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 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 -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "../util.hpp" - -using namespace std; - -std::vector index2coord(int i, std::vector 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 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 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 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 << "" << endl << endl; - fout << "" << endl; - fout << " " << D << "-dimension hypercubic lattice" << endl; - fout << "" << endl << endl; - - fout << " " << Ntau << " " << endl; - fout << " " << NumberOfElements - << " " << endl; - fout << " " << CutoffOfNtau << " " - << endl; - fout << " " << KMAX << " " - << endl; - - fout << endl; - fout << "" - << endl - << endl; - - int NB = 0; // 3 * N ; // number of bonds - - for (int i = 0; i < N; i++) { - std::vector 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 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 << " " << COSrk << " " << SINrk << " " << i << " " << q - << " " << endl; - } - } - - fout << endl; - fout << "" << 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 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"); - }; - }; -};