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

« back to all changes in this revision

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