/
build_Z_RHF.cc
103 lines (90 loc) · 3.13 KB
/
build_Z_RHF.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
/*
* @BEGIN LICENSE
*
* Psi4: an open-source quantum chemistry software package
*
* Copyright (c) 2007-2018 The Psi4 Developers.
*
* The copyrights for code used from other parties are included in
* the corresponding files.
*
* This file is part of Psi4.
*
* Psi4 is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, version 3.
*
* Psi4 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License along
* with Psi4; if not, write to the Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
* @END LICENSE
*/
/*! \file
\ingroup CCDENSITY
\brief Enter brief description of file here
*/
#include <cstdlib>
#include "psi4/libdpd/dpd.h"
#include "psi4/libqt/qt.h"
#include "psi4/libciomr/libciomr.h"
#include <cmath>
#include "MOInfo.h"
#include "Params.h"
#include "Frozen.h"
#define EXTERN
#include "globals.h"
namespace psi {
namespace ccdensity {
/* build_Z_RHF(): Solve the orbital Z-vector equations for RHF refs:
**
** sum E,M A(AI,EM) D(orb)(E,M) = -X(A,I)
**
** where A(AI,EM) is the orbital Hessian computed in build_A(), X(A,I)
** is the orbital rotation gradient computed in build_X(), and
** D(orb)(E,M) is the final Z-vector we want.
**
*/
void build_Z_RHF() {
dpdbuf4 A;
dpdfile2 X1, D;
double *X;
int h, nirreps, a, i, count;
nirreps = moinfo.nirreps;
/* Grab only irrep 0 of the orbital Hessian */
global_dpd_->buf4_init(&A, PSIF_CC_MISC, 0, 11, 11, 11, 11, 0, "A(EM,AI)");
global_dpd_->buf4_mat_irrep_init(&A, 0);
global_dpd_->buf4_mat_irrep_rd(&A, 0);
/* Place all the elements of the orbital rotation gradient, X into a
linear array, Z */
global_dpd_->file2_init(&X1, PSIF_CC_OEI, 0, 1, 0, "XAI");
global_dpd_->file2_mat_init(&X1);
global_dpd_->file2_mat_rd(&X1);
X = init_array(A.params->rowtot[0]);
for (h = 0, count = 0; h < nirreps; h++)
for (a = 0; a < X1.params->rowtot[h]; a++)
for (i = 0; i < X1.params->coltot[h]; i++) X[count++] = -X1.matrix[h][a][i];
global_dpd_->file2_mat_close(&X1);
global_dpd_->file2_close(&X1);
/* Trying out Matt's Pople code --- way to go, Matt! */
pople(A.matrix[0], X, A.params->rowtot[0], 1, 1e-12, "outfile", 0);
/* Build the orbital component of Dai */
global_dpd_->file2_init(&D, PSIF_CC_OEI, 0, 1, 0, "D(orb)(A,I)");
global_dpd_->file2_mat_init(&D);
for (h = 0, count = 0; h < nirreps; h++)
for (a = 0; a < D.params->rowtot[h]; a++)
for (i = 0; i < D.params->coltot[h]; i++) D.matrix[h][a][i] = X[count++];
global_dpd_->file2_mat_wrt(&D);
global_dpd_->file2_mat_close(&D);
global_dpd_->file2_close(&D);
free(X);
global_dpd_->buf4_mat_irrep_close(&A, 0);
global_dpd_->buf4_close(&A);
}
} // namespace ccdensity
} // namespace psi