From 3a7c10d169ed7bf98160ac0753f4510a07b3e6e6 Mon Sep 17 00:00:00 2001 From: Zhi Wang Date: Wed, 31 Aug 2022 14:56:10 -0500 Subject: [PATCH] add dcd support --- include/ff/atom.h | 4 - include/ff/rwcrd.h | 38 ++++ include/md/misc.h | 1 - include/tinker9.h | 4 +- src/atom.cpp | 163 -------------- src/cmakesrc.txt | 1 + src/mdintg.cpp | 31 ++- src/mdsave.cpp | 83 ------- src/rwcrd.cpp | 540 +++++++++++++++++++++++++++++++++++++++++++++ src/xanalyze.cpp | 22 +- src/xbar.cpp | 32 ++- src/xtestgrad.cpp | 9 +- test/bounds.cpp | 7 +- 13 files changed, 638 insertions(+), 297 deletions(-) create mode 100644 include/ff/rwcrd.h create mode 100644 src/rwcrd.cpp diff --git a/include/ff/atom.h b/include/ff/atom.h index 4f541cce0..cbbfd42ed 100644 --- a/include/ff/atom.h +++ b/include/ff/atom.h @@ -29,10 +29,6 @@ void copyPosToXyz(bool refreshNBList); /// - Tinker uses centers of mass. void bounds(); -void readFrameOpen(const std::string& filename, std::ifstream& input); -void readFrameCopyinToXyz(std::ifstream& input, bool& done); -void readFrameClose(std::ifstream& input); - //====================================================================// // // // Global Variables // diff --git a/include/ff/rwcrd.h b/include/ff/rwcrd.h new file mode 100644 index 000000000..80d55d3fe --- /dev/null +++ b/include/ff/rwcrd.h @@ -0,0 +1,38 @@ +#pragma once +#include + +namespace tinker { +enum class CrdFormat +{ + NONE, + TXYZ1, + TXYZ2_PBC, + DCD +}; + +class CrdR; +class CrdReader +{ +protected: + CrdR* m_impl; + +public: + ~CrdReader(); + CrdReader(std::string crdfile, CrdFormat crdformat = CrdFormat::NONE); + int readCurrent(); +}; + +class CrdW; +class CrdWriter +{ +protected: + CrdW* m_impl; + const double *qx, *qy, *qz; + +public: + ~CrdWriter(); + CrdWriter(const double* xx, const double* yy, const double* zz, // + std::string crdfile, CrdFormat crdformat = CrdFormat::NONE); + int writeCurrent(); +}; +} diff --git a/include/md/misc.h b/include/md/misc.h index 206ba90f3..9be994b55 100644 --- a/include/md/misc.h +++ b/include/md/misc.h @@ -68,7 +68,6 @@ const TimeScaleConfig& respaTSConfig(); /// \ingroup md void mdsaveAsync(int istep, time_prec dt); -void mdDebugSaveSync(); /// \ingroup md void mdsaveSynchronize(); /// \ingroup md diff --git a/include/tinker9.h b/include/tinker9.h index 15f393226..a92d3e3b1 100644 --- a/include/tinker9.h +++ b/include/tinker9.h @@ -8,7 +8,7 @@ /// \def TINKER9_VERSION_MINOR /// \def TINKER9_VERSION_PATCH #define TINKER9_VERSION_MAJOR 1 -#define TINKER9_VERSION_MINOR 1 +#define TINKER9_VERSION_MINOR 2 #define TINKER9_VERSION_PATCH 0 /// \} @@ -31,7 +31,7 @@ " ### ### ""\n" \ " ### Tinker9 -- Software Tools for Molecular Design ###""\n" \ " ## ##""\n" \ -" ## Version 1.1.0 August 2022 ##""\n" \ +" ## Version 1.2.0 september 2022 ##""\n" \ " ## ##""\n" \ " ## Copyright (c) Zhi Wang & the Ponder Lab ##""\n" \ " ### All Rights Reserved ###""\n" \ diff --git a/src/atom.cpp b/src/atom.cpp index 15964871e..56150f158 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -154,167 +154,4 @@ void bounds() TINKER_FCALL2(acc1, cu1, bounds); copyPosToXyz(); } - -inline namespace v1 { -enum -{ - DCD_HEADER = 0, - DCD_TDELTA = 10, - DCD_USEBOX = 11, - DCD_CTRL_LEN = 21, - - DCD_TITLE_NCHARS = 80, - - DCD_AX = 0, - DCD_COS_G = 1, - DCD_BX = 2, - DCD_COS_B = 3, - DCD_COS_A = 4, - DCD_CX = 5, - DCD_XTAL_LEN = 6, -}; - -enum class Archive -{ - NONE = 0, - XYZ = 1, - DCD = 2, -}; -} - -static int dcdControl[DCD_CTRL_LEN] = {0}; -static std::vector dcdx, dcdy, dcdz; -static Archive archive = Archive::NONE; - -static void dcdReadIntoBuffer(void* buffer, int nbyte, std::ifstream& ipt) -{ - int size1, size2; - ipt.read((char*)&size1, sizeof(int)); - if (nbyte > 0) assert(nbyte == size1); - ipt.read((char*)buffer, size1); - ipt.read((char*)&size2, sizeof(int)); -} - -void readFrameOpen(const std::string& filename, std::ifstream& ipt) -{ - // get file format type by inspection of first character - char a1; - ipt.open(filename); - ipt >> a1; - auto arc = Archive::NONE; - if (a1 == ' ') - arc = Archive::XYZ; - else if ('0' <= a1 and a1 <= '9') - arc = Archive::XYZ; - else - arc = Archive::DCD; - - if (arc == Archive::DCD) { - ipt.close(); - ipt.open(filename, std::ios::in | std::ios::binary); - - // read header info along with title and number of atoms - dcdReadIntoBuffer(dcdControl, sizeof(int) * DCD_CTRL_LEN, ipt); - - int dcdTitleRecordLen; - ipt.read((char*)&dcdTitleRecordLen, sizeof(int)); - std::vector titlebuf; - titlebuf.resize(dcdTitleRecordLen + sizeof(int)); - ipt.read(titlebuf.data(), dcdTitleRecordLen + sizeof(int)); - - int dcdNAtom; - dcdReadIntoBuffer(&dcdNAtom, sizeof(int), ipt); - assert(n == dcdNAtom); - dcdx.resize(n); - dcdy.resize(n); - dcdz.resize(n); - } - - archive = arc; -} - -static void readFrameDCD(std::ifstream& ipt) -{ - if (dcdControl[DCD_USEBOX]) { - double dcdXtal[DCD_XTAL_LEN]; - dcdReadIntoBuffer(dcdXtal, sizeof(double) * DCD_XTAL_LEN, ipt); - double ax = dcdXtal[DCD_AX], bx = dcdXtal[DCD_BX], cx = dcdXtal[DCD_CX]; - double al = 90., be = 90., ga = 90.; - if (dcdXtal[DCD_COS_A] != 0.0) al = std::acos(dcdXtal[DCD_COS_A]) * radian; - if (dcdXtal[DCD_COS_B] != 0.0) be = std::acos(dcdXtal[DCD_COS_B]) * radian; - if (dcdXtal[DCD_COS_G] != 0.0) ga = std::acos(dcdXtal[DCD_COS_G]) * radian; - Box p; - boxLattice(p, box_shape, ax, bx, cx, al, be, ga); - boxSetCurrent(p); - } - - dcdReadIntoBuffer(dcdx.data(), sizeof(float) * n, ipt); - dcdReadIntoBuffer(dcdy.data(), sizeof(float) * n, ipt); - dcdReadIntoBuffer(dcdz.data(), sizeof(float) * n, ipt); - for (int i = 0; i < n; ++i) { - atoms::x[i] = dcdx[i]; - atoms::y[i] = dcdy[i]; - atoms::z[i] = dcdz[i]; - } -} - -static void readFrameXYZ(std::ifstream& ipt) -{ - std::string line; - std::getline(ipt, line); // n and title - std::getline(ipt, line); // either box size or first atom - // 18.643000 18.643000 18.643000 90.000000 90.000000 90.000000 - // 1 O 8.733783 7.084710 -0.688468 1 2 3 - double l1, l2, l3, a1, a2, a3; - int matched = std::sscanf(line.data(), "%lf%lf%lf%lf%lf%lf", &l1, &l2, &l3, &a1, &a2, &a3); - int row = 0; - int index; - char name[32]; - double xr, yr, zr; - if (matched == 6) { - Box p; - boxLattice(p, box_shape, l1, l2, l3, a1, a2, a3); - boxSetCurrent(p); - } else { - std::sscanf(line.data(), "%d%s%lf%lf%lf", &index, name, &xr, &yr, &zr); - index -= 1; - atoms::x[index] = xr; - atoms::y[index] = yr; - atoms::z[index] = zr; - row = 1; - } - - for (int ir = row; ir < n; ++ir) { - std::getline(ipt, line); - std::sscanf(line.data(), "%d%s%lf%lf%lf", &index, name, &xr, &yr, &zr); - index -= 1; - atoms::x[index] = xr; - atoms::y[index] = yr; - atoms::z[index] = zr; - } -} - -void readFrameCopyinToXyz(std::ifstream& ipt, bool& done) -{ - if (!ipt) done = true; - if (done) return; - - if (archive == Archive::XYZ) - readFrameXYZ(ipt); - else if (archive == Archive::DCD) - readFrameDCD(ipt); - - xyzData(RcOp::INIT); - - if (!ipt.good() or ipt.peek() == EOF) done = true; -} - -void readFrameClose(std::ifstream& ipt) -{ - ipt.close(); - dcdx.clear(); - dcdy.clear(); - dcdz.clear(); - archive = Archive::NONE; -} } diff --git a/src/cmakesrc.txt b/src/cmakesrc.txt index 751fcd664..67d6df9ff 100644 --- a/src/cmakesrc.txt +++ b/src/cmakesrc.txt @@ -72,6 +72,7 @@ pme.cpp potent.cpp rattle.cpp rcman.cpp +rwcrd.cpp spatial.cpp switch.cpp tinkersuppl.cpp diff --git a/src/mdintg.cpp b/src/mdintg.cpp index f54079dfb..6e117163c 100644 --- a/src/mdintg.cpp +++ b/src/mdintg.cpp @@ -1,20 +1,25 @@ +#include "ff/rwcrd.h" #include "math/pow2.h" #include "md/integrator.h" #include "md/misc.h" #include "md/pq.h" +#include "tool/argkey.h" +#include "tool/darray.h" #include "tool/error.h" #include "tool/externfunc.h" #include "tool/iofortstr.h" -#include +#include "tool/tinkersuppl.h" +#include #include #include #include +#include + namespace tinker { void mdData(RcOp op) { - if (not(calc::md & rc_flag)) - return; + if (not(calc::md & rc_flag)) return; RcMan intg42{mdIntegrateData, op}; RcMan save42{mdsaveData, op}; @@ -147,11 +152,31 @@ void mdrestPrintP1(bool prints, double vtot1, double vtot2, double vtot3, double void mdPropagate(int nsteps, time_prec dt_ps) { + bool useDebugArcFile = false; + getKV("T9-DBG-ARCHIVE", useDebugArcFile, false); + std::string dbgfile = ""; + std::vector qx, qy, qz; + if (useDebugArcFile) { + dbgfile = FstrView(files::filename)(1, files::leng).trim() + ".dbg"; + dbgfile = tinker_f_version(dbgfile, "old"); + qx.resize(n); + qy.resize(n); + qz.resize(n); + } + auto iarcdbg = CrdWriter(qx.data(), qy.data(), qz.data(), dbgfile); + for (int istep = 1; istep <= nsteps; ++istep) { intg->dynamic(istep, dt_ps); // mdstat bool save = (istep % inform::iwrite == 0); + if (save and useDebugArcFile) { + darray::copyout(g::q0, n, qx.data(), xpos); + darray::copyout(g::q0, n, qy.data(), ypos); + darray::copyout(g::q0, n, qz.data(), zpos); + waitFor(g::q0); + iarcdbg.writeCurrent(); + } if (save || (istep % BOUNDS_EVERY_X_STEPS) == 0) bounds(); if (save) { diff --git a/src/mdsave.cpp b/src/mdsave.cpp index 6d6a5bdff..79dc50a5d 100644 --- a/src/mdsave.cpp +++ b/src/mdsave.cpp @@ -178,89 +178,6 @@ void mdsaveAsync(int istep, time_prec dt) idle_dup = false; } -void mdDebugSaveSync() -{ - auto DOT3 = [](const double* a, const double* b) -> double { - return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; - }; - - static bool first = true; - static std::vector qx, qy, qz; - static std::string fname, title; - if (first) { - first = false; - qx.resize(n); - qy.resize(n); - qz.resize(n); - FstrView fstr = files::filename; - fname = fstr.trim(); - auto pos = fname.find_last_of('.'); - fname = fname.substr(0, pos) + ".dbg"; - FstrView ftitl = titles::title; - title = ftitl.trim(); - FILE* ftmp = fopen(fname.c_str(), "w"); - fclose(ftmp); - } - - FILE* fout = fopen(fname.c_str(), "a"); - - darray::copyout(g::q0, n, qx.data(), xpos); - darray::copyout(g::q0, n, qy.data(), ypos); - darray::copyout(g::q0, n, qz.data(), zpos); - waitFor(g::q0); - - bool bign = n > 999999; - if (bign) - fprintf(fout, "%8d %s\n", n, title.c_str()); - else - fprintf(fout, "%6d %s\n", n, title.c_str()); - - if (box_shape != BoxShape::UNBOUND) { - double ax[3] = {lvec1.x, lvec2.x, lvec3.x}; - double bx[3] = {lvec1.y, lvec2.y, lvec3.y}; - double cx[3] = {lvec1.z, lvec2.z, lvec3.z}; - - double xb = std::sqrt(DOT3(ax, ax)); - double yb = std::sqrt(DOT3(bx, bx)); - double zb = std::sqrt(DOT3(cx, cx)); - - double cos_a = DOT3(bx, cx) / (yb * zb); - double cos_b = DOT3(cx, ax) / (zb * xb); - double cos_c = DOT3(ax, bx) / (xb * yb); - - double al = 90.0; - double be = 90.0; - double ga = 90.0; - if (cos_a != 0.0) al = (180 / M_PI) * std::acos(cos_a); - if (cos_b != 0.0) be = (180 / M_PI) * std::acos(cos_b); - if (cos_c != 0.0) ga = (180 / M_PI) * std::acos(cos_c); - - if (bign) - fprintf(fout, " %14.6lf%14.6lf%14.6lf%14.6lf%14.6lf%14.6lf\n", xb, yb, zb, al, be, ga); - else - fprintf(fout, " %12.6lf%12.6lf%12.6lf%12.6lf%12.6lf%12.6lf\n", xb, yb, zb, al, be, ga); - } - - const char* fcord; - const char* fcoup; - if (bign) { - fcord = "%8d %c%c%c%14.6lf%14.6lf%14.6lf%8d"; - fcoup = "%8d"; - } else { - fcord = "%6d %c%c%c%12.6lf%12.6lf%12.6lf%6d"; - fcoup = "%6d"; - } - for (int i = 0; i < n; ++i) { - const char* nm = atomid::name[i]; - fprintf(fout, fcord, i + 1, nm[0], nm[1], nm[2], qx[i], qy[i], qz[i], atoms::type[i]); - for (int k = 0; k < couple::n12[i]; ++k) - fprintf(fout, fcoup, couple::i12[i][k]); - fprintf(fout, "%s", "\n"); - } - - fclose(fout); -} - void mdsaveSynchronize() { if (fut_dup_then_write.valid()) diff --git a/src/rwcrd.cpp b/src/rwcrd.cpp new file mode 100644 index 000000000..92f3fe707 --- /dev/null +++ b/src/rwcrd.cpp @@ -0,0 +1,540 @@ +#include "ff/rwcrd.h" +#include "ff/atom.h" +#include "ff/box.h" +#include "math/maxmin.h" +#include "tool/iofortstr.h" +#include "tool/ioprint.h" +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace tinker { +class CrdR +{ +protected: + std::ifstream ipt; + virtual void read() = 0; // read xyz from file and save them in atoms::x,y,z + +public: + static CrdR* create(std::string crdfile, CrdFormat crdformat); + virtual ~CrdR() {} + + int readCurrent() + { + this->read(); + xyzData(RcOp::INIT); + if (ipt.peek() == std::char_traits::eof()) + return 1; + else + return 0; + } +}; + +CrdReader::~CrdReader() +{ + delete m_impl; +} + +CrdReader::CrdReader(std::string crdfile, CrdFormat crdformat) + : m_impl(CrdR::create(crdfile, crdformat)) +{} + +int CrdReader::readCurrent() +{ + return m_impl->readCurrent(); +} +} + +namespace tinker { +class CrdW +{ +protected: + std::string m_name; + virtual void write(const double* qx, const double* qy, const double* qz) = 0; + +public: + static CrdW* create(std::string crdfile, CrdFormat crdformat); + virtual ~CrdW() {} + + CrdW(std::string crdfile) + : m_name(crdfile) + { + assert(crdfile != ""); + } + + int writeCurrent(const double* qx, const double* qy, const double* qz) + { + this->write(qx, qy, qz); + return 0; + } +}; + +CrdWriter::~CrdWriter() +{ + delete m_impl; +} + +CrdWriter::CrdWriter( + const double* xx, const double* yy, const double* zz, std::string crdfile, CrdFormat crdformat) + : m_impl(CrdW::create(crdfile, crdformat)) + , qx(xx) + , qy(yy) + , qz(zz) +{} + +int CrdWriter::writeCurrent() +{ + if (m_impl) + return m_impl->writeCurrent(qx, qy, qz); + else + return 0; +} +} + +//====================================================================// + +namespace tinker { +inline namespace v1 { +enum +{ + DCD_HEADER = 0, + DCD_TDELTA = 10, + DCD_USEBOX = 11, + DCD_CHARMM_VER = 20, + DCD_CTRL_LEN = 21, + + DCD_TITLE_NCHARS = 80, + + DCD_AX = 0, + DCD_COS_G = 1, + DCD_BX = 2, + DCD_COS_B = 3, + DCD_COS_A = 4, + DCD_CX = 5, + DCD_XTAL_LEN = 6, +}; +} + +class CrdRDcd : public CrdR +{ +private: + std::vector dcdx, dcdy, dcdz; + int dcdControl[DCD_CTRL_LEN] = {0}; + + static void readIntoBuffer(void* buffer, int nbyte, std::ifstream& ipt) + { + int size1, size2; + ipt.read((char*)&size1, sizeof(int)); + if (nbyte > 0) assert(nbyte == size1); + ipt.read((char*)buffer, size1); + ipt.read((char*)&size2, sizeof(int)); + } + + void read() override + { + if (dcdControl[DCD_USEBOX]) { + double dcdXtal[DCD_XTAL_LEN]; + readIntoBuffer(dcdXtal, sizeof(double) * DCD_XTAL_LEN, ipt); + double ax = dcdXtal[DCD_AX], bx = dcdXtal[DCD_BX], cx = dcdXtal[DCD_CX]; + double al = 90., be = 90., ga = 90.; + if (dcdXtal[DCD_COS_A] != 0.0) al = std::acos(dcdXtal[DCD_COS_A]) * (180 / M_PI); + if (dcdXtal[DCD_COS_B] != 0.0) be = std::acos(dcdXtal[DCD_COS_B]) * (180 / M_PI); + if (dcdXtal[DCD_COS_G] != 0.0) ga = std::acos(dcdXtal[DCD_COS_G]) * (180 / M_PI); + Box p; + boxLattice(p, box_shape, ax, bx, cx, al, be, ga); + boxSetCurrent(p); + } + + readIntoBuffer(dcdx.data(), sizeof(float) * n, ipt); + readIntoBuffer(dcdy.data(), sizeof(float) * n, ipt); + readIntoBuffer(dcdz.data(), sizeof(float) * n, ipt); + for (int i = 0; i < n; ++i) { + atoms::x[i] = dcdx[i]; + atoms::y[i] = dcdy[i]; + atoms::z[i] = dcdz[i]; + } + } + +public: + ~CrdRDcd() + { + ipt.close(); + } + + CrdRDcd(std::string crdfile) + : CrdR() + { + ipt.open(crdfile, std::ios::in | std::ios::binary); + + // read header info along with title and number of atoms + readIntoBuffer(dcdControl, sizeof(int) * DCD_CTRL_LEN, ipt); + + int dcdTitleRecordLen; + ipt.read((char*)&dcdTitleRecordLen, sizeof(int)); + std::vector titlebuf; + titlebuf.resize(dcdTitleRecordLen + sizeof(int)); + ipt.read(titlebuf.data(), dcdTitleRecordLen + sizeof(int)); + + int dcdNAtom; + readIntoBuffer(&dcdNAtom, sizeof(int), ipt); + assert(n == dcdNAtom); + dcdx.resize(n); + dcdy.resize(n); + dcdz.resize(n); + } +}; + +//====================================================================// + +class CrdRTxyz : public CrdR +{ +private: + bool m_usepbc; + + void read() override + { + std::string line; + std::getline(ipt, line); // n and title + if (m_usepbc) { + std::getline(ipt, line); // boxsize + double l1, l2, l3, a1, a2, a3; + std::sscanf(line.data(), "%lf%lf%lf%lf%lf%lf", &l1, &l2, &l3, &a1, &a2, &a3); + Box p; + boxLattice(p, box_shape, l1, l2, l3, a1, a2, a3); + boxSetCurrent(p); + } + + int index; + char name[32]; + double xr, yr, zr; + for (int ir = 0; ir < n; ++ir) { + std::getline(ipt, line); + std::sscanf(line.data(), "%d%s%lf%lf%lf", &index, name, &xr, &yr, &zr); + index -= 1; + atoms::x[index] = xr; + atoms::y[index] = yr; + atoms::z[index] = zr; + } + } + +public: + ~CrdRTxyz() + { + ipt.close(); + } + + CrdRTxyz(std::string crdfile, bool usepbc) + : CrdR() + , m_usepbc(usepbc) + { + ipt.open(crdfile, std::ios::in); + } +}; + +//====================================================================// + +CrdR* CrdR::create(std::string crdfile, CrdFormat crdformat) +{ + auto fmt = crdformat; + if (fmt == CrdFormat::NONE) { + // get file format type by inspection of first character + char a1; + std::ifstream fipt(crdfile); + fipt >> a1; + + if (a1 == ' ') + fmt = CrdFormat::TXYZ1; + else if ('0' <= a1 and a1 <= '9') + fmt = CrdFormat::TXYZ1; + else if (a1 == 'T') + fmt = CrdFormat::DCD; + + if (fmt == CrdFormat::TXYZ1) { + std::string line; + // rewind + fipt.clear(); + fipt.seekg(0); + std::getline(fipt, line); // n and title + std::getline(fipt, line); // box size or first atom + // 18.643000 18.643000 18.643000 90.000000 90.000000 90.000000 + // 1 O 8.733783 7.084710 -0.688468 1 2 3 + double l1, l2, l3, a1, a2, a3; + int matched = std::sscanf(line.data(), "%lf%lf%lf%lf%lf%lf", &l1, &l2, &l3, &a1, &a2, &a3); + if (matched == 6) fmt = CrdFormat::TXYZ2_PBC; + } + } + + CrdR* p = nullptr; + switch (fmt) { + case CrdFormat::DCD: + p = new CrdRDcd(crdfile); + break; + case CrdFormat::TXYZ2_PBC: + p = new CrdRTxyz(crdfile, true); + break; + default /* CrdFormat::TXYZ1 */: + p = new CrdRTxyz(crdfile, false); + break; + } + return p; +} +} + +//====================================================================// + +namespace tinker { +static auto d3__ = [](const double* a, const double* b) -> double { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +}; + +class CrdWDcd : public CrdW +{ +private: + std::vector dcdx, dcdy, dcdz; + int dcdControl[DCD_CTRL_LEN] = {0}; + + static void writeBufferToExtFile(const void* buffer, int nbyte, FILE* fout) + { + int pad = nbyte; + std::fwrite(&pad, sizeof(int), 1, fout); + std::fwrite(buffer, nbyte, 1, fout); + std::fwrite(&pad, sizeof(int), 1, fout); + } + + void write(const double* qx, const double* qy, const double* qz) override + { + FILE* fout = fopen(m_name.c_str(), "ab"); + if (dcdControl[DCD_USEBOX]) { + double ax[3] = {lvec1.x, lvec2.x, lvec3.x}; + double bx[3] = {lvec1.y, lvec2.y, lvec3.y}; + double cx[3] = {lvec1.z, lvec2.z, lvec3.z}; + + double xb = std::sqrt(d3__(ax, ax)); + double yb = std::sqrt(d3__(bx, bx)); + double zb = std::sqrt(d3__(cx, cx)); + + double cos_a = d3__(bx, cx) / (yb * zb); + double cos_b = d3__(cx, ax) / (zb * xb); + double cos_c = d3__(ax, bx) / (xb * yb); + + double dcdXtal[DCD_XTAL_LEN]; + dcdXtal[DCD_AX] = xb; + dcdXtal[DCD_COS_G] = cos_c; + dcdXtal[DCD_BX] = yb; + dcdXtal[DCD_COS_B] = cos_b; + dcdXtal[DCD_COS_A] = cos_a; + dcdXtal[DCD_CX] = zb; + + writeBufferToExtFile(&dcdXtal[0], sizeof(double) * DCD_XTAL_LEN, fout); + } + for (int i = 0; i < n; ++i) { + dcdx[i] = qx[i]; + dcdy[i] = qy[i]; + dcdz[i] = qz[i]; + } + writeBufferToExtFile(dcdx.data(), sizeof(float) * n, fout); + writeBufferToExtFile(dcdy.data(), sizeof(float) * n, fout); + writeBufferToExtFile(dcdz.data(), sizeof(float) * n, fout); + fclose(fout); + } + +public: + CrdWDcd(std::string crdfile) + : CrdW(crdfile) + { + for (int i = 0; i < DCD_CTRL_LEN; ++i) + dcdControl[i] = 0; + const char* cord = "CORD"; + std::memcpy(&dcdControl[DCD_HEADER], cord, 4); + if (bound::use_bounds) dcdControl[DCD_USEBOX] = 1; + dcdControl[DCD_CHARMM_VER] = 24; + + int dcdTitle[21] = {0}; + char dcd4Space[4] = {' ', ' ', ' ', ' '}; + static_assert(4 == sizeof(int), ""); + int dcd4SpaceInt; + std::memcpy(&dcd4SpaceInt, dcd4Space, sizeof(int)); + for (int i = 1; i < 21; ++i) + dcdTitle[i] = dcd4SpaceInt; + if (titles::ltitle == 0) { + dcdTitle[0] = 0; + } else { + dcdTitle[0] = 1; + std::memcpy(&dcdTitle[1], titles::title, std::min(80, titles::ltitle)); + } + + int dcdNAtom = n; + dcdx.resize(dcdNAtom); + dcdy.resize(dcdNAtom); + dcdz.resize(dcdNAtom); + + FILE* fout = fopen(m_name.c_str(), "ab"); + writeBufferToExtFile(&dcdControl[0], sizeof(int) * DCD_CTRL_LEN, fout); + writeBufferToExtFile(&dcdTitle[0], sizeof(int) * 21, fout); + writeBufferToExtFile(&dcdNAtom, sizeof(int), fout); + fclose(fout); + } +}; + +//====================================================================// + +class CrdWTxyz : public CrdW +{ +private: + std::string title_line; + int ilen, digc; + bool m_usepbc; + + static const char* formatTitleLine(int ilen, bool hasTitle) + { + static const char* f[2][3] = {{"%6d", "%7d", "%8d"}, {"%6d %s", "%7d %s", "%8d %s"}}; + int i = (hasTitle ? 1 : 0); + assert(6 <= ilen and ilen <= 8); + return f[i][ilen - 6]; + } + + static const char* formatInt(int ilen) + { + static const char* f[3] = {"%6d", "%7d", "%8d"}; + assert(6 <= ilen and ilen <= 8); + return f[ilen - 6]; + } + + static const char* formatFlt(int crdsiz, int digc) + { + static const char* f[3][3] = { + {"%12.6lf", "%14.8lf", "%16.10lf"}, // (6,6) (6,8) (6,10) + {"%13.6lf", "%15.8lf", "%17.10lf"}, // (7,6) (7,8) (7,10) + {"%14.6lf", "%16.8lf", "%18.10lf"}, // (8,6) (8,8) (8,10) + }; + assert(crdsiz == 6 or crdsiz == 7 or crdsiz == 8); + assert(digc == 6 or digc == 8 or digc == 10); + return f[crdsiz - 6][digc / 2 - 3]; + } + + void write(const double* qx, const double* qy, const double* qz) override + { + int crdsiz = 6; + double crdmin = 0., crdmax = 0.; + for (int i = 0; i < n; ++i) { + double xi = qx[i], yi = qy[i], zi = qz[i]; + crdmin = minOf(crdmin, xi, yi, zi); + crdmax = maxOf(crdmax, xi, yi, zi); + } + if (crdmin <= -1000.) crdsiz = 7; + if (crdmax >= 10000.) crdsiz = 7; + if (crdmin <= -10000.) crdsiz = 8; + if (crdmax >= 100000.) crdsiz = 8; + const char* fflt = formatFlt(crdsiz, digc); + const char* fint = formatInt(ilen); + + FILE* fout = fopen(m_name.c_str(), "a"); + std::fprintf(fout, "%s\n", title_line.c_str()); + if (m_usepbc) { + double ax[3] = {lvec1.x, lvec2.x, lvec3.x}; + double bx[3] = {lvec1.y, lvec2.y, lvec3.y}; + double cx[3] = {lvec1.z, lvec2.z, lvec3.z}; + + double xb = std::sqrt(d3__(ax, ax)); + double yb = std::sqrt(d3__(bx, bx)); + double zb = std::sqrt(d3__(cx, cx)); + + double cos_a = d3__(bx, cx) / (yb * zb); + double cos_b = d3__(cx, ax) / (zb * xb); + double cos_c = d3__(ax, bx) / (xb * yb); + + double al = 90.0; + double be = 90.0; + double ga = 90.0; + if (cos_a != 0.0) al = (180 / M_PI) * std::acos(cos_a); + if (cos_b != 0.0) be = (180 / M_PI) * std::acos(cos_b); + if (cos_c != 0.0) ga = (180 / M_PI) * std::acos(cos_c); + + std::string fmt_box = " "; + fmt_box = fmt_box + fflt + fflt + fflt + fflt + fflt + fflt + "\n"; + std::fprintf(fout, fmt_box.c_str(), xb, yb, zb, al, be, ga); + } + + std::string fmt_atom = fint; + fmt_atom = fmt_atom + " %c%c%c" + fflt + fflt + fflt + "%6d"; + for (int i = 0; i < n; ++i) { + const char* nm = atomid::name[i]; + std::fprintf(fout, fmt_atom.c_str(), i + 1, nm[0], nm[1], nm[2], qx[i], qy[i], qz[i], + atoms::type[i]); + for (int k = 0; k < couple::n12[i]; ++k) + std::fprintf(fout, fint, couple::i12[i][k]); + std::fprintf(fout, "%s", "\n"); + } + fclose(fout); + } + +public: + CrdWTxyz(std::string crdfile, bool usepbc) + : CrdW(crdfile) + , m_usepbc(usepbc) + { + ilen = 6; + if (n >= 100000) ilen = 7; + if (n >= 1000000) ilen = 8; + const char* f; + + if (titles::ltitle > 0) { + f = formatTitleLine(ilen, true); + FstrView ftitl = titles::title; + title_line = format(f, n, ftitl.trim()); + } else { + f = formatTitleLine(ilen, true); + title_line = format(f, n); + } + + if (inform::digits <= 6) + digc = 6; + else if (inform::digits <= 8) + digc = 8; + else + digc = 10; + } +}; + +//====================================================================// + +CrdW* CrdW::create(std::string crdfile, CrdFormat crdformat) +{ + if (crdfile == "") return nullptr; + + auto fmt = crdformat; + if (fmt == CrdFormat::NONE) { + if (output::dcdsave) { + fmt = CrdFormat::DCD; + } else if (output::arcsave) { + if (bound::use_bounds) + fmt = CrdFormat::TXYZ2_PBC; + else + fmt = CrdFormat::TXYZ1; + } + } + + CrdW* p = nullptr; + switch (fmt) { + case CrdFormat::DCD: + p = new CrdWDcd(crdfile); + break; + case CrdFormat::TXYZ2_PBC: + p = new CrdWTxyz(crdfile, true); + break; + default /* CrdFormat::TXYZ1 */: + p = new CrdWTxyz(crdfile, false); + break; + } + return p; +} +} diff --git a/src/xanalyze.cpp b/src/xanalyze.cpp index 6ffb98a7e..66055b6bc 100644 --- a/src/xanalyze.cpp +++ b/src/xanalyze.cpp @@ -11,6 +11,7 @@ #include "ff/hippomod.h" #include "ff/nblist.h" #include "ff/potent.h" +#include "ff/rwcrd.h" #include "md/osrw.h" #include "tool/argkey.h" #include "tool/iofortstr.h" @@ -426,9 +427,7 @@ void xAnalyze(int, char**) bool exist = false; std::string opt; nextarg(string, exist); - if (exist) { - ioReadString(opt, string); - } + if (exist) ioReadString(opt, string); std::string prompt = R"( The Tinker Energy Analysis Utility Can : @@ -456,24 +455,19 @@ void xAnalyze(int, char**) auto out = stdout; FstrView fsw = files::filename; std::string fname = fsw.trim(); - std::ifstream ipt; - readFrameOpen(fname, ipt); - bool done = false; int nframe_processed = 0; + int done = 0; + auto ipt = CrdReader(fname); do { - readFrameCopyinToXyz(ipt, done); + done = ipt.readCurrent(); nblistRefresh(); nframe_processed++; if (nframe_processed > 1) print(out, "\n Analysis for Archive Structure :%16d\n", nframe_processed); - if (opt.find("E") != failed) - xAnalyzeE(); - if (opt.find("M") != failed) - xAnalyzeM(); - if (opt.find("V") != failed) - xAnalyzeV(); + if (opt.find("E") != failed) xAnalyzeE(); + if (opt.find("M") != failed) xAnalyzeM(); + if (opt.find("V") != failed) xAnalyzeV(); } while (not done); - readFrameClose(ipt); finish(); tinker_f_final(); diff --git a/src/xbar.cpp b/src/xbar.cpp index af8903a8f..8e250eb8d 100644 --- a/src/xbar.cpp +++ b/src/xbar.cpp @@ -1,5 +1,6 @@ #include "ff/energy.h" #include "ff/nblist.h" +#include "ff/rwcrd.h" #include "math/random.h" #include "tool/argkey.h" #include "tool/iofortstr.h" @@ -112,7 +113,7 @@ static void xBarMake() //====================================================================// - bool done; + int done; const char* log_fmt = " Current Potential %lf\n"; const char* process_fmt = " Completed%8d Coordinate Frames\n"; std::memcpy(string, filea, lenga); @@ -152,7 +153,6 @@ static void xBarMake() tinker_f_open(&iarc, str, "old"); else if (output::binary) tinker_f_open(&iarc, str, "unformatted", "old"); - std::ifstream a_arc; if (ua0.size() == 0) { // reset trajectory A using the parameters for state 0 @@ -166,10 +166,10 @@ static void xBarMake() // find potential energies for trajectory A in state 0 initialize(); - readFrameOpen(str, a_arc); - done = false; + done = 0; + auto a_arc = CrdReader(str); do { - readFrameCopyinToXyz(a_arc, done); + done = a_arc.readCurrent(); nblistRefresh(); energy(calc::energy); ua0.push_back(esum); @@ -179,7 +179,6 @@ static void xBarMake() std::fflush(out); } } while (not done); - readFrameClose(a_arc); finish(); } else { int ii = ua0.size(); @@ -208,10 +207,10 @@ static void xBarMake() " Frame State 0 State 1 Delta\n" "\n"); initialize(); - readFrameOpen(str, a_arc); - done = false; + done = 0; + auto a_ar2 = CrdReader(str); do { - readFrameCopyinToXyz(a_arc, done); + done = a_ar2.readCurrent(); nblistRefresh(); energy(calc::energy); double vol = boxVolume(); @@ -222,7 +221,6 @@ static void xBarMake() print(out, "%11d %16.4lf%16.4lf%16.4lf\n", i + 1, ua0[i], ua1[i], ua1[i] - ua0[i]); } } while (not done); - readFrameClose(a_arc); finish(); // save potential energies and volumes for trajectory A @@ -286,10 +284,10 @@ static void xBarMake() // find potential energies for trajectory B in state 1 initialize(); - readFrameOpen(str, b_arc); - done = false; + done = 0; + auto b_arc = CrdReader(str); do { - readFrameCopyinToXyz(b_arc, done); + done = b_arc.readCurrent(); nblistRefresh(); energy(calc::energy); ub1.push_back(esum); @@ -299,7 +297,6 @@ static void xBarMake() std::fflush(out); } } while (not done); - readFrameClose(b_arc); finish(); } else { int ii = ub1.size(); @@ -328,10 +325,10 @@ static void xBarMake() " Frame State 0 State 1 Delta\n" "\n"); initialize(); - readFrameOpen(str, b_arc); - done = false; + done = 0; + auto b_ar2 = CrdReader(str); do { - readFrameCopyinToXyz(b_arc, done); + done = b_ar2.readCurrent(); nblistRefresh(); energy(calc::energy); double vol = boxVolume(); @@ -342,7 +339,6 @@ static void xBarMake() print(out, "%11d %16.4lf%16.4lf%16.4lf\n", i + 1, ub0[i], ub1[i], ub0[i] - ub1[i]); } } while (not done); - readFrameClose(b_arc); finish(); // save potential energies and volumes for trajectory B diff --git a/src/xtestgrad.cpp b/src/xtestgrad.cpp index 3058d51ba..af5cd6f21 100644 --- a/src/xtestgrad.cpp +++ b/src/xtestgrad.cpp @@ -1,5 +1,6 @@ #include "ff/energy.h" #include "ff/nblist.h" +#include "ff/rwcrd.h" #include "tool/iofortstr.h" #include "tool/ioprint.h" #include @@ -49,12 +50,11 @@ void xTestgrad(int, char**) FstrView fsw = files::filename; std::string fname = fsw.trim(); - std::ifstream ipt; - readFrameOpen(fname, ipt); - bool done = false; int nframe_processed = 0; + int done = 0; + auto ipt = CrdReader(fname); do { - readFrameCopyinToXyz(ipt, done); + done = ipt.readCurrent(); nblistRefresh(); nframe_processed++; if (nframe_processed > 1) @@ -112,7 +112,6 @@ void xTestgrad(int, char**) print(out, fmt3, "Total Gradient Norm Value", totnorm, len3, digits); print(out, fmt3, "RMS Gradient over All Atoms", rms, len3, digits); } while (not done); - readFrameClose(ipt); finish(); tinker_f_final(); diff --git a/test/bounds.cpp b/test/bounds.cpp index 0a19e3df5..c79e6e379 100644 --- a/test/bounds.cpp +++ b/test/bounds.cpp @@ -1,6 +1,6 @@ #include "ff/box.h" +#include "ff/rwcrd.h" #include "tool/iofortstr.h" -#include #include #include "test.h" @@ -39,9 +39,8 @@ TEST_CASE("Bounds", "[ff][box]") FstrView fsw = files::filename; std::string fname = fsw.trim(); - std::ifstream ipt(fname); - bool done = false; - readFrameCopyinToXyz(ipt, done); + auto ipt = CrdReader(fname); + ipt.readCurrent(); bounds(); double eps = 1.0e-5;