forked from lijiext/lammps
1240 lines
25 KiB
C
Executable File
1240 lines
25 KiB
C
Executable File
/* **************************************************************** */
|
|
/* functions to read data from files */
|
|
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <unistd.h>
|
|
#include <math.h>
|
|
|
|
/* commented this out for SGI gcc compiler
|
|
#ifdef INCL_FLOAT
|
|
double strtod(char *s, char **t);
|
|
long strtol(char *s, char **t, int base);
|
|
#endif
|
|
*/
|
|
|
|
#include <X11/Intrinsic.h>
|
|
|
|
extern FILE *popen(const char *, const char *);
|
|
extern int pclose(FILE *);
|
|
|
|
#include "xmovie.h"
|
|
|
|
#define LINELEN 256
|
|
#define WHITESPACE " \f\n\r\t\v"
|
|
#define SEPARATORS ",;:" WHITESPACE
|
|
#define NPRINT 500
|
|
#define NREAD 50
|
|
#define NBREAD 2000
|
|
|
|
#define ERROR -1
|
|
#define END_OF_FILE 1
|
|
#define BEING_NICE 2
|
|
#define READ_OK 3
|
|
#define LOOK_FOR_ITEM 10
|
|
#define LOOK_FOR_BOUNDS 11
|
|
#define LOOK_FOR_POSITIONS 12
|
|
#define LOOK_FOR_TIME 13
|
|
#define LOOK_FOR_BONDS 14
|
|
|
|
#define BAD_STATE -2
|
|
|
|
#define ALL_DONE TRUE
|
|
#define MORE_TO_READ FALSE
|
|
|
|
/* **************************************************************** */
|
|
/* local typedefs */
|
|
|
|
typedef struct {
|
|
char *string;
|
|
int state;
|
|
} PARSETABLE;
|
|
|
|
#define NPARSE(x) (sizeof(x)/sizeof(PARSETABLE))
|
|
|
|
typedef struct {
|
|
int nfiles;
|
|
char **fnames;
|
|
FILE *file;
|
|
int in_read;
|
|
int is_compressed;
|
|
int is_binary;
|
|
} READDATA;
|
|
|
|
typedef int (*PFI)();
|
|
|
|
typedef struct {
|
|
RECORD record_id;
|
|
PFI reader;
|
|
} READER;
|
|
|
|
/* **************************************************************** */
|
|
/* local prototypes */
|
|
|
|
PRIVATE int DoSomeReading(void);
|
|
PRIVATE int ReadSomeBinary(void);
|
|
PRIVATE int parser(char *token, PARSETABLE *table, int cnt);
|
|
PRIVATE int get_bounds(FILE *file, char *line);
|
|
PRIVATE int get_positions(FILE *file, char *line);
|
|
PRIVATE int get_time(FILE *file, char *line);
|
|
PRIVATE int get_bonds(FILE *file, char *line);
|
|
PRIVATE int copy_bonds(void);
|
|
PRIVATE void add_point(int index, int type, float x, float y, float z);
|
|
PRIVATE void add_bond(int type, int index1, int index2);
|
|
PRIVATE void NewData(void);
|
|
PRIVATE void ShowStatus(void);
|
|
PRIVATE char *RemoveNL(char *s);
|
|
PRIVATE int getline(FILE *file, char *line);
|
|
PRIVATE int get_float(char **s, float *f);
|
|
PRIVATE int get_int(char **s, int *i);
|
|
PRIVATE void ungetline(void);
|
|
PRIVATE int LineIsBlank(char *s);
|
|
PRIVATE void SortData(DATA *dptr);
|
|
PRIVATE int compare_atypes(POSITION *a, POSITION *b);
|
|
PRIVATE int compare_btypes(BOND *a, BOND *b);
|
|
PRIVATE int IsCompressed(char *s);
|
|
PRIVATE int IsBinary(char *s);
|
|
PRIVATE int TestBinary(FILE *f, int rewind);
|
|
PRIVATE FILE *GetCompressed(char *s);
|
|
PRIVATE FILE *OpenFile(char *s, int compressed, int binary);
|
|
PRIVATE READER *GetReader(RECORD record_id);
|
|
PRIVATE int ReadTime(FILE *f, LENGTH length);
|
|
PRIVATE int ReadBound(FILE *f, LENGTH length);
|
|
PRIVATE int ReadPosition(FILE *f, LENGTH length);
|
|
PRIVATE int ReadBond(FILE *f, LENGTH length);
|
|
PRIVATE int CopyBond(FILE *f, LENGTH length);
|
|
PRIVATE int ReadRecordHeader(FILE *f, RECORD *record, LENGTH *length);
|
|
|
|
/* **************************************************************** */
|
|
/* local data */
|
|
|
|
/* include last 3 out-of-date item names in PARSETABLE for compatability */
|
|
|
|
static PARSETABLE parse_table[] = {
|
|
{ "TIMESTEP", LOOK_FOR_TIME },
|
|
{ "BOX BOUNDS", LOOK_FOR_BOUNDS },
|
|
{ "ATOMS", LOOK_FOR_POSITIONS },
|
|
{ "BONDS", LOOK_FOR_BONDS },
|
|
{ "TIME", LOOK_FOR_TIME },
|
|
{ "BOUNDS", LOOK_FOR_BOUNDS },
|
|
{ "POSITIONS", LOOK_FOR_POSITIONS },
|
|
};
|
|
|
|
static DATA *current_data = (DATA *) NULL;
|
|
static int seen_time = 0;
|
|
static int npositions = 0;
|
|
static int nlines = 0;
|
|
static int nbonds = 0;
|
|
static int read_this_time = 0;
|
|
|
|
static READDATA ReadData = {
|
|
0, /* nfiles */
|
|
NULL, /* filenames */
|
|
NULL, /* file pointer */
|
|
0, /* not in read */
|
|
0, /* not compressed */
|
|
0, /* not binary */
|
|
};
|
|
|
|
static READER read_table[] = {
|
|
{ TimeID, ReadTime },
|
|
{ BoundID, ReadBound },
|
|
{ PositionID, ReadPosition },
|
|
{ BondID, ReadBond },
|
|
{ CopyBondID, CopyBond },
|
|
};
|
|
|
|
/* **************************************************************** */
|
|
|
|
void InitRead(int nfiles, char **fnames)
|
|
{
|
|
static char *Stdin = "stdin";
|
|
|
|
if (nfiles > 0) {
|
|
ReadData.nfiles = nfiles;
|
|
ReadData.fnames = fnames;
|
|
return;
|
|
}
|
|
|
|
if (isatty(fileno(stdin))) {
|
|
Usage();
|
|
exit(0);
|
|
}
|
|
|
|
ReadData.nfiles = 1;
|
|
ReadData.fnames = &Stdin;
|
|
ReadData.file = stdin;
|
|
ReadData.in_read = 1;
|
|
ReadData.is_binary = TestBinary(stdin, TRUE);
|
|
|
|
if (ReadData.is_binary)
|
|
fseek(stdin, 2*sizeof(INT4), SEEK_CUR); /* skip header */
|
|
}
|
|
|
|
Boolean ReadProc(XtPointer client_data)
|
|
{
|
|
char err_msg[512];
|
|
char *file_type;
|
|
|
|
if (ReadData.in_read) goto do_read;
|
|
|
|
next_file:
|
|
|
|
if (ReadData.nfiles <= 0) {
|
|
ShowStatus();
|
|
SortData(current_data);
|
|
return(ALL_DONE);
|
|
}
|
|
|
|
ReadData.is_compressed = FALSE;
|
|
ReadData.is_binary = FALSE;
|
|
|
|
ReadData.is_compressed = IsCompressed(*ReadData.fnames);
|
|
|
|
if (ReadData.is_compressed) {
|
|
file_type = "compressed";
|
|
goto open_file;
|
|
}
|
|
|
|
ReadData.is_binary = IsBinary(*ReadData.fnames);
|
|
if (ReadData.is_binary) {
|
|
file_type = "binary";
|
|
goto open_file;
|
|
}
|
|
|
|
|
|
file_type = "ascii";
|
|
goto open_file;
|
|
|
|
open_file:
|
|
|
|
ReadData.file = OpenFile(*ReadData.fnames,
|
|
ReadData.is_compressed, ReadData.is_binary);
|
|
if (ReadData.file == (FILE *) NULL) {
|
|
sprintf(err_msg,
|
|
"Frames: %i Positions: %i Bonds: %i\n"
|
|
"Error: unable to open %s file \"%s\".\n"
|
|
"Read aborted.",
|
|
Common.ndata, npositions, nbonds,
|
|
file_type, *ReadData.fnames);
|
|
SetReadString(err_msg);
|
|
SortData(current_data);
|
|
return(ALL_DONE);
|
|
}
|
|
|
|
ReadData.in_read = 1;
|
|
nlines = 0;
|
|
ShowStatus();
|
|
|
|
do_read:
|
|
|
|
switch(DoSomeReading()) {
|
|
case END_OF_FILE:
|
|
if (ReadData.is_compressed)
|
|
pclose(ReadData.file);
|
|
else
|
|
if (strcmp("stdin",*ReadData.fnames))
|
|
fclose(ReadData.file);
|
|
ReadData.in_read = 0;
|
|
ReadData.nfiles--;
|
|
ReadData.fnames++;
|
|
goto next_file;
|
|
case ERROR:
|
|
SortData(current_data);
|
|
return(ALL_DONE);
|
|
case BEING_NICE:
|
|
break;
|
|
default:
|
|
fprintf(stderr,"Don't know how I got here: ReadProc.\n");
|
|
exit(0);
|
|
}
|
|
|
|
return(MORE_TO_READ);
|
|
}
|
|
|
|
PRIVATE int IsCompressed(char *s)
|
|
{
|
|
register char *t;
|
|
|
|
t = s + strlen(s) - 3;
|
|
if (t < s) return(FALSE);
|
|
return(!strcmp(t,".gz"));
|
|
}
|
|
|
|
PRIVATE int IsBinary(char *s)
|
|
{
|
|
FILE *f;
|
|
int ok;
|
|
|
|
if ((FILE *) NULL == (f = fopen(s, "r"))) return(FALSE);
|
|
|
|
ok = TestBinary(f, FALSE);
|
|
fclose(f);
|
|
|
|
return(ok);
|
|
}
|
|
|
|
PRIVATE int TestBinary(FILE *f, int rewind)
|
|
{
|
|
union {
|
|
INT4 magic;
|
|
char c[sizeof(INT4)];
|
|
} u;
|
|
int i;
|
|
|
|
if (1 != fread(&u.magic, sizeof(u.magic), 1, f)) return(FALSE);
|
|
|
|
if (rewind)
|
|
for(i = sizeof(u); i; i--)
|
|
ungetc(u.c[i-1], f);
|
|
|
|
return(u.magic == MagicNumber);
|
|
}
|
|
|
|
PRIVATE FILE *GetCompressed(char *name)
|
|
{
|
|
char cmd[256];
|
|
|
|
sprintf(cmd, "gunzip -c %s", name);
|
|
return(popen(cmd,"r"));
|
|
}
|
|
|
|
PRIVATE FILE *OpenFile(char *name, int compressed, int binary)
|
|
{
|
|
FILE *f;
|
|
|
|
if (compressed) return(GetCompressed(name));
|
|
|
|
f = fopen(name, "r");
|
|
if (f == (FILE *) NULL) return(f);
|
|
|
|
if (binary) fseek(f, 2*sizeof(INT4), SEEK_CUR);
|
|
|
|
return(f);
|
|
}
|
|
|
|
PRIVATE int DoSomeReading(void)
|
|
{
|
|
static int state = LOOK_FOR_ITEM;
|
|
static char line[LINELEN];
|
|
static char token[LINELEN];
|
|
static char *t;
|
|
int ntries;
|
|
int newstate;
|
|
char err_msg[512];
|
|
|
|
if (ReadData.is_binary) return(ReadSomeBinary());
|
|
|
|
read_this_time = 0;
|
|
|
|
begin:
|
|
|
|
switch(state) {
|
|
case LOOK_FOR_ITEM:
|
|
ntries = 0;
|
|
while(1) {
|
|
ntries++;
|
|
if (ntries > 100) return(BEING_NICE);
|
|
|
|
if (!getline(ReadData.file, line))
|
|
return(END_OF_FILE);
|
|
t = strtok(line, WHITESPACE);
|
|
if (!t) continue;
|
|
if (strcmp(t, "ITEM:")) continue;
|
|
|
|
t = strtok(NULL, WHITESPACE);
|
|
if (!t) continue;
|
|
sprintf(token,"%s",t);
|
|
|
|
t = strtok(NULL, WHITESPACE);
|
|
while (t != NULL) {
|
|
sprintf(token,"%s %s",token,t);
|
|
t = strtok(NULL, WHITESPACE);
|
|
}
|
|
|
|
newstate = parser(token,
|
|
parse_table,NPARSE(parse_table));
|
|
|
|
if (newstate == BAD_STATE) continue;
|
|
state = newstate;
|
|
goto begin;
|
|
}
|
|
case LOOK_FOR_TIME:
|
|
if (!get_time(ReadData.file, line)) goto print_error;
|
|
state = LOOK_FOR_ITEM;
|
|
goto begin;
|
|
case LOOK_FOR_BOUNDS:
|
|
if (!get_bounds(ReadData.file, line)) goto print_error;
|
|
state = LOOK_FOR_ITEM;
|
|
goto begin;
|
|
case LOOK_FOR_POSITIONS:
|
|
case LOOK_FOR_BONDS:
|
|
switch((state == LOOK_FOR_POSITIONS) ?
|
|
get_positions(ReadData.file, line) :
|
|
get_bonds(ReadData.file, line)) {
|
|
case END_OF_FILE:
|
|
state = LOOK_FOR_ITEM;
|
|
return(END_OF_FILE);
|
|
case BEING_NICE:
|
|
return(BEING_NICE);
|
|
case LOOK_FOR_ITEM:
|
|
state = LOOK_FOR_ITEM;
|
|
return(BEING_NICE);
|
|
case ERROR:
|
|
default:
|
|
goto print_error;
|
|
}
|
|
|
|
case BAD_STATE:
|
|
default:
|
|
fprintf(stderr,"Error: read state corrupted: %i\n", state);
|
|
exit(0);
|
|
}
|
|
|
|
print_error:
|
|
|
|
sprintf(err_msg,"Frames: %i Positions: %i Bonds: %i\n"
|
|
"Error occurred at line %i\n"
|
|
"Processing ITEM: %s\n"
|
|
"Line = \"%s\"\n"
|
|
"Read aborted.",
|
|
Common.ndata, npositions, nbonds, nlines, token,
|
|
RemoveNL(line));
|
|
|
|
SetReadString(err_msg);
|
|
|
|
return(ERROR);
|
|
|
|
/* exit(0); */
|
|
|
|
}
|
|
|
|
/* **************************************************************** */
|
|
|
|
PRIVATE int parser(char *token, PARSETABLE *table, int cnt)
|
|
{
|
|
PARSETABLE *t;
|
|
int i;
|
|
|
|
for(i = cnt, t = table; i; i--, t++)
|
|
if (strstr(token,t->string) == token)
|
|
return t->state;
|
|
return(BAD_STATE);
|
|
}
|
|
|
|
/* **************************************************************** */
|
|
|
|
PRIVATE int get_bounds(FILE *file, char *line)
|
|
{
|
|
int i, n, changed;
|
|
char *t;
|
|
BOUND *b, *cb;
|
|
|
|
if (!current_data) NewData();
|
|
|
|
n = (Common.two_d) ? 2 : 3;
|
|
|
|
for(i = 0; i < n; i++){
|
|
if (!getline(file, line)) return(0);
|
|
|
|
t = line;
|
|
|
|
if (!get_float(&t, &(current_data->bounds[i].low)))
|
|
return(0);
|
|
|
|
if (!get_float(&t, &(current_data->bounds[i].high)))
|
|
return(0);
|
|
}
|
|
|
|
if (Common.two_d) {
|
|
current_data->bounds[2].low = -0.01;
|
|
current_data->bounds[2].high = 0.01;
|
|
}
|
|
|
|
changed = 0;
|
|
for(b=current_data->bounds,cb=Common.bounds,i=3; i;i--,b++,cb++){
|
|
if (b->low < cb->low) {
|
|
cb->low = b->low;
|
|
changed = 1;
|
|
}
|
|
if (b->high > cb->high) {
|
|
cb->high = b->high;
|
|
changed = 1;
|
|
}
|
|
}
|
|
|
|
if (changed) NewDataSetup();
|
|
|
|
return(1);
|
|
}
|
|
|
|
PRIVATE int get_positions(FILE *file, char *line)
|
|
{
|
|
char *t;
|
|
float x, y, z;
|
|
float xlow,xhigh,ylow,yhigh,zlow,zhigh,dx,dy,dz;
|
|
int index, type;
|
|
|
|
if (!current_data) NewData();
|
|
|
|
if (Common.copy_bond) copy_bonds();
|
|
|
|
if (Common.remap || Common.scaleflag) {
|
|
xlow = current_data->bounds[0].low;
|
|
xhigh = current_data->bounds[0].high;
|
|
dx = xhigh - xlow;
|
|
ylow = current_data->bounds[1].low;
|
|
yhigh = current_data->bounds[1].high;
|
|
dy = yhigh - ylow;
|
|
if (!Common.two_d) {
|
|
zlow = current_data->bounds[2].low;
|
|
zhigh = current_data->bounds[2].high;
|
|
dz = zhigh - zlow;
|
|
}
|
|
}
|
|
|
|
z = 0.0;
|
|
|
|
while(1) {
|
|
|
|
if (read_this_time++ > NREAD) return(BEING_NICE);
|
|
|
|
if (!getline(file, line)) return(END_OF_FILE);
|
|
|
|
t = line;
|
|
|
|
if (!get_int(&t, &index)){
|
|
ungetline();
|
|
return(LOOK_FOR_ITEM);
|
|
}
|
|
|
|
if (!get_int(&t, &type)) return(ERROR);
|
|
|
|
if (!get_float(&t, &x)) return(ERROR);
|
|
if (!get_float(&t, &y)) return(ERROR);
|
|
|
|
if (!Common.two_d)
|
|
if (!get_float(&t, &z)) return(ERROR);
|
|
|
|
if (type < 1 || type > Common.natomcolors) type = 1;
|
|
|
|
if (Common.remap) {
|
|
if (x-xlow >= 0.0)
|
|
x = fmod(x-xlow,dx) + xlow;
|
|
else
|
|
x = xhigh + fmod(x-xlow,dx);
|
|
if (y-ylow >= 0.0)
|
|
y = fmod(y-ylow,dy) + ylow;
|
|
else
|
|
y = yhigh + fmod(y-ylow,dy);
|
|
if (!Common.two_d) {
|
|
if (z-zlow >= 0.0)
|
|
z = fmod(z-zlow,dz) + zlow;
|
|
else
|
|
z = zhigh + fmod(z-zlow,dz);
|
|
}
|
|
}
|
|
|
|
if (Common.scaleflag) {
|
|
x = xlow + x*dx;
|
|
y = ylow + y*dy;
|
|
if (!Common.two_d) z = zlow + z*dz;
|
|
}
|
|
|
|
add_point(index, type, x, y, z);
|
|
|
|
npositions++;
|
|
|
|
if ((nbonds + npositions) % NPRINT == 0) ShowStatus();
|
|
|
|
}
|
|
|
|
return(1);
|
|
}
|
|
|
|
PRIVATE int get_bonds(FILE *file, char *line)
|
|
{
|
|
char *t;
|
|
int type, index1, index2;
|
|
|
|
if (!current_data) NewData();
|
|
|
|
while(1) {
|
|
|
|
if (read_this_time++ > NREAD) return(BEING_NICE);
|
|
|
|
if (!getline(file, line)) return(END_OF_FILE);
|
|
|
|
t = line;
|
|
|
|
if (!get_int(&t, &type)){
|
|
ungetline();
|
|
return(LOOK_FOR_ITEM);
|
|
}
|
|
|
|
if (!get_int(&t, &index1)) return(ERROR);
|
|
if (!get_int(&t, &index2)) return(ERROR);
|
|
|
|
if (type < 1 || type > Common.nbondcolors) type = 1;
|
|
|
|
add_bond(type, index1, index2);
|
|
|
|
nbonds++;
|
|
|
|
if ((nbonds + npositions) % NPRINT == 0) ShowStatus();
|
|
}
|
|
|
|
return(1);
|
|
}
|
|
|
|
PRIVATE int get_time(FILE *file, char *line)
|
|
{
|
|
char *t;
|
|
|
|
if (!current_data || seen_time) NewData();
|
|
seen_time = 1;
|
|
|
|
if (!getline(file, line)) return(0);
|
|
|
|
t = line;
|
|
|
|
if (!get_float(&t, &(current_data->time))) return(0);
|
|
|
|
return(1);
|
|
}
|
|
|
|
PRIVATE int copy_bonds()
|
|
{
|
|
DATA *last_data;
|
|
int i;
|
|
BOND *b1, *b2;
|
|
|
|
if (Common.ndata < 2) return(0);
|
|
|
|
last_data = current_data-1;
|
|
|
|
current_data->nbonds = last_data->nbonds;
|
|
current_data->maxbonds = last_data->maxbonds;
|
|
|
|
current_data->bonds =
|
|
(BOND *) XtMalloc(current_data->maxbonds * sizeof(BOND));
|
|
|
|
b1 = current_data->bonds;
|
|
b2 = last_data->bonds;
|
|
for(i = current_data->nbonds; i ; i--, b1++, b2++)
|
|
*b1 = *b2;
|
|
|
|
return(1);
|
|
}
|
|
|
|
/* **************************************************************** */
|
|
|
|
|
|
PRIVATE char *RemoveNL(char *s)
|
|
{
|
|
register char *t;
|
|
|
|
if (!s) return(s);
|
|
|
|
if ((t = s + strlen(s) - 1) < s) return(s);
|
|
if (*t == '\n') *t = '\0';
|
|
|
|
return(s);
|
|
}
|
|
|
|
static int got_line = 0;
|
|
|
|
PRIVATE int getline(FILE *file, char *line)
|
|
{
|
|
int result;
|
|
|
|
if (got_line) {
|
|
got_line = 0;
|
|
return(1);
|
|
}
|
|
|
|
do {
|
|
nlines++;
|
|
result = (fgets(line, LINELEN, file) != NULL);
|
|
} while( result && LineIsBlank(line) );
|
|
|
|
return(result);
|
|
}
|
|
|
|
PRIVATE void ungetline(void)
|
|
{
|
|
got_line = 1;
|
|
}
|
|
|
|
PRIVATE int LineIsBlank(char *s)
|
|
{
|
|
register char *t;
|
|
|
|
t = s + strspn(s, WHITESPACE);
|
|
return( *t == '\0' );
|
|
}
|
|
|
|
PRIVATE int get_int(char **s, int *i)
|
|
{
|
|
register char *t;
|
|
|
|
t = *s + strspn(*s, SEPARATORS);
|
|
*i = strtol(t, s, 10);
|
|
|
|
return(t != *s);
|
|
}
|
|
|
|
PRIVATE int get_float(char **s, float *f)
|
|
{
|
|
register char *t;
|
|
|
|
t = *s + strspn(*s, SEPARATORS);
|
|
*f = strtod(t, s);
|
|
|
|
return(t != *s);
|
|
}
|
|
|
|
/* **************************************************************** */
|
|
|
|
PRIVATE void NewData(void)
|
|
{
|
|
int i;
|
|
DATA *last_data,*olddataptr;
|
|
|
|
olddataptr = Common.dataptr;
|
|
last_data = current_data;
|
|
|
|
if (Common.ndata >= Common.maxdata) { /* need to realloc */
|
|
Common.maxdata += 256;
|
|
Common.dataptr = (DATA *)
|
|
XtMalloc(Common.maxdata * sizeof(DATA));
|
|
memcpy(Common.dataptr,olddataptr,Common.ndata*sizeof(DATA));
|
|
}
|
|
|
|
current_data = Common.dataptr + Common.ndata;
|
|
Common.ndata++;
|
|
|
|
current_data->natypes =
|
|
(INT4 *) XtMalloc(Common.natomcolors * sizeof(INT4));
|
|
|
|
for(i = 0; i < Common.natomcolors; i++)
|
|
current_data->natypes[i] = 0;
|
|
|
|
current_data->nbtypes =
|
|
(INT4 *) XtMalloc(Common.nbondcolors * sizeof(INT4));
|
|
|
|
for(i = 0; i < Common.nbondcolors; i++)
|
|
current_data->nbtypes[i] = 0;
|
|
|
|
if (Common.ndata > 1) {
|
|
current_data->time = last_data->time;
|
|
|
|
for(i = 0; i < 3; i++)
|
|
current_data->bounds[i] = last_data->bounds[i];
|
|
|
|
current_data->maxatoms = last_data->natoms;
|
|
current_data->positions = (POSITION *)
|
|
XtMalloc(last_data->natoms * sizeof(POSITION));
|
|
|
|
current_data->natoms = 0;
|
|
|
|
current_data->maxbonds = last_data->nbonds;
|
|
current_data->bonds = (BOND *)
|
|
XtMalloc(last_data->nbonds * sizeof(BOND));
|
|
|
|
current_data->nbonds = 0;
|
|
|
|
SortData(last_data);
|
|
}
|
|
else {
|
|
current_data->time = 0.0;
|
|
|
|
for(i = 0; i < 3; i++) {
|
|
current_data->bounds[i].low = -1.0;
|
|
current_data->bounds[i].high = 1.0;
|
|
}
|
|
|
|
current_data->natoms = 0;
|
|
current_data->positions = NULL;
|
|
current_data->maxatoms = 0;
|
|
|
|
current_data->nbonds = 0;
|
|
current_data->bonds = NULL;
|
|
current_data->maxbonds = 0;
|
|
}
|
|
|
|
if (olddataptr != Common.dataptr) XtFree((char *) olddataptr);
|
|
}
|
|
|
|
PRIVATE void add_point(int index, int type, float x, float y, float z)
|
|
{
|
|
POSITION *p;
|
|
|
|
if (current_data->natoms >= current_data->maxatoms) {
|
|
current_data->maxatoms += 16;
|
|
current_data->positions = (POSITION *)
|
|
XtRealloc((char *) current_data->positions,
|
|
current_data->maxatoms * sizeof(POSITION));
|
|
}
|
|
|
|
p = &(current_data->positions[current_data->natoms]);
|
|
|
|
p->index = index;
|
|
p->type = type;
|
|
p->x = x;
|
|
p->y = y;
|
|
p->z = z;
|
|
|
|
current_data->natoms++;
|
|
}
|
|
|
|
|
|
PRIVATE void add_bond(int type, int index1, int index2)
|
|
{
|
|
BOND *b;
|
|
|
|
if (current_data->nbonds >= current_data->maxbonds) {
|
|
current_data->maxbonds += 16;
|
|
current_data->bonds = (BOND *)
|
|
XtRealloc((char *) current_data->bonds,
|
|
current_data->maxbonds * sizeof(BOND));
|
|
}
|
|
|
|
b = &(current_data->bonds[current_data->nbonds]);
|
|
|
|
b->type = type;
|
|
b->index1 = index1;
|
|
b->index2 = index2;
|
|
b->atom1 = NULL;
|
|
b->atom2 = NULL;
|
|
|
|
current_data->nbonds++;
|
|
}
|
|
|
|
PRIVATE void ShowStatus(void)
|
|
{
|
|
static char buffer[256];
|
|
|
|
if (ReadData.in_read)
|
|
sprintf(buffer, "Reading file: %s\n"
|
|
"Frames: %i Positions: %i Bonds: %i",
|
|
*ReadData.fnames, Common.ndata, npositions,
|
|
nbonds);
|
|
else
|
|
sprintf(buffer, "Reading done.\n"
|
|
"Frames: %i Positions: %i Bonds %i",
|
|
Common.ndata, npositions, nbonds);
|
|
|
|
SetReadString(buffer);
|
|
}
|
|
|
|
/* **************************************************************** */
|
|
|
|
PRIVATE void SortData(DATA *dptr)
|
|
{
|
|
POSITION *p, **p_by_index, **pp;
|
|
register int i, j;
|
|
int maxindex, minindex, nindex, n;
|
|
register BOND *b, *b2;
|
|
float bond_length, shortest, dx, dy, dz;
|
|
|
|
/* add hpsort in place of standard C qsort, which is very slow
|
|
when atoms or bonds are already in order (N^2 instead of NlnN)
|
|
calling syntax does not change except for first argument
|
|
becomes arg-1, qsort -> hpsort in 2 places
|
|
*/
|
|
void hpsort(POSITION *, int, int, int (*)(POSITION *, POSITION *));
|
|
|
|
if (!dptr) return;
|
|
|
|
/* sort by type */
|
|
|
|
hpsort((dptr->positions)-1, dptr->natoms, sizeof(POSITION),
|
|
compare_atypes);
|
|
|
|
/* find out how many of each type */
|
|
|
|
for(i = dptr->natoms, p = dptr->positions; i ;i--, p++)
|
|
dptr->natypes[p->type-1]++;
|
|
|
|
if (dptr->nbonds < 1) return;
|
|
|
|
/* find max, min indicies */
|
|
|
|
minindex = 1000000000;
|
|
maxindex = -minindex;
|
|
|
|
for(i = dptr->natoms, p = dptr->positions; i; i--, p++) {
|
|
if (p->index > maxindex) maxindex = p->index;
|
|
if (p->index < minindex) minindex = p->index;
|
|
}
|
|
|
|
nindex = maxindex - minindex + 1;
|
|
|
|
if (nindex < 0) {
|
|
dptr->nbonds = 0;
|
|
return;
|
|
}
|
|
|
|
|
|
if (nindex < 1000000) {
|
|
|
|
p_by_index = (POSITION **) XtMalloc(nindex*sizeof(POSITION *));
|
|
|
|
for(pp = p_by_index, i = nindex; i; i--, pp++)
|
|
*pp = (POSITION *) NULL;
|
|
|
|
for(p = dptr->positions, i = dptr->natoms; i; i--, p++){
|
|
pp = p_by_index + p->index - minindex;
|
|
if (*pp != (POSITION *) NULL) {
|
|
fprintf(stderr,
|
|
"Error: atom %i has more than one position at time %g.\n",
|
|
p->index, dptr->time);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
*pp = p;
|
|
}
|
|
|
|
for(b = dptr->bonds, i = dptr->nbonds; i ; i--, b++) {
|
|
b->atom1 = p_by_index[b->index1 - minindex];
|
|
b->atom2 = p_by_index[b->index2 - minindex];
|
|
if (b->atom1 == (POSITION *) NULL) {
|
|
fprintf(stderr,
|
|
"Error: atom %d has no position at time %g.\n",
|
|
b->index1, dptr->time);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
if (b->atom2 == (POSITION *) NULL) {
|
|
fprintf(stderr,
|
|
"Error: atom %d has no position at time %g.\n",
|
|
b->index2, dptr->time);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
}
|
|
|
|
XtFree((char *) p_by_index);
|
|
}
|
|
else {
|
|
for(b = dptr->bonds, i = dptr->nbonds; i; i--, b++) {
|
|
for(p = dptr->positions, j = dptr->natoms; j; j--, p++)
|
|
if (b->index1 == p->index) {
|
|
b->atom1 = p;
|
|
break;
|
|
}
|
|
if (j < 1) {
|
|
fprintf(stderr,
|
|
"Error: atom %i has no position at time %g.\n",
|
|
b->index1, dptr->time);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
for(p = dptr->positions, j = dptr->natoms; j; j--, p++)
|
|
if (b->index2 == p->index) {
|
|
b->atom2 = p;
|
|
break;
|
|
}
|
|
if (j < 1) {
|
|
fprintf(stderr,
|
|
"Error: atom %i has no position at time %g.\n",
|
|
b->index2, dptr->time);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* throw out bonds longer than half the shortest cell size */
|
|
/* only do if periodic-bond switch is set */
|
|
|
|
if (Common.pbc_bond) {
|
|
|
|
n = (Common.two_d) ? 2 : 3;
|
|
|
|
shortest = 1e30;
|
|
for(i = 0; i < n; i++){
|
|
bond_length = dptr->bounds[i].high - dptr->bounds[i].low;
|
|
if (bond_length < shortest) shortest = bond_length;
|
|
}
|
|
|
|
shortest *= 0.5;
|
|
shortest *= shortest;
|
|
|
|
b = dptr->bonds;
|
|
b2 = b;
|
|
for(i = dptr->nbonds; i; i--, b++) {
|
|
dx = b->atom1->x - b->atom2->x;
|
|
dy = b->atom1->y - b->atom2->y;
|
|
dz = b->atom1->z - b->atom2->z;
|
|
bond_length = dx*dx + dy*dy + dz*dz;
|
|
if (bond_length < shortest) {
|
|
if (b2 != b) *b2 = *b;
|
|
b2++;
|
|
}
|
|
}
|
|
|
|
dptr->nbonds -= b - b2;
|
|
|
|
}
|
|
|
|
/* sort by type */
|
|
|
|
hpsort((dptr->bonds)-1, dptr->nbonds, sizeof(BOND),
|
|
compare_btypes);
|
|
|
|
/* find out how many of each type */
|
|
|
|
for(b = dptr->bonds, i = dptr->nbonds; i; i--, b++)
|
|
dptr->nbtypes[b->type-1]++;
|
|
}
|
|
|
|
PRIVATE int compare_atypes(POSITION *a, POSITION *b)
|
|
{
|
|
if (a->type > b->type) return(1);
|
|
if (a->type == b->type) return(0);
|
|
if (a->type < b->type) return(-1);
|
|
}
|
|
|
|
PRIVATE int compare_btypes(BOND *a, BOND *b)
|
|
{
|
|
if (a->type > b->type) return(1);
|
|
if (a->type == b->type) return(0);
|
|
if (a->type < b->type) return(-1);
|
|
}
|
|
|
|
/* routine to print out all data for debugging */
|
|
|
|
#if 0
|
|
|
|
print_them(void)
|
|
{
|
|
int i;
|
|
|
|
for(i = 0; i < Common.ndata; i++)
|
|
print_data(Common.dataptr[i]);
|
|
}
|
|
|
|
print_data(DATA *data)
|
|
{
|
|
int i;
|
|
POSITION *p;
|
|
|
|
fprintf(stderr,"Time: %f\n", data->time);
|
|
|
|
for(i = 0; i < 3; i++)
|
|
fprintf(stderr,"Bounds %i: %f %f\n", i, data->bounds[i].low,
|
|
data->bounds[i].high);
|
|
|
|
for(i = 0; i < data->natoms; i++) {
|
|
p = &(data->positions[i]);
|
|
fprintf(stderr,"Atom: %i %i %f %f %f\n",
|
|
p->index, p->type, p->x, p->y, p->z);
|
|
}
|
|
}
|
|
|
|
#endif
|
|
|
|
/* **************************************************************** */
|
|
|
|
PRIVATE int ReadSomeBinary(void)
|
|
{
|
|
RECORD record_id;
|
|
LENGTH length;
|
|
READER *r;
|
|
char err_msg[256];
|
|
int state;
|
|
int old_positions, old_bonds;
|
|
|
|
old_positions = npositions;
|
|
old_bonds = nbonds;
|
|
|
|
while(1) {
|
|
if (!ReadRecordHeader(ReadData.file, &record_id, &length))
|
|
return(END_OF_FILE);
|
|
|
|
r = GetReader(record_id);
|
|
if (r == (READER *) NULL) { /* dont recognize */
|
|
fseek(ReadData.file, length, SEEK_CUR);
|
|
continue;
|
|
}
|
|
|
|
switch(state = (*r->reader)(ReadData.file, length)) {
|
|
case BEING_NICE:
|
|
if (old_positions + old_bonds + NBREAD >
|
|
nbonds + npositions)
|
|
continue;
|
|
return(state);
|
|
case END_OF_FILE:
|
|
return(state);
|
|
case READ_OK:
|
|
continue;
|
|
case ERROR:
|
|
default:
|
|
sprintf(err_msg,"Frames: %i Positions: %i Bonds: %i\n"
|
|
"Error occurred during binary read.\n"
|
|
"Read aborted.",
|
|
Common.ndata, npositions, nbonds);
|
|
SetReadString(err_msg);
|
|
return(ERROR);
|
|
}
|
|
}
|
|
}
|
|
|
|
PRIVATE READER *GetReader(RECORD record_id)
|
|
{
|
|
register int i;
|
|
register READER *r;
|
|
|
|
i = sizeof(read_table)/sizeof(read_table[0]);
|
|
for(r = read_table; i; i--, r++)
|
|
if (r->record_id == record_id) return(r);
|
|
|
|
return((READER *) NULL);
|
|
}
|
|
|
|
PRIVATE int ReadRecordHeader(FILE *f, RECORD *record, LENGTH *length)
|
|
{
|
|
if (1 != fread(record, sizeof(*record), 1, f)) return(FALSE);
|
|
if (1 != fread(length, sizeof(*length), 1, f)) return(FALSE);
|
|
|
|
return(TRUE);
|
|
}
|
|
|
|
/*ARGSUSED*/
|
|
PRIVATE int ReadBound(FILE *file, LENGTH length)
|
|
{
|
|
int i, changed;
|
|
BOUND *b, *cb;
|
|
|
|
if (!current_data) NewData();
|
|
|
|
|
|
|
|
for(b = current_data->bounds, i = 3; i; i--, b++){
|
|
if (1 != fread(&(b->low), sizeof(b->low), 1, file))
|
|
return(ERROR);
|
|
if (1 != fread(&(b->high), sizeof(b->high), 1, file))
|
|
return(ERROR);
|
|
}
|
|
|
|
changed = 0;
|
|
for(b=current_data->bounds,cb=Common.bounds,i=3; i;i--,b++,cb++){
|
|
if (b->low < cb->low) {
|
|
cb->low = b->low;
|
|
changed = 1;
|
|
}
|
|
if (b->high > cb->high) {
|
|
cb->high = b->high;
|
|
changed = 1;
|
|
}
|
|
}
|
|
|
|
if (changed) NewDataSetup();
|
|
|
|
return(READ_OK);
|
|
}
|
|
|
|
PRIVATE int ReadPosition(FILE *file, LENGTH length)
|
|
{
|
|
POSITION *p;
|
|
int n;
|
|
int i;
|
|
|
|
if (!current_data) NewData();
|
|
|
|
n = length / (sizeof(p->index) + sizeof(p->type) + 3*sizeof(p->x));
|
|
|
|
if (current_data->maxatoms < n) {
|
|
free(current_data->positions);
|
|
current_data->maxatoms = n;
|
|
current_data->positions =
|
|
(POSITION *) XtMalloc(n * sizeof(POSITION));
|
|
}
|
|
|
|
current_data->natoms = 0;
|
|
p = current_data->positions;
|
|
i = 0;
|
|
|
|
for(; i < n; i++, p++) {
|
|
if (1 != fread(&(p->index), sizeof(p->index), 1, file))
|
|
return(ERROR);
|
|
if (1 != fread(&(p->type), sizeof(p->type), 1, file))
|
|
return(ERROR);
|
|
if (1 != fread(&(p->x), sizeof(p->x), 1, file))
|
|
return(ERROR);
|
|
if (1 != fread(&(p->y), sizeof(p->y), 1, file))
|
|
return(ERROR);
|
|
if (1 != fread(&(p->z), sizeof(p->z), 1, file))
|
|
return(ERROR);
|
|
|
|
if (p->type < 1 || p->type > Common.natomcolors)
|
|
p->type = 1;
|
|
|
|
npositions++;
|
|
current_data->natoms++;
|
|
|
|
if ((npositions + nbonds)% NPRINT == 0) ShowStatus();
|
|
}
|
|
|
|
return(BEING_NICE);
|
|
}
|
|
|
|
PRIVATE int ReadBond(FILE *file, LENGTH length)
|
|
{
|
|
BOND *b;
|
|
int n;
|
|
int i;
|
|
|
|
if (!current_data) NewData();
|
|
|
|
n = length / (sizeof(b->type) + sizeof(b->index1) + sizeof(b->index2));
|
|
|
|
if (current_data->maxbonds < n) {
|
|
free(current_data->bonds);
|
|
current_data->maxbonds = n;
|
|
current_data->bonds =
|
|
(BOND *) XtMalloc(n * sizeof(BOND));
|
|
}
|
|
|
|
current_data->nbonds = 0;
|
|
b = current_data->bonds;
|
|
i = 0;
|
|
|
|
for(; i < n; i++, b++) {
|
|
if (1 != fread(&(b->type), sizeof(b->type), 1, file))
|
|
return(ERROR);
|
|
if (1 != fread(&(b->index1), sizeof(b->index1), 1, file))
|
|
return(ERROR);
|
|
if (1 != fread(&(b->index2), sizeof(b->index2), 1, file))
|
|
return(ERROR);
|
|
|
|
if (b->type < 1 || b->type > Common.nbondcolors)
|
|
b->type = 1;
|
|
|
|
nbonds++;
|
|
current_data->nbonds++;
|
|
|
|
if ((npositions + nbonds)% NPRINT == 0) ShowStatus();
|
|
}
|
|
|
|
return(BEING_NICE);
|
|
}
|
|
|
|
/*ARGSUSED*/
|
|
PRIVATE int CopyBond(FILE *file, LENGTH length)
|
|
{
|
|
return( (copy_bonds()) ? READ_OK : ERROR );
|
|
}
|
|
|
|
|
|
|
|
/*ARGSUSED*/
|
|
PRIVATE int ReadTime(FILE *file, LENGTH length)
|
|
{
|
|
if (!current_data || seen_time) NewData();
|
|
seen_time = 1;
|
|
|
|
if (1 != fread(&(current_data->time), sizeof(current_data->time),
|
|
1, file))
|
|
return(ERROR);
|
|
return(READ_OK);
|
|
}
|