00001 /*************************************************************************** 00002 *cr 00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 00004 *cr University of Illinois 00005 *cr All Rights Reserved 00006 *cr 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * RCS INFORMATION: 00011 * 00012 * $RCSfile: QMData.C,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.108 $ $Date: 2020年02月24日 17:48:30 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * 00019 * The QMData class, which stores all QM calculation data that 00020 * are not dependent on the timestep (the latter are handled by 00021 * the QMTimestep class). 00022 * These data are, for instance, basis set and calculation input 00023 * parameters such as method and multiplicity. 00024 * Note that the wavefunctions are stored in QMTimestep but that 00025 * a signature of all wavefunctions that occur in the trajectory 00026 * is kept here in QMData. The signature is needed for setting up 00027 * the Orbital representation GUI. 00028 * 00029 * The basis set data are stored in hierarchical data structures 00030 * for convenient access and better readability of the non 00031 * performance critical code. However, the actual orbital computation 00032 * (performed by the Orbital class) needs simple contiguous 00033 * arrays which is why we keep those too. 00034 * 00035 ***************************************************************************/ 00036 00037 #include <math.h> 00038 #include <stdio.h> 00039 #include "Inform.h" 00040 #include "Timestep.h" 00041 #include "QMData.h" 00042 #include "QMTimestep.h" 00043 #include "Orbital.h" 00044 #include "Molecule.h" 00045 00046 //#define DEBUGGING 1 00047 00049 QMData::QMData(int natoms, int nbasis, int nshells, int nwave) : 00050 num_wave_f(nwave), 00051 num_basis(nbasis), 00052 num_atoms(natoms), 00053 num_shells(nshells) 00054 { 00055 num_wavef_signa = 0; 00056 wavef_signa = NULL; 00057 num_shells_per_atom = NULL; 00058 num_prim_per_shell = NULL; 00059 wave_offset = NULL; 00060 atom_types = NULL; 00061 atom_basis = NULL; 00062 basis_array = NULL; 00063 basis_set = NULL; 00064 shell_types = NULL; 00065 angular_momentum = NULL; 00066 norm_factors = NULL; 00067 carthessian = NULL; 00068 inthessian = NULL; 00069 wavenumbers = NULL; 00070 intensities = NULL; 00071 normalmodes = NULL; 00072 imagmodes = NULL; 00073 runtype = RUNTYPE_UNKNOWN; 00074 scftype = SCFTYPE_UNKNOWN; 00075 status = QMSTATUS_UNKNOWN; 00076 }; 00077 00078 00079 QMData::~QMData() { 00080 int i; 00081 for (i=0; i<num_wavef_signa; i++) { 00082 free(wavef_signa[i].orbids); 00083 free(wavef_signa[i].orbocc); 00084 } 00085 free(wavef_signa); 00086 delete [] basis_array; 00087 delete [] shell_types; 00088 delete [] num_shells_per_atom; 00089 delete [] num_prim_per_shell; 00090 delete [] atom_types; 00091 delete [] wave_offset; 00092 delete [] angular_momentum; 00093 delete [] carthessian; 00094 delete [] inthessian; 00095 delete [] wavenumbers; 00096 delete [] intensities; 00097 delete [] normalmodes; 00098 delete [] imagmodes; 00099 delete [] basis_string; 00100 delete [] atom_basis; 00101 if (norm_factors) { 00102 for (i=0; i<=highest_shell; i++) { 00103 if (norm_factors[i]) delete [] norm_factors[i]; 00104 } 00105 delete [] norm_factors; 00106 } 00107 if (basis_set) 00108 delete_basis_set(); 00109 } 00110 00111 // Free memory of the basis set 00112 void QMData::delete_basis_set() { 00113 int i, j; 00114 for (i=0; i<num_types; i++) { 00115 for (j=0; j<basis_set[i].numshells; j++) { 00116 delete [] basis_set[i].shell[j].prim; 00117 } 00118 delete [] basis_set[i].shell; 00119 } 00120 delete [] basis_set; 00121 00122 basis_set = NULL; 00123 } 00124 00125 00126 00130 void QMData::init_electrons(Molecule *mol, int totcharge) { 00131 00132 int i, nuclear_charge = 0; 00133 for (i=0; i<num_atoms; i++) { 00134 nuclear_charge += mol->atom(i)->atomicnumber; 00135 } 00136 00137 totalcharge = totcharge; 00138 num_electrons = nuclear_charge - totalcharge; 00139 //multiplicity = mult; 00140 00141 #if 0 00142 if (scftype == SCFTYPE_RHF) { 00143 if (mult!=1) { 00144 msgErr << "For RHF calculations the multiplicity has to be 1, but it is " 00145 << multiplicity << "!" 00146 << sendmsg; 00147 } 00148 if (num_electrons%2) { 00149 msgErr << "Unpaired electron(s) in RHF calculation!" 00150 << sendmsg; 00151 } 00152 num_orbitals_A = num_orbitals_B = num_electrons/2; 00153 } 00154 else if ( (scftype == SCFTYPE_ROHF) || 00155 (scftype == SCFTYPE_UHF) ) { 00156 num_orbitals_B = (num_electrons-multiplicity+1)/2; 00157 num_orbitals_A = num_electrons-num_orbitals_B; 00158 } 00159 #endif 00160 } 00161 00162 00163 00164 // ====================================== 00165 // Functions for basis set initialization 00166 // ====================================== 00167 00168 00169 // Populate basis set data and organize them into 00170 // hierarcical data structures. 00171 int QMData::init_basis(Molecule *mol, int num_basis_atoms, 00172 const char *bstring, 00173 const float *basis, 00174 const int *atomic_numbers, 00175 const int *nshells, 00176 const int *nprims, 00177 const int *shelltypes) { 00178 num_types = num_basis_atoms; 00179 00180 basis_string = new char[1+strlen(bstring)]; 00181 strcpy(basis_string, bstring); 00182 00183 if (!basis && (!strcmp(basis_string, "MNDO") || 00184 !strcmp(basis_string, "AM1") || 00185 !strcmp(basis_string, "PM3"))) { 00186 // Semiempirical methods are based on STOs. 00187 // The only parameter we need for orbital rendering 00188 // are the exponents zeta for S, P, D,... shells for 00189 // each atom. Since most QM packages don't print these 00190 // values we have to generate the basis set here using 00191 // hardcoded table values. 00192 00193 // generate_sto_basis(basis_string); 00194 00195 return 1; 00196 } 00197 00198 int i, j; 00199 00200 00201 // Copy the basis set arrays over. 00202 if (!basis || !num_basis) return 1; 00203 basis_array = new float[2*num_basis]; 00204 memcpy(basis_array, basis, 2*num_basis*sizeof(float)); 00205 00206 if (!nshells || !num_basis_atoms) return 0; 00207 num_shells_per_atom = new int[num_basis_atoms]; 00208 memcpy(num_shells_per_atom, nshells, num_basis_atoms*sizeof(int)); 00209 00210 if (!nprims || !num_shells) return 0; 00211 num_prim_per_shell = new int[num_shells]; 00212 memcpy(num_prim_per_shell, nprims, num_shells*sizeof(int)); 00213 00214 if (!shelltypes || !num_shells) return 0; 00215 shell_types = new int[num_shells]; 00216 highest_shell = 0; 00217 for (i=0; i<num_shells; i++) { 00218 // copy shell types ({0, 1, 2, ...} meaning {S, P, D, ...}) 00219 shell_types[i] = shelltypes[i]; 00220 00221 // Translate the combined shell types that have negative 00222 // codes into their corresponding basic types in order 00223 // to be able to determine the highest shell. 00224 // The highest shell is needed by init_angular_norm_factors(). 00225 int basictype = shell_types[i]; 00226 switch (basictype) { 00227 case SP_S_SHELL: basictype = S_SHELL; break; 00228 case SP_P_SHELL: basictype = P_SHELL; break; 00229 case SPD_S_SHELL: basictype = S_SHELL; break; 00230 case SPD_P_SHELL: basictype = P_SHELL; break; 00231 case SPD_D_SHELL: basictype = D_SHELL; break; 00232 } 00233 if (basictype>highest_shell) highest_shell = basictype; 00234 } 00235 #ifdef DEBUGGING 00236 printf("highest shell = %d\n", highest_shell); 00237 #endif 00238 00239 // Create table of angular normalization constants 00240 init_angular_norm_factors(); 00241 00242 // Organize basis set data hierarchically 00243 int boffset = 0; 00244 int shell_counter = 0; 00245 int numcartpershell[14] = {1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 93, 107}; 00246 basis_set = new basis_atom_t[num_basis_atoms]; 00247 00248 for (i=0; i<num_basis_atoms; i++) { 00249 basis_set[i].atomicnum = atomic_numbers[i]; 00250 basis_set[i].numshells = num_shells_per_atom[i]; 00251 basis_set[i].shell = new shell_t[basis_set[i].numshells]; 00252 00253 for (j=0; j<basis_set[i].numshells; j++) { 00254 // We keep the info about the origin of a shell from a 00255 // combined shell (e.g. SP shell having common exponents 00256 // for S and P) in an extra flag in the basis_set structure. 00257 // In shell_types we just store S, P, D ... (0, 1, 2, ...) 00258 // because this is the relevant info. 00259 switch (shell_types[shell_counter]) { 00260 case SP_S_SHELL: 00261 shell_types[shell_counter] = S_SHELL; 00262 basis_set[i].shell[j].combo = 1; 00263 break; 00264 case SP_P_SHELL: 00265 shell_types[shell_counter] = P_SHELL; 00266 basis_set[i].shell[j].combo = 1; 00267 break; 00268 case SPD_S_SHELL: 00269 shell_types[shell_counter] = S_SHELL; 00270 basis_set[i].shell[j].combo = 2; 00271 break; 00272 case SPD_P_SHELL: 00273 shell_types[shell_counter] = P_SHELL; 00274 basis_set[i].shell[j].combo = 2; 00275 break; 00276 case SPD_D_SHELL: 00277 shell_types[shell_counter] = D_SHELL; 00278 basis_set[i].shell[j].combo = 2; 00279 break; 00280 } 00281 00282 basis_set[i].shell[j].type = shell_types[shell_counter]; 00283 00284 int shelltype = shell_types[shell_counter]; 00285 basis_set[i].shell[j].num_cart_func = numcartpershell[shelltype]; 00286 basis_set[i].shell[j].basis = basis_array+2*boffset; 00287 basis_set[i].shell[j].norm_fac = norm_factors[shelltype]; 00288 basis_set[i].shell[j].numprims = num_prim_per_shell[shell_counter]; 00289 00290 basis_set[i].shell[j].prim = new prim_t[basis_set[i].shell[j].numprims]; 00291 #ifdef DEBUGGING 00292 //printf("atom %i shell %i %s\n", i, j, get_shell_type_str(&basis_set[i].shell[j])); 00293 #endif 00294 00295 int k; 00296 for (k=0; k<basis_set[i].shell[j].numprims; k++) { 00297 float expon = basis_array[2*(boffset+k) ]; 00298 float coeff = basis_array[2*(boffset+k)+1]; 00299 basis_set[i].shell[j].prim[k].expon = expon; 00300 basis_set[i].shell[j].prim[k].coeff = coeff; 00301 } 00302 00303 // Offsets to get to this shell in the basis array. 00304 boffset += basis_set[i].shell[j].numprims; 00305 00306 shell_counter++; 00307 } 00308 } 00309 00310 00311 00312 // Collapse basis set so that we have one basis set 00313 // per atom type. 00314 if (!create_unique_basis(mol, num_basis_atoms)) { 00315 return 0; 00316 } 00317 00318 // Multiply the contraction coefficients with 00319 // the shell dependent part of the normalization factor. 00320 normalize_basis(); 00321 00322 return 1; 00323 } 00324 00325 00326 // ================================================= 00327 // Helper functions for building the list of unique 00328 // basis set atoms 00329 // ================================================= 00330 00331 // Return 1 if the two given shell basis sets are identical, 00332 // otherwise return 0. 00333 static int compare_shells(const shell_t *s1, const shell_t *s2) { 00334 if (s1->type != s2->type) return 0; 00335 if (s1->numprims != s2->numprims) return 0; 00336 int i; 00337 for (i=0; i<s1->numprims; i++) { 00338 if (s1->prim[i].expon != s2->prim[i].expon) return 0; 00339 if (s1->prim[i].coeff != s2->prim[i].coeff) return 0; 00340 } 00341 return 1; 00342 } 00343 00344 // Return 1 if the two given atomic basis sets are identical, 00345 // otherwise return 0. 00346 static int compare_atomic_basis(const basis_atom_t *a1, const basis_atom_t *a2) { 00347 if (a2->atomicnum != a1->atomicnum) return 0; 00348 if (a1->numshells != a2->numshells) return 0; 00349 int i; 00350 for (i=0; i<a1->numshells; i++) { 00351 if (!compare_shells(&a1->shell[i], &a2->shell[i])) return 0; 00352 } 00353 return 1; 00354 } 00355 00356 static void copy_shell_basis(const shell_t *s1, shell_t *s2) { 00357 s2->numprims = s1->numprims; 00358 s2->type = s1->type; 00359 s2->combo = s1->combo; 00360 s2->norm_fac = s1->norm_fac; 00361 s2->num_cart_func = s1->num_cart_func; 00362 s2->prim = new prim_t[s2->numprims]; 00363 int i; 00364 for (i=0; i<s2->numprims; i++) { 00365 s2->prim[i].expon = s1->prim[i].expon; 00366 s2->prim[i].coeff = s1->prim[i].coeff; 00367 } 00368 } 00369 00370 static void copy_atomic_basis(const basis_atom_t *a1, basis_atom_t *a2) { 00371 a2->atomicnum = a1->atomicnum; 00372 a2->numshells = a1->numshells; 00373 a2->shell = new shell_t[a2->numshells]; 00374 int i; 00375 for (i=0; i<a2->numshells; i++) { 00376 copy_shell_basis(&a1->shell[i], &a2->shell[i]); 00377 } 00378 } 00379 00380 // Collapse basis set so that we have one basis set per 00381 // atom type rather that per atom. In most cases an atom 00382 // type is a chemical element. Create an array that maps 00383 // individual atoms to their corresponding atomic basis. 00384 int QMData::create_unique_basis(Molecule *mol, int num_basis_atoms) { 00385 basis_atom_t *unique_basis = new basis_atom_t[num_basis_atoms]; 00386 copy_atomic_basis(&basis_set[0], &unique_basis[0]); 00387 int num_unique_atoms = 1; 00388 int i, j, k; 00389 for (i=1; i<num_basis_atoms; i++) { 00390 int found = 0; 00391 for (j=0; j<num_unique_atoms; j++) { 00392 if (compare_atomic_basis(&basis_set[i], &unique_basis[j])) { 00393 found = 1; 00394 break; 00395 } 00396 } 00397 if (!found) { 00398 copy_atomic_basis(&basis_set[i], &unique_basis[j]); 00399 num_unique_atoms++; 00400 } 00401 } 00402 00403 msgInfo << "Number of unique atomic basis sets = " 00404 << num_unique_atoms <<"/"<< num_atoms << sendmsg; 00405 00406 00407 // Free memory of the basis set 00408 delete_basis_set(); 00409 delete [] basis_array; 00410 00411 num_types = num_unique_atoms; 00412 basis_set = unique_basis; 00413 00414 // Count the new number of basis functions 00415 num_basis = 0; 00416 for (i=0; i<num_types; i++) { 00417 for (j=0; j<basis_set[i].numshells; j++) { 00418 num_basis += basis_set[i].shell[j].numprims; 00419 } 00420 } 00421 00422 basis_array = new float[2*num_basis]; 00423 int *basis_offset = new int[num_types]; 00424 00425 int ishell = 0; 00426 int iprim = 0; 00427 for (i=0; i<num_types; i++) { 00428 basis_offset[i] = iprim; 00429 00430 for (j=0; j<basis_set[i].numshells; j++) { 00431 basis_set[i].shell[j].basis = basis_array+iprim; 00432 #ifdef DEBUGGING 00433 printf("atom type %i shell %i %s\n", i, j, get_shell_type_str(&basis_set[i].shell[j])); 00434 #endif 00435 for (k=0; k<basis_set[i].shell[j].numprims; k++) { 00436 basis_array[iprim ] = basis_set[i].shell[j].prim[k].expon; 00437 basis_array[iprim+1] = basis_set[i].shell[j].prim[k].coeff; 00438 #ifdef DEBUGGING 00439 printf("prim %i: % 9.2f % 9.6f \n", k, basis_array[iprim], basis_array[iprim+1]); 00440 #endif 00441 iprim += 2; 00442 } 00443 ishell++; 00444 } 00445 } 00446 00447 atom_types = new int[num_atoms]; 00448 00449 // Assign basis set type to each atom and 00450 // create array of offsets into basis_array. 00451 for (i=0; i<num_atoms; i++) { 00452 int found = 0; 00453 for (j=0; j<num_types; j++) { 00454 //printf("atomicnum %d--%d\n", basis_set[j].atomicnum, mol->atom(i)->atomicnumber); 00455 if (basis_set[j].atomicnum == mol->atom(i)->atomicnumber) { 00456 found = 1; 00457 break; 00458 } 00459 } 00460 if (!found) { 00461 msgErr << "Error reading QM data: Could not assign basis set type to atom " 00462 << i << "." << sendmsg; 00463 delete_basis_set(); 00464 delete [] basis_offset; 00465 return 0; 00466 } 00467 atom_types[i] = j; 00468 #ifdef DEBUGGING 00469 printf("atom_types[%d]=%d\n", i, j); 00470 #endif 00471 } 00472 00473 // Count the new number of shells 00474 num_shells = 0; 00475 for (i=0; i<num_atoms; i++) { 00476 num_shells += basis_set[atom_types[i]].numshells; 00477 } 00478 00479 // Reallocate symmetry expanded arrays 00480 delete [] shell_types; 00481 delete [] num_prim_per_shell; 00482 delete [] num_shells_per_atom; 00483 shell_types = new int[num_shells]; 00484 num_prim_per_shell = new int[num_shells]; 00485 num_shells_per_atom = new int[num_atoms]; 00486 atom_basis = new int[num_atoms]; 00487 wave_offset = new int[num_atoms]; 00488 int shell_counter = 0; 00489 int woffset = 0; 00490 00491 // Populate the arrays again. 00492 for (i=0; i<num_atoms; i++) { 00493 int type = atom_types[i]; 00494 00495 // Offsets into wavefunction array 00496 wave_offset[i] = woffset; 00497 00498 for (j=0; j<basis_set[type].numshells; j++) { 00499 shell_t *shell = &basis_set[type].shell[j]; 00500 00501 woffset += shell->num_cart_func; 00502 00503 shell_types[shell_counter] = shell->type; 00504 num_prim_per_shell[shell_counter] = shell->numprims; 00505 shell_counter++; 00506 } 00507 00508 num_shells_per_atom[i] = basis_set[type].numshells; 00509 00510 // Offsets into basis_array 00511 atom_basis[i] = basis_offset[type]; 00512 } 00513 00514 delete [] basis_offset; 00515 00516 return 1; 00517 } 00518 00519 00520 00521 // Multiply the contraction coefficients with 00522 // the shell dependent part of the normalization factor 00523 // N = (2a/pi)^(3/4) * sqrt[(8a)^L]. 00524 // Here L denotes the shell type with L={0,1,2,3,...} 00525 // for {S,P,D,F,...} shells. 00526 // 00527 // Note that the angular momentum dependent part of the 00528 // normalization factor is 00529 // n = sqrt[i!j!k!/((2i)!(2j)!(2k)!)] 00530 // These values are stored in a table by 00531 // init_angular_norm_factors(). 00532 // 00533 void QMData::normalize_basis() { 00534 int i, j, k; 00535 for (i=0; i<num_types; i++) { 00536 for (j=0; j<basis_set[i].numshells; j++) { 00537 shell_t *shell = &basis_set[i].shell[j]; 00538 int shelltype = shell->type; 00539 for (k=0; k<shell->numprims; k++) { 00540 float expon = shell->prim[k].expon; 00541 float norm = (float) (pow(2.0*expon/VMD_PI, 0.75)*sqrt(pow(8*expon, shelltype))); 00542 #ifdef DEBUGGING 00543 //printf("prim %i: % 9.2f % 9.6f norm=%f\n", k, expon, coeff, norm); 00544 #endif 00545 shell->basis[2*k+1] = norm*shell->prim[k].coeff; 00546 } 00547 } 00548 } 00549 } 00550 00551 // Computes the factorial of n 00552 static int fac(int n) { 00553 if (n==0) return 1; 00554 int i, x=1; 00555 for (i=1; i<=n; i++) x*=i; 00556 return x; 00557 } 00558 00559 // Computes the factorial of n 00560 // Caution: Recursive function! Don't use with large n. 00561 // (For the overlap integrals we need only very small n.) 00562 static int doublefac(int n) { 00563 if (n<=1) return 1; 00564 return n*doublefac(n-2); 00565 } 00566 00567 // Initialize table of angular momentum dependent normalization 00568 // factors containing different factors for each shell and its 00569 // cartesian functions. 00570 void QMData::init_angular_norm_factors() { 00571 int shell; 00572 norm_factors = new float*[highest_shell+1]; 00573 for (shell=0; shell<=highest_shell; shell++) { 00574 int i, j, k; 00575 int numcart = 0; 00576 for (i=0; i<=shell; i++) numcart += i+1; 00577 00578 norm_factors[shell] = new float[numcart]; 00579 int count = 0; 00580 for (k=0; k<=shell; k++) { 00581 for (j=0; j<=shell; j++) { 00582 for (i=0; i<=shell; i++) { 00583 if (i+j+k==shell) { 00584 #ifdef DEBUGGING 00585 printf("count=%i (%i%i%i) %f\n", count, i, j, k, sqrt(((float)(fac(i)*fac(j)*fac(k))) / (fac(2*i)*fac(2*j)*fac(2*k)))); 00586 #endif 00587 norm_factors[shell][count++] = (float) sqrt(((float)(fac(i)*fac(j)*fac(k))) / (fac(2*i)*fac(2*j)*fac(2*k))); 00588 } 00589 } 00590 } 00591 } 00592 } 00593 } 00594 00595 00596 00597 // ================= 00598 // Basis set acccess 00599 // ================= 00600 00601 00602 // Get basis set for an atom 00603 const basis_atom_t* QMData::get_basis(int atom) const { 00604 if (!basis_set || !num_types || atom<0 || atom>=num_atoms) 00605 return NULL; 00606 return &(basis_set[atom_types[atom]]); 00607 } 00608 00609 00610 // Get basis set for a shell 00611 const shell_t* QMData::get_basis(int atom, int shell) const { 00612 if (!basis_set || !num_types || atom<0 || atom>=num_atoms || 00613 shell<0 || shell>=basis_set[atom_types[atom]].numshells) 00614 return NULL; 00615 return &(basis_set[atom_types[atom]].shell[shell]); 00616 } 00617 00618 00619 int QMData::get_num_wavecoeff_per_atom(int atom) const { 00620 if (atom<0 || atom>num_atoms) { 00621 msgErr << "Atom "<<atom<<" does not exist!"<<sendmsg; 00622 return -1; 00623 } 00624 int i; 00625 int a = atom_types[atom]; 00626 int num_cart_func = 0; 00627 for (i=0; i<basis_set[a].numshells; i++) { 00628 num_cart_func += basis_set[a].shell[i].num_cart_func; 00629 } 00630 return num_cart_func; 00631 } 00632 00633 // Get the offset in the wavefunction array for a specified 00634 // shell in an atom. 00635 int QMData::get_wave_offset(int atom, int shell) const { 00636 if (atom<0 || atom>num_atoms) { 00637 msgErr << "Atom "<<atom<<" does not exist!"<<sendmsg; 00638 return -1; 00639 } 00640 if (shell<0 || shell>=basis_set[atom_types[atom]].numshells) { 00641 msgErr << "Shell "<<shell<<" in atom "<<atom 00642 << " does not exist!"<<sendmsg; 00643 return -1; 00644 } 00645 int i; 00646 int numcart = 0; 00647 for (i=0; i<shell; i++) { 00648 numcart += basis_set[atom_types[atom]].shell[i].num_cart_func; 00649 } 00650 return wave_offset[atom]+numcart; 00651 } 00652 00653 00655 const char* QMData::get_shell_type_str(const shell_t *shell) { 00656 const char* map[14] = {"S0円", "P0円", "D0円", "F0円", "G0円", "H0円", 00657 "I0円", "K0円", "L0円", "M0円", "N0円", "O0円", "Q0円", "R0円"}; 00658 00659 return map[shell->type]; 00660 } 00661 00662 00663 00664 int QMData::set_angular_momenta(const int *angmom) { 00665 if (!angmom || !num_wave_f) return 0; 00666 angular_momentum = new int[3*num_wave_f]; 00667 memcpy(angular_momentum, angmom, 3*num_wave_f*sizeof(int)); 00668 return 1; 00669 } 00670 00671 void QMData::set_angular_momentum(int atom, int shell, int mom, 00672 int *array) { 00673 if (!array || !angular_momentum) return; 00674 int offset = get_wave_offset(atom, shell); 00675 if (offset<0) return; 00676 memcpy(&angular_momentum[3*(offset+mom)], array, 3*sizeof(int)); 00677 } 00678 00679 00680 // For a certain atom and shell return the exponent of the 00681 // requested cartesian component of the angular momentum 00682 // (specified by comp=0,1,2 for x,y,z resp.). 00683 // Example: 00684 // For XYYYZZ the exponents of the angular momentum are 00685 // X (comp 0): 1 00686 // Y (comp 1): 3 00687 // Y (comp 2): 2 00688 int QMData::get_angular_momentum(int atom, int shell, int mom, int comp) { 00689 if (!angular_momentum) return -1; 00690 int offset = get_wave_offset(atom, shell); 00691 if (offset<0 || 00692 mom>=get_basis(atom, shell)->num_cart_func) return -1; 00693 //printf("atom=%d, shell=%d, mom=%d, comp=%d\n", atom, shell, mom, comp); 00694 return angular_momentum[3*(offset+mom)+comp]; 00695 } 00696 00697 00698 // Set the angular momentum from a string 00699 void QMData::set_angular_momentum_str(int atom, int shell, int mom, 00700 const char *tag) { 00701 unsigned int j; 00702 int offset = get_wave_offset(atom, shell); 00703 if (offset<0) return; 00704 00705 int xexp=0, yexp=0, zexp=0; 00706 00707 for (j=0; j<strlen(tag); j++) { 00708 switch (tag[j]) { 00709 case 'X': 00710 xexp++; 00711 break; 00712 case 'Y': 00713 yexp++; 00714 break; 00715 case 'Z': 00716 zexp++; 00717 break; 00718 } 00719 } 00720 angular_momentum[3*(offset+mom) ] = xexp; 00721 angular_momentum[3*(offset+mom)+1] = yexp; 00722 angular_momentum[3*(offset+mom)+2] = zexp; 00723 } 00724 00725 00726 // Returns a pointer to a string representing the angular 00727 // momentum of a certain cartesian basis function. 00728 // The strings for an F-shell would be for instance 00729 // XX YY ZZ XY XZ YZ. 00730 // The necessary memory is automatically allocated. 00731 // Caller is responsible delete the string! 00732 char* QMData::get_angular_momentum_str(int atom, int shell, int mom) const { 00733 int offset = get_wave_offset(atom, shell); 00734 if (offset<0) return NULL; 00735 00736 char *s = new char[2+basis_set[atom_types[atom]].shell[shell].type]; 00737 int i, j=0; 00738 for (i=0; i<angular_momentum[3*(offset+mom) ]; i++) s[j++]='X'; 00739 for (i=0; i<angular_momentum[3*(offset+mom)+1]; i++) s[j++]='Y'; 00740 for (i=0; i<angular_momentum[3*(offset+mom)+2]; i++) s[j++]='Z'; 00741 s[j] = '0円'; 00742 if (!strlen(s)) strcpy(s, "S"); 00743 00744 return s; 00745 } 00746 00747 00748 00749 // ======================== 00750 // Hessian and normal modes 00751 // ======================== 00752 00753 00754 void QMData::set_carthessian(int numcart, double *array) { 00755 if (!array || !numcart || numcart!=3*num_atoms) return; 00756 carthessian = new double[numcart*numcart]; 00757 memcpy(carthessian, array, numcart*numcart*sizeof(double)); 00758 } 00759 00760 void QMData::set_inthessian(int numint, double *array) { 00761 if (!array || !numint) return; 00762 nintcoords = numint; 00763 inthessian = new double[numint*numint]; 00764 memcpy(inthessian, array, numint*numint*sizeof(double)); 00765 } 00766 00767 void QMData::set_normalmodes(int numcart, float *array) { 00768 if (!array || !numcart || numcart!=3*num_atoms) return; 00769 normalmodes = new float[numcart*numcart]; 00770 memcpy(normalmodes, array, numcart*numcart*sizeof(float)); 00771 } 00772 00773 void QMData::set_wavenumbers(int numcart, float *array) { 00774 if (!array || !numcart || numcart!=3*num_atoms) return; 00775 wavenumbers = new float[numcart]; 00776 memcpy(wavenumbers, array, numcart*sizeof(float)); 00777 } 00778 00779 void QMData::set_intensities(int numcart, float *array) { 00780 if (!array || !numcart || numcart!=3*num_atoms) return; 00781 intensities = new float[numcart]; 00782 memcpy(intensities, array, numcart*sizeof(float)); 00783 } 00784 00785 void QMData::set_imagmodes(int numimag, int *array) { 00786 if (!array || !numimag) return; 00787 nimag = numimag; 00788 imagmodes = new int[nimag]; 00789 memcpy(imagmodes, array, nimag*sizeof(int)); 00790 } 00791 00792 00793 00794 // ===================== 00795 // Calculation meta data 00796 // ===================== 00797 00798 00799 // Translate the runtype constant into a string 00800 const char *QMData::get_runtype_string(void) const 00801 { 00802 switch (runtype) { 00803 case RUNTYPE_ENERGY: return "single point energy"; break; 00804 case RUNTYPE_OPTIMIZE: return "geometry optimization"; break; 00805 case RUNTYPE_SADPOINT: return "saddle point search"; break; 00806 case RUNTYPE_HESSIAN: return "Hessian/frequency"; break; 00807 case RUNTYPE_SURFACE: return "potential surface scan"; break; 00808 case RUNTYPE_GRADIENT: return "energy gradient"; break; 00809 case RUNTYPE_MEX: return "minimum energy crossing"; break; 00810 case RUNTYPE_DYNAMICS: return "molecular dynamics"; break; 00811 case RUNTYPE_PROPERTIES: return "properties"; break; 00812 default: return "(unknown)"; break; 00813 } 00814 } 00815 00816 00817 // Translate the scftype constant into a string 00818 const char *QMData::get_scftype_string(void) const 00819 { 00820 switch (scftype) { 00821 case SCFTYPE_NONE: return "NONE"; break; 00822 case SCFTYPE_RHF: return "RHF"; break; 00823 case SCFTYPE_UHF: return "UHF"; break; 00824 case SCFTYPE_ROHF: return "ROHF"; break; 00825 case SCFTYPE_GVB: return "GVB"; break; 00826 case SCFTYPE_MCSCF: return "MCSCF"; break; 00827 case SCFTYPE_FF: return "force field"; break; 00828 default: return "(unknown)"; break; 00829 } 00830 } 00831 00832 00833 // Get status of SCF and optimization convergence 00834 const char* QMData::get_status_string() { 00835 if (status==QMSTATUS_OPT_CONV) 00836 return "Optimization converged"; 00837 else if (status==QMSTATUS_OPT_NOT_CONV) 00838 return "Optimization not converged"; 00839 else if (status==QMSTATUS_SCF_NOT_CONV) 00840 return "SCF not converged"; 00841 else if (status==QMSTATUS_FILE_TRUNCATED) 00842 return "File truncated"; 00843 else 00844 return "Unknown"; 00845 } 00846 00847 00848 00849 // ======================= 00850 // Wavefunction signatures 00851 // ======================= 00852 00853 00864 int QMData::assign_wavef_id(int type, int spin, int exci, char *info, 00865 wavef_signa_t *(&signa_curts), 00866 int &num_signa_curts) { 00867 int j, idtag=-1; 00868 00869 // Check if a wavefunction with the same signature exists 00870 // already in the global signature list. In this case we 00871 // return the corresponding idtag (which will cause 00872 // QMTimestep::add_wavefunction to overwrite the existing 00873 // wavefunction with that idtag) . 00874 for (j=0; j<num_wavef_signa; j++) { 00875 if (wavef_signa[j].type==type && 00876 wavef_signa[j].spin==spin && 00877 wavef_signa[j].exci==exci && 00878 (info && !strncmp(wavef_signa[j].info, info, QMDATA_BUFSIZ))) { 00879 idtag = j; 00880 } 00881 } 00882 00883 // Check if we have the same signature in the current timestep 00884 int duplicate = 0; 00885 for (j=0; j<num_signa_curts; j++) { 00886 if (signa_curts[j].type==type && 00887 signa_curts[j].spin==spin && 00888 signa_curts[j].exci==exci && 00889 (info && !strncmp(signa_curts[j].info, info, QMDATA_BUFSIZ))) { 00890 duplicate = 1; 00891 } 00892 } 00893 00894 // Add a new signature for the current timestep 00895 if (!signa_curts) { 00896 signa_curts = (wavef_signa_t *)calloc(1, sizeof(wavef_signa_t)); 00897 } else { 00898 signa_curts = (wavef_signa_t *)realloc(signa_curts, 00899 (num_signa_curts+1)*sizeof(wavef_signa_t)); 00900 } 00901 signa_curts[num_signa_curts].type = type; 00902 signa_curts[num_signa_curts].spin = spin; 00903 signa_curts[num_signa_curts].exci = exci; 00904 if (!info) 00905 signa_curts[num_signa_curts].info[0] = '0円'; 00906 else 00907 strncpy(signa_curts[num_signa_curts].info, info, QMDATA_BUFSIZ-1); 00908 num_signa_curts++; 00909 00910 // Add new wavefunction ID tag in case this signature wasn't 00911 // found at all or in case we have a duplicate in the current 00912 // timestep. 00913 // If in a single timestep two wavefunction with the same type 00914 // are sent by the molfile_plugin then we assume that the plugin 00915 // considers them different wavefunctions. Our categories are 00916 // just not sufficient to distinguish them. 00917 if (idtag<0 || duplicate) { 00918 if (!wavef_signa) { 00919 wavef_signa = (wavef_signa_t *)calloc(1, sizeof(wavef_signa_t)); 00920 } else { 00921 wavef_signa = (wavef_signa_t *)realloc(wavef_signa, 00922 (num_wavef_signa+1)*sizeof(wavef_signa_t)); 00923 } 00924 wavef_signa[num_wavef_signa].type = type; 00925 wavef_signa[num_wavef_signa].spin = spin; 00926 wavef_signa[num_wavef_signa].exci = exci; 00927 wavef_signa[num_wavef_signa].max_avail_orbs = 0; 00928 wavef_signa[num_wavef_signa].orbids = NULL; 00929 if (!info) 00930 wavef_signa[num_wavef_signa].info[0] = '0円'; 00931 else 00932 strncpy(wavef_signa[num_wavef_signa].info, info, QMDATA_BUFSIZ-1); 00933 idtag = num_wavef_signa; 00934 num_wavef_signa++; 00935 } 00936 00937 //printf("idtag=%d (%d, %d, %d, %s)\n", idtag, type, spin, exci, info); 00938 00939 return idtag; 00940 } 00941 00942 00943 // Find the wavefunction ID tag by comparing 00944 // type, spin, and excitation with the signatures 00945 // of existing wavefunctions 00946 // Returns -1 if no such wavefunction exists. 00947 int QMData::find_wavef_id_from_gui_specs(int guitype, int spin, int exci) { 00948 int i, idtag = -1; 00949 for (i=0; i<num_wavef_signa; i++) { 00950 if (spin==wavef_signa[i].spin && 00951 exci==wavef_signa[i].exci) { 00952 if (compare_wavef_guitype_to_type(guitype, wavef_signa[i].type)) { 00953 idtag = i; 00954 } 00955 00956 } 00957 } 00958 //if (idtag<0) { printf("Couldn't find_wavef_id_from_gui_specs: guitype=%d \n", guitype); } 00959 return idtag; 00960 } 00961 00962 00964 int QMData::has_wavef_guitype(int guitype) { 00965 int i; 00966 for (i=0; i<num_wavef_signa; i++) { 00967 if (compare_wavef_guitype_to_type(guitype, wavef_signa[i].type)) { 00968 return 1; 00969 } 00970 } 00971 return 0; 00972 } 00973 00974 int QMData::compare_wavef_guitype_to_type(int guitype, int type) { 00975 if ((guitype==GUI_WAVEF_TYPE_CANON && type==WAVE_CANON) || 00976 (guitype==GUI_WAVEF_TYPE_GEMINAL && type==WAVE_GEMINAL) || 00977 (guitype==GUI_WAVEF_TYPE_MCSCFNAT && type==WAVE_MCSCFNAT) || 00978 (guitype==GUI_WAVEF_TYPE_MCSCFOPT && type==WAVE_MCSCFOPT) || 00979 (guitype==GUI_WAVEF_TYPE_CINAT && type==WAVE_CINATUR) || 00980 (guitype==GUI_WAVEF_TYPE_OTHER && type==WAVE_UNKNOWN) || 00981 (guitype==GUI_WAVEF_TYPE_LOCAL && 00982 (type==WAVE_BOYS || type==WAVE_RUEDEN || type==WAVE_PIPEK))) { 00983 return 1; 00984 } 00985 return 0; 00986 } 00987 00988 00990 int QMData::has_wavef_spin(int spin) { 00991 int i; 00992 for (i=0; i<num_wavef_signa; i++) { 00993 if (wavef_signa[i].spin==spin) return 1; 00994 } 00995 return 0; 00996 } 00997 00998 01000 int QMData::has_wavef_exci(int exci) { 01001 int i; 01002 for (i=0; i<num_wavef_signa; i++) { 01003 if (wavef_signa[i].exci==exci) return 1; 01004 } 01005 return 0; 01006 } 01007 01008 01011 int QMData::has_wavef_signa(int type, int spin, int exci) { 01012 int i; 01013 for (i=0; i<num_wavef_signa; i++) { 01014 if (wavef_signa[i].type==type && 01015 wavef_signa[i].exci==exci && 01016 wavef_signa[i].spin==spin) return 1; 01017 } 01018 return 0; 01019 } 01020 01021 01024 int QMData::get_highest_excitation(int guitype) { 01025 int i, highest=0; 01026 for (i=0; i<num_wavef_signa; i++) { 01027 if (wavef_signa[i].exci>highest && 01028 compare_wavef_guitype_to_type(guitype, wavef_signa[i].type)) { 01029 highest = wavef_signa[i].exci; 01030 } 01031 } 01032 return highest; 01033 } 01034 01035 01036 // ===================================================== 01037 // Functions dealing with the list of orbitals available 01038 // for a given wavefunction. Needed for the GUI. 01039 // ===================================================== 01040 01041 01042 // Merge the provided list of orbital IDs with the existing 01043 // list of available orbitals. Available orbitals are the union 01044 // of all orbital IDs for the wavefunction with ID iwavesig 01045 // occuring throughout the trajectory. 01046 void QMData::update_avail_orbs(int iwavesig, int norbitals, 01047 const int *orbids, const float *orbocc) { 01048 int i, j; 01049 01050 // Signature of wavefunction 01051 wavef_signa_t *cursig = &wavef_signa[iwavesig]; 01052 01053 for (i=0; i<norbitals; i++) { 01054 int found = 0; 01055 for (j=0; j<cursig->max_avail_orbs; j++) { 01056 if (cursig->orbids[j]==orbids[i]) { 01057 found = 1; 01058 break; 01059 } 01060 } 01061 if (!found) { 01062 if (!cursig->orbids) { 01063 cursig->orbids = (int *)calloc(1, sizeof(int)); 01064 cursig->orbocc = (float*)calloc(1, sizeof(float)); 01065 } else { 01066 cursig->orbids = (int *)realloc(cursig->orbids, 01067 (cursig->max_avail_orbs+1)*sizeof(int)); 01068 cursig->orbocc = (float*)realloc(cursig->orbocc, 01069 (cursig->max_avail_orbs+1)*sizeof(float)); 01070 } 01071 cursig->orbids[cursig->max_avail_orbs] = orbids[i]; 01072 cursig->orbocc[cursig->max_avail_orbs] = orbocc[i]; 01073 cursig->max_avail_orbs++; 01074 } 01075 } 01076 // printf("iwavesig=%d, ", iwavesig); 01077 // for (j=0; j<cursig->max_avail_orbs; j++) { 01078 // printf("%d %.2f\n",cursig->orbids[j], cursig->orbocc[j]); 01079 // } 01080 // printf("\n"); 01081 } 01082 01083 01089 int QMData::get_max_avail_orbitals(int iwavesig) { 01090 if (iwavesig<0 || iwavesig>=num_wavef_signa) return -1; 01091 return wavef_signa[iwavesig].max_avail_orbs; 01092 } 01093 01094 01098 int QMData::get_avail_orbitals(int iwavesig, int *(&orbids)) { 01099 if (iwavesig<0 || iwavesig>=num_wavef_signa) return 0; 01100 01101 int i; 01102 for (i=0; i<wavef_signa[iwavesig].max_avail_orbs; i++) { 01103 orbids[i] = wavef_signa[iwavesig].orbids[i]; 01104 } 01105 return 1; 01106 } 01107 01108 01112 int QMData::get_avail_occupancies(int iwavesig, float *(&orbocc)) { 01113 if (iwavesig<0 || iwavesig>=num_wavef_signa) return 0; 01114 01115 int i; 01116 for (i=0; i<wavef_signa[iwavesig].max_avail_orbs; i++) { 01117 orbocc[i] = wavef_signa[iwavesig].orbocc[i]; 01118 } 01119 return 1; 01120 } 01121 01122 01128 int QMData::get_orbital_label_from_gui_index(int iwavesig, int iorb) { 01129 if (iwavesig<0 || iwavesig>=num_wavef_signa || 01130 iorb<0 ||iorb>=wavef_signa[iwavesig].max_avail_orbs) 01131 return -1; 01132 return wavef_signa[iwavesig].orbids[iorb]; 01133 } 01134 01137 int QMData::has_orbital(int iwavesig, int orbid) { 01138 if (iwavesig<0 || iwavesig>=num_wavef_signa) return 0; 01139 01140 int i; 01141 for (i=0; i<wavef_signa[iwavesig].max_avail_orbs; i++) { 01142 if (orbid==wavef_signa[iwavesig].orbids[i]) return 1; 01143 } 01144 return 0; 01145 01146 } 01147 01148 int QMData::expand_atompos(const float *atompos, 01149 float *(&expandedpos)) { 01150 int i, at; 01151 expandedpos = new float[3*num_wave_f]; 01152 01153 int t = 0; 01154 // loop over all the QM atoms 01155 for (at=0; at<num_atoms; at++) { 01156 int a = atom_types[at]; 01157 float x = atompos[3*at ]*ANGS_TO_BOHR; 01158 float y = atompos[3*at+1]*ANGS_TO_BOHR; 01159 float z = atompos[3*at+2]*ANGS_TO_BOHR; 01160 printf("{%.2f %.2f %.2f}\n", x, y, z); 01161 for (i=0; i<basis_set[a].numshells; i++) { 01162 // Loop over the Gaussian primitives of this contracted 01163 // basis function to build the atomic orbital 01164 shell_t *shell = &basis_set[a].shell[i]; 01165 int shelltype = shell->type; 01166 printf("shelltype = %d\n", shelltype); 01167 int l, m, n; 01168 for (n=0; n<=shelltype; n++) { 01169 int mmax = shelltype - n; 01170 for (m=0, l=mmax; m<=mmax; m++, l--) { 01171 expandedpos[3*t ] = x; 01172 expandedpos[3*t+1] = y; 01173 expandedpos[3*t+2] = z; 01174 t++; 01175 } 01176 } 01177 } 01178 } 01179 return 0; 01180 } 01181 01182 01183 int QMData::expand_basis_array(float *(&expandedbasis), int *(&numprims)) { 01184 int i, at; 01185 int num_prim_total = 0; 01186 for (at=0; at<num_atoms; at++) { 01187 int a = atom_types[at]; 01188 for (i=0; i<basis_set[a].numshells; i++) { 01189 num_prim_total += basis_set[a].shell[i].numprims * 01190 basis_set[a].shell[i].num_cart_func; 01191 } 01192 } 01193 01194 numprims = new int[num_wave_f]; 01195 expandedbasis = new float[2*num_prim_total]; 01196 int t=0, ifunc=0; 01197 // loop over all the QM atoms 01198 for (at=0; at<num_atoms; at++) { 01199 int a = atom_types[at]; 01200 printf("atom %d\n", at); 01201 01202 for (i=0; i<basis_set[a].numshells; i++) { 01203 // Loop over the Gaussian primitives of this contracted 01204 // basis function to build the atomic orbital 01205 shell_t *shell = &basis_set[a].shell[i]; 01206 int maxprim = shell->numprims; 01207 int shelltype = shell->type; 01208 printf("shelltype = %d\n", shelltype); 01209 int l, m, n, icart=0; 01210 for (n=0; n<=shelltype; n++) { 01211 int mmax = shelltype - n; 01212 for (m=0, l=mmax; m<=mmax; m++, l--) { 01213 printf("lmn=%d%d%d %d%d%d\n", l, m, n, 01214 angular_momentum[3*ifunc], 01215 angular_momentum[3*ifunc+1], 01216 angular_momentum[3*ifunc+2]); 01217 numprims[ifunc++] = maxprim; 01218 float normfac = shell->norm_fac[icart]; 01219 icart++; 01220 int prim; 01221 for (prim=0; prim<maxprim; prim++) { 01222 expandedbasis[2*t ] = shell->prim[prim].expon; 01223 //expandedbasis[2*t+1] = normfac*shell->prim[prim].coeff; 01224 expandedbasis[2*t+1] = normfac*shell->basis[2*prim+1]; 01225 printf("expon=%f coeff=%f numprims=%d\n", expandedbasis[2*t], expandedbasis[2*t+1], numprims[ifunc-1]); 01226 t++; 01227 } 01228 } 01229 } 01230 } 01231 } 01232 return 1; 01233 } 01234 01235 01236 #define MIN(X,Y) (((X)<(Y))? (X) : (Y)) 01237 #define MAX(X,Y) (((X)>(Y))? (X) : (Y)) 01238 01239 01240 // 1-electron integral evaluation 01241 // ============================== 01242 // 01243 // Below are function for the evaluation of the overlap 01244 // integrals which can be used as a template for developing 01245 // other 1e-integral such as the kinetic energy integrals 01246 // or the nuclear attraction integrals. 01247 // The implemetation for the overlap matrix closely follows 01248 // "Fundamentals of Molecular Integrals Evaluation" by 01249 // Fermann & Valeev (lecture notes): 01250 // http://www.files.chem.vt.edu/chem-dept/valeev/docs/ints.pdf 01251 // See also "Molecular Integrals", Huzinaga 1967, p.68. 01252 // and "Modern Quantum Chemistry", Szabo & Ostlund 1989, 01253 // Appendix A. 01254 // 01255 // The molecular integral scheme used below represents the 01256 // direct computation of the analytical integrals and is not 01257 // optimized. 01258 // Improved algorithms for molecular integrals include 01259 // * the Obara-Saika scheme (exploiting the translational 01260 // and horizontal recurrence relation for cartesian overlap 01261 // integrals) 01262 // * the McMurchie-Davidson scheme (expanding the overlap 01263 // distribution in Hermite Gaussians). 01264 01265 01266 // Binomial coefficient 01267 static int binomial(int n, int k) { 01268 return fac(n)/(fac(k)*fac(n-k)); 01269 } 01270 01271 01272 // Expression (2.46) for f_k on page 13 of Fermann & Valeev. 01273 float overlap_f(int k, int l1, int l2, float PAx, float PBx) { 01274 int q; 01275 float f = 0.f; 01276 for (q=MAX(-k, k-2*l2); q<=MIN(k, 2*l1-k); q+=2) { 01277 int i = (k+q)/2; 01278 int j = (k-q)/2; 01279 f += (float) (binomial(l1, i)*binomial(l2, j)*pow(PAx, l1-i)*pow(PBx, l2-j)); 01280 } 01281 return f; 01282 } 01283 01284 01285 // Expression (3.15) for Ix on page 16 of Fermann & Valeev. 01286 float overlap_I(int l1, int l2, float PAx, float PBx, float gamma) { 01287 int i; 01288 float Ix = 0.f; 01289 for (i=0; i<=(l1+l2)/2; i++) { 01290 Ix += (float) (overlap_f(2*i, l1, l2, PAx, PBx) * doublefac(2*i-1)/pow(2*gamma,i) * sqrtf(float(VMD_PI)/gamma)); 01291 } 01292 return Ix; 01293 } 01294 01295 01296 // Expression (3.12) for S12 on page 16 of Fermann & Valeev: 01297 // Overlap of primitive functions with arbitrary angular momentum 01298 float overlap_S12(int l1, int m1, int n1, int l2, int m2, int n2, 01299 float alpha, float beta, float rAB2, 01300 const float *PA, const float *PB) { 01301 float gamma = alpha+beta; 01302 01303 float Ix = overlap_I(l1, l2, PA[0], PB[0], gamma); 01304 float Iy = overlap_I(m1, m2, PA[1], PB[1], gamma); 01305 float Iz = overlap_I(n1, n2, PA[2], PB[2], gamma); 01306 //printf(" I = {%f %f %f}\n", Ix, Iy, Iz); 01307 return (float) exp(-alpha*beta*rAB2/gamma)*Ix*Iy*Iz; 01308 } 01309 01310 01311 // Overlap of contracted functions with arbitrary angular momentum 01312 // Summing up contributions from oververlap integrals between all 01313 // combinations of primitives of the two contrated functions. 01314 float overlap_S12_contracted(const float *basis1, const float *basis2, 01315 int numprim1, int numprim2, 01316 int l1, int m1, int n1, int l2, int m2, int n2, 01317 const float *A, const float *B) { 01318 int i, j; 01319 float S12_contracted = 0.f; 01320 float sqrtpigamma3 = 1.0f; 01321 float rAB2 = distance2(A, B); 01322 int SS = 0; 01323 01324 if (l1+m1+n1==0 && l2+m2+n2==0) SS = 1; 01325 01326 for (i=0; i<numprim1; i++) { 01327 float alpha = basis1[2*i]; 01328 for (j=0; j<numprim2; j++) { 01329 float beta = basis2[2*j]; 01330 float gamma = alpha+beta; 01331 01332 // P = (alpha*A + beta*B)/gamma 01333 float P[3], PA[3], PB[3]; 01334 vec_scale(P, alpha, A); // P = alpha*A 01335 vec_scaled_add(P, beta, B); // P += beta*B 01336 vec_scale(P, 1/gamma, P); // P = P/gamma 01337 vec_sub(PA, P, A); 01338 vec_sub(PB, P, B); 01339 01340 switch (SS) { 01341 case 1: 01342 // We can use a simpler formula for S-S overlap 01343 sqrtpigamma3 = sqrtf(float(VMD_PI)/gamma); 01344 sqrtpigamma3 = sqrtpigamma3*sqrtpigamma3*sqrtpigamma3; 01345 S12_contracted += float(basis1[2*i+1]*basis2[2*j+1]*exp(-alpha*beta*rAB2/gamma)*sqrtpigamma3); 01346 break; 01347 default: 01348 float S = basis1[2*i+1]*basis2[2*j+1]* 01349 overlap_S12(l1, m1, n1, l2, m2, n2, 01350 alpha, beta, rAB2, PA, PB); 01351 //printf(" prim %d,%d: co %f %f ex %f %f\n", 01352 // i+1, j+1, 01353 // basis1[2*i+1], basis2[2*j+1], 01354 // basis1[2*i], basis2[2*j]); 01355 S12_contracted += S; 01356 } 01357 } 01358 } 01359 return S12_contracted; 01360 } 01361 01362 int get_overlap_matrix(int num_wave_f, const float *expandedbasis, 01363 const int *numprims, 01364 const float *atompos, 01365 const int *lmn, 01366 float *overlap_matrix) { 01367 int i, j; 01368 int t1 = 0; 01369 for (i=0; i<num_wave_f; i++) { 01370 const float *basis1 = &expandedbasis[2*t1]; 01371 int numprim1 = numprims[i]; 01372 const float *pos1 = &atompos[3*i]; 01373 int l1 = lmn[3*i ]; 01374 int m1 = lmn[3*i+1]; 01375 int n1 = lmn[3*i+2]; 01376 int t2 = t1; 01377 for (j=i; j<num_wave_f; j++) { 01378 int l2 = lmn[3*j ]; 01379 int m2 = lmn[3*j+1]; 01380 int n2 = lmn[3*j+2]; 01381 const float *basis2 = &expandedbasis[2*t2]; 01382 int numprim2 = numprims[j]; 01383 const float *pos2 = &atompos[3*j]; 01384 printf("%d,%d %d%d%d--%d%d%d %d-%d {%.2f %.2f %.2f} {%.2f %.2f %.2f}\n", 01385 i, j, l1, m1, n1, l2, m2, n2, numprim1, numprim2, 01386 pos1[0], pos1[1], pos1[2], pos2[0], pos2[1], pos2[2]); 01387 01388 float Sij = overlap_S12_contracted(basis1, basis2, numprim1, numprim2, 01389 l1, m1, n1, l2, m2, n2, pos1, pos2); 01390 overlap_matrix[i*num_wave_f+j] = Sij; 01391 overlap_matrix[j*num_wave_f+i] = Sij; 01392 printf(" S12 = %f\n", Sij); 01393 t2 += numprim2; 01394 } 01395 t1 += numprim1; 01396 } 01397 return 0; 01398 } 01399 01400 #if 0 01401 01402 float* QMData::get_overlap_integrals(const float *atompos) { 01403 float *ao_overlap_integrals=NULL; 01404 if (!ao_overlap_integrals) { 01405 ao_overlap_integrals = new float[num_wave_f*num_wave_f]; 01406 int i; 01407 for (i=0; i<num_wave_f*num_wave_f; i++) { 01408 ao_overlap_integrals[i] = 1.f; 01409 } 01410 01411 int at1; 01412 int shell_counter = 0; 01413 int prim_counter = 0; 01414 // loop over all the QM atoms 01415 for (at1=0; at1<num_atoms; at1++) { 01416 int maxshell1 = num_shells_per_atom[at1]; 01417 float x1 = atompos[3*at1 ]*ANGS_TO_BOHR; 01418 float y1 = atompos[3*at1+1]*ANGS_TO_BOHR; 01419 float z1 = atompos[3*at1+2]*ANGS_TO_BOHR; 01420 int shell1; 01421 for (shell1=0; shell1 < maxshell1; shell1++) { 01422 // Loop over the Gaussian primitives of this contracted 01423 // basis function to build the atomic orbital 01424 int numprims = num_prim_per_shell[shell_counter]; 01425 int shelltype = shell_types[shell_counter]; 01426 int prim1; 01427 for (prim1=0; prim1 < numprims; prim1++) { 01428 float exponent1 = basis_array[prim_counter ]; 01429 float contract_coeff1 = basis_array[prim_counter + 1]; 01430 int at2; 01431 for (at2=0; at2<num_atoms; at2++) { 01432 int maxshell2 = num_shells_per_atom[at2]; 01433 float x2 = atompos[3*at2 ]*ANGS_TO_BOHR; 01434 float y2 = atompos[3*at2+1]*ANGS_TO_BOHR; 01435 float z2 = atompos[3*at2+2]*ANGS_TO_BOHR; 01436 float dx = x2-x1; 01437 float dy = y2-y1; 01438 float dz = z2-z1; 01439 float dist2 = dx*dx + dy*dy + dz*dz; 01440 int shell2; 01441 for (shell2=0; shell2 < maxshell2; shell2++) { 01442 int numprims2 = num_prim_per_shell[shell_counter]; 01443 int shelltype = shell_types[shell_counter]; 01444 int prim2; 01445 for (prim2=0; prim2 < numprims; prim2++) { 01446 float exponent2 = basis_array[prim_counter ]; 01447 float contract_coeff2 = basis_array[prim_counter + 1]; 01448 prim_counter += 2; 01449 } 01450 } 01451 } 01452 } 01453 } 01454 } 01455 } 01456 return ao_overlap_integrals; 01457 } 01458 #endif 01459 01460 01463 void QMData::compute_overlap_integrals(Timestep *ts, 01464 const float *expandedbasis, 01465 const int *numprims, 01466 float *(&overlap_matrix)) { 01467 float *expandedpos = NULL; 01468 expand_atompos(ts->pos, expandedpos); 01469 01470 01471 overlap_matrix = new float[num_wave_f*num_wave_f]; 01472 memset(overlap_matrix, 0, num_wave_f*num_wave_f*sizeof(float)); 01473 01474 get_overlap_matrix(num_wave_f, expandedbasis, numprims, 01475 expandedpos, 01476 angular_momentum, overlap_matrix); 01477 delete [] expandedpos; 01478 } 01479 01480 01481 // matrix multiplication 01482 static void mm_mul(const float *a, int awidth, int aheight, 01483 const float *b, int bwidth, int bheight, 01484 float *(&c)) { 01485 if (awidth!=bheight) 01486 printf("mm_mul(): incompatible sizes %d,%d\n", awidth, bheight); 01487 c = new float[aheight*bwidth]; 01488 for (int i=0; i<aheight; i++) { 01489 for (int j=0; j<bwidth; j++) { 01490 float cc = 0.f; 01491 for (int k=0; k<awidth; k++) 01492 cc += a[i*awidth+k]*b[k*bwidth+j]; 01493 c[i*bwidth+j] = cc; 01494 } 01495 } 01496 } 01497 01498 01499 01500 int QMData::mullikenpop(Timestep *ts, int iwavesig, 01501 const float *expandedbasis, 01502 const int *numprims) { 01503 if (iwavesig<0 || iwavesig>=num_wavef_signa || !ts) return 0; 01504 01505 int iwave = ts->qm_timestep->get_wavef_index(iwavesig); 01506 const Wavefunction *wave = ts->qm_timestep->get_wavefunction(iwave); 01507 01508 01509 float *S; 01510 compute_overlap_integrals(ts, expandedbasis, numprims, S); 01511 01512 int numcoeffs = wave->get_num_coeffs(); 01513 01514 float *P; 01515 wave->population_matrix(S, P); 01516 01517 int i,j; 01518 for (i=0; i<numcoeffs; i++) { 01519 for (j=0; j<numcoeffs; j++) { 01520 printf("P[%d,%d]=%f\n", i, j, P[i*numcoeffs+j]); 01521 } 01522 } 01523 01524 float *PA; 01525 wave->population_matrix(this, 2, S, PA); 01526 //wave->density_matrix(this, 0, PA); 01527 01528 for (i=0; i<get_num_wavecoeff_per_atom(2); i++) { 01529 for (j=0; j<numcoeffs; j++) { 01530 printf("PA[%d,%d]=%f\n", i, j, PA[i*numcoeffs+j]); 01531 } 01532 } 01533 delete [] PA; 01534 01535 float *GOP = new float[numcoeffs]; 01536 for (i=0; i<numcoeffs; i++) { 01537 GOP[i] = 0.f; 01538 for (j=0; j<numcoeffs; j++) { 01539 GOP[i] += P[i*numcoeffs+j]; 01540 } 01541 printf("GOP[%d] = %f\n", i, GOP[i]); 01542 } 01543 01544 float *GAP = new float[num_atoms]; 01545 int coeff_count = 0; 01546 int a; 01547 for (a=0; a<num_atoms; a++) { 01548 int num_cart_func = get_num_wavecoeff_per_atom(a); 01549 GAP[a] = 0.f; 01550 for (i=0; i<num_cart_func; i++) { 01551 GAP[a] += GOP[coeff_count++]; 01552 } 01553 printf("GAP[%d] = %f\n", a, GAP[a]); 01554 } 01555 01556 float *D; 01557 wave->density_matrix(D); 01558 float *DS = new float[numcoeffs*numcoeffs]; 01559 01560 mm_mul(D, numcoeffs, numcoeffs, S, numcoeffs, numcoeffs, DS); 01561 01562 float Nelec = 0.f; 01563 for (i=0; i<numcoeffs; i++) { 01564 printf("DS[%d,%d] = %f\n", i, i, DS[i*numcoeffs+i]); 01565 Nelec += DS[i*numcoeffs+i]; 01566 } 01567 01568 printf("Nelec=%f\n", Nelec); 01569 01570 delete [] S; 01571 delete [] P; 01572 delete [] D; 01573 delete [] DS; 01574 01575 return 1; 01576 } 01577 01578 01579 // ======================================================== 01580 // Orbital localization 01581 // ======================================================== 01582 01601 int QMData::orblocalize(Timestep *ts, int iwavesig, 01602 const float *expandedbasis, 01603 const int *numprims) { 01604 if (iwavesig<0 || iwavesig>=num_wavef_signa || !ts) return 0; 01605 01606 int iwave = ts->qm_timestep->get_wavef_index(iwavesig); 01607 Wavefunction *wave = ts->qm_timestep->get_wavefunction(iwave); 01608 01609 01610 float *S; 01611 compute_overlap_integrals(ts, expandedbasis, numprims, S); 01612 01613 int i, j; 01614 for (i=0; i<num_wave_f; i++) { 01615 for (j=i; j<num_wave_f; j++) { 01616 printf("S12[%d,%d] = %f\n", i, j, S[i*num_wave_f+j]); 01617 } 01618 } 01619 int numoccorbs = wave->get_num_occupied_double(); 01620 // int numorbs = wave->get_num_orbitals(); 01621 float *C = new float[numoccorbs*num_wave_f]; 01622 const float *Ccanon = wave->get_coeffs(); 01623 // for (i=0; i<num_wave_f; i++) { 01624 // memcpy(&C[i*numoccorbs], &Ccanon[i*numorbs], 01625 // numoccorbs*sizeof(float)); 01626 // } 01627 memcpy(C, Ccanon, numoccorbs*num_wave_f*sizeof(float)); 01628 double D = mean_localization_measure(numoccorbs, C, S); 01629 printf("Delocalization D=%f \n", D); 01630 01631 double Dold = D; 01632 int iter; 01633 for (iter=0; iter<20; iter++) { 01634 // Find that orbital pair for which 2x2 maximization of D 01635 // yields the greatest increase in total delocalization D: 01636 // deltaD = Dmax(ui,uj) - D(phii,phij) 01637 // = Aij + sqrt(Aij^2 + Bij^2) 01638 // where 01639 // ui, uj = rotated orbital pair 01640 // phii, phij = non-rotated orbital pair 01641 // Aij = gross Mulliken population of orbital pair i,j 01642 // Bij = see localization_rotation_angle() 01643 01644 double deloc, maxchange = 0.0; 01645 int maxdelocorb1 = 0; 01646 int maxdelocorb2 = 0; 01647 for (i=0; i<numoccorbs; i++) { 01648 for (j=i+1; j<numoccorbs; j++) { 01649 deloc = pair_localization_measure(numoccorbs, i, j, C, S); 01650 double change = localization_orbital_change(i, j, C, S); 01651 printf("deloc[%d,%d] = %f change = %.7f\n", i, j, deloc, change); 01652 if (change>maxchange) { 01653 maxchange = change; 01654 maxdelocorb1 = i; 01655 maxdelocorb2 = j; 01656 } 01657 } 01658 } 01659 if (maxchange<0.000001) { 01660 printf("maxchange = %f\n",maxchange); 01661 break; 01662 } 01663 01664 double gamma = localization_rotation_angle(C, S, maxdelocorb1, 01665 maxdelocorb2); 01666 01667 printf("Rotating orbitals %d,%d by %f\n", maxdelocorb1, maxdelocorb2, gamma); 01668 01669 rotate_2x2_orbitals(C, maxdelocorb1, maxdelocorb2, gamma); 01670 01671 D = mean_localization_measure(numoccorbs, C, S); 01672 printf("Delocalization after rot[%d] D=%f \n", iter, D); 01673 01674 if (fabs(D-Dold)<0.000001) break; 01675 01676 Dold = D; 01677 } 01678 01679 int b; 01680 int ncol = 5; 01681 for (b=0; b<=(numoccorbs-1)/ncol; b++) { 01682 for (j=0; j<num_wave_f; j++) { 01683 printf("%4d ", j); 01684 for (i=0; i<ncol && b*ncol+i<numoccorbs; i++) { 01685 printf("% f ", C[(b*ncol+i)*num_wave_f+j]); 01686 } 01687 printf("\n"); 01688 } 01689 printf("\n"); 01690 } 01691 01692 // XXX dead code? 01693 // float occupancies[numoccorbs]; 01694 // for (i=0; i<numoccorbs; i++) occupancies[i] = 2; 01695 01696 int num_signa_ts = 0; 01697 wavef_signa_t *signa_ts = NULL; 01698 ts->qm_timestep->add_wavefunction(this, num_wave_f, numoccorbs, 01699 C, NULL, NULL, NULL, 0.0, 01700 WAVE_PIPEK, wave->get_spin(), 01701 wave->get_excitation(), 01702 wave->get_multiplicity(), 01703 "Pipek-Mezey localization", 01704 signa_ts, num_signa_ts); 01705 free(signa_ts); 01706 01707 delete [] S; 01708 delete [] C; 01709 01710 return 1; 01711 } 01712 01713 double QMData::localization_orbital_change(int orbid1, int orbid2, 01714 const float *C, 01715 const float *S) { 01716 double Ast = 0.0; 01717 double Bst = 0.0; 01718 int a; 01719 for (a=0; a<num_atoms; a++) { 01720 double QAs = gross_atom_mulliken_pop(C, S, a, orbid1); 01721 double QAt = gross_atom_mulliken_pop(C, S, a, orbid2); 01722 double QAst = orb_pair_gross_atom_mulliken_pop(C, S, a, orbid1, orbid2); 01723 double QAdiff = (QAs-QAt); 01724 Ast += QAst*QAst - 0.25*QAdiff*QAdiff; 01725 Bst += QAst*QAdiff; 01726 } 01727 01728 return Ast + sqrtf(float(Ast*Ast+Bst*Bst)); 01729 } 01730 01731 // Return the delocalization of an orbital pair. 01732 float QMData::pair_localization_measure(int num_orbitals, 01733 int orbid1, int orbid2, 01734 const float *C, 01735 const float *S) { 01736 int a; 01737 float deloc1 = 0.0; 01738 float deloc2 = 0.0; 01739 for (a=0; a<num_atoms; a++) { 01740 // gross Mulliken population of orbital <orbid1> on atom <a> 01741 float mullpop = float(gross_atom_mulliken_pop(C, S, a, orbid1)); 01742 //printf("mullpop[%d,%d] = %f\n", a, orbid1, mullpop); 01743 deloc1 += mullpop*mullpop; 01744 } 01745 for (a=0; a<num_atoms; a++) { 01746 // gross Mulliken population of orbital <orbid2> on atom <a> 01747 float mullpop = float(gross_atom_mulliken_pop(C, S, a, orbid2)); 01748 //printf("mullpop[%d,%d] = %f\n", a, orbid2, mullpop); 01749 deloc2 += mullpop*mullpop; 01750 } 01751 return deloc1+deloc2; 01752 } 01753 01754 01755 // Return the mean over-all localization. 01756 double QMData::mean_localization_measure(int num_orbitals, 01757 const float *C, 01758 const float *S) { 01759 double Dinv = 0.0; 01760 int a, orbid=0; 01761 for (orbid=0; orbid<num_orbitals; orbid++) { 01762 double deloc = 0.0; 01763 01764 for (a=0; a<num_atoms; a++) { 01765 // gross Mulliken population of orbital <orbid> on atom <a> 01766 float mullpop = float(gross_atom_mulliken_pop(C, S, a, orbid)); 01767 01768 //printf("mullpop[%d,%d] = %f\n", a, orbid, 2*mullpop); 01769 deloc += mullpop*mullpop; 01770 } 01771 //deloc[orbid] = deloc; 01772 //printf("deloc[%d] = %f\n", orbid, deloc); 01773 Dinv += deloc; 01774 } 01775 //Dinv /= num_orbitals; 01776 01777 return Dinv; 01778 } 01779 01782 void QMData::rotate_2x2_orbitals(float *C, int orbid1, int orbid2, 01783 double gamma) { 01784 int num_coeffs = num_wave_f; 01785 float *Corb1 = &C[orbid1*num_coeffs]; 01786 float *Corb2 = &C[orbid2*num_coeffs]; 01787 double singamma, cosgamma; 01788 sincos(gamma, &singamma, &cosgamma); 01789 int i; 01790 for (i=0; i<num_coeffs; i++) { 01791 double tmp = cosgamma*Corb1[i] + singamma*Corb2[i]; 01792 Corb2[i] = float(-singamma*Corb1[i] + cosgamma*Corb2[i]); 01793 Corb1[i] = float(tmp); 01794 } 01795 } 01796 01803 double QMData::localization_rotation_angle(const float *C, const float *S, 01804 int orbid1, int orbid2) { 01805 double Ast = 0.0; 01806 double Bst = 0.0; 01807 int a; 01808 for (a=0; a<num_atoms; a++) { 01809 double QAs = gross_atom_mulliken_pop(C, S, a, orbid1); 01810 double QAt = gross_atom_mulliken_pop(C, S, a, orbid2); 01811 double QAst = orb_pair_gross_atom_mulliken_pop(C, S, a, orbid1, orbid2); 01812 double QAdiff = (QAs-QAt); 01813 Ast += QAst*QAst - 0.25*QAdiff*QAdiff; 01814 Bst += QAst*QAdiff; 01815 } 01816 #if 0 01817 double T = 0.25*atan2(4.0*Bst, 4.0*Ast); 01818 while (T>0) { 01819 double sig = 1.0; 01820 if (T>0) sig = -1.0; 01821 T += sig*0.25*VMD_PI; 01822 printf("atan(Bst,Ast) = %f\n", T); 01823 if (fabs(T)<0.00001 || fabs(fabs(T)-0.25*VMD_PI)<0.00001) break; 01824 } 01825 return T; 01826 #endif 01827 01828 double sign = 1.0; 01829 if (Bst<0.f) sign = -1.0; 01830 01831 return sign*0.25*acos(-Ast/sqrt(Ast*Ast+Bst*Bst)); 01832 } 01833 01834 01835 // Return the gross Mulliken population of orbital <orbid> 01836 // on atom <atomid>. 01837 // 01838 // QA = sum_i(sum_j(C_orb,i * C_orb,j * S_ij)) 01839 // 01840 // where i runs over all coefficients related to the atom <atomid> 01841 // and j runs over all coefficients. 01842 01843 // Input arrays: 01844 // C = wavefunction coefficients 01845 // S = overlap matrix 01846 double QMData::gross_atom_mulliken_pop(const float *C, const float *S, 01847 int atomid, int orbid) const { 01848 int num_coeffs = num_wave_f; 01849 int first_coeff = get_wave_offset(atomid, 0); 01850 const float *Corb = &C[orbid*num_coeffs]; 01851 int num_atom_coeffs = get_num_wavecoeff_per_atom(atomid); 01852 double QA = 0.0; 01853 int i, j; 01854 for (i=first_coeff; i<first_coeff+num_atom_coeffs; i++) { 01855 double Corbi = Corb[i]; 01856 for (j=0; j<num_coeffs; j++) { 01857 QA += Corbi*Corb[j]*S[i*num_coeffs+j]; 01858 } 01859 } 01860 01861 return QA; 01862 } 01863 01864 01865 // Return the gross Mulliken population of orbital pair 01866 //<orbid1>/<orbid2> on atom <atomid>. 01867 // Input arrays: 01868 // C = wavefunction coefficients 01869 // S = overlap matrix 01870 double QMData::orb_pair_gross_atom_mulliken_pop(const float *C, 01871 const float *S, 01872 int atomid, 01873 int orbid1, int orbid2) { 01874 int num_coeffs = num_wave_f; 01875 int first_coeff = get_wave_offset(atomid, 0); 01876 const float *Corb1 = &C[orbid1*num_coeffs]; 01877 const float *Corb2 = &C[orbid2*num_coeffs]; 01878 int num_atom_coeffs = get_num_wavecoeff_per_atom(atomid); 01879 double QA = 0.f; 01880 int i, j; 01881 for (i=first_coeff; i<first_coeff+num_atom_coeffs; i++) { 01882 double C1mu = Corb1[i]; 01883 double C2mu = Corb2[i]; 01884 //printf("Cmu[%d]=%f\n", i, Cmu); 01885 for (j=0; j<num_coeffs; j++) { 01886 double C1nu = Corb1[j]; 01887 double C2nu = Corb2[j]; 01888 //printf(" Cnu=%f S[%d,%d]=%f\n", orbital[j], i, j, S[wave_offset+i*num_coeffs+j]); 01889 QA += (C1nu*C2mu+C1mu*C2nu)*S[i*num_coeffs+j]; 01890 } 01891 } 01892 return 0.5*QA; 01893 } 01894 01895 #if 0 01896 void QMData::gross_atom_mulliken_pop_matrix(const float *C, 01897 const float *S, 01898 int atomid) { 01899 01900 } 01901 #endif 01902 01903 // ======================================================== 01904 // Orbital rendering 01905 // ======================================================== 01906 01907 // Create a new Orbital object and return the pointer. 01908 // User is responsible for deleting. 01909 // The orbital object can be used to compute the wavefunction 01910 // amplitude or electron density for rendering. 01911 Orbital* QMData::create_orbital(int iwave, int orbid, float *pos, 01912 QMTimestep *qmts) { 01913 Orbital *orbital = new Orbital(pos, 01914 qmts->get_wavecoeffs(iwave), 01915 basis_array, basis_set, atom_types, 01916 atom_sort, atom_basis, 01917 (const float**)norm_factors, 01918 num_shells_per_atom, 01919 num_prim_per_shell, shell_types, 01920 num_atoms, num_types, num_wave_f, num_basis, 01921 orbid); 01922 return orbital; 01923 } 01924 01925 01926 01927 // ========================================= 01928 // Currently unused stuff I might need later 01929 // ========================================= 01930 01931 01932 #if 0 // XXX: unused 01933 // XXX these quicksort routines are duplicates of the ones 01934 // in QMTimestep. 01935 static void quicksort(const int *tag, int *A, int p, int r); 01936 static int quickpart(const int *tag, int *A, int p, int r); 01937 01938 // Create an index array *atom_sort that sorts the atoms 01939 // by basis set type (usually that means by atomic number). 01940 void QMData::sort_atoms_by_type() { 01941 int i; 01942 if (atom_sort) delete [] atom_sort; 01943 01944 // Initialize the index array; 01945 atom_sort = new int[num_atoms]; 01946 for (i=0; i<num_atoms; i++) { 01947 atom_sort[i] = i; 01948 } 01949 01950 // Sort index array according to the atom_types 01951 quicksort(atom_types, atom_sort, 0, num_atoms-1); 01952 01953 //int *sorted_types = new int[num_atoms]; 01954 01955 // Copy data over into sorted arrays 01956 //for (i=0; i<num_atoms; i++) { 01957 // sorted_types[i] = atom_types[atom_sort[i]]; 01958 //} 01959 } 01960 01961 // The standard quicksort algorithm except for it doesn't 01962 // sort the data itself but rather sorts array of ints *A 01963 // in the same order as it would sort the integers in array 01964 // *tag. Array *A can then be used to reorder any array 01965 // according to the string tags. 01966 // Example: 01967 // tag: BC DD BB AA --> AA BB BC DD 01968 // index: 0 1 2 3 --> 3 2 0 1 01969 // 01970 static void quicksort(const int* tag, int *idx, int p, int r) { 01971 int q; 01972 if (p < r) { 01973 q=quickpart(tag, idx, p, r); 01974 quicksort(tag, idx, p, q); 01975 quicksort(tag, idx, q+1, r); 01976 } 01977 } 01978 01979 // Partitioning for quicksort. 01980 static int quickpart(const int *tag, int *idx, int p, int r) { 01981 int i, j; 01982 int tmp; 01983 int x = tag[idx[p]]; 01984 i = p-1; 01985 j = r+1; 01986 01987 while (1) { 01988 // Find highest element smaller than idx[p] 01989 do j--; while (tag[idx[j]] > x); 01990 01991 // Find lowest element larger than idx[p] 01992 do i++; while (tag[idx[i]] < x); 01993 01994 if (i < j) { 01995 tmp = idx[i]; 01996 idx[i] = idx[j]; 01997 idx[j] = tmp; 01998 } 01999 else { 02000 return j; 02001 } 02002 } 02003 } 02004 #endif