From 2c08dae1f41fc292756c157f0419c0f923722a47 Mon Sep 17 00:00:00 2001 From: Christoph Niethammer Date: Thu, 19 Dec 2024 15:43:52 +0100 Subject: [PATCH] Remove standalone-generators The standalone generators are outdated and cannot produce input files that are compatible with the current ls1 version. Use the internal generators in ls1 instead. Signed-off-by: Christoph Niethammer --- tools/standalone-generators/CMakeLists.txt | 21 - tools/standalone-generators/README | 5 - .../animake/CMakeLists.txt | 9 - .../standalone-generators/animake/Domain.cpp | 1090 ------------- tools/standalone-generators/animake/Domain.h | 184 --- .../standalone-generators/animake/Random.cpp | 30 - tools/standalone-generators/animake/Random.h | 14 - tools/standalone-generators/animake/main.cpp | 265 ---- tools/standalone-generators/atomic_units.txt | 21 - tools/standalone-generators/clean | 3 - tools/standalone-generators/install | 2 - .../standalone-generators/mkcp/CMakeLists.txt | 9 - tools/standalone-generators/mkcp/Domain.cpp | 1356 ----------------- tools/standalone-generators/mkcp/Domain.h | 179 --- tools/standalone-generators/mkcp/Graphit.cpp | 237 --- tools/standalone-generators/mkcp/Graphit.h | 60 - tools/standalone-generators/mkcp/Random.cpp | 30 - tools/standalone-generators/mkcp/Random.h | 14 - tools/standalone-generators/mkcp/main.cpp | 594 -------- .../mkesfera/CMakeLists.txt | 9 - .../standalone-generators/mkesfera/Domain.cpp | 197 --- tools/standalone-generators/mkesfera/Domain.h | 25 - .../standalone-generators/mkesfera/Random.cpp | 30 - tools/standalone-generators/mkesfera/Random.h | 15 - tools/standalone-generators/mkesfera/main.cpp | 106 -- .../mktcts/CMakeLists.txt | 9 - tools/standalone-generators/mktcts/Domain.cpp | 304 ---- tools/standalone-generators/mktcts/Domain.h | 36 - tools/standalone-generators/mktcts/Random.cpp | 30 - tools/standalone-generators/mktcts/Random.h | 15 - tools/standalone-generators/mktcts/main.cpp | 158 -- tools/standalone-generators/produce | 2 - tools/standalone-generators/uninstall | 2 - 33 files changed, 5061 deletions(-) delete mode 100644 tools/standalone-generators/CMakeLists.txt delete mode 100644 tools/standalone-generators/README delete mode 100644 tools/standalone-generators/animake/CMakeLists.txt delete mode 100644 tools/standalone-generators/animake/Domain.cpp delete mode 100644 tools/standalone-generators/animake/Domain.h delete mode 100644 tools/standalone-generators/animake/Random.cpp delete mode 100644 tools/standalone-generators/animake/Random.h delete mode 100644 tools/standalone-generators/animake/main.cpp delete mode 100644 tools/standalone-generators/atomic_units.txt delete mode 100755 tools/standalone-generators/clean delete mode 100755 tools/standalone-generators/install delete mode 100644 tools/standalone-generators/mkcp/CMakeLists.txt delete mode 100644 tools/standalone-generators/mkcp/Domain.cpp delete mode 100644 tools/standalone-generators/mkcp/Domain.h delete mode 100644 tools/standalone-generators/mkcp/Graphit.cpp delete mode 100644 tools/standalone-generators/mkcp/Graphit.h delete mode 100644 tools/standalone-generators/mkcp/Random.cpp delete mode 100644 tools/standalone-generators/mkcp/Random.h delete mode 100644 tools/standalone-generators/mkcp/main.cpp delete mode 100644 tools/standalone-generators/mkesfera/CMakeLists.txt delete mode 100644 tools/standalone-generators/mkesfera/Domain.cpp delete mode 100644 tools/standalone-generators/mkesfera/Domain.h delete mode 100644 tools/standalone-generators/mkesfera/Random.cpp delete mode 100644 tools/standalone-generators/mkesfera/Random.h delete mode 100644 tools/standalone-generators/mkesfera/main.cpp delete mode 100644 tools/standalone-generators/mktcts/CMakeLists.txt delete mode 100644 tools/standalone-generators/mktcts/Domain.cpp delete mode 100644 tools/standalone-generators/mktcts/Domain.h delete mode 100644 tools/standalone-generators/mktcts/Random.cpp delete mode 100644 tools/standalone-generators/mktcts/Random.h delete mode 100644 tools/standalone-generators/mktcts/main.cpp delete mode 100755 tools/standalone-generators/produce delete mode 100755 tools/standalone-generators/uninstall diff --git a/tools/standalone-generators/CMakeLists.txt b/tools/standalone-generators/CMakeLists.txt deleted file mode 100644 index 578ffd48b7..0000000000 --- a/tools/standalone-generators/CMakeLists.txt +++ /dev/null @@ -1,21 +0,0 @@ -cmake_minimum_required(VERSION 3.3) - -project(ls1-standalone-generators LANGUAGES CXX VERSION 1.2.1) -set(AUTHOR "Martin Thomas Horsch") -set(AUTHOR_DETAILS "martin.horsch@stfc.ac.uk") -message(STATUS "cmake building ${PROJECT_NAME}") - -add_subdirectory(animake) -add_subdirectory(mkcp) -add_subdirectory(mkesfera) -add_subdirectory(mktcts) - -set(CPACK_GENERATOR "DEB") -set(CPACK_PACKAGE_DESCRIPTION_SUMMARY "standalone external scenario generators for the MD code ls1 mardyn") -set(CPACK_PACKAGE_VERSION_MAJOR "1") -set(CPACK_PACKAGE_VERSION_MINOR "2") -set(CPACK_PACKAGE_VERSION_PATCH "1") -set(CPACK_PACKAGE_INSTALL_DIRECTORY "/usr/bin") -set(CPACK_DEBIAN_PACKAGE_MAINTAINER "Martin Thomas Horsch") -set(CPACK_PACKAGE_EXECUTABLES mktcts) -include(CPack) diff --git a/tools/standalone-generators/README b/tools/standalone-generators/README deleted file mode 100644 index 353ed71624..0000000000 --- a/tools/standalone-generators/README +++ /dev/null @@ -1,5 +0,0 @@ -auto-installation for Debian-like systems with cmake: - -./install -./clean - diff --git a/tools/standalone-generators/animake/CMakeLists.txt b/tools/standalone-generators/animake/CMakeLists.txt deleted file mode 100644 index a755e7476b..0000000000 --- a/tools/standalone-generators/animake/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -cmake_minimum_required(VERSION 3.3) -message(STATUS "cmake building animake") - -add_executable(animake main.cpp Domain.cpp Random.cpp) -target_include_directories(animake PRIVATE ${PROJECT_SOURCE_DIR}) - -include(GNUInstallDirs) -install(TARGETS animake RUNTIME DESTINATION /usr/bin/) - diff --git a/tools/standalone-generators/animake/Domain.cpp b/tools/standalone-generators/animake/Domain.cpp deleted file mode 100644 index 1f999e8729..0000000000 --- a/tools/standalone-generators/animake/Domain.cpp +++ /dev/null @@ -1,1090 +0,0 @@ -#include "Domain.h" -#include "Random.h" -#include - -#define DT 0.061240 -#define PRECISION 8 -#define BINS 1024 - -Domain::Domain( - int sp_fluid, int sp_fluid2, double* sp_box, double sp_SIG_REF, - double sp_EPS_REF, double sp_REFMASS, bool sp_muVT, - unsigned sp_N, double sp_T, double sp_ETAF, double sp_XIF -) { - this->fluid = sp_fluid; - this->fluid2 = sp_fluid2; - for(int d=0; d < 3; d++) this->box[d] = sp_box[d]; - this->SIG_REF = sp_SIG_REF; - this->EPS_REF = sp_EPS_REF; - this->REFMASS = sp_REFMASS; - this->muVT = sp_muVT; - this->N = sp_N; - this->T = sp_T; - this->ETAF = sp_ETAF; - this->XIF = sp_XIF; -} - -void Domain::write(char* prefix, int format, double mu, double x) -{ - std::ofstream xdr, txt, buchholz; - std::stringstream strstrm, txtstrstrm, buchholzstrstrm; - if(format == FORMAT_BRANCH) - { - strstrm << prefix << ".xdr"; - } - if((format == FORMAT_BERNREUTHER) || (format == FORMAT_BUCHHOLZ)) - { - strstrm << prefix << ".inp"; - } - xdr.open(strstrm.str().c_str(), std::ios::trunc); - if(format == FORMAT_BRANCH) - { - txtstrstrm << prefix << "_1R.txt"; - } - if(format == FORMAT_BUCHHOLZ) - { - txtstrstrm << prefix << "_1R.cfg"; - } - if(format == FORMAT_BERNREUTHER) - { - txtstrstrm << prefix << "_1R.xml"; - } - txt.open(txtstrstrm.str().c_str(), std::ios::trunc); - if(format == FORMAT_BERNREUTHER) - { - txt << "\n\n"; - } - if(format == FORMAT_BUCHHOLZ) - { - buchholzstrstrm << prefix << "_1R.xml"; - buchholz.open(buchholzstrstrm.str().c_str(), std::ios::trunc); - - /* - * Gesamter Inhalt der Buchholz-Datei - */ - buchholz << "\n\n \n " - << prefix << "_1R.cfg\n \n"; - buchholz.close(); - } - - Random* r = new Random(); - r->init( - (int)(3162.3*box[0]) + (int)(31623.0*box[1]) - (int)(316.23*box[2]) - ); - double REFTIME = SIG_REF * sqrt(REFMASS / EPS_REF); - double VEL_REF = SIG_REF / REFTIME; - std::cout << "Velocity unit 1 = " << VEL_REF << " * 1620.34 m/s = " - << 1620.34 * VEL_REF << " m/s.\n"; - double REFCARG = sqrt(EPS_REF * SIG_REF); - std::cout << "Charge unit 1 = " << REFCARG << " e.\n"; - double DIP_REF = SIG_REF*REFCARG; - double QDR_REF = SIG_REF*DIP_REF; - double REFOMGA = 1.0 / REFTIME; - if(format == FORMAT_BERNREUTHER) - { - txt << "" << SIG_REF * 0.0529177 << "" << REFMASS * 1000.0 << "" << EPS_REF * 27.2126 << "\n"; - } - - unsigned fl_units[3]; - double fl_unit[3]; - double N_boxes = N / 3.0; - fl_units[1] = round( - pow( - (N_boxes * box[1]*box[1]) - / (this->box[0] * this->box[2]), 1.0/3.0 - ) - ); - if(fl_units[1] == 0) fl_units[1] = 1; - double bxbz_id = N_boxes / fl_units[1]; - fl_units[0] = round(sqrt(this->box[0] * bxbz_id / this->box[2])); - if(fl_units[0] == 0) fl_units[0] = 1; - fl_units[2] = ceil(bxbz_id / fl_units[0]); - for(int d=0; d < 3; d++) fl_unit[d] = box[d] / (double)fl_units[d]; - std::cout << "Unit cell dimensions: " << fl_unit[0] << " x " << fl_unit[1] << " x " << fl_unit[2] << ".\n"; - bool fill[fl_units[0]][fl_units[1]][fl_units[2]][3]; - unsigned slots = 3 * fl_units[0] * fl_units[1] * fl_units[2]; - /* - double pfill = (double)N / slots; - unsigned N1 = 0; - for(unsigned i=0; i < fl_units[0]; i++) - for(unsigned j=0; j < fl_units[1]; j++) - for(unsigned k=0; k < fl_units[2]; k++) - for(int d=0; d < 3; d++) - { - bool tfill = (pfill >= r->rnd()); - fill[i][j][k][d] = tfill; - if(tfill) N1++; - } - */ - unsigned N1 = slots; - for(unsigned i=0; i < fl_units[0]; i++) - for(unsigned j=0; j < fl_units[1]; j++) - for(unsigned k=0; k < fl_units[2]; k++) - for(int d=0; d < 3; d++) fill[i][j][k][d] = true; - - bool tswap; - double pswap; - for(int m=0; m < PRECISION; m++) - { - tswap = (N1 < N); - pswap = ((double)N - (double)N1) / ((tswap? slots: 0) - (double)N1); - // std::cout << "(N = " << N << ", N1 = " << N1 << ", tswap = " << tswap << ", pswap = " << pswap << ")\n"; - for(unsigned i=0; i < fl_units[0]; i++) - for(unsigned j=0; j < fl_units[1]; j++) - for(unsigned k=0; k < fl_units[2]; k++) - for(int d=0; d < 3; d++) - if(pswap >= r->rnd()) - { - if(fill[i][j][k][d]) N1--; - fill[i][j][k][d] = tswap; - if(tswap) N1++; - } - } - std::cout << "Filling " << N1 << " of 3*" - << fl_units[0] << "*" << fl_units[1] << "*" << fl_units[2] - << " = " << 3*fl_units[0]*fl_units[1]*fl_units[2] - << " slots (ideally " << N << ").\n"; - - double LJ_CUTOFF; - double EL_CUTOFF; - double FLUIDMASS, EPS_FLUID, SIG_FLUID, FLUIDLONG, QDR_FLUID; - if(fluid == FLUID_AR) - { - FLUIDMASS = ARMASS; - EPS_FLUID = EPS_AR; - SIG_FLUID = SIG_AR; - LJ_CUTOFF = CUT_AR; - } - else if(fluid == FLUID_CH4) - { - FLUIDMASS = CH4MASS; - EPS_FLUID = EPS_CH4; - SIG_FLUID = SIG_CH4; - LJ_CUTOFF = CUT_CH4; - } - else if(fluid == FLUID_C2H6) - { - FLUIDMASS = C2H6MASS; - EPS_FLUID = EPS_C2H6; - SIG_FLUID = SIG_C2H6; - FLUIDLONG = C2H6LONG; - QDR_FLUID = QDR_C2H6; - LJ_CUTOFF = CUT_C2H6; - } - else if(fluid == FLUID_N2) - { - FLUIDMASS = N2MASS; - EPS_FLUID = EPS_N2; - SIG_FLUID = SIG_N2; - FLUIDLONG = N2LONG; - QDR_FLUID = QDR_N2; - LJ_CUTOFF = CUT_N2; - } - else if(fluid == FLUID_CO2) - { - FLUIDMASS = CO2MASS; - EPS_FLUID = EPS_CO2; - SIG_FLUID = SIG_CO2; - FLUIDLONG = CO2LONG; - QDR_FLUID = QDR_CO2; - LJ_CUTOFF = CUT_CO2; - } - else if(fluid == FLUID_EOX) - { - FLUIDMASS = 2*CEOXMASS + OEOXMASS; - LJ_CUTOFF = CUTLJEOX; - } - else if(fluid == FLUID_JES) - { - FLUIDMASS = OJESMASS + 2.0*HJESMASS; - EPS_FLUID = EPS_OJES; - SIG_FLUID = SIG_OJES; - FLUIDLONG = JES_LONG; - LJ_CUTOFF = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid == FLUID_MER) - { - FLUIDMASS = CMERMASS + 2.0*OMERMASS; - FLUIDLONG = MER_LONG; - QDR_FLUID = QDR_CMER; - LJ_CUTOFF = 4.0*SIG_OMER + 0.5*FLUIDLONG; - } - else if(fluid == FLUID_TOL) - { - FLUIDMASS = CH3TOLMASS + CTRTOLMASS + 5.0*CH_TOLMASS; - FLUIDLONG = TOL_LONG; - LJ_CUTOFF = 4.0*SIG_CH_TOL + 0.5*FLUIDLONG; - } - else if(fluid == FLUID_VEG) - { - FLUIDMASS = OVEGMASS + 2.0*HVEGMASS; - EPS_FLUID = EPS_OVEG; - SIG_FLUID = SIG_OVEG; - FLUIDLONG = VEG_LONG; - LJ_CUTOFF = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else - { - std::cout << "Unavailable fluid ID " << fluid << ".\n"; - exit(20); - } - - if((fluid == FLUID_AR) || (fluid == FLUID_CH4)) - { - EL_CUTOFF = LJ_CUTOFF; - } - else EL_CUTOFF = 1.2 * LJ_CUTOFF; - - double LJ_CUTOFF2 = 0.0; - double EL_CUTOFF2 = 0.0; - double FLUIDMASS2, EPS_FLUID2, SIG_FLUID2, FLUIDLONG2, QDR_FLUID2; - if(fluid2 == FLUID_AR) - { - FLUIDMASS2 = ARMASS; - EPS_FLUID2 = EPS_AR; - SIG_FLUID2 = SIG_AR; - LJ_CUTOFF2 = 2.5*SIG_FLUID2; - } - else if(fluid2 == FLUID_CH4) - { - FLUIDMASS2 = CH4MASS; - EPS_FLUID2 = EPS_CH4; - SIG_FLUID2 = SIG_CH4; - LJ_CUTOFF2 = 2.5*SIG_FLUID2; - } - else if(fluid2 == FLUID_C2H6) - { - FLUIDMASS2 = C2H6MASS; - EPS_FLUID2 = EPS_C2H6; - SIG_FLUID2 = SIG_C2H6; - FLUIDLONG2 = C2H6LONG; - QDR_FLUID2 = QDR_C2H6; - LJ_CUTOFF2 = 4.0*SIG_FLUID2 + 0.5*FLUIDLONG2; - } - else if(fluid2 == FLUID_N2) - { - FLUIDMASS2 = N2MASS; - EPS_FLUID2 = EPS_N2; - SIG_FLUID2 = SIG_N2; - FLUIDLONG2 = N2LONG; - QDR_FLUID2 = QDR_N2; - LJ_CUTOFF2 = 4.0*SIG_FLUID2 + 0.5*FLUIDLONG2; - } - else if(fluid2 == FLUID_CO2) - { - FLUIDMASS2 = CO2MASS; - EPS_FLUID2 = EPS_CO2; - SIG_FLUID2 = SIG_CO2; - FLUIDLONG2 = CO2LONG; - QDR_FLUID2 = QDR_CO2; - LJ_CUTOFF2 = 4.0*SIG_FLUID2 + 0.5*FLUIDLONG2; - } - else if(fluid2 == FLUID_EOX) - { - FLUIDMASS2 = 2*CEOXMASS + OEOXMASS; - LJ_CUTOFF2 = CUTLJEOX; - } - else if(fluid2 == FLUID_JES) - { - FLUIDMASS2 = OJESMASS + 2.0*HJESMASS; - EPS_FLUID2 = EPS_OJES; - SIG_FLUID2 = SIG_OJES; - FLUIDLONG2 = JES_LONG; - LJ_CUTOFF2 = 4.0*SIG_FLUID2 + 0.5*FLUIDLONG2; - } - else if(fluid2 == FLUID_MER) - { - FLUIDMASS2 = CMERMASS + 2.0*OMERMASS; - FLUIDLONG2 = MER_LONG; - QDR_FLUID2 = QDR_CMER; - LJ_CUTOFF2 = 4.0*SIG_OMER + 0.5*FLUIDLONG2; - } - else if(fluid2 == FLUID_TOL) - { - FLUIDMASS2 = CH3TOLMASS + CTRTOLMASS + 5.0*CH_TOLMASS; - FLUIDLONG2 = TOL_LONG; - LJ_CUTOFF2 = 4.0*SIG_CH_TOL + 0.5*FLUIDLONG2; - } - else if(fluid2 == FLUID_VEG) - { - FLUIDMASS2 = OVEGMASS + 2.0*HVEGMASS; - EPS_FLUID2 = EPS_OVEG; - SIG_FLUID2 = SIG_OVEG; - FLUIDLONG2 = VEG_LONG; - LJ_CUTOFF2 = 4.0*SIG_FLUID2 + 0.5*FLUIDLONG2; - } - else fluid2 = FLUID_NIL; - - if((fluid2 == FLUID_AR) || (fluid2 == FLUID_CH4)) - { - EL_CUTOFF2 = LJ_CUTOFF2; - } - else if(fluid2 != FLUID_NIL) EL_CUTOFF2 = 1.2*LJ_CUTOFF2; - - if(LJ_CUTOFF2 > LJ_CUTOFF) LJ_CUTOFF = LJ_CUTOFF2; - if(EL_CUTOFF2 > EL_CUTOFF) EL_CUTOFF = EL_CUTOFF2; - for(int dim=0; dim < 3; dim++) - { - if(0.5*box[dim] < EL_CUTOFF) EL_CUTOFF = 0.5*box[dim]; - if(0.5*box[dim] < LJ_CUTOFF) LJ_CUTOFF = 0.5*box[dim]; - } - - unsigned fluidcomp = (fluid2 != FLUID_NIL)? 2: 1; - - xdr.precision(9); - if(format == FORMAT_BRANCH) - { - xdr << "mardyn " << TIME << " tersoff\n" - << "# mardyn input file, ls1 project\n" - << "# written by animake, the mesh generator\n"; - } - if(format == FORMAT_BUCHHOLZ) - { - xdr << "mardyn trunk " << TIME << "\n"; - } - - if(format == FORMAT_BERNREUTHER) - { - txt << "\n"; - - txt << "" - << DT/REFTIME << "\n" - << "0100000000\n" - << "1" - << LJ_CUTOFF/SIG_REF << "1.0e+10\n" - << "\n" - << "100" - << prefix << "_1R\n" - << "10000" - << prefix << "_1R\n" - << "\n" - << "" - << T/EPS_REF << "\n" - << "" << box[0]/SIG_REF - << "" << box[1]/SIG_REF << "" - << box[2]/SIG_REF << "\n" - << "\n"; - - txt << "\n"; - if((fluid == FLUID_AR) || (fluid == FLUID_CH4)) - { - txt << "000" - << FLUIDMASS/REFMASS << "" - << SIG_FLUID/SIG_REF << "" - << EPS_FLUID/EPS_REF - << "false\n"; - } - else if(fluid == FLUID_EOX) - { - txt << "" - << R0_C1EOX/SIG_REF << "" << R1_C1EOX/SIG_REF - << "" << R2_C1EOX/SIG_REF << "" - << CEOXMASS/REFMASS << "" - << SIG_CEOX/SIG_REF << "" - << EPS_CEOX/EPS_REF - << "false\n" - << "" - << R0_C2EOX/SIG_REF << "" << R1_C2EOX/SIG_REF - << "" << R2_C2EOX/SIG_REF << "" - << CEOXMASS/REFMASS << "" - << SIG_CEOX/SIG_REF << "" - << EPS_CEOX/EPS_REF - << "false\n" - << "" - << R0_O_EOX/SIG_REF << "" << R1_O_EOX/SIG_REF - << "" << R2_O_EOX/SIG_REF << "" - << OEOXMASS/REFMASS << "" - << SIG_OEOX/SIG_REF << "" - << EPS_OEOX/EPS_REF - << "false\n" - << "" - << R0DIPEOX/SIG_REF << "" << R1DIPEOX/SIG_REF - << "" << R2DIPEOX/SIG_REF - << "001" - << DIPOLEOX/DIP_REF << "\n"; - } - txt << "\n"; - - txt << "\n" - << "" << prefix - << ".inp\n" - << "\n"; - } - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << "T\t" << T/EPS_REF << "\n"; - xdr << "t\t0.0\nL\t" << box[0]/SIG_REF << "\t" - << box[1]/SIG_REF << "\t" << box[2]/SIG_REF - << "\nC\t" << fluidcomp << "\n"; - } - if(format == FORMAT_BRANCH) - { - if((fluid == FLUID_AR) || (fluid == FLUID_CH4)) - { - xdr << "1 0 0 0 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid == FLUID_EOX) - { - xdr << "3 0 0 1 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid == FLUID_JES) - { - xdr << "1 3 0 0 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid == FLUID_MER) - { - xdr << "3 0 1 0 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid == FLUID_TOL) - { - xdr << "7 0 5 1 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid == FLUID_VEG) - { - xdr << "1 3 0 0 0\n"; // LJ, C, Q, D, Tersoff - } - else - { - xdr << "2 0 1 0 0\n"; // LJ, C, Q, D, Tersoff - } - } - if(format == FORMAT_BUCHHOLZ) - { - if((fluid == FLUID_AR) || (fluid == FLUID_CH4)) - { - xdr << "1 0 0 0 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid == FLUID_EOX) - { - xdr << "3 0 1 0 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid == FLUID_JES) - { - xdr << "1 3 0 0 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid == FLUID_MER) - { - xdr << "3 0 0 1 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid == FLUID_TOL) - { - xdr << "7 0 1 5 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid == FLUID_VEG) - { - xdr << "1 3 0 0 0\n"; // LJ, C, D, Q, Tersoff - } - else - { - xdr << "2 0 0 1 0\n"; // LJ, C, D, Q, Tersoff - } - } - - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - if((fluid == FLUID_AR) || (fluid == FLUID_CH4)) - { - xdr << "0.0 0.0 0.0\t" - << FLUIDMASS/REFMASS << " " << EPS_FLUID/EPS_REF << " " - << SIG_FLUID/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\t0.0 0.0 0.0\n"; - } - else if(fluid == FLUID_EOX) - { - xdr << R0_C1EOX/SIG_REF << " " << R1_C1EOX/SIG_REF << " " - << R2_C1EOX/SIG_REF << "\t" << CEOXMASS/REFMASS << " " - << EPS_CEOX/EPS_REF << " " << SIG_CEOX/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << R0_C2EOX/SIG_REF << " " << R1_C2EOX/SIG_REF << " " - << R2_C2EOX/SIG_REF << "\t" << CEOXMASS/REFMASS << " " - << EPS_CEOX/EPS_REF << " " << SIG_CEOX/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << R0_O_EOX/SIG_REF << " " << R1_O_EOX/SIG_REF << " " - << R2_O_EOX/SIG_REF << "\t" << OEOXMASS/REFMASS << " " - << EPS_OEOX/EPS_REF << " " << SIG_OEOX/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n"; - - xdr << R0DIPEOX/SIG_REF << " " << R1DIPEOX/SIG_REF << " " - << R2DIPEOX/SIG_REF << "\t0.0 0.0 1.0 " - << DIPOLEOX/DIP_REF; - - xdr << "\n0.0 0.0 0.0\n"; - } - else if(fluid == FLUID_JES) - { - xdr << R0_O_JES/SIG_REF << " " << R1_O_JES/SIG_REF << " " << R2_O_JES/SIG_REF << "\t" - << OJESMASS/REFMASS << " " << EPS_OJES/EPS_REF << " " << SIG_OJES/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << R0_H1JES/SIG_REF << " " << R1_H1JES/SIG_REF << " " << R2_H1JES/SIG_REF << "\t" - << HJESMASS/REFMASS << " " << CHG_HJES/REFCARG << "\n"; - xdr << R0_H2JES/SIG_REF << " " << R1_H2JES/SIG_REF << " " << R2_H2JES/SIG_REF << "\t" - << HJESMASS/REFMASS << " " << CHG_HJES/REFCARG << "\n"; - xdr << R0_E_JES/SIG_REF << " " << R1_E_JES/SIG_REF << " " << R2_E_JES/SIG_REF << "\t" - << "0.0 " << CHG_EJES/REFCARG << "\n"; - - xdr << "0.0 0.0 0.0\n"; - } - else if(fluid == FLUID_MER) - { - xdr << "0.0 0.0 " << -0.5*MER_LONG/SIG_REF << "\t" - << OMERMASS/REFMASS << " " << EPS_OMER/EPS_REF - << " " << SIG_OMER/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << "0.0 0.0 " << +0.5*MER_LONG/SIG_REF << "\t" - << OMERMASS/REFMASS << " " << EPS_OMER/EPS_REF - << " " << SIG_OMER/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << "0.0 0.0 0.0\t" << CMERMASS/REFMASS << " " - << EPS_CMER/EPS_REF << " " << SIG_CMER/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n"; - - xdr << "0.0 0.0 0.0\t0.0 0.0 1.0\t" << QDR_CMER/QDR_REF; - - xdr << "\n0.0 0.0 0.0\n"; - } - else if(fluid == FLUID_TOL) - { - xdr << "0.0 " << R1_CH3_TOL/SIG_REF << " " - << R2_CH3_TOL/SIG_REF << "\t" << CH3TOLMASS/REFMASS - << " " << EPS_CH3TOL/EPS_REF << " " - << SIG_CH3TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << "0.0 " << R1_CTR_TOL/SIG_REF << " " - << R2_CTR_TOL/SIG_REF << "\t" << CTRTOLMASS/REFMASS - << " " << EPS_CTRTOL/EPS_REF << " " - << SIG_CTRTOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << "0.0 " << R1_CHA_TOL/SIG_REF << " " - << R2_CHA_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << "0.0 " << R1_CHB_TOL/SIG_REF << " " - << R2_CHB_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << "0.0 " << R1_CHC_TOL/SIG_REF << " " - << R2_CHC_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << "0.0 " << R1_CHD_TOL/SIG_REF << " " - << R2_CHD_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << "0.0 " << R1_CHE_TOL/SIG_REF << " " - << R2_CHE_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - if(format == FORMAT_BUCHHOLZ) - { - xdr << "0.0 " << R1_CTR_TOL/SIG_REF << " " - << R2_CTR_TOL/SIG_REF << "\t0.0 0.0 -1.0 " - << DIP_CTRTOL/DIP_REF << "\n"; - } - - xdr << "0.0 " << R1_CHA_TOL/SIG_REF << " " - << R2_CHA_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - xdr << "0.0 " << R1_CHB_TOL/SIG_REF << " " - << R2_CHB_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - xdr << "0.0 " << R1_CHC_TOL/SIG_REF << " " - << R2_CHC_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - xdr << "0.0 " << R1_CHD_TOL/SIG_REF << " " - << R2_CHD_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - xdr << "0.0 " << R1_CHE_TOL/SIG_REF << " " - << R2_CHE_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - - if(format == FORMAT_BRANCH) - { - xdr << "0.0 " << R1_CTR_TOL/SIG_REF << " " - << R2_CTR_TOL/SIG_REF << "\t0.0 0.0 -1.0 " - << DIP_CTRTOL/DIP_REF << "\n"; - } - - xdr << "0.0 0.0 0.0\n"; - } - else if(fluid == FLUID_VEG) - { - xdr << R0_O_VEG/SIG_REF << " " << R1_O_VEG/SIG_REF << " " << R2_O_VEG/SIG_REF << "\t" - << OVEGMASS/REFMASS << " " << EPS_OVEG/EPS_REF << " " << SIG_OVEG/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << R0_H1VEG/SIG_REF << " " << R1_H1VEG/SIG_REF << " " << R2_H1VEG/SIG_REF << "\t" - << HVEGMASS/REFMASS << " " << CHG_HVEG/REFCARG << "\n"; - xdr << R0_H2VEG/SIG_REF << " " << R1_H2VEG/SIG_REF << " " << R2_H2VEG/SIG_REF << "\t" - << HVEGMASS/REFMASS << " " << CHG_HVEG/REFCARG << "\n"; - xdr << R0_E_VEG/SIG_REF << " " << R1_E_VEG/SIG_REF << " " << R2_E_VEG/SIG_REF << "\t" - << "0.0 " << CHG_EVEG/REFCARG << "\n"; - - xdr << "0.0 0.0 0.0\n"; - } - else - { - xdr << "0.0 0.0 " << -0.5*FLUIDLONG/SIG_REF << "\t" - << 0.5*FLUIDMASS/REFMASS << " " << EPS_FLUID/EPS_REF - << " " << SIG_FLUID/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << "0.0 0.0 " << +0.5*FLUIDLONG/SIG_REF << "\t" - << 0.5*FLUIDMASS/REFMASS << " " << EPS_FLUID/EPS_REF - << " " << SIG_FLUID/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n"; - - xdr << "0.0 0.0 0.0\t0.0 0.0 1.0\t" << QDR_FLUID/QDR_REF; - - xdr << "\n0.0 0.0 0.0\n"; - } - } - - if(format == FORMAT_BRANCH) - { - if((fluid2 == FLUID_AR) || (fluid2 == FLUID_CH4)) - { - xdr << "1 0 0 0 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid2 == FLUID_EOX) - { - xdr << "3 0 0 1 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid2 == FLUID_JES) - { - xdr << "1 3 0 0 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid2 == FLUID_MER) - { - xdr << "3 0 1 0 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid2 == FLUID_TOL) - { - xdr << "7 0 5 1 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid2 == FLUID_VEG) - { - xdr << "1 3 0 0 0\n"; // LJ, C, Q, D, Tersoff - } - else if(fluid2 != FLUID_NIL) - { - xdr << "2 0 1 0 0\n"; // LJ, C, Q, D, Tersoff - } - } - if(format == FORMAT_BUCHHOLZ) - { - if((fluid2 == FLUID_AR) || (fluid2 == FLUID_CH4)) - { - xdr << "1 0 0 0 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid2 == FLUID_EOX) - { - xdr << "3 0 1 0 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid2 == FLUID_JES) - { - xdr << "1 3 0 0 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid2 == FLUID_MER) - { - xdr << "3 0 0 1 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid2 == FLUID_TOL) - { - xdr << "7 0 1 5 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid2 == FLUID_VEG) - { - xdr << "1 3 0 0 0\n"; // LJ, C, D, Q, Tersoff - } - else if(fluid2 != FLUID_NIL) - { - xdr << "2 0 0 1 0\n"; // LJ, C, D, Q, Tersoff - } - } - - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - if((fluid2 == FLUID_AR) || (fluid2 == FLUID_CH4)) - { - xdr << "0.0 0.0 0.0\t" - << FLUIDMASS2/REFMASS << " " << EPS_FLUID2/EPS_REF << " " - << SIG_FLUID2/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\t0.0 0.0 0.0\n"; - } - else if(fluid2 == FLUID_EOX) - { - xdr << R0_C1EOX/SIG_REF << " " << R1_C1EOX/SIG_REF << " " - << R2_C1EOX/SIG_REF << "\t" << CEOXMASS/REFMASS << " " - << EPS_CEOX/EPS_REF << " " << SIG_CEOX/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << R0_C2EOX/SIG_REF << " " << R1_C2EOX/SIG_REF << " " - << R2_C2EOX/SIG_REF << "\t" << CEOXMASS/REFMASS << " " - << EPS_CEOX/EPS_REF << " " << SIG_CEOX/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << R0_O_EOX/SIG_REF << " " << R1_O_EOX/SIG_REF << " " - << R2_O_EOX/SIG_REF << "\t" << OEOXMASS/REFMASS << " " - << EPS_OEOX/EPS_REF << " " << SIG_OEOX/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n"; - - xdr << R0DIPEOX/SIG_REF << " " << R1DIPEOX/SIG_REF << " " - << R2DIPEOX/SIG_REF << "\t0.0 0.0 1.0 " - << DIPOLEOX/DIP_REF; - - xdr << "\n0.0 0.0 0.0\n"; - } - else if(fluid2 == FLUID_JES) - { - xdr << R0_O_JES/SIG_REF << " " << R1_O_JES/SIG_REF << " " << R2_O_JES/SIG_REF << "\t" - << OJESMASS/REFMASS << " " << EPS_OJES/EPS_REF << " " << SIG_OJES/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << R0_H1JES/SIG_REF << " " << R1_H1JES/SIG_REF << " " << R2_H1JES/SIG_REF << "\t" - << HJESMASS/REFMASS << " " << CHG_HJES/REFCARG << "\n"; - xdr << R0_H2JES/SIG_REF << " " << R1_H2JES/SIG_REF << " " << R2_H2JES/SIG_REF << "\t" - << HJESMASS/REFMASS << " " << CHG_HJES/REFCARG << "\n"; - xdr << R0_E_JES/SIG_REF << " " << R1_E_JES/SIG_REF << " " << R2_E_JES/SIG_REF << "\t" - << "0.0 " << CHG_EJES/REFCARG << "\n"; - - xdr << "0.0 0.0 0.0\n"; - } - else if(fluid2 == FLUID_MER) - { - xdr << "0.0 0.0 " << -0.5*MER_LONG/SIG_REF << "\t" - << OMERMASS/REFMASS << " " << EPS_OMER/EPS_REF - << " " << SIG_OMER/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << "0.0 0.0 " << +0.5*MER_LONG/SIG_REF << "\t" - << OMERMASS/REFMASS << " " << EPS_OMER/EPS_REF - << " " << SIG_OMER/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << "0.0 0.0 0.0\t" << CMERMASS/REFMASS << " " - << EPS_CMER/EPS_REF << " " << SIG_CMER/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n"; - - xdr << "0.0 0.0 0.0\t0.0 0.0 1.0\t" << QDR_CMER/QDR_REF; - - xdr << "\n0.0 0.0 0.0\n"; - } - else if(fluid2 == FLUID_TOL) - { - xdr << "0.0 " << R1_CH3_TOL/SIG_REF << " " - << R2_CH3_TOL/SIG_REF << "\t" << CH3TOLMASS/REFMASS - << " " << EPS_CH3TOL/EPS_REF << " " - << SIG_CH3TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << "0.0 " << R1_CTR_TOL/SIG_REF << " " - << R2_CTR_TOL/SIG_REF << "\t" << CTRTOLMASS/REFMASS - << " " << EPS_CTRTOL/EPS_REF << " " - << SIG_CTRTOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << "0.0 " << R1_CHA_TOL/SIG_REF << " " - << R2_CHA_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << "0.0 " << R1_CHB_TOL/SIG_REF << " " - << R2_CHB_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << "0.0 " << R1_CHC_TOL/SIG_REF << " " - << R2_CHC_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << "0.0 " << R1_CHD_TOL/SIG_REF << " " - << R2_CHD_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << "0.0 " << R1_CHE_TOL/SIG_REF << " " - << R2_CHE_TOL/SIG_REF << "\t" << CH_TOLMASS/REFMASS - << " " << EPS_CH_TOL/EPS_REF << " " - << SIG_CH_TOL/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - if(format == FORMAT_BUCHHOLZ) - { - xdr << "0.0 " << R1_CTR_TOL/SIG_REF << " " - << R2_CTR_TOL/SIG_REF << "\t0.0 0.0 -1.0 " - << DIP_CTRTOL/DIP_REF << "\n"; - } - - xdr << "0.0 " << R1_CHA_TOL/SIG_REF << " " - << R2_CHA_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - xdr << "0.0 " << R1_CHB_TOL/SIG_REF << " " - << R2_CHB_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - xdr << "0.0 " << R1_CHC_TOL/SIG_REF << " " - << R2_CHC_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - xdr << "0.0 " << R1_CHD_TOL/SIG_REF << " " - << R2_CHD_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - xdr << "0.0 " << R1_CHE_TOL/SIG_REF << " " - << R2_CHE_TOL/SIG_REF << "\t" << "1.0 0.0 0.0\t" - << QDR_CH_TOL/QDR_REF << "\n"; - - if(format == FORMAT_BRANCH) - { - xdr << "0.0 " << R1_CTR_TOL/SIG_REF << " " - << R2_CTR_TOL/SIG_REF << "\t0.0 0.0 -1.0 " - << DIP_CTRTOL/DIP_REF << "\n"; - } - - xdr << "0.0 0.0 0.0\n"; - } - else if(fluid2 == FLUID_VEG) - { - xdr << R0_O_VEG/SIG_REF << " " << R1_O_VEG/SIG_REF << " " << R2_O_VEG/SIG_REF << "\t" - << OVEGMASS/REFMASS << " " << EPS_OVEG/EPS_REF << " " << SIG_OVEG/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << R0_H1VEG/SIG_REF << " " << R1_H1VEG/SIG_REF << " " << R2_H1VEG/SIG_REF << "\t" - << HVEGMASS/REFMASS << " " << CHG_HVEG/REFCARG << "\n"; - xdr << R0_H2VEG/SIG_REF << " " << R1_H2VEG/SIG_REF << " " << R2_H2VEG/SIG_REF << "\t" - << HVEGMASS/REFMASS << " " << CHG_HVEG/REFCARG << "\n"; - xdr << R0_E_VEG/SIG_REF << " " << R1_E_VEG/SIG_REF << " " << R2_E_VEG/SIG_REF << "\t" - << "0.0 " << CHG_EVEG/REFCARG << "\n"; - - xdr << "0.0 0.0 0.0\n"; - } - else if(fluid2 != FLUID_NIL) - { - xdr << "0.0 0.0 " << -0.5*FLUIDLONG2/SIG_REF << "\t" - << 0.5*FLUIDMASS2/REFMASS << " " << EPS_FLUID2/EPS_REF - << " " << SIG_FLUID2/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n" << "0.0 0.0 " << +0.5*FLUIDLONG2/SIG_REF << "\t" - << 0.5*FLUIDMASS2/REFMASS << " " << EPS_FLUID2/EPS_REF - << " " << SIG_FLUID2/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF/SIG_REF << " 0"; - xdr << "\n"; - - xdr << "0.0 0.0 0.0\t0.0 0.0 1.0\t" << QDR_FLUID2/QDR_REF; - - xdr << "\n0.0 0.0 0.0\n"; - } - - if(fluid2 != FLUID_NIL) - { - xdr << XIF << "\t" << ETAF << "\n"; - } - xdr << "1.0e+10\n"; - - txt.precision(6); - txt << "mardynconfig\n# \ntimestepLength\t" << DT/REFTIME - << "\ncutoffRadius\t" << EL_CUTOFF/SIG_REF - << "\nLJCutoffRadius\t" << LJ_CUTOFF/SIG_REF; - } - if((format == FORMAT_BERNREUTHER) || (format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << "N" << "\t" << N1 << "\nM" << "\t" << "ICRVQD\n"; - } - if(format == FORMAT_BRANCH) - { - txt << "\ntersoffCutoffRadius\t" - << LJ_CUTOFF/(4.0*SIG_REF); - } - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt << "\ninitCanonical\t50000\n"; - if(muVT) - { - txt.precision(9); - txt << "chemicalPotential " << mu/EPS_REF - << " component 1 control 0.0 0.0 0.0 to " - << this->box[0]/SIG_REF << " " << this->box[1]/SIG_REF - << " " << this->box[2]/SIG_REF << " conduct " - << 1 + (int)round(0.003 * N) - << " tests every 3 steps\n"; - txt << "planckConstant\t" - << sqrt(6.28319 * T/EPS_REF) << "\n"; // sqrt(2 pi kT) - txt << "initGrandCanonical\t100000\n"; - } - txt.precision(5); - txt << "initStatistics\t150000\n"; - } - if(format == FORMAT_BRANCH) - { - txt << "phaseSpaceFile\t" << prefix << ".xdr\n"; - } - if(format == FORMAT_BUCHHOLZ) - { - txt << "phaseSpaceFile\tOldStyle\t" << prefix << ".inp\nparallelization\tDomainDecomposition\n"; - } - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt << "datastructure\tLinkedCells\t1\noutput\t" - << "ResultWriter\t100\t" << prefix - << "_1R\noutput\tXyzWriter\t10000\t" << prefix - << "_1R.buxyz\n" - << "RDF\t" << (EL_CUTOFF - 4.0*FLUIDLONG)/(SIG_REF * (double)BINS) << " " - << BINS << "\nRDFOutputTimesteps\t150000\nRDFOutputPrefix\t" << prefix << "_1R\n"; - } - if(format == FORMAT_BERNREUTHER) - { - txt << "\n"; - } - txt.close(); - - double I[3]; - for(int d=0; d < 3; d++) I[d] = 0.0; - if(fluid == FLUID_EOX) - { - I[0] = I_XX_EOX; - I[1] = I_YY_EOX; - I[2] = I_ZZ_EOX; - } - else if(fluid == FLUID_JES) - { - I[0] = I_XX_JES; - I[1] = I_YY_JES; - I[2] = I_ZZ_JES; - } - else if(fluid == FLUID_TOL) - { - I[0] = I_XX_TOL; - I[1] = I_YY_TOL; - I[2] = I_ZZ_TOL; - } - else if(fluid == FLUID_VEG) - { - I[0] = I_XX_VEG; - I[1] = I_YY_VEG; - I[2] = I_ZZ_VEG; - } - else if(!((fluid == FLUID_AR) || (fluid == FLUID_CH4))) - { - I[0] = 0.25 * FLUIDMASS * FLUIDLONG * FLUIDLONG; - I[1] = I[0]; - } - - double I2[3]; - for(int k=0; k < 3; k++) I2[k] = 0.0; - if(fluid2 == FLUID_EOX) - { - I2[0] = I_XX_EOX; - I2[1] = I_YY_EOX; - I2[2] = I_ZZ_EOX; - } - else if(fluid2 == FLUID_JES) - { - I2[0] = I_XX_JES; - I2[1] = I_YY_JES; - I2[2] = I_ZZ_JES; - } - else if(fluid2 == FLUID_TOL) - { - I2[0] = I_XX_TOL; - I2[1] = I_YY_TOL; - I2[2] = I_ZZ_TOL; - } - else if(fluid2 == FLUID_VEG) - { - I2[0] = I_XX_VEG; - I2[1] = I_YY_VEG; - I2[2] = I_ZZ_VEG; - } - else if(!((fluid2 == FLUID_AR) || (fluid2 == FLUID_CH4))) - { - I2[0] = 0.25 * FLUIDMASS2 * FLUIDLONG2 * FLUIDLONG2; - I2[1] = I2[0]; - } - - unsigned id = 1; - double tr[3]; - unsigned ii[3]; - double Nf[2]; - Nf[0] = (1.0 - x) * (double)N1 + 0.001; - Nf[1] = x * (double)N1 + 0.001; - for(ii[0]=0; ii[0] < fl_units[0]; (ii[0]) ++) - for(ii[1]=0; ii[1] < fl_units[1]; (ii[1]) ++) - for(ii[2]=0; ii[2] < fl_units[2]; (ii[2]) ++) - for(int d=0; d < 3; d++) - { - if(fill[ ii[0] ][ ii[1] ][ ii[2] ][ d ]) - { - for(int k=0; k < 3; k++) - { - tr[k] = fl_unit[k] * ( - ii[k] + 0.02*r->rnd() + ((k == d)? 0.24: 0.74) - ); - } - for(int k=0; k < 3; k++) - { - if(tr[k] > box[k]) tr[k] -= box[k]; - else if(tr[k] < 0.0) tr[k] += box[k]; - } - double tv = sqrt(3.0*T / FLUIDMASS); - double phi = 6.283185 * r->rnd(); - double omega = 6.283185 * r->rnd(); - unsigned cid; - if(fluidcomp == 1) cid = 1; - else - { - cid = (r->rnd() > (Nf[0] / (Nf[0] + Nf[1])))? 2: 1; - --Nf[cid - 1]; - } - double w[3]; - for(int d=0; d < 3; d++) - if(cid == 1) w[d] = (I[d] == 0)? 0.0: ((r->rnd() > 0.5)? 1: -1) * sqrt(2.0*r->rnd()*T / I[d]); - else w[d] = (I2[d] == 0)? 0.0: ((r->rnd() > 0.5)? 1: -1) * sqrt(2.0*r->rnd()*T / I2[d]); - // xdr << "(" << ii[0] << "/" << ii[1] << "/" << ii[2] << "), j = " << j << ", d = " << d << ":\t"; - if((format == FORMAT_BERNREUTHER) || (format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << id << " " << cid << "\t" << tr[0]/SIG_REF - << " " << tr[1]/SIG_REF << " " << tr[2]/SIG_REF - << "\t" << tv*cos(phi)*cos(omega)/VEL_REF << " " - << tv*cos(phi)*sin(omega)/VEL_REF << " " - << tv*sin(phi)/VEL_REF << "\t1.0 0.0 0.0 0.0\t" - << w[0]/REFOMGA << " " << w[1]/REFOMGA << " " - << w[2]/REFOMGA << "\n"; - } - id++; - } - else xdr << "\n"; - } - xdr << "\n"; - - xdr.close(); -} - diff --git a/tools/standalone-generators/animake/Domain.h b/tools/standalone-generators/animake/Domain.h deleted file mode 100644 index c4852f8832..0000000000 --- a/tools/standalone-generators/animake/Domain.h +++ /dev/null @@ -1,184 +0,0 @@ -#include -#include -#include -#include - - -#define TIME 20130802 - -#define FLUID_NIL -1 -#define FLUID_CH4 0 -#define FLUID_AR 1 -#define FLUID_C2H6 2 -#define FLUID_N2 3 -#define FLUID_CO2 4 -#define FLUID_EOX 5 -#define FLUID_JES 6 -#define FLUID_MER 7 -#define FLUID_TOL 8 -#define FLUID_VEG 9 - -#define EPS_AR 3.69853e-04 -#define SIG_AR 6.41600 -#define ARMASS 0.039948 -#define CUT_AR 38.496 -#define EPS_CH4 4.70431e-04 -#define SIG_CH4 7.04599 -#define CH4MASS 0.016042 -#define CUT_CH4 42.276 - -#define EPS_C2H6 4.33822e-04 -#define SIG_C2H6 6.59439 -#define C2H6MASS 0.03007 -#define C2H6LONG 4.49037 -#define QDR_C2H6 -0.61537 -#define CUT_C2H6 40.689 -#define EPS_N2 1.10512e-04 -#define SIG_N2 6.27597 -#define N2MASS 0.0280134 -#define N2LONG 1.97741 -#define QDR_N2 -1.07038 -#define CUT_N2 38.150 -#define EPS_CO2 4.21883e-04 -#define SIG_CO2 5.64027 -#define CO2MASS 0.0440095 -#define CO2LONG 4.56860 -#define QDR_CO2 -2.82059 -#define CUT_CO2 34.984 - -#define EPS_CEOX 2.68353e-04 -#define SIG_CEOX 6.66431 -#define CEOXMASS 0.014027 -#define EPS_OEOX 1.96742e-04 -#define SIG_OEOX 5.84473 -#define OEOXMASS 0.015999 -#define DIPOLEOX -0.96744 -#define R0_C1EOX +1.47399 -#define R1_C1EOX +0.0 -#define R2_C1EOX -0.83729 -#define R0_C2EOX -1.47399 -#define R1_C2EOX +0.0 -#define R2_C2EOX -0.83729 -#define R0_O_EOX +0.0 -#define R1_O_EOX +0.0 -#define R2_O_EOX +1.46818 -#define R0DIPEOX +0.0 -#define R1DIPEOX +0.0 -#define R2DIPEOX +0.07792 -#define CUTLJEOX 40.833 -#define I_XX_EOX 0.054154 -#define I_YY_EOX 0.11511 -#define I_ZZ_EOX 0.60951 - -#define EPS_OJES 6.58951e-04 -#define SIG_OJES 5.89277 -#define OJESMASS 0.0159994 -#define HJESMASS 0.0010079 -#define CHG_HJES +0.419548 -#define CHG_EJES -0.839096 -#define R0_O_JES +0.0 -#define R1_O_JES -0.14948 -#define R2_O_JES +0.0 -#define R0_H1JES -1.72579 -#define R1_H1JES +1.18644 -#define R2_H1JES +0.0 -#define R0_H2JES +1.72579 -#define R1_H2JES +1.18644 -#define R2_H2JES +0.0 -#define R0_E_JES +0.0 -#define R1_E_JES +0.23765 -#define R2_E_JES +0.0 -#define I_XX_JES 0.0031950 -#define I_YY_JES 0.0060038 -#define I_ZZ_JES 0.0091988 -#define JES_LONG 1.0 - -#define EPS_OMER 3.18243e-04 -#define SIG_OMER 5.62288 -#define OMERMASS 0.0159994 -#define EPS_CMER 3.9181e-05 -#define SIG_CMER 5.31712 -#define QDR_CMER -3.02883 -#define CMERMASS 0.012011 -#define MER_LONG 4.8638 - -#define R1_CH3_TOL +0.0 -#define R2_CH3_TOL -5.1999 -#define R1_CTR_TOL +0.0 -#define R2_CTR_TOL -1.8130 -#define R1_CHA_TOL +2.9707 -#define R2_CHA_TOL -0.8715 -#define R1_CHB_TOL -2.9707 -#define R2_CHB_TOL -0.8715 -#define R1_CHC_TOL +2.9767 -#define R2_CHC_TOL +2.5625 -#define R1_CHD_TOL -2.9767 -#define R2_CHD_TOL +2.5625 -#define R1_CHE_TOL +0.0 -#define R2_CHE_TOL +4.2958 -#define CH3TOLMASS 0.015035 -#define CTRTOLMASS 0.012011 -#define CH_TOLMASS 0.013019 -#define EPS_CH3TOL 3.9107e-04 -#define EPS_CTRTOL 3.4645e-05 -#define EPS_CH_TOL 3.18328e-04 -#define SIG_CH3TOL 6.77656 -#define SIG_CTRTOL 5.2799 -#define SIG_CH_TOL 6.1907 -#define DIP_CTRTOL 0.173144 -#define QDR_CH_TOL -1.2548 -#define I_XX_TOL 1.33752 -#define I_YY_TOL 0.87701 -#define I_ZZ_TOL 0.46050 -#define TOL_LONG 9.4957 - -#define EPS_OVEG 2.9515e-04 -#define SIG_OVEG 5.9695 -#define OVEGMASS 0.0159994 -#define HVEGMASS 0.0010079 -#define CHG_HVEG +0.5564 -#define CHG_EVEG -1.1128 -#define R0_O_VEG +0.0 -#define R1_O_VEG -0.12389 -#define R2_O_VEG +0.0 -#define R0_H1VEG -1.43052 -#define R1_H1VEG +0.98330 -#define R2_H1VEG +0.0 -#define R0_H2VEG +1.43052 -#define R1_H2VEG +0.98330 -#define R2_H2VEG +0.0 -#define R0_E_VEG +0.0 -#define R1_E_VEG +0.16826 -#define R2_E_VEG +0.0 -#define I_XX_VEG 0.0021946 -#define I_YY_VEG 0.0041251 -#define I_ZZ_VEG 0.0063197 -#define VEG_LONG 1.0 - -#define FORMAT_BUCHHOLZ 0 -#define FORMAT_BRANCH 1 -#define FORMAT_BERNREUTHER 2 - -#define FLOW_NONE 0 -#define FLOW_COUETTE 1 -#define FLOW_POISEUILLE 2 - -class Domain -{ - public: - Domain( - int sp_fluid, int sp_fluid2, double* sp_box, double sp_SIG_REF, - double sp_EPS_REF, double sp_REFMASS, bool sp_muVT, - unsigned sp_N, double sp_T, double sp_ETAF, double sp_XIF - ); - void write(char* prefix, int format, double mu, double x); - - private: - bool muVT; - int fluid, fluid2; - unsigned N; - double SIG_REF, EPS_REF, REFMASS, T, ETAF, XIF; - - double box[3]; // offset coordinates for the fluid -}; - diff --git a/tools/standalone-generators/animake/Random.cpp b/tools/standalone-generators/animake/Random.cpp deleted file mode 100644 index 2a96b04a08..0000000000 --- a/tools/standalone-generators/animake/Random.cpp +++ /dev/null @@ -1,30 +0,0 @@ -#include "Random.h" - -Random::Random() { this->init(8624); } - -void Random::init(int seed) -{ - this->ix_muVT = (seed ^ (int)888889999) | (int)1; - this->iy_muVT = seed ^ (int)777755555; - // Calculate normalization factor - this->am_muVT = 2.0 / (1.0 + (unsigned)((int)-1)); -} - -float Random::rnd() -{ - float rnd; - const int IA = 16807; - const int IM = 2147483647; - const int IQ = 127773; - const int IR = 2836; - - ix_muVT ^= (ix_muVT >> 13); - ix_muVT ^= (ix_muVT << 17); - ix_muVT ^= (ix_muVT >> 5); - - int k = iy_muVT / IQ; - iy_muVT = IA * (iy_muVT - k*IQ) - IR*k; - if(iy_muVT < 0) iy_muVT += IM; - rnd = am_muVT * ((IM & (ix_muVT ^ iy_muVT)) | (int)1); - return rnd; -} diff --git a/tools/standalone-generators/animake/Random.h b/tools/standalone-generators/animake/Random.h deleted file mode 100644 index 073f88e2e3..0000000000 --- a/tools/standalone-generators/animake/Random.h +++ /dev/null @@ -1,14 +0,0 @@ -class Random -{ - public: - Random(); - void init(int seed); - - float rnd(); - - int getIX() { return this->ix_muVT; } - - private: - int ix_muVT, iy_muVT; - float am_muVT; -}; diff --git a/tools/standalone-generators/animake/main.cpp b/tools/standalone-generators/animake/main.cpp deleted file mode 100644 index ff9bf06e82..0000000000 --- a/tools/standalone-generators/animake/main.cpp +++ /dev/null @@ -1,265 +0,0 @@ -#include "Domain.h" - -#include -#include -#include -#include - - -int main(int argc, char** argv) -{ - const char* usage = "usage: animake -c [-e] [-f ] [-g ] [-J ] [-m ] -N [-r] [-s ] -T [-u] [-W ] [-x <2nd comp. mole fract.>] [-Y ] [-y [-z ]] [-5 ]\n\n-e\tuse B-e-rnreuther format\n-f\tCH4 (default), Ar, C2H6, N2, CO2, EOX, JES, MER, TOL or VEG\n-r\tuse b-r-anch format\n-s\tgiven in units of Angstrom; default: 1 = 0.5291772 A\n-u\tuse B-u-chholz format (active by default)\n-W\tgiven in units of K; default value: 1 = 315774.5 K\n-Y\tgiven in units of g/mol; default value: 1 = 1000 g/mol\n\n"; - if((argc < 8) || (argc > 23)) - { - std::cout << "There are " << argc - << " arguments where 8 to 23 should be given.\n\n"; - std::cout << usage; - return 1; - } - - unsigned N; - double x, dimx, dimy, dimz, T, rho, XIF, ETAF; - double mu = 0.0; - int fluid; - int fluid2 = FLUID_NIL; - int format = FORMAT_BUCHHOLZ; - char* prefix = argv[1]; - bool muVT = false; - - double SIG_REF = 1.0; - double EPS_REF = 1.0; - double REFMASS = 1.0; - - bool in_N = false; - bool in_rho = false; - bool in_fluid = false; - bool in_fluid2 = false; - bool in_x = false; - bool in_dimy = false; - bool in_dimz = false; - bool in_T = false; - bool in_ETAF = false; - bool in_XIF = false; - - for(int i=2; i < argc; i++) - { - if(*argv[i] != '-') - { - std::cout << "Flag expected where '" << argv[i] - << "' was given.\n\n"; - std::cout << usage; - return 2; - } - for(int j=1; argv[i][j]; j++) - { - if(argv[i][j] == 'c') - { - in_rho = true; - i++; - rho = atof(argv[i]); - break; - } - else if(argv[i][j] == 'e') format = FORMAT_BERNREUTHER; - else if(argv[i][j] == 'f') - { - in_fluid = true; - i++; - if(!strcmp(argv[i], "CH4")) fluid = FLUID_CH4; - else if(!strcmp(argv[i], "Ar")) fluid = FLUID_AR; - else if(!strcmp(argv[i], "C2H6")) fluid = FLUID_C2H6; - else if(!strcmp(argv[i], "N2")) fluid = FLUID_N2; - else if(!strcmp(argv[i], "CO2")) fluid = FLUID_CO2; - else if(!strcmp(argv[i], "EOX")) fluid = FLUID_EOX; - else if(!strcmp(argv[i], "JES")) fluid = FLUID_JES; - else if(!strcmp(argv[i], "MER")) fluid = FLUID_MER; - else if(!strcmp(argv[i], "TOL")) fluid = FLUID_TOL; - else if(!strcmp(argv[i], "VEG")) fluid = FLUID_VEG; - else - { - std::cout << "Fluid '" << argv[i] - << "' is not available.\n\n" << usage; - return 9; - } - break; - } - else if(argv[i][j] == 'g') - { - in_fluid2 = true; - i++; - if(!strcmp(argv[i], "Ar")) fluid2 = FLUID_AR; - else if(!strcmp(argv[i], "CH4")) fluid2 = FLUID_CH4; - else if(!strcmp(argv[i], "C2H6")) fluid2 = FLUID_C2H6; - else if(!strcmp(argv[i], "N2")) fluid2 = FLUID_N2; - else if(!strcmp(argv[i], "CO2")) fluid2 = FLUID_CO2; - else if(!strcmp(argv[i], "EOX")) fluid2 = FLUID_EOX; - else if(!strcmp(argv[i], "JES")) fluid2 = FLUID_JES; - else if(!strcmp(argv[i], "MER")) fluid2 = FLUID_MER; - else if(!strcmp(argv[i], "TOL")) fluid2 = FLUID_TOL; - else if(!strcmp(argv[i], "VEG")) fluid2 = FLUID_VEG; - else - { - std::cout << "(Secondary) fluid '" << argv[i] - << "' is not available.\n\n" << usage; - return 99; - } - break; - } - else if(argv[i][j] == 'J') - { - in_ETAF = true; - i++; - ETAF = atof(argv[i]); - break; - } - else if(argv[i][j] == 'm') - { - muVT = true; - i++; - mu = atof(argv[i]); - break; - } - else if(argv[i][j] == 'N') - { - in_N = true; - i++; - N = atoi(argv[i]); - break; - } - else if(argv[i][j] == 'r') format = FORMAT_BRANCH; - else if(argv[i][j] == 's') - { - i++; - SIG_REF = atof(argv[i]) / 0.5291772; - break; - } - else if(argv[i][j] == 'T') - { - in_T = true; - i++; - T = atof(argv[i]); - break; - } - else if(argv[i][j] == 'u') format = FORMAT_BUCHHOLZ; - else if(argv[i][j] == 'W') - { - i++; - EPS_REF = atof(argv[i]) / 315774.5; - break; - } - else if(argv[i][j] == 'x') - { - in_x = true; - i++; - x = atof(argv[i]); - break; - } - else if(argv[i][j] == 'Y') - { - i++; - REFMASS = atof(argv[i]) / 1000.0; - break; - } - else if(argv[i][j] == 'y') - { - in_dimy = true; - i++; - dimy = atof(argv[i]); - break; - } - else if(argv[i][j] == 'z') - { - in_dimz = true; - i++; - dimz = atof(argv[i]); - break; - } - else if(argv[i][j] == '5') - { - in_XIF = true; - i++; - XIF = atof(argv[i]); - break; - } - } - } - - if(!in_ETAF) ETAF = 1.0; - if(!in_XIF) XIF = 1.0; - if(!in_XIF && in_fluid2 && (((fluid == FLUID_MER) && (fluid2 == FLUID_TOL)) || ((fluid == FLUID_TOL) && (fluid2 == FLUID_MER)))) - { - std::cerr << "Default binary parameter for the modified Berthelot mixing rule with CO2 (Merker) and toluene xi = 0.95.\n"; - XIF = 0.95; - } - if(!in_rho) - { - std::cout << "Fatal error: no fluid density was specified.\n\n"; - std::cout << usage; - return 16; - } - if(in_dimz && !in_dimy) - { - std::cout << "Fatal error: z is specified while y is unknown. Please replace the -z option by -y.\n\n"; - std::cout << usage; - return 17; - } - if(!in_fluid) fluid = FLUID_CH4; - if(!in_N) - { - std::cout << "The essential input parameter " - << "N (number of fluid molecules) is missing.\n\n" << usage; - return 20; - } - if(!in_T) - { - std::cout << "Unspecified temperature: aborting.\n\n"; - std::cout << usage; - return 21; - } - if(in_x && ((x < 0.0) || (x > 1.0))) - { - std::cout << "Invalid mole fraction x = " << x << ".\n\n" << usage; - return 15; - } - if((fluid2 != FLUID_NIL) && !in_x) - { - std::cout << "Unspecified composition: aborting.\n\n"; - std::cout << usage; - return 24; - } - if((in_fluid2 == false) || (fluid2 == FLUID_NIL)) - { - x = 0.0; - } - if(x == 1.0) - { - fluid = fluid2; - x = 0.0; - } - if(x == 0.0) - { - in_fluid2 = false; - fluid2 = FLUID_NIL; - } - - double V = (double)N/rho; - if(!in_dimy) dimy = pow(V, 1.0/3.0); - if(!in_dimz) dimz = sqrt(V/dimy); - dimx = V / (dimy*dimz); - - double box[3]; - box[0] = dimx * SIG_REF; - box[1] = dimy * SIG_REF; - box[2] = dimz * SIG_REF; - T *= EPS_REF; - mu *= EPS_REF; - - Domain* delta = new Domain( - fluid, fluid2, box, SIG_REF, EPS_REF, REFMASS, muVT, N, T, ETAF, XIF - ); - delta->write( - prefix, format, mu, x - ); - - return 0; -} - diff --git a/tools/standalone-generators/atomic_units.txt b/tools/standalone-generators/atomic_units.txt deleted file mode 100644 index ae233e8d81..0000000000 --- a/tools/standalone-generators/atomic_units.txt +++ /dev/null @@ -1,21 +0,0 @@ -Atomic system of units for mkcp (without -L flag): - -Unit length 1 = 1 a0 = 0.0529177 nm -Elementary charge 1 = 1 e = 96485.3 C/mol -Coulomb constant 1 = 1/(4 pi eps_0) = 1.43996e-09 eV m/e^2 -Boltzmann constant 1 = 1.38065e-23 J/K -Unit mass 1 = 1 kg/mol = 1000 u - -Derived units: - -Unit temperature 1 = 315774.5 K -Unit density 1 = 11205.9 mol/l -Unit velocity 1 = 1620.34 m/s -Unit time 1 = 32.6584 fs -Unit energy 1 = 27.2126 eV -Unit acceleration 1 = 49614800 nm/ns^2 -Unit pressure 1 = 2.94210e+10 kPa - -Unit dipole moment 1 = 2.54175 Debye -Unit quadrupole moment 1 = 1.34504 Buckingham (i.e. Debye Ang) - diff --git a/tools/standalone-generators/clean b/tools/standalone-generators/clean deleted file mode 100755 index 1bcafe04c7..0000000000 --- a/tools/standalone-generators/clean +++ /dev/null @@ -1,3 +0,0 @@ -make clean; rm -rf CMakeFiles */CMakeFiles _CPack_Packages install_manifest.txt *cmake */*cmake Makefile */Makefile CMakeCache.txt - - diff --git a/tools/standalone-generators/install b/tools/standalone-generators/install deleted file mode 100755 index 671a37062e..0000000000 --- a/tools/standalone-generators/install +++ /dev/null @@ -1,2 +0,0 @@ -./produce; sudo dpkg -i *.deb - diff --git a/tools/standalone-generators/mkcp/CMakeLists.txt b/tools/standalone-generators/mkcp/CMakeLists.txt deleted file mode 100644 index 58c8f15b48..0000000000 --- a/tools/standalone-generators/mkcp/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -cmake_minimum_required(VERSION 3.3) -message(STATUS "cmake building mkcp") - -add_executable(mkcp main.cpp Domain.cpp Graphit.cpp Random.cpp) -target_include_directories(mkcp PRIVATE ${PROJECT_SOURCE_DIR}) - -include(GNUInstallDirs) -install(TARGETS mkcp RUNTIME DESTINATION /usr/bin/) - diff --git a/tools/standalone-generators/mkcp/Domain.cpp b/tools/standalone-generators/mkcp/Domain.cpp deleted file mode 100644 index af516d2633..0000000000 --- a/tools/standalone-generators/mkcp/Domain.cpp +++ /dev/null @@ -1,1356 +0,0 @@ -#include "Domain.h" -#include "Random.h" -#include "Graphit.h" -#include - -#define DT 0.030620 -#define TAU_ZERO 1.0e+10 -#define ZETA 0.0036143 -#define TAUPRIME 61.240 -#define PRECISION 4 - -Domain::Domain( - int sp_flow, double sp_bondlength, double sp_rho, int sp_d, - int sp_fluid, int sp_fluid2, double sp_h, double sp_ETA, double sp_ETA2, double sp_ETAF, double sp_SIG_REF, - double sp_EPS_REF, double sp_REFMASS, bool sp_muVT, - bool sp_nanotube, double sp_m_per_n, unsigned sp_N, double sp_T, - double sp_XI, double sp_XI2, double sp_XIF, double sp_wo_wall -) { - this->flow = sp_flow; - this->bondlength = sp_bondlength; - this->d = sp_d; - this->fluid = sp_fluid; - this->fluid2 = sp_fluid2; - this->h = sp_h; - this->ETA = sp_ETA; - this->ETA2 = sp_ETA2; - this->ETAF = sp_ETAF; - this->SIG_REF = sp_SIG_REF; - this->EPS_REF = sp_EPS_REF; - this->REFMASS = sp_REFMASS; - this->muVT = sp_muVT; - this->nanotube = sp_nanotube; - this->pfill = 0.0; - this->T = sp_T; - this->XI = sp_XI; - this->XI2 = sp_XI2; - this->XIF = sp_XIF; - this->wo_wall = sp_wo_wall; - - if(fluid == FLUID_AR) this->shielding = this->ETA * SIG_AR; - else if(fluid == FLUID_CH4) this->shielding = ETA * SIG_CH4; - /* - * 2CLJQ fluids: shielding = elongation/4 + eta*sigma - */ - else if(fluid == FLUID_C2H6) - { - this->shielding = 0.5*C2H6LONG + this->ETA*SIG_C2H6; - } - else if(fluid == FLUID_N2) - { - this->shielding = 0.5*N2LONG + this->ETA*SIG_N2; - } - else if(fluid == FLUID_CO2) - { - this->shielding = 0.5*CO2LONG + this->ETA*SIG_CO2; - } - else if(fluid == FLUID_AVE) this->shielding = ETA * SIG_AVE; - else if(fluid == FLUID_H2O) - { - this->shielding = 0.5*H2O_LONG + this->ETA*SIG_OH2O; - } - else if(fluid == FLUID_CH3OH) - { - this->shielding = 0.5*CH3OH_LONG + this->ETA*SIG_CCH3OH; - } - else if(fluid == FLUID_C6H14) - { - this->shielding = 0.5*C6H14_LONG + this->ETA*SIG_MC6H14; - } - - double shielding2 = 0.0; - if(fluid2 == FLUID_AR) shielding2 = this->ETA * SIG_AR; - else if(fluid2 == FLUID_CH4) shielding2 = ETA * SIG_CH4; - /* - * 2CLJQ fluids: shielding = elongation/4 + eta*sigma - */ - else if(fluid2 == FLUID_C2H6) - { - shielding2 = 0.5*C2H6LONG + this->ETA*SIG_C2H6; - } - else if(fluid2 == FLUID_N2) - { - shielding2 = 0.5*N2LONG + this->ETA*SIG_N2; - } - else if(fluid2 == FLUID_CO2) - { - shielding2 = 0.5*CO2LONG + this->ETA*SIG_CO2; - } - else if(fluid2 == FLUID_AVE) shielding2 = ETA * SIG_AVE; - else if(fluid2 == FLUID_H2O) - { - shielding2 = 0.5*H2O_LONG + this->ETA*SIG_OH2O; - } - else if(fluid2 == FLUID_CH3OH) - { - shielding2 = 0.5*CH3OH_LONG + this->ETA*SIG_CCH3OH; - } - else if(fluid2 == FLUID_C6H14) - { - shielding2 = 0.5*C6H14_LONG + this->ETA*SIG_MC6H14; - } - if(shielding2 > this->shielding) this->shielding = shielding2; - - if(nanotube) this->specifyNanotube(sp_rho, sp_m_per_n, sp_N); - else this->specifyGraphite(sp_rho, sp_N); -} - -void Domain::write( - char* prefix, double a, bool empty, int format, double mu, - double TAU, double U, bool original, double wo_acceleration, - double polarity, bool WLJ, bool symmetric, bool widom, double x -) { - - if(this->nanotube) this->writeNanotube( - prefix, a, empty, format, mu, TAU, U, original, - wo_acceleration, polarity, WLJ, symmetric, widom, x - ); - else this->writeGraphite( - prefix, a, empty, format, mu, TAU, U, original, - wo_acceleration, polarity, WLJ, symmetric, widom, x - ); -} - -void Domain::writeGraphite( - char* prefix, double a, bool empty, int format, double mu, - double TAU, double U, bool original, double wo_acceleration, - double polarity, bool WLJ, bool symmetric, bool widom, double x -) { - std::ofstream xdr, txt, buchholz; - std::stringstream strstrm, txtstrstrm, buchholzstrstrm; - if(format == FORMAT_BRANCH) - { - strstrm << prefix << ".xdr"; - } - if(format == FORMAT_BUCHHOLZ) - { - strstrm << prefix << ".inp"; - } - xdr.open(strstrm.str().c_str(), std::ios::trunc); - if(format == FORMAT_BRANCH) - { - txtstrstrm << prefix << "_1R.txt"; - } - if(format == FORMAT_BUCHHOLZ) - { - txtstrstrm << prefix << "_1R.cfg"; - } - txt.open(txtstrstrm.str().c_str(), std::ios::trunc); - if(format == FORMAT_BUCHHOLZ) - { - buchholzstrstrm << prefix << "_1R.xml"; - buchholz.open(buchholzstrstrm.str().c_str(), std::ios::trunc); - - /* - * Gesamter Inhalt der Buchholz-Datei - */ - buchholz << "\n\n \n " - << prefix << "_1R.cfg\n \n"; - } - - Random* r = new Random(); - r->init( - d + (int)(1000000.0*eff[2]) + (int)(10000.0*U) - + (int)(100.0*off[1]) - (int)(10000.0*TAU) + (int)(31623.0*a) - - (int)(316.23*box[1]) - ); - double REFTIME = SIG_REF * sqrt(REFMASS / EPS_REF); - double VEL_REF = SIG_REF / REFTIME; - std::cout << "Velocity unit 1 = " << VEL_REF << " * 1620.34 m/s = " - << 1620.34 * VEL_REF << " m/s.\n"; - double ACC_REF = VEL_REF / REFTIME; - double REFCARG = sqrt(EPS_REF * SIG_REF); - std::cout << "Charge unit 1 = " << REFCARG << " e.\n"; - double QDR_REF = SIG_REF*SIG_REF * REFCARG; - double REFOMGA = 1.0 / REFTIME; - - int repl = (this->flow == FLOW_COUETTE)? 2: 1; - bool tswap; - double pswap, N_id; - - bool fill[fl_units[0]][fl_units[1]][fl_units[2]][repl][3]; - for(unsigned i=0; i < fl_units[0]; i++) - for(unsigned j=0; j < fl_units[1]; j++) - for(unsigned k=0; k < fl_units[2]; k++) - for(int l=0; l < repl; l++) - for(int d=0; d < 3; d++) - fill[i][j][k][l][d] = true; - unsigned N1a = 3*repl * fl_units[0] * fl_units[1] * fl_units[2]; - if(!empty) - { - double slots = N1a; - N_id = pfill * slots; - for(int m=0; m < PRECISION; m++) - { - tswap = (N1a < N_id); - pswap = (N_id - (double)N1a) / ((tswap? slots: 0.0) - (double)N1a); - std::cout << "(N_id = " << N_id << ", N1a = " << N1a << ", tswap = " << tswap << ", pswap = " << pswap << ")\n"; - for(unsigned i=0; i < fl_units[0]; i++) - for(unsigned j=0; j < fl_units[1]; j++) - for(unsigned k=0; k < fl_units[2]; k++) - for(int l=0; l < repl; l++) - for(int d=0; d < 3; d++) - if(pswap >= r->rnd()) - { - if(fill[i][j][k][l][d]) N1a--; - fill[i][j][k][l][d] = tswap; - if(tswap) N1a++; - } - } - std::cout << "Filling " << N1a << " of " << repl << " x 3*" - << fl_units[0] << "*" << fl_units[1] << "*" << fl_units[2] - << " = " << slots << " slots (ideally " << N_id << ").\n\n"; - } - else N1a = 0; - - bool fill_ext[fl_units_ext[0]][fl_units_ext[1]][fl_units_ext[2]][repl][3]; - for(unsigned i=0; i < fl_units_ext[0]; i++) - for(unsigned j=0; j < fl_units_ext[1]; j++) - for(unsigned k=0; k < fl_units_ext[2]; k++) - for(int l=0; l < repl; l++) - for(int d=0; d < 3; d++) - fill_ext[i][j][k][l][d] = true; - unsigned N1b = 3*repl * fl_units_ext[0] * fl_units_ext[1] * fl_units_ext[2]; - if(do_fill_ext && !empty) - { - double slots = N1b; - N_id = pfill_ext * slots; - for(int m=0; m < PRECISION; m++) - { - tswap = (N1b < N_id); - pswap = (N_id - (double)N1b) / ((tswap? slots: 0.0) - (double)N1b); - std::cout << "(N_id = " << N_id << ", N1b = " << N1b << ", tswap = " << tswap << ", pswap = " << pswap << ")\n"; - for(unsigned i=0; i < fl_units_ext[0]; i++) - for(unsigned j=0; j < fl_units_ext[1]; j++) - for(unsigned k=0; k < fl_units_ext[2]; k++) - for(int l=0; l < repl; l++) - for(int d=0; d < 3; d++) - if(pswap >= r->rnd()) - { - if(fill_ext[i][j][k][l][d]) N1b--; - fill_ext[i][j][k][l][d] = tswap; - if(tswap) N1b++; - } - } - std::cout << "Filling " << N1b << " of " << repl << " x 3*" - << fl_units_ext[0] << "*" << fl_units_ext[1] << "*" << fl_units_ext[2] - << " = " << slots << " extension slots (ideally " << N_id << ").\n\n"; - } - else N1b = 0; - - unsigned N1 = N1a + N1b; - bool ignore_first; - if(symmetric && (N1 % 2)) - { - ignore_first = true; - N1a--; - N1--; - } - else ignore_first = false; - - Graphit gra; - unsigned Ngraphene = 0; - unsigned Ntotal = N1; - if(wo_wall < 1.0) - { - gra.calculateCoordinatesOfAtoms( - 1, this->box[0], this->box[2], this->bondlength, this->wo_wall - ); - Ngraphene = gra.getNumberOfAtoms(); - std::cout << "Inserting " << repl*d << " x " << Ngraphene - << " carbon atoms.\n"; - Ntotal += repl * d * Ngraphene; - } - - double LJ_CUTOFF; - double EL_CUTOFF; - double FLUIDMASS, EPS_FLUID, SIG_FLUID, FLUIDLONG, QDR_FLUID; - if(fluid == FLUID_AR) - { - FLUIDMASS = ARMASS; - EPS_FLUID = EPS_AR; - SIG_FLUID = SIG_AR; - LJ_CUTOFF = 2.5*SIG_FLUID; - } - else if(fluid == FLUID_CH4) - { - FLUIDMASS = CH4MASS; - EPS_FLUID = EPS_CH4; - SIG_FLUID = SIG_CH4; - LJ_CUTOFF = 2.5*SIG_FLUID; - } - else if(fluid == FLUID_C2H6) - { - FLUIDMASS = C2H6MASS; - EPS_FLUID = EPS_C2H6; - SIG_FLUID = SIG_C2H6; - FLUIDLONG = C2H6LONG; - QDR_FLUID = QDR_C2H6; - LJ_CUTOFF = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid == FLUID_N2) - { - FLUIDMASS = N2MASS; - EPS_FLUID = EPS_N2; - SIG_FLUID = SIG_N2; - FLUIDLONG = N2LONG; - QDR_FLUID = QDR_N2; - LJ_CUTOFF = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid == FLUID_CO2) - { - FLUIDMASS = CO2MASS; - EPS_FLUID = EPS_CO2; - SIG_FLUID = SIG_CO2; - FLUIDLONG = CO2LONG; - QDR_FLUID = QDR_CO2; - LJ_CUTOFF = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid == FLUID_AVE) - { - FLUIDMASS = AVEMASS; - EPS_FLUID = EPS_AVE; - SIG_FLUID = SIG_AVE; - LJ_CUTOFF = 2.5*SIG_FLUID; - } - else if(fluid == FLUID_H2O) - { - FLUIDMASS = OH2OMASS + 2.0*HH2OMASS; - EPS_FLUID = EPS_OH2O; - SIG_FLUID = SIG_OH2O; - FLUIDLONG = H2O_LONG; - LJ_CUTOFF = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid == FLUID_CH3OH) - { - FLUIDMASS = CCH3OHMASS + OCH3OHMASS + HCH3OHMASS; - FLUIDLONG = CH3OH_LONG; - EPS_FLUID = EPS_WANG; // used for the wall centres only - SIG_FLUID = SIG_WANG; // used for the wall centres only - LJ_CUTOFF = 4.0*SIG_CCH3OH + 0.5*FLUIDLONG; - } - else if(fluid == FLUID_C6H14) - { - FLUIDMASS = 2.0*FC6H14MASS + 4.0*MC6H14MASS; - FLUIDLONG = C6H14_LONG; - EPS_FLUID = EPS_WANG; // used for the wall centres only - SIG_FLUID = SIG_WANG; // used for the wall centres only - LJ_CUTOFF = 4.0*SIG_MC6H14 + 0.5*FLUIDLONG; - } - else - { - std::cout << "Unavailable fluid ID " << fluid << ".\n"; - exit(20); - } - - if((fluid == FLUID_AR) || (fluid == FLUID_CH4) || (fluid == FLUID_AVE)) - { - EL_CUTOFF = LJ_CUTOFF; - } - else if((fluid == FLUID_C2H6) || (fluid == FLUID_N2) || (fluid == FLUID_CO2)) - { - EL_CUTOFF = 1.2*LJ_CUTOFF; - } - else EL_CUTOFF = 1.44*LJ_CUTOFF; - - double LJ_CUTOFF2 = 0.0; - double EL_CUTOFF2 = 0.0; - double FLUIDMASS2, EPS_FLUID2, SIG_FLUID2, FLUIDLONG2, QDR_FLUID2; - if(fluid2 == FLUID_AR) - { - FLUIDMASS2 = ARMASS; - EPS_FLUID2 = EPS_AR; - SIG_FLUID2 = SIG_AR; - LJ_CUTOFF2 = 2.5*SIG_FLUID; - } - else if(fluid2 == FLUID_CH4) - { - FLUIDMASS2 = CH4MASS; - EPS_FLUID2 = EPS_CH4; - SIG_FLUID2 = SIG_CH4; - LJ_CUTOFF2 = 2.5*SIG_FLUID; - } - else if(fluid2 == FLUID_C2H6) - { - FLUIDMASS2 = C2H6MASS; - EPS_FLUID2 = EPS_C2H6; - SIG_FLUID2 = SIG_C2H6; - FLUIDLONG2 = C2H6LONG; - QDR_FLUID2 = QDR_C2H6; - LJ_CUTOFF2 = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid2 == FLUID_N2) - { - FLUIDMASS2 = N2MASS; - EPS_FLUID2 = EPS_N2; - SIG_FLUID2 = SIG_N2; - FLUIDLONG2 = N2LONG; - QDR_FLUID2 = QDR_N2; - LJ_CUTOFF2 = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid2 == FLUID_CO2) - { - FLUIDMASS2 = CO2MASS; - EPS_FLUID2 = EPS_CO2; - SIG_FLUID2 = SIG_CO2; - FLUIDLONG2 = CO2LONG; - QDR_FLUID2 = QDR_CO2; - LJ_CUTOFF2 = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid2 == FLUID_AVE) - { - FLUIDMASS2 = AVEMASS; - EPS_FLUID2 = EPS_AVE; - SIG_FLUID2 = SIG_AVE; - LJ_CUTOFF2 = 2.5*SIG_FLUID; - } - else if(fluid2 == FLUID_H2O) - { - FLUIDMASS2 = OH2OMASS + 2.0*HH2OMASS; - EPS_FLUID2 = EPS_OH2O; - SIG_FLUID2 = SIG_OH2O; - FLUIDLONG2 = H2O_LONG; - LJ_CUTOFF2 = 4.0*SIG_FLUID + 0.5*FLUIDLONG; - } - else if(fluid2 == FLUID_CH3OH) - { - FLUIDMASS2 = CCH3OHMASS + OCH3OHMASS + HCH3OHMASS; - FLUIDLONG2 = CH3OH_LONG; - EPS_FLUID2 = EPS_WANG; // used for the wall centres only - SIG_FLUID2 = SIG_WANG; // used for the wall centres only - LJ_CUTOFF2 = 4.0*SIG_CCH3OH + 0.5*FLUIDLONG; - } - else if(fluid2 == FLUID_C6H14) - { - FLUIDMASS2 = 2.0*FC6H14MASS + 4.0*MC6H14MASS; - FLUIDLONG2 = C6H14_LONG; - EPS_FLUID2 = EPS_WANG; // used for the wall centres only - SIG_FLUID2 = SIG_WANG; // used for the wall centres only - LJ_CUTOFF2 = 4.0*SIG_MC6H14 + 0.5*FLUIDLONG; - } - - if((fluid2 == FLUID_AR) || (fluid2 == FLUID_CH4) || (fluid2 == FLUID_AVE)) - { - EL_CUTOFF2 = LJ_CUTOFF2; - } - else if((fluid2 == FLUID_C2H6) || (fluid2 == FLUID_N2) || (fluid2 == FLUID_CO2)) - { - EL_CUTOFF2 = 1.2*LJ_CUTOFF2; - } - else if(fluid2 != FLUID_NIL) EL_CUTOFF2 = 1.44*LJ_CUTOFF2; - - if(LJ_CUTOFF2 > LJ_CUTOFF) LJ_CUTOFF = LJ_CUTOFF2; - if(EL_CUTOFF2 > EL_CUTOFF) EL_CUTOFF = EL_CUTOFF2; - - unsigned wallcomp_sng = (wo_wall == 1.0)? 0: ((polarity == 0.0)? d: NCOMP_POLAR); - unsigned wallcomp = repl*wallcomp_sng; - unsigned fluidcomp = (symmetric || (fluid2 != FLUID_NIL))? 2: 1; - - xdr.precision(7); - if(format == FORMAT_BRANCH) - { - xdr << "mardyn " << TIME << " tersoff\n" - << "# mardyn input file, ls1 project\n" - << "# generated by the mkcp tool\n"; - } - if(format == FORMAT_BUCHHOLZ) - { - xdr << "mardyn trunk " << TIME << "\n"; - } - - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << "t\t0.0\n" - << "L" << "\t" << box[0]/SIG_REF << "\t" - << repl*box[1]/SIG_REF << "\t" << box[2]/SIG_REF - << "\n"; - } - - unsigned wallthermostats = (polarity == 0.0)? wallcomp: repl; - if(format == FORMAT_BUCHHOLZ) - { - if(wallcomp == 0) - { - xdr << "T\t" << T/EPS_REF << "\n"; - } - else - { - for(unsigned i = fluidcomp + 1; wallcomp + fluidcomp >= i; i++) - { - xdr << "CT\t" << i << " " << ((polarity == 0.0)? i-fluidcomp: ((wallcomp_sng+fluidcomp >= i)? 1: 2)) << "\n"; - } - for(unsigned i=1; wallthermostats >= i; i++) - { - xdr << "ThT\t" << i << " " << T/EPS_REF << "\n"; - if(flow == FLOW_COUETTE) xdr << "U\t" << i << "\n"; - } - for(unsigned i=1; fluidcomp >= i; i++) - xdr << "CT\t" << i << " " << wallthermostats + 1 << "\n"; - xdr << "ThT\t" << wallthermostats + 1 << " " << T/EPS_REF << "\n"; - } - if(flow == FLOW_COUETTE) - { - if(polarity == 0.0) - { - xdr << "S\t" << fluidcomp + 1 << " 1\nS\t" << fluidcomp + this->d << " 2\n" - << "A\t1\t0.0 0.0 " << U/VEL_REF << "\t" << TAU/REFTIME - << "\t0.0 0.0 " << a/ACC_REF << "\n" - << "A\t2\t0.0 0.0 " << U/VEL_REF << "\t" << TAU/REFTIME - << "\t0.0 0.0 " << a/ACC_REF << "\n" - << "S\t" << fluidcomp + d + 1 << " 3\n" - << "S\t" << fluidcomp + 2*d << " 4\n" - << "A\t3\t0.0 0.0 0.0\t" << TAU_ZERO/REFTIME - << "\t0.0 0.0 " << -1.0*a/ACC_REF << "\n" - << "A\t4\t0.0 0.0 0.0\t" << TAU_ZERO/REFTIME - << "\t0.0 0.0 " << -1.0*a/ACC_REF << "\n"; - for(unsigned i = fluidcomp + 2; i < fluidcomp + d; i++) - xdr << "S\t" << i << " 5\nS\t" << d+i << " 6\n"; - xdr << "A\t" << 5 << "\t0.0 0.0 " << U/VEL_REF << "\t" - << TAU/REFTIME << "\t0.0 0.0 0.0\n"; - xdr << "A\t" << 6 << "\t0.0 0.0 0.0\t" << TAU_ZERO/REFTIME - << "\t0.0 0.0 0.0\n"; - } - else - { - for(unsigned i=1; i < wallcomp_sng; i++) - { - xdr << "S\t" << fluidcomp + i << " 1\nS\t" << wallcomp_sng + fluidcomp + i << "\t 2\n"; - } - xdr << "A\t1\t0.0 0.0 " << U/VEL_REF << "\t" << TAU/REFTIME - << "\t0.0 0.0 " << a/ACC_REF << "\n" - << "A\t2\t0.0 0.0 0.0\t" << TAU/REFTIME - << "\t0.0 0.0 " << -1.0*a/ACC_REF << "\n"; - } - } - else if(flow == FLOW_POISEUILLE) - { - xdr << "S\t1 1\n"; - if(fluid2 != FLUID_NIL) xdr << "S\t2 1\n"; - xdr << "A\t1\t0.0 0.0 " << U/VEL_REF << "\t" - << TAU/REFTIME << "\t0.0 0.0 " << a/ACC_REF << "\n"; - if(symmetric) - xdr << "S\t2 2\nA\t2\t0.0 0.0 " << -U/VEL_REF << "\t" - << TAU/REFTIME << "\t0.0 0.0 " << -a/ACC_REF << "\n"; - } - } - - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << "C" << "\t" << wallcomp + fluidcomp << "\n"; - - int tfluid = fluid; - double TFLUIDMASS = FLUIDMASS; - double TEPS_FLUID = EPS_FLUID; - double TSIG_FLUID = SIG_FLUID; - double TFLUIDLONG = FLUIDLONG; - double TQDR_FLUID = QDR_FLUID; - for(unsigned flit = 0; flit < fluidcomp; flit++) - { - if((tfluid == FLUID_AR) || (tfluid == FLUID_CH4) || (tfluid == FLUID_AVE)) - { - xdr << "1 0 0 0 0\n" // LJ, C, Q, D, Tersoff - << "0.0 0.0 0.0\t" - << TFLUIDMASS/REFMASS << " " << TEPS_FLUID/EPS_REF << " " - << TSIG_FLUID/SIG_REF << "\t0.0 0.0 0.0"; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 1"; - xdr << "\n"; - } - else if((tfluid == FLUID_C2H6) || (tfluid == FLUID_N2) || (tfluid == FLUID_CO2)) - { - xdr << "2 0 1 0 0\n" // LJ, C, Q, D, Tersoff - << "0.0 0.0 " << -0.5*TFLUIDLONG/SIG_REF << "\t" - << 0.5*TFLUIDMASS/REFMASS << " " << TEPS_FLUID/EPS_REF - << " " << TSIG_FLUID/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n" - << "0.0 0.0 " << +0.5*TFLUIDLONG/SIG_REF << "\t" - << 0.5*TFLUIDMASS/REFMASS << " " << TEPS_FLUID/EPS_REF - << " " << TSIG_FLUID/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n" - << "0.0 0.0 0.0\t0.0 0.0 1.0\t" << TQDR_FLUID/QDR_REF - << "\n0.0 0.0 0.0\n"; - } - else if(tfluid == FLUID_H2O) - { - xdr << "1 3 0 0 0\n"; // LJ, C, Q, D, Tersoff - - xdr << R0_O_H2O/SIG_REF << " " << R1_O_H2O/SIG_REF << " " << R2_O_H2O/SIG_REF << "\t" - << OH2OMASS/REFMASS << " " << EPS_OH2O/EPS_REF << " " << SIG_OH2O/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << R0_H1H2O/SIG_REF << " " << R1_H1H2O/SIG_REF << " " << R2_H1H2O/SIG_REF << "\t" - << HH2OMASS/REFMASS << " " << CHG_HH2O/REFCARG << "\n"; - xdr << R0_H2H2O/SIG_REF << " " << R1_H2H2O/SIG_REF << " " << R2_H2H2O/SIG_REF << "\t" - << HH2OMASS/REFMASS << " " << CHG_HH2O/REFCARG << "\n"; - xdr << R0_E_H2O/SIG_REF << " " << R1_E_H2O/SIG_REF << " " << R2_E_H2O/SIG_REF << "\t" - << "0.0 " << CHG_EH2O/REFCARG << "\n"; - - xdr << "0.0 0.0 0.0\n"; - } - else if(tfluid == FLUID_CH3OH) - { - xdr << "2 3 0 0 0\n"; // LJ, C, Q, D, Tersoff - - xdr << R0_C_CH3OH/SIG_REF << " " << R1_C_CH3OH/SIG_REF << " " << R2_C_CH3OH/SIG_REF << "\t" - << CCH3OHMASS/REFMASS << " " << EPS_CCH3OH/EPS_REF << " " << SIG_CCH3OH/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << R0_O_CH3OH/SIG_REF << " " << R1_O_CH3OH/SIG_REF << " " << R2_O_CH3OH/SIG_REF << "\t" - << OCH3OHMASS/REFMASS << " " << EPS_OCH3OH/EPS_REF << " " << SIG_OCH3OH/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << R0_C_CH3OH/SIG_REF << " " << R1_C_CH3OH/SIG_REF << " " << R2_C_CH3OH/SIG_REF << "\t" - << "0.0 " << CHG_CCH3OH/REFCARG << "\n"; - xdr << R0_O_CH3OH/SIG_REF << " " << R1_O_CH3OH/SIG_REF << " " << R2_O_CH3OH/SIG_REF << "\t" - << "0.0 " << CHG_OCH3OH/REFCARG << "\n"; - xdr << R0_H_CH3OH/SIG_REF << " " << R1_H_CH3OH/SIG_REF << " " << R2_H_CH3OH/SIG_REF << "\t" - << HCH3OHMASS/REFMASS << " " << CHG_HCH3OH/REFCARG << "\n"; - - xdr << "0.0 0.0 0.0\n"; - } - else if(tfluid == FLUID_C6H14) - { - xdr << "6 0 0 0 0\n"; // LJ, C, Q, D, Tersoff - - xdr << R0_F1C6H14/SIG_REF << " " << R1_F1C6H14/SIG_REF << " " << R2_F1C6H14/SIG_REF << "\t" - << FC6H14MASS/REFMASS << " " << EPS_FC6H14/EPS_REF << " " << SIG_FC6H14/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << R0_M1C6H14/SIG_REF << " " << R1_M1C6H14/SIG_REF << " " << R2_M1C6H14/SIG_REF << "\t" - << MC6H14MASS/REFMASS << " " << EPS_MC6H14/EPS_REF << " " << SIG_MC6H14/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << R0_M2C6H14/SIG_REF << " " << R1_M2C6H14/SIG_REF << " " << R2_M2C6H14/SIG_REF << "\t" - << MC6H14MASS/REFMASS << " " << EPS_MC6H14/EPS_REF << " " << SIG_MC6H14/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << R0_M3C6H14/SIG_REF << " " << R1_M3C6H14/SIG_REF << " " << R2_M3C6H14/SIG_REF << "\t" - << MC6H14MASS/REFMASS << " " << EPS_MC6H14/EPS_REF << " " << SIG_MC6H14/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << R0_M4C6H14/SIG_REF << " " << R1_M4C6H14/SIG_REF << " " << R2_M4C6H14/SIG_REF << "\t" - << MC6H14MASS/REFMASS << " " << EPS_MC6H14/EPS_REF << " " << SIG_MC6H14/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - xdr << R0_F2C6H14/SIG_REF << " " << R1_F2C6H14/SIG_REF << " " << R2_F2C6H14/SIG_REF << "\t" - << FC6H14MASS/REFMASS << " " << EPS_FC6H14/EPS_REF << " " << SIG_FC6H14/SIG_REF; - if(format == FORMAT_BUCHHOLZ) xdr << "\t" << LJ_CUTOFF << " 0"; - xdr << "\n"; - - xdr << "0.0 0.0 0.0\n"; - } - else - { - std::cout << "Fluid code " << tfluid << ": Not yet implemented.\n"; - exit(1000+tfluid); - } - - if(!symmetric) - { - tfluid = fluid2; - TFLUIDMASS = FLUIDMASS2; - TEPS_FLUID = EPS_FLUID2; - TSIG_FLUID = SIG_FLUID2; - TFLUIDLONG = FLUIDLONG2; - TQDR_FLUID = QDR_FLUID2; - } - } - - double crga[repl*NCOMP_POLAR]; - for(unsigned i=0; i < wallcomp; i++) crga[i] = 0.0; - for(int r=0; r < repl; r++) - { - crga[1 + r*NCOMP_POLAR] = polarity/3.0; - crga[2 + r*NCOMP_POLAR] = -2.0*polarity/3.0; - crga[3 + r*NCOMP_POLAR] = -1.0*polarity/3.0; - } - - for(unsigned i=0; i < wallcomp; i++) - { - if(polarity == 0.0) xdr << "1 0 0 0 1\n"; // LJ, C, Q, D, Tersoff - else if(crga[i] == 0.0) xdr << "1 0 0 0 1\n"; - else xdr << "1 1 0 0 1\n"; - - if(polarity == 0.0) - { - xdr << "0.0 0.0 0.0\t" << 0.5*ATOMIC_MASS_C/REFMASS << " " - << EPS_FLUID/EPS_REF << " " << SIG_FLUID/SIG_REF; - if(format == FORMAT_BUCHHOLZ) - xdr << "\t" << LJ_CUTOFF - << (((fluid == FLUID_AR) || (fluid == FLUID_CH4) || (fluid == FLUID_AVE))? " 1" - : " 0"); - xdr << "\n"; - } - else if(crga[i] != 0.0) - { - xdr << "0.0 0.0 0.0\t0.0 " - << EPS_FLUID/EPS_REF << " " << SIG_FLUID/SIG_REF; - if(format == FORMAT_BUCHHOLZ) - xdr << "\t" << LJ_CUTOFF - << (((fluid == FLUID_AR) || (fluid == FLUID_CH4) || (fluid == FLUID_AVE))? " 1" - : " 0"); - xdr << "\n" - << "0.0 0.0 0.0\t" << 0.5*ATOMIC_MASS_C/REFMASS - << " " << crga[i]/REFCARG << "\n"; - } - else - { - xdr << "0.0 0.0 0.0\t0.0 " - << EPS_FLUID/EPS_REF << " " << SIG_FLUID/SIG_REF; - if(format == FORMAT_BUCHHOLZ) - xdr << "\t" << LJ_CUTOFF - << (((fluid == FLUID_AR) || (fluid == FLUID_CH4) || (fluid == FLUID_AVE))? " 1" - : " 0"); - xdr << "\n"; - } - - if((polarity != 0.0) && (crga[i] == 0.0)) - xdr << "0.0 0.0 0.0\t" << ATOMIC_MASS_C/REFMASS << " "; - else - xdr << "0.0 0.0 0.0\t" << 0.5*ATOMIC_MASS_C/REFMASS << " "; - - xdr << TERSOFF_A/EPS_REF << " " << TERSOFF_B/EPS_REF << "\t" - << (original? TERSOFF_LAMBDA_ORIG - : TERSOFF_LAMBDA) * SIG_REF << " " - << (original? TERSOFF_MU_ORIG: TERSOFF_MU) * SIG_REF - << " " << (original? TERSOFF_R_ORIG: TERSOFF_R)/SIG_REF - << " " << (original? TERSOFF_S_ORIG: TERSOFF_S)/SIG_REF - << "\t" << TERSOFF_C << " " << TERSOFF_D << " " - << TERSOFF_H << " " << TERSOFF_N << " " << TERSOFF_BETA - << "\n0.0 0.0 0.0\n"; - } - if(symmetric) - { - xdr << "1.0 1.0 "; - for(unsigned j=3; wallcomp + 2 >= j; j++) - xdr << this->XI << " " << this->ETA << " "; - xdr << "\n"; - for(unsigned j=3; wallcomp + 2 >= j; j++) - xdr << this->XI << " " << this->ETA << " "; - xdr << "\n"; - } - else if(fluid2 != FLUID_NIL) - { - xdr << this->XIF << " " << this->ETAF << " "; - for(unsigned j=3; wallcomp + 2 >= j; j++) - xdr << this->XI << " " << this->ETA << " "; - xdr << "\n"; - for(unsigned j=3; wallcomp + 2 >= j; j++) - xdr << this->XI2 << " " << this->ETA2 << " "; - xdr << "\n"; - } - else - { - for(unsigned j=2; wallcomp + 1 >= j; j++) - xdr << this->XI << " " << this->ETA << " "; - xdr << "\n"; - } - for(unsigned i=2; wallcomp >= i; i++) - { - for(unsigned j = i+1; wallcomp + 1 >= j; j++) - { - if(WLJ) xdr << "1.0 1.0 "; - else xdr << "0.0 1.0 "; - } - xdr << "\n"; - } - xdr << "1.0e+10\n"; - } - - if(format == FORMAT_BRANCH) - { - if(wallcomp == 0) - { - xdr << "T\t" << T/EPS_REF << "\n"; - } - else - { - for(unsigned i = fluidcomp + 1; wallcomp + fluidcomp >= i; i++) - { - xdr << "CT\t" << i << " " << ((polarity == 0.0)? i-fluidcomp: ((wallcomp_sng+fluidcomp >= i)? 1: 2)) << "\n"; - } - for(unsigned i=1; wallthermostats >= i; i++) - { - xdr << "ThT\t" << i << " " << T/EPS_REF << "\n"; - if(flow == FLOW_COUETTE) xdr << "U\t" << i << "\n"; - } - for(unsigned i=1; fluidcomp >= i; i++) - xdr << "CT\t" << i << " " << wallthermostats + 1 << "\n"; - xdr << "ThT\t" << wallthermostats + 1 << " " << T/EPS_REF << "\n"; - } - if(flow == FLOW_COUETTE) - { - if(polarity == 0.0) - { - xdr << "S\t" << fluidcomp + 1 << " 1\nS\t" << fluidcomp + this->d << " 2\n" - << "A\t1\t0.0 0.0 " << U/VEL_REF << "\t" << TAU/REFTIME - << "\t0.0 0.0 " << a/ACC_REF << "\n" - << "A\t2\t0.0 0.0 " << U/VEL_REF << "\t" << TAU/REFTIME - << "\t0.0 0.0 " << a/ACC_REF << "\n" - << "S\t" << fluidcomp + d + 1 << " 3\n" - << "S\t" << fluidcomp + 2*d << " 4\n" - << "A\t3\t0.0 0.0 0.0\t" << TAU_ZERO/REFTIME - << "\t0.0 0.0 " << -1.0*a/ACC_REF << "\n" - << "A\t4\t0.0 0.0 0.0\t" << TAU_ZERO/REFTIME - << "\t0.0 0.0 " << -1.0*a/ACC_REF << "\n"; - for(unsigned i = fluidcomp + 2; i < fluidcomp + d; i++) - xdr << "S\t" << i << " 5\nS\t" << d+i << " 6\n"; - xdr << "A\t" << 5 << "\t0.0 0.0 " << U/VEL_REF << "\t" - << TAU/REFTIME << "\t0.0 0.0 0.0\n"; - xdr << "A\t" << 6 << "\t0.0 0.0 0.0\t" << TAU_ZERO/REFTIME - << "\t0.0 0.0 0.0\n"; - } - else - { - for(unsigned i=1; i < wallcomp_sng; i++) - { - xdr << "S\t" << fluidcomp + i << " 1\nS\t" << wallcomp_sng + fluidcomp + i << "\t 2\n"; - } - xdr << "A\t1\t0.0 0.0 " << U/VEL_REF << "\t" << TAU/REFTIME - << "\t0.0 0.0 " << a/ACC_REF << "\n" - << "A\t2\t0.0 0.0 0.0\t" << TAU/REFTIME - << "\t0.0 0.0 " << -1.0*a/ACC_REF << "\n"; - } - } - else if(flow == FLOW_POISEUILLE) - { - xdr << "S\t1 1\n"; - xdr << "A\t1\t0.0 0.0 " << U/VEL_REF << "\t" - << TAU/REFTIME << "\t0.0 0.0 " << a/ACC_REF << "\n"; - if(symmetric) - xdr << "S\t2 2\nA\t2\t0.0 0.0 " << -U/VEL_REF << "\t" - << TAU/REFTIME << "\t0.0 0.0 " << -a/ACC_REF << "\n"; - else if(fluid2 != FLUID_NIL) - xdr << "S\t2 2\nA\t2\t0.0 0.0 " << U/VEL_REF << "\t" - << TAU/REFTIME << "\t0.0 0.0 " << a/ACC_REF << "\n"; - } - } - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << "N" << "\t" << Ntotal << "\nM" << "\t" << "ICRVQD\n\n"; - } - - txt.precision(6); - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt << "mardynconfig\ntimestepLength\t" << DT/REFTIME - << "\ncutoffRadius\t" << EL_CUTOFF/SIG_REF - << "\nLJCutoffRadius\t" << LJ_CUTOFF/SIG_REF << "\n" - << "tersoffCutoffRadius\t" - << 1.0001*(original? TERSOFF_S_ORIG: TERSOFF_S) / SIG_REF << "\n" - << "initCanonical\t20001\n"; - if(muVT || widom) - { - txt.precision(9); - txt << "chemicalPotential " << (muVT? (mu/EPS_REF): 0.0) - << " component 1 control 0.0 0.0 0.0 to " - << this->box[0]/SIG_REF << " " << 0.5*this->h/SIG_REF - << " " << this->box[2]/SIG_REF << " conduct " - << 1 + (int)round(((flow == FLOW_COUETTE)? 0.0005 - : 0.001) * N1) - << " tests every 2 steps\n"; - if(symmetric && widom && !muVT) - { - txt << "chemicalPotential " << 0.0 - << " component 2 control 0.0 0.0 0.0 to " - << this->box[0]/SIG_REF << " " << 0.5*this->h/SIG_REF - << " " << this->box[2]/SIG_REF << " conduct " - << 1 + (int)round(((flow == FLOW_COUETTE)? 0.0005 - : 0.001) * N1) - << " tests every 2 steps\n"; - } - if(flow == FLOW_COUETTE) - { - txt << "chemicalPotential " << (muVT? (mu/EPS_REF): 0.0) - << " component 1 control 0.0 " << this->box[1]/SIG_REF - << " 0.0 to " << this->box[0]/SIG_REF << " " - << (this->box[1] + 0.5*this->h)/SIG_REF << " " - << box[2]/SIG_REF << " conduct " - << 1 + (int)round(((flow == FLOW_COUETTE)? 0.0005 - : 0.001) * N1) - << " tests every 2 steps\n"; - if(symmetric && widom && !muVT) - { - txt << "chemicalPotential " << 0.0 - << " component 2 control 0.0 " << this->box[1]/SIG_REF - << " 0.0 to " << this->box[0]/SIG_REF << " " - << (this->box[1] + 0.5*this->h)/SIG_REF << " " - << box[2]/SIG_REF << " conduct " - << 1 + (int)round(((flow == FLOW_COUETTE)? 0.0005 - : 0.001) * N1) - << " tests every 2 steps\n"; - } - } - txt << "planckConstant\t" - << sqrt(6.28319 * T/EPS_REF) << "\n"; // sqrt(2 pi kT) - txt << "initGrandCanonical\t" << (muVT? 40001: 60001) << "\n"; - } - if(widom) txt << "Widom\n"; - txt << "initStatistics\t60001\n"; - if(flow != FLOW_NONE) - { - txt << "constantAccelerationTimesteps\t32\n" - << "zetaFlow\t" << ZETA*sqrt(REFTIME) << "\n" - << "tauPrimeFlow\t" << TAUPRIME/REFTIME << "\n"; - } - } - - txt.precision(5); - if(format == FORMAT_BRANCH) - { - txt << "phaseSpaceFile\t" << prefix - << ".xdr\n# for LinkedCells, the cellsInCutoffRadius has to" - << " be provided\ndatastructure\tLinkedCells\t1\noutput\t"; - } - if(format == FORMAT_BUCHHOLZ) - { - txt << "phaseSpaceFile\tOldStyle\t" << prefix << ".inp\nparallelization\tDomainDecomposition\n" - << "datastructure\tLinkedCells\t1\noutput\t"; - } - - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt << "ResultWriter\t40\t" << prefix - << "_1R\noutput\tXyzWriter\t8000\t" << prefix - << "_1R.buxyz\noutput\tVisittWriter\t10000000\t" << prefix - << "_1R\nprofile\t1 362 362\nprofileRecordingTimesteps\t1\n" - << "profileOutputTimesteps\t60000" - << "\nprofiledComponent\t1\n"; - if(fluidcomp >= 2) txt << "profiledComponent\t2\n"; - txt << "profileOutputPrefix\t" << prefix << "_1R\n"; - if(wo_wall < 1.0) - { - if(WLJ) - { - txt << "wallLJ\ton\n"; - if(flow == FLOW_COUETTE) - { - txt << "zOscillator 1536\n"; - } - else - { - txt << "oscillator 1536\n"; - } - } - else - { - txt << "wallLJ\toff\n"; - if(flow == FLOW_COUETTE) - { - txt << "zOscillator 384\n"; - } - else - { - txt << "oscillator 384\n"; - } - } - } - if(symmetric || (flow == FLOW_NONE)) - txt << "nomomentum 1024\n"; - if(wo_acceleration != 0.0) - { - txt << "flowControl\t0 0 " << 0.5*wo_acceleration*box[2]/SIG_REF - << " to " << box[0]/SIG_REF << " " << repl*box[1]/SIG_REF << " " - << (1.0 - 0.5*wo_acceleration)*box[2]/SIG_REF << "\n"; - } - } - - double I[3]; - for(int k=0; k < 3; k++) I[k] = 0.0; - if((fluid == FLUID_C2H6) || (fluid == FLUID_N2) || (fluid == FLUID_CO2)) - { - I[0] = 0.25 * FLUIDMASS * FLUIDLONG * FLUIDLONG; - I[1] = I[0]; - } - else if(fluid == FLUID_H2O) - { - I[0] = I_XX_H2O; - I[1] = I_YY_H2O; - I[2] = I_ZZ_H2O; - } - else if(fluid == FLUID_CH3OH) - { - I[0] = I_XX_CH3OH; - I[1] = I_YY_CH3OH; - I[2] = I_ZZ_CH3OH; - } - else if(fluid == FLUID_C6H14) - { - I[0] = I_XX_C6H14; - I[1] = I_YY_C6H14; - I[2] = I_ZZ_C6H14; - } - - double I2[3]; - if(symmetric) - { - for(int k=0; k < 3; k++) I2[k] = I[k]; - } - else - { - for(int k=0; k < 3; k++) I2[k] = 0.0; - if((fluid2 == FLUID_C2H6) || (fluid2 == FLUID_N2) || (fluid2 == FLUID_CO2)) - { - I2[0] = 0.25 * FLUIDMASS * FLUIDLONG * FLUIDLONG; - I2[1] = I2[0]; - } - else if(fluid2 == FLUID_H2O) - { - I2[0] = I_XX_H2O; - I2[1] = I_YY_H2O; - I2[2] = I_ZZ_H2O; - } - else if(fluid2 == FLUID_CH3OH) - { - I2[0] = I_XX_CH3OH; - I2[1] = I_YY_CH3OH; - I2[2] = I_ZZ_CH3OH; - } - else if(fluid2 == FLUID_C6H14) - { - I2[0] = I_XX_C6H14; - I2[1] = I_YY_C6H14; - I2[2] = I_ZZ_C6H14; - } - } - - /* - * main fluid volume - */ - unsigned id = 1; - double tr[3]; - unsigned ii[3]; - double Nf[2]; - Nf[0] = (1.0 - x) * (double)N1 + 0.00001; - Nf[1] = x * (double)N1 + 0.00001; - if(!empty) for(ii[0]=0; ii[0] < this->fl_units[0]; (ii[0]) ++) - for(ii[1]=0; ii[1] < this->fl_units[1]; (ii[1]) ++) - for(ii[2]=0; ii[2] < this->fl_units[2]; (ii[2]) ++) - for(int j=0; j < repl; j++) - { - for(int d=0; d < 3; d++) - { - if(fill[ ii[0] ][ ii[1] ][ ii[2] ][ j ][ d ]) - { - if(!ignore_first) - { - for(int k=0; k < 3; k++) - { - tr[k] = off[k] + fl_unit[k] * ( - ii[k] + 0.004*r->rnd() + ((k == d)? 0.248: 0.748) - ); - } - tr[1] += j*box[1]; - box[1] *= repl; - for(int k=0; k < 3; k++) - { - if(tr[k] > box[k]) tr[k] -= box[k]; - else if(tr[k] < 0.0) tr[k] += box[k]; - } - box[1] /= repl; - double tv = sqrt(3.0*T / FLUIDMASS); - double phi = 6.283185 * r->rnd(); - double omega = 6.283185 * r->rnd(); - double w[3]; - for(int k=0; k < 3; k++) - w[k] = (I[k] == 0)? 0.0: ((r->rnd() > 0.5)? 1: -1) * sqrt(2.0*r->rnd()*T / I[k]); - unsigned cid; - if(fluidcomp == 1) cid = 1; - else - { - cid = (r->rnd() > (Nf[0] / (Nf[0] + Nf[1])))? 2: 1; - --Nf[cid - 1]; - } - xdr << id << " " << cid << "\t" << tr[0]/SIG_REF - << " " << tr[1]/SIG_REF << " " << tr[2]/SIG_REF - << "\t" << tv*cos(phi)*cos(omega)/VEL_REF << " " - << tv*cos(phi)*sin(omega)/VEL_REF << " " - << tv*sin(phi)/VEL_REF << "\t1.0 0.0 0.0 0.0\t" - << w[0]/REFOMGA << " " << w[1]/REFOMGA << " " - << w[2]/REFOMGA << "\n"; - id++; - } - else ignore_first = false; - } - else xdr << "\n"; - } - } - xdr << "\n\n"; - - /* - * fluid volume: extension - */ - if(do_fill_ext && !empty) for(ii[0]=0; ii[0] < this->fl_units_ext[0]; (ii[0]) ++) - for(ii[1]=0; ii[1] < this->fl_units_ext[1]; (ii[1]) ++) - for(ii[2]=0; ii[2] < this->fl_units_ext[2]; (ii[2]) ++) - for(int j=0; j < repl; j++) - { - for(int d=0; d < 3; d++) - { - if(fill_ext[ ii[0] ][ ii[1] ][ ii[2] ][ j ][ d ]) - { - for(int k=0; k < 3; k++) - { - tr[k] = off_ext[k] + fl_unit_ext[k] * ( - ii[k] + 0.004*r->rnd() + ((k == d)? 0.248: 0.748) - ); - } - tr[1] += j*box[1]; - box[1] *= repl; - for(int k=0; k < 3; k++) - { - if(tr[k] > box[k]) tr[k] -= box[k]; - else if(tr[k] < 0.0) tr[k] += box[k]; - } - box[1] /= repl; - double tv = sqrt(3.0*T / FLUIDMASS); - double phi = 6.283185 * r->rnd(); - double omega = 6.283185 * r->rnd(); - double w[3]; - for(int k=0; k < 3; k++) - w[k] = (I[k] == 0)? 0.0: ((r->rnd() > 0.5)? 1: -1) * sqrt(2.0*r->rnd()*T / I[k]); - unsigned cid; - if(fluidcomp == 1) cid = 1; - else - { - cid = (r->rnd() > (Nf[0] / (Nf[0] + Nf[1])))? 2: 1; - --Nf[cid - 1]; - } - xdr << id << " " << cid << "\t" << tr[0]/SIG_REF - << " " << tr[1]/SIG_REF << " " << tr[2]/SIG_REF - << "\t" << tv*cos(phi)*cos(omega)/VEL_REF << " " - << tv*cos(phi)*sin(omega)/VEL_REF << " " - << tv*sin(phi)/VEL_REF << "\t1.0 0.0 0.0 0.0\t" - << w[0]/REFOMGA << " " << w[1]/REFOMGA << " " - << w[2]/REFOMGA << "\n"; - id++; - } - else xdr << "\n"; - } - } - xdr << "\n\n"; - - for(int j=0; j < repl; j++) - { - double layerU = ((flow == FLOW_COUETTE) && (j == 0)) - ? U - : 0.0; - gra.calculateVelocities(T, layerU); - for(unsigned k = fluidcomp+1; fluidcomp+d >= k; k++) - { - unsigned layerid = j*d + k; - double yoffset = 0.5*h + (k - fluidcomp - 1.0)*Z; - for(unsigned l=0; l < Ngraphene; l++) - { - tr[0] = gra.getX(l); - tr[1] = gra.getY(l) + yoffset; - tr[2] = gra.getZ(l); - unsigned tcid = gra.getComponent(l) + fluidcomp + j*NCOMP_POLAR; - for(int m=0; m < 3; m++) - { - tr[m] += (0.004*r->rnd() - 0.002) * bondlength; - } - if(j == 1) tr[1] += box[1]; - box[1] *= repl; - for(int m=0; m < 3; m++) - { - if(tr[m] > box[m]) tr[m] -= box[m]; - else if(tr[m] < 0.0) tr[m] += box[m]; - } - box[1] /= repl; - - xdr << id << " " << ((polarity == 0.0)? layerid: tcid) - << "\t" << tr[0]/SIG_REF - << " " << tr[1]/SIG_REF << " " << tr[2]/SIG_REF - << "\t" << gra.getVelocityX(l)/VEL_REF << " " - << gra.getVelocityY(l)/VEL_REF << " " - << gra.getVelocityZ(l)/VEL_REF - << "\t1.0 0.0 0.0 0.0\t0.0 0.0 0.0\n"; - - id++; - } - xdr << "\n"; - } - xdr << "\n"; - } - - xdr.close(); - txt.close(); - if(format == FORMAT_BUCHHOLZ) - { - buchholz.close(); - } -} - -void Domain::writeNanotube( - char* prefix, double a, bool empty, int format, double mu, - double TAU, double U, bool original, double wo_acceleration, - double polarity, bool WLJ, bool symmetric, bool widom, double x -) { - std::cout << "Cannot create the nanotube - implementation missing.\n"; - exit(19); -} - -void Domain::specifyGraphite(double rho, unsigned N) -{ - /* - * y direction - */ - if(wo_wall == 1.0) - { - this->box[1] = this->h; - this->eff[1] = this->h; - this->off[1] = 0.01 * this->h; - - this->shielding = 0.0; - } - else - { - this->box[1] = this->h + ((double)(this->d)-1.0)*Z; - if(1.1 * this->shielding > 0.5 * this->h) - { - std::cout << "Warning: shielding = " << shielding - << " versus h = " << h << ", "; - this->shielding = 0.5*this->h - 0.1 * this->shielding; - if(this->shielding < 0.0) this->shielding = 0.0; - std::cout << "corrected to shielding = " << shielding << ".\n"; - } - this->eff[1] = this->h - 2.0*this->shielding; - this->off[1] = 0.5*this->h + ((double)(this->d)-1.0)*Z - + this->shielding; - } - - double V_id = (double)(N/rho) / (wo_wall + (1.0 - wo_wall)*eff[1]/box[1]); - if(this->flow == FLOW_COUETTE) V_id *= 0.5; - std::cout << "Carbon-carbon bond length: " << bondlength - << " * 0.05291772 nm.\n"; - std::cout << "Total volume should approach " << V_id - << " * 1.4818e-04 nm^3.\n"; - - /* - * x and z direction - */ - double A = V_id / this->box[1]; - if(wo_wall == 1.0) - { - this->box[0] = sqrt(A); - this->box[2] = sqrt(A); - } - else - { - double tX = this->bondlength * 3.0; - double tZ = this->bondlength * 12.124356 / (1.0 - wo_wall); // effektiv 14 * sin(pi/3) - int zeta = round(sqrt(2.0*A) / tZ); // bewirkt, dass etwa = 2 wird - if(zeta == 0) zeta = 1; - int xi = round(A / (tX*tZ*zeta)); - if(xi == 0) - { - xi = 1; - zeta = round(A / (tX*tZ*xi)); - if(zeta == 0) - { - std::cout << "Warning: The generated box will be larger than specified, due to technical reasons.\n\n"; - zeta = 1; - } - } - this->box[0] = xi*tX; - this->box[2] = zeta*tZ; - } - this->eff[0] = this->box[0]; - this->off[0] = 0.0; - this->eff[2] = this->box[2]; - this->off[2] = 0.0; - - /* - * fluid unit box dimensions - */ - double V_eff = this->eff[0] * this->eff[1] * this->eff[2]; - std::cout << "Symmetry volume " << box[0]*box[1]*box[2] - << " * 1.4818e-04 nm^3, effectively " << V_eff - << " * 1.4818e-04 nm^3.\n"; - double N_id = rho*V_eff; - double N_boxes = N_id / 3.0; - fl_units[1] = round( - pow( - (N_boxes * this->eff[1] * this->eff[1]) - / (this->eff[0] * this->eff[2]), 1.0/3.0 - ) - ); - if(fl_units[1] == 0) fl_units[1] = 1; - double bxbz_id = N_boxes / fl_units[1]; - fl_units[0] = round(sqrt(this->eff[0] * bxbz_id / this->eff[2])); - if(fl_units[0] == 0) this->fl_units[0] = 1; - fl_units[2] = ceil(bxbz_id / fl_units[0]); - std::cout << "Elementary cell: " << this->eff[0]/fl_units[0] << " a0 x " << this->eff[1]/fl_units[1] << " a0 x " << this->eff[2]/fl_units[2] << " a0.\n\n"; - for(int i=0; i < 3; i++) - this->fl_unit[i] = this->eff[i] / (double)fl_units[i]; - this->pfill = N_boxes / ((double)fl_units[0]*fl_units[1]*fl_units[2]); - - /* - * additional free volume (outside the pore) - */ - if((wo_wall*box[2] > 2.0*shielding) && (wo_wall != 1.0)) - { - this->do_fill_ext = true; - - this->ext[0] = this->eff[0]; - this->off_ext[0] = this->off[0]; - this->ext[1] = this->box[1] - this->eff[1]; // Wand + Shielding - this->off_ext[1] = this->off[1] + this->eff[1]; - if(off_ext[1] < 0.0) off_ext[1] += box[1]; - else if(off_ext[1] > box[1]) off_ext[1] -= box[1]; - this->ext[2] = this->box[2]*wo_wall - 2.0*shielding; - this->off_ext[2] = this->off[2] + 0.5*(1.0 - wo_wall)*box[2] + shielding; // Shielding (Aussenseite der Wand) - double V_ext = this->ext[0] * this->ext[1] * this->ext[2]; - std::cout << "Additionally available: " << V_ext << " * 1.4818e-04 nm^3," - << "i.e. " << this->ext[0] << " (off " << this->off_ext[0] - << ") x " << this->ext[1] << " (off " << this->off_ext[1] - << ") x " << this->ext[2] << " (off " << this->off_ext[2] << ").\n"; - - double N_id_ext = rho*V_ext; - double N_boxes_ext = N_id_ext / 3.0; - fl_units_ext[1] = round( - pow( - (N_boxes_ext * this->ext[1] * this->ext[1]) - / (this->ext[0] * this->ext[2]), 1.0/3.0 - ) - ); - if(fl_units_ext[1] == 0) fl_units_ext[1] = 1; - double bxbz_id_ext = N_boxes_ext / fl_units_ext[1]; - fl_units_ext[0] = round(sqrt(this->ext[0] * bxbz_id_ext / this->ext[2])); - if(fl_units_ext[0] == 0) this->fl_units_ext[0] = 1; - fl_units_ext[2] = ceil(bxbz_id_ext / fl_units_ext[0]); - std::cout << "Elementary cell (extension): " << this->ext[0]/fl_units_ext[0] - << " a0 x " << this->ext[1]/fl_units_ext[1] << " a0 x " - << this->ext[2]/fl_units_ext[2] << " a0.\n\n"; - for(int i=0; i < 3; i++) - this->fl_unit_ext[i] = this->ext[i] / (double)fl_units_ext[i]; - this->pfill_ext = N_boxes_ext / ((double)fl_units_ext[0]*fl_units_ext[1]*fl_units_ext[2]); - } - else - { - this->do_fill_ext = false; - this->pfill_ext = 0.0; - } -} - -void Domain::specifyNanotube(double rho, double m_per_n, unsigned N) -{ - std::cout << "Nanotubes are not yet implemented.\n"; - exit(16); -} - diff --git a/tools/standalone-generators/mkcp/Domain.h b/tools/standalone-generators/mkcp/Domain.h deleted file mode 100644 index 2974e358b2..0000000000 --- a/tools/standalone-generators/mkcp/Domain.h +++ /dev/null @@ -1,179 +0,0 @@ -#include -#include -#include -#include - - -#define TIME 20111116 - -#define FLUID_NIL -1 -#define FLUID_AR 0 -#define FLUID_CH4 1 -#define FLUID_C2H6 2 -#define FLUID_N2 3 -#define FLUID_CO2 4 -#define FLUID_AVE 5 -#define FLUID_H2O 6 -#define FLUID_CH3OH 7 -#define FLUID_C6H14 8 - -#define EPS_AR 4.36704e-04 -#define SIG_AR 6.40920 -#define ARMASS 0.039948 -#define EPS_CH4 5.54383e-04 -#define SIG_CH4 7.03753 -#define CH4MASS 0.016042 -#define EPS_C2H6 4.33822e-04 -#define SIG_C2H6 6.59439 -#define C2H6MASS 0.03007 -#define C2H6LONG 4.49037 -#define QDR_C2H6 -0.61537 -#define EPS_N2 1.10512e-04 -#define SIG_N2 6.27597 -#define N2MASS 0.0280134 -#define N2LONG 1.97741 -#define QDR_N2 -1.07038 -#define EPS_CO2 4.21883e-04 -#define SIG_CO2 5.64027 -#define CO2MASS 0.0440095 -#define CO2LONG 4.56860 -#define QDR_CO2 -2.82059 -#define SIG_AVE 5.89406 -#define EPS_AVE 0.00189803 -#define AVEMASS 0.018015 - -#define EPS_OH2O 6.58951e-04 -#define SIG_OH2O 5.89277 -#define OH2OMASS 0.0159994 -#define HH2OMASS 0.0010079 -#define CHG_HH2O +0.419548 -#define CHG_EH2O -0.839096 -#define R0_O_H2O +0.0 -#define R1_O_H2O -0.14948 -#define R2_O_H2O +0.0 -#define R0_H1H2O -1.72579 -#define R1_H1H2O +1.18644 -#define R2_H1H2O +0.0 -#define R0_H2H2O +1.72579 -#define R1_H2H2O +1.18644 -#define R2_H2H2O +0.0 -#define R0_E_H2O +0.0 -#define R1_E_H2O +0.23765 -#define R2_E_H2O +0.0 -#define I_XX_H2O 0.0031950 -#define I_YY_H2O 0.0060038 -#define I_ZZ_H2O 0.0091988 -#define H2O_LONG 0.29896 - -#define EPS_CCH3OH 3.81893e-04 -#define SIG_CCH3OH 7.09460 -#define CCH3OHMASS 0.0150348 -#define CHG_CCH3OH +0.24746 -#define EPS_OCH3OH 2.78297e-04 -#define SIG_OCH3OH 5.72587 -#define OCH3OHMASS 0.0159994 -#define CHG_OCH3OH -0.67874 -#define HCH3OHMASS 0.0010079 -#define CHG_HCH3OH +0.43128 -#define R0_C_CH3OH -1.44677 -#define R1_C_CH3OH -0.05327 -#define R2_C_CH3OH +0.0 -#define R0_O_CH3OH +1.24534 -#define R1_O_CH3OH -0.05327 -#define R2_O_CH3OH +0.0 -#define R0_H_CH3OH +1.81292 -#define R1_H_CH3OH +1.64012 -#define R2_H_CH3OH +0.0 -#define I_XX_CH3OH 0.0027993 -#define I_YY_CH3OH 0.0595954 -#define I_ZZ_CH3OH 0.0623947 -#define CH3OH_LONG 2.6942 - -#define EPS_FC6H14 3.10e-04 -#define SIG_FC6H14 7.086 -#define FC6H14MASS 0.01503 -#define EPS_MC6H14 1.46e-04 -#define SIG_MC6H14 7.464 -#define MC6H14MASS 0.01403 -#define R0_F1C6H14 -6.139 -#define R1_F1C6H14 +0.418 -#define R2_F1C6H14 +0.0 -#define R0_M1C6H14 -3.606 -#define R1_M1C6H14 -1.015 -#define R2_M1C6H14 +0.0 -#define R0_M2C6H14 -1.266 -#define R1_M2C6H14 +0.716 -#define R2_M2C6H14 +0.0 -#define R0_M3C6H14 +1.266 -#define R1_M3C6H14 -0.716 -#define R2_M3C6H14 +0.0 -#define R0_M4C6H14 +3.606 -#define R1_M4C6H14 +1.015 -#define R2_M4C6H14 +0.0 -#define R0_F2C6H14 +6.139 -#define R1_F2C6H14 -0.418 -#define R2_F2C6H14 +0.0 -#define I_XX_C6H14 0.04855 -#define I_YY_C6H14 1.54288 -#define I_ZZ_C6H14 1.59144 -#define C6H14_LONG 12.31 - -#define FORMAT_BUCHHOLZ 0 -#define FORMAT_BRANCH 1 -#define FORMAT_BERNREUTHER 2 - -#define FLOW_NONE 0 -#define FLOW_COUETTE 1 -#define FLOW_POISEUILLE 2 - -class Domain -{ - public: - Domain( - int sp_flow, double sp_bondlength, double sp_rho, int sp_d, - int sp_fluid, int sp_fluid2, double sp_h, double sp_ETA, double sp_ETA2, double sp_ETAF, double sp_SIG_REF, - double sp_EPS_REF, double sp_REFMASS, bool sp_muVT, - bool sp_nanotube, double sp_m_per_n, unsigned sp_N, double sp_T, - double sp_XI, double sp_XI2, double sp_XIF, double sp_wo_wall - ); - void write( - char* prefix, double a, bool empty, int format, double mu, - double TAU, double U, bool original, double wo_acceleration, - double polarity, bool WLJ, bool symmetric, bool widom, double x - ); - - private: - void specifyGraphite(double rho, unsigned N); - void specifyNanotube(double rho, double m_per_n, unsigned N); - void writeGraphite( - char* prefix, double a, bool empty, int format, double mu, - double TAU, double U, bool original, double wo_acceleration, - double polarity, bool WLJ, bool symmetric, bool widom, double x - ); - void writeNanotube( - char* prefix, double a, bool empty, int format, double mu, - double TAU, double U, bool original, double wo_acceleration, - double polarity, bool WLJ, bool symmetric, bool widom, double x - ); - - bool muVT, nanotube; - int flow, fluid, fluid2; - unsigned d, m, n; - double bondlength, h, ETA, ETA2, ETAF, SIG_REF, EPS_REF, REFMASS, T, XI, XI2, XIF, wo_wall; - - double box[3], // box dimensions (Couette == two boxes) - eff[3], // effective space for the fluid - off[3]; // offset coordinates for the fluid - double shielding; // shielding of the wall due to ETA*SIGMA - - double fl_unit[3]; - unsigned fl_units[3]; - double pfill; // probability of placing a molecule in a fluid unit - - bool do_fill_ext; - double ext[3], off_ext[3]; - double fl_unit_ext[3]; - unsigned fl_units_ext[3]; - double pfill_ext; // probability of placing a molecule in a fluid unit -}; - diff --git a/tools/standalone-generators/mkcp/Graphit.cpp b/tools/standalone-generators/mkcp/Graphit.cpp deleted file mode 100644 index 387e394143..0000000000 --- a/tools/standalone-generators/mkcp/Graphit.cpp +++ /dev/null @@ -1,237 +0,0 @@ -#include "Graphit.h" -#include "Random.h" - -#include -#include - -int Graphit::getNumberOfAtoms() -{ - return numberOfAtoms; -} - -void Graphit::calculateCoordinatesOfAtoms( - int numberOfLayers, double xLength, double zLength, double A, double wo_wall -) { - double B = 2.0 * A; - double C = 0.86602540378 * A; // sin(pi/3) - - double xCoor = 0.05*A; - double yCoor = 0.0; - double zCoor = 0.05*C; - numberOfAtoms=0; - for(int Index=1; Index<=numberOfLayers; Index++) - { - /******************* - * Creates the odd layer (1,3,5,..) - *******************/ - if(Index%2==1) - { - xCoor = 0.1*A; - zCoor = 0.1*C; - int i = 0; - int j = 0; - int temp = 0; - if((xCoor < xLength && zCoor < zLength) && ((zCoor < (0.5 - 0.5*wo_wall)*zLength) || (zCoor > (0.5 + 0.5*wo_wall)*zLength))) - { - coordinatesOfAtoms[0][numberOfAtoms] = xCoor; - coordinatesOfAtoms[1][numberOfAtoms] = yCoor; - coordinatesOfAtoms[2][numberOfAtoms] = zCoor; - componentsOfAtoms[numberOfAtoms] = this->comp(i, j); - numberOfAtoms++; - } - while(zCoor < zLength) - { - while(xCoor < xLength) - { - if (i%2==0) - { - xCoor += B; - i=1; - } - else - { - xCoor += A; - i=0; - temp=0; - } - if((xCoor < xLength && zCoor < zLength) && ((zCoor < (0.5 - 0.5*wo_wall)*zLength) || (zCoor > (0.5 + 0.5*wo_wall)*zLength))) - { - coordinatesOfAtoms[0][numberOfAtoms] = xCoor; - coordinatesOfAtoms[1][numberOfAtoms] = yCoor; - coordinatesOfAtoms[2][numberOfAtoms] = zCoor; - componentsOfAtoms[numberOfAtoms] = this->comp(i, j); - numberOfAtoms++; - } - } - zCoor += C; - j++; - if(j == 14) j = 0; - if(j%2==1) - { - i = 1; - xCoor = 0.6*A; - } - else - { - i = 0; - xCoor = 0.1*A; - } - if((xCoor < xLength && zCoor < zLength) && ((zCoor < (0.5 - 0.5*wo_wall)*zLength) || (zCoor > (0.5 + 0.5*wo_wall)*zLength))) - { - coordinatesOfAtoms[0][numberOfAtoms] = xCoor; - coordinatesOfAtoms[1][numberOfAtoms] = yCoor; - coordinatesOfAtoms[2][numberOfAtoms] = zCoor; - componentsOfAtoms[numberOfAtoms] = this->comp(i, j); - numberOfAtoms++; - } - temp = j; - } - } - /******************** - * Creates the even layer - *********************/ - else - { - // Odd Layer has a different starting point - xCoor = B-0.4*A; - zCoor = 0.1*C; - int i = 0; //Makes it possible to change the different lengths in one direction - int j = 7; //Necessary for changing the different startpoints of the xCoor - - // Code for writing the Coordinates together - if((xCoor < xLength && zCoor < zLength) && ((zCoor < (0.5 - 0.5*wo_wall)*zLength) || (zCoor > (0.5 + 0.5*wo_wall)*zLength))) - { - coordinatesOfAtoms[0][numberOfAtoms] = xCoor; - coordinatesOfAtoms[1][numberOfAtoms] = yCoor; - coordinatesOfAtoms[2][numberOfAtoms] = zCoor; - componentsOfAtoms[numberOfAtoms] = this->comp(i, j); - numberOfAtoms++; - } - while(zCoor < zLength) - { - while(xCoor < xLength) - { - if (i%2==0) - { - xCoor += A; - i=0; - } - else - { - xCoor += B; - i=1; - } - if((xCoor < xLength && zCoor < zLength) && ((zCoor < (0.5 - 0.5*wo_wall)*zLength) || (zCoor > (0.5 + 0.5*wo_wall)*zLength))) - { - coordinatesOfAtoms[0][numberOfAtoms] = xCoor; - coordinatesOfAtoms[1][numberOfAtoms] = yCoor; - coordinatesOfAtoms[2][numberOfAtoms] = zCoor; - componentsOfAtoms[numberOfAtoms] = this->comp(i, j); - numberOfAtoms++; - } - - } - zCoor += C; - j++; - if(j == 14) j = 0; - if(j%2==1) - { - i = 1; - xCoor = 0.1*A; - } - else - { - i = 0; - xCoor = B-0.4*A; - } - if((xCoor < xLength && zCoor < zLength) && ((zCoor < (0.5 - 0.5*wo_wall)*zLength) || (zCoor > (0.5 + 0.5*wo_wall)*zLength))) - { - coordinatesOfAtoms[0][numberOfAtoms] = xCoor; - coordinatesOfAtoms[1][numberOfAtoms] = yCoor; - coordinatesOfAtoms[2][numberOfAtoms] = zCoor; - componentsOfAtoms[numberOfAtoms] = this->comp(i, j); - numberOfAtoms++; - } - } - } - yCoor = yCoor + Z; - } -} - -double Graphit::getX(int number) -{ - return coordinatesOfAtoms[0][number]; -} - -double Graphit::getY(int number) -{ - return coordinatesOfAtoms[1][number]; -} - -double Graphit::getZ(int number) -{ - return coordinatesOfAtoms[2][number]; -} - -void Graphit::calculateVelocities(double T, double U) -{ - double absoluteVelocity = sqrt(3.0*T / ATOMIC_MASS_C); - double phi, omega; - - Random* r = new Random(); - r->init((int)(10000000.0*T) - (int)(3162300.0*U)); - - for(int i=0; i < this->numberOfAtoms; i++) - { - phi = 6.283185 * r->rnd(); - omega = 6.283185 * r->rnd(); - velocitiesOfAtoms[0][i] = absoluteVelocity*cos(phi)*cos(omega); - velocitiesOfAtoms[1][i] = absoluteVelocity*cos(phi)*sin(omega); - velocitiesOfAtoms[2][i] = U + absoluteVelocity*sin(phi); - } -} - -double Graphit::getVelocityX(int number) -{ - return velocitiesOfAtoms[0][number]; -} - -double Graphit::getVelocityY(int number) -{ - return velocitiesOfAtoms[1][number]; -} - -double Graphit::getVelocityZ(int number) -{ - return velocitiesOfAtoms[2][number]; -} - -void Graphit::reset() -{ - this->numberOfAtoms = 0; - for(int d=0; d<4; d++) - { - this->coordinatesOfAtoms[d].clear(); - this->velocitiesOfAtoms[d].clear(); - } -} - -unsigned Graphit::comp(int ti, int tj) -{ - int i = ti % 2; - int j = tj % 14; - - // std::cout << "\t\t(i, j) \t=\t (" << i << ", " << j << ")\n"; - // std::cout.flush(); - - if((j == 2) && (i == 1)) return CID_I; - else if((j == 3) && (i == 0)) return CID_I; - else if((j == 3) && (i == 1)) return CID_ZI; - else if((j == 7) && (i == 0)) return CID_I; - else if((j == 7) && (i == 1)) return CID_Z; - else if((j == 11) && (i == 0)) return CID_I; - else if((j == 11) && (i == 1)) return CID_ZI; - else if((j == 12) && (i == 1)) return CID_I; - else return CID_C; -} - diff --git a/tools/standalone-generators/mkcp/Graphit.h b/tools/standalone-generators/mkcp/Graphit.h deleted file mode 100644 index 504651045c..0000000000 --- a/tools/standalone-generators/mkcp/Graphit.h +++ /dev/null @@ -1,60 +0,0 @@ -#include - -#define SIG_WANG 6.2860 -#define EPS_WANG 6.909e-05 -#define ATOMIC_MASS_C 0.0120107 -#define DEFAULT_A 2.6853 -#define DEFAULT_A_ORIGINAL 2.7666 -#define Z 6.33 - -#define NCOMP_POLAR 4 -#define CID_I 2 -#define CID_C 1 -#define CID_Z 4 -#define CID_ZI 3 - -#define TERSOFF_A 51.214 -#define TERSOFF_B 12.740 -#define TERSOFF_LAMBDA 1.8982 -#define TERSOFF_LAMBDA_ORIG 1.8457 -#define TERSOFF_MU 1.2039 -#define TERSOFF_MU_ORIG 1.1705 -#define TERSOFF_R 3.78 -#define TERSOFF_R_ORIG 3.40 -#define TERSOFF_S 4.44 -#define TERSOFF_S_ORIG 3.97 -#define TERSOFF_C 38049 -#define TERSOFF_D 4.384 -#define TERSOFF_H -0.57058 -#define TERSOFF_N 0.72751 -#define TERSOFF_BETA 1.5724e-07 - -class Graphit -{ - public: - int getNumberOfAtoms(); - - void calculateCoordinatesOfAtoms( - int numberOfLayers, double xLength, double zLength, double A, double wo_wall - ); - double getX(int number); - double getY(int number); - double getZ(int number); - - void calculateVelocities(double T, double U); - double getVelocityX(int number); - double getVelocityY(int number); - double getVelocityZ(int number); - - unsigned getComponent(int number) { return this->componentsOfAtoms[number]; }; - - void reset(); - - private: - unsigned comp(int ti, int tj); - - int numberOfAtoms; - std::map coordinatesOfAtoms[3]; - std::map velocitiesOfAtoms[3]; - std::map componentsOfAtoms; -}; diff --git a/tools/standalone-generators/mkcp/Random.cpp b/tools/standalone-generators/mkcp/Random.cpp deleted file mode 100644 index 2a96b04a08..0000000000 --- a/tools/standalone-generators/mkcp/Random.cpp +++ /dev/null @@ -1,30 +0,0 @@ -#include "Random.h" - -Random::Random() { this->init(8624); } - -void Random::init(int seed) -{ - this->ix_muVT = (seed ^ (int)888889999) | (int)1; - this->iy_muVT = seed ^ (int)777755555; - // Calculate normalization factor - this->am_muVT = 2.0 / (1.0 + (unsigned)((int)-1)); -} - -float Random::rnd() -{ - float rnd; - const int IA = 16807; - const int IM = 2147483647; - const int IQ = 127773; - const int IR = 2836; - - ix_muVT ^= (ix_muVT >> 13); - ix_muVT ^= (ix_muVT << 17); - ix_muVT ^= (ix_muVT >> 5); - - int k = iy_muVT / IQ; - iy_muVT = IA * (iy_muVT - k*IQ) - IR*k; - if(iy_muVT < 0) iy_muVT += IM; - rnd = am_muVT * ((IM & (ix_muVT ^ iy_muVT)) | (int)1); - return rnd; -} diff --git a/tools/standalone-generators/mkcp/Random.h b/tools/standalone-generators/mkcp/Random.h deleted file mode 100644 index 073f88e2e3..0000000000 --- a/tools/standalone-generators/mkcp/Random.h +++ /dev/null @@ -1,14 +0,0 @@ -class Random -{ - public: - Random(); - void init(int seed); - - float rnd(); - - int getIX() { return this->ix_muVT; } - - private: - int ix_muVT, iy_muVT; - float am_muVT; -}; diff --git a/tools/standalone-generators/mkcp/main.cpp b/tools/standalone-generators/mkcp/main.cpp deleted file mode 100644 index 0fd3132ac5..0000000000 --- a/tools/standalone-generators/mkcp/main.cpp +++ /dev/null @@ -1,594 +0,0 @@ -#include "Domain.h" -#include "Graphit.h" - -#include -#include -#include -#include - - -#define DEFAULT_TAU 1.0e+10 - -#define DEFAULT_CC_BONDLENGTH 2.6853 -#define ORIGINAL_CC_BONDLENGTH 2.7609 - -// interlayer spacing graphite: 3.35 A -// interlayer spacing MWCNT: 3.40 A -// (http://www.wag.caltech.edu/foresight/foresight_2.html) - -int main(int argc, char** argv) -{ - const char* usage = "usage: mkcp (-C|-P|-0) [-a ] [-A ] -c -d [-e] [-E] [-f ] [-g ] -h [-H ] [-I ] [-J ] [-k] [-l] [-L] [-m ] [-M ] -N [-p ] [-r] [-s ] [-S] [-t ] -T [-u] -U [-v ] [-V ] [-x <2nd comp. mole fract.>] [-Y ] [-3 ] [-4 ] [-5 ] [-8]\n\n-A\treduced C-C bond length; default: 2.6853 a0 = 0.1421 nm, original Tersoff: 2.7609 a0\n-C\tCouette flow (flag followed by output prefix)\n-e\tuse B-e-rnreuther format\n-E\tgenerate an empty channel\n-f\tAr (default), Ave, CH4, C2H6, N2, CO2, H2O, CH3OH, or C6H14\n-H, -I\tdefault: eta according to Wang et al.\n-J\tdefault: eta = 1\n-k\tonly harmonic potentials active within the wall (default for polar walls)\n-l\tLJ interaction within the wall as well (default for unpolar walls)\n-L\tuse Lennard-Jones units instead of atomic units (cf. atomic_units.txt)\n-M\tgenerate a carbon nanotube (only for Poiseuille flow)\n-p\tdefault polarity coefficient: 0 (unpolar walls)\n-r\tuse b-r-anch format (active by default)\n-s\tgiven in units of Angstrom; default: 1 = 0.5291772 A 0\n-S\tsymmetric system (with two identical fluid components)\n-t\tdefault: tau extremely large (about 30 us)\n-u\tuse B-u-chholz format\n-w\tWidom\n-W\tgiven in units of K; default value: 1 = 315774.5 K\n-x\tdefault value: 0 (i.e. pure first comp. fluid)\n-Y\tgiven in units of g/mol; default value: 1 = 1000 g/mol\n-0\tstatic scenario without flow (followed by output prefix)\n-3, -4\tdefault: xi according to Wang et al.\n-5\tdefault: xi = 1\n-8\toriginal Tersoff potential as published in the 80s\n"; - if((argc < 15) || (argc > 69)) - { - std::cout << "There are " << argc - << " arguments where 15 to 61 should be given.\n\n"; - std::cout << usage; - return 1; - } - - unsigned d, N; - double h, T, rho, U, XI, ETA, XI2, ETA2, XIF, ETAF, TAU, bondlength, a, m_per_n; - double mu = 0.0; - double x = 0.0; - int fluid, flow; - int fluid2 = FLUID_NIL; - int format = FORMAT_BRANCH; - char* prefix = (char*)0; - - double SIG_REF = 1.0; - double EPS_REF = 1.0; - double REFMASS = 1.0; - - double wo_acceleration = 0.0; - double wo_wall = 0.0; - double polarity = 0.0; - - bool in_a = false; - bool in_bondlength = false; - bool empty = false; - bool in_ETA = false; - bool in_ETA2 = false; - bool in_ETAF = false; - bool LJunits = false; - bool in_TAU = false; - bool in_N = false; - bool in_XI = false; - bool in_XI2 = false; - bool in_XIF = false; - bool in_rho = false; - bool in_d = false; - bool in_fluid = false; - bool in_fluid2 = false; - bool in_h = false; - bool in_T = false; - bool in_U = false; - bool in_WLJ = false; - bool in_x = false; - bool original = false; - bool muVT = false; - bool nanotube = false; - bool WLJ = true; - bool symmetric = false; - bool widom = false; - - for(int i=1; i < argc; i++) - { - if(*argv[i] != '-') - { - std::cout << "Flag expected where '" << argv[i] - << "' was given.\n\n"; - std::cout << usage; - return 2; - } - for(int j=1; argv[i][j]; j++) - { - // std::cout << argv[i][j] << "\n"; - if(argv[i][j] == 'a') - { - in_a = true; - i++; - a = atof(argv[i]); - break; - } - else if(argv[i][j] == 'A') - { - in_bondlength = true; - i++; - bondlength = atof(argv[i]); - break; - } - else if(argv[i][j] == 'c') - { - in_rho = true; - i++; - rho = atof(argv[i]); - break; - } - else if(argv[i][j] == 'C') - { - flow = FLOW_COUETTE; - i++; - prefix = argv[i]; - break; - } - else if(argv[i][j] == 'd') - { - in_d = true; - i++; - d = atoi(argv[i]); - break; - } - else if(argv[i][j] == 'e') format = FORMAT_BERNREUTHER; - else if(argv[i][j] == 'E') empty = true; - else if(argv[i][j] == 'f') - { - in_fluid = true; - i++; - if(!strcmp(argv[i], "Ar")) fluid = FLUID_AR; - else if(!strcmp(argv[i], "Ave")) fluid = FLUID_AVE; - else if(!strcmp(argv[i], "CH4")) fluid = FLUID_CH4; - else if(!strcmp(argv[i], "C2H6")) fluid = FLUID_C2H6; - else if(!strcmp(argv[i], "N2")) fluid = FLUID_N2; - else if(!strcmp(argv[i], "CO2")) fluid = FLUID_CO2; - else if(!strcmp(argv[i], "H2O")) fluid = FLUID_H2O; - else if(!strcmp(argv[i], "CH3OH")) fluid = FLUID_CH3OH; - else if(!strcmp(argv[i], "C6H14")) fluid = FLUID_C6H14; - else - { - std::cout << "Fluid '" << argv[i] - << "' is not available.\n\n" << usage; - return 9; - } - break; - } - else if(argv[i][j] == 'g') - { - in_fluid2 = true; - i++; - if(!strcmp(argv[i], "Ar")) fluid2 = FLUID_AR; - else if(!strcmp(argv[i], "Ave")) fluid2 = FLUID_AVE; - else if(!strcmp(argv[i], "CH4")) fluid2 = FLUID_CH4; - else if(!strcmp(argv[i], "C2H6")) fluid2 = FLUID_C2H6; - else if(!strcmp(argv[i], "N2")) fluid2 = FLUID_N2; - else if(!strcmp(argv[i], "CO2")) fluid2 = FLUID_CO2; - else if(!strcmp(argv[i], "H2O")) fluid2 = FLUID_H2O; - else if(!strcmp(argv[i], "CH3OH")) fluid2 = FLUID_CH3OH; - else if(!strcmp(argv[i], "C6H14")) fluid2 = FLUID_C6H14; - else - { - std::cout << "(Secondary) fluid '" << argv[i] - << "' is not available.\n\n" << usage; - return 99; - } - break; - } - else if(argv[i][j] == 'h') - { - in_h = true; - i++; - h = atof(argv[i]); - break; - } - else if(argv[i][j] == 'H') - { - in_ETA = true; - i++; - ETA = atof(argv[i]); - break; - } - else if(argv[i][j] == 'I') - { - in_ETA2 = true; - i++; - ETA2 = atof(argv[i]); - break; - } - else if(argv[i][j] == 'J') - { - in_ETAF = true; - i++; - ETAF = atof(argv[i]); - break; - } - else if(argv[i][j] == 'k') - { - in_WLJ = true; - WLJ = false; - } - else if(argv[i][j] == 'l') - { - in_WLJ = true; - WLJ = true; - } - else if(argv[i][j] == 'L') LJunits = true; - else if(argv[i][j] == 'm') - { - muVT = true; - i++; - mu = atof(argv[i]); - break; - } - else if(argv[i][j] == 'M') - { - nanotube = true; - i++; - m_per_n = atof(argv[i]); - break; - } - else if(argv[i][j] == 'N') - { - in_N = true; - i++; - N = atoi(argv[i]); - break; - } - else if(argv[i][j] == 'p') - { - i++; - polarity = atof(argv[i]); - if(!in_WLJ) - { - WLJ = (polarity == 0.0); - } - break; - } - else if(argv[i][j] == 'P') - { - flow = FLOW_POISEUILLE; - i++; - prefix = argv[i]; - break; - } - else if(argv[i][j] == 'r') format = FORMAT_BRANCH; - else if(argv[i][j] == 'S') symmetric = true; - else if(argv[i][j] == 's') - { - i++; - SIG_REF = atof(argv[i]) / 0.5291772; - break; - } - else if(argv[i][j] == 't') - { - in_TAU = true; - i++; - TAU = atof(argv[i]); - break; - } - else if(argv[i][j] == 'T') - { - in_T = true; - i++; - T = atof(argv[i]); - break; - } - else if(argv[i][j] == 'u') format = FORMAT_BUCHHOLZ; - else if(argv[i][j] == 'U') - { - in_U = true; - i++; - U = atof(argv[i]); - break; - } - else if(argv[i][j] == 'v') - { - i++; - wo_acceleration = atof(argv[i]); - break; - } - else if(argv[i][j] == 'V') - { - i++; - wo_wall = atof(argv[i]); - break; - } - else if(argv[i][j] == 'w') widom = true; - else if(argv[i][j] == 'W') - { - i++; - EPS_REF = atof(argv[i]) / 315774.5; - break; - } - else if(argv[i][j] == 'x') - { - in_x = true; - i++; - x = atof(argv[i]); - break; - } - else if(argv[i][j] == 'Y') - { - i++; - REFMASS = atof(argv[i]) / 1000.0; - break; - } - else if(argv[i][j] == '0') - { - flow = FLOW_NONE; - i++; - prefix = argv[i]; - break; - } - else if(argv[i][j] == '3') - { - in_XI = true; - i++; - XI = atof(argv[i]); - break; - } - else if(argv[i][j] == '4') - { - in_XI2 = true; - i++; - XI2 = atof(argv[i]); - break; - } - else if(argv[i][j] == '5') - { - in_XIF = true; - i++; - XIF = atof(argv[i]); - break; - } - else if(argv[i][j] == '8') original = true; - else - { - std::cout << "Invalid flag '-" << argv[i][j] - << "' was detected.\n\n" << usage; - return 2; - } - } - } - - if(wo_wall < 0.003) wo_wall = 0.0; - else if(wo_wall > 0.997) wo_wall = 1.0; - - if(nanotube) - { - if(flow == FLOW_COUETTE) - { - std::cout << "Couette flow is not well-defined for nanotubes.\n\n" - << usage; - return 18; - } - - std::cout << "Carbon nanotubes are not yet implemented.\n\n" - << usage; - return 11; - } - if(WLJ && (polarity != 0.0)) - { - std::cout << "Severe error: Non-zero polarity is incompatible with the LJ interlayer interaction.\n\n" - << usage; - return 181; - } - if(!prefix) - { - std::cout << "You have to specify an output prefix " - << "via the -C, -P, or -0 option.\n\n" << usage; - return 5; - } - if(!in_rho) - { - std::cout << "Fatal error: no fluid density was specified.\n\n"; - std::cout << usage; - return 6; - } - if(!in_d) - { - std::cout << "Wall thickness unknown. Reconsider parameter set.\n\n"; - std::cout << usage; - return 7; - } - if(format == FORMAT_BERNREUTHER) - { - std::cout << "B-e-rnreuther format (flag -e) " - << "is unavailable at present.\n\n" << usage; - return 14; - } - - if(fluid == fluid2) - { - symmetric = true; - fluid2 = FLUID_NIL; - } - if((x < 0.0) || (x > 1.0)) - { - std::cout << "Invalid mole fraction x = " << x << ".\n\n" << usage; - return 15; - } - if((in_fluid2 == false) || (fluid2 == FLUID_NIL)) - { - x = 0.0; - } - if(x == 1.0) - { - fluid = fluid2; - x = 0.0; - } - if(x == 0.0) - { - in_fluid2 = false; - fluid2 = FLUID_NIL; - } - - if(symmetric && (fluid2 != FLUID_NIL)) - { - std::cout << "The symmetric scenario cannot contain an (actual) fluid mixture.\n\n" << usage; - return 16; - } - if(symmetric) x = 0.5; - - if(!in_fluid) fluid = FLUID_AR; - double SIG_FLUID, EPS_FLUID, FLUIDMASS; - double SIG_FLUID2, EPS_FLUID2, FLUIDMASS2; - if(fluid == FLUID_CH4) - { - SIG_FLUID = SIG_CH4; - EPS_FLUID = EPS_CH4; - FLUIDMASS = CH4MASS; - } - else if(fluid == FLUID_C2H6) - { - SIG_FLUID = SIG_C2H6; - EPS_FLUID = EPS_C2H6; - FLUIDMASS = C2H6MASS; - } - else if(fluid == FLUID_N2) - { - SIG_FLUID = SIG_N2; - EPS_FLUID = EPS_N2; - FLUIDMASS = N2MASS; - } - else if(fluid == FLUID_CO2) - { - SIG_FLUID = SIG_CO2; - EPS_FLUID = EPS_CO2; - FLUIDMASS = CO2MASS; - } - else if(fluid == FLUID_AVE) // Avendaño mode - { - SIG_FLUID = SIG_AVE; - EPS_FLUID = EPS_AVE; - FLUIDMASS = AVEMASS; - } - else if(fluid == FLUID_AR) - { - SIG_FLUID = SIG_AR; - EPS_FLUID = EPS_AR; - FLUIDMASS = ARMASS; - } - else // effectively sets ETA = XI = 1 in all other cases - { - SIG_FLUID = SIG_WANG; - EPS_FLUID = EPS_WANG; - if(fluid == FLUID_H2O) FLUIDMASS = OH2OMASS + 2.0*HH2OMASS; - else if(fluid == FLUID_CH3OH) FLUIDMASS = CCH3OHMASS + OCH3OHMASS + HCH3OHMASS; - else FLUIDMASS = 2.0*FC6H14MASS + 4.0*MC6H14MASS; - } - if(fluid2 == FLUID_CH4) - { - SIG_FLUID2 = SIG_CH4; - EPS_FLUID2 = EPS_CH4; - FLUIDMASS2 = CH4MASS; - } - else if(fluid2 == FLUID_C2H6) - { - SIG_FLUID2 = SIG_C2H6; - EPS_FLUID2 = EPS_C2H6; - FLUIDMASS2 = C2H6MASS; - } - else if(fluid2 == FLUID_N2) - { - SIG_FLUID2 = SIG_N2; - EPS_FLUID2 = EPS_N2; - FLUIDMASS2 = N2MASS; - } - else if(fluid2 == FLUID_CO2) - { - SIG_FLUID2 = SIG_CO2; - EPS_FLUID2 = EPS_CO2; - FLUIDMASS2 = CO2MASS; - } - else if(fluid2 == FLUID_AVE) // Avendaño mode - { - SIG_FLUID2 = SIG_AVE; - EPS_FLUID2 = EPS_AVE; - FLUIDMASS2 = AVEMASS; - } - else if(fluid2 == FLUID_AR) - { - SIG_FLUID2 = SIG_AR; - EPS_FLUID2 = EPS_AR; - FLUIDMASS2 = ARMASS; - } - else // effectively sets ETA2 = XI2 = 1 in all other cases - { - SIG_FLUID2 = SIG_WANG; - EPS_FLUID2 = EPS_WANG; - if(fluid2 == FLUID_H2O) FLUIDMASS2 = OH2OMASS + 2.0*HH2OMASS; - else if(fluid2 == FLUID_CH3OH) FLUIDMASS2 = CCH3OHMASS + OCH3OHMASS + HCH3OHMASS; - else FLUIDMASS2 = 2.0*FC6H14MASS + 4.0*MC6H14MASS; - } - if(!in_N) - { - std::cout << "Missing essential input parameter " - << "N (number of fluid molecules).\n\n" << usage; - return 12; - } - if(!in_T) - { - std::cout << "No temperature specified: aborting.\n\n"; - std::cout << usage; - return 13; - } - if(!in_U) U = 0.0; - else if(FLOW_NONE && (U != 0.0)) - { - std::cout << "Specified flow velocity U = " << U - << " is incompatible with -0 option (static scenario).\n\n" - << usage; - return 15; - } - if(muVT && symmetric) - { - std::cout << "The combination between the uVT (-m) and symmetry (-S) options remains unsupported.\n\n" - << usage; - return 115; - } - - if(!in_h) h = pow((double)N/rho, 1.0/3.0); - if(!in_ETA) ETA = 0.5*(1.0 + SIG_WANG/SIG_FLUID); - if(!in_XI) XI = sqrt(EPS_WANG / EPS_FLUID); - if(!in_ETA2) ETA2 = (SIG_WANG + SIG_FLUID2) / (SIG_FLUID + SIG_FLUID2); - if(!in_XI2) XI2 = sqrt(EPS_WANG / EPS_FLUID); - if(!in_ETAF) ETAF = 1.0; - if(!in_XIF) XIF = 1.0; - - SIG_REF = LJunits? SIG_FLUID: 1.0; - EPS_REF = LJunits? EPS_FLUID: 1.0; - REFMASS = LJunits? FLUIDMASS: 1.0; - if(LJunits) - { - SIG_REF = SIG_FLUID; - EPS_REF = EPS_FLUID; - REFMASS = FLUIDMASS; - } - double REFTIME = SIG_REF * sqrt(REFMASS / EPS_REF); - - h *= SIG_REF; - T *= EPS_REF; - rho /= (SIG_REF * SIG_REF * SIG_REF); - U *= SIG_REF / REFTIME; - TAU *= REFTIME; - bondlength *= SIG_REF; - a *= SIG_REF / (REFTIME * REFTIME); - mu *= EPS_REF; - - if(!in_bondlength) - bondlength = original? ORIGINAL_CC_BONDLENGTH - : DEFAULT_CC_BONDLENGTH; - if(!in_TAU) TAU = DEFAULT_TAU; - if(!in_a) a = 0.01 * U / TAU; - - Domain* dalet = new Domain( - flow, bondlength, rho, d, fluid, fluid2, h, ETA, ETA2, ETAF, SIG_REF, EPS_REF, - REFMASS, muVT, nanotube, m_per_n, N, T, XI, XI2, XIF, wo_wall - ); - if(nanotube) - { - /* - * not yet implemented - */ - } - else - { - dalet->write( - prefix, a, empty, format, mu, TAU, U, original, - wo_acceleration, polarity, WLJ, symmetric, widom, x - ); - } - - return 0; -} - diff --git a/tools/standalone-generators/mkesfera/CMakeLists.txt b/tools/standalone-generators/mkesfera/CMakeLists.txt deleted file mode 100644 index 6cb95256c8..0000000000 --- a/tools/standalone-generators/mkesfera/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -cmake_minimum_required(VERSION 3.3) -message(STATUS "cmake building mkesfera") - -add_executable(mkesfera main.cpp Domain.cpp Random.cpp) -target_include_directories(mkesfera PRIVATE ${PROJECT_SOURCE_DIR}) - -include(GNUInstallDirs) -install(TARGETS mkesfera RUNTIME DESTINATION /usr/bin/) - diff --git a/tools/standalone-generators/mkesfera/Domain.cpp b/tools/standalone-generators/mkesfera/Domain.cpp deleted file mode 100644 index 712f363e41..0000000000 --- a/tools/standalone-generators/mkesfera/Domain.cpp +++ /dev/null @@ -1,197 +0,0 @@ -#include "Domain.h" -#include "Random.h" -#include - -#define DT 0.002 -#define TIME 20111102 -#define VARFRACTION 0.07 -#define BOXOVERLOAD 1.3333 -#define AVGBIN 0.025 - -void Domain::write(char* prefix, double cutoff, double T, bool do_shift, int format) -{ - std::ofstream xdr, txt, buchholz; - std::stringstream strstrm, txtstrstrm, buchholzstrstrm; - if(format == FORMAT_BRANCH) - { - strstrm << prefix << ".xdr"; - } - if(format == FORMAT_BUCHHOLZ) - { - strstrm << prefix << ".inp"; - } - xdr.open(strstrm.str().c_str(), std::ios::trunc); - if(format == FORMAT_BRANCH) - { - txtstrstrm << prefix << "_1R.txt"; - } - if(format == FORMAT_BUCHHOLZ) - { - txtstrstrm << prefix << "_1R.cfg"; - } - txt.open(txtstrstrm.str().c_str(), std::ios::trunc); - if(format == FORMAT_BUCHHOLZ) - { - buchholzstrstrm << prefix << "_1R.xml"; - buchholz.open(buchholzstrstrm.str().c_str(), std::ios::trunc); - - /* - * Gesamter Inhalt der Buchholz-Datei - */ - buchholz << "\n\n \n " - << prefix << "_1R.cfg\n \n"; - } - - unsigned fl_units; - double rhomax = (rho_i > rho_o)? rho_i: rho_o; - double N_boxes = 8.0*R_o*R_o*R_o * (BOXOVERLOAD*rhomax) / 3.0; - fl_units = ceil(pow(N_boxes, 1.0/3.0)); - double fl_unit = 2.0*R_o / (double)fl_units; - - Random* r = new Random(); - r->init( - (int)(10000.0*R_o) - (int)(3162.3*cutoff) - + (int)(1000.0*T) - (int)(316.23*rho_i) - + (int)(100.0*R_i) - ); - - bool**** fill = new bool***[fl_units]; - for (unsigned int i = 0; i < fl_units; i++) { - fill[i] = new bool**[fl_units]; - for (unsigned int j = 0; j < fl_units; j++) { - fill[i][j] = new bool*[fl_units]; - for (unsigned int k = 0; k < fl_units; k++) { - fill[i][j][k] = new bool[3]; - } - } - } - unsigned slots = 3.0*fl_units*fl_units*fl_units; - double boxdensity = (double)slots / (8.0*R_o*R_o*R_o); - std::cerr << "Box density: " << boxdensity << " (unit cell: " << fl_unit << ").\n"; - double P_in = rho_i / boxdensity; - double P_out = rho_o / boxdensity; - std::cerr << "Insertion probability: " << P_in << " inside, " << P_out << " outside.\n"; - - double goffset[3][3]; - goffset[0][0] = 0.0; goffset[1][0] = 0.5; goffset[2][0] = 0.5; - goffset[0][1] = 0.5; goffset[1][1] = 0.0; goffset[2][1] = 0.5; - goffset[0][2] = 0.5; goffset[1][2] = 0.5; goffset[2][2] = 0.0; - - unsigned N = 0; - unsigned idx[3]; - for(idx[0]=0; idx[0] < fl_units; idx[0]++) - for(idx[1]=0; idx[1] < fl_units; idx[1]++) - for(idx[2]=0; idx[2] < fl_units; idx[2]++) - for(int p=0; p < 3; p++) - { - double qq = 0.0; - for(int d = 0; d < 3; d++) - { - double qrel = (idx[d] + goffset[d][p])*fl_unit - R_o; - qq += qrel*qrel; - } - double tP = (qq > R_i*R_i)? P_out: P_in; - bool tfill = (tP >= r->rnd()); - fill[idx[0]][idx[1]][idx[2]][p] = tfill; - if(tfill) N++; - } - std::cerr << "Filling " << N << " out of " << slots << " slots (total density: " << N / (8.0*R_o*R_o*R_o) << ").\n"; - - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt.precision(6); - txt << "mardynconfig\n# \ntimestepLength\t" << DT << "\ncutoffRadius\t" << cutoff << "\nLJCutoffRadius\t" << cutoff << "\ntersoffCutoffRadius\t0.5\n"; - } - - if(format == FORMAT_BRANCH) txt << "phaseSpaceFile\t" << prefix << ".xdr\n"; - if(format == FORMAT_BUCHHOLZ) txt << "phaseSpaceFile\tOldStyle\t" << prefix << ".inp\nparallelization\tDomainDecomposition\n"; - - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt << "datastructure\tLinkedCells 1\n"; - - txt << "output\tResultWriter 500\t" << prefix << "_1R\nresultOutputTimesteps\t500\noutput\tXyzWriter 60000\t" << prefix << "_1R\noutput\tXdrWriter 60000\t" << prefix << "_1R\ninitCanonical\t10\ninitStatistics\t3003003\n"; - txt << "profile\t1 " << (unsigned)round(R_o / AVGBIN) << " 1\nprofileRecordingTimesteps\t1\nprofileOutputTimesteps\t500000\nprofiledComponent\t1\nprofileOutputPrefix\t" << prefix << "_1R\nesfera\nnomomentum 16384\nAlignCentre 333 0.003\nchemicalPotential 0 component 1 conduct " << (unsigned)round(N/3.0) << " tests every 1 steps\nWidom\n"; - } - - xdr.precision(8); - if(format == FORMAT_BRANCH) - { - xdr << "mardyn " << TIME << " tersoff\n" - << "# mardyn input file, ls1 project\n" - << "# generated by the mkesfera tool\n# \n"; - } - if(format == FORMAT_BUCHHOLZ) - { - xdr << "mardyn trunk " << TIME << "\n"; - } - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << "t\t0\nL\t" << 2.0*R_o << " " << 2.0*R_o << " " << 2.0*R_o << "\n"; - } - if(format == FORMAT_BUCHHOLZ) - { - xdr << "T\t" << T << "\n"; - xdr << "C\t1\t1 0 0 0 0\t0 0 0 1 1 1\t" << cutoff << " " << (do_shift? " 1": " 0") << "\t0 0 0 1e+10\n"; - } - if(format == FORMAT_BRANCH) - { - xdr << "C\t1\t1 0 0 0 0\t0 0 0 1 1 1\t0 0 0 1e+10\n"; - xdr << "T\t" << T << "\n"; - } - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << "N\t" << N << "\nM\tICRVQD\n"; - } - - double v = sqrt(3.0 * T); - - unsigned ID = 1; - for(idx[0]=0; idx[0] < fl_units; idx[0]++) - for(idx[1]=0; idx[1] < fl_units; idx[1]++) - for(idx[2]=0; idx[2] < fl_units; idx[2]++) - for(unsigned p=0; p < 3; p++) - { - // std::cerr << idx[0] << "\t" << idx[1] << "\t" << idx[2] << "\t|\t" << fill[idx[0]][idx[1]][idx[2]][p] << "\n"; - if(fill[idx[0]][idx[1]][idx[2]][p]) - { - double q[3]; - for(int d=0; d < 3; d++) - { - q[d] = (idx[d] + VARFRACTION*(r->rnd() - 0.5) + goffset[d][p])*fl_unit; - if(q[d] < 0.0) q[d] += 2.0*R_o; - else if(q[d] > 2.0*R_o) q[d] -= 2.0*R_o; - } - double phi = 6.283185 * r->rnd(); - double omega = 6.283185 * r->rnd(); - - xdr << ID << " " << 1 << "\t" << q[0] - << " " << q[1] << " " << q[2] - << "\t" << v*cos(phi)*cos(omega) << " " - << v*cos(phi)*sin(omega) << " " - << v*sin(phi) << "\t1 0 0 0 0 0 0\n"; - - ID++; - } - } - - xdr.close(); - txt.close(); - if(format == FORMAT_BUCHHOLZ) - { - buchholz.close(); - } - - for (unsigned int i = 0; i < fl_units; i++) { - for (unsigned int j = 0; j < fl_units; j++) { - for (unsigned int k = 0; k < fl_units; k++) { - delete[] fill[i][j][k]; - } - delete[] fill[i][j]; - } - delete[] fill[i]; - } - delete[] fill; -} - diff --git a/tools/standalone-generators/mkesfera/Domain.h b/tools/standalone-generators/mkesfera/Domain.h deleted file mode 100644 index 977366e2de..0000000000 --- a/tools/standalone-generators/mkesfera/Domain.h +++ /dev/null @@ -1,25 +0,0 @@ -#include -#include -#include -#include - - -#define FORMAT_BUCHHOLZ 0 -#define FORMAT_BRANCH 1 -#define FORMAT_BERNREUTHER 2 - -class Domain -{ - public: - Domain(double tRin, double tRout, double trhoin, double trhoout) { - R_i = tRin; R_o = tRout; rho_i = trhoin; rho_o = trhoout; - }; - - void write(char* prefix, double cutoff, double T, bool do_shift, int format); - - private: - double box[3]; - - double R_i, R_o, rho_i, rho_o; -}; - diff --git a/tools/standalone-generators/mkesfera/Random.cpp b/tools/standalone-generators/mkesfera/Random.cpp deleted file mode 100644 index 2a96b04a08..0000000000 --- a/tools/standalone-generators/mkesfera/Random.cpp +++ /dev/null @@ -1,30 +0,0 @@ -#include "Random.h" - -Random::Random() { this->init(8624); } - -void Random::init(int seed) -{ - this->ix_muVT = (seed ^ (int)888889999) | (int)1; - this->iy_muVT = seed ^ (int)777755555; - // Calculate normalization factor - this->am_muVT = 2.0 / (1.0 + (unsigned)((int)-1)); -} - -float Random::rnd() -{ - float rnd; - const int IA = 16807; - const int IM = 2147483647; - const int IQ = 127773; - const int IR = 2836; - - ix_muVT ^= (ix_muVT >> 13); - ix_muVT ^= (ix_muVT << 17); - ix_muVT ^= (ix_muVT >> 5); - - int k = iy_muVT / IQ; - iy_muVT = IA * (iy_muVT - k*IQ) - IR*k; - if(iy_muVT < 0) iy_muVT += IM; - rnd = am_muVT * ((IM & (ix_muVT ^ iy_muVT)) | (int)1); - return rnd; -} diff --git a/tools/standalone-generators/mkesfera/Random.h b/tools/standalone-generators/mkesfera/Random.h deleted file mode 100644 index 30cfe95573..0000000000 --- a/tools/standalone-generators/mkesfera/Random.h +++ /dev/null @@ -1,15 +0,0 @@ -class Random -{ - public: - Random(); - void init(int seed); - - float rnd(); - - int getIX() { return this->ix_muVT; } - - private: - int ix_muVT, iy_muVT; - float am_muVT; -}; - diff --git a/tools/standalone-generators/mkesfera/main.cpp b/tools/standalone-generators/mkesfera/main.cpp deleted file mode 100644 index 2a3128fb6e..0000000000 --- a/tools/standalone-generators/mkesfera/main.cpp +++ /dev/null @@ -1,106 +0,0 @@ -#include "Domain.h" - -#include -#include -#include -#include - - -int main(int argc, char** argv) -{ - const char* usage = "usage: mkesfera [-e] -I -i -O -o [-R ] [-r] [-S] -T [-U] [-u]\n\n-e\tuse B-e-rnreuther format\n-r\tuse b-r-anch format (active by default)\n-S\tshift (active by default)\n-U\tunshift\n-u\tuse B-u-chholz format\n"; - if((argc < 12) || (argc > 16)) - { - std::cout << "There are " << argc - << " arguments where 12 to 16 should be given.\n\n"; - std::cout << usage; - return 1; - } - - bool do_shift = true; - unsigned format = FORMAT_BRANCH; - - double cutoff = 2.5; - double rho = 0.319; - double rho2 = 0.319; - double R = 7.0; - double R2 = 14.0; - double T = 1.0; - - char* prefix = argv[1]; - - for(int i=2; i < argc; i++) - { - if(*argv[i] != '-') - { - std::cout << "Flag expected where '" << argv[i] - << "' was given.\n\n"; - std::cout << usage; - return 2; - } - for(int j=1; argv[i][j]; j++) - { - if(argv[i][j] == 'e') format = FORMAT_BERNREUTHER; - else if(argv[i][j] == 'I') - { - i++; - R = atof(argv[i]); - break; - } - else if(argv[i][j] == 'i') - { - i++; - rho = atof(argv[i]); - break; - } - else if(argv[i][j] == 'O') - { - i++; - R2 = atof(argv[i]); - break; - } - else if(argv[i][j] == 'o') - { - i++; - rho2 = atof(argv[i]); - break; - } - else if(argv[i][j] == 'R') - { - i++; - cutoff = atof(argv[i]); - break; - } - else if(argv[i][j] == 'r') format = FORMAT_BRANCH; - else if(argv[i][j] == 'S') do_shift = true; - else if(argv[i][j] == 'T') - { - i++; - T = atof(argv[i]); - break; - } - else if(argv[i][j] == 'U') do_shift = false; - else if(argv[i][j] == 'u') format = FORMAT_BUCHHOLZ; - else - { - std::cout << "Invalid flag '-" << argv[i][j] - << "' was detected.\n\n" << usage; - return 2; - } - } - } - - if(format == FORMAT_BERNREUTHER) - { - std::cout << "B-e-rnreuther format (flag -e) " - << "is unavailable at present.\n\n" << usage; - return 3; - } - - Domain* dalet; - dalet = new Domain(R, R2, rho, rho2); - dalet->write(prefix, cutoff, T, do_shift, format); - - return 0; -} - diff --git a/tools/standalone-generators/mktcts/CMakeLists.txt b/tools/standalone-generators/mktcts/CMakeLists.txt deleted file mode 100644 index 99e83b76f8..0000000000 --- a/tools/standalone-generators/mktcts/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -cmake_minimum_required(VERSION 3.3) -message(STATUS "cmake building mktcts") - -add_executable(mktcts main.cpp Domain.cpp Random.cpp) -target_include_directories(mktcts PRIVATE ${PROJECT_SOURCE_DIR}) - -include(GNUInstallDirs) -install(TARGETS mktcts RUNTIME DESTINATION /usr/bin/) - diff --git a/tools/standalone-generators/mktcts/Domain.cpp b/tools/standalone-generators/mktcts/Domain.cpp deleted file mode 100644 index 8ff4a6de08..0000000000 --- a/tools/standalone-generators/mktcts/Domain.cpp +++ /dev/null @@ -1,304 +0,0 @@ -#include "Domain.h" -#include "Random.h" -#include - -#define BINS 512 -#define BINS_AUTOCORR 16 -#define DT 0.0025 -#define PRECISION 5 -#define TIME 20150730 -#define VARFRACTION 0.125 - -Domain::Domain(double t_h, unsigned t_N, double t_rho, double t_rho2) -{ - this->gradient = true; - this->rho = t_rho; - this->rho2 = t_rho2; - this->RDF = 0.0; - - double V = 2.0 * t_N / (rho + rho2); - this->box[0] = sqrt(V / t_h); - this->box[1] = t_h; - this->box[2] = sqrt(V / t_h); - - this->use_hato = false; -} - -Domain::Domain(unsigned t_N, double t_rho, double t_RDF) -{ - this->gradient = false; - this->rho = t_rho; - this->rho2 = t_rho; - this->RDF = t_RDF; - - this->box[0] = pow((double)t_N/rho, 1.0/3.0); - this->box[1] = this->box[0]; - this->box[2] = this->box[0]; - - this->use_hato = false; -} - -void Domain::write(char* prefix, double cutoff, double mu, double T, bool do_shift, bool use_mu, bool compute_autocorr, int format) -{ - std::ofstream xdr, txt, buchholz; - std::stringstream strstrm, txtstrstrm, buchholzstrstrm; - if(format == FORMAT_BRANCH) - { - strstrm << prefix << ".xdr"; - } - if(format == FORMAT_BUCHHOLZ) - { - strstrm << prefix << ".inp"; - } - xdr.open(strstrm.str().c_str(), std::ios::trunc); - if(format == FORMAT_BRANCH) - { - txtstrstrm << prefix << "_1R.txt"; - } - if(format == FORMAT_BUCHHOLZ) - { - txtstrstrm << prefix << "_1R.cfg"; - } - txt.open(txtstrstrm.str().c_str(), std::ios::trunc); - if(format == FORMAT_BUCHHOLZ) - { - buchholzstrstrm << prefix << "_1R.xml"; - buchholz.open(buchholzstrstrm.str().c_str(), std::ios::trunc); - - /* - * Gesamter Inhalt der Buchholz-Datei - */ - buchholz << "\n\n \n " - << prefix << "_1R.cfg\n \n"; - } - - unsigned fl_units[3][2]; - double fl_unit[3][2]; - double N_id[2]; - for(int i=0; i < 2; i++) - { - N_id[i] = box[0]*(0.5*box[1])*box[2] * ((i == 0)? rho: rho2); - double N_boxes = N_id[i] / 3.0; - fl_units[1][i] = round( - pow( - (N_boxes * (0.5*box[1]) * (0.5*box[1])) - / (this->box[0] * this->box[2]), 1.0/3.0 - ) - ); - if(fl_units[1][i] == 0) fl_units[1][i] = 1; - double bxbz_id = N_boxes / fl_units[1][i]; - fl_units[0][i] = round(sqrt(this->box[0] * bxbz_id / this->box[2])); - if(fl_units[0][i] == 0) fl_units[0][i] = 1; - fl_units[2][i] = ceil(bxbz_id / fl_units[0][i]); - for(int d=0; d < 3; d++) fl_unit[d][i] = ((d == 1)? 0.5: 1.0) * box[d] / (double)fl_units[d][i]; - std::cout << "Elementary cell " << i << ": " << fl_unit[0][i] << " x " << fl_unit[1][i] << " x " << fl_unit[2][i] << ".\n"; - } - - Random* r = new Random(); - r->init( - (int)(10000.0*box[0]) - (int)(3162.3*cutoff) - + (int)(1000.0*T) - (int)(316.23*mu) - + (int)(100.0*box[1]) - ); - - bool fill0[fl_units[0][0]][fl_units[1][0]][fl_units[2][0]][3]; - bool fill1[fl_units[0][1]][fl_units[1][1]][fl_units[2][1]][3]; - unsigned N[2]; - unsigned slots[2]; - for(unsigned l=0; l < 2; l++) - { - for(unsigned i=0; i < fl_units[0][l]; i++) - for(unsigned j=0; j < fl_units[1][l]; j++) - for(unsigned k=0; k < fl_units[2][l]; k++) - for(unsigned d=0; d < 3; d++) - { - if(l == 0) fill0[i][j][k][d] = true; - else fill1[i][j][k][d] = true; - } - slots[l] = 3 * fl_units[0][l] * fl_units[1][l] * fl_units[2][l]; - N[l] = slots[l]; - } - bool tswap; - double pswap; - for(unsigned l=0; l < 2; l++) - { - for(unsigned m=0; m < PRECISION; m++) - { - tswap = (N[l] < N_id[l]); - pswap = (N_id[l] - (double)N[l]) / ((tswap? slots[l]: 0) - (double)N[l]); - // std::cout << "N = " << N[l] << ", N_id = " << N_id[l] << " => tswap = " << tswap << ", pswap = " << pswap << "\n"; - for(unsigned i=0; i < fl_units[0][l]; i++) - for(unsigned j=0; j < fl_units[1][l]; j++) - for(unsigned k=0; k < fl_units[2][l]; k++) - for(unsigned d=0; d < 3; d++) - if(pswap >= r->rnd()) - { - if(((l == 0) && fill0[i][j][k][d]) || ((l == 1) && fill1[i][j][k][d])) N[l] --; - if(l == 0) fill0[i][j][k][d] = tswap; - else fill1[i][j][k][d] = tswap; - if(tswap) N[l] ++; - } - } - std::cout << "Filling " << N[l] << " of 3*" - << fl_units[0][l] << "*" << fl_units[1][l] << "*" << fl_units[2][l] - << " = " << slots[l] << " slots (ideally " << N_id[l] << ").\n"; - } - - txt.precision(6); - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt << "mardynconfig\ntimestepLength\t" << DT << "\ncutoffRadius\t" - << ((RDF > cutoff)? 1.031623*RDF: cutoff) << "\nLJCutoffRadius\t" << cutoff << "\n"; - } - if(format == FORMAT_BRANCH) - { - txt << "tersoffCutoffRadius\t0.5\nphaseSpaceFile\t" << prefix << ".xdr\n"; - } - if(format == FORMAT_BUCHHOLZ) - { - txt << "phaseSpaceFile\tOldStyle\t" << prefix << ".inp\nparallelization\tDomainDecomposition\n"; - } - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt << "datastructure\tLinkedCells\t1\n"; - - txt << "output\tResultWriter\t1500\t" << prefix << "_1R\n"; - } - if(format == FORMAT_BRANCH) txt << "resultOutputTimesteps\t1500\n"; - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - txt << "initCanonical\t10\ninitStatistics\t100000\n"; - - if(compute_autocorr) - { - double rho_max = (rho > rho2)? rho: rho2; - double rho_min = (rho > rho2)? rho2: rho; - txt << "activatePsi\t" << rho_min * (0.001 + rho_min/(rho_max + 9.0*rho_min)) + 0.000001 << "\nautocorrProfileUnits\t" << ((gradient || this->use_hato)? BINS_AUTOCORR: 1) << "\nautocorrMinStep\t2\nautocorrOutReset\t32768 131072\nautocorrArraySize\t32\nautocorrLevels\t3 16\n"; - } - if(gradient || this->use_hato) - { - txt << "output\tXyzWriter 400\t" << prefix << "_1R\nprofile\t1 " << BINS << " 1\nprofileRecordingTimesteps\t1\nprofileOutputTimesteps\t200000\nprofiledComponent\t1\nprofileOutputPrefix\t" << prefix << "_1R\n"; - if(!this->use_hato && !compute_autocorr) txt << "AlignCentre\t100 0.01\n"; - } - else - { - txt << "output\tXyzWriter 40000\t" << prefix << "_1R\nRDF\t" << RDF/(double)BINS << " " << BINS << "\nRDFOutputTimesteps\t500000\nRDFOutputPrefix\t" << prefix << "_1R\n"; - } - txt << "nomomentum\t1024\n"; - - if(this->use_hato) - { - txt << "planckConstant\t" << sqrt(6.2831853 * T) << "\ninitGrandCanonical\t124816\n"; - double p_high = (p1 > p2)? p1: p2; - double p_low = (p1 > p2)? p2: p1; - double base_pressure = (p1 > 0)? p1: -p1; - if(p2 > base_pressure) base_pressure = p2; - else if(-p2 > base_pressure) base_pressure = -p2; - double coupling = 4.0 * (mu_high - mu_low) / (p_high - p_low + base_pressure); - txt << "Hatonian target PTOTAL " << p1 << " chemicalPotential " << mu_low << " " << mu_high << " coupling " << coupling << " " << coupling << " component 1 control 0.0 " << box[1]/6.0 << " 0.0 to " << box[0] << " " << box[1]/3.0 << " " << box[2] << " conduct " << (int)round(N[0] / 500.0) << " tests every 8 steps\n"; - txt << "Hatonian target PTOTAL " << p2 << " chemicalPotential " << mu_low << " " << mu_high << " coupling " << coupling << " " << coupling << " component 1 control 0.0 " << 2.0*box[1]/3.0 << " 0.0 to " << box[0] << " " << 5.0*box[1]/6.0 << " " << box[2] << " conduct " << (int)round(N[1] / 500.0) << " tests every 8 steps\n"; - txt << "AccumulatorSize\t512\n"; - } - else if(use_mu) - { - if(gradient) - { - txt << "chemicalPotential\t" << mu << " component 1\tcontrol 0 0 0 to " << box[0] << " " << 0.1*box[1] << " " << box[2] << "\tconduct " << (int)round(N[0] / 2000.0) << " tests every 2 steps\n"; - txt << "chemicalPotential\t" << mu << " component 1\tcontrol 0 " << 0.5*box[1] << " 0 to " << box[0] << " " << 0.6*box[1] << " " << box[2] << "\tconduct " << (int)round(N[1] / 2000.0) << " tests every 2 steps\n"; - } - else - { - std::cout << "N[0] = " << N[0] << ", N[1] = " << N[1] << ", conducting " << round((double)(N[0]+N[1]) / 5000.0) << " test actions.\n"; - txt << "chemicalPotential\t" << mu << " component 1\tconduct " << (int)round((double)(N[0]+N[1]) / 5000.0) << " tests every 2 steps\n"; - } - txt << "planckConstant\t" << sqrt(6.2831853 * T) << "\ninitGrandCanonical\t60000\n"; - } - } - - xdr.precision(8); - if(format == FORMAT_BRANCH) - { - xdr << "mardyn " << TIME << " tersoff\n" - << "# mardyn input file, ls1 project\n" - << "# generated by the mkTcTS tool\n# \n" - << "t\t0\nT\t" << T << "\nL\t" << box[0] << " " << box[1] - << " " << box[2] << "\nC\t1\t1 0 0 0 0\t0 0 0 1 1 1\t0 0 0 1e+10\nN\t" - << N[0]+N[1] << "\nM\tICRVQD\n\n"; - } - if(format == FORMAT_BUCHHOLZ) - { - xdr << "mardyn trunk " << TIME << "\n" - << "t\t0\nT\t" << T << "\nL\t" << box[0] << " " << box[1] - << " " << box[2] << "\nC\t1\t1 0 0 0 0\t0 0 0 1 1 1\t" << cutoff << " " - << (do_shift? "1": "0") << "\t0 0 0 1e+10\nN\t" - << N[0]+N[1] << "\nM\tICRVQD\n\n"; - } - - double v = sqrt(3.0 * T); - double loffset[3][2]; - if(this->use_hato) - { - loffset[0][0] = 0.1; loffset[1][0] = 0.0; loffset[2][0] = 0.1; - loffset[0][1] = 0.1; loffset[1][1] = 0.5; loffset[2][1] = 0.1; - } - else - { - loffset[0][0] = 0.1; loffset[1][0] = 0.3; loffset[2][0] = 0.1; - loffset[0][1] = 0.1; loffset[1][1] = 0.8; loffset[2][1] = 0.1; - } - double goffset[3][3]; - goffset[0][0] = 0.0; goffset[1][0] = 0.5; goffset[2][0] = 0.5; - goffset[0][1] = 0.5; goffset[1][1] = 0.0; goffset[2][1] = 0.5; - goffset[0][2] = 0.5; goffset[1][2] = 0.5; goffset[2][2] = 0.0; - - unsigned ID = 1; - for(unsigned l=0; l < 2; l++) - for(unsigned i=0; i < fl_units[0][l]; i++) - for(unsigned j=0; j < fl_units[1][l]; j++) - for(unsigned k=0; k < fl_units[2][l]; k++) - for(unsigned d=0; d < 3; d++) - { - if(((l == 0) && fill0[i][j][k][d]) || ((l == 1) && fill1[i][j][k][d])) - { - double q[3]; - q[0] = i * fl_unit[0][l]; - q[1] = j * fl_unit[1][l]; - q[2] = k * fl_unit[2][l]; - for(int m=0; m < 3; m++) - { - q[m] += box[m]*loffset[m][l] + fl_unit[m][l]*goffset[m][d]; - q[m] += VARFRACTION * fl_unit[m][l] * (r->rnd() - 0.5); - if(q[m] > box[m]) q[m] -= box[m]; - else if(q[m] < 0.0) q[m] += box[m]; - } - double phi = 6.283185 * r->rnd(); - double omega = 6.283185 * r->rnd(); - - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << ID << " " << 1 << "\t" << q[0] - << " " << q[1] << " " << q[2] - << "\t" << v*cos(phi)*cos(omega) << " " - << v*cos(phi)*sin(omega) << " " - << v*sin(phi) << "\t1 0 0 0 0 0 0\n"; - } - ID++; - } - else - { - if((format == FORMAT_BRANCH) || (format == FORMAT_BUCHHOLZ)) - { - xdr << "\n"; - } - } - } - - xdr.close(); - txt.close(); - if(format == FORMAT_BUCHHOLZ) - { - buchholz.close(); - } -} - diff --git a/tools/standalone-generators/mktcts/Domain.h b/tools/standalone-generators/mktcts/Domain.h deleted file mode 100644 index aa97e8bb69..0000000000 --- a/tools/standalone-generators/mktcts/Domain.h +++ /dev/null @@ -1,36 +0,0 @@ -#include -#include -#include -#include - - -#define FORMAT_BUCHHOLZ 0 -#define FORMAT_BRANCH 1 -#define FORMAT_BERNREUTHER 2 - -class Domain -{ - public: - Domain( - double t_h, unsigned t_N, double t_rho, double t_rho2 - ); - Domain( - unsigned t_N, double t_rho, double t_RDF - ); - void write( - char* prefix, double cutoff, double mu, double T, bool do_shift, bool use_mu, bool compute_autocorr, int format - ); - void hato(double in1, double in2, double inc1, double inc2) { use_hato = true; p1 = in1; p2 = in2; mu_low = inc1; mu_high = inc2; } - - private: - double box[3]; - - bool gradient; - double rho; - double rho2; - double RDF; - - bool use_hato; - double p1, p2; - double mu_low, mu_high; -}; diff --git a/tools/standalone-generators/mktcts/Random.cpp b/tools/standalone-generators/mktcts/Random.cpp deleted file mode 100644 index 2a96b04a08..0000000000 --- a/tools/standalone-generators/mktcts/Random.cpp +++ /dev/null @@ -1,30 +0,0 @@ -#include "Random.h" - -Random::Random() { this->init(8624); } - -void Random::init(int seed) -{ - this->ix_muVT = (seed ^ (int)888889999) | (int)1; - this->iy_muVT = seed ^ (int)777755555; - // Calculate normalization factor - this->am_muVT = 2.0 / (1.0 + (unsigned)((int)-1)); -} - -float Random::rnd() -{ - float rnd; - const int IA = 16807; - const int IM = 2147483647; - const int IQ = 127773; - const int IR = 2836; - - ix_muVT ^= (ix_muVT >> 13); - ix_muVT ^= (ix_muVT << 17); - ix_muVT ^= (ix_muVT >> 5); - - int k = iy_muVT / IQ; - iy_muVT = IA * (iy_muVT - k*IQ) - IR*k; - if(iy_muVT < 0) iy_muVT += IM; - rnd = am_muVT * ((IM & (ix_muVT ^ iy_muVT)) | (int)1); - return rnd; -} diff --git a/tools/standalone-generators/mktcts/Random.h b/tools/standalone-generators/mktcts/Random.h deleted file mode 100644 index 30cfe95573..0000000000 --- a/tools/standalone-generators/mktcts/Random.h +++ /dev/null @@ -1,15 +0,0 @@ -class Random -{ - public: - Random(); - void init(int seed); - - float rnd(); - - int getIX() { return this->ix_muVT; } - - private: - int ix_muVT, iy_muVT; - float am_muVT; -}; - diff --git a/tools/standalone-generators/mktcts/main.cpp b/tools/standalone-generators/mktcts/main.cpp deleted file mode 100644 index eb6302d22f..0000000000 --- a/tools/standalone-generators/mktcts/main.cpp +++ /dev/null @@ -1,158 +0,0 @@ -#include "Domain.h" - -#include -#include -#include -#include - - -int main(int argc, char** argv) -{ - const char* usage = "usage: mkTcTS -c [-a] [-d ] [-e] [-H ] [-h ] [-m ] -N [-p ] [-R ] [-r] [-S] -T [-U] [-u]\n\n-a\tcompute velocity autocorrelation\n-e\tuse B-e-rnreuther format\n-P\tcompute pseudomomenta\n-r\tuse b-r-anch format\n-S\tshift (active by default)\n-U\tunshift\n-u\tuse B-u-chholz format (active by default)\n"; - if((argc < 8) || (argc > 28)) - { - std::cout << "There are " << argc - << " arguments where 8 to 28 should be given.\n\n"; - std::cout << usage; - return 1; - } - - bool compute_autocorr = false; - bool do_shift = true; - bool in_h = false; - bool use_mu = false; - bool gradient = false; - unsigned format = FORMAT_BUCHHOLZ; - - double cutoff = 2.5; - double RDF = 12.0; - double h = 0.0; - double mu = 0.0; - unsigned N = 131072; - double rho = 0.319; - double rho2 = 0.319; - double T = 1.0779; - - bool use_hato = false; - double p1 = 0.0; - double p2 = 0.0; - double mu_low = 0.0; - double mu_high = 0.0; - - char* prefix = argv[1]; - - for(int i=2; i < argc; i++) - { - if(*argv[i] != '-') - { - std::cout << "Flag expected where '" << argv[i] - << "' was given.\n\n"; - std::cout << usage; - return 2; - } - for(int j=1; argv[i][j]; j++) - { - if(argv[i][j] == 'a') compute_autocorr = true; - else if(argv[i][j] == 'c') - { - i++; - rho = atof(argv[i]); - break; - } - else if(argv[i][j] == 'd') - { - gradient = true; - i++; - rho2 = atof(argv[i]); - break; - } - else if(argv[i][j] == 'e') format = FORMAT_BERNREUTHER; - else if(argv[i][j] == 'H') - { - use_hato = true; - i++; - p1 = atof(argv[i]); - i++; - p2 = atof(argv[i]); - i++; - mu_low = atof(argv[i]); - i++; - mu_high = atof(argv[i]); - break; - } - else if(argv[i][j] == 'h') - { - in_h = true; - i++; - h = atof(argv[i]); - break; - } - else if(argv[i][j] == 'm') - { - use_mu = true; - i++; - mu = atof(argv[i]); - break; - } - else if(argv[i][j] == 'N') - { - i++; - N = atoi(argv[i]); - break; - } - else if(argv[i][j] == 'p') - { - i++; - RDF = atof(argv[i]); - break; - } - else if(argv[i][j] == 'R') - { - i++; - cutoff = atof(argv[i]); - break; - } - else if(argv[i][j] == 'r') format = FORMAT_BRANCH; - else if(argv[i][j] == 'S') do_shift = true; - else if(argv[i][j] == 'T') - { - i++; - T = atof(argv[i]); - break; - } - else if(argv[i][j] == 'U') do_shift = false; - else if(argv[i][j] == 'u') format = FORMAT_BUCHHOLZ; - else - { - std::cout << "Invalid flag '-" << argv[i][j] - << "' was detected.\n\n" << usage; - return 2; - } - } - } - - if(in_h && !gradient) - { - std::cout << "The box dimension can only be specified for " - << "systems with a density gradient.\n\n" << usage; - return 5; - } - - if(format == FORMAT_BERNREUTHER) - { - std::cout << "B-e-rnreuther format (flag -e) " - << "is unavailable at present.\n\n" << usage; - return 3; - } - - if(!in_h) h = pow((double)N/rho, 1.0/3.0); - - Domain* dalet; - if(gradient) dalet = new Domain(h, N, rho, rho2); - else dalet = new Domain(N, rho, RDF); - if(use_hato) dalet->hato(p1, p2, mu_low, mu_high); - dalet->write(prefix, cutoff, mu, T, do_shift, use_mu, compute_autocorr, format); - - return 0; -} - diff --git a/tools/standalone-generators/produce b/tools/standalone-generators/produce deleted file mode 100755 index 67d182c2ae..0000000000 --- a/tools/standalone-generators/produce +++ /dev/null @@ -1,2 +0,0 @@ -cmake .; make package; make - diff --git a/tools/standalone-generators/uninstall b/tools/standalone-generators/uninstall deleted file mode 100755 index 801a352d9f..0000000000 --- a/tools/standalone-generators/uninstall +++ /dev/null @@ -1,2 +0,0 @@ -sudo dpkg -r ls1-standalone-generators -