1 /********************************************************************** 2 XTCFormat - portable compressed trajectory format (gromacs) 3 4 Copyright (C) 2008 by Tim Vandermeersch 5 6 Parts of the code are from libxdrf: 7 frans van hoesel hoesel@chem.rug.nl 8 http://hpcv100.rc.rug.nl/xdrf.html 9 10 This program is free software; you can redistribute it and/or modify 11 it under the terms of the GNU General Public License as published by 12 the Free Software Foundation version 2 of the License. 13 14 This program is distributed in the hope that it will be useful, 15 but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 GNU General Public License for more details. 18 ***********************************************************************/ 19 20 #include <openbabel/babelconfig.h> 21 #include <openbabel/obmolecformat.h> 22 #include <openbabel/mol.h> 23 24 // This include is not necessary for VS2017 25 #if defined(_MSC_VER) && _MSC_VER <= 1910 26 #include <rpc/types.h> 27 #endif 28 #include <rpc/xdr.h> 29 #include <vector> 30 31 #define MAXID 20 32 #define MAXABS INT_MAX-2 33 34 #ifndef MIN 35 #define MIN(x,y) ((x) < (y) ? (x):(y)) 36 #endif 37 #ifndef MAX 38 #define MAX(x,y) ((x) > (y) ? (x):(y)) 39 #endif 40 #ifndef SQR 41 #define SQR(x) ((x)*(x)) 42 #endif 43 44 #define FIRSTIDX 9 45 /* note that magicints[FIRSTIDX-1] == 0 */ 46 #define LASTIDX (sizeof(magicints) / sizeof(*magicints)) 47 48 namespace OpenBabel 49 { 50 51 class XTCFormat : public OBMoleculeFormat 52 { 53 private: 54 XDR xd; 55 FILE* xdrfiles[MAXID]; 56 XDR* xdridptr[MAXID]; 57 char xdrmodes[MAXID]; 58 unsigned int cnt; 59 60 static int magicints[]; // defined below 61 int xdropen(XDR *xdrs, const char *filename, const char *type); 62 int xdrclose(XDR *xdrs); 63 void sendbits(int buf[], int num_of_bits, int num); 64 int sizeofint(const int size); 65 int sizeofints( const int num_of_ints, unsigned int sizes[]); 66 void sendints(int buf[], const int num_of_ints, const int num_of_bits, 67 unsigned int sizes[], unsigned int nums[]); 68 int receivebits(int buf[], int num_of_bits); 69 void receiveints(int buf[], const int num_of_ints, int num_of_bits, 70 unsigned int sizes[], int nums[]); 71 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision); 72 73 public: 74 //Register this format type ID XTCFormat()75 XTCFormat() 76 { 77 OBConversion::RegisterFormat("xtc",this); 78 } 79 Description()80 virtual const char* Description() //required 81 { 82 return 83 "XTC format\n" 84 "A portable format for trajectories (gromacs)\n"; 85 }; 86 SpecificationURL()87 virtual const char* SpecificationURL() 88 { 89 return "http://manual.gromacs.org/documentation/current/reference-manual/file-formats.html#xtc"; 90 } 91 92 //Flags() can return be any the following combined by | or be omitted if none apply 93 // NOTREADABLE READONEONLY NOTWRITABLE WRITEONEONLY Flags()94 virtual unsigned int Flags() 95 { 96 return NOTWRITABLE; 97 }; 98 99 //*** This section identical for most OBMol conversions *** 100 //////////////////////////////////////////////////// 101 /// The "API" interface functions 102 virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv); 103 //virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv); 104 }; 105 //*** 106 107 // define the static magicints 108 int XTCFormat::magicints[] = { 109 0, 0, 0, 0, 0, 0, 0, 0, 0, 110 8, 10, 12, 16, 20, 25, 32, 40, 50, 64, 111 80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 112 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501, 113 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536, 114 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561, 115 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042, 116 8388607, 10568983, 13316085, 16777216 }; 117 118 //Make an instance of the format class 119 XTCFormat theXTCFormat; 120 121 ///////////////////////////////////////////////////////////////// ReadMolecule(OBBase * pOb,OBConversion * pConv)122 bool XTCFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv) 123 { 124 //OBMol* pmol = pOb->CastAndClear<OBMol>(); 125 OBMol* pmol = dynamic_cast<OBMol*>(pOb); 126 if (pmol == nullptr) 127 return false; 128 129 OBMol &mol = *pmol; 130 std::string filename = pConv->GetInFilename(); 131 132 if (xdropen(&xd, filename.c_str(),"r") == 0) { 133 std::stringstream errorMsg; 134 errorMsg << "Error while opening " << filename << " for reading."; 135 obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning); 136 return false; 137 } 138 139 int magic, natoms, step; 140 float prec = 1000.0; 141 float time, box[3][3]; 142 std::vector<float> floatCoord; 143 std::vector<double*> vconf; 144 145 146 while (true) { 147 // Check the magic int starting each frame 148 xdr_int(&xd, &magic); 149 if (magic != 1995) { 150 std::stringstream errorMsg; 151 errorMsg << "Error: magic int is " << magic << ", should be 1995."; 152 obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning); 153 return false; 154 } 155 156 // Get the number of atoms 157 xdr_int(&xd, &natoms); 158 if (natoms != mol.NumAtoms()) { 159 std::stringstream errorMsg; 160 errorMsg << "Error: number of atoms in the trajectory (" << natoms 161 << ") doesn't match the number of atoms in the supplied " 162 << "molecule (" << mol.NumAtoms() << ")."; 163 obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning); 164 return false; 165 } 166 167 // Get the frame step number 168 xdr_int(&xd, &step); 169 // Get the frame time 170 xdr_float(&xd, &time); 171 // Get the frame box 172 xdr_float(&xd, &box[0][0]); 173 xdr_float(&xd, &box[0][1]); 174 xdr_float(&xd, &box[0][2]); 175 xdr_float(&xd, &box[1][0]); 176 xdr_float(&xd, &box[1][1]); 177 xdr_float(&xd, &box[1][2]); 178 xdr_float(&xd, &box[2][0]); 179 xdr_float(&xd, &box[2][1]); 180 xdr_float(&xd, &box[2][2]); 181 182 if (!floatCoord.size()) 183 floatCoord.resize(natoms * 3); 184 185 // Read the positions 186 if (xdr3dfcoord(&xd, &floatCoord[0], &natoms, &prec) == 0) { 187 // printf("end reached...\n"); 188 break; 189 } 190 191 // Convert positions from single to double precision and convert from 192 // nm to A 193 double *confs = new double[natoms *3]; 194 for (int i=0; i < natoms * 3; ++i) // unroll?? 195 confs[i] = static_cast<double>(10.0 * floatCoord.at(i)); 196 197 vconf.push_back(confs); 198 } 199 200 // Close the XDR file 201 xdrclose(&xd); 202 203 // Set the conformers in the mol object 204 mol.SetConformers(vconf); 205 206 return(true); 207 } 208 209 /*____________________________________________________________________________ 210 | 211 | libxdrf - portable fortran interface to xdr. some xdr routines 212 | are C routines for compressed coordinates 213 | 214 | version 1.1 215 | 216 | This collection of routines is intended to write and read 217 | data in a portable way to a file, so data written on one type 218 | of machine can be read back on a different type. 219 | ________________________________________________________________________ 220 | 221 | Below are the routines to be used by C programmers. Use the 'normal' 222 | xdr routines to write integers, floats, etc (see man xdr) 223 | 224 | int xdropen(XDR *xdrs, const char *filename, const char *type) 225 | This will open the file with the given filename and the 226 | given mode. You should pass it an allocated XDR struct 227 | in xdrs, to be used in all other calls to xdr routines. 228 | Mode is 'w' to create, or update an file, and for all 229 | other values of mode the file is opened for reading. 230 | You need to call xdrclose to flush the output and close 231 | the file. 232 | 233 | Note that you should not use xdrstdio_create, which 234 | comes with the standard xdr library. 235 | 236 | int xdrclose(XDR *xdrs) 237 | Flush the data to the file, and close the file; 238 | You should not use xdr_destroy (which comes standard 239 | with the xdr libraries). 240 | 241 | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) 242 | This is \fInot\fR a standard xdr routine. I named it this 243 | way, because it invites people to use the other xdr 244 | routines. 245 | 246 | frans van hoesel hoesel@chem.rug.nl 247 */ 248 249 250 /*___________________________________________________________________________ 251 | 252 | what follows are the C routines for opening, closing xdr streams 253 | and the routine to read/write compressed coordinates together 254 | with some routines to assist in this task (those are marked 255 | static and cannot be called from user programs) 256 */ 257 258 /*__________________________________________________________________________ 259 | 260 | xdropen - open xdr file 261 | 262 | This versions differs from xdrstdio_create, because I need to know 263 | the state of the file (read or write) so I can use xdr3dfcoord 264 | in eigther read or write mode, and the file descriptor 265 | so I can close the file (something xdr_destroy doesn't do). 266 | 267 */ 268 xdropen(XDR * xdrs,const char * filename,const char * type)269 int XTCFormat::xdropen(XDR *xdrs, const char *filename, const char *type) { 270 int init_done = 0; 271 enum xdr_op lmode; 272 int xdrid; 273 274 if (init_done == 0) { 275 for (xdrid = 1; xdrid < MAXID; xdrid++) { 276 xdridptr[xdrid] = nullptr; 277 } 278 init_done = 1; 279 } 280 xdrid = 1; 281 while (xdrid < MAXID && xdridptr[xdrid] != nullptr) { 282 xdrid++; 283 } 284 if (xdrid == MAXID) { 285 return 0; 286 } 287 if (*type == 'w' || *type == 'W') { 288 type = "w+"; 289 lmode = XDR_ENCODE; 290 } else { 291 type = "r"; 292 lmode = XDR_DECODE; 293 } 294 xdrfiles[xdrid] = fopen(filename, type); 295 if (xdrfiles[xdrid] == nullptr) { 296 xdrs = nullptr; 297 return 0; 298 } 299 xdrmodes[xdrid] = *type; 300 /* next test isn't useful in the case of C language 301 * but is used for the Fortran interface 302 * (C users are expected to pass the address of an already allocated 303 * XDR staructure) 304 */ 305 if (xdrs == nullptr) { 306 xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR)); 307 xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode); 308 } else { 309 xdridptr[xdrid] = xdrs; 310 xdrstdio_create(xdrs, xdrfiles[xdrid], lmode); 311 } 312 return xdrid; 313 } 314 315 /*_________________________________________________________________________ 316 | 317 | xdrclose - close a xdr file 318 | 319 | This will flush the xdr buffers, and destroy the xdr stream. 320 | It also closes the associated file descriptor (this is *not* 321 | done by xdr_destroy). 322 | 323 */ 324 xdrclose(XDR * xdrs)325 int XTCFormat::xdrclose(XDR *xdrs) { 326 int xdrid; 327 328 if (xdrs == nullptr) { 329 fprintf(stderr, "xdrclose: passed a NULL pointer\n"); 330 return 0; 331 } 332 for (xdrid = 1; xdrid < MAXID; xdrid++) { 333 if (xdridptr[xdrid] == xdrs) { 334 335 xdr_destroy(xdrs); 336 fclose(xdrfiles[xdrid]); 337 xdridptr[xdrid] = nullptr; 338 return 1; 339 } 340 } 341 fprintf(stderr, "xdrclose: no such open xdr file\n"); 342 return 0; 343 344 } 345 346 /*____________________________________________________________________________ 347 | 348 | sendbits - encode num into buf using the specified number of bits 349 | 350 | This routines appends the value of num to the bits already present in 351 | the array buf. You need to give it the number of bits to use and you 352 | better make sure that this number of bits is enough to hold the value 353 | Also num must be positive. 354 | 355 */ 356 sendbits(int buf[],int num_of_bits,int num)357 void XTCFormat::sendbits(int buf[], int num_of_bits, int num) { 358 359 unsigned int cnt, lastbyte; 360 int lastbits; 361 unsigned char * cbuf; 362 363 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf); 364 cnt = (unsigned int) buf[0]; 365 lastbits = buf[1]; 366 lastbyte =(unsigned int) buf[2]; 367 while (num_of_bits >= 8) { 368 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/); 369 cbuf[cnt++] = lastbyte >> lastbits; 370 num_of_bits -= 8; 371 } 372 if (num_of_bits > 0) { 373 lastbyte = (lastbyte << num_of_bits) | num; 374 lastbits += num_of_bits; 375 if (lastbits >= 8) { 376 lastbits -= 8; 377 cbuf[cnt++] = lastbyte >> lastbits; 378 } 379 } 380 buf[0] = cnt; 381 buf[1] = lastbits; 382 buf[2] = lastbyte; 383 if (lastbits>0) { 384 cbuf[cnt] = lastbyte << (8 - lastbits); 385 } 386 } 387 388 /*_________________________________________________________________________ 389 | 390 | sizeofint - calculate bitsize of an integer 391 | 392 | return the number of bits needed to store an integer with given max size 393 | 394 */ 395 sizeofint(const int size)396 int XTCFormat::sizeofint(const int size) { 397 unsigned int num = 1; 398 int num_of_bits = 0; 399 400 while (size >= num && num_of_bits < 32) { 401 num_of_bits++; 402 num <<= 1; 403 } 404 return num_of_bits; 405 } 406 407 /*___________________________________________________________________________ 408 | 409 | sizeofints - calculate 'bitsize' of compressed ints 410 | 411 | given the number of small unsigned integers and the maximum value 412 | return the number of bits needed to read or write them with the 413 | routines receiveints and sendints. You need this parameter when 414 | calling these routines. Note that for many calls I can use 415 | the variable 'smallidx' which is exactly the number of bits, and 416 | So I don't need to call 'sizeofints for those calls. 417 */ 418 sizeofints(const int num_of_ints,unsigned int sizes[])419 int XTCFormat::sizeofints( const int num_of_ints, unsigned int sizes[]) { 420 int i, num; 421 unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp; 422 num_of_bytes = 1; 423 bytes[0] = 1; 424 num_of_bits = 0; 425 for (i=0; i < num_of_ints; i++) { 426 tmp = 0; 427 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) { 428 tmp = bytes[bytecnt] * sizes[i] + tmp; 429 bytes[bytecnt] = tmp & 0xff; 430 tmp >>= 8; 431 } 432 while (tmp != 0) { 433 bytes[bytecnt++] = tmp & 0xff; 434 tmp >>= 8; 435 } 436 num_of_bytes = bytecnt; 437 } 438 num = 1; 439 num_of_bytes--; 440 while (bytes[num_of_bytes] >= num) { 441 num_of_bits++; 442 num *= 2; 443 } 444 return num_of_bits + num_of_bytes * 8; 445 446 } 447 448 /*____________________________________________________________________________ 449 | 450 | sendints - send a small set of small integers in compressed format 451 | 452 | this routine is used internally by xdr3dfcoord, to send a set of 453 | small integers to the buffer. 454 | Multiplication with fixed (specified maximum ) sizes is used to get 455 | to one big, multibyte integer. Allthough the routine could be 456 | modified to handle sizes bigger than 16777216, or more than just 457 | a few integers, this is not done, because the gain in compression 458 | isn't worth the effort. Note that overflowing the multiplication 459 | or the byte buffer (32 bytes) is unchecked and causes bad results. 460 | 461 */ 462 sendints(int buf[],const int num_of_ints,const int num_of_bits,unsigned int sizes[],unsigned int nums[])463 void XTCFormat::sendints(int buf[], const int num_of_ints, const int num_of_bits, 464 unsigned int sizes[], unsigned int nums[]) { 465 466 int i; 467 unsigned int bytes[32], num_of_bytes, bytecnt, tmp; 468 469 tmp = nums[0]; 470 num_of_bytes = 0; 471 do { 472 bytes[num_of_bytes++] = tmp & 0xff; 473 tmp >>= 8; 474 } while (tmp != 0); 475 476 for (i = 1; i < num_of_ints; i++) { 477 if (nums[i] >= sizes[i]) { 478 fprintf(stderr,"major breakdown in sendints num %d doesn't " 479 "match size %d\n", nums[i], sizes[i]); 480 return; 481 } 482 /* use one step multiply */ 483 tmp = nums[i]; 484 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) { 485 tmp = bytes[bytecnt] * sizes[i] + tmp; 486 bytes[bytecnt] = tmp & 0xff; 487 tmp >>= 8; 488 } 489 while (tmp != 0) { 490 bytes[bytecnt++] = tmp & 0xff; 491 tmp >>= 8; 492 } 493 num_of_bytes = bytecnt; 494 } 495 if (num_of_bits >= num_of_bytes * 8) { 496 for (i = 0; i < num_of_bytes; i++) { 497 sendbits(buf, 8, bytes[i]); 498 } 499 sendbits(buf, num_of_bits - num_of_bytes * 8, 0); 500 } else { 501 for (i = 0; i < num_of_bytes-1; i++) { 502 sendbits(buf, 8, bytes[i]); 503 } 504 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]); 505 } 506 } 507 508 509 /*___________________________________________________________________________ 510 | 511 | receivebits - decode number from buf using specified number of bits 512 | 513 | extract the number of bits from the array buf and construct an integer 514 | from it. Return that value. 515 | 516 */ 517 receivebits(int buf[],int num_of_bits)518 int XTCFormat::receivebits(int buf[], int num_of_bits) { 519 int cnt, num; 520 unsigned int lastbits, lastbyte; 521 unsigned char * cbuf; 522 int mask = (1 << num_of_bits) -1; 523 524 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf); 525 cnt = buf[0]; 526 lastbits = (unsigned int) buf[1]; 527 lastbyte = (unsigned int) buf[2]; 528 529 num = 0; 530 while (num_of_bits >= 8) { 531 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++]; 532 num |= (lastbyte >> lastbits) << (num_of_bits - 8); 533 num_of_bits -=8; 534 } 535 if (num_of_bits > 0) { 536 if (lastbits < num_of_bits) { 537 lastbits += 8; 538 lastbyte = (lastbyte << 8) | cbuf[cnt++]; 539 } 540 lastbits -= num_of_bits; 541 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1); 542 } 543 num &= mask; 544 buf[0] = cnt; 545 buf[1] = lastbits; 546 buf[2] = lastbyte; 547 return num; 548 } 549 550 /*____________________________________________________________________________ 551 | 552 | receiveints - decode 'small' integers from the buf array 553 | 554 | this routine is the inverse from sendints() and decodes the small integers 555 | written to buf by calculating the remainder and doing divisions with 556 | the given sizes[]. You need to specify the total number of bits to be 557 | used from buf in num_of_bits. 558 | 559 */ 560 receiveints(int buf[],const int num_of_ints,int num_of_bits,unsigned int sizes[],int nums[])561 void XTCFormat::receiveints(int buf[], const int num_of_ints, int num_of_bits, 562 unsigned int sizes[], int nums[]) { 563 int bytes[32]; 564 int i, j, num_of_bytes, p, num; 565 566 bytes[1] = bytes[2] = bytes[3] = 0; 567 num_of_bytes = 0; 568 while (num_of_bits > 8) { 569 bytes[num_of_bytes++] = receivebits(buf, 8); 570 num_of_bits -= 8; 571 } 572 if (num_of_bits > 0) { 573 bytes[num_of_bytes++] = receivebits(buf, num_of_bits); 574 } 575 for (i = num_of_ints-1; i > 0; i--) { 576 num = 0; 577 for (j = num_of_bytes-1; j >=0; j--) { 578 num = (num << 8) | bytes[j]; 579 p = num / sizes[i]; 580 bytes[j] = p; 581 num = num - p * sizes[i]; 582 } 583 nums[i] = num; 584 } 585 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24); 586 } 587 588 /*____________________________________________________________________________ 589 | 590 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file. 591 | 592 | this routine reads or writes (depending on how you opened the file with 593 | xdropen() ) a large number of 3d coordinates (stored in *fp). 594 | The number of coordinates triplets to write is given by *size. On 595 | read this number may be zero, in which case it reads as many as were written 596 | or it may specify the number if triplets to read (which should match the 597 | number written). 598 | Compression is achieved by first converting all floating numbers to integer 599 | using multiplication by *precision and rounding to the nearest integer. 600 | Then the minimum and maximum value are calculated to determine the range. 601 | The limited range of integers so found, is used to compress the coordinates. 602 | In addition the differences between successive coordinates is calculated. 603 | If the difference happens to be 'small' then only the difference is saved, 604 | compressing the data even more. The notion of 'small' is changed dynamically 605 | and is enlarged or reduced whenever needed or possible. 606 | Extra compression is achieved in the case of GROMOS and coordinates of 607 | water molecules. GROMOS first writes out the Oxygen position, followed by 608 | the two hydrogens. In order to make the differences smaller (and thereby 609 | compression the data better) the order is changed into first one hydrogen 610 | then the oxygen, followed by the other hydrogen. This is rather special, but 611 | it shouldn't harm in the general case. 612 | 613 */ 614 xdr3dfcoord(XDR * xdrs,float * fp,int * size,float * precision)615 int XTCFormat::xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) { 616 int *ip = nullptr; 617 int oldsize = 0; 618 int *buf = nullptr; 619 620 int minint[3], maxint[3], mindiff, *lip, diff; 621 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx; 622 int minidx, maxidx; 623 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip; 624 int flag, k; 625 int small_, smaller, larger, i; 626 int is_small, is_smaller, run, prevrun; 627 float *lfp, lf; 628 int tmp, *thiscoord, prevcoord[3]; 629 unsigned int tmpcoord[30]; 630 631 int bufsize, xdrid, lsize; 632 unsigned int bitsize; 633 float inv_precision; 634 int errval = 1; 635 636 /* find out if xdrs is opened for reading or for writing */ 637 xdrid = 0; 638 while (xdridptr[xdrid] != xdrs) { 639 xdrid++; 640 if (xdrid >= MAXID) { 641 fprintf(stderr, "xdr error. no open xdr stream\n"); 642 return 0; 643 } 644 } 645 if (xdrmodes[xdrid] == 'w') { 646 647 /* xdrs is open for writing */ 648 649 if (xdr_int(xdrs, size) == 0) 650 return 0; 651 size3 = *size * 3; 652 /* when the number of coordinates is small, don't try to compress; just 653 * write them as floats using xdr_vector 654 */ 655 if (*size <= 9 ) { 656 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp), 657 (xdrproc_t)xdr_float)); 658 } 659 660 xdr_float(xdrs, precision); 661 if (ip == nullptr) { 662 ip = (int *)malloc(size3 * sizeof(*ip)); 663 if (ip == nullptr) { 664 fprintf(stderr,"malloc failed\n"); 665 return 0; 666 } 667 bufsize = static_cast<int> (size3 * 1.2); 668 buf = (int *)malloc(bufsize * sizeof(*buf)); 669 if (buf == nullptr) { 670 fprintf(stderr,"malloc failed\n"); 671 return 0; 672 } 673 oldsize = *size; 674 } else if (*size > oldsize) { 675 ip = (int *)realloc(ip, size3 * sizeof(*ip)); 676 if (ip == nullptr) { 677 fprintf(stderr,"malloc failed\n"); 678 return 0; 679 } 680 bufsize = static_cast<int> (size3 * 1.2); 681 buf = (int *)realloc(buf, bufsize * sizeof(*buf)); 682 if (buf == nullptr) { 683 fprintf(stderr,"malloc failed\n"); 684 return 0; 685 } 686 oldsize = *size; 687 } 688 /* buf[0-2] are special and do not contain actual data */ 689 buf[0] = buf[1] = buf[2] = 0; 690 minint[0] = minint[1] = minint[2] = INT_MAX; 691 maxint[0] = maxint[1] = maxint[2] = INT_MIN; 692 prevrun = -1; 693 lfp = fp; 694 lip = ip; 695 mindiff = INT_MAX; 696 oldlint1 = oldlint2 = oldlint3 = 0; 697 while(lfp < fp + size3 ) { 698 /* find nearest integer */ 699 if (*lfp >= 0.0) 700 lf = *lfp * *precision + 0.5f; 701 else 702 lf = *lfp * *precision - 0.5f; 703 if (fabs(lf) > MAXABS) { 704 /* scaling would cause overflow */ 705 errval = 0; 706 } 707 lint1 = static_cast<int> (lf); 708 if (lint1 < minint[0]) minint[0] = lint1; 709 if (lint1 > maxint[0]) maxint[0] = lint1; 710 *lip++ = lint1; 711 lfp++; 712 if (*lfp >= 0.0) 713 lf = *lfp * *precision + 0.5f; 714 else 715 lf = *lfp * *precision - 0.5f; 716 if (fabs(lf) > MAXABS) { 717 /* scaling would cause overflow */ 718 errval = 0; 719 } 720 lint2 = static_cast<int> (lf); 721 if (lint2 < minint[1]) minint[1] = lint2; 722 if (lint2 > maxint[1]) maxint[1] = lint2; 723 *lip++ = lint2; 724 lfp++; 725 if (*lfp >= 0.0) 726 lf = *lfp * *precision + 0.5f; 727 else 728 lf = *lfp * *precision - 0.5f; 729 if (fabs(lf) > MAXABS) { 730 /* scaling would cause overflow */ 731 errval = 0; 732 } 733 lint3 = static_cast<int> (lf); 734 if (lint3 < minint[2]) minint[2] = lint3; 735 if (lint3 > maxint[2]) maxint[2] = lint3; 736 *lip++ = lint3; 737 lfp++; 738 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3); 739 if (diff < mindiff && lfp > fp + 3) 740 mindiff = diff; 741 oldlint1 = lint1; 742 oldlint2 = lint2; 743 oldlint3 = lint3; 744 } 745 xdr_int(xdrs, &(minint[0])); 746 xdr_int(xdrs, &(minint[1])); 747 xdr_int(xdrs, &(minint[2])); 748 749 xdr_int(xdrs, &(maxint[0])); 750 xdr_int(xdrs, &(maxint[1])); 751 xdr_int(xdrs, &(maxint[2])); 752 753 if ((float)maxint[0] - (float)minint[0] >= MAXABS || 754 (float)maxint[1] - (float)minint[1] >= MAXABS || 755 (float)maxint[2] - (float)minint[2] >= MAXABS) { 756 /* turning value in unsigned by subtracting minint 757 * would cause overflow 758 */ 759 errval = 0; 760 } 761 sizeint[0] = maxint[0] - minint[0]+1; 762 sizeint[1] = maxint[1] - minint[1]+1; 763 sizeint[2] = maxint[2] - minint[2]+1; 764 765 /* check if one of the sizes is to big to be multiplied */ 766 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) { 767 bitsizeint[0] = sizeofint(sizeint[0]); 768 bitsizeint[1] = sizeofint(sizeint[1]); 769 bitsizeint[2] = sizeofint(sizeint[2]); 770 bitsize = 0; /* flag the use of large sizes */ 771 } else { 772 bitsize = sizeofints(3, sizeint); 773 } 774 lip = ip; 775 luip = (unsigned int *) ip; 776 smallidx = FIRSTIDX; 777 while (smallidx < LASTIDX && magicints[smallidx] < mindiff) { 778 smallidx++; 779 } 780 xdr_int(xdrs, &smallidx); 781 maxidx = MIN(LASTIDX, smallidx + 8) ; 782 minidx = maxidx - 8; /* often this equal smallidx */ 783 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2; 784 small_ = magicints[smallidx] / 2; 785 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; 786 larger = magicints[maxidx] / 2; 787 i = 0; 788 while (i < *size) { 789 is_small = 0; 790 thiscoord = (int *)(luip) + i * 3; 791 if (smallidx < maxidx && i >= 1 && 792 abs(thiscoord[0] - prevcoord[0]) < larger && 793 abs(thiscoord[1] - prevcoord[1]) < larger && 794 abs(thiscoord[2] - prevcoord[2]) < larger) { 795 is_smaller = 1; 796 } else if (smallidx > minidx) { 797 is_smaller = -1; 798 } else { 799 is_smaller = 0; 800 } 801 if (i + 1 < *size) { 802 if (abs(thiscoord[0] - thiscoord[3]) < small_ && 803 abs(thiscoord[1] - thiscoord[4]) < small_ && 804 abs(thiscoord[2] - thiscoord[5]) < small_) { 805 /* interchange first with second atom for better 806 * compression of water molecules 807 */ 808 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3]; 809 thiscoord[3] = tmp; 810 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4]; 811 thiscoord[4] = tmp; 812 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5]; 813 thiscoord[5] = tmp; 814 is_small = 1; 815 } 816 817 } 818 tmpcoord[0] = thiscoord[0] - minint[0]; 819 tmpcoord[1] = thiscoord[1] - minint[1]; 820 tmpcoord[2] = thiscoord[2] - minint[2]; 821 if (bitsize == 0) { 822 sendbits(buf, bitsizeint[0], tmpcoord[0]); 823 sendbits(buf, bitsizeint[1], tmpcoord[1]); 824 sendbits(buf, bitsizeint[2], tmpcoord[2]); 825 } else { 826 sendints(buf, 3, bitsize, sizeint, tmpcoord); 827 } 828 prevcoord[0] = thiscoord[0]; 829 prevcoord[1] = thiscoord[1]; 830 prevcoord[2] = thiscoord[2]; 831 thiscoord = thiscoord + 3; 832 i++; 833 834 run = 0; 835 if (is_small == 0 && is_smaller == -1) 836 is_smaller = 0; 837 while (is_small && run < 8*3) { 838 if (is_smaller == -1 && ( 839 SQR(thiscoord[0] - prevcoord[0]) + 840 SQR(thiscoord[1] - prevcoord[1]) + 841 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) { 842 is_smaller = 0; 843 } 844 845 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small_; 846 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small_; 847 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small_; 848 849 prevcoord[0] = thiscoord[0]; 850 prevcoord[1] = thiscoord[1]; 851 prevcoord[2] = thiscoord[2]; 852 853 i++; 854 thiscoord = thiscoord + 3; 855 is_small = 0; 856 if (i < *size && 857 abs(thiscoord[0] - prevcoord[0]) < small_ && 858 abs(thiscoord[1] - prevcoord[1]) < small_ && 859 abs(thiscoord[2] - prevcoord[2]) < small_) { 860 is_small = 1; 861 } 862 } 863 if (run != prevrun || is_smaller != 0) { 864 prevrun = run; 865 sendbits(buf, 1, 1); /* flag the change in run-length */ 866 sendbits(buf, 5, run+is_smaller+1); 867 } else { 868 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */ 869 } 870 for (k=0; k < run; k+=3) { 871 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]); 872 } 873 if (is_smaller != 0) { 874 smallidx += is_smaller; 875 if (is_smaller < 0) { 876 small_ = smaller; 877 smaller = magicints[smallidx-1] / 2; 878 } else { 879 smaller = small_; 880 small_ = magicints[smallidx] / 2; 881 } 882 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; 883 } 884 } 885 if (buf[1] != 0) buf[0]++;; 886 xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */ 887 return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0])); 888 } else { 889 890 /* xdrs is open for reading */ 891 892 if (xdr_int(xdrs, &lsize) == 0) 893 return 0; 894 if (*size != 0 && lsize != *size) { 895 fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; " 896 "%d arg vs %d in file", *size, lsize); 897 } 898 *size = lsize; 899 size3 = *size * 3; 900 if (*size <= 9) { 901 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp), 902 (xdrproc_t)xdr_float)); 903 } 904 xdr_float(xdrs, precision); 905 if (ip == nullptr) { 906 ip = (int *)malloc(size3 * sizeof(*ip)); 907 if (ip == nullptr) { 908 fprintf(stderr,"malloc failed\n"); 909 return 0; 910 } 911 bufsize = static_cast<int> (size3 * 1.2); 912 buf = (int *)malloc(bufsize * sizeof(*buf)); 913 if (buf == nullptr) { 914 fprintf(stderr,"malloc failed\n"); 915 return 0; 916 } 917 oldsize = *size; 918 } else if (*size > oldsize) { 919 ip = (int *)realloc(ip, size3 * sizeof(*ip)); 920 if (ip == nullptr) { 921 fprintf(stderr,"malloc failed\n"); 922 return 0; 923 } 924 bufsize = static_cast<int> (size3 * 1.2); 925 buf = (int *)realloc(buf, bufsize * sizeof(*buf)); 926 if (buf == nullptr) { 927 fprintf(stderr,"malloc failed\n"); 928 return 0; 929 } 930 oldsize = *size; 931 } 932 buf[0] = buf[1] = buf[2] = 0; 933 934 xdr_int(xdrs, &(minint[0])); 935 xdr_int(xdrs, &(minint[1])); 936 xdr_int(xdrs, &(minint[2])); 937 938 xdr_int(xdrs, &(maxint[0])); 939 xdr_int(xdrs, &(maxint[1])); 940 xdr_int(xdrs, &(maxint[2])); 941 942 sizeint[0] = maxint[0] - minint[0]+1; 943 sizeint[1] = maxint[1] - minint[1]+1; 944 sizeint[2] = maxint[2] - minint[2]+1; 945 946 /* check if one of the sizes is to big to be multiplied */ 947 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) { 948 bitsizeint[0] = sizeofint(sizeint[0]); 949 bitsizeint[1] = sizeofint(sizeint[1]); 950 bitsizeint[2] = sizeofint(sizeint[2]); 951 bitsize = 0; /* flag the use of large sizes */ 952 } else { 953 bitsize = sizeofints(3, sizeint); 954 } 955 956 xdr_int(xdrs, &smallidx); 957 maxidx = MIN(LASTIDX, smallidx + 8) ; 958 minidx = maxidx - 8; /* often this equal smallidx */ 959 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2; 960 small_ = magicints[smallidx] / 2; 961 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; 962 larger = magicints[maxidx]; 963 964 /* buf[0] holds the length in bytes */ 965 966 if (xdr_int(xdrs, &(buf[0])) == 0) 967 return 0; 968 if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0) 969 return 0; 970 buf[0] = buf[1] = buf[2] = 0; 971 972 lfp = fp; 973 inv_precision = 1.0f / * precision; 974 run = 0; 975 i = 0; 976 lip = ip; 977 while ( i < lsize ) { 978 thiscoord = (int *)(lip) + i * 3; 979 980 if (bitsize == 0) { 981 thiscoord[0] = receivebits(buf, bitsizeint[0]); 982 thiscoord[1] = receivebits(buf, bitsizeint[1]); 983 thiscoord[2] = receivebits(buf, bitsizeint[2]); 984 } else { 985 receiveints(buf, 3, bitsize, sizeint, thiscoord); 986 } 987 988 i++; 989 thiscoord[0] += minint[0]; 990 thiscoord[1] += minint[1]; 991 thiscoord[2] += minint[2]; 992 993 prevcoord[0] = thiscoord[0]; 994 prevcoord[1] = thiscoord[1]; 995 prevcoord[2] = thiscoord[2]; 996 997 998 flag = receivebits(buf, 1); 999 is_smaller = 0; 1000 if (flag == 1) { 1001 run = receivebits(buf, 5); 1002 is_smaller = run % 3; 1003 run -= is_smaller; 1004 is_smaller--; 1005 } 1006 if (run > 0) { 1007 thiscoord += 3; 1008 for (k = 0; k < run; k+=3) { 1009 receiveints(buf, 3, smallidx, sizesmall, thiscoord); 1010 i++; 1011 thiscoord[0] += prevcoord[0] - small_; 1012 thiscoord[1] += prevcoord[1] - small_; 1013 thiscoord[2] += prevcoord[2] - small_; 1014 if (k == 0) { 1015 /* interchange first with second atom for better 1016 * compression of water molecules 1017 */ 1018 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0]; 1019 prevcoord[0] = tmp; 1020 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1]; 1021 prevcoord[1] = tmp; 1022 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2]; 1023 prevcoord[2] = tmp; 1024 *lfp++ = prevcoord[0] * inv_precision; 1025 *lfp++ = prevcoord[1] * inv_precision; 1026 *lfp++ = prevcoord[2] * inv_precision; 1027 } else { 1028 prevcoord[0] = thiscoord[0]; 1029 prevcoord[1] = thiscoord[1]; 1030 prevcoord[2] = thiscoord[2]; 1031 } 1032 *lfp++ = thiscoord[0] * inv_precision; 1033 *lfp++ = thiscoord[1] * inv_precision; 1034 *lfp++ = thiscoord[2] * inv_precision; 1035 } 1036 } else { 1037 *lfp++ = thiscoord[0] * inv_precision; 1038 *lfp++ = thiscoord[1] * inv_precision; 1039 *lfp++ = thiscoord[2] * inv_precision; 1040 } 1041 smallidx += is_smaller; 1042 if (is_smaller < 0) { 1043 small_ = smaller; 1044 if (smallidx > FIRSTIDX) { 1045 smaller = magicints[smallidx - 1] /2; 1046 } else { 1047 smaller = 0; 1048 } 1049 } else if (is_smaller > 0) { 1050 smaller = small_; 1051 small_ = magicints[smallidx] / 2; 1052 } 1053 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; 1054 } 1055 } 1056 return 1; 1057 } 1058 1059 } //namespace OpenBabel 1060