← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/oeprop/grid_unitvec.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/oeprop/grid_unitvec.cc
 
1
/*! \file
 
2
    \ingroup OEPROP
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#define EXTERN
 
6
#include "includes.h"
 
7
#include "prototypes.h"
 
8
#include "globals.h"
 
9
 
 
10
namespace {
 
11
  int* fock_to_pitzer(int, int*);
 
12
}
 
13
 
 
14
namespace psi { namespace oeprop {
 
15
 
 
16
void grid_parse()
 
17
{
 
18
  int i, errcod, atom;
 
19
  int mo, *fock_mos;
 
20
  char *tmpstring;
 
21
  double xmin, ymin, zmin, xmax, ymax, zmax;    /* molecular dimensions */
 
22
 
 
23
  errcod = ip_data("GRID","%d",&grid,0);
 
24
  if (grid < 0)
 
25
    punt("GRID type must be positive");
 
26
  
 
27
  if (grid == 5 || grid == 6)
 
28
    grid3d = 1;
 
29
  
 
30
  /* which MOs to plot? */
 
31
  if (grid == 5) {
 
32
    
 
33
    read_opdm = 0;
 
34
    
 
35
    if (ip_exist("MO_TO_PLOT",0)) {
 
36
      errcod = ip_count("MO_TO_PLOT",&num_mos_to_plot,0);
 
37
      if (errcod == IPE_NOT_AN_ARRAY)
 
38
        punt("MO_TO_PLOT must be an array. See the manual.");
 
39
      
 
40
      errcod = ip_string("MO_TO_PLOT",&tmpstring,1,0);
 
41
      /* If signed integers are used - Fock ordering is used, convert indices back to Pitzer ... */
 
42
      if (tmpstring[0] == '+' || tmpstring[0] == '-') {
 
43
        fock_mos = init_int_array(num_mos_to_plot);
 
44
        for(mo=0;mo<num_mos_to_plot;mo++) {
 
45
          errcod = ip_data("MO_TO_PLOT","%d",&fock_mos[mo],1,mo);
 
46
          if (errcod != IPE_OK)
 
47
            punt("Problem reading elements of MO_TO_PLOT. See the manual.");
 
48
        }
 
49
        mos_to_plot = fock_to_pitzer(num_mos_to_plot,fock_mos);
 
50
        free(fock_mos);
 
51
      }
 
52
      /* ... else indices are already Pitzer indices */
 
53
      else {
 
54
        mos_to_plot = init_int_array(num_mos_to_plot);
 
55
        for(mo=0;mo<num_mos_to_plot;mo++) {
 
56
          errcod = ip_data("MO_TO_PLOT","%d",&mos_to_plot[mo],1,mo);
 
57
          if (mos_to_plot[mo] <= 0 || mos_to_plot[mo] > nmo)
 
58
            punt("One of the elements of MO_TO_PLOT out of range");
 
59
          mos_to_plot[mo]--;
 
60
          if (errcod != IPE_OK)
 
61
            punt("Problem reading elements of MO_TO_PLOT. See the manual.");
 
62
        }
 
63
      }
 
64
      free(tmpstring);
 
65
    }
 
66
    /* If MO_TO_PLOT is not specified just use HOMO and LUMO */
 
67
    else {
 
68
      num_mos_to_plot = 2;
 
69
      fock_mos = init_int_array(num_mos_to_plot);
 
70
      fock_mos[0] = -1;
 
71
      fock_mos[1] = 1;
 
72
      mos_to_plot = fock_to_pitzer(num_mos_to_plot,fock_mos);
 
73
      free(fock_mos);
 
74
    }
 
75
  }
 
76
  
 
77
  /*--------------------------------------------------------------------
 
78
    If GRID_ORIGIN is present read in user's specification for the grid
 
79
    else create a default grid using molecular dimensions
 
80
   --------------------------------------------------------------------*/
 
