3
#include <libipv1/ip_lib.h>
4
#include <libciomr/libciomr.h>
5
#include <libchkpt/chkpt.h>
8
#include "globaldefs.h"
11
#define MO_HESS_MIN 1.0E-2
17
** This function calculates the step in theta space for the orbitals
18
** given the orbital gradient and an approximate orbital Hessian
23
void calc_orb_step(int npairs, double *grad, double *hess_diag, double *theta)
29
for (pair=0; pair<npairs; pair++) {
31
denom = hess_diag[pair];
33
fprintf(outfile, "Warning: MO Hessian denominator negative\n");
36
if (denom < MO_HESS_MIN) {
37
fprintf(outfile, "Warning: MO Hessian denominator too small\n");
40
theta[pair] = - numer / denom;
47
** calc_orb_step_full()
49
** This function calculates the step in theta space for the orbitals
50
** given the orbital gradient and a square matrix orbital Hessian
55
void calc_orb_step_full(int npairs, double *grad, double **hess, double *theta)
58
double **hess_copy; /* for testing! */
64
double hess_det = 1.0;
68
hess_copy = block_matrix(npairs, npairs);
69
indx = init_int_array(npairs);
71
for (i=0; i<npairs; i++) {
72
for (j=0; j<npairs; j++) {
73
hess_copy[i][j] = hess[i][j];
77
ludcmp(hess_copy,npairs,indx,&hess_det);
78
for (j=0;j<npairs;j++){
79
hess_det *= hess_copy[j][j];
81
fprintf(outfile,"The determinant of the hessian is %8.3E\n",hess_det);
85
if the orbital Hessian is not positive definite, we may have some
86
trouble converging the orbitals. Guarantee it's positive definite
89
if (Params.level_shift) {
90
while (hess_det < Params.determ_min) {
91
fprintf(outfile,"Level shifting the hessian by %8.3E\n",Params.shift);
92
for (i=0;i<npairs;i++) {
93
hess[i][i] += Params.shift;
95
for (i=0;i<npairs;i++) {
96
for (j=0;j<npairs;j++) {
97
hess_copy[i][j]=hess[i][j];
100
ludcmp(hess_copy,npairs,indx,&hess_det);
101
for (j=0;j<npairs;j++){
102
hess_det *= hess_copy[j][j];
104
fprintf(outfile,"The determinant of the hessian is %8.3E\n",hess_det);
106
fprintf(outfile,"Determinant of the hessian is greater than %8.3E\n",
111
/* let's re-copy hess into hess_copy because we ludcmp'd hess_copy */
112
for (i=0;i<npairs;i++) {
113
for (j=0;j<npairs;j++) {
114
hess_copy[i][j]=hess[i][j];
118
if (!Params.invert_hessian) { /* solve H delta = - g */
119
fprintf(outfile,"Solving system of linear equations for orbital step...");
120
BVector = init_array(npairs);
121
pivots = init_int_array(0.5 * npairs * (npairs - 1));
122
for(i=0;i<npairs;i++){
123
BVector[i] = -grad[i];
126
solved = C_DGESV(npairs,1,&(hess_copy[0][0]),npairs,pivots,
129
fprintf(outfile,"equations solved!\n");
130
for(i=0;i<npairs;i++) {
131
theta[i] = BVector[i];
135
fprintf(outfile,"FAILED TO SOLVE FOR THETA VALUES\n");
136
fprintf(outfile,"DGESV returned error %5d \n",solved);
141
} /* end solution of linear equations H delta = -g */
143
else { /* direct inversion of orbital Hessian */
144
fprintf(outfile,"Attempting to directly invert the Hessian matrix\n");
145
hess_inv = block_matrix(npairs,npairs);
147
/* note: this will destroy hessian matrix; don't use it again later! */
148
invert_matrix(hess_copy,hess_inv,npairs,outfile);
151
mmult(hess_inv,0,hess,0,hess_copy,0,npairs,npairs,npairs,0);
152
fprintf(outfile, "Hessian * Hessian inverse = \n");
153
print_mat(hess_copy,npairs,npairs,outfile);
154
fprintf(outfile, "\n");
156
/* step = - B^{-1} * g */
157
zero_arr(theta,npairs);
158
/* the line below seems to have trouble unless I take out the -1
159
which should be there, and even then it's not really working */
161
C_DGEMV('n',npairs,npairs,-1.0,hess_inv[0],npairs,grad,1,0.0,theta,1);
164
for (i=0; i<npairs; i++) {
166
for (j=0; j<npairs; j++) {
167
tval += hess_inv[i][j] * grad[j];
171
free_block(hess_inv);
172
} /* end direct inversion of Hessian */
174
/* make sure the step is not too big */
176
for (i=0; i<npairs; i++) {
178
if (fabs(tval) > biggest_step) biggest_step = fabs(tval);
180
fprintf(outfile,"\nLargest step in theta space is %12.6lf \n", biggest_step);
181
if (biggest_step > Params.step_max) {
182
fprintf(outfile, "Scaling the step\n");
183
for (i=0;i<npairs;i++) {
184
theta[i] = theta[i] * Params.step_max / biggest_step;
188
free_block(hess_copy);
194
** calc_orb_step_bfgs()
196
** This function calculates the step in theta space for the orbitals
197
** given the orbital gradient and a square matrix orbital Hessian INVERSE.
198
** With the inverse already available, this is very straightforward.
203
void calc_orb_step_bfgs(int npairs, double *grad, double **hess, double *theta)
207
double tval, biggest_step;
209
for (i=0; i<npairs; i++) {
211
for (j=0; j<npairs; j++) {
212
tval += hess[i][j] * grad[j];
217
/* make sure the step is not too big */
219
for (i=0; i<npairs; i++) {
221
if (fabs(tval) > biggest_step) biggest_step = fabs(tval);
223
fprintf(outfile,"\nLargest step in theta space is %12.6lf \n", biggest_step);
224
if (biggest_step > Params.step_max) {
225
fprintf(outfile, "Largest allowed step %12.6lf --- scaling the step\n",
227
for (i=0;i<npairs;i++) {
228
theta[i] = theta[i] * Params.step_max / biggest_step;
238
** This function prints out the information for a given orbital iteration
240
int print_step(int npairs, int steptype)
243
char sumfile_name[] = "file14.dat";
244
int i, entries, iter, *nind;
245
double *rmsgrad, *scaled_rmsgrad, *energies, energy;
248
/* open ascii file, get number of entries already in it */
250
ffile_noexit(&sumfile,sumfile_name,2);
251
if (sumfile == NULL) { /* the file doesn't exist yet */
253
if (Params.print_lvl)
254
fprintf(outfile, "\nPreparing new file %s\n", sumfile_name);
257
if (fscanf(sumfile, "%d", &entries) != 1) {
258
fprintf(outfile,"(print_step): Trouble reading num entries in file %s\n",
265
rmsgrad = init_array(entries+1);
266
scaled_rmsgrad = init_array(entries+1);
267
energies= init_array(entries+1);
268
nind = init_int_array(entries+1);
269
comments = (char **) malloc ((entries+1) * sizeof (char *));
270
for (i=0; i<entries+1; i++) {
271
comments[i] = (char *) malloc (MAX_COMMENT * sizeof(char));
274
for (i=0; i<entries; i++) {
275
fscanf(sumfile, "%d %d %lf %lf %lf %s", &iter, &(nind[i]),
276
&(scaled_rmsgrad[i]), &(rmsgrad[i]), &(energies[i]), comments[i]);
279
chkpt_init(PSIO_OPEN_OLD);
280
if (chkpt_exist("State averaged energy")) {
281
energy = chkpt_rd_e_labeled("State averged energy");
284
energy = chkpt_rd_etot();
287
scaled_rmsgrad[entries] = CalcInfo.scaled_mo_grad_rms;
288
rmsgrad[entries] = CalcInfo.mo_grad_rms;
289
energies[entries] = energy;
290
nind[entries] = npairs;
293
strcpy(comments[entries], "CONV");
294
else if (steptype == 1)
295
strcpy(comments[entries], "NR");
296
else if (steptype == 2)
297
strcpy(comments[entries], "DIIS");
299
fprintf(outfile, "(print_step): Unrecognized steptype %d\n", steptype);
300
strcpy(comments[entries], "?");
303
if (entries) fclose(sumfile);
305
/* now open file for writing, write out old info plus new */
306
ffile_noexit(&sumfile,"file14.dat",0);
307
if (sumfile == NULL) {
308
fprintf(outfile, "(print_step): Unable to open file %s\n", sumfile_name);
312
fprintf(sumfile, "%5d\n", entries);
313
for (i=0; i<entries; i++) {
314
fprintf(sumfile, "%5d %5d %14.9lf %14.9lf %20.12lf %9s\n", i+1, nind[i],
315
scaled_rmsgrad[i], rmsgrad[i], energies[i], comments[i]);
320
free(scaled_rmsgrad);
324
for (i=0; i<entries; i++)
added
removed
src/bin/detcas/step.c
