00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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;
00083 if ((rc = plugin->read_structure(rv, &optflags, atoms))) {
00084 free(atoms);
00085 return rc;
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;
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
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
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
00131
00132
00133
00134 mass[i] = (optflags & MOLFILE_MASS) ?
00135 atom.mass : ((atomicnumber > 0) ?
00136 get_pte_mass(atomicnumber) : m->default_mass(atom.name));
00137
00138
00139
00140
00141
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
00159
00160
00161 msgErr << "MolFilePlugin: file load aborted" << sendmsg;
00162 free(atoms);
00163 return MOLFILE_ERROR;
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
00174
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
00192
00193
00194 if (nbonds == 0 && from == NULL && to == NULL) {
00195 if (autobonds)
00196 m->find_bonds_from_timestep();
00197 } else if (nbonds > 0) {
00198
00199
00200 if (bondtypename != NULL)
00201 for (i=0; i < nbondtypes; i++)
00202 m->bondTypeNames.add_name(bondtypename[i], i);
00203
00204
00205
00206
00207
00208 if (optflags & MOLFILE_BONDSSPECIAL) {
00209 if (autobonds) {
00210 m->find_unique_bonds_from_timestep();
00211 }
00212
00213
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
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
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
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
00261
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
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
00297
00298
00299 if ( (angletypes != NULL) || (dihedraltypes != NULL)
00300 || (impropertypes != NULL) )
00301 m->set_dataset_flag(BaseMolecule::ANGLETYPES);
00302
00303
00304
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
00318
00319
00320
00321
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
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;
00404 if ((rc = plugin->read_structure(rv, &optflags, atoms))) {
00405 free(atoms);
00406 return rc;
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;
00413 }
00414
00415
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
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
00466 if (!can_read_bonds())
00467 return MOLFILE_SUCCESS;
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 if (!(optflags & MOLFILE_BONDSSPECIAL) && filebonds) {
00478 int nbonds, *from, *to, nbondtypes, *bondtype;
00479 float *bondorder;
00480 char **bondtypename;
00481
00482
00483
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
00498
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);
00530 }
00531 } else {
00532
00533 return MOLFILE_SUCCESS;
00534 }
00535
00536
00537
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
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
00573
00574
00575 if ( (angletypes != NULL) || (dihedraltypes != NULL)
00576 || (impropertypes != NULL) )
00577 m->set_dataset_flag(BaseMolecule::ANGLETYPES);
00578
00579
00580
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
00594
00595
00596
00597
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
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
00668
00669 m->analyze();
00670
00671
00672 m->force_recalc(DrawMolItem::COL_REGEN | DrawMolItem::SEL_REGEN);
00673
00674
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
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;
00727 #else
00728
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(×tep, 0, sizeof(molfile_timestep_t));
00741
00742
00743
00744
00745 float *velocities = NULL;
00746
00747
00748
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
00759
00760
00761 molfile_qm_timestep_metadata_t qmmeta;
00762 if (can_read_qm_timestep_metadata()) {
00763 memset(&qmmeta, 0, sizeof(molfile_qm_timestep_metadata));
00764
00765
00766
00767 plugin->read_qm_timestep_metadata(rv, &qmmeta);
00768 }
00769
00770
00771
00772
00773
00774 timestep.A = timestep.B = timestep.C = 0.0f;
00775
00776
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;
00786 molfile_qm_timestep_t qm_timestep;
00787 memset(&qm_timestep, 0, sizeof(molfile_qm_timestep_t));
00788
00789
00790
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, ×tep, qm_metadata, &qm_timestep);
00816 } else {
00817 rc = plugin->read_next_timestep(rv, numatoms, ×tep);
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
00876
00877 int num_signa_ts = 0;
00878 wavef_signa_t *signa_ts = NULL;
00879 for (i=0; i<qmmeta.num_wavef; i++) {
00880
00881
00882
00883
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
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
00947
00948
00949
00950
00951 if (m->numframes()>0) {
00952
00953 Timestep *prevts = m->get_frame(m->numframes()-1);
00954
00955 msgInfo << "sort frame " << m->numframes() << sendmsg;
00956
00957
00958
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
01000
01001
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
01021
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
01037 for (i=0, j=0; i<m->nAtoms; i++) {
01038
01039 if (on && !on[i])
01040 continue;
01041
01042
01043
01044
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
01059
01060
01061 char name[6], *nameptr;
01062 name[0] = ' ';
01063 strncpy(name+1, (m->atomNames).name(atom->nameindex), 4);
01064 name[5] = '\0';
01065
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;
01096
01097 j++;
01098 }
01099
01100
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
01110 int optflags = MOLFILE_NOOPTIONS;
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
01137
01138
01139
01140
01141
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;
01159
01160 for (k=0; k<atom->bonds; k++) {
01161 int bto = atom->bondTo[k];
01162 int btmap = atomindexmap[bto] + 1;
01163
01164
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
01175
01176 if (m->test_dataset_flag(BaseMolecule::BONDORDERS))
01177 bondorderptr=&bondorder[0];
01178 else
01179 bondorderptr=NULL;
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
01191
01192 if (m->test_dataset_flag(BaseMolecule::BONDTYPES))
01193 bondtypeptr= &bondtype[0];
01194 else
01195 bondtypeptr=NULL;
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
01212
01213
01214
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
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);
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);
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);
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
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
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
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
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
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
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
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
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
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;
01528
01529 int setsinfile = 0;
01530 plugin->read_volumetric_metadata(rv, &setsinfile, &metadata);
01531
01532
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
01558
01559
01560
01561
01562
01563
01564
01565
01566
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++;
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;
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
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
01669
01670
01671
01672 molfile_qm_metadata_t metadata;
01673
01674
01675
01676 if (plugin->read_qm_metadata(rv, &metadata) != MOLFILE_SUCCESS)
01677 return MOLFILE_ERROR;
01678
01679
01680
01681
01682
01683
01684
01685
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
01692
01693
01694
01695
01696 qm_data = mol->qm_data;
01697
01698 molfile_qm_t *qmdata = (molfile_qm_t *) calloc(1, sizeof(molfile_qm_t));
01699
01700
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
01726
01727
01728
01729
01730
01731
01732
01733 plugin->read_qm_rundata(rv, qmdata);
01734
01735
01736
01737 if (metadata.have_sysinfo) {
01738
01739
01740 mol->qm_data->nproc = qmdata->run.nproc;
01741 mol->qm_data->memory = qmdata->run.memory;
01742
01743
01744
01745
01746
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
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
01830
01831
01832
01833 mol->qm_data->init_electrons(mol, qmdata->run.totalcharge);
01834 }
01835
01836
01837
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
01850 if (metadata.wavef_size) {
01851 mol->qm_data->set_angular_momenta(qmdata->basis.angular_momentum);
01852 }
01853
01854
01855
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
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
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