81
  if (ip_exist("GRID_ORIGIN",0)) {
 
82
    
 
83
    ip_count("GRID_ORIGIN",&i,0);
 
84
    if (i != 3)
 
85
      punt("GRID_ORIGIN must have 3 components");
 
86
    for (i=0;i<3;i++) {
 
87
      errcod = ip_data("GRID_ORIGIN","%lf",&grid_origin[i],1,i);
 
88
      if (errcod != IPE_OK)
 
89
        punt("Error in the definition of GRID_ORIGIN");
 
90
    }
 
91
    
 
92
    if (ip_exist("GRID_UNIT_X",0)) {
 
93
      ip_count("GRID_UNIT_X",&i,0);
 
94
      if (i != 3)
 
95
        punt("GRID_UNIT_X must have 3 components");
 
96
      for (i=0;i<3;i++) {
 
97
        errcod = ip_data("GRID_UNIT_X","%lf",&grid_unit_x[i],1,i);
 
98
        if (errcod != IPE_OK)
 
99
          punt("Problem in parsing GRID_UNIT_X");
 
100
      }
 
101
    }
 
102
    else
 
103
      punt("GRID_UNIT_X must be defined when GRID_ORIGIN is given");
 
104
    
 
105
    if (ip_exist("GRID_UNIT_Y",0)) {
 
106
      ip_count("GRID_UNIT_Y",&i,0);
 
107
      if (i != 3)
 
108
        punt("GRID_UNIT_Y must have 3 components");
 
109
      for (i=0;i<3;i++) {
 
110
        errcod = ip_data("GRID_UNIT_Y","%lf",&grid_unit_y[i],1,i);
 
111
        if (errcod != IPE_OK)
 
112
          punt("Problem in parsing GRID_UNIT_Y");
 
113
      }
 
114
    }
 
115
    else
 
116
      punt("GRID_UNIT_Y must be defined when GRID_ORIGIN is given");
 
117
 
 
118
    if (grid3d == 0) {
 
119
      if (ip_exist("GRID_XY0",0)) {
 
120
        ip_count("GRID_XY0",&i,0);
 
121
        if (i != 2)
 
122
          punt("GRID_XY0 must have 2 components");
 
123
        for (i=0;i<2;i++) {
 
124
          errcod = ip_data("GRID_XY0","%lf",&grid_xy0[i],1,i);
 
125
          if (errcod != IPE_OK)
 
126
            punt("Error in the definition of GRID_XY0");
 
127
        }
 
128
      }
 
129
      else
 
130
        punt("GRID_XY0 is not defined");
 
131
      
 
132
      if (ip_exist("GRID_XY1",0)) {
 
133
        ip_count("GRID_XY1",&i,0);
 
134
        if (i != 2)
 
135
          punt("GRID_XY1 must have 2 components");
 
136
        for (i=0;i<2;i++) {
 
137
          errcod = ip_data("GRID_XY1","%lf",&grid_xy1[i],1,i);
 
138
          if (errcod != IPE_OK)
 
139
            punt("Error in the definition of GRID_XY1");
 
140
          else if (grid_xy1[i] <= grid_xy0[i])
 
141
            punt("GRID_XY1 must point to the upper right corner of the grid");
 
142
        }
 
143
      }
 
144
      else
 
145
        punt("GRID_XY1 is not defined");
 
146
    }
 
147
    else {
 
148
      if (ip_exist("GRID_XYZ0",0)) {
 
149
        ip_count("GRID_XYZ0",&i,0);
 
150
        if (i != 3)
 
151
          punt("GRID_XYZ0 must have 3 components");
 
152
        for (i=0;i<3;i++) {
 
153
          errcod = ip_data("GRID_XYZ0","%lf",&grid_xyz0[i],1,i);
 
154
          if (errcod != IPE_OK)
 
155
            punt("Error in the definition of GRID_XYZ0");
 
156
        }
 
157
      }
 
158
      else
 
159
        punt("GRID_XYZ0 is not defined");
 
160
      
 
161
      if (ip_exist("GRID_XYZ1",0)) {
 
162
        ip_count("GRID_XYZ1",&i,0);
 
163
        if (i != 3)
 
164
          punt("GRID_XYZ1 must have 3 components");
 
165
        for (i=0;i<3;i++) {
 
166
          errcod = ip_data("GRID_XYZ1","%lf",&grid_xyz1[i],1,i);
 
167
          if (errcod != IPE_OK)
 
168
            punt("Error in the definition of GRID_XYZ1");
 
169
          else if (grid_xyz1[i] <= grid_xyz0[i])
 
170
            punt("GRID_XYZ1 must point to the upper right corner of the grid parallelepiped");
 
171
        }
 
172
      }
 
173
      else
 
174
        punt("GRID_XYZ1 is not defined");
 
175
      
 
176
    }
 
177
  }
 
