Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80967907
molfile_interface.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Sep 3, 12:49
Size
33 KB
Mime Type
text/x-c
Expires
Thu, Sep 5, 12:49 (2 d)
Engine
blob
Format
Raw Data
Handle
20482871
Attached To
rLAMMPS lammps
molfile_interface.cpp
View Options
This document is not UTF8. It was detected as Shift JIS and converted to UTF8 for display.
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple)
------------------------------------------------------------------------- */
#include "molfile_interface.h"
#include <sys/types.h>
#include <stdio.h>
#include <dirent.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>
#if defined(_WIN32)
#include <windows.h>
#else
#include <dirent.h>
#include <dlfcn.h>
#endif
#include "molfile_plugin.h"
#if vmdplugin_ABIVERSION < 16
#error "unsupported VMD molfile plugin ABI version"
#endif
#define DEBUG 0
extern "C" {
typedef int (*initfunc)(void);
typedef int (*regfunc)(void *, vmdplugin_register_cb);
typedef int (*finifunc)(void);
typedef struct {
void *p;
const char *name;
} plugin_reginfo_t;
// callback function for plugin registration.
static int plugin_register_cb(void *v, vmdplugin_t *p)
{
plugin_reginfo_t *r = static_cast<plugin_reginfo_t *>(v);
// make sure we have the proper plugin type (native reader)
// for the desired file type (called "name" at this level)
if ((strcmp(MOLFILE_PLUGIN_TYPE,p->type) == 0)
&& (strcmp(r->name, p->name) == 0) ) {
r->p = static_cast<void *>(p);
}
return 0;
}
/* periodic table of elements for translation of ordinal to atom type */
static const char *pte_label[] = {
"X", "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne",
"Na", "Mg", "Al", "Si", "P" , "S", "Cl", "Ar", "K", "Ca", "Sc",
"Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge",
"As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc",
"Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe",
"Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb",
"Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os",
"Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr",
"Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf",
"Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt",
"Ds", "Rg"
};
static const int nr_pte_entries = sizeof(pte_label) / sizeof(char *);
/* corresponding table of masses. */
static const float pte_mass[] = {
/* X */ 0.00000, 1.00794, 4.00260, 6.941, 9.012182, 10.811,
/* C */ 12.0107, 14.0067, 15.9994, 18.9984032, 20.1797,
/* Na */ 22.989770, 24.3050, 26.981538, 28.0855, 30.973761,
/* S */ 32.065, 35.453, 39.948, 39.0983, 40.078, 44.955910,
/* Ti */ 47.867, 50.9415, 51.9961, 54.938049, 55.845, 58.9332,
/* Ni */ 58.6934, 63.546, 65.409, 69.723, 72.64, 74.92160,
/* Se */ 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585,
/* Zr */ 91.224, 92.90638, 95.94, 98.0, 101.07, 102.90550,
/* Pd */ 106.42, 107.8682, 112.411, 114.818, 118.710, 121.760,
/* Te */ 127.60, 126.90447, 131.293, 132.90545, 137.327,
/* La */ 138.9055, 140.116, 140.90765, 144.24, 145.0, 150.36,
/* Eu */ 151.964, 157.25, 158.92534, 162.500, 164.93032,
/* Er */ 167.259, 168.93421, 173.04, 174.967, 178.49, 180.9479,
/* W */ 183.84, 186.207, 190.23, 192.217, 195.078, 196.96655,
/* Hg */ 200.59, 204.3833, 207.2, 208.98038, 209.0, 210.0, 222.0,
/* Fr */ 223.0, 226.0, 227.0, 232.0381, 231.03588, 238.02891,
/* Np */ 237.0, 244.0, 243.0, 247.0, 247.0, 251.0, 252.0, 257.0,
/* Md */ 258.0, 259.0, 262.0, 261.0, 262.0, 266.0, 264.0, 269.0,
/* Mt */ 268.0, 271.0, 272.0
};
/*
* corresponding table of VDW radii.
* van der Waals radii are taken from A. Bondi,
* J. Phys. Chem., 68, 441 - 452, 1964,
* except the value for H, which is taken from R.S. Rowland & R. Taylor,
* J.Phys.Chem., 100, 7384 - 7391, 1996. Radii that are not available in
* either of these publications have RvdW = 2.00 ナ.
* The radii for Ions (Na, K, Cl, Ca, Mg, and Cs are based on the CHARMM27
* Rmin/2 parameters for (SOD, POT, CLA, CAL, MG, CES) by default.
*/
static const float pte_vdw_radius[] = {
/* X */ 1.5, 1.2, 1.4, 1.82, 2.0, 2.0,
/* C */ 1.7, 1.55, 1.52, 1.47, 1.54,
/* Na */ 1.36, 1.18, 2.0, 2.1, 1.8,
/* S */ 1.8, 2.27, 1.88, 1.76, 1.37, 2.0,
/* Ti */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
/* Ni */ 1.63, 1.4, 1.39, 1.07, 2.0, 1.85,
/* Se */ 1.9, 1.85, 2.02, 2.0, 2.0, 2.0,
/* Zr */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
/* Pd */ 1.63, 1.72, 1.58, 1.93, 2.17, 2.0,
/* Te */ 2.06, 1.98, 2.16, 2.1, 2.0,
/* La */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
/* Eu */ 2.0, 2.0, 2.0, 2.0, 2.0,
/* Er */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
/* W */ 2.0, 2.0, 2.0, 2.0, 1.72, 1.66,
/* Hg */ 1.55, 1.96, 2.02, 2.0, 2.0, 2.0, 2.0,
/* Fr */ 2.0, 2.0, 2.0, 2.0, 2.0, 1.86,
/* Np */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
/* Md */ 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
/* Mt */ 2.0, 2.0, 2.0
};
/* lookup functions */
static const char *get_pte_label(const int idx)
{
if ((idx < 1) || (idx >= nr_pte_entries)) return pte_label[0];
return pte_label[idx];
}
static float get_pte_mass(const int idx)
{
if ((idx < 1) || (idx >= nr_pte_entries)) return pte_mass[0];
return pte_mass[idx];
}
static float get_pte_vdw_radius(const int idx)
{
if ((idx < 1) || (idx >= nr_pte_entries)) return pte_vdw_radius[0];
#if 1
/* Replace with Hydrogen radius with an "all-atom" radius */
if (idx == 1)
return 1.0; /* H */
#else
/* Replace with old VMD atom radii values */
switch (idx) {
case 1: return 1.0; /* H */
case 6: return 1.5; /* C */
case 7: return 1.4; /* N */
case 8: return 1.3; /* O */
case 9: return 1.2; /* F */
case 15: return 1.5; /* P */
case 16: return 1.9; /* S */
}
#endif
return pte_vdw_radius[idx];
}
static int get_pte_idx(const char *label)
{
int i;
char atom[3];
/* zap string */
atom[0] = (char) 0;
atom[1] = (char) 0;
atom[2] = (char) 0;
/* if we don't have a null-pointer, there must be at least two
* chars, which is all we need. we convert to the capitalization
* convention of the table above during assignment. */
if (label != NULL) {
atom[0] = (char) toupper((int) label[0]);
atom[1] = (char) tolower((int) label[1]);
}
/* discard numbers in atom label */
if (isdigit(atom[1])) atom[1] = (char) 0;
for (i=0; i < nr_pte_entries; ++i) {
if ( (pte_label[i][0] == atom[0])
&& (pte_label[i][1] == atom[1]) ) return i;
}
return 0;
}
static int get_pte_idx_from_string(const char *label) {
int i, ind;
char atom[3];
if (label != NULL) {
/* zap string */
atom[0] = atom[1] = atom[2] = '\0';
for (ind=0,i=0; (ind<2) && (label[i]!='\0'); i++) {
if (label[i] != ' ') {
atom[ind] = toupper(label[i]);
ind++;
}
}
if (ind < 1)
return 0; /* no non-whitespace characters */
for (i=0; i < nr_pte_entries; ++i) {
if ((toupper(pte_label[i][0]) == atom[0]) && (toupper(pte_label[i][1]) == atom[1]))
return i;
}
}
return 0;
}
// directory traversal helper functions
#if defined(_WIN32)
// Win32 directory traversal handle
typedef struct {
HANDLE h;
WIN32_FIND_DATA fd;
char *name;
char *searchname;
int dlen;
} dirhandle_t;
// open a directory handle
static dirhandle_t *my_opendir(const char *dirname)
{
dirhandle_t *d;
int len;
if (dirname == NULL)
return NULL;
d = new dirhandle_t;
len = 2 + strlen(dirname);
d->name = new char[len];
strcpy(d->name, dirname);
strcat(d->name, "\\");
d->dlen = len;
len += 1;
d->searchname = new char[len];
strcpy(d->searchname, dirname);
strcat(d->searchname, "\\*");
d->h = FindFirstFile(d->searchname, &(d->fd));
if (d->h == ((HANDLE)(-1))) {
delete[] d->searchname;
delete[] d->name;
delete d;
return NULL;
}
return d;
}
// get next file name from directory handle
static char *my_readdir(dirhandle_t *d)
{
if (FindNextFile(d->h, &(d->fd))) {
return d->fd.cFileName;
}
return NULL;
}
// close directory handle
static void my_closedir(dirhandle_t *d)
{
if (d->h != NULL) {
FindClose(d->h);
}
delete[] d->searchname;
delete[] d->name;
delete d;
}
// open a shared object file
static void *my_dlopen(const char *fname) {
return (void *)LoadLibrary(fname);
}
// report error message from dlopen
static const char *my_dlerror(void) {
static CHAR szBuf[80];
DWORD dw = GetLastError();
sprintf(szBuf, "my_dlopen failed: GetLastError returned %u\n", dw);
return szBuf;
}
// resolve a symbol in shared object
static void *my_dlsym(void *h, const char *sym) {
return (void *)GetProcAddress((HINSTANCE)h, sym);
}
// close a shared object
static int my_dlclose(void *h) {
/* FreeLibrary returns nonzero on success */
return !FreeLibrary((HINSTANCE)h);
}
#else
// Unix directory traversal handle
typedef struct {
DIR *d;
char *name;
int dlen;
} dirhandle_t;
// open a directory handle
static dirhandle_t *my_opendir(const char *dirname)
{
dirhandle_t *d;
int len;
if (dirname == NULL) return NULL;
d = new dirhandle_t;
len = 2 + strlen(dirname);
d->name = new char[len];
strcpy(d->name,dirname);
strcat(d->name,"/");
d->dlen = len;
d->d = opendir(d->name);
if (d->d == NULL) {
delete[] d->name;
delete d;
return NULL;
}
return d;
}
// get next file name from directory handle
static char *my_readdir(dirhandle_t *d)
{
struct dirent *p;
if ((p = readdir(d->d)) != NULL) {
return p->d_name;
}
return NULL;
}
// close directory handle
static void my_closedir(dirhandle_t *d)
{
if (d->d != NULL) {
closedir(d->d);
}
delete[] d->name;
delete d;
return;
}
// open a shared object file
static void *my_dlopen(const char *fname) {
return dlopen(fname, RTLD_NOW);
}
// report error message from dlopen
static const char *my_dlerror(void) {
return dlerror();
}
// resolve a symbol in shared object
static void *my_dlsym(void *h, const char *sym) {
return dlsym(h, sym);
}
// close a shared object
static int my_dlclose(void *h) {
return dlclose(h);
}
#endif
} // end of extern "C" region
using namespace LAMMPS_NS;
// constructor.
MolfileInterface::MolfileInterface(const char *type, const int mode)
: _plugin(0), _dso(0), _ptr(0), _info(0), _natoms(0),
_mode(mode), _caps(M_NONE)
{
_name = new char[5];
strcpy(_name,"none");
_type = new char[1+strlen(type)];
strcpy(_type,type);
}
// destructor.
MolfileInterface::~MolfileInterface()
{
forget_plugin();
if (_info) {
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
delete[] a;
_info = NULL;
}
delete[] _name;
delete[] _type;
}
// register the best matching plugin in a given directory
int MolfileInterface::find_plugin(const char *pluginpath)
{
dirhandle_t *dir;
char *filename, *ext, *next, *path, *plugindir;
int retval = E_NONE;
#if defined(_WIN32)
#define MY_PATHSEP ';'
#else
#define MY_PATHSEP ':'
#endif
if (pluginpath == NULL) return E_DIR;
plugindir = path = strdup(pluginpath);
while (plugindir) {
// check if this a single directory or path.
next = strchr(plugindir,MY_PATHSEP);
if (next) {
*next = '\0';
++next;
}
dir = my_opendir(plugindir);
if (!dir)
retval = (retval > E_DIR) ? retval : E_DIR;
// search for suitable file names and try to inspect them
while(dir) {
char *fullname;
int len;
filename = my_readdir(dir);
if (filename == NULL) break;
// only look at .so files
ext = strrchr(filename, '.');
if (ext == NULL) continue;
if (strcasecmp(ext,".so") != 0) continue;
// construct full pathname of potential DSO
len = dir->dlen;
len += strlen(filename);
fullname = new char[len];
strcpy(fullname,dir->name);
strcat(fullname,filename);
// try to register plugin at file name.
int rv = load_plugin(fullname);
if (rv > retval) retval = rv;
delete[] fullname;
}
if (dir)
my_closedir(dir);
plugindir = next;
}
free(path);
return retval;
}
// register the best matching plugin in a given directory
int MolfileInterface::load_plugin(const char *filename)
{
void *dso;
int len, retval = E_NONE;
// access shared object
dso = my_dlopen(filename);
if (dso == NULL)
return E_FILE;
// check for required plugin symbols
void *ifunc = my_dlsym(dso,"vmdplugin_init");
void *rfunc = my_dlsym(dso,"vmdplugin_register");
void *ffunc = my_dlsym(dso,"vmdplugin_fini");
if (ifunc == NULL || rfunc == NULL || ffunc == NULL) {
my_dlclose(dso);
return E_SYMBOL;
}
// intialize plugin. skip plugin if it fails.
if (((initfunc)(ifunc))()) {
my_dlclose(dso);
return E_SYMBOL;
}
// pre-register plugin.
// the callback will be called for each plugin in the DSO and
// check the file type. plugin->name will change if successful.
plugin_reginfo_t reginfo;
reginfo.p = NULL;
reginfo.name=_type;
((regfunc)rfunc)(®info, plugin_register_cb);
// make some checks to see if the plugin is suitable or not.
molfile_plugin_t *plugin = static_cast<molfile_plugin_t *>(reginfo.p);
// if the callback found a matching plugin and copied the struct,
// its name element will point to a different location now.
if (plugin == NULL) {
retval = E_TYPE;
// check if the ABI matches the one used to compile this code
} else if (plugin->abiversion != vmdplugin_ABIVERSION) {
retval = E_ABI;
// check if (basic) reading is supported
} else if ((_mode & M_READ) &&
( (plugin->open_file_read == NULL) ||
(plugin->read_next_timestep == NULL) ||
(plugin->close_file_read == NULL) )) {
retval = E_MODE;
// check if (basic) writing is supported
} else if ( (_mode & M_WRITE) &&
( (plugin->open_file_write == NULL) ||
(plugin->write_timestep == NULL) ||
(plugin->close_file_write == NULL) )) {
retval = E_MODE;
// make some additional check, if we
// already have a plugin registered.
// NOTE: this has to come last.
} else if (_dso && _plugin) {
molfile_plugin_t *p;
p = static_cast<molfile_plugin_t *>(_plugin);
// check if the new plugin is of a newer major version
if (p->majorv > plugin->majorv) {
retval = E_VERSION;
// check if the new plugin is of a newer minor version
} else if ( (p->majorv == plugin->majorv) &&
(p->minorv >= plugin->minorv) ) {
retval = E_VERSION;
}
}
// bingo! this one is a keeper.
if (retval == E_NONE) {
// make sure any existing plugin is wiped out
forget_plugin();
delete[] _name;
len = 16;
len += strlen(plugin->prettyname);
len += strlen(plugin->author);
_name = new char[len];
sprintf(_name,"%s v%d.%d by %s",plugin->prettyname,
plugin->majorv, plugin->minorv, plugin->author);
// determine plugin capabilities
_caps = M_NONE;
if (plugin->read_next_timestep) _caps |= M_READ;
if (plugin->write_timestep) _caps |= M_WRITE;
#if vmdplugin_ABIVERSION > 10
// required to tell if velocities are present
if (plugin->read_timestep_metadata) _caps |= M_RVELS;
// we can always offer velocities. we may not know if
// they will be written by the plugin though.
if (plugin->write_timestep) _caps |= M_WVELS;
#endif
if (plugin->read_structure) _caps |= M_RSTRUCT;
if (plugin->write_structure) _caps |= M_WSTRUCT;
if (plugin->read_bonds) _caps |= M_RBONDS;
if (plugin->write_bonds) _caps |= M_WBONDS;
if (plugin->read_angles) _caps |= M_RANGLES;
if (plugin->write_angles) _caps |= M_WANGLES;
if (plugin->read_volumetric_data) _caps |= M_RVOL;
if (plugin->write_volumetric_data) _caps |= M_WVOL;
if (_mode & M_WRITE) {
_mode |= (_caps & M_WSTRUCT);
_mode |= (_caps & M_WVELS);
} else if (_mode & M_READ) {
_mode |= (_caps & M_RSTRUCT);
_mode |= (_caps & M_RVELS);
}
_plugin = plugin;
_dso = dso;
return E_MATCH;
}
// better luck next time. clean up and return.
my_dlclose(dso);
return retval;
}
// deregister a plugin and close or reset all associated objects.
void MolfileInterface::forget_plugin()
{
if (_ptr)
close();
if (_plugin)
_plugin = NULL;
if (_dso) {
void *ffunc = my_dlsym(_dso,"vmdplugin_fini");
if (ffunc)
((finifunc)ffunc)();
my_dlclose(_dso);
}
_dso = NULL;
delete[] _name;
_name = new char[5];
strcpy(_name,"none");
_caps = M_NONE;
}
// open file for reading or writing
int MolfileInterface::open(const char *name, int *natoms)
{
if (!_plugin || !_dso || !natoms)
return E_FILE;
molfile_plugin_t *p = static_cast<molfile_plugin_t *>(_plugin);
if (_mode & M_WRITE)
_ptr = p->open_file_write(name,_type,*natoms);
else if (_mode & M_READ)
_ptr = p->open_file_read(name,_type,natoms);
if (_ptr == NULL)
return E_FILE;
_natoms = *natoms;
// we need to deal with structure information,
// so we allocate and initialize storage for it.
if (_mode & (M_RSTRUCT|M_WSTRUCT)) {
molfile_atom_t *a = new molfile_atom_t[_natoms];
_info = a;
memset(_info,0,_natoms*sizeof(molfile_atom_t));
for (int i=0; i < _natoms; ++i) {
a[i].name[0] = 'X';
a[i].type[0] = a[i].resname[0] = a[i].segid[0] = 'U';
a[i].type[1] = a[i].resname[1] = a[i].segid[1] = 'N';
a[i].type[2] = a[i].resname[2] = a[i].segid[2] = 'K';
a[i].chain[0] = 'X';
}
}
return E_NONE;
}
// get of set atom structure information
int MolfileInterface::structure()
{
if (!_plugin || !_dso)
return E_FILE;
molfile_plugin_t *p = static_cast<molfile_plugin_t *>(_plugin);
int optflags = MOLFILE_NOOPTIONS;
if (_mode & M_WSTRUCT) {
optflags |= (_props & P_BFAC) ? MOLFILE_BFACTOR : 0;
optflags |= (_props & P_OCCP) ? MOLFILE_OCCUPANCY : 0;
optflags |= (_props & P_MASS) ? MOLFILE_MASS : 0;
optflags |= (_props & P_CHRG) ? MOLFILE_CHARGE : 0;
optflags |= (_props & P_RADS) ? MOLFILE_RADIUS : 0;
optflags |= (_props & P_ATMN) ? MOLFILE_ATOMICNUMBER : 0;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
p->write_structure(_ptr,optflags,a);
} else if (_mode & M_RSTRUCT) {
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
p->read_structure(_ptr,&optflags,a);
// mandatory properties
_props = P_NAME|P_TYPE|P_RESN|P_RESI|P_SEGN|P_CHAI;
// optional properties
_props |= (optflags & MOLFILE_BFACTOR) ? P_BFAC : 0;
_props |= (optflags & MOLFILE_OCCUPANCY) ? P_OCCP : 0;
_props |= (optflags & MOLFILE_MASS) ? P_MASS : 0;
_props |= (optflags & MOLFILE_CHARGE) ? P_CHRG : 0;
_props |= (optflags & MOLFILE_RADIUS) ? P_RADS : 0;
_props |= (optflags & MOLFILE_ATOMICNUMBER) ? P_ATMN : 0;
}
return 0;
}
// safely close file
int MolfileInterface::close()
{
if (!_plugin || !_dso || !_ptr)
return E_FILE;
molfile_plugin_t *p = static_cast<molfile_plugin_t *>(_plugin);
if (_mode & M_WRITE) {
p->close_file_write(_ptr);
} else if (_mode & M_READ) {
p->close_file_read(_ptr);
}
if (_info) {
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
delete[] a;
_info = NULL;
}
_ptr = NULL;
_natoms = 0;
return E_NONE;
}
// read or write timestep
int MolfileInterface::timestep(float *coords, float *vels,
float *cell, double *simtime)
{
if (!_plugin || !_dso || !_ptr)
return 1;
molfile_plugin_t *p = static_cast<molfile_plugin_t *>(_plugin);
molfile_timestep_t *t = new molfile_timestep_t;
int rv;
if (_mode & M_WRITE) {
t->coords = coords;
t->velocities = vels;
if (cell != NULL) {
t->A = cell[0];
t->B = cell[1];
t->C = cell[2];
t->alpha = cell[3];
t->beta = cell[4];
t->gamma = cell[5];
} else {
t->A = 0.0f;
t->B = 0.0f;
t->C = 0.0f;
t->alpha = 90.0f;
t->beta = 90.0f;
t->gamma = 90.0f;
}
if (simtime)
t->physical_time = *simtime;
else
t->physical_time = 0.0;
rv = p->write_timestep(_ptr,t);
} else {
// no coordinate storage => skip step
if (coords == NULL) {
rv = p->read_next_timestep(_ptr, _natoms, NULL);
} else {
t->coords = coords;
t->velocities = vels;
t->A = 0.0f;
t->B = 0.0f;
t->C = 0.0f;
t->alpha = 90.0f;
t->beta = 90.0f;
t->gamma = 90.0f;
t->physical_time = 0.0;
rv = p->read_next_timestep(_ptr, _natoms, t);
if (cell != NULL) {
cell[0] = t->A;
cell[1] = t->B;
cell[2] = t->C;
cell[3] = t->alpha;
cell[4] = t->beta;
cell[5] = t->gamma;
}
if (simtime)
*simtime = t->physical_time;
}
if (rv == MOLFILE_EOF)
return 1;
}
return 0;
}
// functions to read properties from molfile structure
#define PROPUPDATE(PROP,ENTRY,VAL) \
if (propid == PROP) { VAL = a.ENTRY; }
#define PROPSTRCPY(PROP,ENTRY,VAL) \
if (propid == PROP) { strcpy(VAL,a.ENTRY); }
// single precision floating point props
static float read_float_property(molfile_atom_t &a, const int propid)
{
float prop = 0.0f;
int iprop = 0;
PROPUPDATE(MolfileInterface::P_OCCP,occupancy,prop);
PROPUPDATE(MolfileInterface::P_BFAC,bfactor,prop);
PROPUPDATE(MolfileInterface::P_MASS,mass,prop);
PROPUPDATE(MolfileInterface::P_CHRG,charge,prop);
PROPUPDATE(MolfileInterface::P_RADS,radius,prop);
PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_MASS),
atomicnumber,iprop);
PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_RADS),
atomicnumber,iprop);
if (propid & MolfileInterface::P_ATMN) {
if (propid & MolfileInterface::P_MASS)
prop = get_pte_mass(iprop);
if (propid & MolfileInterface::P_RADS)
prop = get_pte_vdw_radius(iprop);
}
return prop;
}
// integer and derived props
static int read_int_property(molfile_atom_t &a, const int propid)
{
int prop = 0;
const char * sprop;
PROPUPDATE(MolfileInterface::P_RESI,resid,prop);
PROPUPDATE(MolfileInterface::P_ATMN,atomicnumber,prop);
PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_NAME),
name,sprop);
PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_TYPE),
type,sprop);
if (propid & MolfileInterface::P_ATMN) {
if (propid & (MolfileInterface::P_NAME|MolfileInterface::P_TYPE))
prop = get_pte_idx_from_string(sprop);
}
return prop;
}
// string and derived props
static const char *read_string_property(molfile_atom_t &a,
const int propid)
{
const char *prop = NULL;
int iprop = 0;
PROPUPDATE(MolfileInterface::P_NAME,name,prop);
PROPUPDATE(MolfileInterface::P_TYPE,type,prop);
PROPUPDATE(MolfileInterface::P_RESN,resname,prop);
PROPUPDATE(MolfileInterface::P_SEGN,segid,prop);
PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_NAME),
atomicnumber,iprop);
PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_TYPE),
atomicnumber,iprop);
if (propid & MolfileInterface::P_ATMN) {
if (propid & (MolfileInterface::P_NAME|MolfileInterface::P_TYPE))
prop = get_pte_label(iprop);
}
return prop;
}
#undef PROPUPDATE
#undef PROPSTRCPY
// functions to store properties into molfile structure
#define PROPUPDATE(PROP,ENTRY,VAL) \
if ((propid & PROP) == PROP) { a.ENTRY = VAL; plist |= PROP; }
#define PROPSTRCPY(PROP,ENTRY,VAL) \
if ((propid & PROP) == PROP) { strcpy(a.ENTRY,VAL); plist |= PROP; }
// floating point props
static int write_atom_property(molfile_atom_t &a,
const int propid,
const float prop)
{
int plist = MolfileInterface::P_NONE;
PROPUPDATE(MolfileInterface::P_OCCP,occupancy,prop);
PROPUPDATE(MolfileInterface::P_BFAC,bfactor,prop);
PROPUPDATE(MolfileInterface::P_MASS,mass,prop);
PROPUPDATE(MolfileInterface::P_CHRG,charge,prop);
PROPUPDATE(MolfileInterface::P_RADS,radius,prop);
return plist;
}
// double precision floating point props
static int write_atom_property(molfile_atom_t &a,
const int propid,
const double prop)
{
return write_atom_property(a,propid,static_cast<float>(prop));
}
// integer and derived props
static int write_atom_property(molfile_atom_t &a,
const int propid,
const int prop)
{
int plist = MolfileInterface::P_NONE;
PROPUPDATE(MolfileInterface::P_RESI,resid,prop);
PROPUPDATE(MolfileInterface::P_ATMN,atomicnumber,prop);
PROPUPDATE((MolfileInterface::P_ATMN|MolfileInterface::P_MASS),
mass,get_pte_mass(prop));
PROPSTRCPY((MolfileInterface::P_ATMN|MolfileInterface::P_NAME),
name,get_pte_label(prop));
PROPSTRCPY((MolfileInterface::P_ATMN|MolfileInterface::P_TYPE),
type,get_pte_label(prop));
return plist;
}
// integer and derived props
static int write_atom_property(molfile_atom_t &a,
const int propid,
const char *prop)
{
int plist = MolfileInterface::P_NONE;
PROPSTRCPY(MolfileInterface::P_NAME,name,prop);
PROPSTRCPY(MolfileInterface::P_TYPE,type,prop);
PROPSTRCPY(MolfileInterface::P_RESN,resname,prop);
PROPSTRCPY(MolfileInterface::P_SEGN,segid,prop);
return plist;
}
#undef PROPUPDATE
#undef PROPSTRCPY
// set/get atom floating point property
int MolfileInterface::property(int propid, int idx, float *prop)
{
if ((_info == NULL) || (prop == NULL) || (idx < 0) || (idx >= _natoms))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT)
_props |= write_atom_property(a[idx], propid, *prop);
if (_mode & M_RSTRUCT)
*prop = read_float_property(a[idx], propid);
return _props;
}
// set/get per type floating point property
int MolfileInterface::property(int propid, int *types, float *prop)
{
if ((_info == NULL) || (types == NULL) || (prop == NULL))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
for (int i=0; i < _natoms; ++i)
_props |= write_atom_property(a[i], propid, prop[types[i]]);
}
// useless for reading.
if (_mode & M_RSTRUCT)
return P_NONE;
return _props;
}
// set/get per atom floating point property
int MolfileInterface::property(int propid, float *prop)
{
if ((_info == NULL) || (prop == NULL))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
for (int i=0; i < _natoms; ++i)
_props |= write_atom_property(a[i], propid, prop[i]);
}
if (_mode & M_RSTRUCT) {
for (int i=0; i < _natoms; ++i)
prop[i] = read_float_property(a[i], propid);
}
return _props;
}
// set/get atom floating point property
int MolfileInterface::property(int propid, int idx, double *prop)
{
if ((_info == NULL) || (prop == NULL) || (idx < 0) || (idx >= _natoms))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT)
return write_atom_property(a[idx], propid, *prop);
if (_mode & M_RSTRUCT)
*prop = static_cast<double>(read_float_property(a[idx], propid));
return _props;
}
// set/get per type floating point property
int MolfileInterface::property(int propid, int *types, double *prop)
{
if ((_info == NULL) || (types == NULL) || (prop == NULL))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
for (int i=0; i < _natoms; ++i)
_props |= write_atom_property(a[i], propid, prop[types[i]]);
}
// useless for reading
if (_mode & M_RSTRUCT)
return P_NONE;
return _props;
}
// set/get per atom floating point property
int MolfileInterface::property(int propid, double *prop)
{
if ((_info == NULL) || (prop == NULL))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
for (int i=0; i < _natoms; ++i)
_props |= write_atom_property(a[i], propid, prop[i]);
}
if (_mode & M_RSTRUCT) {
for (int i=0; i < _natoms; ++i)
prop[i] = static_cast<double>(read_float_property(a[i], propid));
}
return _props;
}
#define INT_TO_STRING_BODY(IDX) \
buf[15] = 0; \
if (propid & P_NAME) \
_props |= write_atom_property(a[IDX],P_NAME,buf); \
if (propid & P_TYPE) \
_props |= write_atom_property(a[IDX],P_TYPE,buf); \
buf[7] = 0; \
if (propid & P_RESN) \
_props |= write_atom_property(a[IDX],P_RESN,buf); \
if (propid & P_SEGN) \
_props |= write_atom_property(a[IDX],P_SEGN,buf); \
buf[1] = 0; \
if (propid & P_CHAI) \
_props |= write_atom_property(a[IDX],P_CHAI,buf)
// set/get atom integer property
int MolfileInterface::property(int propid, int idx, int *prop)
{
if ((_info == NULL) || (prop == NULL) || (idx < 0) || (idx >= _natoms))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
char buf[64];
_props |= write_atom_property(a[idx], propid, *prop);
if (propid & (P_NAME|P_TYPE|P_RESN|P_SEGN|P_CHAI)) {
sprintf(buf,"%d",*prop);
INT_TO_STRING_BODY(idx);
}
}
if (_mode & M_RSTRUCT)
*prop = read_int_property(a[idx], propid);
return _props;
}
// set/get per type integer property
int MolfileInterface::property(int propid, int *types, int *prop)
{
if ((_info == NULL) || (types == NULL) || (prop == NULL))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
char buf[64];
for (int i=0; i < _natoms; ++i)
_props |= write_atom_property(a[i], propid, prop[types[i]]);
if (propid & (P_NAME|P_TYPE|P_RESN|P_SEGN|P_CHAI)) {
for (int i=0; i < _natoms; ++i) {
sprintf(buf,"%d",prop[types[i]]);
INT_TO_STRING_BODY(i);
}
}
}
// useless when reading
if (_mode & M_RSTRUCT)
return P_NONE;
return _props;
}
// set/get per atom integer property
int MolfileInterface::property(int propid, int *prop)
{
if ((_info == NULL) || (prop == NULL))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
char buf[64];
for (int i=0; i < _natoms; ++i)
_props |= write_atom_property(a[i],propid,prop[i]);
if (propid & (P_NAME|P_TYPE|P_RESN|P_SEGN|P_CHAI)) {
for (int i=0; i < _natoms; ++i) {
sprintf(buf,"%d",prop[i]);
INT_TO_STRING_BODY(i);
}
}
}
if (_mode & M_RSTRUCT) {
for (int i=0; i < _natoms; ++i)
prop[i] = read_int_property(a[i], propid);
}
return _props;
}
#undef INT_TO_STRING_BODY
// set/get atom string property
int MolfileInterface::property(int propid, int idx, char *prop)
{
if ((_info == NULL) || (prop == NULL) || (idx < 0) || (idx >= _natoms))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
_props |= write_atom_property(a[idx], propid, prop);
}
if (_mode & M_RSTRUCT)
strcpy(prop,read_string_property(a[idx], propid));
return _props;
}
// set/get per type string property
int MolfileInterface::property(int propid, int *types, char **prop)
{
if ((_info == NULL) || (types == NULL) || (prop == NULL))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
for (int i=0; i < _natoms; ++i) {
_props |= write_atom_property(a[i], propid, prop[types[i]]);
}
}
// useless when reading
if (_mode & M_RSTRUCT)
return P_NONE;
return _props;
}
// set/get per atom string property
int MolfileInterface::property(int propid, char **prop)
{
if ((_info == NULL) || (prop == NULL))
return P_NONE;
molfile_atom_t *a = static_cast<molfile_atom_t *>(_info);
if (_mode & M_WSTRUCT) {
for (int i=0; i < _natoms; ++i) {
_props |= write_atom_property(a[i], propid, prop[i]);
}
}
// not supported right now. XXX: should we use strdup() here?
if (_mode & M_RSTRUCT)
return P_NONE;
return _props;
}
Event Timeline
Log In to Comment