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

tprplugin.C

Go to the documentation of this file.
00001 #define VMDPLUGIN_STATIC
00002 #include "molfile_plugin.h"
00003 #include <math.h>
00004 #include <stdio.h>
00005 #include <stdlib.h>
00006 
00007 
00008 /* I am making a very concious decision here. I'm only going to TRY supporting gromacs TPX files that were written by
00009 version 4.0 or later. This means TPX format version 58 and later.
00010 */
00011 
00012 #include "Gromacs.h"
00013 
00014 #define STRLEN 4096
00015 #define SAVELEN 8
00016 #define TRUE (1)
00017 #define FALSE (0)
00018 
00019 #include "tprplugin.h"
00020 
00021 int readff(md_file *mf, int version) {
00022     int atnr, ntypes, i, j;
00023     int *functype;
00024     double reppow =12.0;
00025     float fudge;
00026     if (trx_int(mf, &atnr)) return MOLFILE_ERROR;
00027     #ifdef TPRDEBUG
00028     printf("Atnr: %d\n", atnr);
00029     #endif
00030     if (trx_int(mf, &ntypes)) return MOLFILE_ERROR;
00031     #ifdef TPRDEBUG
00032     printf("Ntypes: %d\n", ntypes);
00033     #endif
00034     functype = (int *) malloc(sizeof(int) * ntypes);
00035     for (i=0; i < ntypes; i++)
00036         if (trx_int(mf, &functype[i])) return MOLFILE_ERROR;
00037     if (version >= 66)
00038         if (trx_double(mf, &reppow)) return MOLFILE_ERROR;
00039         if (trx_real(mf, &fudge)) return MOLFILE_ERROR;
00040     for (i = 0; i < ntypes; i++) {
00041         for (j = 0; (j < NFTUPD); j++) {
00042             if (version < ftupd[j].fvnr && functype[i] >= ftupd[j].ftype)
00043                 functype[i] += 1;
00044         }
00045         readparams(mf, version, functype[i]);
00046     }
00047     free(functype);
00048     return MOLFILE_SUCCESS;
00049 }
00050 //do_atomtypes
00051 int read_atomtypes(md_file *mf, int version) {
00052     int i;
00053     int err = 0;
00054     int numtypes;
00055     err |= trx_int(mf, &numtypes);
00056 #ifdef TPRDEBUG
00057     printf("Num atom types:%d\n", numtypes);
00058 #endif
00059     //Read a bunch of stuff that is set optionally.
00060     for (i = 0; i < numtypes; i++) {
00061         if (version < 113) {
00062             err |= trx_real(mf, NULL);//Radius?
00063             err |= trx_real(mf, NULL);//Volume?
00064             err |= trx_real(mf, NULL);//Surface tension?
00065         }
00066         if (version >= 40) {
00067             err |= trx_real(mf, NULL);
00068         }
00069         if (version >= 60 && version < 113) {
00070             err |= trx_real(mf, NULL);
00071             err |= trx_real(mf, NULL);
00072         }
00073     }
00074     return err;
00075 }
00076 //Reading in this case means discarding.
00077 int read_cmap(md_file *mf) {
00078     int i, ngrid, gs;
00079     int err = 0;
00080     err |= trx_int(mf, &ngrid);
00081     err |= trx_int(mf, &gs);
00082     #ifdef TPRDEBUG
00083     printf("ngrid: %d, gs: %d\n", ngrid, gs);
00084     printf("This many reals should be read: %d\n", ngrid * gs * gs);
00085     #endif
00086     for (i = 0; i < ngrid * gs * gs; i++) {
00087         err |= trx_real(mf, NULL);
00088         err |= trx_real(mf, NULL);
00089         err |= trx_real(mf, NULL);
00090         err |= trx_real(mf, NULL);
00091     }
00092     return err;
00093 }
00094 //do_groups
00095 int read_groups(md_file *mf, int ngrp, tprdata *tpr) {
00096         int g, i, j, tmp;
00097         int err = 0;
00098         //do_grps
00099         #ifdef TPRDEBUG
00100         printf("ngrp = %d\n", ngrp);
00101         #endif
00102         for(i = 0; i < ngrp; i++) {
00103                 err |= trx_int(mf, &tmp);
00104                 #ifdef TPRDEBUG
00105                 printf("Number of reads = %d\n", tmp);
00106                 #endif
00107                 for (j = 0; j < tmp; j++) {
00108                         trx_int(mf, NULL);
00109                 }
00110         }
00111         //end do_grps
00112         //ngrpname
00113         err |= trx_int(mf, &tmp);
00114         #ifdef TPRDEBUG
00115         printf("ngrpname = %d\n", tmp);
00116         #endif
00117         for (i = 0; i < tmp; i++) {
00118                 err |= trx_int(mf, &g);
00119                 #ifdef TPRDEBUG
00120                 printf("symtabid = %d\n", g);
00121                 for (j = 0; j < SAVELEN; j++)
00122                         printf("%c", tpr->symtab[SAVELEN*g+j]);
00123                 printf("\n");
00124                 #endif
00125         }
00126         for (g = 0; g < ngrp; g++) {
00127                 err |= trx_int(mf, &tmp);
00128                 #ifdef TPRDEBUG
00129                 printf("This many chars: %d\n", tmp);
00130                 #endif
00131                 if (tmp) {
00132                         //For reasons that aren't clear to me, writer version 27 has a different way of reading/writing the chars.
00133                         if (tpr->wversion >=27) {
00134                                 // Advance the pointer by tmp bytes.
00135                                 fseek(mf->f, tmp, SEEK_CUR);
00136                         }
00137                         else {
00138                                 fseek(mf->f, 4*tmp, SEEK_CUR);
00139                         }
00140                 }
00141         }
00142         return err;
00143 }
00144 
00145 #define MIN(a, b) ((a)<(b)? (a):(b))
00146 
00147 void tpr_save_string(md_file* mf, char* saveloc, int genversion) {
00148     #ifdef _WIN32
00149         long long int len = 0;
00150     #else
00151         long len = 0;
00152     #endif
00153     int i, j;
00154     fpos_t pos;
00155     char buf[STRLEN];
00156     trx_long(mf, &len);
00157     fread(buf, 1, int(len), mf->f);
00158     // GROMACS is weird. Before writer version 27, the reads were always aligned to 4 bytes.
00159     // In subsequent versions, they are not. So to maintain backwards compatability, add an
00160     // extra seek.
00161     if (genversion < 27 && len % 4) {
00162         fseek(mf->f, 4 - (len % 4), SEEK_CUR);
00163     }
00164 
00165     for (i = 0; i < MIN(int(len), (SAVELEN-1)); i++) {
00166         saveloc[i] = buf[i];
00167     }
00168     saveloc[i] = '\0';
00169 
00170 }
00171 
00172 int readtprAfterPrecision (tprdata *tpr) {
00173         //Start at do_tpx, trace down to do_tpxheader
00174         md_file *mf = tpr->mf;
00175         char buf[STRLEN];
00176         tpr->readcoordinates = 0;
00177         long fsize;
00178         float dummy;
00179         int tmp, i, j,k;
00180         unsigned int ui;
00181         int hasIR, hasCoord, hasV, hasF, hasTop, hasDim, hasIntermoleculeBonds, hasGBSA;
00182         int numcmap;
00183         float boxoffset[9];
00184         float boxv[9];
00185         unsigned short sdummy;
00186         bool bClear;
00187         int sum = 0;
00188 
00189         if (trx_int(mf, &(tpr->version))) return MOLFILE_ERROR;
00190         printf("File Format Version: %d\n", tpr->version);
00191 
00192         if (tpr->version >= 77 && tpr->version <= 79) {
00193                 trx_string(mf, buf, STRLEN);
00194         }
00195 
00196         if (trx_int(mf, &(tpr->wversion))) return MOLFILE_ERROR;
00197         printf("Generator Version: %d\n", tpr->wversion);
00198 
00199         //Bailouts if things are too new/we can't guarantee accurately reading them.
00200         if (tpr->wversion > 28 || tpr->version <= 57) {
00201                 printf("Your file cannot be read, as it has version %d, but we can read from version 57 to at least 128.\n", tpr->version);
00202                 printf("The generator version for your file is %d, but we can only read up to 28\n", tpr->wversion);
00203                 return MOLFILE_ERROR;
00204         }
00205         
00206         if (tpr->version >= 81) {
00207                 j = tpr_string(mf, buf, STRLEN);
00208                 fread(buf, 4, 1, mf->f);
00209                 // //I dunno why this is wrong. Should say "release", but it totally doesn't.
00210 
00211         }
00212 
00213         //If we were dealing with TPA files, we'd need to do something here. Not a clue what.
00214         //do_section line in the gromacs source.
00215         if (trx_int(mf, &(tpr->natoms))) return MOLFILE_ERROR;
00216         #ifdef TPRDEBUG
00217         printf("Natoms: %d\n", tpr->natoms);
00218         #endif
00219 
00220         if (trx_int(mf, &(tpr->ngtc))) return MOLFILE_ERROR;
00221         #ifdef TPRDEBUG
00222         printf("Ngtc: %d\n", tpr->ngtc);
00223         #endif
00224         if (tpr->version < 62) {
00225                 if (trx_int(mf, &tmp)) return MOLFILE_ERROR;
00226                 if (trx_real(mf, &dummy)) return MOLFILE_ERROR;
00227         }
00228 
00229         if (tpr->version >=79) {
00230                 if (trx_int(mf, &tmp)) return MOLFILE_ERROR;
00231                 #ifdef TPRDEBUG
00232                 printf("FEP state: %d\n", tmp);
00233                 #endif
00234         }
00235         if (trx_real(mf, &dummy)) return MOLFILE_ERROR;
00236         #ifdef TPRDEBUG
00237         printf("Lambda: %f\n", dummy);
00238         #endif
00239         if (trx_int(mf, &hasIR)) return MOLFILE_ERROR;
00240         #ifdef TPRDEBUG
00241         printf("IR state?: %d\n", hasIR);
00242         #endif
00243         if (trx_int(mf, &hasTop)) return MOLFILE_ERROR;
00244         #ifdef TPRDEBUG
00245         printf("Contains topology: %d\n", hasTop);
00246         #endif
00247         if (trx_int(mf, &hasCoord)) return MOLFILE_ERROR;
00248         #ifdef TPRDEBUG
00249         printf("Contains coordinates: %d\n", hasCoord);
00250         #endif
00251         if (trx_int(mf, &hasV)) return MOLFILE_ERROR;
00252         #ifdef TPRDEBUG
00253         printf("Contains velocities: %d\n", hasV);
00254         #endif
00255         tpr->has_velocities = 0;
00256         if(hasV) tpr->has_velocities = 1;
00257         if (trx_int(mf, &hasF)) return MOLFILE_ERROR;
00258         #ifdef TPRDEBUG
00259         printf("Contains forces: %d\n", hasF);
00260         #endif
00261         if (trx_int(mf, &hasDim)) return MOLFILE_ERROR;
00262         #ifdef TPRDEBUG
00263         printf("Contains dimensions: %d\n", hasDim);
00264         #endif
00265         if (tpr->wversion >= 27) {
00266                 if (trx_long(mf, &fsize)) return MOLFILE_ERROR;
00267                 #ifdef TPRDEBUG
00268                 printf("Filesize: %d\n", fsize);
00269                 #endif
00270         }
00271 
00272         //End of do_tpxheader
00273         if (hasDim) {
00274                 if (trx_rvector(mf, &(tpr->boxdims[0]))) return MOLFILE_ERROR;
00275                 if (trx_rvector(mf, &(tpr->boxdims[3]))) return MOLFILE_ERROR;
00276                 if (trx_rvector(mf, &(tpr->boxdims[6]))) return MOLFILE_ERROR;
00277                 #ifdef TPRDEBUG
00278                 printf("Box dim: (%f %f %f), (%f %f %f), (%f %f %f)\n", tpr->boxdims[0],tpr->boxdims[1],tpr->boxdims[2],
00279                         tpr->boxdims[3],tpr->boxdims[4],tpr->boxdims[5],tpr->boxdims[6],tpr->boxdims[7],tpr->boxdims[8]);
00280                 #endif
00281                 if (tpr->version >= 51) {
00282                         if (trx_rvector(mf, &(boxoffset[0]))) return MOLFILE_ERROR;
00283                         if (trx_rvector(mf, &(boxoffset[3]))) return MOLFILE_ERROR;
00284                         if (trx_rvector(mf, &(boxoffset[6]))) return MOLFILE_ERROR;
00285                 }
00286                 else {
00287                         boxoffset[0] = boxoffset[1] = boxoffset[2] = 0;
00288                 }
00289                 #ifdef TPRDEBUG
00290                 printf("Box offset: (%f %f %f), (%f %f %f), (%f %f %f)\n", boxoffset[0],boxoffset[1],boxoffset[2],
00291                         boxoffset[3],boxoffset[4],boxoffset[5],boxoffset[6],boxoffset[7],boxoffset[8]);
00292                 #endif
00293                 if (trx_rvector(mf, &(boxv[0]))) return MOLFILE_ERROR;
00294                 if (trx_rvector(mf, &(boxv[3]))) return MOLFILE_ERROR;
00295                 if (trx_rvector(mf, &(boxv[6]))) return MOLFILE_ERROR;
00296                 #ifdef TPRDEBUG
00297                 printf("Box vel: (%f %f %f), (%f %f %f), (%f %f %f)\n", boxv[0],boxv[1],boxv[2],
00298                         boxv[3],boxv[4],boxv[5],boxv[6],boxv[7],boxv[8]);
00299                 #endif
00300         }
00301         //Dump the thermostat storage stuff. VMD wouldn't know what to do with it anyway.
00302         if (tpr->ngtc) {
00303                 if (tpr->version < 69) {
00304                         for (i = 0; i < tpr->ngtc; i++) {
00305                                 if (trx_real(mf, &dummy)) return MOLFILE_ERROR;
00306                                 //printf("NGTC %d: %f\n", i, dummy);
00307                         }
00308                 }
00309                 for (i = 0; i < tpr->ngtc; i++) {
00310                         if (trx_real(mf, &dummy)) return MOLFILE_ERROR;
00311                         //printf("NGTC %d: %f\n", i, dummy);
00312                 }
00313         }
00314 
00315         //do_mtop starts here, which starts by reading the symtab (do_symtab)
00316         if (trx_int(mf, &(tpr->symtablen))) return MOLFILE_ERROR;
00317         #ifdef TPRDEBUG
00318         printf("Symtab length: %d\n", tpr->symtablen);
00319         #endif
00320         tpr->symtab = (char*) malloc(SAVELEN * tpr->symtablen * sizeof(char));
00321         for (i = 0; i < tpr->symtablen; i++) {
00322                 #ifdef TPREXTRADEBUG
00323                 printf("i=%d\n", i);
00324                 #endif
00325                 tpr_save_string(mf, &(tpr->symtab[SAVELEN * i]),tpr->wversion);
00326         }
00327         #ifdef TPRDEBUG
00328         for (i = 0; i < tpr->symtablen; i++) {
00329                 for (j = 0; j < SAVELEN; j++)
00330                         printf("%c", tpr->symtab[SAVELEN*i+j]);
00331                 printf("\n");
00332         }
00333         #endif
00334 
00335         if (trx_int(mf, &tmp)) return MOLFILE_ERROR;
00336         /*printf ("System name: ");
00337         for (j = 0; j < SAVELEN; j++)
00338                 printf("%c", tpr->symtab[SAVELEN*tmp+j]);
00339         printf("\n");*/
00340         //Now "read" in the forcefield. This SHOULD be do_ffparams
00341         readff(mf, tpr->version);
00342 
00343         //Read in the type of molecules.
00344         if (trx_int(mf, &(tpr->nmoltypes))) return MOLFILE_ERROR;
00345         #ifdef TPRDEBUG
00346         printf("nmoltypes: %d\n", tpr->nmoltypes);
00347         #endif
00348         tpr->atomsinmol = (int*) malloc(sizeof(int) * tpr->nmoltypes);
00349         //printf("Malloced a thing\n");
00350         tpr->molnames = (int*) malloc(sizeof(int) * tpr->nmoltypes);
00351         tpr->resinmol = (int*) malloc(sizeof(int) * tpr->nmoltypes);
00352         tpr->charges = (float**) malloc(sizeof(float*) * tpr->nmoltypes);
00353         tpr->masses = (float**) malloc(sizeof(float*) * tpr->nmoltypes);
00354         tpr->types = (unsigned short**) malloc(sizeof(unsigned short*) * tpr->nmoltypes);
00355         tpr->ptypes = (int**) malloc(sizeof(int*) * tpr->nmoltypes);
00356         tpr->resids = (int**) malloc(sizeof(int*) * tpr->nmoltypes);
00357         tpr->atomnameids = (int**) malloc(sizeof(int*) * tpr->nmoltypes);
00358         tpr->atomtypeids = (int**) malloc(sizeof(int*) * tpr->nmoltypes);
00359         tpr->atomicnumbers = (int**) malloc(sizeof(int*) * tpr->nmoltypes);
00360         tpr->resnames = (int**) malloc(sizeof(int*) * tpr->nmoltypes);
00361         for (j = 0; j < F_NRE; j++) {
00362                 tpr->interactionlist[j] = (int**) malloc(sizeof(int*) * tpr->nmoltypes);
00363                 tpr->nr[j] = (int*) malloc(sizeof(int) * tpr->nmoltypes);
00364         }
00365         
00366         //printf("Malloced things\n");
00367         //do_moltype
00368         for (i = 0; i < tpr->nmoltypes; i++) {
00369                 if (trx_int(mf, &(tpr->molnames[i]))) return MOLFILE_ERROR;
00370                 #ifdef TPRDEBUG
00371                 printf ("Mol name %d: ", i);
00372                 for (j = 0; j < SAVELEN; j++)
00373                         printf("%c", tpr->symtab[SAVELEN*tpr->molnames[i]+j]);
00374                 printf("\n");
00375                 #endif
00376                 if (trx_int(mf, &(tpr->atomsinmol[i]))) return MOLFILE_ERROR;
00377                 if (trx_int(mf, &(tpr->resinmol[i]))) return MOLFILE_ERROR;
00378                 tpr->charges[i] = (float*) malloc(sizeof(float) * tpr->atomsinmol[i]);
00379                 tpr->masses[i] = (float*) malloc(sizeof(float) * tpr->atomsinmol[i]);
00380                 tpr->types[i] = (unsigned short*) malloc(sizeof(unsigned short) * tpr->atomsinmol[i]);
00381                 tpr->ptypes[i] = (int*) malloc(sizeof(int) * tpr->atomsinmol[i]);
00382                 tpr->resids[i] = (int*) malloc(sizeof(int) * tpr->atomsinmol[i]);
00383                 tpr->atomicnumbers[i] = (int*) malloc(sizeof(int) * tpr->atomsinmol[i]);
00384                 #ifdef TPRDEBUG
00385                 printf("%d atoms, %d residues in molecule %d\n", tpr->atomsinmol[i], tpr->resinmol[i], i);
00386                 #endif
00387                 for (j=0; j < tpr->atomsinmol[i]; j++) {
00388                         if (trx_real(mf, &(tpr->masses[i][j]))) return MOLFILE_ERROR;
00389                         if (trx_real(mf, &(tpr->charges[i][j]))) return MOLFILE_ERROR;
00390                         if (trx_real(mf, NULL)) return MOLFILE_ERROR; //mB
00391                         if (trx_real(mf, NULL)) return MOLFILE_ERROR; //cB
00392 
00393                         if (trx_int(mf, &tmp)) return MOLFILE_ERROR;
00394                         tpr->types[i][j] = (unsigned short) tmp;
00395                         if (trx_int(mf, &tmp)) return MOLFILE_ERROR;//typeB
00396                         if (tpr->wversion >=27) {
00397                                 fseek(mf->f, -4, SEEK_CUR);
00398                         }
00399                         if (trx_int(mf, &(tpr->ptypes[i][j]))) return MOLFILE_ERROR;
00400                         if (trx_int(mf, &(tpr->resids[i][j]))) return MOLFILE_ERROR;
00401                         if (tpr->version >= 52) {
00402                                 if (trx_int(mf, &(tpr->atomicnumbers[i][j]))) return MOLFILE_ERROR;
00403                         }
00404                         #ifdef TPREXTRADEBUG
00405                         if (i == 0)
00406                                 printf("%d: mass %f charge %f type %d ptype %d resid %d, periodic table number %d\n", j,
00407                                         tpr->masses[i][j], tpr->charges[i][j], tpr->types[i][j], tpr->ptypes[i][j],tpr->resids[i][j], tmp);
00408                         #endif
00409                 }
00410                 tpr->atomnameids[i] = (int*) malloc(sizeof(int) * tpr->atomsinmol[i]);
00411                 if (tpr_ivector(mf, tpr->atomnameids[i], tpr->atomsinmol[i])) return MOLFILE_ERROR;
00412                 tpr->atomtypeids[i] = (int*) malloc(sizeof(int) * tpr->atomsinmol[i]);
00413                 if (tpr_ivector(mf, tpr->atomtypeids[i], tpr->atomsinmol[i])) return MOLFILE_ERROR;
00414                 for (j = 0; j < tpr->atomsinmol[i]; j++) {
00415                         if (trx_real(mf, NULL)) return MOLFILE_ERROR;
00416                 }
00417 
00418                 //Read residues
00419                 tpr->resnames[i] = (int*) malloc(sizeof(int) * tpr->resinmol[i]);
00420                 #ifdef TPRDEBUG
00421                 printf("%d residues in mol %d\n", tpr->resinmol[i], i);
00422                 #endif
00423                 for (j = 0; j < tpr->resinmol[i]; j++) {
00424                         if (trx_int(mf, &(tpr->resnames[i][j]))) return MOLFILE_ERROR;
00425                         if (tpr->version >=63) {
00426                                 if (tpr->wversion < 27)
00427                                         fseek(mf->f, 8, SEEK_CUR);
00428                                 else
00429                                         fseek(mf->f, 5, SEEK_CUR);
00430                         }
00431                         else
00432                                 tpr->resnames[i][j] += 1;
00433                         #ifdef TPREXTRADEBUG
00434                         printf("%d\n", tpr->resnames[i][j]);
00435                         if ( i == 0) {
00436                                 for (k = 0; k < SAVELEN; k++)
00437                                         printf("%c", tpr->symtab[SAVELEN*tpr->resnames[i][j]+k]);
00438                                 printf("\n");
00439                         }
00440                         #endif
00441                 }
00442                 //do_ilists
00443                 //printf ("Atoms in mol %d: %d\n", i, tpr->atomsinmol[i]);
00444                 //Here is where we would read the i-lists, however we only need a subset.
00445                 //Based on the stuff in src/gmxlib/ifunc.c, we only need bonds (index 0),
00446                 //angles (index 10), dihedrals (index 18 and 19), impropers (index 21 and 22), and CMAP cross terms (index 25)
00447                 for (j = 0; j < F_NRE; j++) {
00448                         bClear = FALSE;
00449                         for (k = 0; k < NFTUPD; k++) {
00450                                 if (tpr->version < ftupd[k].fvnr && j == ftupd[k].ftype) {
00451                                         bClear = TRUE;
00452                                 }
00453                         }
00454                         if (bClear) {
00455                                 tpr->nr[j][i] = 0;
00456                         }
00457                         else {
00458                                 if (trx_int(mf, &(tpr->nr[j][i]))) return MOLFILE_ERROR;
00459                                 #ifdef TPREXTRADEBUG
00460                                 printf("j, k, interactions: %d, %d, %d\n", j, k, tpr->nr[j][i]);
00461                                 #endif
00462                         }
00463                         if (tpr->nr[j][i] != 0) {
00464                                 tpr->interactionlist[j][i] = (int*) malloc(sizeof(int) * tpr->nr[j][i]);
00465                                 #ifdef TPRDEBUG
00466                                 printf("Interaction id: %d number of interactants %d\n", j, tpr->nr[j][i]);
00467                                 #endif
00468                                 tpr_ivector(mf, tpr->interactionlist[j][i], tpr->nr[j][i]);
00469                                 #ifdef TPREXTRADEBUG
00470                                 for (k=0; k < tpr->nr[j][i]; k++) {
00471                                         printf("%d ", tpr->interactionlist[j][i][k]);
00472                                 }
00473                                 printf("\n");
00474                                 #endif
00475                         }
00476                         else {
00477                                 tpr->interactionlist[j][i] = NULL;
00478                         }
00479                 }
00480                 //end do_ilists
00481 
00482 
00483                 //do_block and do_blocka. Remove these. VMD doesn't know what they are.
00484                 //They refer to atomgroups and exclusions. Useful for dynamics, not useful for visualization.
00485                 if (trx_int(mf, &k)) return MOLFILE_ERROR;
00486                 for (j = 0; j <= k; j++) {
00487                         if (trx_int(mf, NULL)) return MOLFILE_ERROR;
00488                 }
00489                 if (trx_int(mf, &k)) return MOLFILE_ERROR;
00490                 if (trx_int(mf, &tmp)) return MOLFILE_ERROR;
00491                 for (j = 0; j <= k; j++) {
00492                         if (trx_int(mf, NULL)) return MOLFILE_ERROR;
00493                 }
00494                 for (j = 0; j < tmp; j++) {
00495                         if (trx_int(mf, NULL)) return MOLFILE_ERROR;
00496                 }
00497         }
00498         //end of do_moltype
00499 
00500         if (trx_int(mf, &(tpr->nmolblock))) return MOLFILE_ERROR;
00501         #ifdef TPRDEBUG
00502         printf("nmolblock: %d\n", tpr->nmolblock);
00503         #endif
00504         //do_molblock
00505         tpr->molbtype = (int*) malloc(tpr->nmolblock * sizeof(int));
00506         tpr->molbnmol = (int*) malloc(tpr->nmolblock * sizeof(int));
00507         tpr->molbnatoms = (int*) malloc(tpr->nmolblock * sizeof(int));
00508         for (i = 0; i < tpr->nmolblock; i++) {
00509                 if (trx_int(mf, &(tpr->molbtype[i]))) return MOLFILE_ERROR;
00510                 if (trx_int(mf, &(tpr->molbnmol[i]))) return MOLFILE_ERROR;
00511                 if (trx_int(mf, &(tpr->molbnatoms[i]))) return MOLFILE_ERROR;
00512                 if (trx_int(mf, &k)) return MOLFILE_ERROR;
00513                 #ifdef TPRDEBUG
00514                 printf("posresXA: %d\n", k);
00515                 #endif
00516                 //Position restraints have a float for every coordinate, hence the multiplier by 3.
00517                 for (j = 0; j < 3 * k; j++) {
00518                         if (trx_real(mf, NULL)) return MOLFILE_ERROR;
00519                 }
00520                 if (trx_int(mf, &k)) return MOLFILE_ERROR;
00521                 #ifdef TPRDEBUG
00522                 printf("posresXB: %d\n", k);
00523                 #endif
00524                 for (j = 0; j < 3 * k; j++) {
00525                         if (trx_real(mf, NULL)) return MOLFILE_ERROR;
00526                 }
00527                 #ifdef TPRDEBUG
00528                 printf("Segname: %d ", tpr->molbtype[i]);
00529                 for (j = 0; j < SAVELEN; j++) {
00530                         printf("%c", tpr->symtab[SAVELEN*tpr->molnames[tpr->molbtype[i]]+j]);
00531                 }
00532                 printf("\n\tNumatoms: %d\n", tpr->molbnatoms[i]);
00533                 printf("\tNumcopies: %d\n", tpr->molbnmol[i]);
00534                 #endif
00535                 //sum += tpr->molbnatoms[i] * tpr->molbnmol[i];
00536         }
00537         //Burn off the number of atoms after the do_molblock
00538         if (trx_int(mf, &tmp)) return MOLFILE_ERROR;
00539 #ifdef TPRDEBUG
00540     printf("What is this? %d. It should be the number of atoms.\n", tmp);
00541 #endif
00542     if (tpr->version >= 103) { //103 is tpxv_IntermolecularBondeds
00543         if (trx_int(mf, &hasIntermoleculeBonds)) return MOLFILE_ERROR;
00544 #ifdef TPRDEBUG
00545         printf("intermolecularbondeds %d\n", hasIntermoleculeBonds);
00546 #endif
00547         if (hasIntermoleculeBonds) {
00548             printf("Systems with intermolecular bonds are not supported. The file reports %d intermolecular bonds.\n", hasIntermoleculeBonds);
00549             return MOLFILE_ERROR;
00550         }
00551         if (tpr->wversion >= 27) {
00552                 fseek(mf->f, -3, SEEK_CUR);
00553         }
00554     }
00555 
00556     if (tpr->version < 128) //128 is tpxv_RemoveAtomtypes
00557     {
00558             //do_atomtypes
00559             if (read_atomtypes(mf, tpr->version)) return MOLFILE_ERROR;
00560         }
00561 
00562 #ifdef TPRDEBUG
00563     printf("Reading cmap terms\n");
00564 #endif
00565     //do_cmap
00566     if (tpr->version >= 65) {
00567         if (read_cmap(mf)) return MOLFILE_ERROR;
00568     }
00569     //do_groups
00570 #ifdef TPRDEBUG
00571     printf("Reading groups\n");
00572 #endif
00573     read_groups(mf, egcNR, tpr);
00574     if (tpr->version >= 120) {
00575         long len;
00576         if (trx_long(mf, &len)) return MOLFILE_ERROR;
00577         int* jj = new int[len];
00578         #ifdef TPRDEBUG
00579         printf("Intermolecular Exclusions: %d\n", len);
00580         printf("%d\n", ftell(mf->f));
00581         #endif
00582         tpr_ivector(mf, jj, len);
00583     }
00584 #ifdef TPRDEBUG
00585     printf("This is my file position: %d\n", ftell(mf->f));
00586     printf("Returning control\n");
00587 #endif
00588     return MOLFILE_SUCCESS;
00589 }
00590 
00591 static int read_tpr_structure(void *mydata, int *optflags, molfile_atom_t *atoms) {
00592     tprdata *tpr = (tprdata *) mydata;
00593     //printf("Structure stuff\n");
00594     *optflags = MOLFILE_MASS | MOLFILE_CHARGE; // | MOLFILE_ATOMICNUMBER; //For the coarse grained residues I'm testing, the atomic number is -1. That is bad.
00595     int idx = 0;
00596     int i, j, k;
00597     for (i = 0; i < tpr->nmolblock; i++) {
00598         for (j = 0; j < tpr->molbnmol[i]; j++) {
00599             for (k = 0; k < tpr->molbnatoms[i]; k++) {
00600                 //printf("%d %d %d\n", i, j, k);
00601                 //printf("%s\n",&(tpr->symtab[SAVELEN*tpr->molnames[tpr->molbtype[i]]]));
00602                 strcpy(atoms[idx].segid, &(tpr->symtab[SAVELEN*tpr->molnames[tpr->molbtype[i]]]));
00603                 //printf("%s\n",&(tpr->symtab[SAVELEN*tpr->atomnameids[tpr->molbtype[i]][k]]));
00604                 strcpy(atoms[idx].name, &(tpr->symtab[SAVELEN*tpr->atomnameids[tpr->molbtype[i]][k]]));
00605                 //printf("%s\n",&(tpr->symtab[SAVELEN*tpr->atomtypeids[tpr->molbtype[i]][k]]));
00606                 strcpy(atoms[idx].type, &(tpr->symtab[SAVELEN*tpr->atomtypeids[tpr->molbtype[i]][k]]));
00607                 //printf("%d\n",tpr->resids[tpr->molbtype[i]][k]);
00608                 atoms[idx].resid = tpr->resids[tpr->molbtype[i]][k];
00609                 //printf("%s\n",&(tpr->symtab[SAVELEN*tpr->resnames[tpr->molbtype[i]][atoms[idx].resid]]));
00610                 strcpy(atoms[idx].resname, &(tpr->symtab[SAVELEN*tpr->resnames[tpr->molbtype[i]][atoms[idx].resid]]));
00611                 //printf("%f\n",tpr->masses[tpr->molbtype[i]][k]);
00612                 atoms[idx].mass = tpr->masses[tpr->molbtype[i]][k];
00613                 //printf("%f\n",tpr->charges[tpr->molbtype[i]][k]);
00614                 atoms[idx].charge = tpr->charges[tpr->molbtype[i]][k];
00615                 //atoms[idx].atomicnumber = tpr->atomicnumbers[tpr->molbtype[i]][k];
00616                 idx++;
00617             }
00618         }
00619     }
00620     //printf("I'm done here\n");
00621     return MOLFILE_SUCCESS;
00622 }
00623 
00624 static void *open_tpr_read(const char *filename, const char *,
00625     int *natoms) {
00626     tprdata *tprdat = NULL;
00627     FILE *fin;
00628     char buf[STRLEN];
00629     md_file *mf = new md_file;
00630 
00631     int i, j;
00632     if (!(fin = fopen(filename, "rb"))) {
00633         fprintf(stderr, "tprplugin) Cannot open tpr file '%s'\n", filename);
00634         return NULL;
00635     }
00636     mf->f = fin;
00637     if (trx_int(mf, &i)) {
00638         fprintf(stderr, "tprplugin) Could not read initial integer from file.\n");
00639         return NULL;
00640     }
00641         if (i > STRLEN) {//If i value is large, everything in the file should be endian swapped.
00642                 mf->rev = 1;
00643         }
00644         j = tpr_string(mf, buf, STRLEN);
00645         if (trx_int(mf, &(mf->prec))) {
00646                 fprintf(stderr, "tprplugin) Could not read precision from file.\n");
00647         return NULL;
00648         }
00649     if (mf->prec == 4) {
00650         tprdat = new tprdata;
00651         memset(tprdat, 0, sizeof(tprdata));
00652         tprdat->mf = mf;
00653         if (readtprAfterPrecision(tprdat) != MOLFILE_SUCCESS) {
00654             delete tprdat;
00655             return NULL;
00656         }
00657         *natoms = tprdat->natoms;
00658         return tprdat;
00659     }
00660     else {
00661         fprintf(stderr, "tprplugin) Illegal precision (requires single)\n");
00662         return NULL;
00663     }
00664     fclose(fin);
00665     return NULL;
00666 }
00667 
00668 static int read_tpr_timestep(void *v, int natoms, molfile_timestep_t *ts) {
00669     tprdata *tpr = (tprdata *)v;
00670     md_file *mf = tpr->mf;
00671     if (tpr->readcoordinates) {
00672         return MOLFILE_ERROR;
00673     }
00674     //Get the positions.
00675     if (ts != NULL) {
00676         tpr_rvector(mf, ts->coords, 3 * tpr->natoms);
00677         for (int i = 0; i < 3 * natoms; i++) {
00678             ts->coords[i] = 10 * ts->coords[i]; //A
00679         }
00680         #ifdef TPRDEBUG
00681         printf("\nAtoms Coordinates: (A)\n");
00682         for (int i = 0; i < natoms; i++) {
00683             printf("x[%d]: %f %f %f\n", i,   ts->coords[3*i+0], ts->coords[3*i+1],ts->coords[3*i+2]);
00684         }
00685         printf("coordinate end position: %d\n", ftell(mf->f));
00686 
00687         for (int i = 0; i < 3*natoms; i++) {
00688             if (! isfinite(ts->coords[i])) {
00689                 printf("The %d coordinate of atom %d is not finite!\n", (i % 3), i / 3);
00690             }
00691         }
00692         #endif
00693 
00694         if(tpr->has_velocities)
00695         {
00696             tpr_rvector(mf, ts->velocities, 3 * tpr->natoms);
00697             for (int i = 0; i < 3 * natoms; i++)
00698             {
00699                 ts->velocities[i] *= 10; //A/ps
00700                 //fprintf(stderr, "%f\n", ts->velocities[i]);
00701             }
00702         }
00703 
00704         #ifdef TPRDEBUG
00705         if(tpr->has_velocities)
00706         {
00707             printf("\nAtoms Velocities: (A/ps)\n");
00708             for (int i = 0; i < natoms; i++)
00709             {
00710                 printf("v[%d]: %f %f %f\n", i,   ts->velocities[3*i+0], ts->velocities[3*i+1],ts->velocities[3*i+2]);
00711             }
00712             printf("velocity end position: %d\n", ftell(mf->f));
00713         }
00714         #endif
00715 
00716         ts->A = sqrt(tpr->boxdims[0]*tpr->boxdims[0] + tpr->boxdims[1]*tpr->boxdims[1] + tpr->boxdims[2] * tpr->boxdims[2]) * 10;
00717         ts->B = sqrt(tpr->boxdims[3]*tpr->boxdims[3] + tpr->boxdims[4]*tpr->boxdims[4] + tpr->boxdims[5] * tpr->boxdims[5]) * 10;
00718         ts->C = sqrt(tpr->boxdims[6]*tpr->boxdims[6] + tpr->boxdims[7]*tpr->boxdims[7] + tpr->boxdims[8] * tpr->boxdims[8]) * 10;
00719 
00720         if(ts->A <= 0 || ts->B <= 0 || ts->C <=0)
00721         {
00722             ts->A     = ts->B    = ts->C     = 0;
00723             ts->alpha = ts->beta = ts->gamma = 0;
00724         }
00725         else
00726         {
00727                 ts->alpha = acos((tpr->boxdims[3] * tpr->boxdims[6] + tpr->boxdims[4] * tpr->boxdims[7] + tpr->boxdims[5] * tpr->boxdims[8])*100/(ts->A*ts->B)) * 90.0 / M_PI_2;
00728                 ts->beta  = acos((tpr->boxdims[0] * tpr->boxdims[6] + tpr->boxdims[1] * tpr->boxdims[7] + tpr->boxdims[2] * tpr->boxdims[8])*100/(ts->A*ts->C)) * 90.0 / M_PI_2;
00729                 ts->gamma = acos((tpr->boxdims[0] * tpr->boxdims[3] + tpr->boxdims[1] * tpr->boxdims[4] + tpr->boxdims[2] * tpr->boxdims[5])*100/(ts->B*ts->C)) * 90.0 / M_PI_2;
00730         }
00731         //printf("%f %f %f %f %f %f\n", ts->A, ts->B, ts->C, ts->alpha, ts->beta, ts->gamma);
00732     }
00733     tpr->readcoordinates = 1;
00734     return MOLFILE_SUCCESS;
00735 }
00736 static int read_tpr_bonds(void *v, int *nbonds, int **fromptr, int **toptr,
00737                          float **bondorder, int **bondtype,
00738                          int *nbondtypes, char ***bondtypename) {
00739     tprdata *tpr = (tprdata *)v;
00740     int i, j, k, l, itraction, mtype;
00741     const int bondinteractions = 4;
00742     int bondinteraction[bondinteractions] = {F_BONDS,F_G96BONDS,F_CONSTR,F_SETTLE};
00743     int aoffset = 0;
00744     int boffset = 0;
00745     *nbonds = 0;
00746     *fromptr = NULL;
00747     *toptr = NULL;
00748     *bondorder = NULL;
00749     *bondtype = NULL;
00750     *nbondtypes = 0;
00751     *bondtypename = NULL;
00752     //Gromacs made this complicated. I'm doing to try and do this the simple way, checking one at a time to see if the interaction types aren't empty.
00753     for (i = 0; i < tpr->nmolblock; i++) {
00754         mtype = tpr->molbtype[i];
00755         //printf("%d atoms in block %d\n", tpr->atomsinmol[i], i);
00756         //printf("Type %d\n", tpr->molbtype[i]);
00757         for (j = 0; j < tpr->molbnmol[i]; j++) {
00758             for (k = 0; k < bondinteractions; k++) {
00759                 itraction = bondinteraction[k];
00760                 if (k == 3) {
00761                     //printf("%d elements in settle\n", tpr->nr[itraction][mtype]);
00762                     *nbonds += 2 * tpr->nr[itraction][mtype] / 4;
00763                 } else {
00764                     *nbonds += tpr->nr[itraction][mtype] / 3;
00765                 }
00766                 //printf("nbonds = %d\n", *nbonds);
00767             }
00768         }
00769     }
00770     int *fromlist = (int *) malloc(*nbonds * sizeof(int));
00771     int *tolist = (int *) malloc(*nbonds * sizeof(int));
00772     int *bolist = (int *) malloc(*nbonds * sizeof(int));
00773     //printf("I Malloced\n");
00774     for (i = 0; i < tpr->nmolblock; i++) {
00775         mtype = tpr->molbtype[i];
00776         for (j = 0; j < tpr->molbnmol[i]; j++) {
00777             for (k = 0; k < bondinteractions; k++) {
00778                 itraction = bondinteraction[k];
00779                 if (k==3) {
00780                     for (l = 0; l < 2 * tpr->nr[itraction][mtype] / 4; l++) {
00781                         fromlist[boffset + l] = 1 + tpr->interactionlist[itraction][mtype][4*(l/2)+1] + aoffset;
00782                         tolist[boffset + l] = 1 + tpr->interactionlist[itraction][mtype][4*(l/2)+2+l] + aoffset;
00783                         bolist[boffset + l] = k;
00784                     }
00785                     boffset += 2 * tpr->nr[itraction][mtype] / 4;
00786                     /*for (l = 0; l < tpr->nr[itraction][mtype]; l++) {
00787                         printf("%d ", tpr->interactionlist[itraction][mtype][l]);
00788                     }
00789                     printf("\n");*/
00790                 }
00791                 else {
00792                     for (l = 0; l < tpr->nr[itraction][mtype] / 3; l++) {
00793                         //The 1+ comes because the from and to lists are 1-indexed, since psfs are 1 indexed.
00794                         fromlist[boffset + l] = 1 + tpr->interactionlist[itraction][mtype][3*l+1] + aoffset;
00795                         tolist[boffset + l] = 1 + tpr->interactionlist[itraction][mtype][3*l+2] + aoffset;
00796                         bolist[boffset + l] = k;
00797                     }
00798                     boffset += tpr->nr[itraction][mtype] / 3;
00799                 }
00800             }
00801             aoffset += tpr->atomsinmol[mtype];
00802         }
00803     }
00804     /*printf("Bondlist:\n");
00805     for (i = 0; i < *nbonds; i++) {
00806         printf("%d to %d\n", fromlist[i], tolist[i]);
00807     }*/
00808     *fromptr = fromlist;
00809     *toptr = tolist;
00810     return MOLFILE_SUCCESS;
00811 }
00812 
00813 static int read_tpr_angles(void *v, int *numangles, int **angles,
00814                          int **angletypes, int *numangletypes,
00815                          char ***angletypenames, int *numdihedrals,
00816                          int **dihedrals, int **dihedraltypes,
00817                          int *numdihedraltypes, char ***dihedraltypenames,
00818                          int *numimpropers, int **impropers,
00819                          int **impropertypes, int *numimpropertypes,
00820                          char ***impropertypenames, int *numcterms,
00821                          int **cterms, int *ctermcols, int *ctermrows) {
00822     tprdata *tpr = (tprdata *)v;
00823     int i, j, k, l, itraction, mtype;
00824     int aoffset = 0;
00825     int boffset = 0;
00826     /* initialize data to zero */
00827     *numangles         = 0;
00828     *angles            = NULL;
00829     *angletypes        = NULL;
00830     *numangletypes     = 0;
00831     *angletypenames    = NULL;
00832     *numdihedrals      = 0;
00833     *dihedrals         = NULL;
00834     *dihedraltypes     = NULL;
00835     *numdihedraltypes  = 0;
00836     *dihedraltypenames = NULL;
00837     *numimpropers      = 0;
00838     *impropers         = NULL;
00839     *impropertypes     = NULL;
00840     *numimpropertypes  = 0;
00841     *impropertypenames = NULL;
00842     *numcterms         = 0;
00843     *cterms            = NULL;
00844     *ctermrows         = 0;
00845     *ctermcols         = 0;
00846     const int angleinteractions = 6;
00847     int angleinteraction[angleinteractions] = {F_ANGLES, F_G96ANGLES, F_RESTRANGLES, F_LINEAR_ANGLES, F_UREY_BRADLEY, F_SETTLE};
00848     for (i = 0; i < tpr->nmolblock; i++) {
00849         mtype = tpr->molbtype[i];
00850         for (j = 0; j < tpr->molbnmol[i]; j++) {
00851             for (k = 0; k < angleinteractions; k++) {
00852                 itraction = angleinteraction[k];
00853                 *numangles += tpr->nr[itraction][mtype] / 4;
00854             }
00855         }
00856     }
00857     //printf("Numangles read in: %d\n", *numangles);
00858     int *anglelist = (int *) malloc(3 * (*numangles) * sizeof(int));
00859     //printf("I Malloced\n");
00860     for (i = 0; i < tpr->nmolblock; i++) {
00861         mtype = tpr->molbtype[i];
00862         //printf("%d atoms in block %d\n", tpr->atomsinmol[i], i);
00863         for (j = 0; j < tpr->molbnmol[i]; j++) {
00864             for (k = 0; k < angleinteractions; k++) {
00865                 itraction = angleinteraction[k];
00866                 for (l = 0; l < tpr->nr[itraction][mtype] / 4; l++) {
00867                     //The 1+ comes because the angle lists are 1-indexed, since psfs are 1 indexed.
00868                     if (k == 5) {
00869                         anglelist[boffset + 3*l + 0] = 1 + tpr->interactionlist[itraction][mtype][4*l+2] + aoffset;
00870                         anglelist[boffset + 3*l + 1] = 1 + tpr->interactionlist[itraction][mtype][4*l+1] + aoffset;
00871                         anglelist[boffset + 3*l + 2] = 1 + tpr->interactionlist[itraction][mtype][4*l+3] + aoffset;
00872                     }
00873                     else {
00874                         anglelist[boffset + 3*l + 0] = 1 + tpr->interactionlist[itraction][mtype][4*l+1] + aoffset;
00875                         anglelist[boffset + 3*l + 1] = 1 + tpr->interactionlist[itraction][mtype][4*l+2] + aoffset;
00876                         anglelist[boffset + 3*l + 2] = 1 + tpr->interactionlist[itraction][mtype][4*l+3] + aoffset;
00877                     }
00878                 }
00879                 boffset += 3*(tpr->nr[itraction][mtype] / 4);
00880             }
00881             aoffset += tpr->atomsinmol[mtype];
00882         }
00883     }
00884     *angles = anglelist;
00885 
00886     aoffset = 0;
00887     boffset = 0;
00888     const int pdihsinteractions = 1;
00889     int pdihsinteraction[pdihsinteractions] = {F_PDIHS};
00890     for (i = 0; i < tpr->nmolblock; i++) {
00891         mtype = tpr->molbtype[i];
00892         for (j = 0; j < tpr->molbnmol[i]; j++) {
00893             for (k = 0; k < pdihsinteractions; k++) {
00894                 itraction = pdihsinteraction[k];
00895                 *numdihedrals += tpr->nr[itraction][mtype] / 5;
00896             }
00897         }
00898     }
00899     //printf("Numdihedrals read in: %d\n", *numdihedrals);
00900     int *pdihlist = (int *) malloc(4 * (*numdihedrals) * sizeof(int));
00901     //printf("I Malloced\n");
00902     for (i = 0; i < tpr->nmolblock; i++) {
00903         mtype = tpr->molbtype[i];
00904         //printf("%d atoms in block %d\n", tpr->atomsinmol[i], i);
00905         for (j = 0; j < tpr->molbnmol[i]; j++) {
00906             for (k = 0; k < pdihsinteractions; k++) {
00907                 itraction = pdihsinteraction[k];
00908                 for (l = 0; l < tpr->nr[itraction][mtype] / 5; l++) {
00909                     //The 1+ comes because the angle lists are 1-indexed, since psfs are 1 indexed.
00910                     pdihlist[boffset + 4*l + 0] = 1 + tpr->interactionlist[itraction][mtype][5*l+1] + aoffset;
00911                     pdihlist[boffset + 4*l + 1] = 1 + tpr->interactionlist[itraction][mtype][5*l+2] + aoffset;
00912                     pdihlist[boffset + 4*l + 2] = 1 + tpr->interactionlist[itraction][mtype][5*l+3] + aoffset;
00913                     pdihlist[boffset + 4*l + 3] = 1 + tpr->interactionlist[itraction][mtype][5*l+4] + aoffset;
00914                 }
00915                 boffset += 4*(tpr->nr[itraction][mtype] / 5);
00916             }
00917             aoffset += tpr->atomsinmol[mtype];
00918         }
00919     }
00920     *dihedrals = pdihlist;
00921 
00922     aoffset = 0;
00923     boffset = 0;
00924     const int idihsinteractions = 1;
00925     int idihsinteraction[idihsinteractions] = {F_IDIHS};
00926     for (i = 0; i < tpr->nmolblock; i++) {
00927         mtype = tpr->molbtype[i];
00928         for (j = 0; j < tpr->molbnmol[i]; j++) {
00929             for (k = 0; k < idihsinteractions; k++) {
00930                 itraction = idihsinteraction[k];
00931                 *numimpropers += tpr->nr[itraction][mtype] / 5;
00932             }
00933         }
00934     }
00935     //printf("Numimpropers read in: %d\n", *numimpropers);
00936     int *idihlist = (int *) malloc(4 * (*numimpropers) * sizeof(int));
00937     //printf("I Malloced\n");
00938     for (i = 0; i < tpr->nmolblock; i++) {
00939         mtype = tpr->molbtype[i];
00940         //printf("%d atoms in block %d\n", tpr->atomsinmol[i], i);
00941         for (j = 0; j < tpr->molbnmol[i]; j++) {
00942             for (k = 0; k < idihsinteractions; k++) {
00943                 itraction = idihsinteraction[k];
00944                 for (l = 0; l < tpr->nr[itraction][mtype] / 5; l++) {
00945                     //The 1+ comes because the angle lists are 1-indexed, since psfs are 1 indexed.
00946                     idihlist[boffset + 4*l + 0] = 1 + tpr->interactionlist[itraction][mtype][5*l+1] + aoffset;
00947                     idihlist[boffset + 4*l + 1] = 1 + tpr->interactionlist[itraction][mtype][5*l+2] + aoffset;
00948                     idihlist[boffset + 4*l + 2] = 1 + tpr->interactionlist[itraction][mtype][5*l+3] + aoffset;
00949                     idihlist[boffset + 4*l + 3] = 1 + tpr->interactionlist[itraction][mtype][5*l+4] + aoffset;
00950                 }
00951                 boffset += 4*(tpr->nr[itraction][mtype] / 5);
00952             }
00953             aoffset += tpr->atomsinmol[mtype];
00954         }
00955     }
00956     *impropers = idihlist;
00957 
00958     aoffset = 0;
00959     boffset = 0;
00960     const int cmapinteraction = 1;
00961     int cmapinteractions[cmapinteraction] = {F_CMAP};
00962     for (i = 0; i < tpr->nmolblock; i++) {
00963         mtype = tpr->molbtype[i];
00964         for (j = 0; j < tpr->molbnmol[i]; j++) {
00965             for (k = 0; k < cmapinteraction; k++) {
00966                 itraction = cmapinteractions[k];
00967                 *numcterms += tpr->nr[itraction][mtype] / 6;
00968             }
00969         }
00970     }
00971     //printf("Numcmap terms read in: %d\n", *numcterms);
00972     int *cmaplist = (int *) malloc(8 * (*numcterms) * sizeof(int));
00973     //printf("I Malloced\n");
00974     for (i = 0; i < tpr->nmolblock; i++) {
00975         mtype = tpr->molbtype[i];
00976         //printf("%d atoms in block %d\n", tpr->atomsinmol[i], i);
00977         for (j = 0; j < tpr->molbnmol[i]; j++) {
00978             for (k = 0; k < cmapinteraction; k++) {
00979                 itraction = cmapinteractions[k];
00980                 for (l = 0; l < tpr->nr[itraction][mtype] / 6; l++) {
00981                     //The 1+ comes because the angle lists are 1-indexed, since psfs are 1 indexed.
00982                     cmaplist[boffset + 8*l + 0] = 1 + tpr->interactionlist[itraction][mtype][6*l+1] + aoffset;
00983                     cmaplist[boffset + 8*l + 1] = 1 + tpr->interactionlist[itraction][mtype][6*l+2] + aoffset;
00984                     cmaplist[boffset + 8*l + 2] = 1 + tpr->interactionlist[itraction][mtype][6*l+3] + aoffset;
00985                     cmaplist[boffset + 8*l + 3] = 1 + tpr->interactionlist[itraction][mtype][6*l+4] + aoffset;
00986                     cmaplist[boffset + 8*l + 4] = 1 + tpr->interactionlist[itraction][mtype][6*l+2] + aoffset;
00987                     cmaplist[boffset + 8*l + 5] = 1 + tpr->interactionlist[itraction][mtype][6*l+3] + aoffset;
00988                     cmaplist[boffset + 8*l + 6] = 1 + tpr->interactionlist[itraction][mtype][6*l+4] + aoffset;
00989                     cmaplist[boffset + 8*l + 7] = 1 + tpr->interactionlist[itraction][mtype][6*l+5] + aoffset;
00990                 }
00991                 boffset += 8*(tpr->nr[itraction][mtype] / 6);
00992             }
00993             aoffset += tpr->atomsinmol[mtype];
00994         }
00995     }
00996     *cterms = cmaplist;
00997 
00998     *ctermcols = 0;
00999     *ctermrows = 0;
01000     return MOLFILE_SUCCESS;
01001 }
01002 
01003 static int read_tpr_timestep_metadata(void *v, molfile_timestep_metadata_t *metadata)
01004 {
01005     tprdata *tpr = (tprdata *)v;
01006     if(tpr->has_velocities == 1)
01007     {
01008         metadata->has_velocities = 1;
01009     }
01010     else
01011     {
01012         metadata->has_velocities = 0;
01013     }
01014     return MOLFILE_SUCCESS;
01015 }
01016 
01017 static void close_tpr_read(void *mydata) {
01018         int i, j;
01019         //printf("Closing TPR file\n");
01020         tprdata *tpr = (tprdata *)mydata;
01021         fclose(tpr->mf->f);
01022         for (i = 0; i < tpr->nmoltypes; i++) {
01023                 //printf("%d of %d (%d atoms)\n", i+1, tpr->nmoltypes, tpr->atomsinmol[i]);
01024                 free(tpr->charges[i]);
01025                 free(tpr->masses[i]);
01026                 free(tpr->types[i]);
01027                 free(tpr->ptypes[i]);
01028                 free(tpr->resids[i]);
01029                 free(tpr->atomnameids[i]);
01030                 free(tpr->atomtypeids[i]);
01031                 free(tpr->resnames[i]);
01032                 free(tpr->atomicnumbers[i]);
01033         }
01034         for (i = 0; i < F_NRE; i++) {
01035                 if (tpr->interactionlist[i] != NULL) {
01036                         for (j = 0; j < tpr->nmoltypes; j++) {
01037                                 if (tpr->interactionlist[i][j] != NULL)
01038                                         free(tpr->interactionlist[i][j]);
01039                         }
01040                         free(tpr->interactionlist[i]);
01041                 }
01042                 if (tpr->nr[i] != NULL) {
01043                         free(tpr->nr[i]);
01044                 }
01045         }
01046         free(tpr->atomicnumbers);
01047         free(tpr->molnames);
01048         free(tpr->molbtype);
01049         free(tpr->molbnmol);
01050         free(tpr->molbnatoms);
01051         free(tpr->resnames);
01052         free(tpr->atomtypeids);
01053         free(tpr->atomnameids);
01054         free(tpr->types);
01055         free(tpr->ptypes);
01056         free(tpr->resids);
01057         free(tpr->charges);
01058         free(tpr->masses);
01059         free(tpr->symtab);
01060         free(tpr->atomsinmol);
01061         free(tpr->resinmol);
01062         delete tpr->mf;
01063         delete tpr;
01064         //printf("TPR file read completely\n");
01065 }
01066 
01067 
01068 /*
01069  * Initialization stuff down here
01070  */
01071 
01072 #ifndef TPRTEST
01073 
01074 static molfile_plugin_t tpr_plugin;
01075 
01076 VMDPLUGIN_API int VMDPLUGIN_init(void) {
01077     memset(&tpr_plugin, 0, sizeof(molfile_plugin_t));
01078     tpr_plugin.abiversion = vmdplugin_ABIVERSION;
01079     tpr_plugin.type = MOLFILE_PLUGIN_TYPE;
01080     tpr_plugin.name = "tpr";
01081     tpr_plugin.prettyname = "Gromacs Binary Topology";
01082     tpr_plugin.author = "Josh Vermaas";
01083     tpr_plugin.majorv = 2023;
01084     tpr_plugin.minorv = 0;//Corresponds to the Gromacs version I was basing this on.
01085     tpr_plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
01086     tpr_plugin.filename_extension = "tpr";
01087     tpr_plugin.open_file_read = open_tpr_read;
01088     tpr_plugin.read_timestep_metadata = read_tpr_timestep_metadata;
01089     tpr_plugin.read_structure = read_tpr_structure;
01090     tpr_plugin.read_next_timestep = read_tpr_timestep;
01091     tpr_plugin.read_bonds = read_tpr_bonds;
01092     tpr_plugin.read_angles = read_tpr_angles;
01093     tpr_plugin.close_file_read = close_tpr_read;
01094     return VMDPLUGIN_SUCCESS;
01095 }
01096 
01097 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
01098     (*cb)(v, (vmdplugin_t *)&tpr_plugin);
01099     return 0;
01100 }
01101 
01102 VMDPLUGIN_API int VMDPLUGIN_fini(void) { return VMDPLUGIN_SUCCESS; }
01103 
01104 #else
01105 
01106 int main (int argc, char *argv[]) {
01107         int natoms;
01108         tprdata *tprdat = NULL;
01109         for (int fcount=1; fcount < argc; fcount++) {
01110                 printf("Attempting to read %s\n", argv[fcount]);
01111                 tprdat = (tprdata*)open_tpr_read(argv[fcount], NULL, &natoms);
01112                 if (tprdat != NULL){
01113                 printf("Total number of atoms: %d, %d\n", natoms, tprdat->natoms);
01114                         printf("Finished initial reading\n");
01115                         molfile_timestep_t *ts = new molfile_timestep_t;
01116                         ts->coords = new float[3*tprdat->natoms];
01117                         ts->velocities = new float[3*tprdat->natoms];
01118                         read_tpr_timestep(tprdat, tprdat->natoms, ts);
01119             }
01120             else {
01121                 fprintf(stderr, "open_tpr_read failed\n");
01122                 return NULL;
01123             }
01124             FILE *fin = tprdat->mf->f;
01125                 fseek(fin, 0L, SEEK_END);
01126                 long length = ftell(fin);
01127         printf("END : %ld\n", length);
01128                 fclose(fin);
01129         }
01130         return 0;
01131 }
01132 
01133 #endif
01134 

Generated on Mon Oct 14 03:09:44 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002