178
  /* GRID_ORIGIN not specified -- use molecular dimensions in 3D case */
 
179
  else if (grid3d) {
 
180
 
 
181
    xmax = ymax = zmax = 0.0;
 
182
    xmin = ymin = zmin = 0.0;
 
183
    
 
184
    for(atom=0; atom<natom; atom++) {
 
185
      if (geom[atom][0] > xmax)
 
186
        xmax = geom[atom][0];
 
187
      if (geom[atom][0] < xmin)
 
188
        xmin = geom[atom][0];
 
189
      if (geom[atom][1] > ymax)
 
190
        ymax = geom[atom][1];
 
191
      if (geom[atom][1] < ymin)
 
192
        ymin = geom[atom][1];
 
193
      if (geom[atom][2] > zmax)
 
194
        zmax = geom[atom][2];
 
195
      if (geom[atom][2] < zmin)
 
196
        zmin = geom[atom][2];
 
197
    }
 
198
 
 
199
    grid_origin[0] = 0.0;
 
200
    grid_origin[1] = 0.0;
 
201
    grid_origin[2] = 0.0;
 
202
 
 
203
    grid_unit_x[0] = 1.0;
 
204
    grid_unit_x[1] = 0.0;
 
205
    grid_unit_x[2] = 0.0;
 
206
    grid_unit_y[0] = 0.0;
 
207
    grid_unit_y[1] = 1.0;
 
208
    grid_unit_y[2] = 0.0;
 
209
 
 
210
    grid_xyz0[0] = xmin - 2.0;
 
211
    grid_xyz0[1] = ymin - 2.0;
 
212
    grid_xyz0[2] = zmin - 2.0;
 
213
    grid_xyz1[0] = xmax + 2.0;
 
214
    grid_xyz1[1] = ymax + 2.0;
 
215
    grid_xyz1[2] = zmax + 2.0;
 
216
  }
 
217
  /* It's a 2D grid and no GRID_ORIGIN is given -- fail */
 
218
  else
 
219
    punt("GRID_ORIGIN must be specified for two-dimentional grids");
 
220
 
 
221
  i = 0; errcod = ip_data("NIX","%d",&i,0);
 
222
  nix = i - 1;
 
223
  if (i <= 1)
 
224
    punt("NIX must be greater than 1");
 
225
  i = 0; errcod = ip_data("NIY","%d",&i,0);
 
226
  niy = i - 1;
 
227
  if (i <= 1)
 
228
    punt("NIY must be greater than 1");
 
229
  if (grid3d) {
 
230
    i = 0; errcod = ip_data("NIZ","%d",&i,0);
 
231
    niz = i - 1;
 
232
    if (i <= 1)
 
233
      punt("NIZ must be greater than 1");
 
234
  }
 
235
 
 
236
  /* orthonormalize grid unit vectors */
 
237
  grid_unitvec();
 
238
  
 
239
  /* Which format to use for the 3D grid */
 
240
  grid_format = NULL;
 
241
  errcod = ip_string("GRID_FORMAT", &grid_format, 0);
 
