3
\brief Enter brief description of file here
8
#include <libciomr/libciomr.h>
10
#include <libiwl/iwl.h>
18
namespace psi { namespace ccresponse {
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.
28
** Notes on specific integrals (all produced by "cints --oeprop") used
31
** (1) Mu: Length-gauge electric dipole moment integrals = -r. These
32
** already include the electronic charge, and they are symmetric wrt
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
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).
53
void transpert(const char *pert)
55
int nao, nso, nmo, noei_ao;
58
double *scratch, **TMP, **X, **target;
60
double prefactor, anti, sign;
65
noei_ao = moinfo.noei_ao;
67
TMP = block_matrix(nao, nao);
68
X = block_matrix(nao, nao);
69
scratch = init_array(noei_ao);
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; }
77
for(alpha=0; alpha < 3; alpha++) {
79
target = block_matrix(nmo,nmo);
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; }
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; }
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; }
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];
104
C_DGEMM('n','t',nao,nso,nao,1,&(TMP[0][0]),nao,&(moinfo.usotao[0][0]),nao,
106
C_DGEMM('n','n',nso,nso,nao,1,&(moinfo.usotao[0][0]),nao,&(X[0][0]),nao,
109
C_DGEMM('n','n',nso,nmo,nso,1,&(TMP[0][0]),nao,&(moinfo.scf[0][0]),nmo,
111
C_DGEMM('t','n',nmo,nmo,nso,1,&(moinfo.scf[0][0]),nmo,&(X[0][0]),nao,
112
0,&(target[0][0]),nmo);
114
zero_arr(scratch,noei_ao);
123
}} // namespace psi::ccresponse