1
/** FREQ_GRAD_CART computes frequencies from gradients and cartesian disps */
11
#include <libchkpt/chkpt.h>
15
#include <libciomr/libciomr.h>
17
#include <libipv1/ip_lib.h>
18
#include <physconst.h>
19
#include <libpsio/psio.h>
26
#include "cartesians.h"
27
#include "internals.h"
29
#include "bond_lengths.h"
32
extern double **compute_B(internals &simples, salc_set &salcs);
33
extern double *compute_q(internals &simples, salc_set &symm);
34
extern double **compute_G(double **B, int num_intcos, cartesians &carts);
35
extern int get_irrep_xyz(double **cartrep, int xyz);
36
void sort_evals_all(int nsalc_all, double *evals_all, int *evals_all_irrep);
39
/**** FREQ_GRAD_CART compute frequencies from gradients for cartesian
42
void freq_grad_cart(cartesians &carts) {
43
int i,j,k,l,a,b, ii, cnt, dim, natom, xyz, cnt_eval = -1, *evals_all_irrep,op_disp;
44
int xyzA, xyzB, atomA, atomB, row, col,h,nirreps,xyz_irr,I,cnt_all,match,sign;
45
int nsalcs,ncoord,start_disp, start_salc, atom, atom2, op, loner, natom_unique;
46
double **B, **G, **G_inv, *masses, **u, *geom, *forces, **force_constants;
47
double energy, *energies, **displacements, cm_convert, k_convert;
48
double *f, *f_q, *temp_arr, *temp_arr2, tval, **geom2D, **B_inv;
49
double **evects, *evals, **FG, tmp, **force_constants_x, **tmp_mat_inv;
50
double *micro_e, **micro_geom, **disp_grad, *grad, **grads_adapted, **tmp_mat;
51
double *evals_all, **cartrep, **disp_grad_all;
52
int *nsalc, *ndisp, ndisp_all, nsalc_all, **ict;
53
double ***salc, ***disp;
56
print = optinfo.print_cartesians;
58
nirreps = syminfo.nirreps;
59
natom = carts.get_natom();
60
masses = carts.get_mass();
62
ndisp = init_int_array(nirreps);
63
nsalc = init_int_array(nirreps);
65
chkpt_init(PSIO_OPEN_OLD);
66
cartrep = chkpt_rd_cartrep();
68
natom_unique = chkpt_rd_num_unique_atom();
71
/* Read in data from PSIF_OPTKING */
73
psio_read_entry(PSIF_OPTKING, "OPT: Num. of disp.",
74
(char *) &(ndisp_all), sizeof(int));
75
psio_read_entry(PSIF_OPTKING, "OPT: Num. of coord.",
76
(char *) &(nsalc_all), sizeof(int));
77
psio_read_entry(PSIF_OPTKING, "OPT: Disp. per irrep",
78
(char *) &(ndisp[0]), nirreps*sizeof(int));
79
psio_read_entry(PSIF_OPTKING, "OPT: Coord. per irrep",
80
(char *) &(nsalc[0]), nirreps*sizeof(int));
82
fprintf(outfile,"\n total nsalc: %d, total ndisp: %d\n", nsalc_all, ndisp_all);
83
fprintf(outfile,"nsalc per irrep: "); for (h=0; h<nirreps; ++h) fprintf(outfile,"%d ",nsalc[h]);
84
fprintf(outfile,"\n");
85
fprintf(outfile,"ndisp per irrep: "); for (h=0; h<nirreps; ++h) fprintf(outfile,"%d ",ndisp[h]);
86
fprintf(outfile,"\n\n");
88
disp_grad = block_matrix(ndisp_all,3*natom);
89
if (optinfo.grad_dat) { /* read in gradients from "file11.dat" */
90
ffile(&fp11, "file11.dat", 2);
92
line1 = new char[MAX_LINE+1];
93
for (i=0;i<ndisp_all;++i) {
94
/* read in 2 header lines */
95
fgets(line1, MAX_LINE, fp11);
96
fgets(line1, MAX_LINE, fp11);
98
for (j=0; j<natom; ++j) {
99
fgets(line1, MAX_LINE, fp11);
101
/* read in xyz for N atoms */
102
for (j=0; j<natom; ++j) {
103
fgets(line1, MAX_LINE, fp11);
104
sscanf(line1, "%lf %lf %lf",
105
&(disp_grad[i][3*j]), &(disp_grad[i][3*j+1]), &(disp_grad[i][3*j+2]) );
111
psio_read_entry(PSIF_OPTKING, "OPT: Displaced gradients",
112
(char *) &(disp_grad[0][0]), ndisp_all*3*natom*sizeof(double));
115
fprintf(outfile,"Gradients of displaced geometries\n");
116
print_mat(disp_grad,ndisp_all,3*natom,outfile);
119
/**** construct full gradient list using unique ones ****/
120
if (optinfo.points == 3)
121
disp_grad_all = block_matrix(2*nsalc_all, 3*natom);
122
else if (optinfo.points == 5)
123
disp_grad_all = block_matrix(4*nsalc_all, 3*natom);
124
/* just copy symmetric gradients */
125
for (j=0; j<ndisp[0]; ++j) {
126
for (I=0; I<3*natom; ++I)
127
disp_grad_all[j][I] = disp_grad[j][I];
130
/* add gradients to list for asymmetric displacements */
131
cnt_all = cnt = ndisp[0] - 1;
132
for (h=1; h<nirreps; ++h) {
134
if (!ndisp[h]) continue;
136
/* copy computed displacement in */
137
for (j=0; j<ndisp[h]; ++j) {
140
for (I=0; I<3*natom; ++I)
141
disp_grad_all[cnt_all][I] = disp_grad[cnt][I];
143
/* determine which operation takes minus displacement into plus displacement */
144
for (op_disp=0; op_disp<syminfo.nirreps; ++op_disp) {
145
if ( syminfo.ct[h][op_disp] == -1 )
148
fprintf(outfile,"\tOperation that takes plus displacement %d to minus is %s.\n",
149
cnt+1, syminfo.op_lbls[op_disp]);
151
/* what is the parity of the irrep for this special operation ? */
152
//sign = syminfo.ct[h][op_disp];
153
//fprintf(outfile,"\tParity of irrep %s for this operation is %d.\n",
154
//syminfo.clean_irrep_lbls[h], sign);
157
for (atom=0; atom<natom; ++atom) {
158
for (xyz=0; xyz<3; ++xyz) {
159
atom2 = ict[op_disp][atom] - 1;
160
disp_grad_all[cnt_all][3*atom2+xyz] = disp_grad[cnt][3*atom+xyz] * cartrep[op_disp][3*xyz+xyz];
164
ndisp[h] += ndisp[h];
166
free_block(disp_grad);
168
/* fix number of displacements to match full list of gradients */
170
for (h=0; h<nirreps; ++h) {
171
ndisp_all += ndisp[h];
174
B = block_matrix(nsalc_all,3*natom);
175
psio_read_entry(PSIF_OPTKING, "OPT: Adapted cartesians",
176
(char *) &(B[0][0]), nsalc_all*3*natom*sizeof(double));
178
fprintf(outfile,"B matrix\n");
179
print_mat(B,nsalc_all,3*natom,outfile);
182
evals_all = init_array(nsalc_all);
183
evals_all_irrep = init_int_array(nsalc_all);
185
/**** loop over irreps of salcs and vibrations ****/
188
for (h=0; h<nirreps; ++h) {
190
if (!nsalc[h]) continue;
192
/* construct B_inv = (BuBt)-1 B u = G-1 B u */
193
fprintf(outfile,"Computing B^-1: ");
194
G = compute_G(&(B[start_salc]),nsalc[h],carts);
195
G_inv = symm_matrix_invert(G,nsalc[h],1,0);
198
/* mass-weight gradients; g_xm = (1/sqrt(m)) * g_x */
199
/* works for H2 what is more general? */
200
for (k=0; k<ndisp[h]; ++k)
201
for (j=0;j<3*natom;++j)
202
disp_grad_all[k+start_disp][j] /= sqrt(masses[j]);
205
fprintf(outfile,"Displaced gradients/sqrt(masses[j])\n");
206
print_mat(&(disp_grad_all[start_disp]),ndisp[h],3*natom,outfile);
210
// compute forces in internal coordinates, f_q = G_inv B u f_x
211
f_q = init_array(nsalc[h]);
212
temp_arr = init_array(nsalc[h]);
213
temp_arr2 = init_array(3*natom);
214
grads_adapted = block_matrix(ndisp[h],3*natom);
216
if (print) fprintf(outfile,"Gradients recomputed from internal coordinate gradients\n");
217
for (k=0; k<ndisp[h]; ++k) {
219
if (print) fprintf(outfile,"ndisp[h]: %d, start_disp: %d\n", ndisp[h], start_disp);
220
mmult(u,0,&(disp_grad_all[k+start_disp]),1,&temp_arr2,1,3*natom,3*natom,1,0);
221
mmult(&(B[start_salc]),0,&temp_arr2,1,&temp_arr,1,nsalc[h],3*natom,1,0);
222
mmult(G_inv,0,&temp_arr,1,&f_q,1,nsalc[h],nsalc[h],1,0);
224
for (j=0;j<3*natom;++j)
225
grads_adapted[k][j] = f_q[j];
227
/* test by transforming f_q back to cartesian forces and compare */
228
mmult(B,1,&(grads_adapted[k]),1,&temp_arr2,1,3*natom,nsalc[h],1,0);
229
if (print) print_mat2(&temp_arr2, 1, 3*natom, outfile);
237
/* fprintf(outfile,"Adapted gradients\n");
238
print_mat(grads_adapted,ndisps,ncoord,outfile);
241
force_constants = block_matrix(nsalc[h],nsalc[h]);
243
/*** Construct force constant matrix from finite differences of forces ***/
244
fprintf(outfile," ** Using %d-point formula.\n",optinfo.points);
245
if (optinfo.points == 3) {
246
for (i=0; i<nsalc[h]; ++i) {
247
for (j=0; j<nsalc[h]; ++j) {
248
force_constants[i][j] =
249
(grads_adapted[2*i+1][j] - grads_adapted[2*i][j])
250
/ (2.0 * optinfo.disp_size);
254
else if (optinfo.points == 5) {
255
for (i=0; i<nsalc[h]; ++i) {
256
for (j=0; j<nsalc[h]; ++j) {
257
force_constants[i][j] =
258
/* 1f(-2) - 1f(+2) - 8f(-1) + 8f(+1) / 12h */
259
( 1.0 * grads_adapted[4*i][j] - 1.0 * grads_adapted[4*i+1][j]
260
- 8.0 * grads_adapted[4*i+2][j] + 8.0 * grads_adapted[4*i+3][j] )
261
/ (12.0 * optinfo.disp_size);
266
fprintf(outfile,"Force Constants\n");
267
print_mat(force_constants, nsalc[h], nsalc[h], outfile);
270
start_salc += nsalc[h];
271
start_disp += ndisp[h];
275
/* masses = carts.get_mass();
276
for (i=0;i<3*natom;++i)
277
for (j=0;j<3*natom;++j) {
278
force_constants[i][j] *= 1.0 / sqrt(masses[i] * masses[j]);
280
fprintf(outfile,"Mass Weighted Force Constants 3N by 3N\n");
281
print_mat(force_constants, dim, dim, outfile);
284
/** Find eigenvalues of force constant matrix **/
285
evals = init_array(dim);
286
evects = block_matrix(dim, dim);
287
dgeev_optking(dim, force_constants, evals, evects);
288
free_block(force_constants);
289
sort(evals, evects, dim);
292
for (i=0; i<dim; ++i) {
294
evals_all[cnt_eval] = evals[i];
295
evals_all_irrep[cnt_eval] = h;
301
free_block(disp_grad_all);
303
sort_evals_all(nsalc_all,evals_all, evals_all_irrep);
305
/* convert evals from H/(kg bohr^2) to J/(kg m^2) = 1/s^2 */
306
/* v = 1/(2 pi c) sqrt( eval ) */
307
fprintf(outfile, "\n\t Harmonic Vibrational Frequencies in cm^(-1)\n");
308
fprintf(outfile, "\t--------------------------------------------\n");
309
k_convert = _hartree2J/(_bohr2m * _bohr2m * _amu2kg);
310
cm_convert = 1.0/(2.0 * _pi * _c * 100.0);
311
for(i=nsalc_all-1; i>-1; --i) {
312
if(evals_all[i] < 0.0)
313
fprintf(outfile, "\t %5s %10.3fi\n", syminfo.irrep_lbls[evals_all_irrep[i]],
314
cm_convert * sqrt(-k_convert * evals_all[i]));
316
fprintf(outfile, "\t %5s %10.3f\n", syminfo.irrep_lbls[evals_all_irrep[i]],
317
cm_convert * sqrt(k_convert * evals_all[i]));
319
fprintf(outfile, "\t----------------------------\n");
322
free(evals_all_irrep);
324
optinfo.disp_num = 0;
325
psio_write_entry(PSIF_OPTKING, "OPT: Current disp_num",
326
(char *) &(optinfo.disp_num),sizeof(int));
331
void sort_evals_all(int nsalc_all, double *evals_all, int *evals_all_irrep) {
332
int i,j,tmp_irrep,min_index;
335
for (j=0; j<nsalc_all; ++j) {
338
for (i=j; i<nsalc_all; ++i) {
339
if (evals_all[i] < min) {
344
tmp_eval = evals_all[j];
345
evals_all[j]=evals_all[min_index];
346
evals_all[min_index] = tmp_eval;
348
tmp_irrep = evals_all_irrep[j];
349
evals_all_irrep[j] = evals_all_irrep[min_index];
350
evals_all_irrep[min_index] = tmp_irrep;
added
removed
src/bin/optking/freq_grad_cart.cc