242
  if (grid_format == NULL) {
 
243
    if (grid3d)
 
244
      grid_format = strdup("GAUSSCUBE");
 
245
    else
 
246
      grid_format = strdup("PLOTMTV");
 
247
  }
 
248
  else if (strcmp(grid_format,"GAUSSCUBE") &&
 
249
           strcmp(grid_format,"PLOTMTV") &&
 
250
           strcmp(grid_format,"MEGAPOVPLUS"))
 
251
    punt("Invalid value for GRID_FORMAT");
 
252
  
 
253
  /* Check if the grid output format is compatibe with the type of grid */
 
254
  if (!strcmp(grid_format,"PLOTMTV") && grid3d)
 
255
    punt("GRID_FORMAT=PLOTMTV can only be used for 2-d grids");
 
256
  if (strcmp(grid_format,"PLOTMTV") && !grid3d)
 
257
    punt("Only GRID_FORMAT=PLOTMTV can be used for 2-d grids");
 
258
  
 
259
  if (grid3d == 0) {
 
260
    errcod = ip_data("GRID_ZMIN","%lf",&grid_zmin,0);
 
261
    errcod = ip_data("GRID_ZMAX","%lf",&grid_zmax,0);
 
262
    if (grid_zmin >= grid_zmax)
 
263
      punt("GRID_ZMIN must be less than GRID_ZMAX");
 
264
    errcod = ip_data("EDGRAD_LOGSCALE","%d",&edgrad_logscale,0);
 
265
  }
 
266
}
 
267
 
 
268
}} //namespace psi::oeprop
 
