3
\brief Enter brief description of file here
8
#include <libciomr/libciomr.h>
10
#include <libiwl/iwl.h>
16
namespace psi { namespace ccsort {
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.
26
** Notes on specific integrals (all produced by "cints --oeprop") used
29
** (1) Mu: Length-gauge electric dipole moment integrals = -r. These
30
** already include the electronic charge, and they are symmetric wrt
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
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).
51
void transpert(const char *pert)
53
int nao, nso, nmo, noei_ao;
56
double *scratch, **TMP, **X, **target;
58
double prefactor, anti, sign;
63
noei_ao = nao * (nao+1)/2;
65
TMP = block_matrix(nao, nao);
66
X = block_matrix(nao, nao);
67
scratch = init_array(noei_ao);
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; }
75
for(alpha=0; alpha < 3; alpha++) {
77
target = block_matrix(nmo,nmo);
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; }
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; }
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; }
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];
102
C_DGEMM('n','t',nao,nso,nao,1,&(TMP[0][0]),nao,&(moinfo.usotao[0][0]),nao,
104
C_DGEMM('n','n',nso,nso,nao,1,&(moinfo.usotao[0][0]),nao,&(X[0][0]),nao,
107
C_DGEMM('n','n',nso,nmo,nso,1,&(TMP[0][0]),nao,&(moinfo.scf[0][0]),nmo,
109
C_DGEMM('t','n',nmo,nmo,nso,1,&(moinfo.scf[0][0]),nmo,&(X[0][0]),nao,
110
0,&(target[0][0]),nmo);
112
zero_arr(scratch,noei_ao);
121
}} // namespace psi::ccsort