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: DrawMolItemOrbital.C,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.53 $ $Date: 2021年10月28日 21:12:15 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * A continuation of rendering types from DrawMolItem 00019 * 00020 * This file only contains representations for visualizing QM orbital data 00021 ***************************************************************************/ 00022 #include <math.h> 00023 #include <stdlib.h> 00024 #include <stdio.h> 00025 00026 #include "DrawMolItem.h" 00027 #include "DrawMolecule.h" 00028 #include "BaseMolecule.h" // need volume data definitions 00029 #include "Inform.h" 00030 #include "Isosurface.h" 00031 #include "MoleculeList.h" 00032 #include "Scene.h" 00033 #include "Orbital.h" 00034 #include "QMData.h" 00035 #include "VolumetricData.h" 00036 #include "VMDApp.h" 00037 #include "utilities.h" 00038 #include "WKFUtils.h" // timers 00039 00040 #define MYSGN(a) (((a) > 0) ? 1 : -1) 00041 00042 void DrawMolItem::draw_orbital(int density, int wavefnctype, int wavefncspin, 00043 int wavefncexcitation, int orbid, 00044 float isovalue, 00045 int drawbox, int style, 00046 float gridspacing, int stepsize, int thickness) { 00047 if (!mol->numframes() || gridspacing <= 0.0f) 00048 return; 00049 00050 // only recalculate the orbital grid if necessary 00051 int regenorbital=1; 00052 int useorbgridfromrep = -1; // XXX repID with an existing grid we can reuse 00053 00054 if (density != orbgridisdensity || 00055 wavefnctype != waveftype || 00056 wavefncspin != wavefspin || 00057 wavefncexcitation != wavefexcitation || 00058 orbid != gridorbid || 00059 gridspacing != orbgridspacing || 00060 orbvol == NULL || 00061 needRegenerate & MOL_REGEN || 00062 needRegenerate & SEL_REGEN) { 00063 00064 #if defined(VMDENABLEORBITALGRIDBACKDOOR) 00065 // 00066 // XXX hack to look for an existing orbital grid we can reuse as-is... 00067 // 00068 // This search loop allows the molecular orbital representations 00069 // within the same molecule to reuse any existing 00070 // rep's molecular orbital grid if the orbital ID and various 00071 // grid-specific parameters are all compatible. This optimization 00072 // short-circuits the need for a rep to compute its own grid if 00073 // any other rep already has what it needs. For large QM/MM scenes, 00074 // this optimization can be worth as much as a 2X speedup when 00075 // orbital computation dominates animation performance. 00076 int repcnt = mol->repList.num(); 00077 int r; 00078 for (r=0; r<repcnt && useorbgridfromrep < 0; r++) { 00079 DrawMolItem *dmi = mol->repList[r]; 00080 if (dmi->repNumber != repNumber) { 00081 AtomRep *ar = dmi->atomRep; 00082 if (ar->method() == AtomRep::ORBITAL) { 00083 if (orbid == dmi->gridorbid && 00084 wavefnctype == dmi->waveftype && 00085 wavefncspin == dmi->wavefspin && 00086 wavefncexcitation == dmi->wavefexcitation && 00087 gridspacing == dmi->orbgridspacing && 00088 density == dmi->orbgridisdensity && 00089 dmi->orbvol != NULL) { 00090 // printf("Rep[%d]: orbid %d, rep[%d]: %d\n", repNumber, orbid, r, dmi->gridorbid); 00091 useorbgridfromrep=r; 00092 delete orbvol; 00093 orbvol = dmi->orbvol; // XXX watch out when borrowing orbvol! 00094 } 00095 } 00096 } 00097 } 00098 00099 if (useorbgridfromrep >= 0) 00100 regenorbital=0; 00101 else 00102 #endif 00103 regenorbital=1; 00104 } 00105 00106 double motime=0, voltime=0, gradtime=0; 00107 wkf_timerhandle timer = wkf_timer_create(); 00108 wkf_timer_start(timer); 00109 00110 if (regenorbital) { 00111 // XXX this needs to be fixed so that things like the 00112 // draw multiple frames feature will work correctly for Orbitals 00113 int frame = mol->frame(); // draw currently active frame 00114 const Timestep *ts = mol->get_frame(frame); 00115 00116 if (!ts->qm_timestep || !mol->qm_data || 00117 !mol->qm_data->num_basis || orbid < 1) { 00118 wkf_timer_destroy(timer); 00119 return; 00120 } 00121 00122 // Find the timestep independent wavefunction ID tag 00123 // by comparing type, spin, and excitation with the 00124 // signatures of existing wavefunctions. 00125 int waveid = mol->qm_data->find_wavef_id_from_gui_specs( 00126 wavefnctype, wavefncspin, wavefncexcitation); 00127 00128 // Translate the wavefunction ID into the index the 00129 // wavefunction has in this timestep 00130 int iwave = ts->qm_timestep->get_wavef_index(waveid); 00131 00132 if (iwave<0 || 00133 !ts->qm_timestep->get_wavecoeffs(iwave) || 00134 !ts->qm_timestep->get_num_orbitals(iwave) || 00135 orbid > ts->qm_timestep->get_num_orbitals(iwave)) { 00136 wkf_timer_destroy(timer); 00137 return; 00138 } 00139 00140 // Get the orbital index for this timestep from the orbital ID. 00141 int orbindex = ts->qm_timestep->get_orbital_index_from_id(iwave, orbid); 00142 00143 // Build an Orbital object and prepare to calculate a grid 00144 Orbital *orbital = mol->qm_data->create_orbital(iwave, orbindex, 00145 ts->pos, ts->qm_timestep); 00146 00147 // Set the bounding box of the atom coordinates as the grid dimensions 00148 orbital->set_grid_to_bbox(ts->pos, 3.0, gridspacing); 00149 00150 // XXX needs more testing, can get stuck for certain orbitals 00151 #if 0 00152 // XXX for GPU, we need to only optimize to a stepsize of 4 or more, as 00153 // otherwise doing this actually slows us down rather than speeding up 00154 // orbital.find_optimal_grid(0.01, 4, 8); 00155 // 00156 // optimize: minstep 2, maxstep 8, threshold 0.01 00157 orbital->find_optimal_grid(0.01, 2, 8); 00158 #endif 00159 00160 // Calculate the molecular orbital 00161 orbital->calculate_mo(mol, density); 00162 00163 motime = wkf_timer_timenow(timer); 00164 00165 // query orbital grid origin, dimensions, and axes 00166 const int *numvoxels = orbital->get_numvoxels(); 00167 const float *origin = orbital->get_origin(); 00168 00169 float xaxis[3], yaxis[3], zaxis[3]; 00170 orbital->get_grid_axes(xaxis, yaxis, zaxis); 00171 00172 // build a VolumetricData object for rendering 00173 char dataname[64]; 00174 sprintf(dataname, "molecular orbital %i", orbid); 00175 00176 // update attributes of cached orbital grid 00177 orbgridisdensity = density; 00178 waveftype = wavefnctype; 00179 wavefspin = wavefncspin; 00180 wavefexcitation = wavefncexcitation; 00181 gridorbid = orbid; 00182 orbgridspacing = gridspacing; 00183 delete orbvol; 00184 orbvol = new VolumetricData(dataname, origin, 00185 xaxis, yaxis, zaxis, 00186 numvoxels[0], numvoxels[1], numvoxels[2], 00187 orbital->get_grid_data()); 00188 delete orbital; 00189 00190 voltime = wkf_timer_timenow(timer); 00191 00192 orbvol->compute_volume_gradient(); // calc gradients: smooth vertex normals 00193 00194 gradtime = wkf_timer_timenow(timer); 00195 } // regen the orbital grid... 00196 00197 00198 // draw the newly created VolumetricData object 00199 sprintf(commentBuffer, "Mol[%d] Rep[%d] Orbital", mol->id(), repNumber); 00200 cmdCommentX.putdata(commentBuffer, cmdList); 00201 00202 if (drawbox > 0) { 00203 // don't texture the box if color by volume is active 00204 if (atomColor->method() == AtomColor::VOLUME) { 00205 append(DVOLTEXOFF); 00206 } 00207 // wireframe only? or solid? 00208 if (style > 0 || drawbox == 2) { 00209 draw_volume_box_lines(orbvol); 00210 } else { 00211 draw_volume_box_solid(orbvol); 00212 } 00213 if (atomColor->method() == AtomColor::VOLUME) { 00214 append(DVOLTEXON); 00215 } 00216 } 00217 00218 if ((drawbox == 2) || (drawbox == 0)) { 00219 switch (style) { 00220 case 3: 00221 // shaded points isosurface looping over X-axis, 1 point per voxel 00222 draw_volume_isosurface_lit_points(orbvol, isovalue, stepsize, thickness); 00223 break; 00224 00225 case 2: 00226 // points isosurface looping over X-axis, max of 1 point per voxel 00227 draw_volume_isosurface_points(orbvol, isovalue, stepsize, thickness); 00228 break; 00229 00230 case 1: 00231 // lines implementation, max of 18 line per voxel (3-per triangle) 00232 draw_volume_isosurface_lines(orbvol, isovalue, stepsize, thickness); 00233 break; 00234 00235 case 0: 00236 default: 00237 // trimesh polygonalized surface, max of 6 triangles per voxel 00238 draw_volume_isosurface_trimesh(orbvol, isovalue, stepsize); 00239 break; 00240 } 00241 } 00242 00243 // XXX if we reused the orbital grid from another rep, we have to 00244 // null out orbvol so we don't try and free another reps memory later... 00245 if (useorbgridfromrep >= 0) { 00246 orbvol = NULL; // XXX watch out, un-copy the pointer to the borrowed grid 00247 } 00248 00249 if (regenorbital) { 00250 double surftime = wkf_timer_timenow(timer); 00251 if (surftime > 5) { 00252 char strmsg[1024]; 00253 sprintf(strmsg, "Total MO rep time: %.3f [MO: %.3f vol: %.3f grad: %.3f surf: %.2f]", 00254 surftime, motime, voltime - motime, gradtime - motime, surftime - gradtime); 00255 00256 msgInfo << strmsg << sendmsg; 00257 } 00258 } 00259 00260 wkf_timer_destroy(timer); 00261 } 00262