MolFilePlugin.C

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