~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/ccresponse/transpert.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:
 
1
/*! \file
 
2
    \ingroup CCRESPONSE
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <cstring>
 
8
#include <libciomr/libciomr.h>
 
9
#include <libqt/qt.h>
 
10
#include <libiwl/iwl.h>
 
11
#include <psifiles.h>
 
12
#include "MOInfo.h"
 
13
#include "Params.h"
 
14
#include "Local.h"
 
15
#define EXTERN
 
16
#include "globals.h"
 
17
 
 
18
namespace psi { namespace ccresponse {
 
19
 
 
20
/* transpert(): Transform various one-electron property integrals from
 
21
** the AO to the MO basis.  In some cases, we must also add
 
22
** appropriate prefactors and signs.  The only argument is a
 
23
** character string indicating the type of integral we want: "Mu",
 
24
** "L", "L*", "P", or "P*".  The cints code produces only lower
 
25
** triangles, so we must unpack the integrals and keep up with
 
26
** symmetric vs. antisymmetric cases.
 
27
**
 
28
** Notes on specific integrals (all produced by "cints --oeprop") used
 
29
** here:
 
30
**
 
31
** (1) Mu: Length-gauge electric dipole moment integrals = -r.  These
 
32
** already include the electronic charge, and they are symmetric wrt
 
33
** index permutation.
 
34
 
 
35
** (2) L: Magnetic dipole integrals = -1/2 (r x p).  OK, actually, the
 
36
** input integrals are really just angular momentum integrals (r x p),
 
37
** but we multiply these by -0.5 to account for both the sign of the
 
38
** electronic charge and the definition of the magnetic dipole.  These
 
39
** integrals are antisymmetric wrt index permutation.  Use "L*" to use
 
40
** the complex conjugate of the operator (i.e., this multiplies by
 
41
** -1).
 
42
**
 
43
** (3) P: Velocity-gauge electric dipole moment integrals = -p.  OK,
 
44
** cints actually produces -del integrals, which already include the
 
45
** electronic charge, so we must multiply by -1 for the definition of
 
46
** the linear momentum operator.  They are antisymmetric wrt to index
 
47
** permutation. Use "P*" to use the complex conjugate of the operator
 
48
** (i.e., this multiplies by -1).
 
49
**
 
50
** -TDC, 11/05
 
51
*/
 
52
 
 
53
void transpert(const char *pert)
 
54
{
 
55
  int nao, nso, nmo, noei_ao;
 
56
  int alpha;
 
57
  int i, j, ij;
 
58
  double *scratch, **TMP, **X, **target;
 
59
  const char *name;
 
60
  double prefactor, anti, sign;
 
61
 
 
62
  nao = moinfo.nao;
 
63
  nso = moinfo.nso;
 
64
  nmo = moinfo.nmo;
 
65
  noei_ao = moinfo.noei_ao;
 
66
 
 
67
  TMP = block_matrix(nao, nao);
 
68
  X = block_matrix(nao, nao);
 
69
  scratch = init_array(noei_ao);
 
70
 
 
71
  if(!strcmp(pert,"Mu")) { prefactor = 1.0; anti = 1.0; sign = 1.0; }
 
72
  else if(!strcmp(pert, "L")) { prefactor = -0.5; anti = -1.0; sign = 1.0; }
 
73
  else if(!strcmp(pert, "L*")) { prefactor = -0.5; anti = -1.0; sign = -1.0; }
 
74
  else if(!strcmp(pert, "P")) { prefactor = -1.0; anti = -1.0; sign = 1.0; }
 
75
  else if(!strcmp(pert, "P*")) { prefactor = -1.0; anti = -1.0; sign = -1.0; }
 
76
 
 
77
  for(alpha=0; alpha < 3; alpha++) {
 
78
 
 
79
    target = block_matrix(nmo,nmo);
 
80
 
 
81
    if(!strcmp(pert,"Mu")) {
 
82
      if(alpha == 0) { name = PSIF_AO_MX; moinfo.MUX = target; }
 
83
      else if(alpha == 1) { name = PSIF_AO_MY; moinfo.MUY = target; }
 
84
      else if(alpha == 2) { name = PSIF_AO_MZ; moinfo.MUZ = target; }
 
85
    }
 
86
    else if(!strcmp(pert,"L") || !strcmp(pert, "L*")) {
 
87
      if(alpha == 0) { name = PSIF_AO_LX; moinfo.LX = target; }
 
88
      else if(alpha == 1) { name = PSIF_AO_LY; moinfo.LY = target; }
 
89
      else if(alpha == 2) { name = PSIF_AO_LZ; moinfo.LZ = target; }
 
90
    }
 
91
    else if(!strcmp(pert,"P") || !strcmp(pert, "P*")) {
 
92
      if(alpha == 0) { name = PSIF_AO_NablaX; moinfo.PX = target; }
 
93
      else if(alpha == 1) { name = PSIF_AO_NablaY; moinfo.PY = target; }
 
94
      else if(alpha == 2) { name = PSIF_AO_NablaZ; moinfo.PZ = target; }
 
95
    }
 
96
 
 
97
    iwl_rdone(PSIF_OEI, name, scratch, noei_ao, 0, 0, outfile);
 
98
    for(i=0,ij=0; i < nao; i++)
 
99
      for(j=0; j <= i; j++,ij++) {
 
100
        TMP[i][j] = prefactor * sign * scratch[ij];
 
101
        TMP[j][i] = anti * prefactor * sign * scratch[ij];
 
102
      }
 
103
 
 
104
    C_DGEMM('n','t',nao,nso,nao,1,&(TMP[0][0]),nao,&(moinfo.usotao[0][0]),nao,
 
105
            0,&(X[0][0]),nao);
 
106
    C_DGEMM('n','n',nso,nso,nao,1,&(moinfo.usotao[0][0]),nao,&(X[0][0]),nao,
 
107
            0,&(TMP[0][0]),nao);
 
108
 
 
109
    C_DGEMM('n','n',nso,nmo,nso,1,&(TMP[0][0]),nao,&(moinfo.scf[0][0]),nmo,
 
110
            0,&(X[0][0]),nao);
 
111
    C_DGEMM('t','n',nmo,nmo,nso,1,&(moinfo.scf[0][0]),nmo,&(X[0][0]),nao,
 
112
            0,&(target[0][0]),nmo);
 
113
 
 
114
    zero_arr(scratch,noei_ao);
 
115
 
 
116
  }
 
117
 
 
118
  free(scratch);
 
119
  free_block(TMP);
 
120
  free_block(X);
 
121
}
 
122
 
 
123
}} // namespace psi::ccresponse