1
/* DISP_FREQ_ENERGY_CART makes symmetry-adapted cartesian coordinates and displaces
2
along them to setup frequencies by energy points */
12
#include <libchkpt/chkpt.h>
16
#include <libciomr/libciomr.h>
18
#include <libipv1/ip_lib.h>
19
#include <physconst.h>
20
#include <libpsio/psio.h>
27
#include "cartesians.h"
28
#include "internals.h"
30
#include "bond_lengths.h"
32
extern int get_irrep_xyz( double **cartrep, int xyz);
33
extern int check_coordinates(int natom, double *coord, double *masses, double *Zvals,
34
int *ndisp_small, double ***disp_small);
36
int disp_freq_energy_cart(cartesians &carts)
38
int i,j,a,b,I,k,dim, ndisp_all, nsalc_all, natom, atom, xyz, cnt, loner, *ndisp;
39
int op, atom2, nirreps, natom_unique, irrep, diag_ind, atom_unique;
40
double *coord, energy, ***disp_orig, **displacements, *character, *disp_e;
41
double *v, *f, *q, *masses, **constraints, normval, dotval, **constraints_ortho, *Zvals;
42
int **ict, print=0,nconstraints, *ua2a;
43
int *nsalc_orig, *nsalc, *irrep_per_disp;
44
double ***geoms, **moi, *v1_A, *v1_B, *big_zeroes;
45
double *com, **X, *evals, tval, tval1, tval2, *disp_grad, **cartrep;
46
double ***salc_orig, ***salc;
50
coord = carts.get_coord();
51
energy = carts.get_energy();
52
natom = carts.get_natom();
53
masses = carts.get_mass();
54
Zvals = carts.get_Zvals();
55
nirreps = syminfo.nirreps;
57
chkpt_init(PSIO_OPEN_OLD);
58
natom_unique = chkpt_rd_num_unique_atom();
60
ua2a = chkpt_rd_ua2a();
61
cartrep = chkpt_rd_cartrep();
65
fprintf(outfile,"Character table from syminfo\n");
66
for (irrep=0; irrep<nirreps; ++irrep) {
67
for (op=0; op<nirreps; ++op)
68
fprintf(outfile,"%d",syminfo.ct[irrep][op]);
69
fprintf(outfile,"\n");
71
fprintf(outfile,"Cartesian rep from chkpt\n");
72
print_mat(cartrep,nirreps,9,outfile);
73
fprintf(outfile,"ICT table\n");
74
for (i=0; i<natom; ++i) {
75
for (op=0; op<nirreps; ++op)
76
fprintf(outfile," %d ", ict[op][i]);
77
fprintf(outfile,"\n");
81
/**** Make sure center of mass is at origin ****/
84
for (i=0; i<natom; ++i) {
85
for (xyz=0; xyz<3 ; ++xyz)
86
com[xyz] += masses[3*i+xyz] * coord[3*i+xyz];
89
for (xyz=0; xyz<3 ; ++xyz)
90
com[xyz] = com[xyz] / tval;
92
for (i=0; i<natom; ++i)
93
for (xyz=0; xyz<3 ; ++xyz)
94
coord[3*i+xyz] -= com[xyz];
97
fprintf(outfile,"\n\nCenter of mass\n");
98
for (xyz=0; xyz<3 ; ++xyz)
99
fprintf(outfile,"%15.10lf", com[xyz]);
100
fprintf(outfile,"\n\nNew Geometry wrt COM\n");
101
for (i=0; i<natom; ++i)
102
fprintf(outfile,"%15.10lf%15.10lf%15.10lf\n",coord[3*i],coord[3*i+1],coord[3*i+2]);
107
/**** Determine the rotational coordinates from MOI tensor ****/
108
moi = block_matrix(3,3);
109
for (i=0; i<natom; ++i) {
110
moi[0][0] += masses[3*i] * (SQR(coord[3*i+1]) + SQR(coord[3*i+2]));
111
moi[1][1] += masses[3*i] * (SQR(coord[3*i+0]) + SQR(coord[3*i+2]));
112
moi[2][2] += masses[3*i] * (SQR(coord[3*i+0]) + SQR(coord[3*i+1]));
113
moi[0][1] -= masses[3*i] * coord[3*i+0] * coord[3*i+1];
114
moi[0][2] -= masses[3*i] * coord[3*i+0] * coord[3*i+2];
115
moi[1][2] -= masses[3*i] * coord[3*i+1] * coord[3*i+2];
117
moi[1][0] = moi[0][1];
118
moi[2][0] = moi[0][2];
119
moi[2][1] = moi[1][2];
121
X = block_matrix(3,3);
122
evals = init_array(3);
123
sq_rsp(3,3,moi,evals,1,X,1.0E-14);
126
fprintf(outfile,"Eigenvectors (X) of MOI tensor\n");
127
print_mat(X,3,3,outfile);
130
/**** Build matrix with COM and rotational constraints as rows ****/
131
constraints = block_matrix(6, 3*natom);
132
/* COM constraints */
133
for (i=0; i<natom; ++i) {
134
constraints[0][3*i] = 1.0 * sqrt(masses[3*i]);
135
constraints[1][3*i+1] = 1.0 * sqrt(masses[3*i+1]);
136
constraints[2][3*i+2] = 1.0 * sqrt(masses[3*i+2]);
138
/* Rotational constraints */
139
for (i=0; i<natom; ++i) {
140
tval1 = (coord[3*i+0] * X[0][1]) + (coord[3*i+1] * X[1][1]) + (coord[3*i+2] * X[2][1]);
141
tval2 = (coord[3*i+0] * X[0][2]) + (coord[3*i+1] * X[1][2]) + (coord[3*i+2] * X[2][2]);
142
for (xyz = 0; xyz < 3; ++xyz)
143
constraints[3][3*i+xyz] = (tval1 * X[xyz][2] - tval2 * X[xyz][1]) * sqrt(masses[3*i]);
145
for (i=0; i<natom; ++i) {
146
tval1 = (coord[3*i+0] * X[0][2]) + (coord[3*i+1] * X[1][2]) + (coord[3*i+2] * X[2][2]);
147
tval2 = (coord[3*i+0] * X[0][0]) + (coord[3*i+1] * X[1][0]) + (coord[3*i+2] * X[2][0]);
148
for (xyz = 0; xyz < 3; ++xyz)
149
constraints[4][3*i+xyz] = (tval1 * X[xyz][0] - tval2 * X[xyz][2]) * sqrt(masses[3*i]);
151
for (i=0; i<natom; ++i) {
152
tval1 = (coord[3*i+0] * X[0][0]) + (coord[3*i+1] * X[1][0]) + (coord[3*i+2] * X[2][0]);
153
tval2 = (coord[3*i+0] * X[0][1]) + (coord[3*i+1] * X[1][1]) + (coord[3*i+2] * X[2][1]);
154
for (xyz = 0; xyz < 3; ++xyz)
155
constraints[5][3*i+xyz] = (tval1 * X[xyz][1] - tval2 * X[xyz][0]) * sqrt(masses[3*i]);
158
fprintf(outfile,"Raw COM/Rotational Constraints\n");
159
print_mat(constraints,6,3*natom,outfile);
163
/* Remove NULL constraint (if present) and normalize rest of them */
165
for (i = 0; i < 6; ++i) {
166
dot_arr(constraints[i], constraints[i], 3*natom, &normval);
167
if (normval > 1.0E-10) {
169
for (j=0; j<3*natom; ++j)
170
constraints[cnt][j] = constraints[i][j] / sqrt(normval) ;
173
nconstraints = cnt+1;
176
fprintf(outfile,"Normalized constraints\n");
177
print_mat(constraints,nconstraints,3*natom,outfile);
181
/* Orthogonalize rotations and translations exactly-is this ever necessary?*/
182
constraints_ortho = block_matrix(nconstraints,3*natom);
184
for (i=0; i<nconstraints; ++i)
185
cnt += schmidt_add(constraints_ortho, cnt, 3*natom, constraints[i]);
186
for (i=0; i<nconstraints; ++i)
187
for (j=0; j<3*natom; ++j)
188
constraints[i][j] = constraints_ortho[i][j];
189
free_block(constraints_ortho);
192
fprintf(outfile,"Orthogonal constraints\n");
193
print_mat(constraints,nconstraints,3*natom,outfile);
197
/**** Form symmetry-adapted cartesian vectors ****/
198
salc_orig = (double ***) malloc(nirreps*sizeof(double **));
199
nsalc_orig = init_int_array(nirreps);
201
/* loop over irrep of coordinates */
202
for (irrep=0; irrep < nirreps; ++irrep) {
203
salc_orig[irrep] = block_matrix(3*natom,3*natom);
206
/* loop over unique atoms and xyz coordinates */
207
for (atom_unique=0; atom_unique < natom_unique; ++atom_unique) {
208
atom = ua2a[atom_unique];
209
/* determine if atom is loner, having no symmetry-related partners */
211
for (op=0; op < nirreps; ++op) {
212
if (atom != ict[op][atom] - 1)
216
for (xyz=0; xyz<3; ++xyz) {
218
/* if a loner, don't try to adapt */
220
if ( get_irrep_xyz( cartrep, xyz ) == irrep) {
221
salc_orig[irrep][cnt][3*atom+xyz] = 1.0;
225
diag_ind = xyz*3 + xyz;
226
for (op=0; op < nirreps; ++op) {
227
atom2 = ict[op][atom] - 1;
228
for (i=0; i<3*natom; ++i)
229
salc_orig[irrep][cnt][3*atom2+xyz] += cartrep[op][diag_ind] / syminfo.ct[irrep][op];
234
nsalc_orig[irrep] = cnt+1;
237
for (irrep=0; irrep < nirreps; ++irrep) {
238
fprintf(outfile,"Cartesian SALCs for Irrep %d\n", irrep);
239
print_mat(salc_orig[irrep], nsalc_orig[irrep], 3*natom, outfile);
243
/**** Schmidt orthogonalize coordinates to remove COM/ROT constraints ****/
244
salc = (double ***) malloc(nirreps*sizeof(double **));
245
nsalc = init_int_array(nirreps);
246
v = init_array(3*natom);
248
for (irrep=0; irrep< nirreps; ++irrep) {
249
salc[irrep] = block_matrix(nsalc_orig[irrep],3*natom);
252
for (i=0; i<nsalc_orig[irrep]; ++i) {
255
for (I=0; I<3*natom; I++)
256
v[I] = salc_orig[irrep][i][I];
258
for (j=0; j<nconstraints; j++) {
259
dot_arr(salc_orig[irrep][i], constraints[j], 3*natom, &dotval) ;
260
for (I=0; I<3*natom; I++)
261
v[I] -= dotval * constraints[j][I] ;
264
dot_arr(v, v, 3*natom, &normval);
265
if (normval > 1.0E-10) {
266
for (j=0; j<3*natom; ++j)
267
v[j] = v[j] / sqrt(normval) ;
269
cnt += schmidt_add(salc[irrep], cnt, 3*natom, v);
275
for (irrep=0; irrep<nirreps; ++irrep) {
276
free_block(salc_orig[irrep]);
278
free_block(constraints);
282
if (optinfo.freq_irrep != -1) {
283
for (irrep=0; irrep < nirreps; ++irrep)
284
if (irrep != (optinfo.freq_irrep - 1) ) {
285
nsalc_all -= nsalc[irrep];
291
for (irrep=0; irrep < nirreps; ++irrep)
293
fprintf(outfile,"Projected Cartesian Displacements, irrep %d\n", irrep);
294
print_mat(salc[irrep], nsalc[irrep], 3*natom, outfile);
298
check_coordinates(natom,coord,masses,Zvals,nsalc,salc);
300
/* determine number of displacements */
301
/* 1/h^2 * (f(1,0) + f(-1,0) + 2f(0,0) ) */
303
ndisp = init_int_array(nirreps);
304
for (irrep=0; irrep<nirreps; ++irrep) {
306
/* diagonal displacements */
308
if (optinfo.points == 3)
309
ndisp[irrep] = 2 * nsalc[irrep];
310
else if (optinfo.points == 5)
311
ndisp[irrep] = 4 * nsalc[irrep];
314
if (optinfo.points == 3)
315
ndisp[irrep] = nsalc[irrep];
316
else if (optinfo.points == 5)
317
ndisp[irrep] = 2* nsalc[irrep];
320
/* off-diagonal displacements */
321
if (optinfo.points == 3)
322
ndisp[irrep] += 4 * nsalc[irrep] * (nsalc[irrep] - 1) / 2;
323
else if (optinfo.points == 5)
324
ndisp[irrep] += 16 * nsalc[irrep] * (nsalc[irrep] - 1) / 2;
326
ndisp_all += ndisp[irrep];
329
for (irrep=0; irrep<nirreps;++irrep)
330
fprintf(outfile,"ndisp[%d]: %d\n",irrep,ndisp[irrep]);
333
/*** construct displaced geometries ***/
334
geoms = (double ***) malloc(nirreps*sizeof(double **));
335
for (irrep=0; irrep<nirreps; ++irrep) {
337
geoms[irrep] = block_matrix(ndisp[irrep],3*natom);
338
for (i=0; i<ndisp[irrep]; ++i)
339
for (j=0; j<3*natom; ++j)
340
geoms[irrep][i][j] = coord[j];
342
/* construct diagonal displacements */
345
for (i=0; i<nsalc[irrep]; ++i) {
346
if (optinfo.points == 3) {
347
for (j=0; j < 3*natom; ++j) {
348
geoms[irrep][cnt+0][j] -= optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
349
geoms[irrep][cnt+1][j] += optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
353
else if (optinfo.points == 5) {
354
for (j=0; j < 3*natom; ++j) {
355
geoms[irrep][cnt+0][j] -= 2.0 * optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
356
geoms[irrep][cnt+1][j] -= 1.0 * optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
357
geoms[irrep][cnt+2][j] += 1.0 * optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
358
geoms[irrep][cnt+3][j] += 2.0 * optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
364
else { /* asymmetric diagonal displacements */
366
for (i=0; i<nsalc[irrep]; ++i) {
367
if (optinfo.points == 3) {
368
for (j=0; j < 3*natom; ++j)
369
geoms[irrep][cnt][j] -= optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
372
else if (optinfo.points == 5) {
373
for (j=0; j < 3*natom; ++j) {
374
geoms[irrep][cnt+0][j] -= 2.0 * optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
375
geoms[irrep][cnt+1][j] -= 1.0 * optinfo.disp_size * salc[irrep][i][j] / sqrt(masses[j]);
381
/* add off-diagonal displacements */
383
for (i=0; i<nsalc[irrep]; ++i) {
384
for (j=i+1; j<nsalc[irrep]; ++j) {
385
if (optinfo.points == 3) {
386
for (k=0; k < 3*natom; ++k) {
387
geoms[irrep][cnt+0][k] += ( + salc[irrep][i][k] + salc[irrep][j][k] )
388
* optinfo.disp_size / sqrt(masses[k]);
389
geoms[irrep][cnt+1][k] += ( + salc[irrep][i][k] - salc[irrep][j][k] )
390
* optinfo.disp_size / sqrt(masses[k]);
391
geoms[irrep][cnt+2][k] += ( - salc[irrep][i][k] + salc[irrep][j][k] )
392
* optinfo.disp_size / sqrt(masses[k]);
393
geoms[irrep][cnt+3][k] += ( - salc[irrep][i][k] - salc[irrep][j][k] )
394
* optinfo.disp_size / sqrt(masses[k]);
398
else if (optinfo.points == 5) {
399
for (k=0; k < 3*natom; ++k) {
400
geoms[irrep][cnt+0][k] += ( - 2.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
401
* optinfo.disp_size / sqrt(masses[k]);
402
geoms[irrep][cnt+1][k] += ( - 1.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
403
* optinfo.disp_size / sqrt(masses[k]);
404
geoms[irrep][cnt+2][k] += ( + 1.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
405
* optinfo.disp_size / sqrt(masses[k]);
406
geoms[irrep][cnt+3][k] += ( + 2.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
407
* optinfo.disp_size / sqrt(masses[k]);
408
geoms[irrep][cnt+4][k] += ( - 2.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
409
* optinfo.disp_size / sqrt(masses[k]);
410
geoms[irrep][cnt+5][k] += ( - 1.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
411
* optinfo.disp_size / sqrt(masses[k]);
412
geoms[irrep][cnt+6][k] += ( + 1.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
413
* optinfo.disp_size / sqrt(masses[k]);
414
geoms[irrep][cnt+7][k] += ( + 2.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
415
* optinfo.disp_size / sqrt(masses[k]);
416
geoms[irrep][cnt+8][k] += ( - 2.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
417
* optinfo.disp_size / sqrt(masses[k]);
418
geoms[irrep][cnt+9][k] += ( - 1.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
419
* optinfo.disp_size / sqrt(masses[k]);
420
geoms[irrep][cnt+10][k] += ( + 1.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
421
* optinfo.disp_size / sqrt(masses[k]);
422
geoms[irrep][cnt+11][k] += ( + 2.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
423
* optinfo.disp_size / sqrt(masses[k]);
424
geoms[irrep][cnt+12][k] += ( - 2.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
425
* optinfo.disp_size / sqrt(masses[k]);
426
geoms[irrep][cnt+13][k] += ( - 1.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
427
* optinfo.disp_size / sqrt(masses[k]);
428
geoms[irrep][cnt+14][k] += ( + 1.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
429
* optinfo.disp_size / sqrt(masses[k]);
430
geoms[irrep][cnt+15][k] += ( + 2.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
431
* optinfo.disp_size / sqrt(masses[k]);
440
/**** write out info to PSIF_OPTKING, geoms and coordinates in irrep order ****/
443
big_zeroes = init_array(ndisp_all*3*natom);
444
psio_write_entry(PSIF_OPTKING, "OPT: Displaced geometries", (char *) &(big_zeroes[0]),
445
ndisp_all*3*natom*sizeof(double));
446
psio_write_entry(PSIF_OPTKING, "OPT: Adapted cartesians", (char *) &(big_zeroes[0]),
447
ndisp_all*3*natom*sizeof(double));
448
psio_write_entry(PSIF_OPTKING, "OPT: Displaced gradients", (char *) &(big_zeroes[0]),
449
ndisp_all*3*natom*sizeof(double));
453
for (irrep=0; irrep<nirreps; ++irrep) {
455
psio_write(PSIF_OPTKING, "OPT: Displaced geometries",
456
(char *) &(geoms[irrep][0][0]), ndisp[irrep]*3*natom*sizeof(double), next, &next);
457
fprintf(outfile,"Displaced Geometries of irrep %d.\n",irrep);
458
print_mat(geoms[irrep],ndisp[irrep],3*natom,outfile);
460
free_block(geoms[irrep]);
465
for (irrep=0; irrep<nirreps; ++irrep) {
467
psio_write(PSIF_OPTKING, "OPT: Adapted cartesians", (char *) &(salc[irrep][0][0]),
468
nsalc[irrep]*3*natom*sizeof(double), next, &next);
470
free_block(salc[irrep]);
474
/* make list of irreps for each disp used to make prefixes */
475
irrep_per_disp = (int *) malloc(ndisp_all*sizeof(int));
477
for (irrep=0; irrep<nirreps; ++irrep) {
478
for (i=0; i<ndisp[irrep]; ++i)
479
irrep_per_disp[++cnt] = irrep;
481
fprintf(outfile,"Irreps per disp: ");
482
for (i=0; i<ndisp_all; ++i)
483
fprintf(outfile,"%d ",irrep_per_disp[i]);
484
fprintf(outfile,"\n");
486
/**** Write parameters to PSIO file ****/
487
optinfo.disp_num = 0;
488
disp_e = init_array(ndisp_all);
489
optinfo.micro_iteration = 0;
502
psio_write_entry(PSIF_OPTKING, "OPT: Displaced energies",
503
(char *) &(disp_e[0]), ndisp_all*sizeof(double));
504
psio_write_entry(PSIF_OPTKING, "OPT: Reference geometry",
505
(char *) &(coord[0]), 3*natom*sizeof(double));
506
psio_write_entry(PSIF_OPTKING, "OPT: Reference energy",
507
(char *) &(energy), sizeof(double));
508
psio_write_entry(PSIF_OPTKING, "OPT: Current disp_num",
509
(char *) &(optinfo.disp_num),sizeof(int));
510
psio_write_entry(PSIF_OPTKING, "OPT: Num. of disp.",
511
(char *) &(ndisp_all), sizeof(int));
512
psio_write_entry(PSIF_OPTKING, "OPT: Disp. per irrep",
513
(char *) &(ndisp[0]), nirreps*sizeof(int));
514
psio_write_entry(PSIF_OPTKING, "OPT: Irrep per disp",
515
(char *) &(irrep_per_disp[0]), ndisp_all*sizeof(int));
516
psio_write_entry(PSIF_OPTKING, "OPT: Num. of coord.",
517
(char *) &(nsalc_all), sizeof(int));
518
psio_write_entry(PSIF_OPTKING, "OPT: Coord. per irrep",
519
(char *) &(nsalc[0]), nirreps*sizeof(int));
520
psio_write_entry(PSIF_OPTKING, "Micro_iteration",
521
(char *) &(optinfo.micro_iteration),sizeof(int));
524
for (irrep=0; irrep<nirreps; ++irrep)
525
fprintf(outfile,"Number of %s displaced geometries is %d.\n",syminfo.irrep_lbls[irrep],ndisp[irrep]);
526
fprintf(outfile,"Total number of displaced geometries is %d.\n",ndisp_all);
added
removed
src/bin/optking/disp_freq_energy_cart.cc
