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

MolFilePlugin.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: MolFilePlugin.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.198 $      $Date: 2020/10/22 03:43:23 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   VMD interface to 'molfile' plugins.  Molfile plugins read coordinate
00019  *   files, structure files, volumetric data, and graphics data.  The data
00020  *   is loaded into a new or potentially preexisting molecule in VMD.
00021  *   Some molefile plugins can also export data from VMD to files.
00022  * 
00023  * LICENSE:
00024  *   UIUC Open Source License
00025  *   http://www.ks.uiuc.edu/Research/vmd/plugins/pluginlicense.html
00026  * 
00027  ***************************************************************************/
00028 
00029 #include <stdio.h>
00030 #include <stddef.h>
00031 
00032 #include "MolFilePlugin.h"
00033 #include "Molecule.h"
00034 #include "AtomSel.h"
00035 #include "Timestep.h"
00036 #include "Inform.h"
00037 #include "Scene.h"
00038 #include "molfile_plugin.h"
00039 #include "PeriodicTable.h"
00040 #include "VolumetricData.h"
00041 #include "MoleculeGraphics.h"
00042 #include "QMData.h"
00043 #include "QMTimestep.h"
00044 
00045 #if defined(VMDTKCON)
00046 #include "vmdconsole.h"
00047 #endif
00048 
00049 #define MINSIZEOF(a, b) \
00050         ((sizeof(a) < sizeof(b)) ? sizeof(a) : sizeof(b))
00051 
00052 #define SAFESTRNCPY(a, b) \
00053         (strncpy(a, b, MINSIZEOF(a, b)))
00054 
00055 MolFilePlugin::MolFilePlugin(vmdplugin_t *p)
00056 : plugin((molfile_plugin_t *)p) {
00057   rv = wv = NULL; 
00058   numatoms = 0;
00059   _filename = NULL;
00060   qm_data = NULL;
00061 #if defined(VMDTKCON)
00062    plugin->cons_fputs = &vmdcon_fputs;
00063 #else
00064    plugin->cons_fputs = NULL;
00065 #endif
00066 }
00067 
00068 MolFilePlugin::~MolFilePlugin() {
00069   close();
00070   delete [] _filename;
00071 }
00072 
00073 int MolFilePlugin::read_structure(Molecule *m, int filebonds, int autobonds) {
00074   if (!rv) return MOLFILE_ERROR; 
00075   if (!can_read_structure()) return MOLFILE_ERROR;
00076   if (!m->init_atoms(numatoms)) return MOLFILE_ERROR;
00077 
00078   molfile_atom_t *atoms = 
00079     (molfile_atom_t *) calloc(1, numatoms*sizeof(molfile_atom_t));
00080 
00081   int rc, i;
00082   int optflags = MOLFILE_BADOPTIONS; /* plugin must reset this correctly */
00083   if ((rc = plugin->read_structure(rv, &optflags, atoms))) {
00084     free(atoms);
00085     return rc; // propagate error to caller
00086   }
00087   if (optflags == int(MOLFILE_BADOPTIONS)) {
00088     free(atoms);
00089     msgErr << "MolFilePlugin: plugin didn't initialize optional data flags" << sendmsg;
00090     msgErr << "MolFilePlugin: file load aborted" << sendmsg;
00091     return MOLFILE_ERROR; /* abort load, data can't be trusted */
00092   }
00093 
00094 
00095   float *charge = m->charge();
00096   float *mass = m->mass();
00097   float *radius = m->radius();
00098   float *beta = m->beta();
00099   float *occupancy = m->occupancy();
00100 
00101   // set molecule dataset flags
00102   if (optflags & MOLFILE_INSERTION)
00103     m->set_dataset_flag(BaseMolecule::INSERTION);
00104   if (optflags & MOLFILE_OCCUPANCY)
00105     m->set_dataset_flag(BaseMolecule::OCCUPANCY);
00106   if (optflags & MOLFILE_BFACTOR)
00107     m->set_dataset_flag(BaseMolecule::BFACTOR);
00108   if (optflags & MOLFILE_MASS)
00109     m->set_dataset_flag(BaseMolecule::MASS);
00110   if (optflags & MOLFILE_CHARGE)
00111     m->set_dataset_flag(BaseMolecule::CHARGE);
00112   if (optflags & MOLFILE_RADIUS)
00113     m->set_dataset_flag(BaseMolecule::RADIUS);
00114   if (optflags & MOLFILE_ALTLOC)
00115     m->set_dataset_flag(BaseMolecule::ALTLOC);
00116   if (optflags & MOLFILE_ATOMICNUMBER)
00117     m->set_dataset_flag(BaseMolecule::ATOMICNUMBER);
00118 
00119   // Force min/max atom radii to be updated since we've loaded new data
00120   if (optflags & MOLFILE_RADIUS)
00121     m->set_radii_changed();
00122 
00123   for (i=0; i<numatoms; i++) {
00124     molfile_atom_t &atom = atoms[i];
00125     int atomicnumber = (optflags & MOLFILE_ATOMICNUMBER) ? 
00126       atom.atomicnumber : -1;
00127     charge[i] = (optflags & MOLFILE_CHARGE) ?
00128       atom.charge : m->default_charge(atom.name);
00129  
00130     // If we're given an explicit mass value, use it.
00131     // Failing that, try doing a periodic table lookup if we have
00132     // a valid atomicnumber.  If that fails also, then we use a crude
00133     // guess based on the atom name string. 
00134     mass[i] = (optflags & MOLFILE_MASS) ?
00135       atom.mass : ((atomicnumber > 0) ?
00136         get_pte_mass(atomicnumber) : m->default_mass(atom.name));
00137 
00138     // If we're given an explicit VDW radius value, use it.
00139     // Failing that, try doing a periodic table lookup if we have
00140     // a valid atomicnumber.  If that fails also, then we use a crude
00141     // guess based on the atom name string. 
00142     radius[i] = (optflags & MOLFILE_RADIUS) ?
00143       atom.radius : ((atomicnumber > 0) ?
00144         get_pte_vdw_radius(atomicnumber) : m->default_radius(atom.name));
00145 
00146     beta[i] = (optflags & MOLFILE_BFACTOR) ?
00147       atom.bfactor : m->default_beta();
00148     occupancy[i] = (optflags & MOLFILE_OCCUPANCY) ?
00149       atom.occupancy : m->default_occup();
00150     const char *insertion = (optflags & MOLFILE_INSERTION) ? 
00151       atom.insertion : " ";
00152     const char *altloc = (optflags & MOLFILE_ALTLOC) ?
00153       atom.altloc : "";
00154     if (0 > m->add_atoms(1,
00155                          atom.name, atom.type, atomicnumber, 
00156                          atom.resname, atom.resid, 
00157                          atom.chain, atom.segid, (char *)insertion, altloc)) {
00158       // if an error occured while adding an atom, we should delete
00159       // the offending molecule since the data is presumably inconsistent,
00160       // or at least not representative of what we tried to load
00161       msgErr << "MolFilePlugin: file load aborted" << sendmsg;
00162       free(atoms);
00163       return MOLFILE_ERROR; // abort load, data can't be trusted
00164     }
00165   }
00166   free(atoms);
00167 
00168   if (filebonds && can_read_bonds()) {
00169     int nbonds, *from, *to, nbondtypes, *bondtype;
00170     float *bondorder;
00171     char **bondtypename;
00172 
00173     // must explicitly set these since they may are otherwise only 
00174     // initialized by the read_bonds() call in the new ABI
00175     nbondtypes = 0;
00176     bondtype = NULL;
00177     bondtypename = NULL;
00178  
00179 #if vmdplugin_ABIVERSION >= 15
00180     if (plugin->read_bonds(rv, &nbonds, &from, &to, &bondorder, &bondtype, 
00181                            &nbondtypes, &bondtypename)) 
00182 #else
00183     if (plugin->read_bonds(rv, &nbonds, &from, &to, &bondorder)) 
00184 #endif
00185     {
00186       
00187       msgErr << "Error reading bond information." << sendmsg;
00188       if (autobonds)
00189         m->find_bonds_from_timestep();
00190     } else {
00191       // if we didn't get an error, but the plugin didn't define bonds
00192       // for some reason (i.e. no usable PDB CONECT records found)
00193       // fall back to automatic bond search
00194       if (nbonds == 0 && from == NULL && to == NULL) {
00195         if (autobonds)
00196           m->find_bonds_from_timestep();
00197       } else if (nbonds > 0) {
00198         // insert bondtype names into our pool of names. this will preserve 
00199         // the indexing if no other names do already exist.
00200         if (bondtypename != NULL)
00201           for (i=0; i < nbondtypes; i++)
00202             m->bondTypeNames.add_name(bondtypename[i], i);
00203 
00204         // If the bonds provided by the plugin were only the for
00205         // non-standard residues (i.e. PDB CONECT records) then we'll
00206         // also do our regular bond search, combining all together, but 
00207         // preventing duplicates
00208         if (optflags & MOLFILE_BONDSSPECIAL) {
00209           if (autobonds) {
00210             m->find_unique_bonds_from_timestep(); // checking for duplicates
00211           }
00212 
00213           // indicate what kind of data will be available.
00214           m->set_dataset_flag(BaseMolecule::BONDS);
00215           if (bondorder != NULL) 
00216             m->set_dataset_flag(BaseMolecule::BONDORDERS);
00217           if (bondtype  != NULL) 
00218             m->set_dataset_flag(BaseMolecule::BONDTYPES);
00219           
00220           // Now add in all of the special bonds provided by the plugin
00221           for (i=0; i<nbonds; i++) {
00222             float bo =  1;
00223             int   bt = -1;
00224             if (bondorder != NULL) 
00225               bo = bondorder[i];
00226             if (bondtype  != NULL) {
00227               if (bondtypename != NULL) {
00228                 bt = m->bondTypeNames.typecode(bondtypename[bondtype[i]]);
00229               } else {
00230                 bt = bondtype[i];
00231               }
00232             }
00233             m->add_bond_dupcheck(from[i]-1, to[i]-1, bo, bt); 
00234           }
00235         } else {
00236 
00237           // indicate what kind of data will be available...
00238           m->set_dataset_flag(BaseMolecule::BONDS);
00239           if (bondorder != NULL) 
00240             m->set_dataset_flag(BaseMolecule::BONDORDERS);
00241           if (bondtype  != NULL) 
00242             m->set_dataset_flag(BaseMolecule::BONDTYPES);
00243 
00244           // ...and add the bonds.
00245           for (i=0; i<nbonds; i++) {
00246             float bo =  1;
00247             int   bt = -1;
00248             if (bondorder != NULL) bo = bondorder[i];
00249             if (bondtype  != NULL) bt = bondtype[i];
00250             m->add_bond(from[i]-1, to[i]-1, bo, bt); 
00251           }
00252         }
00253       }
00254     }
00255   } else {
00256     if (autobonds)
00257       m->find_bonds_from_timestep();
00258   }
00259 
00260   // if the plugin can read angles, dihedrals, impropers, cross-terms, 
00261   // do it here...
00262   if (can_read_angles()) {
00263     int numangles, *angles, *angletypes, numangletypes;
00264     int numdihedrals, *dihedrals, *dihedraltypes, numdihedraltypes;
00265     int numimpropers, *impropers, *impropertypes, numimpropertypes; 
00266     int numcterms, *cterms, ctermcols, ctermrows;
00267     char **angletypenames, **dihedraltypenames, **impropertypenames;
00268 
00269     // initialize to empty so that the rest of the code works well.
00270     angletypes = dihedraltypes = impropertypes = NULL;
00271     numangletypes = numdihedraltypes = numimpropertypes = 0;
00272     angletypenames = dihedraltypenames = impropertypenames = NULL;
00273     
00274 #if vmdplugin_ABIVERSION >= 16
00275     if (plugin->read_angles(rv, &numangles, &angles, &angletypes, 
00276                             &numangletypes, &angletypenames, &numdihedrals,
00277                             &dihedrals,  &dihedraltypes, &numdihedraltypes, 
00278                             &dihedraltypenames, &numimpropers, &impropers,
00279                             &impropertypes, &numimpropertypes, 
00280                             &impropertypenames, &numcterms, &cterms, 
00281                             &ctermcols, &ctermrows))
00282       
00283 #else
00284     double *angleforces, *dihedralforces, *improperforces, *ctermforces;
00285     angleforces = dihedralforces = improperforces = ctermforces = NULL;
00286     if (plugin->read_angles(rv, 
00287                             &numangles, &angles, &angleforces,
00288                             &numdihedrals, &dihedrals, &dihedralforces,
00289                             &numimpropers, &impropers, &improperforces,
00290                             &numcterms, &cterms, &ctermcols, &ctermrows, 
00291                             &ctermforces)) 
00292 #endif
00293     {
00294       msgErr << "Error reading angle and cross-term information." << sendmsg;
00295     } else {
00296       // we consider angles, dihedrals and impropers all as "angles".
00297       // it is hard to come up with a scenario where something else
00298       // as an all-or-nothing setup would happen.
00299       if ( (angletypes != NULL) || (dihedraltypes != NULL) 
00300            || (impropertypes != NULL) )
00301         m->set_dataset_flag(BaseMolecule::ANGLETYPES);
00302 
00303       // insert type names into our pool of names. this will preserve 
00304       // the indexing if no other names do already exist.
00305       if (angletypenames != NULL)
00306         for (i=0; i < numangletypes; i++)
00307           m->angleTypeNames.add_name(angletypenames[i], i);
00308 
00309       if (dihedraltypenames != NULL)
00310         for (i=0; i < numdihedraltypes; i++)
00311           m->dihedralTypeNames.add_name(dihedraltypenames[i], i);
00312 
00313       if (impropertypenames != NULL)
00314         for (i=0; i < numimpropertypes; i++)
00315           m->improperTypeNames.add_name(impropertypenames[i], i);
00316 
00317       // XXX: using a similar interface as with bonds would become
00318       //      quite complicated. not sure if it is really needed.
00319       //      for the moment we forgo it and force a simple all-or-nothing
00320       //      scheme instead. incrementally adding and removing of
00321       //      angles has to be done from (topology building) scripts.
00322       if (numangles > 0 || numdihedrals > 0 || numimpropers > 0) {
00323         m->set_dataset_flag(BaseMolecule::ANGLES);
00324         m->clear_angles();
00325         if (angletypes != NULL) {
00326           int type;
00327 
00328           for (i=0; i<numangles; i++) {
00329             type = angletypes[i];
00330             if (angletypenames != NULL)
00331               type = m->angleTypeNames.typecode(angletypenames[type]);
00332 
00333             m->add_angle(angles[3L*i]-1, angles[3L*i+1]-1, angles[3L*i+2]-1, type);
00334           }
00335         } else {
00336           for (i=0; i<numangles; i++) {
00337             m->add_angle(angles[3L*i]-1, angles[3L*i+1]-1, angles[3L*i+2]-1);
00338           }
00339         }
00340 
00341         m->clear_dihedrals();
00342         if (dihedraltypes != NULL) {
00343           int type;
00344 
00345           for (i=0; i<numdihedrals; i++) {
00346             type = dihedraltypes[i];
00347             if (dihedraltypenames != NULL)
00348               type = m->dihedralTypeNames.typecode(dihedraltypenames[type]);
00349 
00350             m->add_dihedral(dihedrals[4L*i]-1,   dihedrals[4L*i+1]-1, 
00351                             dihedrals[4L*i+2]-1, dihedrals[4L*i+3]-1, type);
00352           }
00353         } else {
00354           for (i=0; i<numdihedrals; i++) {
00355             m->add_dihedral(dihedrals[4L*i]-1,   dihedrals[4L*i+1]-1, 
00356                             dihedrals[4L*i+2]-1, dihedrals[4L*i+3]-1);
00357           }
00358         }
00359         
00360         m->clear_impropers();
00361         if (impropertypes != NULL) {
00362           int type;
00363 
00364           for (i=0; i<numimpropers; i++) {
00365             type = impropertypes[i];
00366             if (impropertypenames != NULL)
00367               type = m->improperTypeNames.typecode(impropertypenames[type]);
00368 
00369             m->add_improper(impropers[4L*i]-1,   impropers[4L*i+1]-1, 
00370                             impropers[4L*i+2]-1, impropers[4L*i+3]-1, type);
00371           }
00372         } else {
00373           for (i=0; i<numimpropers; i++) {
00374             m->add_improper(impropers[4L*i]-1,   impropers[4L*i+1]-1, 
00375                             impropers[4L*i+2]-1, impropers[4L*i+3]-1);
00376           }
00377         }
00378       }
00379 
00380       // cross terms have no types. they are a CHARMM-only thing.
00381       if (numcterms > 0) {
00382         m->set_dataset_flag(BaseMolecule::CTERMS);
00383         for (i=0; i<numcterms; i++)
00384           m->add_cterm(cterms[8L*i]-1,   cterms[8L*i+1]-1, cterms[8L*i+2]-1,
00385                        cterms[8L*i+3]-1, cterms[8L*i+4]-1, cterms[8L*i+5]-1,
00386                        cterms[8L*i+6]-1, cterms[8L*i+7]-1);
00387       }
00388     }
00389   }
00390 
00391   return MOLFILE_SUCCESS;
00392 }
00393 
00394   
00395 int MolFilePlugin::read_optional_structure(Molecule *m, int filebonds) {
00396   if (!rv) return MOLFILE_ERROR; 
00397   if (!can_read_structure()) return MOLFILE_ERROR;
00398   if (numatoms != m->nAtoms) return MOLFILE_ERROR;
00399 
00400   molfile_atom_t *atoms = (molfile_atom_t *) calloc(1, numatoms*sizeof(molfile_atom_t));
00401 
00402   int rc, i;
00403   int optflags = MOLFILE_BADOPTIONS; /* plugin must reset this correctly */
00404   if ((rc = plugin->read_structure(rv, &optflags, atoms))) {
00405     free(atoms);
00406     return rc; // propagate error to caller
00407   }
00408   if (optflags == int(MOLFILE_BADOPTIONS)) {
00409     free(atoms);
00410     msgErr << "MolFilePlugin: plugin didn't initialize optional data flags" << sendmsg;
00411     msgErr << "MolFilePlugin: file load aborted" << sendmsg;
00412     return MOLFILE_ERROR; /* abort load, data can't be trusted */
00413   }
00414 
00415   // set molecule dataset flags
00416   if (optflags & MOLFILE_OCCUPANCY)
00417     m->set_dataset_flag(BaseMolecule::OCCUPANCY);
00418   if (optflags & MOLFILE_BFACTOR)
00419     m->set_dataset_flag(BaseMolecule::BFACTOR);
00420   if (optflags & MOLFILE_MASS)
00421     m->set_dataset_flag(BaseMolecule::MASS);
00422   if (optflags & MOLFILE_CHARGE)
00423     m->set_dataset_flag(BaseMolecule::CHARGE);
00424   if (optflags & MOLFILE_RADIUS)
00425     m->set_dataset_flag(BaseMolecule::RADIUS);
00426   if (optflags & MOLFILE_ATOMICNUMBER)
00427     m->set_dataset_flag(BaseMolecule::ATOMICNUMBER);
00428 
00429   // Force min/max atom radii to be updated since we've loaded new data
00430   if (optflags & MOLFILE_RADIUS)
00431     m->set_radii_changed();
00432 
00433   float *charge = m->charge();
00434   float *mass = m->mass();
00435   float *radius = m->radius();
00436   float *beta = m->beta();
00437   float *occupancy = m->occupancy();
00438   for (i=0; i<numatoms; i++) {
00439     if (optflags & MOLFILE_OCCUPANCY) {
00440       occupancy[i] = atoms[i].occupancy;
00441     }
00442 
00443     if (optflags & MOLFILE_BFACTOR) {
00444       beta[i] = atoms[i].bfactor;
00445     }
00446 
00447     if (optflags & MOLFILE_MASS) {
00448       mass[i] = atoms[i].mass;
00449     }
00450 
00451     if (optflags & MOLFILE_CHARGE) {
00452       charge[i] = atoms[i].charge;
00453     }
00454 
00455     if (optflags & MOLFILE_RADIUS) {
00456       radius[i] = atoms[i].radius;
00457     }
00458 
00459     if (optflags & MOLFILE_ATOMICNUMBER) {
00460       m->atom(i)->atomicnumber = atoms[i].atomicnumber;
00461     }
00462   }
00463   free(atoms);
00464 
00465   // if no bonds are added, then we do not re-analyze the structure
00466   if (!can_read_bonds()) 
00467     return MOLFILE_SUCCESS;
00468 
00469   // When tacking on trajectory frames from PDB files with CONECT records,
00470   // we have to prevent partial CONECT record bonding information from
00471   // blowing away existing connectivity information derived from 
00472   // automatic bond search results or from a previously loaded file with
00473   // complete bond information such as a PSF.  We only accept 
00474   // complete bonding connectivity updates, no partial updates are allowed.
00475   // Also bonding information updates can be disabled by the user with the
00476   // filebonds flag.
00477   if (!(optflags & MOLFILE_BONDSSPECIAL) && filebonds) {
00478     int nbonds, *from, *to, nbondtypes, *bondtype;
00479     float *bondorder;
00480     char **bondtypename;
00481   
00482     // must explicitly set these since they may are otherwise only
00483     // initialized by the read_bonds() call in the new ABI
00484     nbondtypes = 0;
00485     bondtype = NULL;
00486     bondtypename = NULL;
00487 
00488 #if vmdplugin_ABIVERSION >= 15
00489     if (plugin->read_bonds(rv, &nbonds, &from, &to, &bondorder, &bondtype, 
00490                            &nbondtypes, &bondtypename)) 
00491       return MOLFILE_SUCCESS;
00492 #else
00493     if (plugin->read_bonds(rv, &nbonds, &from, &to, &bondorder)) 
00494       return MOLFILE_SUCCESS;
00495 #endif
00496 
00497     // insert bondtype names into our pool of names. this will preserve 
00498     // the indexing if no other names do already exist.
00499     if (bondtypename != NULL)
00500       for (i=0; i < nbondtypes; i++)
00501         m->bondTypeNames.add_name(bondtypename[i], i);
00502       
00503     if (nbonds == 0) 
00504       return MOLFILE_SUCCESS;
00505 
00506     for (i=0; i<numatoms; i++) 
00507       m->atom(i)->bonds = 0;
00508 
00509     m->set_dataset_flag(BaseMolecule::BONDS);
00510     if (bondorder != NULL)
00511       m->set_dataset_flag(BaseMolecule::BONDORDERS);
00512     if (bondtype != NULL)
00513       m->set_dataset_flag(BaseMolecule::BONDTYPES);
00514 
00515     for (i=0; i<nbonds; i++) {
00516       float bo = 1;
00517       int bt = -1;
00518 
00519       if (bondorder != NULL) 
00520         bo = bondorder[i];
00521       if (bondtype != NULL) {
00522         if (bondtypename != NULL) {
00523           bt = m->bondTypeNames.typecode(bondtypename[bondtype[i]]);
00524         } else {
00525           bt = bondtype[i];
00526         }
00527       }
00528         
00529       m->add_bond(from[i]-1, to[i]-1, bo, bt); // real bond order
00530     } 
00531   } else {
00532     // if no bonds are added, then we do not re-analyze the structure
00533     return MOLFILE_SUCCESS;
00534   }
00535 
00536   // if the plugin can read angles, dihedrals, impropers, cross-terms, 
00537   // do it here...
00538   if (can_read_angles()) {
00539     int numangles, *angles, *angletypes, numangletypes;
00540     int numdihedrals, *dihedrals, *dihedraltypes, numdihedraltypes;
00541     int numimpropers, *impropers, *impropertypes, numimpropertypes; 
00542     int numcterms, *cterms, ctermcols, ctermrows;
00543     char **angletypenames, **dihedraltypenames, **impropertypenames;
00544 
00545     // initialize to empty so that the rest of the code works well.
00546     angletypes = dihedraltypes = impropertypes = NULL;
00547     numangletypes = numdihedraltypes = numimpropertypes = 0;
00548     angletypenames = dihedraltypenames = impropertypenames = NULL;
00549     
00550 #if vmdplugin_ABIVERSION >= 16
00551     if (plugin->read_angles(rv, &numangles, &angles, &angletypes, 
00552                             &numangletypes, &angletypenames, &numdihedrals,
00553                             &dihedrals,  &dihedraltypes, &numdihedraltypes, 
00554                             &dihedraltypenames, &numimpropers, &impropers,
00555                             &impropertypes, &numimpropertypes, 
00556                             &impropertypenames, &numcterms, &cterms, 
00557                             &ctermcols, &ctermrows))
00558       
00559 #else
00560     double *angleforces, *dihedralforces, *improperforces, *ctermforces;
00561     angleforces = dihedralforces = improperforces = ctermforces = NULL;
00562     if (plugin->read_angles(rv, 
00563                             &numangles, &angles, &angleforces,
00564                             &numdihedrals, &dihedrals, &dihedralforces,
00565                             &numimpropers, &impropers, &improperforces,
00566                             &numcterms, &cterms, &ctermcols, &ctermrows, 
00567                             &ctermforces)) 
00568 #endif
00569     {
00570       msgErr << "Error reading angle and cross-term information." << sendmsg;
00571     } else {
00572       // we consider angles, dihedrals and impropers all as "angles".
00573       // it is hard to come up with a scenario where something else
00574       // as an all-or-nothing setup would happen.
00575       if ( (angletypes != NULL) || (dihedraltypes != NULL) 
00576            || (impropertypes != NULL) )
00577         m->set_dataset_flag(BaseMolecule::ANGLETYPES);
00578 
00579       // insert type names into our pool of names. this will preserve 
00580       // the indexing if no other names do already exist.
00581       if (angletypenames != NULL)
00582         for (i=0; i < numangletypes; i++)
00583           m->angleTypeNames.add_name(angletypenames[i], i);
00584 
00585       if (dihedraltypenames != NULL)
00586         for (i=0; i < numdihedraltypes; i++)
00587           m->dihedralTypeNames.add_name(dihedraltypenames[i], i);
00588 
00589       if (impropertypenames != NULL)
00590         for (i=0; i < numimpropertypes; i++)
00591           m->improperTypeNames.add_name(impropertypenames[i], i);
00592 
00593       // XXX: using a similar interface as with bonds would become
00594       //      quite complicated. not sure if it is really needed.
00595       //      for the moment we force a simple all-or-nothing
00596       //      scheme and only allow incrementally adding and removing
00597       //      angles from (topology building) scripts.
00598       if (numangles > 0 || numdihedrals > 0 || numimpropers > 0) {
00599         m->set_dataset_flag(BaseMolecule::ANGLES);
00600         m->clear_angles();
00601         if (angletypes != NULL) {
00602           int type;
00603 
00604           for (i=0; i<numangles; i++) {
00605             type = angletypes[i];
00606             if (angletypenames != NULL)
00607               type = m->angleTypeNames.typecode(angletypenames[type]);
00608 
00609             m->add_angle(angles[3L*i]-1, angles[3L*i+1]-1, angles[3L*i+2]-1, type);
00610           }
00611         } else {
00612           for (i=0; i<numangles; i++) {
00613             m->add_angle(angles[3L*i]-1, angles[3L*i+1]-1, angles[3L*i+2]-1);
00614           }
00615         }
00616 
00617         m->clear_dihedrals();
00618         if (dihedraltypes != NULL) {
00619           int type;
00620 
00621           for (i=0; i<numdihedrals; i++) {
00622             type = dihedraltypes[i];
00623             if (dihedraltypenames != NULL)
00624               type = m->dihedralTypeNames.typecode(dihedraltypenames[type]);
00625 
00626             m->add_dihedral(dihedrals[4L*i]-1,   dihedrals[4L*i+1]-1, 
00627                             dihedrals[4L*i+2]-1, dihedrals[4L*i+3]-1, type);
00628           }
00629         } else {
00630           for (i=0; i<numdihedrals; i++) {
00631             m->add_dihedral(dihedrals[4L*i]-1,   dihedrals[4L*i+1]-1, 
00632                             dihedrals[4L*i+2]-1, dihedrals[4L*i+3]-1);
00633           }
00634         }
00635         
00636         m->clear_impropers();
00637         if (impropertypes != NULL) {
00638           int type;
00639 
00640           for (i=0; i<numimpropers; i++) {
00641             type = impropertypes[i];
00642             if (impropertypenames != NULL)
00643               type = m->improperTypeNames.typecode(impropertypenames[type]);
00644 
00645             m->add_improper(impropers[4L*i]-1,   impropers[4L*i+1]-1, 
00646                             impropers[4L*i+2]-1, impropers[4L*i+3]-1, type);
00647           }
00648         } else {
00649           for (i=0; i<numimpropers; i++) {
00650             m->add_improper(impropers[4L*i]-1,   impropers[4L*i+1]-1, 
00651                             impropers[4L*i+2]-1, impropers[4L*i+3]-1);
00652           }
00653         }
00654       }
00655 
00656       // cross terms have no types. they are a CHARMM-only thing.
00657       if (numcterms > 0) {
00658         m->set_dataset_flag(BaseMolecule::CTERMS);
00659         for (i=0; i<numcterms; i++)
00660           m->add_cterm(cterms[8L*i]-1,   cterms[8L*i+1]-1, cterms[8L*i+2]-1,
00661                        cterms[8L*i+3]-1, cterms[8L*i+4]-1, cterms[8L*i+5]-1,
00662                        cterms[8L*i+6]-1, cterms[8L*i+7]-1);
00663       }
00664     }
00665   }
00666 
00667   // (re)analyze the molecular structure, since bonds/angles/etc
00668   // may have been changed
00669   m->analyze();
00670 
00671   // force all reps and selections to be recalculated
00672   m->force_recalc(DrawMolItem::COL_REGEN | DrawMolItem::SEL_REGEN);
00673 
00674   // force secondary structure to be recalculated
00675   m->invalidate_ss(); 
00676  
00677   return MOLFILE_SUCCESS;
00678 }
00679 
00680 void MolFilePlugin::set_natoms(int n) {
00681   numatoms = n;
00682 }
00683 
00684 int MolFilePlugin::init_read(const char *file) {
00685   rv = NULL;
00686   if (can_read_structure() || can_read_timesteps() || can_read_graphics() ||
00687       can_read_volumetric() || can_read_metadata() || can_read_qm()) { 
00688     rv = plugin->open_file_read(file, plugin->name, &numatoms);
00689   }
00690   if (!rv) return MOLFILE_ERROR;
00691   delete [] _filename;
00692   _filename = stringdup(file);
00693   return MOLFILE_SUCCESS;
00694 }
00695 
00696 
00697 int MolFilePlugin::read_timestep_pagealign_size(void) {
00698 #if vmdplugin_ABIVERSION > 17
00699   if (!rv) return 1;
00700   if (can_read_pagealigned_timesteps()) {
00701     int sz;
00702     plugin->read_timestep_pagealign_size(rv, &sz);
00703 
00704     // If the page alignment size passes sanity check, we use it...
00705     if ((sz != 1) && 
00706         ((sz < MOLFILE_DIRECTIO_MIN_BLOCK_SIZE) &&
00707          (sz > MOLFILE_DIRECTIO_MAX_BLOCK_SIZE))) {
00708       msgWarn << "Plugin returned bad page alignment size!: " 
00709               << sz << sendmsg;
00710     }
00711 
00712     if (getenv("VMDPLUGINVERBOSE") != NULL) {
00713       msgInfo << "Plugin returned page alignment size: " << sz << sendmsg;
00714     }
00715 
00716     return sz;
00717   } else {
00718     if (getenv("VMDPLUGINVERBOSE") != NULL) {
00719       msgInfo << "Plugin can't read page aligned timesteps." << sendmsg;
00720     }
00721   }
00722 
00723   return 1;
00724 #else
00725 #if 1
00726   return 1; // assume non-blocked I/O
00727 #else
00728   // Enable VMD to cope with hard-coded revs of jsplugin if we want
00729   return MOLFILE_DIRECTIO_MAX_BLOCK_SIZE;
00730 #endif
00731 #endif
00732 }
00733 
00734 
00735 Timestep *MolFilePlugin::next(Molecule *m, int ts_pagealign_sz) {
00736   if (!rv) return NULL;
00737   if (numatoms <= 0) return NULL;
00738   if (!(can_read_timesteps() || can_read_qm_timestep())) return NULL;
00739   molfile_timestep_t timestep;
00740   memset(&timestep, 0, sizeof(molfile_timestep_t));
00741 
00742   // allocate space for velocities only if 
00743   // 1) the plugin implements read_timestep_metadata;
00744   // 2) metadata->has_velocities is TRUE.
00745   float *velocities = NULL;
00746 
00747   // XXX this needs to be pulled out into the read init code
00748   // rather than being done once per timestep
00749   molfile_timestep_metadata_t meta;
00750   if (can_read_timestep_metadata()) {
00751     memset(&meta, 0, sizeof(molfile_timestep_metadata));
00752     plugin->read_timestep_metadata(rv, &meta);
00753     if (meta.has_velocities) {
00754       velocities = new float[3L*numatoms];
00755     }
00756   }
00757 
00758   // QM timestep metadata can be different on every single timestep
00759   // so we need to query it prior to reading the timestep so we can 
00760   // allocate the right buffer sizes etc. 
00761   molfile_qm_timestep_metadata_t qmmeta;
00762   if (can_read_qm_timestep_metadata()) {
00763     memset(&qmmeta, 0, sizeof(molfile_qm_timestep_metadata));
00764     // XXX need to add timestep parameter or other method to specify
00765     //     which frame this applies to, else keep it the way it is
00766     //     and rename the plugin function appropriately.
00767     plugin->read_qm_timestep_metadata(rv, &qmmeta);
00768   }
00769 
00770   // set useful defaults for unit cell information
00771   // a non-periodic structure has cell lengths of zero
00772   // if it's periodic in less than 3 dimensions, then only the
00773   // periodic directions will be non-zero.
00774   timestep.A = timestep.B = timestep.C = 0.0f;
00775 
00776   // cells are rectangular until told otherwise
00777   timestep.alpha = timestep.beta = timestep.gamma = 90.0f;
00778 
00779   Timestep *ts = new Timestep(numatoms, ts_pagealign_sz);
00780   timestep.coords = ts->pos; 
00781   timestep.velocities = velocities;
00782   ts->vel = velocities;
00783  
00784   int rc = 0;
00785   molfile_qm_metadata_t *qm_metadata = NULL; // this just a dummy
00786   molfile_qm_timestep_t qm_timestep;
00787   memset(&qm_timestep, 0, sizeof(molfile_qm_timestep_t));
00788 
00789   // XXX this needs to be fixed so that a file format that can 
00790   //     optionally contain either QM or non-QM data will work correctly
00791   if (can_read_qm_timestep()) {
00792     qm_timestep.scfenergies = new double[qmmeta.num_scfiter];
00793     qm_timestep.wave = new molfile_qm_wavefunction_t[qmmeta.num_wavef];
00794     memset(qm_timestep.wave, 0, qmmeta.num_wavef*sizeof(molfile_qm_wavefunction_t));
00795     int i;
00796     for (i=0; (i<MOLFILE_MAXWAVEPERTS && i<qmmeta.num_wavef); i++) {
00797       qm_timestep.wave[i].wave_coeffs = 
00798         new float[qmmeta.num_orbitals_per_wavef[i]*qmmeta.wavef_size];
00799       if (qmmeta.has_orben_per_wavef[i]) {
00800         qm_timestep.wave[i].orbital_energies =
00801           new float[qmmeta.num_orbitals_per_wavef[i]];
00802       }
00803       if (qmmeta.has_occup_per_wavef[i]) {
00804         qm_timestep.wave[i].occupancies = 
00805           new float[qmmeta.num_orbitals_per_wavef[i]];
00806       }
00807     }
00808     if (qmmeta.has_gradient) {
00809       qm_timestep.gradient = new float[3L*numatoms];
00810     }
00811     if (qmmeta.num_charge_sets) {
00812       qm_timestep.charges = new double[numatoms*qmmeta.num_charge_sets];
00813       qm_timestep.charge_types = new int[qmmeta.num_charge_sets];
00814     }
00815     rc = plugin->read_timestep(rv, numatoms, &timestep, qm_metadata, &qm_timestep);
00816   } else {
00817     rc = plugin->read_next_timestep(rv, numatoms, &timestep);
00818   }
00819 
00820   if (rc) {
00821     if (can_read_qm_timestep()) {
00822       delete [] qm_timestep.scfenergies;
00823       if (qm_timestep.gradient) delete [] qm_timestep.gradient;
00824       if (qm_timestep.charges)  delete [] qm_timestep.charges;
00825       if (qm_timestep.charge_types) delete [] qm_timestep.charge_types;
00826       int i;
00827       for (i=0; i<qmmeta.num_wavef; i++) {
00828         delete [] qm_timestep.wave[i].wave_coeffs;
00829         delete [] qm_timestep.wave[i].orbital_energies;
00830       }
00831       delete [] qm_timestep.wave;
00832     }
00833     delete ts;
00834     ts = NULL; 
00835   } else {
00836     ts->a_length = timestep.A;
00837     ts->b_length = timestep.B;
00838     ts->c_length = timestep.C;
00839     ts->alpha = timestep.alpha;
00840     ts->beta = timestep.beta;
00841     ts->gamma = timestep.gamma;
00842     ts->physical_time = timestep.physical_time;
00843     if (can_read_qm_timestep()) {
00844       int i;
00845       int *chargetypes = new int[qmmeta.num_charge_sets];
00846       for (i=0; i<qmmeta.num_charge_sets; i++) {
00847         switch (qm_timestep.charge_types[i]) {
00848         case  MOLFILE_QMCHARGE_MULLIKEN:
00849           chargetypes[i] = QMCHARGE_MULLIKEN;
00850           break;
00851         case MOLFILE_QMCHARGE_LOWDIN:
00852           chargetypes[i] = QMCHARGE_LOWDIN;
00853           break;
00854         case MOLFILE_QMCHARGE_ESP:
00855           chargetypes[i] = QMCHARGE_ESP;
00856           break;
00857         case MOLFILE_QMCHARGE_NPA:
00858           chargetypes[i] = QMCHARGE_NPA;
00859           break;
00860         default:
00861           chargetypes[i] = QMCHARGE_UNKNOWN;
00862         }
00863       }
00864 
00865       ts->qm_timestep = new QMTimestep(numatoms);
00866       ts->qm_timestep->set_charges(qm_timestep.charges,
00867                                    chargetypes,
00868                                    numatoms, qmmeta.num_charge_sets);
00869       delete [] chargetypes;
00870 
00871       ts->qm_timestep->set_scfenergies(qm_timestep.scfenergies,
00872                                        qmmeta.num_scfiter);
00873 
00874 
00875       // signa_ts is the list of numsigts signatures for the wavefunctions
00876       // already processed for this timestep. 
00877       int num_signa_ts = 0;
00878       wavef_signa_t *signa_ts = NULL;
00879       for (i=0; i<qmmeta.num_wavef; i++) {
00880         // We need to translate between the macros used in the plugins
00881         // and the one ones in VMD, here so that they can be independent
00882         // and we don't have to include molfile_plugin.h anywhere else
00883         // in VMD.
00884         int wavef_type;
00885         switch (qm_timestep.wave[i].type) {
00886         case MOLFILE_WAVE_CANON:
00887           wavef_type = WAVE_CANON;
00888           break;
00889         case MOLFILE_WAVE_CINATUR:
00890           wavef_type = WAVE_CINATUR;
00891           break;
00892         case MOLFILE_WAVE_GEMINAL:
00893           wavef_type = WAVE_GEMINAL;
00894           break;
00895         case MOLFILE_WAVE_BOYS:
00896           wavef_type = WAVE_BOYS;
00897           break;
00898         case MOLFILE_WAVE_RUEDEN:
00899           wavef_type = WAVE_RUEDEN;
00900           break;
00901         case MOLFILE_WAVE_PIPEK:
00902           wavef_type = WAVE_PIPEK;
00903           break;
00904         case MOLFILE_WAVE_MCSCFOPT:
00905           wavef_type = WAVE_MCSCFOPT;
00906           break;
00907         case MOLFILE_WAVE_MCSCFNAT:
00908           wavef_type = WAVE_MCSCFNAT;
00909           break;
00910         default:
00911           wavef_type = WAVE_UNKNOWN;
00912         }
00913 
00914         // Add new wavefunction to QM timestep
00915         ts->qm_timestep->add_wavefunction(qm_data,
00916                            qmmeta.wavef_size,
00917                            qmmeta.num_orbitals_per_wavef[i],
00918                            qm_timestep.wave[i].wave_coeffs,
00919                            qm_timestep.wave[i].orbital_energies,
00920                            qm_timestep.wave[i].occupancies,
00921                            qm_timestep.wave[i].orbital_ids,
00922                            qm_timestep.wave[i].energy,
00923                            wavef_type,
00924                            qm_timestep.wave[i].spin,
00925                            qm_timestep.wave[i].excitation,
00926                            qm_timestep.wave[i].multiplicity,
00927                            qm_timestep.wave[i].info,
00928                            signa_ts, num_signa_ts);
00929       }
00930       free(signa_ts);
00931 
00932       ts->qm_timestep->set_gradients(qm_timestep.gradient, numatoms);
00933 
00934       delete [] qm_timestep.scfenergies;
00935       if (qm_timestep.gradient)     delete [] qm_timestep.gradient;
00936       if (qm_timestep.charges)      delete [] qm_timestep.charges;
00937       if (qm_timestep.charge_types) delete [] qm_timestep.charge_types;
00938 
00939       for (i=0; i<qmmeta.num_wavef; i++) {
00940         delete [] qm_timestep.wave[i].wave_coeffs;
00941         delete [] qm_timestep.wave[i].orbital_energies;
00942         delete [] qm_timestep.wave[i].occupancies;
00943       }
00944       delete [] qm_timestep.wave;
00945 #if 0
00946       // If we have at least two timesteps then sort the orbitals
00947       // of the current timestep according to the ones from the
00948       // previous timestep.
00949       // Note that the frame counter m->numframes() is updated
00950       // after this function call. 
00951       if (m->numframes()>0) {
00952         // Get the previous Timestep
00953         Timestep *prevts = m->get_frame(m->numframes()-1);
00954 
00955         msgInfo << "sort frame " << m->numframes() << sendmsg;
00956 
00957         // Sort the orbitals by comparing them to the ones from
00958         // the previous timestep.
00959         ts->qm_timestep->sort_orbitals(prevts->qm_timestep);
00960       } else {
00961         msgInfo << "ignore frame " << m->numframes() << sendmsg;
00962       }
00963 #endif 
00964     }
00965   }
00966   return ts;
00967 }
00968 
00969 int MolFilePlugin::skip(Molecule *m) {
00970   if (!rv) return MOLFILE_ERROR;
00971   if (!can_read_timesteps()) return MOLFILE_ERROR;
00972   return plugin->read_next_timestep(rv, numatoms, 0);
00973 }
00974 
00975 void MolFilePlugin::close() {
00976   if (rv && (can_read_structure() || can_read_timesteps() || 
00977              can_read_graphics() || can_read_volumetric() ||
00978              can_read_metadata() || can_read_bonds())) { 
00979     plugin->close_file_read(rv);
00980     rv = NULL;
00981   }
00982   if (wv && (can_write_structure() || can_write_timesteps() ||
00983              can_write_bonds())) { 
00984     plugin->close_file_write(wv);
00985     wv = NULL;
00986   }
00987 }
00988 
00989 int MolFilePlugin::init_write(const char *file, int natoms) {
00990   wv = NULL;
00991   if (can_write_structure() || can_write_timesteps() || 
00992       can_write_volumetric()) { 
00993     wv = plugin->open_file_write(file, plugin->name, natoms);
00994   } 
00995   if (!wv) return MOLFILE_ERROR;
00996   delete [] _filename;
00997   _filename = stringdup(file);
00998 
00999   // Cache the number of atoms to be written.  It's not strictly necessary,
01000   // but it lets us allocate only the necessary space in write_structure
01001   // and write_timestep.
01002   numatoms = natoms;
01003 
01004   return MOLFILE_SUCCESS;
01005 }
01006 
01007 
01008 int MolFilePlugin::write_structure(Molecule *m, const int *on) {
01009   if (!wv) return MOLFILE_ERROR;
01010   if (!can_write_structure()) return MOLFILE_ERROR;
01011   if (!m->has_structure()) {
01012     msgErr << "Molecule's structure has not been initialized." << sendmsg;
01013     return MOLFILE_ERROR;
01014   }
01015   
01016   ptrdiff_t i, j, k;
01017   molfile_atom_t *atoms = (molfile_atom_t *) calloc(1, numatoms*sizeof(molfile_atom_t));
01018   int *atomindexmap = (int *) calloc(1, m->nAtoms*sizeof(int));
01019 
01020   // initialize the atom index map to an invalid atom index value, so that
01021   // we can use this to eliminate bonds to atoms that aren't selected.
01022   for (i=0; i<m->nAtoms; i++) 
01023     atomindexmap[i] = -1;
01024 
01025   const float *charge = m->charge();
01026   const float *mass = m->mass();
01027   const float *radius = m->radius();
01028   const float *beta = m->beta();
01029   const float *occupancy = m->occupancy();
01030 
01031   int mangleatomnames = 1;
01032   if (getenv("VMDNOMANGLEATOMNAMES")) {
01033     mangleatomnames = 0;
01034   }
01035 
01036   // build the array of selected atoms to be written out
01037   for (i=0, j=0; i<m->nAtoms; i++) {
01038     // skip atoms that aren't 'on', if we've been given an 'on' array
01039     if (on && !on[i]) 
01040       continue;
01041 
01042     // Check that the number of atoms specified in init_write is no smaller
01043     // than the number of atoms in the selection; otherwise we would 
01044     // corrupt memory.
01045     if (j >= numatoms) {
01046       msgErr << 
01047         "MolFilePlugin: Internal error, selection size exceeds numatoms ("
01048         << numatoms << ")" << sendmsg;
01049       free(atoms);
01050       free(atomindexmap);
01051       return MOLFILE_ERROR;
01052     }
01053 
01054     const MolAtom *atom = m->atom(i);
01055     molfile_atom_t &atm = atoms[j];
01056 
01057     if (mangleatomnames) {
01058       // Try to restore the spacing on the name since VMD destroys it when it
01059       // reads it in.
01060       // XXX this is PDB-centric thinking, need to reconsider this
01061       char name[6], *nameptr;
01062       name[0] = ' ';
01063       strncpy(name+1, (m->atomNames).name(atom->nameindex), 4);
01064       name[5] = '\0';
01065       // the name must be left-justified
01066       if(strlen(name) == 5) {
01067         nameptr = name + 1;
01068       } else {
01069         nameptr = name;
01070         int p;
01071         while((p = strlen(name)) < 4) {
01072           name[p] = ' ';
01073           name[p+1] = '\0';
01074         }
01075       }
01076       strcpy(atm.name, nameptr);
01077     } else {
01078       strncpy(atm.name, m->atomNames.name(atom->nameindex), sizeof(atm.name)-1);
01079     }
01080 
01081     strcpy(atm.type, m->atomTypes.name(atom->typeindex));
01082     strcpy(atm.resname, m->resNames.name(atom->resnameindex));
01083     atm.resid = atom->resid;
01084     strcpy(atm.chain, m->chainNames.name(atom->chainindex));
01085     strcpy(atm.segid, m->segNames.name(atom->segnameindex));
01086     strcpy(atm.insertion, atom->insertionstr);
01087     strcpy(atm.altloc, m->altlocNames.name(atom->altlocindex));
01088     atm.atomicnumber = atom->atomicnumber;
01089     atm.occupancy = occupancy[i];
01090     atm.bfactor = beta[i];
01091     atm.mass = mass[i]; 
01092     atm.charge = charge[i]; 
01093     atm.radius = radius[i];
01094 
01095     atomindexmap[i] = j; // build index map for bond/angle/dihedral/cterms
01096 
01097     j++;
01098   }
01099 
01100   // check that the selection size matches numatoms
01101   if (j != numatoms) {
01102     msgErr 
01103       << "MolFilePlugin: Internal error, selection size (" << j 
01104       << ") doesn't match numatoms (" << numatoms << ")" << sendmsg;
01105     free (atoms);
01106     return MOLFILE_ERROR;
01107   }
01108 
01109   // set optional data fields we're providing
01110   int optflags = MOLFILE_NOOPTIONS; // initialize optflags
01111 
01112   if (m->test_dataset_flag(BaseMolecule::INSERTION))
01113     optflags |= MOLFILE_INSERTION;
01114 
01115   if (m->test_dataset_flag(BaseMolecule::OCCUPANCY))
01116     optflags |= MOLFILE_OCCUPANCY;
01117 
01118   if (m->test_dataset_flag(BaseMolecule::BFACTOR))
01119     optflags |= MOLFILE_BFACTOR;
01120 
01121   if (m->test_dataset_flag(BaseMolecule::MASS))
01122     optflags |= MOLFILE_MASS;
01123 
01124   if (m->test_dataset_flag(BaseMolecule::CHARGE))
01125     optflags |= MOLFILE_CHARGE;
01126 
01127   if (m->test_dataset_flag(BaseMolecule::RADIUS))
01128     optflags |= MOLFILE_RADIUS;
01129   
01130   if (m->test_dataset_flag(BaseMolecule::ALTLOC))
01131     optflags |= MOLFILE_ALTLOC;
01132 
01133   if (m->test_dataset_flag(BaseMolecule::ATOMICNUMBER))
01134     optflags |= MOLFILE_ATOMICNUMBER;
01135 
01136   // Build and save a bond list if this plugin has a write_bonds() callback.
01137   // Bonds must be specified to the plugin before write_structure() is called.
01138   // Only store bond information if it was either set by the user or 
01139   // by loading from other files.  We don't save auto-generated bond info 
01140   // by default anymore.  We consider auto-generated bonds worth saving
01141   // only if the user has customized other bond properties.
01142   if (can_write_bonds() && 
01143       (m->test_dataset_flag(BaseMolecule::BONDS) ||
01144        m->test_dataset_flag(BaseMolecule::BONDORDERS) ||
01145        m->test_dataset_flag(BaseMolecule::BONDTYPES))) {
01146     ResizeArray<int> bondfrom, bondto; 
01147     ResizeArray<float> bondorder;
01148     ResizeArray<int> bondtype;
01149     ResizeArray<char *>bondtypename;
01150     
01151     float *bondorderptr=NULL;
01152     int *bondtypeptr=NULL;
01153     char **bondtypenameptr=NULL;
01154     int numbondtypes=0;
01155     
01156     for (i=0; i<m->nAtoms; i++) {
01157       const MolAtom *atom = m->atom(i);
01158       ptrdiff_t bfmap = atomindexmap[i] + 1; // 1-based mapped atom index
01159 
01160       for (k=0; k<atom->bonds; k++) {
01161         int bto = atom->bondTo[k];
01162         int btmap = atomindexmap[bto] + 1; // 1-based mapped atom index
01163  
01164         // add 1-based bonds to 'on' atoms to the bond list
01165         if (bfmap > 0 && btmap > bfmap) {
01166           bondfrom.append(bfmap);
01167           bondto.append(btmap);
01168           bondorder.append(m->getbondorder(i, k)); 
01169           bondtype.append(m->getbondtype(i, k));
01170         }
01171       }
01172     } 
01173 
01174     // only store bond orders if they were set by the user or read in 
01175     // from files 
01176     if (m->test_dataset_flag(BaseMolecule::BONDORDERS))
01177       bondorderptr=&bondorder[0];
01178     else 
01179       bondorderptr=NULL; // no bond orders provided
01180 
01181     numbondtypes = m->bondTypeNames.num();
01182     for (i=0; i < numbondtypes; i++) {
01183       bondtypename.append((char *)m->bondTypeNames.name(i));
01184     }
01185     if (numbondtypes > 0)
01186       bondtypenameptr = &bondtypename[0];
01187     else 
01188       bondtypenameptr = NULL;
01189     
01190     // only store bond types if they were set by the user or read in 
01191     // from files 
01192     if (m->test_dataset_flag(BaseMolecule::BONDTYPES))
01193       bondtypeptr= &bondtype[0];
01194     else 
01195       bondtypeptr=NULL; // no bond types provided
01196 
01197 #if vmdplugin_ABIVERSION >= 15
01198     if (plugin->write_bonds(wv, bondfrom.num(), &bondfrom[0], &bondto[0], 
01199                             bondorderptr, bondtypeptr, 
01200                             numbondtypes, bondtypenameptr)) {
01201 #else
01202     if (plugin->write_bonds(wv, bondfrom.num(), &bondfrom[0], &bondto[0], 
01203                             bondorderptr)) {
01204 #endif
01205       free(atoms);
01206       free(atomindexmap);
01207       return MOLFILE_ERROR;
01208     }
01209   }
01210 
01211   // Only write angle info if all atoms are selected.
01212   // It's not clear whether there's a point in trying to preserve this
01213   // kind of information when writing out a sub-structure from an 
01214   // atom selection.
01215   if (can_write_angles() && m->test_dataset_flag(BaseMolecule::ANGLES)) {
01216     ResizeArray<int> angles;
01217     ResizeArray<int> angleTypes;
01218     ResizeArray<const char *>angleTypeNames;
01219     ResizeArray<int> dihedrals;
01220     ResizeArray<int> dihedralTypes;
01221     ResizeArray<const char *>dihedralTypeNames;
01222     ResizeArray<int> impropers;
01223     ResizeArray<int> improperTypes;
01224     ResizeArray<const char *>improperTypeNames;
01225     ResizeArray<int> cterms;
01226 
01227     int numangles = m->num_angles();
01228     int numdihedrals = m->num_dihedrals();
01229     int numimpropers = m->num_impropers();
01230 
01231     // generate packed arrays with 1-based indexing
01232     for (i=0; i<numangles; i++) {
01233       ptrdiff_t i3addr = i*3L;
01234       int idx0 = atomindexmap[m->angles[i3addr    ]];
01235       int idx1 = atomindexmap[m->angles[i3addr + 1]];
01236       int idx2 = atomindexmap[m->angles[i3addr + 2]];
01237       if ((idx0 >= 0) && (idx1 >= 0) && (idx2 >= 0)) {
01238         angles.append3(idx0+1, idx1+1, idx2+1); // 1-based indices
01239 #if vmdplugin_ABIVERSION >= 16
01240         if (m->test_dataset_flag(BaseMolecule::ANGLETYPES)) {
01241           angleTypes.append(m->get_angletype(i));
01242         }
01243 #endif
01244       } 
01245     } 
01246     for (i=0; i<numdihedrals; i++) {
01247       ptrdiff_t i4addr = i*4L;
01248       int idx0 = atomindexmap[m->dihedrals[i4addr    ]];
01249       int idx1 = atomindexmap[m->dihedrals[i4addr + 1]];
01250       int idx2 = atomindexmap[m->dihedrals[i4addr + 2]];
01251       int idx3 = atomindexmap[m->dihedrals[i4addr + 3]];
01252       if ((idx0 >= 0) && (idx1 >= 0) && (idx2 >= 0) && (idx3 >= 0)) {
01253         dihedrals.append4(idx0+1, idx1+1, idx2+1, idx3+1); // 1-based indices
01254 #if vmdplugin_ABIVERSION >= 16
01255         if (m->test_dataset_flag(BaseMolecule::ANGLETYPES)) {
01256           dihedralTypes.append(m->get_dihedraltype(i));
01257         }
01258 #endif
01259       } 
01260     } 
01261     for (i=0; i<numimpropers; i++) {
01262       ptrdiff_t i4addr = i*4L;
01263       int idx0 = atomindexmap[m->impropers[i4addr    ]];
01264       int idx1 = atomindexmap[m->impropers[i4addr + 1]];
01265       int idx2 = atomindexmap[m->impropers[i4addr + 2]];
01266       int idx3 = atomindexmap[m->impropers[i4addr + 3]];
01267       if ((idx0 >= 0) && (idx1 >= 0) && (idx2 >= 0) && (idx3 >= 0)) {
01268         impropers.append4(idx0+1, idx1+1, idx2+1, idx3+1); // 1-based indices
01269 #if vmdplugin_ABIVERSION >= 16
01270         if (m->test_dataset_flag(BaseMolecule::ANGLETYPES)) {
01271           improperTypes.append(m->get_impropertype(i));
01272         }
01273 #endif
01274       } 
01275     } 
01276 
01277     int ctermcnt=0;
01278     int *ctermlist=NULL;
01279     if (m->test_dataset_flag(BaseMolecule::CTERMS)) {
01280       int numcterms = m->num_cterms();
01281       for (i=0; i<numcterms; i++) {
01282         int goodcount=0;
01283         for (j=0; j<8; j++) {
01284           if (atomindexmap[m->cterms[i*8L + j]] >= 0)
01285             goodcount++; 
01286         }
01287         if (goodcount == 8) {
01288           ctermcnt++;
01289           for (j=0; j<8; j++) {
01290             // 1-based atom index map
01291             cterms.append(atomindexmap[m->cterms[i*8L + j]]+1);
01292           }
01293         }
01294       }
01295       if (ctermcnt > 0)
01296         ctermlist = &cterms[0];
01297     }
01298 
01299 #if vmdplugin_ABIVERSION >= 16
01300     // copy all names.
01301     if (m->test_dataset_flag(BaseMolecule::ANGLETYPES)) {
01302       for (i=0; i < m->angleTypeNames.num(); i++)
01303         angleTypeNames.append(m->angleTypeNames.name(i));
01304       for (i=0; i < m->dihedralTypeNames.num(); i++)
01305         dihedralTypeNames.append(m->dihedralTypeNames.name(i));
01306       for (i=0; i < m->improperTypeNames.num(); i++)
01307         improperTypeNames.append(m->improperTypeNames.name(i));
01308     }
01309 #endif
01310     
01311     int anglescnt         = angles.num()/3;
01312     int *anglelist        = (anglescnt    > 0) ? &angles[0] : NULL;
01313 #if vmdplugin_ABIVERSION >= 16
01314     int *angletypelist    = (angleTypes.num() > 0) ? &angleTypes[0] : NULL;
01315     int angletypecnt      = angleTypeNames.num();
01316     const char **anglenmlist = (angletypecnt>0) ? &angleTypeNames[0] : NULL;
01317 #endif
01318 
01319     int dihedralscnt      = dihedrals.num()/4;
01320     int *dihedrallist     = (dihedralscnt > 0) ? &dihedrals[0] : NULL;
01321 #if vmdplugin_ABIVERSION >= 16
01322     int *dihedraltypelist = (dihedralTypes.num() > 0) ? &dihedralTypes[0] : NULL;
01323     int dihedraltypecnt   = dihedralTypeNames.num();
01324     const char **dihedralnmlist = (dihedraltypecnt>0) ? &dihedralTypeNames[0] : NULL;
01325 #endif
01326 
01327     int improperscnt      = impropers.num()/4;
01328     int *improperlist     = (improperscnt > 0) ? &impropers[0] : NULL;
01329 #if vmdplugin_ABIVERSION >= 16
01330     int *impropertypelist = (improperTypes.num() > 0) ? &improperTypes[0] : NULL;
01331     int impropertypecnt   = improperTypeNames.num();
01332     const char **impropernmlist = (impropertypecnt>0) ? &improperTypeNames[0] : NULL;
01333     if (plugin->write_angles(wv, anglescnt, anglelist, angletypelist, angletypecnt, 
01334                              anglenmlist, dihedralscnt, dihedrallist, dihedraltypelist,
01335                              dihedraltypecnt, dihedralnmlist, improperscnt, 
01336                              improperlist, impropertypelist, impropertypecnt, 
01337                              impropernmlist, ctermcnt, ctermlist, 0, 0)) {
01338       free(atoms);
01339       free(atomindexmap);
01340       return MOLFILE_ERROR;
01341     }
01342 #else
01343     if (plugin->write_angles(wv, anglescnt, anglelist, NULL,
01344                              dihedralscnt, dihedrallist, NULL,
01345                              improperscnt, improperlist, NULL,
01346                              ctermcnt, ctermlist, 0, 0, NULL)) {
01347       free(atoms);
01348       free(atomindexmap);
01349       return MOLFILE_ERROR;
01350     }
01351 #endif
01352   }
01353 
01354   // write the structure
01355   if (plugin->write_structure(wv, optflags, atoms)) {
01356     free(atoms);
01357     free(atomindexmap);
01358     return MOLFILE_ERROR;
01359   }
01360 
01361   free(atoms);
01362   free(atomindexmap);
01363 
01364   return MOLFILE_SUCCESS;
01365 }
01366  
01367 int MolFilePlugin::write_timestep(const Timestep *ts, const int *on) {
01368   // it isn't an error if this file format doesn't write timesteps
01369   if (!can_write_timesteps()) 
01370     return MOLFILE_SUCCESS; 
01371 
01372   molfile_timestep_t mol_ts;
01373   memset(&mol_ts, 0, sizeof(molfile_timestep_t));
01374   mol_ts.A = ts->a_length;
01375   mol_ts.B = ts->b_length;
01376   mol_ts.C = ts->c_length;
01377   mol_ts.alpha = ts->alpha;
01378   mol_ts.beta = ts->beta;
01379   mol_ts.gamma = ts->gamma;
01380   mol_ts.physical_time = ts->physical_time;
01381 
01382   if (!on) {
01383     mol_ts.coords = ts->pos;
01384     mol_ts.velocities = ts->vel;
01385     return plugin->write_timestep(wv, &mol_ts);
01386   }
01387 
01388   float *coords = new float[3L*numatoms];
01389   float *vel = NULL;
01390   if (ts->vel) vel = new float[3L*numatoms];
01391   ptrdiff_t j=0;
01392   for (int i=0; i<ts->num; i++) {
01393     if (on[i]) {
01394       if (on && !on[i]) continue;
01395       // check that the selection doesn't contain too many atoms
01396       if (j >= 3L*numatoms) {
01397         msgErr << "MolFilePlugin::write_timestep: Internal error" << sendmsg;
01398         msgErr << "Selection size exceeds numatoms (" << numatoms << ")" 
01399                << sendmsg;
01400         delete [] coords;
01401         delete vel;
01402         return MOLFILE_ERROR;
01403       }
01404       coords[j  ] = ts->pos[3L*i];
01405       coords[j+1] = ts->pos[3L*i+1];
01406       coords[j+2] = ts->pos[3L*i+2];
01407       if (ts->vel) {
01408         vel[j  ] = ts->vel[3L*i  ];
01409         vel[j+1] = ts->vel[3L*i+1];
01410         vel[j+2] = ts->vel[3L*i+2];
01411       }
01412       j += 3;
01413     }
01414   }
01415   // check that the selection size matches numatoms
01416   if (j != 3L*numatoms) {
01417     msgErr << "MolFilePlugin::write_timestep: Internal error" << sendmsg;
01418     msgErr << "selection size (" << j << ") doesn't match numatoms (" 
01419            << numatoms << ")" << sendmsg;
01420     delete [] coords;
01421     return MOLFILE_ERROR;
01422   }
01423   mol_ts.coords = coords;
01424   mol_ts.velocities = vel;
01425   int rc = plugin->write_timestep(wv, &mol_ts);
01426   delete [] coords;
01427   delete vel;
01428   return rc;
01429 }
01430  
01431 int MolFilePlugin::read_rawgraphics(Molecule *m, Scene *sc) {
01432   if (!rv || !can_read_graphics()) return MOLFILE_ERROR;
01433   const molfile_graphics_t *graphics = NULL;
01434   int nelem = -1;
01435   if (plugin->read_rawgraphics(rv, &nelem, &graphics)) return MOLFILE_ERROR;
01436   msgInfo << "Reading " << nelem << " graphics elements..." << sendmsg;
01437   MoleculeGraphics *mg = m->moleculeGraphics();
01438   for (int i=0; i<nelem; i++) {
01439     const float *data = graphics[i].data;
01440     switch (graphics[i].type) {
01441       case MOLFILE_POINT: 
01442         mg->add_point(data); 
01443         break;
01444       case MOLFILE_TRIANGLE:
01445         mg->add_triangle(data, data+3, data+6);
01446         break;
01447       case MOLFILE_TRINORM:
01448         {
01449           const float *ndata;
01450           // next element must be the norms
01451           if (i+1 >= nelem || graphics[i+1].type != MOLFILE_NORMS) {
01452             msgErr << "Invalid rawgraphics: NORMS must follow TRINORM."
01453               << sendmsg;
01454             return MOLFILE_ERROR;
01455           }
01456           ++i;
01457           ndata = graphics[i].data;
01458           mg->add_trinorm(data, data+3, data+6, ndata, ndata+3, ndata+6);
01459         }
01460         break;
01461       case MOLFILE_TRICOLOR: 
01462         {
01463           const float *ndata, *cdata;
01464           // next element must be the norms
01465           if (i+1 >= nelem || graphics[i+1].type != MOLFILE_NORMS) {
01466             msgErr << "Invalid rawgraphics: NORMS must follow TRINORM."
01467               << sendmsg;
01468             return MOLFILE_ERROR;
01469           }
01470           ++i;
01471           ndata = graphics[i].data;
01472           // next element must be the vertex colors
01473           if (i+1 >= nelem || graphics[i+1].type != MOLFILE_COLOR) {
01474             msgErr << "Invalid rawgraphics: NORMS and COLOR must fullow TRICOLOR."
01475               << sendmsg;
01476             return MOLFILE_ERROR;
01477           }
01478           ++i;
01479           cdata = graphics[i].data;
01480           mg->add_tricolor(data, data+3, data+6, ndata, ndata+3, ndata+6,
01481               sc->nearest_index(cdata[0], cdata[1], cdata[2]), 
01482               sc->nearest_index(cdata[3], cdata[4], cdata[5]), 
01483               sc->nearest_index(cdata[6], cdata[7], cdata[8]));
01484         }
01485         break;
01486       case MOLFILE_LINE:
01487         mg->add_line(data, data+3, graphics[i].style, (int)graphics[i].size);
01488         break;
01489       case MOLFILE_CYLINDER:
01490         mg->add_cylinder(data, data+3, graphics[i].size, graphics[i].style, 0);
01491         break;
01492       case MOLFILE_CAPCYL:
01493         mg->add_cylinder(data, data+3, graphics[i].size, graphics[i].style, 1);
01494         break;
01495       case MOLFILE_CONE:
01496         mg->add_cone(data, data+3, graphics[i].size, data[6], graphics[i].style);
01497         break;
01498       case MOLFILE_SPHERE:
01499         mg->add_sphere(data, graphics[i].size, graphics[i].style);
01500         break;
01501       case MOLFILE_TEXT:
01502         {
01503           char text[24];
01504           strncpy(text, (char *)data+3, 24);
01505           text[23] = '\0';
01506           mg->add_text(data, text, graphics[i].size, 1.0f);
01507         }
01508         break;
01509       case MOLFILE_COLOR:
01510         mg->use_color(sc->nearest_index(data[0], data[1], data[2]));
01511         break;
01512       case MOLFILE_NORMS:
01513         msgErr << "Invalid rawgraphics: NORMS must follow TRINORM." << sendmsg;
01514         return MOLFILE_ERROR;
01515         break;
01516       default:
01517         msgErr << "Invalid rawgraphics: unknown type " << graphics[i].type
01518                << sendmsg;
01519     }
01520   }
01521 
01522   return MOLFILE_SUCCESS;
01523 }
01524 
01525 
01526 int MolFilePlugin::read_volumetric(Molecule *m, int nsets, const int *setids) {
01527   molfile_volumetric_t *metadata; // fetch metadata from file
01528 
01529   int setsinfile = 0;
01530   plugin->read_volumetric_metadata(rv, &setsinfile, &metadata);
01531 
01532   // Get datasets specified in setids
01533   int n;
01534   int *sets;
01535   if (nsets < 0) {
01536     n = setsinfile;
01537     sets = new int[n];
01538     for (int i=0; i<n; i++) sets[i] = i;
01539   } else {
01540     n = nsets;
01541     sets = new int [n];
01542     for (int i=0; i<n; i++) sets[i] = setids[i];
01543   }
01544 
01545   for (int i=0; i< n; i++) {
01546     if (sets[i] < 0 || sets[i] >= setsinfile) {
01547       msgErr << "Bogus setid passed to read_volumetric: " << sets[i]
01548              << sendmsg;
01549       continue;
01550     }  
01551 
01552     const molfile_volumetric_t *v = metadata+sets[i];
01553     ptrdiff_t size = ptrdiff_t(v->xsize) * ptrdiff_t(v->ysize) * ptrdiff_t(v->zsize);
01554 
01555     char *dataname = stringdup(v->dataname);
01556     if (_filename) {
01557       // prepend the filename to the dataname; otherwise it's hard to tell
01558       // multiple data sets apart in the GUI.  This should be done here,
01559       // within VMD, rather than within each plugin because otherwise 
01560       // different plugins will end up exhibiting different behavior with
01561       // regard to naming their datasets.
01562       //
01563       // XXX The breakup_filename command uses forward slashes only and
01564       // is therefore Unix-specific.  Also, to avoid super-long dataset
01565       // names I'm going to use just the 'basename' part of the file.
01566       // It's easier just to code a correct version of what I want here.
01567       char sep = 
01568 #ifdef WIN32
01569         '\\'
01570 #else
01571         '/'
01572 #endif
01573         ;
01574       const char *basename = strrchr(_filename, sep);
01575       if (!basename) {
01576         basename = _filename;
01577       } else {
01578         basename++; // skip the separator
01579       }
01580       char *tmp = new char[strlen(dataname)+5+strlen(basename)];
01581       sprintf(tmp, "%s : %s", basename, dataname);
01582       delete [] dataname;
01583       dataname = tmp;
01584     }
01585 
01586 
01587 #if vmdplugin_ABIVERSION > 16
01588     if (plugin->read_volumetric_data_ex != NULL) {
01589 msgInfo << "Loading voumetric data using ABI 17+ ... " << sendmsg;
01590       molfile_volumetric_readwrite_t rwparms;
01591       memset(&rwparms, 0, sizeof(rwparms));
01592  
01593       rwparms.setidx = sets[i];
01594       if (v->has_scalar) 
01595         rwparms.scalar = new float[size];
01596       if (v->has_gradient) 
01597         rwparms.gradient = new float[3L*size];
01598 #if 0
01599       if (v->has_variance) 
01600         rwparms.variance = new float[size];
01601       if (v->has_color == ...) 
01602         rwparms.rgb3f = new float[3L*size];
01603       if (v->has_color == ...) 
01604         rwparms.rgb3u = new unsigned char[3L*size];
01605 #endif
01606 
01607       if (plugin->read_volumetric_data_ex(rv, &rwparms)) {
01608         msgErr << "Error reading volumetric data set " << sets[i]+1 << sendmsg;
01609         delete [] dataname;
01610         delete [] rwparms.scalar;
01611         delete [] rwparms.gradient;
01612         delete [] rwparms.variance;
01613         delete [] rwparms.rgb3f;
01614         delete [] rwparms.rgb3u;
01615         continue;
01616       }
01617 
01618       m->add_volume_data(dataname, v->origin, v->xaxis, v->yaxis, v->zaxis, 
01619                          v->xsize, v->ysize, v->zsize, 
01620                          rwparms.scalar, rwparms.gradient, rwparms.variance);
01621       delete [] dataname;
01622     } else if (plugin->read_volumetric_data != NULL) {
01623 #endif 
01624       float *scalar=NULL, *rgb3f=NULL;
01625       scalar = new float[size];
01626       if (v->has_color) 
01627         rgb3f = new float[3L*size];
01628       if (plugin->read_volumetric_data(rv, sets[i], scalar, rgb3f)) {
01629         msgErr << "Error reading volumetric data set " << sets[i]+1 << sendmsg;
01630         delete [] dataname;
01631         delete [] scalar;
01632         delete [] rgb3f;
01633         continue;
01634       }
01635 
01636       m->add_volume_data(dataname, v->origin, v->xaxis, v->yaxis, v->zaxis, 
01637                          v->xsize, v->ysize, v->zsize, scalar);
01638       delete [] dataname;
01639       delete [] rgb3f; // have to delete if created, since we don't use yet
01640 #if vmdplugin_ABIVERSION > 16
01641     }
01642 #endif 
01643   }
01644 
01645   delete [] sets;
01646 
01647   return MOLFILE_SUCCESS;
01648 }
01649 
01650 
01651 int MolFilePlugin::read_metadata(Molecule *m) {
01652   // Fetch metadata from file
01653   molfile_metadata_t *metadata;
01654 
01655   plugin->read_molecule_metadata(rv, &metadata);
01656 
01657   m->record_database(metadata->database, metadata->accession);
01658   if (metadata->remarks != NULL) 
01659     m->record_remarks(metadata->remarks);
01660   else 
01661     m->record_remarks("");
01662 
01663   return MOLFILE_SUCCESS;
01664 }
01665 
01666 
01667 int MolFilePlugin::read_qm_data(Molecule *mol) {
01668   // Fetch metadata from file.
01669   // It provides us with a bunch of sizes for the arrays
01670   // that we have to allocate and provide to read_qm_rundata().
01671   // Note that this probably should be done 
01672   molfile_qm_metadata_t metadata;
01673 
01674   // check for failures while parsing metadata and bail out if
01675   // an error occurs.
01676   if (plugin->read_qm_metadata(rv, &metadata) != MOLFILE_SUCCESS)
01677     return MOLFILE_ERROR;
01678 
01679   // If the plugin didn't provide the number of atoms
01680   // (e.g. because it only read the basis set) we set it
01681   // to the number of atoms in the molecule.
01682   // XXX This is kind of dangerous: If you load a basis set
01683   // on top of a multimillion atom structure then a basis
01684   // will be added to each of the atoms whic cost a lot of
01685   // memory!
01686   if (!numatoms) numatoms = mol->nAtoms;
01687   mol->qm_data = new QMData(numatoms,
01688                             metadata.num_basis_funcs,
01689                             metadata.num_shells,
01690                             metadata.wavef_size);
01691   //mol->qm_data->nintcoords = metadata.nintcoords;
01692 
01693   // We need to keep a pointer to the QMData object
01694   // to be used in next() when we are sorting the
01695   // wavefunction coefficients for the current timestep.
01696   qm_data = mol->qm_data;
01697 
01698   molfile_qm_t *qmdata = (molfile_qm_t *) calloc(1, sizeof(molfile_qm_t));
01699 
01700   // Allocate memory for the arrays:
01701   if (metadata.num_basis_atoms) {
01702     qmdata->basis.num_shells_per_atom = new int[metadata.num_basis_atoms];
01703     qmdata->basis.atomic_number       = new int[metadata.num_basis_atoms];
01704     qmdata->basis.num_prim_per_shell  = new int[metadata.num_shells];
01705     qmdata->basis.basis            = new float[2L*metadata.num_basis_funcs];
01706     qmdata->basis.shell_types      = new int[metadata.num_shells];
01707     qmdata->basis.angular_momentum = new int[3L*metadata.wavef_size];
01708   }
01709 
01710   if (metadata.nimag)
01711     qmdata->hess.imag_modes = new int[metadata.nimag];
01712 
01713   if (metadata.have_normalmodes) {
01714     qmdata->hess.normalmodes = new float[metadata.ncart*metadata.ncart];
01715     qmdata->hess.wavenumbers = new float[metadata.ncart];
01716     qmdata->hess.intensities = new float[metadata.ncart];
01717   }
01718 
01719   if (metadata.have_carthessian)
01720     qmdata->hess.carthessian = new double[metadata.ncart*metadata.ncart];
01721 
01722   if (metadata.have_inthessian)
01723     qmdata->hess.inthessian = new double[metadata.nintcoords*metadata.nintcoords];
01724 
01725 //   if (metadata.have_esp) {
01726 //     qmdata->run.esp_charges = new double[numatoms];
01727 //   }
01728 //   else
01729 //     qmdata->run.esp_charges = NULL;
01730 
01731   // All necessary arrays are allocated.
01732   // Now get the data from the plugin:
01733   plugin->read_qm_rundata(rv, qmdata);
01734 
01735   // Copy data from molfile_plugin structs into VMD's data structures.
01736 
01737   if (metadata.have_sysinfo) {
01738     //mol->qm_data->num_orbitals_A = qmdata->run.num_orbitals_A;
01739     //mol->qm_data->num_orbitals_B = qmdata->run.num_orbitals_B;
01740     mol->qm_data->nproc  = qmdata->run.nproc;
01741     mol->qm_data->memory = qmdata->run.memory;
01742 
01743     // We need to translate between the macros used in the plugins
01744     // and the one ones in VMD, here so that they can be independent
01745     // and we don't have to include molfile_plugin.h anywhere else
01746     // in VMD.
01747     switch (qmdata->run.scftype) {
01748     case MOLFILE_SCFTYPE_NONE:
01749       mol->qm_data->scftype = SCFTYPE_NONE;
01750       break;
01751     case MOLFILE_SCFTYPE_RHF:
01752       mol->qm_data->scftype = SCFTYPE_RHF;
01753       break;
01754     case MOLFILE_SCFTYPE_UHF:
01755       mol->qm_data->scftype = SCFTYPE_UHF;
01756       break;
01757     case MOLFILE_SCFTYPE_ROHF:
01758       mol->qm_data->scftype = SCFTYPE_ROHF;
01759       break;
01760     case MOLFILE_SCFTYPE_GVB:
01761       mol->qm_data->scftype = SCFTYPE_GVB;
01762       break;
01763     case MOLFILE_SCFTYPE_MCSCF:
01764       mol->qm_data->scftype = SCFTYPE_MCSCF;
01765       break;
01766     case MOLFILE_SCFTYPE_FF:
01767       mol->qm_data->scftype = SCFTYPE_FF;
01768       break;
01769     default:
01770       mol->qm_data->scftype = SCFTYPE_UNKNOWN;
01771     }
01772 
01773     switch (qmdata->run.runtype) {
01774     case MOLFILE_RUNTYPE_ENERGY:
01775       mol->qm_data->runtype = RUNTYPE_ENERGY;
01776       break;
01777     case MOLFILE_RUNTYPE_OPTIMIZE:
01778       mol->qm_data->runtype = RUNTYPE_OPTIMIZE;
01779       break;
01780     case MOLFILE_RUNTYPE_SADPOINT:
01781       mol->qm_data->runtype = RUNTYPE_SADPOINT;
01782       break;
01783     case MOLFILE_RUNTYPE_HESSIAN:
01784       mol->qm_data->runtype = RUNTYPE_HESSIAN;
01785       break;
01786     case MOLFILE_RUNTYPE_SURFACE:
01787       mol->qm_data->runtype = RUNTYPE_SURFACE;
01788       break;
01789     case MOLFILE_RUNTYPE_GRADIENT:
01790       mol->qm_data->runtype = RUNTYPE_GRADIENT;
01791       break;
01792     case MOLFILE_RUNTYPE_MEX:
01793       mol->qm_data->runtype = RUNTYPE_MEX;
01794       break;
01795     case MOLFILE_RUNTYPE_DYNAMICS:
01796       mol->qm_data->runtype = RUNTYPE_DYNAMICS;
01797       break;
01798     case MOLFILE_RUNTYPE_PROPERTIES:
01799       mol->qm_data->runtype = RUNTYPE_PROPERTIES;
01800       break;
01801     default:
01802       mol->qm_data->runtype = RUNTYPE_UNKNOWN;
01803     }
01804 
01805     switch (qmdata->run.status) {
01806     case MOLFILE_QMSTATUS_OPT_CONV:
01807       mol->qm_data->status = QMSTATUS_OPT_CONV;
01808       break;
01809     case MOLFILE_QMSTATUS_OPT_NOT_CONV:
01810       mol->qm_data->status = QMSTATUS_OPT_NOT_CONV;
01811       break;
01812     case MOLFILE_QMSTATUS_SCF_NOT_CONV:
01813       mol->qm_data->status = QMSTATUS_SCF_NOT_CONV;
01814       break;
01815     case MOLFILE_QMSTATUS_FILE_TRUNCATED:
01816       mol->qm_data->status = QMSTATUS_FILE_TRUNCATED;
01817       break;
01818     default:
01819       mol->qm_data->status = QMSTATUS_UNKNOWN;
01820     }
01821 
01822     // Run data:
01823     SAFESTRNCPY(mol->qm_data->version_string, qmdata->run.version_string);
01824     SAFESTRNCPY(mol->qm_data->runtitle, qmdata->run.runtitle);
01825     SAFESTRNCPY(mol->qm_data->geometry, qmdata->run.geometry);
01826   }
01827 
01828   if (metadata.have_sysinfo) {
01829     // Initialize total charge, multiplicity, number of electrons,
01830     // number of occupied orbitals.
01831     // Note that mol->qm_data->scftyp must have been
01832     // assign before.
01833     mol->qm_data->init_electrons(mol, qmdata->run.totalcharge);
01834   }
01835 
01836   // Populate basis set data and organize them into
01837   // hierarcical data structures.
01838   if (!mol->qm_data->init_basis(mol, metadata.num_basis_atoms,
01839                                 qmdata->run.basis_string,
01840                                 qmdata->basis.basis,
01841                                 qmdata->basis.atomic_number,
01842                                 qmdata->basis.num_shells_per_atom,
01843                                 qmdata->basis.num_prim_per_shell,
01844                                 qmdata->basis.shell_types)) {
01845     msgWarn << "Incomplete basis set info in QM data."
01846            << sendmsg;
01847   }
01848 
01849   // Exponents of angular momenta in wave function
01850   if (metadata.wavef_size) {
01851     mol->qm_data->set_angular_momenta(qmdata->basis.angular_momentum);
01852   }
01853 
01854 
01855   // Hessian data:
01856   if (metadata.have_carthessian) {
01857     mol->qm_data->set_carthessian(metadata.ncart, qmdata->hess.carthessian);
01858   }
01859   
01860   if (metadata.have_inthessian) {
01861     mol->qm_data->set_inthessian(metadata.nintcoords, qmdata->hess.inthessian);
01862   }
01863 
01864   if (metadata.have_normalmodes) {
01865     mol->qm_data->set_normalmodes(metadata.ncart, qmdata->hess.normalmodes);
01866     mol->qm_data->set_wavenumbers(metadata.ncart, qmdata->hess.wavenumbers);
01867     mol->qm_data->set_intensities(metadata.ncart, qmdata->hess.intensities);
01868   }
01869 
01870   if (metadata.nimag) {
01871     mol->qm_data->set_imagmodes(metadata.nimag, qmdata->hess.imag_modes);
01872   }
01873 
01874   // Cleanup the arrays we needed to get the data from the plugin.
01875   if (metadata.num_basis_atoms) {
01876     delete [] qmdata->basis.num_shells_per_atom;
01877     delete [] qmdata->basis.atomic_number;
01878     delete [] qmdata->basis.num_prim_per_shell;
01879     delete [] qmdata->basis.basis;
01880     delete [] qmdata->basis.shell_types;
01881     delete [] qmdata->basis.angular_momentum;
01882   }
01883   delete [] qmdata->hess.carthessian;
01884   delete [] qmdata->hess.inthessian;
01885   delete [] qmdata->hess.normalmodes;
01886   delete [] qmdata->hess.wavenumbers;
01887   delete [] qmdata->hess.intensities;
01888   delete [] qmdata->hess.imag_modes;
01889   free(qmdata);
01890 
01891   return MOLFILE_SUCCESS;
01892 }
01893 
01894 
01895 int MolFilePlugin::write_volumetric(Molecule *m, int set) {
01896   if (set < 0 || set > m->num_volume_data()) {
01897     msgErr << "Bogus setid passed to write_volumetric: " << set
01898            << sendmsg;
01899     return MOLFILE_SUCCESS;
01900   } 
01901 
01902   const VolumetricData *v = m->get_volume_data(set); 
01903 
01904   molfile_volumetric_t volmeta;
01905 
01906   // SAFESTRNCPY(volmeta.dataname, v->name); 
01907   int n = ((sizeof(volmeta.dataname) < strlen(v->name)+1) ? 
01908             sizeof(volmeta.dataname) : strlen(v->name)+1);
01909   strncpy(volmeta.dataname, v->name, n);
01910 
01911   for (int i=0; i<3; i++) {
01912       volmeta.origin[i] = (float)v->origin[i];
01913       volmeta.xaxis[i] = (float)v->xaxis[i];
01914       volmeta.yaxis[i] = (float)v->yaxis[i];
01915       volmeta.zaxis[i] = (float)v->zaxis[i];
01916   }
01917   volmeta.xsize = v->xsize;
01918   volmeta.ysize = v->ysize;
01919   volmeta.zsize = v->zsize;
01920   volmeta.has_color = 0;
01921  
01922   float *datablock = v->data;
01923   float *colorblock = NULL;
01924 
01925   plugin->write_volumetric_data(wv, &volmeta, datablock, colorblock);
01926 
01927   return MOLFILE_SUCCESS;
01928 }
01929 
01930 
01931 

Generated on Thu Mar 28 02:43:31 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002