269
namespace {
 
270
  using namespace psi::oeprop;
 
271
 
 
272
int* fock_to_pitzer(int num_mo, int* fock_mos)
 
273
{
 
274
  int mo, irrep, o, v;
 
275
  int nvirt, nocc, ndocc, nsocc;
 
276
  int nv, no, na;
 
277
  int current_highest, current_lowest;
 
278
  int *f2p_virt, *f2p_occ;
 
279
  int *mo_offset, *occ_offset, *virt_offset;
 
280
  int *pitzer_mos;
 
281
  double eval;
 
282
 
 
283
  ndocc = nvirt = nsocc = 0;
 
284
  for(irrep=0;irrep<nirreps;irrep++) {
 
285
    ndocc += clsdpi[irrep];
 
286
    nsocc += openpi[irrep];
 
287
    nvirt += orbspi[irrep] - clsdpi[irrep] - openpi[irrep];
 
288
  }
 
289
  nocc = ndocc + nsocc;
 
290
  f2p_occ = init_int_array(nocc);
 
291
  f2p_virt = init_int_array(nvirt);
 
292
  pitzer_mos = init_int_array(num_mo);
 
293
  virt_offset = init_int_array(nirreps);
 
294
  occ_offset = init_int_array(nirreps);
 
295
  mo_offset = init_int_array(nirreps);
 
296
 
 
297
  virt_offset[0] = clsdpi[0] + openpi[0];
 
298
  occ_offset[0] = clsdpi[0] + openpi[0] - 1;
 
299
  mo_offset[0] = 0;
 
300
  for(irrep=1;irrep<nirreps;irrep++) {
 
301
    na = orbspi[irrep];
 
302
    mo_offset[irrep] = mo_offset[irrep-1] + orbspi[irrep-1];
 
303
    no = clsdpi[irrep] + openpi[irrep];
 
304
    nv = na - no;
 
305
    if (no)
 
306
      occ_offset[irrep] = mo_offset[irrep] + no - 1;
 
307
    else
 
308
      occ_offset[irrep] = -1;
 
309
    if(nv)
 
310
      virt_offset[irrep] = mo_offset[irrep] + no;
 
311
    else
 
312
      virt_offset[irrep] = -1;
 
313
  }
 
314
 
 
315
  
 
316
  /*---------------------------------------------------------------
 
317
    construct Fock->Pitzer mapping for virtual orbitals (f2p_virt)
 
318
   ---------------------------------------------------------------*/
 
319
  for(v=0;v<nvirt;v++) {
 
320
    eval = 1E300;    /* to be greater than any eigenvalue to be encountered */
 
321
    for(irrep=0;irrep<nirreps;irrep++) {
 
322
      if (virt_offset[irrep] >= 0)
 
323
        if (eval > scf_evals[virt_offset[irrep]]) {
 
324
          current_lowest = irrep;
 
325
          eval = scf_evals[virt_offset[irrep]];
 
326
        }
 
327
    }
 
328
    f2p_virt[v] = virt_offset[current_lowest];
 
329
    virt_offset[current_lowest]++;
 
330
    /* If all virtual MOs in this irrep have been exhausted - set virt_offset for
 
331
       that irrep to -1 so that it's skipped on subsequent passes */
 
332
    if (current_lowest == nirreps-1) {
 
333
      if (virt_offset[current_lowest] == nmo)
 
334
        virt_offset[current_lowest] = -1;
 
335
    }
 
336
    else {
 
337
      if (virt_offset[current_lowest] == mo_offset[current_lowest+1])
 
338
        virt_offset[current_lowest] = -1;
 
339
    }
 
340
  }
 
341
 
 
342
  /*----------------------------------------------------------------
 
343
    construct Fock->Pitzer mapping for occupied orbitals (f2o_occ)
 
344
   ----------------------------------------------------------------*/
 
345
  for(o=0;o<nocc;o++) {
 
346
    eval = -1E300;    /* to be smaller than any eigenvalue to be encountered */
 
347
    for(irrep=0;irrep<nirreps;irrep++) {
 
348
      if (occ_offset[irrep] >= 0)
 
349
        if (eval < scf_evals[occ_offset[irrep]]) {
 
350
          current_highest = irrep;
 
351
          eval = scf_evals[occ_offset[irrep]];
 
352
        }
 
353
    }
 
354
    f2p_occ[o] = occ_offset[current_highest];
 
355
    occ_offset[current_highest]--;
 
356
    /* If all occupied MOs in this irrep have been exhausted - set occ_offset for
 
357
       that irrep to -1 so that it's skipped on subsequent passes */
 
358
    if (occ_offset[current_highest] == mo_offset[current_highest]-1)
 
359
      occ_offset[current_highest] = -1;
 
360
  }
 
361
 
 
362
  for(mo=0;mo<num_mo;mo++)
 
363
    if (fock_mos[mo] < 0) {
 
364
      o = -fock_mos[mo] - 1;
 
365
      if (o >= nocc)
 
366
        punt("MOS_TO_PLOT contains occupied indices outside of allowed range"); 
 
367
      else
 
368
        pitzer_mos[mo] = f2p_occ[o];
 
369
    }
 
370
    else {
 
371
      v = fock_mos[mo] - 1;
 
372
      if (v >= nvirt)
 
373
        punt("MOS_TO_PLOT contains unoccupied indices outside of allowed range"); 
 
374
      else
 
375
        pitzer_mos[mo] = f2p_virt[v];
 
376
    }
 
377
  
 
378
  free(virt_offset);
 
379
  free(occ_offset);
 
380
  free(mo_offset);
 
381
  free(f2p_occ);
 
382
  free(f2p_virt);
 
383
 
 
384
  return pitzer_mos;
 
385
}
 
386
 
 
387
} // namespace
 
