Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

lmplugin.C

Go to the documentation of this file.
00001 /*
00002  * University of Illinois Open Source License
00003  * Copyright 2018-2018 Luthey-Schulten Group,
00004  * All rights reserved.
00005  *
00006  * Developed by: Luthey-Schulten Group
00007  *                           University of Illinois at Urbana-Champaign
00008  *                           http://www.scs.uiuc.edu/~schulten
00009  *
00010  * Permission is hereby granted, free of charge, to any person obtaining a copy of
00011  * this software and associated documentation files (the Software), to deal with
00012  * the Software without restriction, including without limitation the rights to
00013  * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
00014  * of the Software, and to permit persons to whom the Software is furnished to
00015  * do so, subject to the following conditions:
00016  *
00017  * - Redistributions of source code must retain the above copyright notice,
00018  * this list of conditions and the following disclaimers.
00019  *
00020  * - Redistributions in binary form must reproduce the above copyright notice,
00021  * this list of conditions and the following disclaimers in the documentation
00022  * and/or other materials provided with the distribution.
00023  *
00024  * - Neither the names of the Luthey-Schulten Group, University of Illinois at
00025  * Urbana-Champaign, nor the names of its contributors may be used to endorse or
00026  * promote products derived from this Software without specific prior written
00027  * permission.
00028  *
00029  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
00030  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
00031  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
00032  * THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
00033  * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
00034  * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
00035  * OTHER DEALINGS WITH THE SOFTWARE.
00036  *
00037  * Author(s): Tyler M. Earnest
00038  */
00039 
00040 #include <cmath>
00041 #include <cstdarg>
00042 #include <cstdio>
00043 #include <cstring>
00044 #include <ctime>
00045 
00046 #include <algorithm>
00047 #include <exception>
00048 #include <numeric>
00049 #include <set>
00050 #include <sstream>
00051 #include <string>
00052 #include <vector>
00053 
00054 #include <hdf5.h>
00055 #include <hdf5_hl.h>
00056 
00057 #include "molfile_plugin.h"
00058 #include "vmdconio.h"
00059 
00060 namespace {
00061 
00062 #define MAX_VMD_NAME 16
00063 #define RDME_RADIUS_SCALE 0.1
00064 #define TO_ANGSTROM 1e10
00065 #define MAX_LMPEXC_MSG 256
00066 
00067 typedef unsigned char particle_int_t;
00068 typedef unsigned char site_int_t;
00069 
00070 class LMPException : public std::exception {
00071 protected:
00072     char msg[MAX_LMPEXC_MSG];
00073 public:
00074     void vmdcon_report() { vmdcon_printf(VMDCON_ERROR, "LMPlugin) Error: %s\n", msg); }
00075 };
00076 
00077 class Exception : public LMPException {
00078 public:
00079     Exception(const char *fmt, ...)
00080     {
00081         va_list ap;
00082         va_start(ap, fmt);
00083         vsnprintf(msg, sizeof(msg), fmt, ap);
00084         va_end(ap);
00085     }
00086 };
00087 
00088 class H5Exception : public LMPException {
00089     static herr_t
00090     walk_cb(unsigned int n,
00091             const H5E_error2_t *err,
00092             void *data)
00093     {
00094         if (n == 0) {
00095             snprintf(static_cast<char*>(data), MAX_LMPEXC_MSG, "HDF5: %s: %s",
00096                     err->func_name, err->desc);
00097         }
00098         return 0;
00099     }
00100 
00101 public:
00102     H5Exception() { H5Ewalk2(H5E_DEFAULT, H5E_WALK_DOWNWARD, H5Exception::walk_cb, msg); }
00103 };
00104 
00105 #define H5_CALL_ASSIGN(val, call)                                                                  \
00106     do {                                                                                           \
00107         val = (call);                                                                              \
00108         if (val < 0) {                                                                             \
00109             throw H5Exception();                                                                   \
00110         }                                                                                          \
00111     } while (0)
00112 
00113 #define H5_CALL(call)                                                                              \
00114     do {                                                                                           \
00115         hid_t val;                                                                                 \
00116         H5_CALL_ASSIGN(val, (call));                                                               \
00117     } while (0)
00118 
00119 
00120 struct Timer {
00121     // using low resolution time() to avoid portability issues
00122     time_t start, end;
00123 
00124     Timer() { reset(); }
00125     void reset() { time(&start); }
00126 
00127     int
00128     elapsed()
00129     {
00130         time(&end);
00131         return end-start;
00132     }
00133 };
00134 
00135 namespace H5File_specialized {
00136 
00137 template <typename T> hid_t get_type_id();
00138 template <> inline hid_t get_type_id<char>() { return H5T_STD_I8LE; }
00139 template <> inline hid_t get_type_id<unsigned char>() { return H5T_STD_U8LE; }
00140 template <> inline hid_t get_type_id<int>() { return H5T_STD_I32LE; }
00141 template <> inline hid_t get_type_id<unsigned int>() { return H5T_STD_U32LE; }
00142 template <> inline hid_t get_type_id<float>() { return H5T_IEEE_F32LE; }
00143 template <> inline hid_t get_type_id<double>() { return H5T_IEEE_F64LE; }
00144 
00145 template <typename T>
00146 void
00147 get_attr(hid_t loc_id,
00148          std::string obj_name,
00149          std::string attr_name,
00150          T &data)
00151 {
00152     H5_CALL(H5LTget_attribute(loc_id, obj_name.c_str(), attr_name.c_str(), get_type_id<T>(),
00153                 &data));
00154 }
00155 
00156 template <>
00157 void
00158 get_attr<std::string>(hid_t loc_id,
00159                       std::string obj_name,
00160                       std::string attr_name,
00161                       std::string &str)
00162 {
00163     hid_t attr, type;
00164     hsize_t size;
00165 
00166     H5_CALL_ASSIGN(attr, H5Aopen_by_name(loc_id, obj_name.c_str(), attr_name.c_str(),
00167                 H5P_DEFAULT, H5P_DEFAULT));
00168     H5_CALL_ASSIGN(type, H5Aget_type(attr));
00169 
00170     if (H5Tget_class(type) == H5T_STRING) {
00171         hid_t space;
00172         H5_CALL_ASSIGN(space, H5Aget_space(attr));
00173         H5_CALL_ASSIGN(size, H5Aget_storage_size(attr));
00174 
00175         hid_t memtype;
00176         H5_CALL_ASSIGN(memtype, H5Tcopy(H5T_C_S1));
00177         H5_CALL(H5Tset_size(memtype, size));
00178 
00179         char *data = new char[size + 1];
00180         H5_CALL(H5Aread(attr, memtype, data));
00181         data[size] = 0;
00182         str = std::string(data);
00183         delete [] data;
00184 
00185         H5_CALL(H5Sclose(space));
00186         H5_CALL(H5Tclose(memtype));
00187     }
00188 
00189     H5_CALL(H5Tclose(type));
00190     H5_CALL(H5Aclose(attr));
00191 }
00192 
00193 template <typename T>
00194 bool
00195 str_conv1(std::string str,
00196           T &val)
00197 {
00198     std::istringstream ss(str);
00199     return ss >> val;
00200 }
00201 
00202 template <typename T>
00203 bool
00204 str_conv_last(std::string str,
00205               std::vector<T> &vals)
00206 {
00207     return str.size() == 0;
00208 }
00209 
00210 // empty strings are valid
00211 template <>
00212 bool
00213 str_conv1(std::string str,
00214           std::string &val)
00215 {
00216     std::istringstream ss(str);
00217     ss >> val;
00218     return true;
00219 }
00220 
00221 template <>
00222 bool
00223 str_conv_last(std::string str,
00224               std::vector<std::string> &vals)
00225 {
00226     if (str.size() > 0 && *(str.end()-1) == ',') {
00227         vals.push_back("");
00228     }
00229     return true;
00230 }
00231 } // end namespace H5File_specialized
00232 
00233 class H5File {
00234     hid_t h5file;
00235     std::string filename;
00236 
00237     static herr_t
00238     get_group_list__proc(hid_t g_id,
00239                          const char *name,
00240                          const H5L_info_t *info,
00241                          void *op_data)
00242     {
00243         std::vector<std::string> *vec = static_cast<std::vector<std::string> *>(op_data);
00244         vec->push_back(name);
00245         return 0;
00246     }
00247 
00248 public:
00249 
00250     H5File () : h5file(-1) {}
00251     ~H5File () { if (h5file>=0) { close(); } }
00252 
00253     void
00254     open(std::string filename_)
00255     {
00256         filename = filename_;
00257         H5_CALL_ASSIGN(h5file, H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT));
00258     }
00259 
00260     void
00261     open_rdwr(std::string filename_)
00262     {
00263         filename = filename_;
00264         H5_CALL_ASSIGN(h5file, H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT));
00265     }
00266 
00267     void
00268     close()
00269     {
00270         H5_CALL(H5Fclose(h5file));
00271         h5file = -1;
00272     }
00273 
00274     void
00275     get_userblock(std::vector<char>& data)
00276     {
00277         hid_t userblock;
00278         hsize_t userblockSize;
00279         H5_CALL_ASSIGN(userblock, H5Fget_create_plist(h5file));
00280         if (H5Pget_userblock(userblock, &userblockSize) < 0) {
00281             data.resize(0);
00282         }
00283         else {
00284             data.resize(userblockSize);
00285             close();
00286             FILE *fp;
00287             fp = fopen(filename.c_str(), "rb");
00288             fread(&data[0], userblockSize, 1, fp);
00289             fclose(fp);
00290             open(filename);
00291         }
00292     }
00293 
00294     template <typename T>
00295     void
00296     get_attr(std::string obj_name,
00297              std::string attr_name,
00298              T &data)
00299     {
00300         H5File_specialized::get_attr(h5file, obj_name, attr_name, data);
00301     }
00302 
00303     template <typename T>
00304     void
00305     get_csv_attr(std::string obj_name,
00306                  std::string attr_name,
00307                  std::vector<T> &data)
00308     {
00309         std::string str;
00310         H5File_specialized::get_attr(h5file, obj_name.c_str(), attr_name.c_str(), str);
00311 
00312         std::istringstream ss;
00313         ss.str(str);
00314         data.clear();
00315 
00316         std::string sval;
00317 
00318         while (std::getline(ss, sval, ',')) {
00319             T val;
00320             if (! H5File_specialized::str_conv1(sval, val) ) {
00321                 throw Exception("String conversion failed for attribute %s on %s: Bad token: \"%s\"",
00322                         attr_name.c_str(), obj_name.c_str(), sval.c_str());
00323             }
00324             data.push_back(val);
00325         }
00326 
00327         if ( !H5File_specialized::str_conv_last(sval, data) ) {
00328             throw Exception("String conversion failed for attribute %s on %s: Bad token: \"%s\"",
00329                     attr_name.c_str(), obj_name.c_str(), sval.c_str());
00330         }
00331     }
00332 
00333     template <typename T>
00334     void
00335     get_dataset(std::string path,
00336                 std::vector<T> &data,
00337                 std::vector<hsize_t> &shape,
00338                 hsize_t* nelemptr=NULL)
00339     {
00340         hid_t datasetHandle, dataspaceHandle, typeHandle, h5type;
00341         hsize_t dims[H5S_MAX_RANK];
00342         hsize_t ndims;
00343 
00344         H5_CALL_ASSIGN(datasetHandle, H5Dopen2(h5file, path.c_str(), H5P_DEFAULT));
00345         H5_CALL_ASSIGN(dataspaceHandle, H5Dget_space(datasetHandle));
00346         H5_CALL(H5Sget_simple_extent_dims(dataspaceHandle, dims, NULL));
00347         H5_CALL_ASSIGN(ndims, H5Sget_simple_extent_ndims(dataspaceHandle));
00348         H5_CALL_ASSIGN(typeHandle, H5Dget_type(datasetHandle));
00349         H5_CALL_ASSIGN(h5type, H5Tget_native_type(typeHandle, H5T_DIR_ASCEND));
00350 
00351         if (H5Tequal(H5File_specialized::get_type_id<T>(), h5type) <= 0) {
00352             throw Exception("Wrong HDF5 data type for ``%s''", path.c_str());
00353         }
00354 
00355         shape.resize(ndims);
00356         std::copy(dims, dims + ndims, shape.begin());
00357 
00358         hsize_t nelem = std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<hsize_t>());
00359 
00360         if (nelemptr) {
00361             *nelemptr = nelem;
00362         }
00363 
00364         if (nelem != data.size()) {
00365             data.resize(nelem);
00366         }
00367 
00368         H5_CALL(H5Dread(datasetHandle, h5type, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(data[0])));
00369 
00370         H5_CALL(H5Tclose(typeHandle));
00371         H5_CALL(H5Sclose(dataspaceHandle));
00372         H5_CALL(H5Dclose(datasetHandle));
00373     }
00374 
00375     bool
00376     group_exists(std::string path)
00377     {
00378         H5E_auto2_t funcOld;
00379         void *dataOld;
00380         bool exists;
00381 
00382         H5Eget_auto2(H5E_DEFAULT, &funcOld, &dataOld);
00383         H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
00384 
00385         hid_t handle = H5Gopen2(h5file, path.c_str(), H5P_DEFAULT);
00386 
00387         exists = handle >= 0;
00388         if (exists) {
00389             H5Gclose(handle);
00390         }
00391         H5Eset_auto2(H5E_DEFAULT, funcOld, dataOld);
00392         return exists;
00393     }
00394 
00395     bool
00396     dataset_exists(std::string path)
00397     {
00398         H5E_auto2_t funcOld;
00399         void *dataOld;
00400         bool exists;
00401 
00402         H5Eget_auto2(H5E_DEFAULT, &funcOld, &dataOld);
00403         H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
00404 
00405         hid_t handle = H5Dopen2(h5file, path.c_str(), H5P_DEFAULT);
00406 
00407         exists = handle >= 0;
00408         if (exists) {
00409             H5Dclose(handle);
00410         }
00411         H5Eset_auto2(H5E_DEFAULT, funcOld, dataOld);
00412         return exists;
00413     }
00414 
00415     bool
00416     attr_exists(std::string path,
00417                 std::string attr)
00418     {
00419         if (group_exists(path)) {
00420             hid_t objHandle = H5Oopen(h5file, path.c_str(), H5P_DEFAULT);
00421             int exists = H5LTfind_attribute(objHandle, attr.c_str());
00422             H5Oclose(objHandle);
00423             return exists;
00424         }
00425         else {
00426             return false;
00427         }
00428     }
00429 
00430     void
00431     get_group_listing(std::string path,
00432                       std::vector<std::string> &listing)
00433     {
00434         hid_t groupHandle;
00435         hsize_t ix = 0;
00436         listing.clear();
00437         H5_CALL_ASSIGN(groupHandle, H5Gopen(h5file, path.c_str(), H5P_DEFAULT));
00438         H5_CALL(H5Literate(groupHandle, H5_INDEX_NAME, H5_ITER_INC, &ix, get_group_list__proc, &listing));
00439         H5_CALL(H5Gclose(groupHandle));
00440     }
00441 
00442     template <typename T>
00443     void
00444     write_ds(std::string path,
00445              std::string name,
00446              const std::vector<hsize_t>& shape,
00447              const std::vector<T>& data)
00448     {
00449         hid_t groupHandle;
00450         hsize_t dims[H5S_MAX_RANK];
00451         std::copy(shape.begin(), shape.end(), dims);
00452 
00453         H5_CALL_ASSIGN(groupHandle, H5Gopen(h5file, path.c_str(), H5P_DEFAULT));
00454         H5_CALL(H5LTmake_dataset(groupHandle, name.c_str(), shape.size(), dims,
00455                     H5File_specialized::get_type_id<T>(), &(data[0])));
00456         H5_CALL(H5Gclose(groupHandle));
00457     }
00458 };
00459 
00460 class LMPlugin {
00461     H5File h5file;
00462     std::string filename;
00463     int replicate;
00464 
00465     std::vector<std::string> speciesNames, siteNames;
00466     std::vector<std::string> dsNames;
00467 
00468     unsigned int totalRdmeAtom, totalSiteAtom, nAtom;
00469     std::vector<unsigned int> nRdmeAtoms, nSiteAtoms;
00470     std::vector<unsigned int> rdmeAtomOffsets, siteAtomOffsets;
00471 
00472     std::vector<double> frameTimes;
00473     std::vector<particle_int_t> pLattice;
00474     std::vector<site_int_t> sLattice;
00475 
00476     std::vector<unsigned int> trajIndices;
00477 
00478     unsigned int nSpecies, nSiteType;
00479     unsigned int nLatticeX, nLatticeY, nLatticeZ, nLatticeP;
00480     float latticeSpacing;
00481     float rdmeRadius;
00482     unsigned int frameIx, nFrame;
00483     bool hasVariableSites, hasTrajData;
00484 
00485     int bondFrom, bondTo;
00486 
00487     static std::string
00488     disambiguate_name(std::string name,
00489                       int dupCt)
00490     {
00491         if (name.size() < MAX_VMD_NAME && dupCt == 0) {
00492             return name;
00493         }
00494         else {
00495             char shortName[MAX_VMD_NAME];
00496             int numlen = dupCt == 0 ? 1 : log10(dupCt) + 1;
00497             int j;
00498 
00499             for (j = 0; j < MAX_VMD_NAME - numlen - 2; j++) { // 2 = 1(separator)+1(null)
00500                 shortName[j] = name[j];
00501             }
00502             snprintf(shortName + j, MAX_VMD_NAME, "~%d", dupCt);
00503 
00504             return std::string(shortName);
00505         }
00506     }
00507 
00508     void
00509     open_lm_file()
00510     {
00511         try {
00512             h5file.open(filename);
00513         }
00514         catch (H5Exception &exc) {
00515             throw Exception("Invalid LM file: Not an HDF5 file.");
00516         }
00517         std::vector<char> userblock;
00518         h5file.get_userblock(userblock);
00519 
00520         if (userblock.size() != 512 || strncmp("LMH5",  &userblock[0], 5) != 0) {
00521             throw Exception("Invalid LM file: Bad userblock");
00522         }
00523     }
00524 
00525     std::string
00526     sim_path(std::string relPath = "")
00527     {
00528         char path[128];
00529         snprintf(path, sizeof(path), "/Simulations/%07d/%s", replicate, relPath.c_str());
00530         return std::string(path);
00531     }
00532 
00533     void
00534     read_times()
00535     {
00536         std::vector<hsize_t> shape;
00537 
00538         h5file.get_dataset(sim_path("LatticeTimes"), frameTimes, shape);
00539         h5file.get_group_listing(sim_path("Lattice"), dsNames);
00540 
00541         // Verify site data consistent with rdme particle data
00542         if (hasVariableSites) {
00543             std::vector<double> times;
00544             std::vector<std::string> names;
00545             h5file.get_dataset(sim_path("SiteTimes"), times, shape);
00546             h5file.get_group_listing(sim_path("Sites"), names);
00547             if (dsNames != names || frameTimes != times) {
00548                 throw Exception("Particle lattice and site lattice data appear to be inconsistent");
00549             }
00550         }
00551 
00552         nFrame = dsNames.size();
00553     }
00554 
00555     void
00556     read_parameters()
00557     {
00558         h5file.get_attr("/Model/Diffusion", "numberSiteTypes", nSiteType);
00559         h5file.get_attr("/Model/Diffusion", "numberSpecies", nSpecies);
00560         h5file.get_attr("/Model/Diffusion", "latticeSpacing", latticeSpacing);
00561         h5file.get_attr("/Model/Diffusion", "latticeXSize", nLatticeX);
00562         h5file.get_attr("/Model/Diffusion", "latticeYSize", nLatticeY);
00563         h5file.get_attr("/Model/Diffusion", "latticeZSize", nLatticeZ);
00564         h5file.get_attr("/Model/Diffusion", "particlesPerSite", nLatticeP);
00565 
00566         hasVariableSites = h5file.dataset_exists(sim_path("Sites/0000000001"));
00567         hasTrajData = h5file.group_exists(sim_path("Lattice"));
00568 
00569         latticeSpacing *= TO_ANGSTROM;
00570         rdmeRadius = RDME_RADIUS_SCALE * latticeSpacing;
00571     }
00572 
00573     template <typename T>
00574     void
00575     count_lattice_objects(std::string path,
00576                           std::vector<unsigned int>& count,
00577                           std::vector<T>& data,
00578                           std::vector<unsigned int>& tmp)
00579     {
00580         std::vector<hsize_t> shape;
00581         hsize_t nData;
00582 
00583         h5file.get_dataset(path, data, shape, &nData);
00584 
00585         std::fill(tmp.begin(), tmp.end(), 0);
00586 
00587         for (unsigned int pIx = 0; pIx < nData; pIx++) {
00588             T type = data[pIx];
00589             if (type != 0) {
00590                 tmp[type]++;
00591             }
00592         }
00593 
00594         for (unsigned int i = 0; i < tmp.size(); i++) {
00595             count[i] = std::max(tmp[i], count[i]);
00596         }
00597     }
00598 
00599     void
00600     write_particle_counts()
00601     {
00602         h5file.close();
00603         h5file.open_rdwr(filename);
00604 
00605         std::vector<hsize_t> shape;
00606         shape.resize(1);
00607 
00608         shape[0] = nSpecies + 1;
00609         h5file.write_ds(sim_path(), "MaxParticleCounts", shape, nRdmeAtoms);
00610 
00611         shape[0] = nSiteType;
00612         h5file.write_ds(sim_path(), "MaxSiteCounts", shape, nSiteAtoms);
00613 
00614         h5file.close();
00615         h5file.open(filename);
00616 
00617         vmdcon_printf(
00618             VMDCON_INFO,
00619             "LMPlugin) Maximum particle counts written to %s:%s\n",
00620             filename.c_str(), sim_path().c_str());
00621     }
00622 
00623     void
00624     calculate_particle_counts()
00625     {
00626         std::vector<unsigned int> tmp;
00627 
00628         Timer timer;
00629         float t;
00630 
00631         // Count RDME max particles. If there is no trajectory stored in the file,
00632         // load the initial conditions. The maximum particle counts will not be
00633         // saved to the HDF5 file.
00634         tmp.resize(nSpecies + 1);
00635 
00636         timer.reset();
00637         if (hasTrajData) {
00638             vmdcon_printf(VMDCON_WARN,
00639                           "LMPlugin) Missing maximum particle count data. Rebuilding it now.\n");
00640 
00641             for (unsigned int fIx = 0, nproc = 0; fIx < dsNames.size(); fIx++, nproc++) {
00642                 if ((t = timer.elapsed()) > 60) {
00643                     vmdcon_printf(VMDCON_INFO,
00644                                   "LMPlugin) %lu/%lu particle lattices processed (%.1f/sec).\n", fIx,
00645                                   dsNames.size(), nproc / t);
00646                     nproc = 0;
00647                     timer.reset();
00648                 }
00649 
00650                 count_lattice_objects(sim_path("Lattice/" + dsNames[fIx]), nRdmeAtoms, pLattice, tmp);
00651             }
00652         }
00653         else {
00654             count_lattice_objects("Model/Diffusion/Lattice", nRdmeAtoms, pLattice, tmp);
00655         }
00656 
00657         vmdcon_printf(VMDCON_INFO, "LMPlugin) RDME particles counted.\n");
00658 
00659         // Count site max particles. Use either the static site lattice or
00660         // variable data if available
00661         tmp.resize(nSiteType + 1);
00662 
00663         timer.reset();
00664         if (hasVariableSites) {
00665             for (unsigned int fIx = 0, nproc = 0; fIx < dsNames.size(); fIx++, nproc++) {
00666                 if ((t = timer.elapsed()) > 60) {
00667                     vmdcon_printf(VMDCON_INFO,
00668                                   "LMPlugin) %lu/%lu site lattices processed (%.1f/sec).\n", fIx,
00669                                   dsNames.size(), nproc / t);
00670                     nproc = 0;
00671                     timer.reset();
00672                 }
00673                 count_lattice_objects(sim_path("Sites/" + dsNames[fIx]), nSiteAtoms, sLattice, tmp);
00674             }
00675             vmdcon_printf(VMDCON_INFO, "LMPlugin) Site particles counted (time dependent).\n");
00676         }
00677         else {
00678             count_lattice_objects("/Model/Diffusion/LatticeSites", nSiteAtoms, sLattice, tmp);
00679             vmdcon_printf(VMDCON_INFO, "LMPlugin) Site particles counted (constant).\n");
00680         }
00681 
00682         // Write histograms to the trajectory file so we don't have to do it again next time.
00683         if (hasTrajData) {
00684             write_particle_counts();
00685         }
00686     }
00687 
00688     void
00689     make_offsets(unsigned int &total,
00690                  std::vector<unsigned int> &maxCts,
00691                  std::vector<unsigned int> &offsets)
00692     {
00693         // The offset vector is constructed such that particle/site ids can be used as indexes.
00694         // Since we ignore the id=0 particles and sites, we explicitly set it to zero in case
00695         // they got counted anyway.
00696         maxCts[0] = 0;
00697         total = std::accumulate(maxCts.begin(), maxCts.end(), 0);
00698 
00701         offsets.resize(maxCts.size());
00702         offsets[0] = 0;
00703         std::partial_sum(maxCts.begin(), maxCts.end() - 1, offsets.begin() + 1);
00704     }
00705 
00706     void
00707     read_particle_counts()
00708     {
00709         std::vector<hsize_t> shape;
00710 
00711         nRdmeAtoms.resize(nSpecies + 1);
00712         std::fill(nRdmeAtoms.begin(), nRdmeAtoms.end(), 0);
00713         nSiteAtoms.resize(nSiteType + 1);
00714         std::fill(nSiteAtoms.begin(), nSiteAtoms.end(), 0);
00715 
00716         if (h5file.dataset_exists(sim_path("MaxParticleCounts")) ) {
00717             h5file.get_dataset(sim_path("MaxParticleCounts"), nRdmeAtoms, shape);
00718             h5file.get_dataset(sim_path("MaxSiteCounts"), nSiteAtoms, shape);
00719         }
00720         else {
00721             calculate_particle_counts();
00722         }
00723 
00724         // make offsets for each particle type into buffer provided by read_timestep
00725         make_offsets(totalRdmeAtom, nRdmeAtoms, rdmeAtomOffsets);
00726         make_offsets(totalSiteAtom, nSiteAtoms, siteAtomOffsets);
00727 
00728         // sites are placed after particles in the read_timestep buffer
00729         for (unsigned int i = 0; i < siteAtomOffsets.size(); i++) {
00730             siteAtomOffsets[i] += totalRdmeAtom;
00731         }
00732 
00733         nAtom = totalRdmeAtom + totalSiteAtom;
00734 
00735         // Extra element since the species counts do not include id=0
00736         trajIndices.resize(std::max(nSpecies + 1, nSiteType));
00737     }
00738 
00739     void
00740     load_and_fix_names(std::vector<std::string>& names,
00741                        const std::vector<std::string> &oldNames,
00742                        const char* type)
00743     {
00744         std::set<std::string> seen;
00745 
00746         for (unsigned int i = 0; i < oldNames.size(); i++) {
00747             std::string name = oldNames[i];
00748             std::string name0 = name;
00749             for (int dupCt = 0; seen.count(name) > 0; dupCt++) {
00750                 name = disambiguate_name(name0, dupCt);
00751             }
00752             seen.insert(name);
00753             names.push_back(name);
00754             if (name != name0) {
00755                 vmdcon_printf(VMDCON_WARN, "LMPlugin) %s ``%s'' renamed to ``%s''\n",
00756                               type, name0.c_str(), name.c_str());
00757             }
00758         }
00759     }
00760 
00761     void
00762     read_names()
00763     {
00764         std::vector<std::string> nameTmp;
00765         h5file.get_csv_attr("/Parameters", "speciesNames", nameTmp);
00766         speciesNames.push_back(""); // blank for id=0
00767 
00768         // ensure all names are unique after truncation to MAX_VMD_NAME characters
00769         load_and_fix_names(speciesNames, nameTmp, "Species");
00770 
00771         if (h5file.attr_exists("/Parameters", "siteTypeNames")) {
00772             h5file.get_csv_attr("/Parameters", "siteTypeNames", nameTmp);
00773         }
00774         else { // older LM files are missing this data
00775             nameTmp.clear();
00776             for (unsigned int i = 0; i < nSiteType; i++) {
00777                 std::stringstream ss;
00778                 ss << "site" << i;
00779                 nameTmp.push_back(ss.str());
00780             }
00781         }
00782 
00783         load_and_fix_names(siteNames, nameTmp, "Site");
00784     }
00785 
00786     void
00787     define_molfile_atom(molfile_atom_t *&atomPtr,
00788                         unsigned int ix,
00789                         const char *segname,
00790                         const char *fmt,
00791                         float radius,
00792                         const std::vector<std::string >& names,
00793                         const std::vector<unsigned int>& counts)
00794     {
00795         unsigned int ncopies = counts[ix];
00796         if (ncopies > 0) {
00797             molfile_atom_t atomType;
00798             memset(&atomType, 0, sizeof(atomType));
00799 
00800             std::copy(names[ix].begin(), names[ix].end(), atomType.name);
00801             std::strcpy(atomType.segid, segname);
00802             snprintf(atomType.type, sizeof(atomType.type), fmt, ix);
00803 
00804             atomType.resid = 1;
00805             atomType.radius = radius;
00806 
00807             std::fill(atomPtr, atomPtr + ncopies, atomType);
00808             atomPtr += ncopies;
00809         }
00810     }
00811 
00812     static inline
00813     float
00814     lattice_jitter(unsigned int co,
00815                    int axis)
00816     {
00817         // We're using this custom prng so that particles which do not transition between subvolumes
00818         // between frames do not appear to move
00819         co = (3*co + axis)*0xabcd1234;
00820 
00821         // xorshift32
00822         co ^= co << 13;
00823         co ^= co >> 17;
00824         co ^= co << 5;
00825         co &= 0xffffffff;
00826 
00827         // Scale/translate the noise vector to keep the particles inside a subvolume.
00828         // Here, 2*RDME_RADIUS_SCALE*latticeSpacing = diameter.
00829         return 2.328306e-10f*(1.0f - 2.0f * RDME_RADIUS_SCALE) * co + RDME_RADIUS_SCALE;
00830     }
00831 
00832   public:
00833     LMPlugin(const char *file)
00834         : filename(file)
00835         , replicate(0)
00836         , totalRdmeAtom(0)
00837         , totalSiteAtom(0)
00838         , frameIx(0)
00839         , nFrame(0)
00840         , bondFrom(0)
00841         , bondTo(0)
00842     {
00843         if (getenv("LM_REPLICATE") != NULL) {
00844             replicate = atoi(getenv("LM_REPLICATE"));
00845         }
00846 
00847         // replicate 0 shouldn't exist, if atoi fails we'll get 0
00848         if (replicate <= 0) {
00849             replicate = 1;
00850         }
00851 
00852         open_lm_file();
00853 
00854         // set nSiteType, nSpecies, latticeSpacing, nLatticeX, nLatticeY, nLatticeZ, nLatticeP,
00855         // hasVariableSites, hasTrajData, rdmeRadius
00856         read_parameters();
00857 
00858         // set frameTimes, dsNames, nFrame
00859         if (hasTrajData) {
00860             read_times();
00861             vmdcon_printf(VMDCON_INFO, "LMPlugin) %s, replicate %d:\n", file, replicate);
00862             vmdcon_printf(VMDCON_INFO, "LMPlugin)    Lattice dimensions: %dx%dx%dx%d\n", nLatticeX, nLatticeY,
00863                           nLatticeZ, nLatticeP);
00864             vmdcon_printf(VMDCON_INFO, "LMPlugin)    %d frames, %.3f to %.3f seconds.\n", nFrame,
00865                           frameTimes[0], frameTimes[nFrame - 1]);
00866             if (hasVariableSites) {
00867                 vmdcon_printf(VMDCON_INFO, "LMPlugin)    Contains site type trajectory\n");
00868             }
00869         }
00870         else {
00871             vmdcon_printf(VMDCON_INFO, "LMPlugin) %s, initial conditions (replicate %d missing)\n",
00872                           file, replicate);
00873             vmdcon_printf(VMDCON_INFO, "LMPlugin)    Lattice dimensions: %dx%dx%dx%d\n", nLatticeX, nLatticeY,
00874                           nLatticeZ, nLatticeP);
00875             frameTimes.push_back(0.0);
00876             nFrame = 1;
00877         }
00878 
00879         // set totalRdmeAtom, nRdmeAtoms, rdmeAtomOffsets,
00880         //     totalSiteAtom, nSiteAtoms, siteAtomOffsets,
00881         //     nAtom, trajIndices
00882         read_particle_counts();
00883 
00884         // set speciesNames, siteNames
00885         read_names();
00886 
00887         vmdcon_printf(VMDCON_INFO, "LMPlugin) Atom usage:\n");
00888         vmdcon_printf(VMDCON_INFO, "LMPlugin)    RDME species types: %3d (%d atoms used)\n",
00889                       nSpecies, totalRdmeAtom);
00890         vmdcon_printf(VMDCON_INFO, "LMPlugin)    Site types:         %3d (%d atoms used)\n",
00891                       nSiteType, totalSiteAtom);
00892     }
00893 
00894     unsigned int get_nframe() const { return nFrame; }
00895     unsigned int get_natoms() const { return nAtom; }
00896 
00897     int
00898     read_structure(int *optflags,
00899                    molfile_atom_t *atoms)
00900     {
00901         *optflags = MOLFILE_RADIUS;
00902         molfile_atom_t *atomPtr = atoms;
00903 
00904         // Segid field is used to distinguish particles (RDME) and sites (SITE).
00905 
00906         for (unsigned int sp = 0; sp < speciesNames.size(); sp++) {
00907             define_molfile_atom(atomPtr, sp, "RDME", "s%u", rdmeRadius, speciesNames, nRdmeAtoms);
00908         }
00909 
00910         for (unsigned int st = 0; st < siteNames.size(); st++) {
00911             define_molfile_atom(atomPtr, st, "SITE", "st%u", 0.5f*latticeSpacing, siteNames, nSiteAtoms);
00912         }
00913         return MOLFILE_SUCCESS;
00914     }
00915 
00916     int
00917     read_bonds(int *nbonds,
00918                int **from,
00919                int **to,
00920                float **bondorder,
00921                int **bondtype,
00922                int *nbondtypes,
00923                char ***bondtypename)
00924     {
00925         // Defined to prevent VMD from trying to create bonds on its own
00926         *nbonds = 0;
00927         *from = &(bondFrom);
00928         *to = &(bondTo);
00929         *bondorder = NULL;
00930         *bondtype = NULL;
00931         *nbondtypes = 0;
00932         *bondtypename = NULL;
00933         return MOLFILE_SUCCESS;
00934     }
00935 
00936     int
00937     read_timestep(molfile_timestep_t *ts)
00938     {
00939         if (frameIx >= nFrame) {
00940             return MOLFILE_EOF;
00941         }
00942 
00943         if (ts) { // ts can be null if we're skipping a frame
00944             ts->physical_time = frameTimes[frameIx];
00945             std::vector<hsize_t> shape;
00946             unsigned int ix;
00947 
00948             // Unused atoms banished to {-lambda, -lambda, -lambda}
00949             std::fill(ts->coords, ts->coords + 3 * nAtom, -1.0f * latticeSpacing);
00950 
00951             // RDME particles
00952             if (hasTrajData) {
00953                 h5file.get_dataset(sim_path("Lattice/" + dsNames[frameIx]), pLattice, shape);
00954             }
00955             else {
00956                 h5file.get_dataset("/Model/Diffusion/Lattice", pLattice, shape);
00957             }
00958 
00959             // Keep track of the offset of to the next empty slot for each particle type
00960             // using a histogram
00961             std::fill(trajIndices.begin(), trajIndices.end(), 0);
00962 
00963             ix = 0;
00964             for (unsigned int i = 0; i < shape[0]; i++) {
00965                 for (unsigned int j = 0; j < shape[1]; j++) {
00966                     for (unsigned int k = 0; k < shape[2]; k++) {
00967                         for (unsigned int p = 0; p < shape[3]; p++, ix++) {
00968                             particle_int_t type = pLattice[ix];
00969                             if (type != 0) {
00970                                 unsigned int offset = 3 * (rdmeAtomOffsets[type] + trajIndices[type]);
00971                                 float *coord = ts->coords + offset;
00972                                 // Semi-deterministic position randomization to break up moire
00973                                 // pattern
00974                                 coord[0] = latticeSpacing * (i + lattice_jitter(ix, 0));
00975                                 coord[1] = latticeSpacing * (j + lattice_jitter(ix, 1));
00976                                 coord[2] = latticeSpacing * (k + lattice_jitter(ix, 2));
00977                                 trajIndices[type]++;
00978                             }
00979                         }
00980                     }
00981                 }
00982             }
00983 
00984             // site types
00985             if (hasVariableSites) {
00986                 h5file.get_dataset(sim_path("Sites/" + dsNames[frameIx]), sLattice, shape);
00987             }
00988             else {
00989                 h5file.get_dataset("/Model/Diffusion/LatticeSites", sLattice, shape);
00990             }
00991 
00992             std::fill(trajIndices.begin(), trajIndices.end(), 0);
00993 
00994             ix = 0;
00995             for (unsigned int i = 0; i < shape[0]; i++) {
00996                 for (unsigned int j = 0; j < shape[1]; j++) {
00997                     for (unsigned int k = 0; k < shape[2]; k++, ix++) {
00998                         site_int_t type = sLattice[ix];
00999                         if (type != 0) {
01000                             unsigned int offset = 3 * (siteAtomOffsets[type] + trajIndices[type]);
01001                             float *coord = ts->coords + offset;
01002                             coord[0] = latticeSpacing * (i + 0.5f);
01003                             coord[1] = latticeSpacing * (j + 0.5f);
01004                             coord[2] = latticeSpacing * (k + 0.5f);
01005                             trajIndices[type]++;
01006                         }
01007                     }
01008                 }
01009             }
01010         }
01011 
01012         frameIx++;
01013         return MOLFILE_SUCCESS;
01014     }
01015 };
01016 } // end anonymous namespace
01017 
01018 static void *
01019 open_read(const char *filename,
01020           const char *filetype,
01021           int *nAtom)
01022 {
01023     H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
01024     LMPlugin *lmd = NULL;
01025     try {
01026         lmd = new LMPlugin(filename);
01027         *nAtom = lmd->get_natoms();
01028         return lmd;
01029     }
01030     catch (LMPException &e) {
01031         e.vmdcon_report();
01032         return NULL;
01033     }
01034 }
01035 
01036 static int
01037 read_structure(void *mydata,
01038                int *optflags,
01039                molfile_atom_t *atoms)
01040 {
01041     LMPlugin *lmd = static_cast<LMPlugin *>(mydata);
01042     try {
01043         return lmd->read_structure(optflags, atoms);
01044     }
01045     catch (LMPException &e) {
01046         e.vmdcon_report();
01047         return MOLFILE_ERROR;
01048     }
01049 }
01050 
01051 static int
01052 read_bonds(void *mydata,
01053            int *nbonds,
01054            int **from,
01055            int **to,
01056            float **bondorder,
01057            int **bondtype,
01058            int *nbondtypes,
01059            char ***bondtypename)
01060 {
01061     LMPlugin *lmd = static_cast<LMPlugin *>(mydata);
01062     try {
01063         return lmd->read_bonds(nbonds, from, to, bondorder, bondtype, nbondtypes, bondtypename);
01064     }
01065     catch (LMPException &e) {
01066         e.vmdcon_report();
01067         return MOLFILE_ERROR;
01068     }
01069 }
01070 
01071 static int
01072 read_timestep(void *mydata,
01073               int nAtom,
01074               molfile_timestep_t *ts)
01075 {
01076     LMPlugin *lmd = static_cast<LMPlugin *>(mydata);
01077     try {
01078         return lmd->read_timestep(ts);
01079     }
01080     catch (LMPException &e) {
01081         e.vmdcon_report();
01082         return MOLFILE_ERROR;
01083     }
01084 }
01085 
01086 static void
01087 close_read(void *mydata)
01088 {
01089     LMPlugin *lmd = static_cast<LMPlugin *>(mydata);
01090     delete lmd;
01091 }
01092 
01093 VMDPLUGIN_API int
01094 VMDPLUGIN_init()
01095 {
01096     memset(&plugin, 0, sizeof(molfile_plugin_t));
01097     plugin.abiversion = vmdplugin_ABIVERSION;
01098     plugin.type = MOLFILE_PLUGIN_TYPE;
01099     plugin.name = "LM";
01100     plugin.prettyname = "Lattice Microbes (RDME)";
01101     plugin.author = "Tyler M. Earnest";
01102     plugin.majorv = 1;
01103     plugin.minorv = 0;
01104     plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
01105     plugin.filename_extension = "lm";
01106     plugin.open_file_read = open_read;
01107     plugin.read_structure = read_structure;
01108     plugin.read_bonds = read_bonds;
01109     plugin.read_next_timestep = read_timestep;
01110     plugin.close_file_read = close_read;
01111     return VMDPLUGIN_SUCCESS;
01112 }
01113 
01114 VMDPLUGIN_API int
01115 VMDPLUGIN_register(void *v,
01116                    vmdplugin_register_cb cb)
01117 {
01118     (*cb)(v, (vmdplugin_t *) &plugin);
01119     return VMDPLUGIN_SUCCESS;
01120 }
01121 
01122 VMDPLUGIN_API int
01123 VMDPLUGIN_fini()
01124 {
01125     return VMDPLUGIN_SUCCESS;
01126 }
01127 
01128 #ifdef LMP_MAIN
01129 int
01130 main(int argc,
01131      char *argv[])
01132 {
01133     if (argc != 2) {
01134         fprintf(stderr, "usage: %s test.lm\n", argv[0]);
01135         return 1;
01136     }
01137     VMDPLUGIN_init();
01138     int nAtom, flags;
01139     LMPlugin *lmp = static_cast<LMPlugin *>(open_read(argv[1], NULL, &nAtom));
01140     if (!lmp) {
01141         fprintf(stderr, "%s: open_read() failed\n", argv[1]);
01142         return 1;
01143     }
01144     molfile_atom_t *atoms = new molfile_atom_t[nAtom];
01145     if (read_structure(lmp, &flags, atoms) != MOLFILE_SUCCESS) {
01146         fprintf(stderr, "%s: read_structure() failed\n", argv[1]);
01147         return 1;
01148     }
01149 
01150     int nbonds;
01151     int *from, *to, *bondtype, nbondtypes;
01152     float *bondorder;
01153     char **bondtypename;
01154     if (read_bonds(lmp, &nbonds, &from, &to, &bondorder, &bondtype, &nbondtypes, &bondtypename)
01155             != MOLFILE_SUCCESS) {
01156         fprintf(stderr, "%s: read_bonds() failed\n", argv[1]);
01157         return 1;
01158     }
01159 
01160     molfile_timestep_t trajdata;
01161     trajdata.coords = new float[3*nAtom];
01162 
01163     int stride = lmp->get_nframe()/3;
01164     stride = stride > 0 ? stride : 1;
01165     int nframes = stride*((lmp->get_nframe() + stride - 1)/stride);
01166 
01167     for (unsigned int i=0; i < nframes; i+=stride) {
01168         if (read_timestep(lmp, nAtom, &trajdata) != MOLFILE_SUCCESS) {
01169             fprintf(stderr, "%s: read_timestep(), frame index %d failed\n", argv[1], i);
01170             return 1;
01171         }
01172         printf("frame %d/%d\n", i+1, lmp->get_nframe());
01173     }
01174 
01175     close_read(lmp);
01176     delete [] trajdata.coords;
01177     delete [] atoms;
01178     return 0;
01179 }
01180 #endif

Generated on Fri Jan 17 02:57:46 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002