add dcd support

Author: zhi-wang

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                          //

include/ff/rwcrd.h

@@ -0,0 +1,38 @@
+#pragma once
+#include <string>
+
+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();
+};
+}

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

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" \

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<float> 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<char> 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;
-}
 }

src/cmakesrc.txt

@@ -72,6 +72,7 @@ pme.cpp
 potent.cpp
 rattle.cpp
 rcman.cpp
+rwcrd.cpp
 spatial.cpp
 switch.cpp
 tinkersuppl.cpp

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 <cassert>
+#include "tool/tinkersuppl.h"
+#include <tinker/detail/files.hh>
 #include <tinker/detail/inform.hh>
 #include <tinker/detail/mdstuf.hh>
 #include <tinker/detail/units.hh>
 
+#include <cassert>
+
 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<pos_prec> 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) {