388
 
 
389
namespace psi { namespace oeprop {
 
390
 
 
391
void grid_unitvec()
 
392
{
 
393
  int i;
 
394
  double sum1, sum2, dot;
 
395
  double step_x, step_y, step_z;
 
396
 
 
397
  /*-------------------------------------------------
 
398
    First orthonormalize grid_unit_x and grid_unit_y
 
399
    and compute grid_unit_z
 
400
   -------------------------------------------------*/
 
401
 
 
402
  /* Normalizing the grid_unit's */
 
403
  dot_arr(grid_unit_x,grid_unit_x,3,&sum1);
 
404
  dot_arr(grid_unit_y,grid_unit_y,3,&sum2);
 
405
  sum1 = sqrt(sum1); sum2 = sqrt(sum2);
 
406
  for(i=0;i<3;i++) {
 
407
    grid_unit_x[i] /= sum1;
 
408
    grid_unit_y[i] /= sum2;
 
409
  }
 
410
 
 
411
  /* Checking if vectors are parallel */
 
412
  dot_arr(grid_unit_x,grid_unit_y,3,&dot);
 
413
  if (1.0 - fabs(dot) < ADOTB_ORTHOGONAL) /* Vectors are parallel - abort */
 
414
    punt("Vectors GRID_UNIT_X and GRID_UNIT_Y are parallel");
 
415
  else if (fabs(dot) > ADOTB_ORTHOGONAL) {  /* Vectors are not orthogonal - orthonormalizing */
 
416
    for(i=0;i<3;i++)
 
417
      grid_unit_y[i] -= dot*grid_unit_x[i];
 
418
    dot_arr(grid_unit_y,grid_unit_y,3,&sum1);
 
419
    sum1 = sqrt(sum1);
 
420
    for(i=0;i<3;i++)
 
421
      grid_unit_y[i] /= sum1;
 
422
  }
 
423
 
 
424
  /* Get grid_unit_z as a vector product */
 
425
  if (grid3d) {
 
426
    grid_unit_z[0] = grid_unit_x[1]*grid_unit_y[2] - grid_unit_x[2]*grid_unit_y[1];
 
427
    grid_unit_z[1] = grid_unit_x[2]*grid_unit_y[0] - grid_unit_x[0]*grid_unit_y[2];
 
428
    grid_unit_z[2] = grid_unit_x[0]*grid_unit_y[1] - grid_unit_x[1]*grid_unit_y[0];
 
429
  }
 
430
 
 
431
  /* grid_origin will now contain the origin of the grid rectangle/box rather
 
432
     than the grid coordinate systems */
 
433
  if (grid3d)
 
434
    for(i=0;i<3;i++)
 
435
      grid_origin[i] += grid_xyz0[0]*grid_unit_x[i] + grid_xyz0[1]*grid_unit_y[i] + grid_xyz0[2]*grid_unit_z[i];
 
436
  else
 
437
    for(i=0;i<3;i++)
 
438
      grid_origin[i] += grid_xy0[0]*grid_unit_x[i] + grid_xy0[1]*grid_unit_y[i];
 
439
 
 
440
 
 
441
  /* Compute grid size and unit cell vectors */
 
442
  if (grid3d) {
 
443
    step_x = (grid_xyz1[0] - grid_xyz0[0])/nix;
 
444
    step_y = (grid_xyz1[1] - grid_xyz0[1])/niy;
 
445
    step_z = (grid_xyz1[2] - grid_xyz0[2])/niz;
 
446
    for(i=0;i<3;i++) {
 
447
      grid_step_x[i] = grid_unit_x[i]*step_x;
 
448
      grid_step_y[i] = grid_unit_y[i]*step_y;
 
449
      grid_step_z[i] = grid_unit_z[i]*step_z;
 
450
    }
 
451
  }
 
452
  else {
 
453
    step_x = (grid_xy1[0] - grid_xy0[0])/nix;
 
454
    step_y = (grid_xy1[1] - grid_xy0[1])/niy;
 
455
    for(i=0;i<3;i++) {
 
456
      grid_step_x[i] = grid_unit_x[i]*step_x;
 
457
      grid_step_y[i] = grid_unit_y[i]*step_y;
 
458
    }
 
459
  }
 
460
 
 
461
  return;
 
462
}
 
463
 
 
464
}} // namespace psi::oeprop