add dcd support in analyze bar testgrad

Author: zhi-wang

include/ff/atom.h

@@ -1,7 +1,7 @@
 #pragma once
 #include "ff/precision.h"
 #include "tool/rcman.h"
-#include <istream>
+#include <fstream>
 
 namespace tinker {
 /// \addtogroup ff
@@ -29,7 +29,9 @@ void copyPosToXyz(bool refreshNBList);
 ///    - Tinker uses centers of mass.
 void bounds();
 
-void readFrameCopyinToXyz(std::istream& input, bool& done);
+void readFrameOpen(const std::string& filename, std::ifstream& input);
+void readFrameCopyinToXyz(std::ifstream& input, bool& done);
+void readFrameClose(std::ifstream& input);
 
 //====================================================================//
 //                                                                    //

include/tool/tinkersuppl.h

@@ -10,6 +10,7 @@ void tinker_f_rewind(int* unit);
 void tinker_f_close(int* unit);
 
 void tinker_f_open(int* unit, std::string file, std::string status);
+void tinker_f_open(int* unit, std::string file, std::string form, std::string status);
 
 // memory
 int tinker_f_allocated(void* p);

src/atom.cpp

@@ -6,12 +6,13 @@
 #include "tool/error.h"
 #include "tool/externfunc.h"
 #include "tool/gpucard.h"
-#include <cassert>
 #include <tinker/detail/atomid.hh>
 #include <tinker/detail/atoms.hh>
 #include <tinker/detail/bound.hh>
 #include <tinker/detail/usage.hh>
 
+#include <cassert>
+
 namespace tinker {
 void nData(RcOp op)
 {
@@ -154,49 +155,166 @@ void bounds()
    copyPosToXyz();
 }
 
-void readFrameCopyinToXyz(std::istream& ipt, bool& done)
+inline namespace v1 {
+enum
 {
-   if (done)
-      return;
+   DCD_HEADER = 0,
+   DCD_TDELTA = 10,
+   DCD_USEBOX = 11,
+   DCD_CTRL_LEN = 21,
 
-   if (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;
-      }
+   DCD_TITLE_NCHARS = 80,
 
-      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;
-      }
+   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;
+   }
 
-      xyzData(RcOp::INIT);
+   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 (ipt.peek() == EOF)
-      done = true;
+   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/f/suppl.f

@@ -27,15 +27,15 @@ subroutine suppl_close (unit)
 
 
 #ifdef TINKER_SUPPL_DECL
-      void suppl_open_(int* unit, const char* file, const char* status,
-         tinker_fchar_len_t file_cap, tinker_fchar_len_t status_cap);
+      void suppl_open_(int* unit, const char* file, const char* form, const char* status,
+         tinker_fchar_len_t file_cap, tinker_fchar_len_t form_cap, tinker_fchar_len_t status_cap);
 #endif
 #ifdef TINKER_SUPPL_IMPL
-      subroutine suppl_open (unit,file,status)
+      subroutine suppl_open (unit,file,form,status)
       implicit none
       integer unit
-      character*(*) file,status
-      open (unit=unit,file=file,status=status)
+      character*(*) file,form,status
+      open (unit=unit,file=file,form=form,status=status)
       return
       end
 #endif

src/tinkersuppl.cpp

@@ -29,7 +29,15 @@ void tinker_f_close(int* unit)
 
 void tinker_f_open(int* unit, std::string file, std::string status)
 {
-   suppl_open_(unit, file.c_str(), status.c_str(), file.length(), status.length());
+   std::string form = "formatted";
+   suppl_open_(unit, file.c_str(), form.c_str(), status.c_str(), file.length(), form.length(),
+      status.length());
+}
+
+void tinker_f_open(int* unit, std::string file, std::string form, std::string status)
+{
+   suppl_open_(unit, file.c_str(), form.c_str(), status.c_str(), file.length(), form.length(),
+      status.length());
 }
 
 int tinker_f_allocated(void* p)

src/xanalyze.cpp

@@ -33,8 +33,6 @@
 
 #include "tinker9.h"
 
-#include <fstream>
-
 namespace tinker {
 static void xAnalyzeE()
 {
@@ -419,7 +417,8 @@ namespace tinker {
 void xAnalyze(int, char**)
 {
    initial();
-   tinker_f_getxyz();
+   int ixyz;
+   tinker_f_getcart(&ixyz);
    tinker_f_mechanic();
    mechanic2();
 
@@ -457,7 +456,8 @@ void xAnalyze(int, char**)
    auto out = stdout;
    FstrView fsw = files::filename;
    std::string fname = fsw.trim();
-   std::ifstream ipt(fname);
+   std::ifstream ipt;
+   readFrameOpen(fname, ipt);
    bool done = false;
    int nframe_processed = 0;
    do {
@@ -473,6 +473,7 @@ void xAnalyze(int, char**)
       if (opt.find("V") != failed)
          xAnalyzeV();
    } while (not done);
+   readFrameClose(ipt);
 
    finish();
    tinker_f_final();