00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include <stdlib.h>
00020 #include <string.h>
00021 #include <ctype.h>
00022 #include <errno.h>
00023 #include <time.h>
00024 #include <math.h>
00025
00026 #include "gaussianplugin.h"
00027 #include "periodic_table.h"
00028 #include "unit_conversion.h"
00029
00030 #define THISPLUGIN plugin
00031 #include "vmdconio.h"
00032
00033
00034
00035
00036 #ifndef GAUSSIAN_DEBUG
00037 #define GAUSSIAN_DEBUG 0
00038 #endif
00039 #define GAUSSIAN_BASIS_DEBUG 1
00040
00041 #if GAUSSIAN_DEBUG
00042 #define PRINTERR vmdcon_printf(VMDCON_ERROR, \
00043 "\n In file %s, line %d: \n %s \n \n", \
00044 __FILE__, __LINE__, strerror(errno))
00045 #else
00046 #define PRINTERR (void)(0)
00047 #endif
00048
00049
00050
00051
00052
00053 #if GAUSSIAN_DEBUG
00054 #define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE ; \
00055 else fprintf(stderr,"%s:%d %s",__FILE__, __LINE__, x)
00056 #else
00057 #define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE
00058 #endif
00059
00060
00061 #define SAFE_FREE(ptr) free(ptr); ptr=NULL
00062
00063 #define SAFE_CALLOC(ptr,type,count) \
00064 ptr = (type *)calloc(count,sizeof(type)); \
00065 if (ptr == NULL) { \
00066 PRINTERR; \
00067 return MOLFILE_ERROR; \
00068 }
00069
00070 #define UNK_SHELL -666
00071 #define SPD_D_SHELL -5
00072 #define SPD_P_SHELL -4
00073 #define SPD_S_SHELL -3
00074 #define SP_S_SHELL -2
00075 #define SP_P_SHELL -1
00076 #define S_SHELL 0
00077 #define P_SHELL 1
00078 #define D_SHELL 2
00079 #define F_SHELL 3
00080 #define G_SHELL 4
00081 #define H_SHELL 5
00082 #define I_SHELL 6
00083
00084 #define SPIN_ALPHA 0
00085 #define SPIN_BETA 1
00086
00087 #define STATUS_UNKNOWN -1
00088 #define STATUS_CONVERGED 0
00089 #define STATUS_SCF_NOT_CONV 1
00090 #define STATUS_TOO_MANY_STEPS 2
00091 #define STATUS_BROKEN_OFF 3
00092
00093 static const char *runtypes[] = {
00094 "(unknown)", "ENERGY", "OPTIMIZE", "SADPOINT", "HESSIAN",
00095 "SURFACE", "DYNAMICS", "PROPERTIES", NULL};
00096
00097 static const char *scftypes[] = {
00098 "(unknown)", "RHF", "UHF", "ROHF", "GVB", "MCSCF", "FF", NULL };
00099
00100
00101
00102
00103
00104
00107 static int parse_static_data(gaussiandata *, int *);
00108
00111 static int have_gaussian(gaussiandata *);
00112
00114 static int get_proc_mem(gaussiandata *);
00115
00117 static int get_basis_options(gaussiandata *);
00118
00120 static int get_runtitle(gaussiandata *);
00121
00123 static int get_input_structure(gaussiandata *);
00125 static int get_coordinates(FILE *file, qm_atom_t *atoms, int numatoms);
00127 static int get_int_coordinates(FILE *file, qm_atom_t *atoms, int numatoms);
00128
00130 static int get_contrl(gaussiandata *);
00131
00134 static int get_final_info(gaussiandata *);
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 static int get_basis(gaussiandata *);
00145
00146
00147
00148
00149
00150 static int get_internal_basis(gaussiandata *);
00151
00152
00153 static int shellsymm_int(char *symm);
00154
00155
00156 static int fill_basis_arrays(gaussiandata *);
00157
00158 static int read_first_frame(gaussiandata *);
00159
00160
00161
00162 static int get_traj_frame(gaussiandata *);
00163
00164
00165 static int find_traj_end(gaussiandata *);
00166
00167
00168
00169 static int get_wavefunction(gaussiandata *, qm_timestep_t *);
00170
00172 static int get_population(gaussiandata *, qm_timestep_t *);
00173
00174
00175 static void fix_fortran_exp(char *string) {
00176 while (*string) {
00177 if ( (*string == 'D') || (*string == 'd')) *string='e';
00178 ++string;
00179 }
00180 }
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198 static void *open_gaussian_read(const char *filename,
00199 const char *filetype, int *natoms) {
00200
00201 FILE *fd;
00202 gaussiandata *data;
00203
00204
00205 fd = fopen(filename, "rb");
00206 if (fd == NULL) {
00207 PRINTERR;
00208 return NULL;
00209 }
00210
00211
00212 data = (gaussiandata *) calloc(1,sizeof(gaussiandata));
00213 if (data == NULL) return NULL;
00214
00215 data->runtyp = RUNTYP_UNKNOWN;
00216 data->scftyp = SCFTYP_UNKNOWN;
00217 data->file = fd;
00218 data->file_name = strdup(filename);
00219
00220
00221
00222 if (have_gaussian(data)==TRUE) {
00223
00224
00225 if ((data->version < 19940000) || (data->version > 20040000)) {
00226 vmdcon_printf(VMDCON_ERROR,
00227 "gaussianplugin) Gaussian version %s is not "
00228 "(yet) supported. Bailing out.\n",
00229 data->version_string);
00230 free(data);
00231 return NULL;
00232 }
00233
00234
00235 if (parse_static_data(data, natoms) == FALSE) {
00236 free(data);
00237 return NULL;
00238 }
00239 }
00240 else {
00241 free(data);
00242 return NULL;
00243 }
00244
00245 return data;
00246 }
00247
00248
00249
00250
00251
00252
00253
00254
00255 static int read_gaussian_structure(void *mydata, int *optflags,
00256 molfile_atom_t *atoms)
00257 {
00258 gaussiandata *data = (gaussiandata *)mydata;
00259 qm_atom_t *cur_atom;
00260 molfile_atom_t *atom;
00261 int i = 0;
00262
00263
00264 *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS;
00265
00266 if (data->have_mulliken)
00267 *optflags |= MOLFILE_CHARGE;
00268
00269
00270
00271
00272
00273
00274
00275 cur_atom = data->initatoms;
00276
00277 for(i=0; i<data->numatoms; i++) {
00278 atom = atoms+i;
00279 strncpy(atom->name, cur_atom->type, sizeof(atom->name));
00280 strncpy(atom->type, get_pte_label(cur_atom->atomicnum), sizeof(atom->type));
00281 strncpy(atom->resname,"QM", sizeof(atom->resname));
00282 atom->resid = 1;
00283 atom->chain[0] = '\0';
00284 atom->segid[0] = '\0';
00285 atom->atomicnumber = cur_atom->atomicnum;
00286 atom->radius = get_pte_vdw_radius(cur_atom->atomicnum);
00287
00288 atom->mass = get_pte_mass(cur_atom->atomicnum);
00289 if (data->have_mulliken)
00290 atom->charge = data->qm_timestep->mulliken_charges[i];
00291
00292 cur_atom++;
00293 }
00294
00295 return MOLFILE_SUCCESS;
00296 }
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 static int read_gaussian_metadata(void *mydata,
00307 molfile_qm_metadata_t *gaussian_metadata) {
00308
00309 gaussiandata *data = (gaussiandata *)mydata;
00310
00311 gaussian_metadata->ncart = 0;
00312 gaussian_metadata->nimag = 0;
00313 gaussian_metadata->nintcoords = 0;
00314
00315
00316 gaussian_metadata->num_basis_funcs = data->num_basis_funcs;
00317 gaussian_metadata->num_basis_atoms = data->num_basis_atoms;
00318 gaussian_metadata->num_shells = data->num_shells;
00319 gaussian_metadata->wavef_size = data->wavef_size;
00320
00321 #if vmdplugin_ABIVERSION > 11
00322 gaussian_metadata->have_sysinfo = 1;
00323
00324
00325 gaussian_metadata->have_esp = 0;
00326 gaussian_metadata->have_carthessian = 0;
00327 gaussian_metadata->have_internals = 0;
00328 gaussian_metadata->have_normalmodes = FALSE;
00329 #endif
00330
00331 return MOLFILE_SUCCESS;
00332 }
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342 static int read_gaussian_rundata(void *mydata, molfile_qm_t *qm_data) {
00343
00344 gaussiandata *data = (gaussiandata *)mydata;
00345
00346 molfile_qm_basis_t *basis_data = &qm_data->basis;
00347 molfile_qm_sysinfo_t *sys_data = &qm_data->run;
00348
00349
00350 sys_data->nproc = data->nproc;
00351 sys_data->memory = data->memory;
00352 sys_data->runtype = data->runtyp;
00353 sys_data->scftype = data->scftyp;
00354 sys_data->totalcharge = data->totalcharge;
00355
00356 sys_data->num_electrons = data->num_electrons;
00357 sys_data->num_occupied_A = data->occ_orbitals_A;
00358 sys_data->num_occupied_B = data->occ_orbitals_B;
00359
00360 strncpy(sys_data->basis_string, data->basis_string,
00361 sizeof(sys_data->basis_string));
00362
00363 strncpy(sys_data->runtitle, data->runtitle, sizeof(sys_data->runtitle));
00364 strncpy(sys_data->geometry, data->geometry, sizeof(sys_data->geometry));
00365 strncpy(sys_data->version_string, data->version_string,
00366 sizeof(sys_data->version_string));
00367
00368 #if vmdplugin_ABIVERSION > 11
00369
00370 if (data->num_basis_funcs) {
00371 memcpy(basis_data->num_shells_per_atom, data->num_shells_per_atom,
00372 data->num_basis_atoms*sizeof(int));
00373 memcpy(basis_data->num_prim_per_shell, data->num_prim_per_shell,
00374 data->num_shells*sizeof(int));
00375 memcpy(basis_data->shell_symmetry, data->shell_symmetry,
00376 data->num_shells*sizeof(int));
00377 memcpy(basis_data->basis, data->basis, 2*data->num_basis_funcs*sizeof(float));
00378 memcpy(basis_data->angular_momentum, data->angular_momentum,
00379 3*data->wavef_size*sizeof(int));
00380 }
00381 #endif
00382
00383 return MOLFILE_SUCCESS;
00384 }
00385
00386
00387 #if vmdplugin_ABIVERSION > 11
00388
00389
00390
00391
00392
00393 static int read_timestep_metadata(void *mydata,
00394 molfile_timestep_metadata_t *meta) {
00395 meta->count = -1;
00396 meta->has_velocities = 0;
00397
00398 return MOLFILE_SUCCESS;
00399 }
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 static int read_qm_timestep_metadata(void *mydata,
00410 molfile_qm_timestep_metadata_t *meta) {
00411 int i, have = 0;
00412 gaussiandata *data = (gaussiandata *)mydata;
00413 #if GAUSSIAN_DEBUG
00414 vmdcon_printf(VMDCON_INFO,
00415 "gaussianplugin) read_qm_timestep_metadata(): %d/%d/%d\n",
00416 data->num_frames,
00417 data->num_frames_read,
00418 data->num_frames_sent);
00419 #endif
00420
00421 meta->count = -1;
00422 meta->has_gradient = 0;
00423
00424 if (data->num_frames_read > data->num_frames_sent) {
00425 have = 1;
00426 } else if (data->num_frames_read < data->num_frames) {
00427 #if GAUSSIAN_DEBUG
00428 vmdcon_printf(VMDCON_INFO,
00429 "gaussianplugin) Probing timestep %d\n",
00430 data->num_frames_read);
00431 #endif
00432 have = get_traj_frame(data);
00433 }
00434
00435 if (have) {
00436
00437 qm_timestep_t *cur_qm_ts = data->qm_timestep+data->num_frames_sent;
00438 #if GAUSSIAN_DEBUG
00439 vmdcon_printf(VMDCON_INFO,
00440 "gaussianplugin) Approved timestep %d\n",
00441 data->num_frames_sent);
00442 #endif
00443 meta->num_scfiter = cur_qm_ts->num_scfiter;
00444 for (i=0; (i<MAX_NUM_WAVE && i<cur_qm_ts->numwave); i++) {
00445 #if GAUSSIAN_DEBUG
00446 vmdcon_printf(VMDCON_INFO,
00447 "gaussianplugin) num_orbitals_per_wavef[%d/%d]=%d\n",
00448 i+1, cur_qm_ts->numwave, cur_qm_ts->wave[i].num_orbitals);
00449 #endif
00450 meta->num_orbitals_per_wavef[i] = cur_qm_ts->wave[i].num_orbitals;
00451 }
00452 meta->num_wavef = cur_qm_ts->numwave;
00453 meta->wavef_size = data->wavef_size;
00454 } else {
00455 meta->num_scfiter = 0;
00456 meta->num_orbitals_per_wavef[0] = 0;
00457 meta->num_wavef = 0;
00458 meta->wavef_size = 0;
00459
00460 data->end_of_trajectory = TRUE;
00461 }
00462
00463 return MOLFILE_SUCCESS;
00464 }
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476 static int read_timestep(void *mydata, int natoms,
00477 molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
00478 molfile_qm_timestep_t *qm_ts)
00479 {
00480 gaussiandata *data = (gaussiandata *)mydata;
00481 qm_atom_t *cur_atom;
00482 int i = 0;
00483 qm_timestep_t *cur_qm_ts;
00484
00485 if (data->end_of_trajectory == TRUE) return MOLFILE_ERROR;
00486
00487 #if GAUSSIAN_DEBUG
00488 vmdcon_printf(VMDCON_INFO,
00489 "gaussianplugin) Sending timestep %d\n",
00490 data->num_frames_sent);
00491 #endif
00492
00493
00494 cur_atom = data->initatoms;
00495
00496
00497 for(i=0; i<natoms; i++) {
00498 ts->coords[3*i ] = cur_atom->x;
00499 ts->coords[3*i+1] = cur_atom->y;
00500 ts->coords[3*i+2] = cur_atom->z;
00501 cur_atom++;
00502 }
00503
00504
00505 cur_qm_ts = data->qm_timestep+data->num_frames_sent;
00506
00507
00508 for (i=0; i<cur_qm_ts->num_scfiter; i++) {
00509 qm_ts->scfenergies[i] = cur_qm_ts->scfenergies[i];
00510 }
00511
00512
00513 if (cur_qm_ts->wave) {
00514 for (i=0; i<cur_qm_ts->numwave; i++) {
00515 qm_wavefunction_t *wave = &cur_qm_ts->wave[i];
00516 if (wave->wave_coeffs && wave->orb_energies) {
00517 memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs,
00518 wave->num_orbitals*data->wavef_size*sizeof(float));
00519 memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
00520 wave->num_orbitals*sizeof(float));
00521 }
00522 }
00523 }
00524
00525 if (data->runtyp == RUNTYP_ENERGY || data->runtyp == RUNTYP_HESSIAN) {
00526
00527 data->end_of_trajectory = TRUE;
00528 }
00529
00530 data->num_frames_sent++;
00531
00532 return MOLFILE_SUCCESS;
00533 }
00534 #endif
00535
00536
00537
00542 static void close_gaussian_read(void *mydata) {
00543
00544 gaussiandata *data = (gaussiandata *)mydata;
00545 int i, j;
00546 fclose(data->file);
00547
00548 free(data->file_name);
00549 free(data->initatoms);
00550 free(data->basis);
00551 free(data->shell_symmetry);
00552 free(data->num_shells_per_atom);
00553 free(data->num_prim_per_shell);
00554 free(data->mulliken_charges);
00555 free(data->internal_coordinates);
00556 free(data->wavenumbers);
00557 free(data->intensities);
00558 free(data->normal_modes);
00559 free(data->angular_momentum);
00560
00561 if (data->basis_set) {
00562 for(i=0; i<data->numatoms; i++) {
00563 for (j=0; j<data->basis_set[i].numshells; j++) {
00564 free(data->basis_set[i].shell[j].prim);
00565 }
00566 free(data->basis_set[i].shell);
00567 }
00568 free(data->basis_set);
00569 }
00570
00571 for (i=0; i<data->num_frames_read; i++) {
00572 free(data->qm_timestep[i].scfenergies);
00573 free(data->qm_timestep[i].gradient);
00574 free(data->qm_timestep[i].mulliken_charges);
00575 for (j=0; j<data->qm_timestep[i].numwave; j++) {
00576 free(data->qm_timestep[i].wave[j].wave_coeffs);
00577 free(data->qm_timestep[i].wave[j].orb_energies);
00578
00579 }
00580 free(data->qm_timestep[i].wave);
00581 }
00582 free(data->qm_timestep);
00583
00584 free(data);
00585 }
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600 static int parse_static_data(gaussiandata *data, int *natoms)
00601 {
00602
00603 if (!get_proc_mem(data)) return FALSE;
00604
00605
00606 if (!get_basis_options(data)) return FALSE;
00607
00608
00609
00610 if (!get_contrl(data)) return FALSE;
00611
00612
00613 if (!get_runtitle(data)) return FALSE;
00614
00615
00616 if (!get_input_structure(data)) return FALSE;
00617
00618 *natoms = data->numatoms;
00619
00620 read_first_frame(data);
00621
00622
00623 get_final_info(data);
00624
00625 vmdcon_printf(VMDCON_INFO, "gaussianplugin) found %d QM data frames.\n", data->num_frames);
00626 #if GAUSSIAN_DEBUG
00627 vmdcon_printf(VMDCON_INFO, "gaussianplugin) num_frames_read = %d\n", data->num_frames_read);
00628 vmdcon_printf(VMDCON_INFO, "gaussianplugin) num_frames_sent = %d\n", data->num_frames_sent);
00629 #endif
00630 return TRUE;
00631 }
00632
00633
00634
00635
00636
00637
00638
00639 static int have_gaussian(gaussiandata *data)
00640 {
00641 char word[4][MOLFILE_BUFSIZ];
00642 char buffer[BUFSIZ];
00643 char *ptr;
00644 int i = 0;
00645
00646 buffer[0] = '\0';
00647 for (i=0; i<3; i++) word[i][0] = '\0';
00648
00649
00650
00651
00652 i=0;
00653 do {
00654 GET_LINE(buffer, data->file);
00655 sscanf(buffer,"%s%s%s",word[0],word[1],word[2]);
00656 ++i;
00657 } while( (strcmp(word[0],"Entering") ||
00658 strcmp(word[1],"Gaussian") ||
00659 strcmp(word[2],"System,")) && (i<100) );
00660 if (i>=100) return FALSE;
00661 vmdcon_printf(VMDCON_INFO, "gaussianplugin) Analyzing Gaussian log file: %s\n",data->file_name);
00662
00663
00664
00665 i=0;
00666 do {
00667 GET_LINE(buffer, data->file);
00668 buffer[20] = '\0';
00669 } while ( (i<100) && strcmp(buffer,
00670 " *******************"));
00671 if (i>=100) return FALSE;
00672
00673
00674 GET_LINE(buffer, data->file);
00675 sscanf(buffer,"%s%s%s%s",word[0],word[1],
00676 word[2],word[3]);
00677
00678 data->version = atoi(word[1]) * 10000;
00679 if (data->version > 700000) {
00680 data->version += 19000000;
00681 } else {
00682 data->version += 20000000;
00683 }
00684 strcpy(data->version_string,word[2]);
00685
00686 ptr=strrchr(word[2],'-');
00687 if (ptr != NULL) {
00688
00689 ptr += 7;
00690 i = (*ptr) - 'A' + 1;
00691 data->version += i*100;
00692 ++ptr; ++ptr;
00693 data->version += atoi(ptr);
00694 }
00695
00696 vmdcon_printf(VMDCON_INFO,
00697 "gaussianplugin) Gaussian version = %s (Version code: %d)\n",
00698 data->version_string, data->version);
00699 vmdcon_printf(VMDCON_INFO,
00700 "gaussianplugin) Compiled on = %s \n", word[3]);
00701
00702
00703 return TRUE;
00704 }
00705
00706
00707
00708
00709
00710
00711
00712
00713 static int get_proc_mem(gaussiandata *data) {
00714
00715 char word[5][MOLFILE_BUFSIZ];
00716 char buffer[BUFSIZ];
00717 int nproc,maxmem;
00718 int i, link;
00719
00720 buffer[0] = '\0';
00721 for (i=0; i<3; i++) word[i][0] = '\0';
00722
00723 rewind(data->file);
00724
00725
00726 nproc = 1;
00727 maxmem = -1;
00728
00729 do {
00730 GET_LINE(buffer, data->file);
00731 sscanf(buffer,"%s%s%s%*s%s%*s%*s%*s%*s%*s%s",
00732 word[0],word[1],word[2],word[3],word[4]);
00733
00734
00735 if ((strcmp(word[0],"Leave") == 0) &&
00736 (strcmp(word[1],"Link") == 0)) {
00737 link = atoi(word[2]);
00738
00739 if (link > 1) maxmem=atoi(word[4])/128/1024;
00740 }
00741
00742
00743 if ( (strcmp(word[0],"Will") == 0) &&
00744 (strcmp(word[1],"use") == 0) &&
00745 (strcmp(word[2],"up") == 0) ) {
00746 nproc = atoi(word[3]);
00747 }
00748
00749
00750 if (((strcmp(word[0],"Standard") == 0) ||
00751 (strcmp(word[0],"Z-Matrix") == 0) ||
00752 (strcmp(word[0],"Input") == 0) ) &&
00753 (strcmp(word[1],"orientation:") == 0)) {
00754
00755 maxmem=0;
00756 }
00757 } while (maxmem < 0);
00758
00759
00760 data->nproc = nproc;
00761 data->memory = maxmem;
00762 if (maxmem)
00763 vmdcon_printf(VMDCON_INFO,
00764 "gaussianplugin) Gaussian used %2d SMP process(es), "
00765 "% 6d Mbytes of memory \n", nproc, maxmem);
00766
00767 return TRUE;
00768 }
00769
00770
00771
00772
00773
00774
00775
00776 static int get_basis_options(gaussiandata *data) {
00777
00778 char word[5][MOLFILE_BUFSIZ];
00779 char buffer[BUFSIZ], *ptr;
00780 int i = 0, nfunc, nprim, nume_a, nume_b;
00781
00782 buffer[0] = '\0';
00783 for (i=0; i<3; i++) word[i][0] = '\0';
00784
00785
00786 rewind(data->file);
00787
00788
00789 nume_a=-1;
00790 nfunc=-1;
00791 nprim=-1;
00792 do {
00793 GET_LINE(buffer, data->file);
00794 i=sscanf(buffer,"%s%s%s",word[0],word[1],word[2]);
00795 if (i==3) {
00796
00797 if ( (strcmp(word[0],"Standard") == 0) &&
00798 (strcmp(word[1],"basis:") == 0) ) {
00799
00800 ptr = &buffer[0] + strlen(buffer) - 1;
00801 while(*ptr==' ') --ptr;
00802 *ptr='\0';
00803 strncpy(data->gbasis, word[2], 10);
00804 strncpy(data->basis_string, buffer+17, MOLFILE_BUFSIZ);
00805
00806
00807 ptr=data->gbasis;
00808 for (;*ptr;++ptr) *ptr=toupper(*ptr);
00809
00810 } else if ( (strcmp(word[0],"General") == 0) &&
00811 (strcmp(word[1],"basis") == 0) ) {
00812
00813
00814
00815 ptr = &buffer[0] + strlen(buffer) - 1;
00816 while(*ptr==' ') --ptr;
00817 *ptr='\0';
00818
00819 strncpy(data->gbasis, "GEN", 4);
00820 strncpy(data->basis_string, buffer, MOLFILE_BUFSIZ);
00821
00822 } else if ( (strcmp(word[0],"AO") == 0)
00823 && (strcmp(word[1],"basis") == 0) ) {
00824
00825
00826 int numcenter, numgauss, numbasis, numshell;
00827 numcenter=numgauss=numbasis=numshell=0;
00828 GET_LINE(buffer, data->file);
00829 do {
00830 int numprim=0;
00831 ++numcenter;
00832 do {
00833
00834
00835 GET_LINE(buffer, data->file);
00836 i=sscanf(buffer,"%s%d",word[0], &numprim);
00837 if (i==2) {
00838 ++numshell;
00839 if (strcmp(word[0],"S") == 0) {
00840 numbasis += 1;
00841 numgauss += numprim;
00842 } else if (strcmp(word[0],"P") == 0) {
00843 numbasis += 3;
00844 numgauss += numprim;
00845 } else if (strcmp(word[0],"SP") == 0) {
00846 numbasis += 1+3;
00847 numgauss += 2*numprim;
00848 numshell += 1;
00849 } else if (strcmp(word[0],"D") == 0) {
00850 numbasis += 6;
00851 numgauss += numprim;
00852 } else if (strcmp(word[0],"SPD") == 0) {
00853 numbasis += 1+3+6;
00854 numgauss += 3*numprim;
00855 numshell += 2;
00856 } else if (strcmp(word[0],"F") == 0) {
00857 numbasis += 10;
00858 numgauss += numprim;
00859 } else if (strcmp(word[0],"G") == 0) {
00860 numbasis += 15;
00861 numgauss += numprim;
00862 } else {
00863 vmdcon_printf(VMDCON_ERROR, "gaussianplugin) support for %s-"
00864 "shells is not yet programmed.\n", word[0]);
00865 return MOLFILE_ERROR;
00866 }
00867 } else if (i==1) {
00868 if (strcmp(word[0],"****") == 0) {
00869 break;
00870 }
00871 } else {
00872 return MOLFILE_ERROR;
00873 }
00874
00875 for (i=0; i < numprim; ++i)
00876 GET_LINE(buffer, data->file);
00877
00878 #if GAUSSIAN_DEBUG
00879 vmdcon_printf(VMDCON_INFO, "numcenter:% 4d numbasis:% 4d numgauss:% 4d "
00880 "numshell:% 4d\n", numcenter, numbasis, numgauss, numshell);
00881 #endif
00882 } while (strcmp(word[0],"****"));
00883 GET_LINE(buffer, data->file);
00884 i=sscanf(buffer,"%s%s",word[0], word[1]);
00885 } while (i == 2);
00886 data->wavef_size = numbasis;
00887 data->num_shells = numshell;
00888 data->num_basis_funcs = numgauss;
00889 data->num_basis_atoms = numcenter;
00890 } else {
00891 i=sscanf(buffer,"%d%s%s%d%s",&nfunc,word[0],word[1],
00892 &nprim,word[2]);
00893 if (i==5) {
00894 if ((strncmp(word[0],"basis",5) == 0) &&
00895 (strncmp(word[1],"functions",9) == 0) &&
00896 (strncmp(word[2],"primitive",9) == 0) ) {
00897 GET_LINE(buffer, data->file);
00898 sscanf(buffer,"%d%s%s%d%s",&nume_a,word[0],word[1],
00899 &nume_b,word[2]);
00900 }
00901 }
00902 }
00903 }
00904 } while (nume_a < 0);
00905
00906 vmdcon_printf(VMDCON_INFO,
00907 "gaussianplugin) %d shells with %d/%d basis functions, "
00908 "%d/%d primitive gaussians\n", data->num_shells,
00909 nfunc, data->wavef_size, nprim, data->num_basis_funcs);
00910 vmdcon_printf(VMDCON_INFO,
00911 "gaussianplugin) %d QM atoms with %d alpha electrons, "
00912 "%d beta electron\n", data->num_basis_atoms, nume_a, nume_b);
00913
00914 data->num_orbitals = nfunc;
00915 data->occ_orbitals_A = nume_a;
00916 data->occ_orbitals_B = nume_b;
00917 data->num_electrons = nume_a + nume_b;
00918 data->multiplicity = (nume_a+nume_b)/2-(nume_a>nume_b?nume_b:nume_a)+1;
00919
00920 return TRUE;
00921 }
00922
00923
00924
00925
00926
00927
00928
00929 static int get_runtitle(gaussiandata *data) {
00930
00931 char buffer[BUFSIZ];
00932 char *temp;
00933 size_t offs;
00934
00935
00936 do {
00937 GET_LINE(buffer, data->file);
00938 } while (strncmp(buffer," -", 2) );
00939
00940 GET_LINE(buffer, data->file);
00941 do {
00942 char s;
00943
00944 offs=strlen(buffer);
00945 temp=&buffer[offs];
00946
00947 while (*temp == '\n' || *temp == '\r' || *temp == '\0') --temp;
00948 s = *temp;
00949 fgets(temp,BUFSIZ-offs,data->file);
00950 *temp = s;
00951 } while ( strncmp(temp+1,"--",2) );
00952
00953 offs=strlen(buffer);
00954 temp=&buffer[offs-1];
00955 while (*temp == '-' || *temp == '\n' || *temp == '\r') --temp;
00956 ++temp; *temp='\0';
00957 strncpy(data->runtitle,buffer,sizeof(data->runtitle));
00958
00959 return TRUE;
00960 }
00961
00962
00963
00964 static int get_input_structure(gaussiandata *data) {
00965 char buffer[BUFSIZ];
00966 char word[4][MOLFILE_BUFSIZ];
00967 int i, numatoms;
00968
00969 buffer[0] = '\0';
00970 for (i=0; i<4; i++) word[i][0] = '\0';
00971
00972 GET_LINE(buffer, data->file);
00973 sscanf(buffer,"%s%s",word[0],word[1]);
00974
00975 if ( ( (strcmp(word[0],"Symbolic") == 0) &&
00976 (strcmp(word[1],"Z-matrix:") == 0) ) ||
00977 ( (strcmp(word[0],"Redundant") == 0) &&
00978 (strcmp(word[1],"internal") == 0) ) ||
00979 ( (strcmp(word[0],"Z-Matrix") == 0) &&
00980 (strcmp(word[1],"taken") == 0) ) ) {
00981
00982
00983 if ( ( (strcmp(word[0],"Redundant") == 0) &&
00984 (strcmp(word[1],"internal") == 0) ) ||
00985 ( (strcmp(word[0],"Z-Matrix") == 0) &&
00986 (strcmp(word[1],"taken") == 0) ) ) {
00987 GET_LINE(buffer, data->file);
00988 }
00989
00990
00991 GET_LINE(buffer, data->file);
00992 sscanf(buffer,"%*s%*s%d%*s%*s%d", &(data->totalcharge),
00993 &(data->multiplicity));
00994
00995
00996 numatoms=0;
00997 data->initatoms=NULL;
00998 i=1;
00999 do {
01000 char *ptr;
01001 qm_atom_t *atm;
01002
01003 i=1;
01004 GET_LINE(buffer, data->file);
01005 ptr = strtok(buffer," ,\t\n");
01006
01007 if (ptr == NULL) break;
01008 if (strcmp(ptr, "Variables:") == 0) break;
01009 if (strcmp(ptr, "Recover") == 0) break;
01010 if (strcmp(ptr, "The") == 0) break;
01011 if (strcmp(ptr, "Leave") == 0) break;
01012
01013
01014 data->initatoms=realloc(data->initatoms,(numatoms+1)*sizeof(qm_atom_t));
01015 atm = data->initatoms + numatoms;
01016 strncpy(atm->type, ptr, sizeof(atm->type));
01017 atm->atomicnum=get_pte_idx(ptr);
01018 ++numatoms;
01019 } while ( i >= 0 );
01020
01021 } else {
01022 vmdcon_printf(VMDCON_ERROR,
01023 "gaussianplugin) ERROR, cannot parse input coordinates.\n");
01024 return FALSE;
01025 }
01026
01027
01028 data->numatoms = numatoms;
01029 vmdcon_printf(VMDCON_INFO,
01030 "gaussianplugin) Atoms: %d Charge: %d Multiplicity: %d\n",
01031 numatoms, data->totalcharge, data->multiplicity);
01032 return TRUE;
01033 }
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046 static int get_contrl(gaussiandata *data) {
01047
01048 char buffer[BUFSIZ];
01049 const char *vmdbasis;
01050 char *temp;
01051 size_t offs;
01052
01053 buffer[0] = '\0';
01054 rewind(data->file);
01055
01056
01057 do {
01058 GET_LINE(buffer, data->file);
01059 } while( strncmp(buffer," #", 2) );
01060
01061
01062
01063
01064
01065
01066 do {
01067 char s;
01068
01069 offs=strlen(buffer);
01070 temp=&buffer[offs];
01071
01072 while (*temp == '\n' || *temp == '\r' || *temp == '\0') --temp;
01073 s = *temp;
01074 fgets(temp,BUFSIZ-offs,data->file);
01075 *temp = s;
01076 } while ( strncmp(temp+1,"--",2) );
01077
01078 offs=strlen(buffer);
01079 temp=&buffer[offs-1];
01080 while (*temp == '-' || *temp == '\n' || *temp == '\r') --temp;
01081
01082
01083
01084
01085
01086
01087 ++temp; *temp=' '; ++temp; *temp='\0';
01088
01089
01090 temp=&buffer[0];
01091 while (*temp++) *temp = toupper(*temp);
01092
01093
01094
01095
01096 if (strstr(buffer," IOP(6/7=3) ")) {
01097 data->have_wavefunction=TRUE;
01098 } else {
01099 data->have_wavefunction=FALSE;
01100 }
01101
01102
01103 if (strstr(buffer," GFINPUT ")) {
01104 data->have_basis=TRUE;
01105 } else {
01106 data->have_basis=FALSE;
01107 }
01108
01109
01110 if (strstr(buffer," 6D ")) {
01111 data->have_cart_basis |= 1;
01112 }
01113
01114 if (strstr(buffer," 10F ")) {
01115 data->have_cart_basis |= 2;
01116 }
01117
01118 if (strstr(buffer," 5D ")) {
01119 data->have_cart_basis &= ~1;
01120 }
01121
01122 if (strstr(buffer," 7F ")) {
01123 data->have_cart_basis &= ~2;
01124 }
01125
01126
01127 if ((strstr(buffer," ROHF/")) ||
01128 (strstr(buffer," ROHF ")) ||
01129 (strstr(buffer," ROMP"))) {
01130 data->scftyp = SCFTYP_ROHF;
01131 } else if (data->multiplicity != 1) {
01132 data->scftyp = SCFTYP_UHF;
01133 } else {
01134 data->scftyp = SCFTYP_RHF;
01135 }
01136
01137
01138 if ((strstr(buffer," AM1/")) ||
01139 (strstr(buffer," AM1 ")) ||
01140 (strstr(buffer," PM3/")) ||
01141 (strstr(buffer," PM3 ")) ||
01142 (strstr(buffer," MNDO/")) ||
01143 (strstr(buffer," MNDO "))) {
01144
01145 vmdbasis = getenv("VMDDEFBASISSET");
01146 if (vmdbasis == NULL)
01147 vmdbasis = "VSTO-3G";
01148
01149 if(strlen(data->gbasis) == 0)
01150 strncpy(data->gbasis, vmdbasis, sizeof(data->gbasis));
01151
01152 if(strlen(data->basis_string) == 0) {
01153 strncpy(data->basis_string, "Internal ", sizeof(data->basis_string));
01154 strncat(data->basis_string, vmdbasis, sizeof(data->basis_string) - 10);
01155
01156 if(data->have_cart_basis & 1)
01157 strcat(data->basis_string," 6D");
01158 else
01159 strcat(data->basis_string," 5D");
01160 if(data->have_cart_basis & 2)
01161 strcat(data->basis_string," 10F");
01162 else
01163 strcat(data->basis_string," 7F");
01164 }
01165 }
01166
01167
01168 data->runtyp = RUNTYP_ENERGY;
01169 if ((strstr(buffer," FOPT ")) ||
01170 (strstr(buffer," FOPT=")) ||
01171 (strstr(buffer," FOPT(")) ||
01172 (strstr(buffer," OPT=")) ||
01173 (strstr(buffer," OPT(")) ||
01174 (strstr(buffer," OPT "))) {
01175 data->runtyp = RUNTYP_OPTIMIZE;
01176 }
01177 if (strstr(buffer," FREQ ")) {
01178 data->runtyp = RUNTYP_HESSIAN;
01179 }
01180 if (strstr(buffer," SCAN ")) {
01181 data->runtyp = RUNTYP_SURFACE;
01182 }
01183
01184
01185
01186 vmdbasis = getenv("VMDDEFBASISSET");
01187 if(strlen(data->gbasis) == 0) {
01188 if (vmdbasis == NULL) {
01189 strncpy(data->gbasis, "(unknown)", sizeof(data->gbasis));
01190 strncpy(data->basis_string, "(unknown)", sizeof(data->basis_string));
01191 } else {
01192 strncpy(data->gbasis, vmdbasis, sizeof(data->gbasis));
01193 strncpy(data->basis_string, "Internal ", sizeof(data->basis_string));
01194 strncat(data->basis_string, vmdbasis, sizeof(data->basis_string) - 10);
01195
01196 if(data->have_cart_basis & 1)
01197 strcat(data->basis_string," 6D");
01198 else
01199 strcat(data->basis_string," 5D");
01200 if(data->have_cart_basis & 2)
01201 strcat(data->basis_string," 10F");
01202 else
01203 strcat(data->basis_string," 7F");
01204 }
01205 }
01206
01207
01208
01209
01210 vmdcon_printf(VMDCON_INFO, "gaussianplugin) Run-type: %s, SCF-type: %s\n",
01211 runtypes[data->runtyp], scftypes[data->scftyp]);
01212 vmdcon_printf(VMDCON_INFO,
01213 "gaussianplugin) using %s basis set.\n", data->basis_string);
01214
01215 return TRUE;
01216 }
01217
01218 static int read_first_frame(gaussiandata *data) {
01219 data->qm_timestep = NULL;
01220
01221
01222
01223
01224
01225 vmdcon_printf(VMDCON_INFO,
01226 "gaussianplugin) preparing for %d atomic basis functions "
01227 "per wavefunction.\n", data->wavef_size);
01228 SAFE_CALLOC(data->angular_momentum,int,3*data->wavef_size);
01229
01230 if (!get_traj_frame(data)) {
01231 return FALSE;
01232 }
01233
01234 data->num_frames = 1;
01235 return TRUE;
01236 }
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249 static int get_final_info(gaussiandata *data) {
01250 long filepos;
01251 filepos = ftell(data->file);
01252
01253 if (data->runtyp == RUNTYP_OPTIMIZE ||
01254 data->runtyp == RUNTYP_SADPOINT ||
01255 data->runtyp == RUNTYP_SURFACE) {
01256
01257
01258
01259 if (!find_traj_end(data)) return FALSE;
01260 }
01261
01262 if (data->runtyp == RUNTYP_HESSIAN || data->runtyp == RUNTYP_SURFACE) {
01263 vmdcon_printf(VMDCON_WARN, "gaussianplugin) this run type is not fully supported\n");
01264 }
01265
01266 fseek(data->file, filepos, SEEK_SET);
01267 return TRUE;
01268 }
01269
01270
01271 static int get_coordinates(FILE *file, qm_atom_t *atoms, int numatoms) {
01272 char buffer[BUFSIZ];
01273 int atomicnum;
01274 float x,y,z;
01275 int i,n;
01276
01277
01278
01279
01280
01281
01282
01283
01284 do {
01285 GET_LINE(buffer, file);
01286 } while ((strstr(buffer,"Input orientation:") == NULL) &&
01287 (strstr(buffer,"Z-Matrix orientation:") == NULL));
01288 GET_LINE(buffer, file);
01289 GET_LINE(buffer, file);
01290 GET_LINE(buffer, file);
01291 GET_LINE(buffer, file);
01292
01293 for (i=0; i < numatoms; ++i) {
01294 GET_LINE(buffer, file);
01295 n = sscanf(buffer,"%*d%d%*d%f%f%f",&atomicnum,&x,&y,&z);
01296 if (n!=4) return FALSE;
01297 atoms[i].x=x;
01298 atoms[i].y=y;
01299 atoms[i].z=z;
01300 }
01301 return TRUE;
01302 }
01303
01306 static int get_int_coordinates(FILE *file, qm_atom_t *atoms, int numatoms) {
01307 char buffer[BUFSIZ];
01308 int atomicnum;
01309 float x,y,z;
01310 int i,n;
01311
01312
01313
01314
01315
01316
01317
01318
01319 do {
01320 GET_LINE(buffer, file);
01321 } while ((strstr(buffer,"Standard orientation:") == NULL) &&
01322 (strstr(buffer,"Z-Matrix orientation:") == NULL));
01323 GET_LINE(buffer, file);
01324 GET_LINE(buffer, file);
01325 GET_LINE(buffer, file);
01326 GET_LINE(buffer, file);
01327
01328 for (i=0; i < numatoms; ++i) {
01329 GET_LINE(buffer, file);
01330 n = sscanf(buffer,"%*d%d%*d%f%f%f",&atomicnum,&x,&y,&z);
01331 if (n!=4) return FALSE;
01332 atoms[i].x=x;
01333 atoms[i].y=y;
01334 atoms[i].z=z;
01335 }
01336 return TRUE;
01337 }
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373 int get_basis(gaussiandata *data) {
01374
01375 char buffer[BUFSIZ];
01376 char word[3][MOLFILE_BUFSIZ];
01377 int i;
01378
01379
01380
01381 if (!data->have_basis) return FALSE;
01382
01383
01384 i=0;
01385 if(data->version < 20030000) {
01386 do {
01387 GET_LINE(buffer, data->file);
01388 sscanf(buffer,"%s%s",word[0],word[1]);
01389 ++i;
01390 } while( (strcmp(word[0],"Basis") ||
01391 strcmp(word[1],"set")) && (i<1000) );
01392 } else {
01393 do {
01394 GET_LINE(buffer, data->file);
01395 sscanf(buffer,"%s%s",word[0],word[1]);
01396 ++i;
01397 } while( (strcmp(word[0],"AO") ||
01398 strcmp(word[1],"basis")) && (i<1000) );
01399 }
01400 if (i>=1000) {
01401
01402 data->num_shells_per_atom=NULL;
01403 data->have_basis=FALSE;
01404 return MOLFILE_ERROR;
01405 }
01406
01407
01408
01409
01410 SAFE_CALLOC(data->basis_set,basis_atom_t,data->num_basis_atoms);
01411
01412 for (i=0; i< data->num_basis_atoms; ++i) {
01413 int numshells, numprim;
01414 int numread, ishell;
01415 float scalef;
01416 shell_t *shell;
01417
01418
01419
01420
01421
01422 GET_LINE(buffer, data->file);
01423 if (atoi(buffer) != (i+1) ) {
01424 vmdcon_printf(VMDCON_WARN,
01425 "gaussianplugin) basis set atom counter mismatch: %d vs %s\n",
01426 atoi(buffer), i+1);
01427 SAFE_FREE(data->basis_set);
01428 data->num_shells_per_atom=NULL;
01429 return MOLFILE_ERROR;
01430 }
01431 strncpy(data->basis_set[i].name,data->gbasis,10);
01432
01433 numshells=0;
01434 shell=NULL;
01435 GET_LINE(buffer, data->file);
01436 while( strncmp(buffer," ****",5) ) {
01437 int n;
01438
01439 numread=sscanf(buffer,"%s%d%f",word[0],&numprim,&scalef);
01440 if (numread == 3) {
01441 #if GAUSSIAN_DEBUG && GAUSSIAN_BASIS_DEBUG
01442 vmdcon_printf(VMDCON_INFO, "gaussianplugin) atom: %d, element: %s, shell: %d "
01443 "%s-type shell, %d primitives, scalefactor %f\n", i,
01444 get_pte_label(data->initatoms[i].atomicnum), numshells+1,
01445 word[0], numprim, scalef);
01446 #endif
01447 ;
01448 } else {
01449 vmdcon_printf(VMDCON_WARN,
01450 "gaussianplugin) basis set parse error: %s",buffer);
01451 free(data->basis_set);
01452 data->have_basis=FALSE;
01453 data->basis_set=NULL;
01454 data->num_shells_per_atom=NULL;
01455 return FALSE;
01456 }
01457 ishell=numshells;
01458 ++numshells;
01459 shell=realloc(shell,numshells*sizeof(shell_t));
01460 shell[ishell].numprims=numprim;
01461 shell[ishell].symmetry=shellsymm_int(word[0]);
01462 shell[ishell].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01463 if (shell[ishell].symmetry == SP_S_SHELL) {
01464 ++numshells;
01465 shell=realloc(shell,numshells*sizeof(shell_t));
01466 shell[ishell+1].numprims=numprim;
01467 shell[ishell+1].symmetry=SP_P_SHELL;
01468 shell[ishell+1].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01469 }
01470
01471 for (n=0; n<numprim; ++n) {
01472 GET_LINE(buffer, data->file);
01473 sscanf(buffer,"%s%s%s", word[0],word[1],word[2]);
01474 fix_fortran_exp(word[0]);
01475 shell[ishell].prim[n].exponent=atof(word[0])*scalef*scalef;
01476 fix_fortran_exp(word[1]);
01477 shell[ishell].prim[n].contraction_coeff=atof(word[1]);
01478 if (shell[ishell].symmetry == SP_S_SHELL) {
01479 shell[ishell+1].prim[n].exponent=shell[ishell].prim[n].exponent;
01480 fix_fortran_exp(word[2]);
01481 shell[ishell+1].prim[n].contraction_coeff=atof(word[2]);
01482 }
01483 }
01484 GET_LINE(buffer, data->file);
01485 }
01486
01487
01488 data->basis_set[i].numshells = numshells;
01489 data->basis_set[i].shell = shell;
01490 }
01491
01492 vmdcon_printf(VMDCON_INFO,
01493 "gaussianplugin) Parsed %d uncontracted basis functions with %d shells.\n",
01494 data->num_basis_funcs, data->num_shells);
01495
01496 #if GAUSSIAN_DEBUG
01497 for (i=0; i<data->num_basis_atoms; i++) {
01498 int j,k,primcount,shellcount;
01499 primcount=0;
01500 shellcount=0;
01501
01502 printf("%-8s (%10s)\n\n", data->initatoms[i].type, data->basis_set[i].name);
01503 printf("\n");
01504
01505 for (j=0; j<data->basis_set[i].numshells; j++) {
01506
01507 for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
01508 printf("%6d %d %7d %22f%22f\n", j,
01509 data->basis_set[i].shell[j].symmetry,
01510 primcount+1,
01511 data->basis_set[i].shell[j].prim[k].exponent,
01512 data->basis_set[i].shell[j].prim[k].contraction_coeff);
01513 primcount++;
01514 }
01515
01516 printf("\n");
01517 shellcount++;
01518 }
01519 }
01520 #endif
01521
01522
01523 return fill_basis_arrays(data);
01524 }
01525
01526
01527
01528
01529
01530
01531
01532
01533 int get_internal_basis(gaussiandata *data) {
01534
01535 char *vmddir = NULL;
01536 FILE *fp;
01537 char buffer[BUFSIZ];
01538 char word[3][MOLFILE_BUFSIZ];
01539 char filepath[256];
01540 int i,n;
01541
01542
01543 if (data->have_basis) return TRUE;
01544
01545
01546 sprintf(filepath,"%s.gbs",data->gbasis);
01547 fp=fopen(filepath,"rb");
01548 if (fp == NULL) {
01549 vmddir=getenv("VMDDIR");
01550 if (vmddir == NULL) {
01551 vmddir="/usr/local/lib/vmd";
01552 }
01553 sprintf(filepath,"%s/basis/%s.gbs",vmddir,data->gbasis);
01554 fp=fopen(filepath,"rb");
01555 }
01556
01557 if (fp == NULL) {
01558 vmdcon_printf(VMDCON_ERROR, "gaussianplugin) failed to read basis set "
01559 "from data base file %s\n", filepath);
01560 data->num_shells_per_atom=NULL;
01561 data->have_basis=FALSE;
01562 return FALSE;
01563 } else {
01564 vmdcon_printf(VMDCON_INFO, "gaussianplugin) reading basis set "
01565 "from data base file %s\n", filepath);
01566 }
01567
01568
01569
01570
01571 SAFE_CALLOC(data->basis_set,basis_atom_t,data->num_basis_atoms);
01572
01573 for (i=0; i < data->num_basis_atoms; ++i) {
01574 int numshells, numprim;
01575 int numread, ishell;
01576 float scalef;
01577 shell_t *shell;
01578
01579
01580 rewind(fp);
01581 do {
01582 fgets(buffer, sizeof(buffer), fp);
01583 sscanf(buffer,"%s%s",word[0],word[1]);
01584 } while(strcmp(word[0],"****"));
01585
01586
01587 do {
01588 fgets(buffer, sizeof(buffer), fp);
01589 if (feof(fp)) {
01590 free(data->basis_set);
01591 data->basis_set=NULL;
01592 data->num_shells_per_atom=NULL;
01593 data->have_basis=FALSE;
01594 vmdcon_printf(VMDCON_ERROR, "gaussianplugin) EOF in data base "
01595 "file %s while looking for element %s.\n", filepath,
01596 get_pte_label(data->initatoms[i].atomicnum), buffer);
01597 fclose(fp);
01598 return FALSE;
01599 }
01600 n=sscanf(buffer,"%s%s",word[0],word[1]);
01601 } while ( (n != 2) || strcmp(word[1],"0") ||
01602 (strcmp(word[0],get_pte_label(data->initatoms[i].atomicnum))) );
01603
01604 strncpy(data->basis_set[i].name,data->gbasis,sizeof(data->gbasis));
01605
01606 numshells=0;
01607 shell=NULL;
01608
01609
01610 do {
01611
01612 fgets(buffer, sizeof(buffer), fp);
01613 if (strstr(buffer,"****")) break;
01614 if (ferror(fp)) {
01615 vmdcon_printf(VMDCON_ERROR, "gaussianplugin) read error in data "
01616 "base file %s while reading basis of element %s.\n",
01617 filepath, get_pte_label(data->initatoms[i].atomicnum));
01618 free(data->basis_set);
01619 data->basis_set=NULL;
01620 return FALSE;
01621 }
01622
01623 numread=sscanf(buffer,"%s%d%f",word[0],&numprim,&scalef);
01624 if (numread == 3) {
01625 #if GAUSSIAN_DEBUG && GAUSSIAN_BASIS_DEBUG
01626 vmdcon_printf(VMDCON_INFO, "gaussianplugin) atom: %d, element: %s, shell: %d "
01627 "%s-type shell, %d primitives, scalefactor %f\n", i,
01628 get_pte_label(data->initatoms[i].atomicnum), numshells+1,
01629 word[0], numprim, scalef);
01630 #endif
01631 ;
01632 } else {
01633 vmdcon_printf(VMDCON_WARN,
01634 "gaussianplugin) basis set parse error: %s",buffer);
01635 free(data->basis_set);
01636 data->basis_set=NULL;
01637 return FALSE;
01638 }
01639
01640 ishell=numshells;
01641 ++numshells;
01642 shell=realloc(shell,numshells*sizeof(shell_t));
01643 shell[ishell].numprims=numprim;
01644 shell[ishell].symmetry=shellsymm_int(word[0]);
01645 shell[ishell].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01646 if (shell[ishell].symmetry == SP_S_SHELL) {
01647 ++numshells;
01648 shell=realloc(shell,numshells*sizeof(shell_t));
01649 shell[ishell+1].numprims=numprim;
01650 shell[ishell+1].symmetry=SP_P_SHELL;
01651 shell[ishell+1].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01652 }
01653
01654 for (n=0; n<numprim; ++n) {
01655 fgets(buffer, sizeof(buffer), fp);
01656 if (ferror(fp)) {
01657 vmdcon_printf(VMDCON_ERROR, "gaussianplugin) read error in data "
01658 "base file %s while reading basis of element %s.\n",
01659 filepath, get_pte_label(data->initatoms[i].atomicnum));
01660 free(data->basis_set);
01661 data->basis_set=NULL;
01662 return FALSE;
01663 }
01664 sscanf(buffer,"%s%s%s", word[0],word[1],word[2]);
01665 fix_fortran_exp(word[0]);
01666 shell[ishell].prim[n].exponent=atof(word[0])*scalef*scalef;
01667 fix_fortran_exp(word[1]);
01668 shell[ishell].prim[n].contraction_coeff=atof(word[1]);
01669 if (shell[ishell].symmetry == SP_S_SHELL) {
01670 shell[ishell+1].prim[n].exponent=shell[ishell].prim[n].exponent;
01671 fix_fortran_exp(word[2]);
01672 shell[ishell+1].prim[n].contraction_coeff=atof(word[2]);
01673 }
01674 }
01675 } while(1);
01676
01677
01678
01679 data->basis_set[i].numshells = numshells;
01680 data->basis_set[i].shell = shell;
01681
01682
01683 data->num_shells += numshells;
01684 }
01685
01686 vmdcon_printf(VMDCON_INFO,
01687 "gaussianplugin) Parsed %d uncontracted basis functions. \n",
01688 data->num_basis_funcs);
01689
01690
01691 data->have_basis = TRUE;
01692 return fill_basis_arrays(data);
01693 }
01694
01695
01696
01697
01698
01699
01700
01701 static int shellsymm_int(char *symm) {
01702 int shell_symmetry;
01703
01704 switch (toupper(symm[0])) {
01705 case 'S':
01706 if (symm[1] == '\0') {
01707 shell_symmetry = S_SHELL;
01708 } else if (toupper(symm[1]) == 'P') {
01709 if (symm[2] == '\0') {
01710 shell_symmetry = SP_S_SHELL;
01711 } else if (toupper(symm[1]) == 'D') {
01712 shell_symmetry = SPD_S_SHELL;
01713 } else {
01714 shell_symmetry = UNK_SHELL;
01715 }
01716 } else {
01717 shell_symmetry = UNK_SHELL;
01718 }
01719 break;
01720 case 'L':
01721 shell_symmetry = SP_S_SHELL;
01722 break;
01723 case 'M':
01724 shell_symmetry = SP_P_SHELL;
01725 break;
01726 case 'P':
01727 shell_symmetry = P_SHELL;
01728 break;
01729 case 'D':
01730 shell_symmetry = D_SHELL;
01731 break;
01732 case 'F':
01733 shell_symmetry = F_SHELL;
01734 break;
01735 case 'G':
01736 shell_symmetry = G_SHELL;
01737 break;
01738 default:
01739 shell_symmetry = UNK_SHELL;
01740 break;
01741 }
01742
01743 return shell_symmetry;
01744 }
01745
01746
01748 static int fill_basis_arrays(gaussiandata *data) {
01749 int i, j, k;
01750 int shellcount, primcount;
01751 float *basis;
01752 int *num_shells_per_atom;
01753 int *num_prim_per_shell;
01754 int *shell_symmetry;
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764 SAFE_CALLOC(basis,float,2*data->num_basis_funcs);
01765 SAFE_CALLOC(num_shells_per_atom,int,data->num_basis_atoms);
01766 SAFE_CALLOC(shell_symmetry,int,data->num_shells);
01767 SAFE_CALLOC(num_prim_per_shell,int,data->num_shells);
01768
01769
01770 data->basis = basis;
01771 data->shell_symmetry = shell_symmetry;
01772 data->num_shells_per_atom = num_shells_per_atom;
01773 data->num_prim_per_shell = num_prim_per_shell;
01774
01775 shellcount=primcount=0;
01776
01777 for(i=0; i<data->num_basis_atoms; i++) {
01778 num_shells_per_atom[i] = data->basis_set[i].numshells;
01779
01780 for (j=0; j<data->basis_set[i].numshells; j++) {
01781 shell_symmetry[shellcount] = data->basis_set[i].shell[j].symmetry;
01782 num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;
01783
01784 for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
01785 basis[2*primcount ] = data->basis_set[i].shell[j].prim[k].exponent;
01786 basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
01787 primcount++;
01788 }
01789 shellcount++;
01790 }
01791 }
01792 return TRUE;
01793 }
01794
01795
01796
01799 static int get_traj_frame(gaussiandata *data) {
01800 qm_timestep_t *cur_qm_ts;
01801 int i;
01802 long fpos;
01803 char buffer[BUFSIZ];
01804 char word[MOLFILE_BUFSIZ];
01805 buffer[0] = '\0';
01806 word[0] = '\0';
01807
01808 #if GAUSSIAN_DEBUG
01809 vmdcon_printf(VMDCON_INFO,
01810 "gaussianplugin) Timestep %d: =======================\n",
01811 data->num_frames_read);
01812 #endif
01813
01814 fpos=ftell(data->file);
01815 if (!get_coordinates(data->file, data->initatoms, data->numatoms)) {
01816 vmdcon_printf(VMDCON_WARN,
01817 "gaussianplugin) Couldn't find input orientation coordinates"
01818 " for timestep %d\n", data->num_frames_read);
01819 vmdcon_printf(VMDCON_WARN,
01820 "gaussianplugin) Trying internal coordinates instead.\n");
01821
01822 fseek(data->file, fpos, SEEK_SET);
01823 if(!get_int_coordinates(data->file, data->initatoms, data->numatoms)) {
01824 vmdcon_printf(VMDCON_ERROR,
01825 "gaussianplugin) Couldn't find any coordinates.\n",
01826 data->num_frames_read);
01827 return FALSE;
01828 }
01829 }
01830
01831
01832 data->qm_timestep =
01833 (qm_timestep_t *)realloc(data->qm_timestep,
01834 (data->num_frames_read+1)*sizeof(qm_timestep_t));
01835
01836
01837 cur_qm_ts = data->qm_timestep+data->num_frames_read;
01838 memset(cur_qm_ts, 0, sizeof(qm_timestep_t));
01839
01840
01841 cur_qm_ts->numwave = 1;
01842 if (data->scftyp == SCFTYP_UHF)
01843 cur_qm_ts->numwave = 2;
01844
01845
01846 if (data->num_frames_read == 0) {
01847 get_basis(data);
01848 get_internal_basis(data);
01849 }
01850
01851 SAFE_CALLOC(cur_qm_ts->wave,qm_wavefunction_t,cur_qm_ts->numwave);
01852 for (i=0; i < cur_qm_ts->numwave; ++i) {
01853 cur_qm_ts->wave[i].idtag = i;
01854 SAFE_CALLOC(cur_qm_ts->wave[i].orb_indices,int,data->wavef_size);
01855 SAFE_CALLOC(cur_qm_ts->wave[i].orb_energies,float,data->wavef_size);
01856 SAFE_CALLOC(cur_qm_ts->wave[i].occupancies,float,data->wavef_size);
01857 SAFE_CALLOC(cur_qm_ts->wave[i].wave_coeffs,float,
01858 data->wavef_size*data->wavef_size);
01859 }
01860
01861
01862 if (get_wavefunction(data, cur_qm_ts) == FALSE) {
01863 vmdcon_printf(VMDCON_WARN, "gaussianplugin) No wavefunction present for timestep %d\n", data->num_frames_read);
01864
01865 for (i=0; i < cur_qm_ts->numwave; ++i) {
01866 free(cur_qm_ts->wave[i].wave_coeffs);
01867 free(cur_qm_ts->wave[i].orb_energies);
01868 free(cur_qm_ts->wave[i].occupancies);
01869 }
01870 free(cur_qm_ts->wave);
01871 cur_qm_ts->wave=NULL;
01872 cur_qm_ts->numwave=0;
01873 } else {
01874 vmdcon_printf(VMDCON_INFO, "gaussianplugin) Wavefunction found for timestep %d\n", data->num_frames_read);
01875 }
01876
01877 #if 0
01878 if (get_population(data, cur_qm_ts)) {
01879 vmdcon_printf(VMDCON_INFO, "gaussianplugin) Mulliken charges found\n");
01880 }
01881 #endif
01882
01883 data->num_frames_read++;
01884
01885 return TRUE;
01886 }
01887
01888
01889
01890
01891
01892
01893
01894
01895 static int find_traj_end(gaussiandata *data) {
01896 char buffer[BUFSIZ];
01897 long filepos;
01898 filepos = ftell(data->file);
01899
01900 while (1) {
01901 if (!fgets(buffer, sizeof(buffer), data->file)) break;
01902
01903 if (strstr(buffer, "Berny optimization.")) {
01904 data->num_frames++;
01905 } else if (strstr(buffer, "Optimization completed.")) {
01906 #if GAUSSIAN_DEBUG
01907 vmdcon_printf(VMDCON_INFO, "gaussianplugin) ==== End of trajectory. ====\n");
01908 #endif
01909 data->opt_status = STATUS_CONVERGED;
01910 return TRUE;
01911 } else if (strstr(buffer, "Optimization stopped.")) {
01912 #if GAUSSIAN_DEBUG
01913 vmdcon_printf(VMDCON_INFO, "gaussianplugin) ==== End of trajectory. ====\n");
01914 #endif
01915 data->opt_status = STATUS_TOO_MANY_STEPS;
01916 return TRUE;
01917 } else if (strstr(buffer, "Convergence failure -- run terminated.")) {
01918 data->opt_status = STATUS_SCF_NOT_CONV;
01919 return FALSE;
01920 }
01921 }
01922
01923
01924
01925
01926 data->opt_status = STATUS_BROKEN_OFF;
01927
01928 fseek(data->file, filepos, SEEK_SET);
01929 return FALSE;
01930 }
01931
01932
01933
01934
01935
01936
01937
01938 static int get_wavefunction(gaussiandata *data, qm_timestep_t *ts)
01939 {
01940 long filepos;
01941 char buffer[BUFSIZ];
01942 char word[6][MOLFILE_BUFSIZ];
01943 float *orb_erg, *orb_occ;
01944 int orbital_counter;
01945 int i, j, num_values, num_orbs, *orb_idx;
01946 qm_wavefunction_t *wf;
01947
01948 i=j=num_values=num_orbs=orbital_counter=0;
01949
01950 buffer[0] = '\0';
01951 for (i=0; i<6; i++) word[i][0] = '\0';
01952
01953 wf = ts->wave;
01954 if (wf == NULL || ts->numwave == 0) {
01955 PRINTERR;
01956 return FALSE;
01957 }
01958
01959
01960 if (!data->have_wavefunction) return FALSE;
01961
01962
01963 wf->type = MOLFILE_WAVE_CANON;
01964 wf->spin = SPIN_ALPHA;
01965 wf->cartesian = 1;
01966 orb_erg = wf->orb_energies;
01967 orb_occ = wf->occupancies;
01968 orb_idx = wf->orb_indices;
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006 filepos = ftell(data->file);
02007
02008 do {
02009 GET_LINE(buffer, data->file);
02010
02011 if (strstr(buffer,"Input orientation") ) {
02012 fseek(data->file, filepos, SEEK_SET);
02013 return FALSE;
02014 }
02015 } while(!strstr(buffer,"Molecular Orbital Coefficients"));
02016
02017 while (orbital_counter < data->num_orbitals) {
02018
02019 GET_LINE(buffer, data->file);
02020 num_orbs = sscanf(buffer,"%s%s%s%s%s",word[0],word[1],word[2],word[3],word[4]);
02021
02022
02023
02024
02025
02026
02027 GET_LINE(buffer, data->file);
02028 sscanf(buffer,"%s%s%s%s%s", word[0],word[1],word[2],word[3],word[4]);
02029 for (i=0; i<num_orbs; i++) {
02030 j=strlen(word[i]);
02031 orb_occ[i] = (word[i][j-1] == 'O') ? 1.0f : 0.0f;
02032 }
02033
02034 GET_LINE(buffer, data->file);
02035 sscanf(buffer,"%*s%*s%s%s%s%s%s",word[0],word[1],word[2],word[3],word[4]);
02036 for (i=0; i<num_orbs; i++)
02037 orb_erg[i] = atof(word[i]);
02038
02039
02040 orb_erg += num_orbs;
02041 orb_occ += num_orbs;
02042
02043
02044 for (i=0; i<data->num_orbitals; i++) {
02045 int xexp=0, yexp=0, zexp=0;
02046
02047
02048
02049 GET_LINE(buffer, data->file);
02050 num_values = sscanf(buffer+12,"%4s%s%s%s%s%s",
02051 word[0], word[1], word[2],
02052 word[3], word[4], word[5]);
02053
02054
02055
02056 for (j=1; j<strlen(word[0]); j++) {
02057 switch (word[0][j]) {
02058 case 'X':
02059 xexp++;
02060 break;
02061 case 'Y':
02062 yexp++;
02063 break;
02064 case 'Z':
02065 zexp++;
02066 break;
02067
02068
02069 case '+':
02070 case '-':
02071 case '0':
02072 case '1':
02073 case '2':
02074 if (wf->cartesian) {
02075 wf->cartesian=0;
02076 vmdcon_printf(VMDCON_ERROR, "gaussianplugin) pure basis function "
02077 "detected: '%s'. those are not supported yet. bailing out...\n", word[0]);
02078 return FALSE;
02079 }
02080 break;
02081 default:
02082
02083 break;
02084 }
02085 }
02086 data->angular_momentum[3*i ] = xexp;
02087 data->angular_momentum[3*i+1] = yexp;
02088 data->angular_momentum[3*i+2] = zexp;
02089 #if GAUSSIAN_DEBUG && GAUSSIAN_BASIS_DEBUG && 0
02090 vmdcon_printf(VMDCON_INFO,"%s:%d orbitals %d/%d/%d/%d shell %d/%d/%s: %d %d %d\n",
02091 __FILE__, __LINE__, i, data->num_orbitals, num_orbs, orbital_counter,
02092 i, data->wavef_size, word[0], xexp, yexp, zexp);
02093 #endif
02094
02095
02096
02097
02098 for (j=0 ; j<num_orbs; j++) {
02099 wf->wave_coeffs[(orbital_counter+j)*data->wavef_size+i] = atof(&word[j+1][0]);
02100 }
02101 }
02102 orbital_counter += num_orbs;
02103 }
02104
02105
02106 wf->num_orbitals = orbital_counter;
02107 vmdcon_printf(VMDCON_INFO,
02108 "gaussianplugin) Number of orbitals scanned: %d \n",
02109 orbital_counter);
02110 return TRUE;
02111 }
02112
02113
02114
02115
02116
02117 static int get_population(gaussiandata *data, qm_timestep_t *ts) {
02118 #if 0
02119 int i;
02120 char buffer[BUFSIZ];
02121 data->have_mulliken = FALSE;
02122
02123
02124
02125 ts->mulliken_charges =
02126 (double *)calloc(data->num_basis_atoms, sizeof(double));
02127
02128 if (!ts->mulliken_charges) {
02129 PRINTERR;
02130 return FALSE;
02131 }
02132
02133 for (i=0; i<data->num_basis_atoms; i++) {
02134 int n;
02135 float mullpop, mullcharge, lowpop, lowcharge;
02136 GET_LINE(buffer, data->file);
02137 n = sscanf(buffer,"%*i%*s%f%f%f%f",
02138 &mullpop, &mullcharge, &lowpop, &lowcharge);
02139 if (n!=4) return FALSE;
02140 ts->mulliken_charges[i] = mullcharge;
02141 }
02142
02143 if (i!=data->numatoms) return FALSE;
02144
02145 data->have_mulliken = TRUE;
02146
02147 #if GAUSSIAN_DEBUG
02148 vmdcon_printf(VMDCON_INFO,
02149 "gaussianplugin) Number of orbitals scanned: %d \n",
02150 orbital_counter);
02151 #endif
02152
02153 #endif
02154 return TRUE;
02155 }
02156
02157
02158
02159
02160
02161
02162
02163
02164 VMDPLUGIN_API int VMDPLUGIN_init(void) {
02165 memset(&plugin, 0, sizeof(molfile_plugin_t));
02166 plugin.abiversion = vmdplugin_ABIVERSION;
02167 plugin.type = MOLFILE_PLUGIN_TYPE;
02168 plugin.name = "gaussian";
02169 plugin.prettyname = "Gaussian Logfile (g94,g98,g03)";
02170 plugin.author = "Axel Kohlmeyer, Markus Dittrich, Jan Saam";
02171 plugin.majorv = 0;
02172 plugin.minorv = 2;
02173 plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
02174 plugin.filename_extension = "log";
02175 plugin.open_file_read = open_gaussian_read;
02176 plugin.read_structure = read_gaussian_structure;
02177 plugin.close_file_read = close_gaussian_read;
02178
02179 plugin.read_qm_metadata = read_gaussian_metadata;
02180 plugin.read_qm_rundata = read_gaussian_rundata;
02181
02182 #if vmdplugin_ABIVERSION > 11
02183 plugin.read_timestep_metadata = read_timestep_metadata;
02184 plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
02185 plugin.read_timestep = read_timestep;
02186 #endif
02187 return VMDPLUGIN_SUCCESS;
02188 }
02189
02190 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
02191 (*cb)(v, (vmdplugin_t *)&plugin);
02192 return VMDPLUGIN_SUCCESS;
02193 }
02194
02195 VMDPLUGIN_API int VMDPLUGIN_fini(void) {
02196 return VMDPLUGIN_SUCCESS;
02197 }
02198
02199
02200 #ifdef TEST_PLUGIN
02201
02202 int main(int argc, char *argv[]) {
02203 int numatoms, i, j, optflags;
02204 molfile_atom_t *atoms;
02205 molfile_timestep_t timestep;
02206 molfile_metadata_t metadata;
02207 molfile_timestep_metadata_t ts_metadata;
02208 molfile_qm_timestep_metadata_t qm_ts_metadata;
02209 molfile_qm_metadata_t qm_metadata;
02210 molfile_qm_timestep_t qm_ts;
02211
02212 void *v;
02213
02214 while (--argc) {
02215 ++argv;
02216 v = open_gaussian_read(*argv, "log", &numatoms);
02217 if (!v) {
02218 fprintf(stderr, "open_gaussian_read failed for file %s\n", *argv);
02219 return 1;
02220 }
02221 fprintf(stderr, "open_gaussian_read succeeded for file %s\n", *argv);
02222 fprintf(stderr, "number of atoms: %d\n", numatoms);
02223 atoms = (molfile_atom_t *)malloc(sizeof(molfile_atom_t)*numatoms);
02224 read_gaussian_structure(v,&optflags, atoms);
02225
02226 i = 0;
02227 timestep.coords = (float *)malloc(3*sizeof(float)*numatoms);
02228
02229 while (1) {
02230 memset(&ts_metadata, 0, sizeof(molfile_timestep_metadata_t));
02231 read_timestep_metadata(v, &ts_metadata);
02232 memset(&qm_metadata, 0, sizeof(molfile_qm_metadata_t));
02233 read_qm_timestep_metadata(v, &qm_ts_metadata);
02234 qm_ts.scfenergies = (double *)malloc(qm_ts_metadata.num_scfiter*sizeof(double));
02235 qm_ts.wave = (molfile_qm_wavefunction_t *)malloc(qm_ts_metadata.num_wavef
02236 *sizeof(molfile_qm_wavefunction_t));
02237 memset(qm_ts.wave, 0, qm_ts_metadata.num_wavef*sizeof(molfile_qm_wavefunction_t));
02238 for (j=0; (j<10 && j<qm_ts_metadata.num_wavef); j++) {
02239 qm_ts.wave[j].wave_coeffs = (float *) malloc(qm_ts_metadata.num_orbitals_per_wavef[j]
02240 * qm_ts_metadata.wavef_size * sizeof(float));
02241 qm_ts.wave[j].orbital_energies = (float *) malloc(qm_ts_metadata.num_orbitals_per_wavef[j]*sizeof(float));
02242 }
02243 qm_ts.gradient = (float *) malloc(3*numatoms*sizeof(float));
02244 if (read_timestep(v, numatoms, ×tep, &qm_metadata, &qm_ts)) break;
02245
02246
02247
02248 free(qm_ts.gradient);
02249 for (j=0; (j<10 && j<qm_ts_metadata.num_wavef); j++) {
02250 free(qm_ts.wave[j].wave_coeffs);
02251 free(qm_ts.wave[j].orbital_energies);
02252 }
02253 free(qm_ts.wave);
02254 free(qm_ts.scfenergies);
02255 i++;
02256 }
02257 fprintf(stderr, "ended read_timestep on frame %d\n", i);
02258 close_gaussian_read(v);
02259 }
02260 return 0;
02261 }
02262 #endif