← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/optking/cartesians.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/optking/cartesians.cc
 
1
/*** cartesians.cc : contains member functions for cartesian class ***/
 
2
 
 
3
#if HAVE_CMATH
 
4
# include <cmath>
 
5
#else
 
6
# include <math.h>
 
7
#endif
 
8
 
 
9
extern "C" {
 
10
#include <stdio.h>
 
11
#include <stdlib.h>
 
12
#include <string.h>
 
13
#include <libciomr/libciomr.h>
 
14
#include <libipv1/ip_lib.h>
 
15
#include <libchkpt/chkpt.h>
 
16
#include <physconst.h>
 
17
#include <masses.h>
 
18
}
 
19
 
 
20
#define EXTERN
 
21
#include "opt.h"
 
22
#include "cartesians.h"
 
23
#undef EXTERN
 
24
 
 
25
 
 
26
/*** FORCES returns forces in cartesian coordinates in aJ/Ang ***/ 
 
27
double *cartesians:: get_forces() {
 
28
  int i;
 
29
  double *f;
 
30
  f = init_array(3*natom);
 
31
  for (i=0;i<3*natom;++i)
 
32
    f[i] = -1.0 * grad[i] * _hartree2J * 1.0E18 / _bohr2angstroms;
 
33
  return f;
 
34
}
 
35
 
 
36
double *cartesians:: get_fforces() {
 
37
  int i;
 
38
  double *f;
 
39
  f = init_array(3*nallatom);
 
40
  for (i=0; i<3*nallatom; ++i) {
 
41
    f[i] = -1.0 * fgrad[i] * _hartree2J * 1.0E18 / _bohr2angstroms;
 
42
  }
 
43
  return f;
 
44
}
 
45
 
 
46
 
 
47
/*** CARTESIANS this is the class constructor for cartesian ***/ 
 
48
cartesians::cartesians() {
 
49
  int i, j, a, count, isotopes_given = 0, masses_given = 0;
 
50
  char label[MAX_LINELENGTH], line1[MAX_LINELENGTH];
 
51
  FILE *fp_11;
 
52
  double an,x,y,z,tval,tval1,tval2,tval3,tval4;
 
53
  double **geom, **fgeom, *zvals;
 
54
 
 
55
  natom = optinfo.natom;
 
56
  nallatom = optinfo.nallatom;
 
57
 
 
58
  /* All data (but maybe masses) now read from chkpt */
 
59
  chkpt_init(PSIO_OPEN_OLD);
 
60
  geom = chkpt_rd_geom();
 
61
  fgeom = chkpt_rd_fgeom();
 
62
 
 
63
  if ((optinfo.mode == MODE_OPT_STEP) || (optinfo.mode == MODE_GRAD_SAVE))
 
64
    grad = chkpt_rd_grad();
 
65
  else grad = init_array(3*natom);
 
66
 
 
67
  if ((optinfo.mode == MODE_OPT_STEP) || (optinfo.mode == MODE_ENERGY_SAVE)
 
68
    ||(optinfo.mode == MODE_DISP_NOSYMM) || (optinfo.mode == MODE_DISP_IRREP)
 
69
    ||(optinfo.mode == MODE_DISP_LOAD) || (optinfo.mode == MODE_DISP_USER))
 
70
      energy = chkpt_rd_etot();
 
71
 
 
72
  zvals = chkpt_rd_zvals();
 
73
  chkpt_close();
 
74
 
 
75
  coord = new double[3*natom];
 
76
  atomic_num = new double[natom];
 
77
  mass = new double[3*natom];
 
78
  fcoord = new double[3*nallatom];
 
79
  fgrad = new double[3*nallatom];
 
80
  fatomic_num = new double[nallatom];
 
81
  fmass = new double[3*nallatom];
 
82
 
 
83
  count = -1;
 
84
  for(i=0; i < natom; i++) {
 
85
    atomic_num[i] = zvals[i];
 
86
    for(j=0; j < 3; j++) {
 
87
      coord[++count] = geom[i][j];
 
88
    }
 
89
  }
 
90
  free_block(geom);
 
91
 
 
92
  zero_arr(fatomic_num,nallatom);
 
93
  for(i=0; i<natom; i++)
 
94
    fatomic_num[optinfo.to_dummy[i]] = zvals[i];
 
95
 
 
96
  count = -1;
 
97
  for(i=0; i<nallatom; i++)
 
98
    for(j=0; j<3; j++)
 
99
      fcoord[++count] = fgeom[i][j];
 
100
 
 
101
  free(zvals);
 
102
  free_block(fgeom);
 
103
 
 
104
  zero_arr(fgrad, 3*nallatom);
 
105
  count = -1;
 
106
  for(i=0; i < natom; i++)
 
107
    for(j=0; j < 3; j++)
 
108
      fgrad[3*optinfo.to_dummy[i]+j]  = grad[3*i+j];
 
109
 
 
110
  /* read masses from input.dat or use default masses */
 
111
  count = -1;
 
112
 
 
113
  if (ip_exist("ISOTOPES",0)) {
 
114
    isotopes_given = 1;
 
115
    a = 0;
 
116
    ip_count("ISOTOPES", &a, 0);
 
117
    if (a != natom) {
 
118
      fprintf(outfile,"ISOTOPES array has wrong dimension.\n");
 
119
      exit(2);
 
120
    }
 
121
    for (i=0;i<natom;++i) {
 
122
      ip_data("ISOTOPES","%s", line1,1,i);
 
123
      for (j=0;j<138;j++) {
 
124
        if (strcmp(line1, mass_labels[j]) == 0) {
 
125
          mass[++count] = atomic_masses[j];
 
126
          mass[++count] = atomic_masses[j];
 
127
          mass[++count] = atomic_masses[j];
 
128
          break;
 
129
        }
 
130
      }
 
131
      if (j == 138) {
 
132
        fprintf(outfile,
 
133
                "Isotope label %s is unidentifiable.\n",mass_labels[j]);
 
134
        exit(2);
 
135
      }
 
136
    }
 
137
  }
 
138
  if (ip_exist("MASSES",0)) {
 
139
    masses_given = 1;
 
140
    if (isotopes_given)
 
141
      fprintf(outfile,"Ignoring ISOTOPES keyword, using given MASSES.\n");
 
142
    a = 0;
 
143
    ip_count("MASSES",&a,0);
 
144
    if (a != natom) {
 
145
      fprintf(outfile,"MASSES array has wrong dimension\n");
 
146
      exit(2);
 
147
    }
 
148
    else {
 
149
      for(i=0;i<natom;++i) {
 
150
        ip_data("MASSES","%lf",&tval,1,i);
 
151
        mass[++count] = tval;
 
152
        mass[++count] = tval;
 
153
        mass[++count] = tval;
 
154
      }
 
155
    }
 
156
  }
 
157
  if ((isotopes_given == 0) && (masses_given == 0)) {
 
158
    for(i=0;i<natom;++i) {
 
159
      a = (int) get_atomic_num(i); // casting to an int for index
 
160
      mass[++count] = an2masses[a];
 
161
      mass[++count] = an2masses[a];
 
162
      mass[++count] = an2masses[a];
 
163
    }
 
164
  }
 
165
 
 
166
  zero_arr(fmass, 3*nallatom);
 
167
  for (i=0; i<natom; ++i) {
 
168
    j = optinfo.to_dummy[i];
 
169
    fmass[3*j+0] = fmass[3*j+1] = fmass[3*j+2] = mass[3*i];
 
170
  }
 
171
 
 
172
/*
 
173
  fprintf(outfile,"masses without dummy atoms\n");
 
174
  for (i=0; i<nallatom; ++i)
 
175
    fprintf(outfile,"%15.10lf%15.10lf%15.10lf\n",
 
176
        fmass[3*i],fmass[3*i+1],fmass[3*i+2]);
 
177
     */
 
178
 
 
179
  return;
 
180
}
 
181
 
 
182
/*** CARTESIANS::PRINT
 
183
 * print_flag                                         to
 
184
 *    0        fcoord                                fp_out
 
185
 *    1        fatomic_num, fcoord                   fp_out
 
186
 *    2        fatomic_num, fmass, fcoord            fp_out
 
187
 *    3        fatomic_num, fmass, fcoord, fgrad     fp_out
 
188
 *    4        coord                                 fp_out
 
189
 *    5        atomic_num, coord                     fp_out
 
190
 *    6        atomic_num, mass, coord               fp_out
 
191
 *    7        atomic_num, mass, coord, grad         fp_out
 
192
 *   10        coord                                 geom.dat
 
193
 *   11        coord, grad                           file11.dat
 
194
 *   32        fcoord (implicitly coord)             chkpt
 
195
 *   12        fatomic_symb, fcoord                  fp_out
 
196
 * disp_label is only used for geom.dat writing ***/
 
197
 
 
198
void cartesians :: print(int print_flag, FILE *fp_out, int new_geom_file,
 
199
                         char *disp_label, int disp_num) {
 
200
  int i,j;
 
201
  double x,y,z;
 
202
  int cnt = -1;
 
203
  char sym[3];
 
204
 
 
205
  if (print_flag == 0) {
 
206
    for (i = 0; i < natom; ++i) {
 
207
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
208
      fprintf(fp_out,"%20.10f%20.10f%20.10f\n",x,y,z);
 
209
    }
 
210
  }
 
211
  else if (print_flag == 1) {
 
212
    for (i = 0; i < natom; ++i) {
 
213
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
214
      fprintf(fp_out,"%5.1lf%15.10f%15.10f%15.10f\n",atomic_num[i],x,y,z);
 
215
    }
 
216
  }
 
217
  else if (print_flag == 2) {
 
218
    cnt = -1;
 
219
    for (i = 0; i < natom; ++i) {
 
220
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
221
      fprintf(fp_out,
 
222
          "%5.1lf%15.8lf%15.10f%15.10f%15.10f\n",atomic_num[i],mass[3*i],x,y,z);
 
223
    }
 
224
  }
 
225
  else if (print_flag == 3) {
 
226
    cnt = -1;
 
227
    for (i = 0; i < natom; ++i) {
 
228
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
229
      fprintf(fp_out,
 
230
          "%5.1lf%15.8lf%15.10f%15.10f%15.10f\n",atomic_num[i],mass[3*i],x,y,z);
 
231
    }
 
232
    cnt = -1;
 
233
    for (i = 0; i < natom; ++i) {
 
234
      x = grad[++cnt]; y = grad[++cnt]; z = grad[++cnt];
 
235
      fprintf(fp_out,
 
236
              "%35.10lf%15.10f%15.10f\n",x,y,z);
 
237
    }
 
238
  }
 
239
  if (print_flag == 4) {
 
240
    for (i = 0; i < natom; ++i) {
 
241
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
242
      fprintf(fp_out,"%20.10f%20.10f%20.10f\n",x,y,z);
 
243
    }
 
244
  } 
 
245
  else if (print_flag == 5) {
 
246
    for (i = 0; i < natom; ++i) {
 
247
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
248
      fprintf(fp_out,"%5.1lf%15.10f%15.10f%15.10f\n",atomic_num[i],x,y,z);
 
249
    }
 
250
  }
 
251
  else if (print_flag == 6) {
 
252
    cnt = -1;
 
253
    for (i = 0; i < natom; ++i) {
 
254
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
255
      fprintf(fp_out,
 
256
          "%5.1lf%15.8lf%15.10f%15.10f%15.10f\n",atomic_num[i],mass[3*i],x,y,z);
 
257
    }
 
258
  }
 
259
  else if (print_flag == 7) {
 
260
    cnt = -1;
 
261
    for (i = 0; i < natom; ++i) {
 
262
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
263
      fprintf(fp_out,
 
264
          "%5.1lf%15.8lf%15.10f%15.10f%15.10f\n",atomic_num[i],mass[3*i],x,y,z);
 
265
    }
 
266
    cnt = -1;
 
267
    for (i = 0; i < natom; ++i) {
 
268
      x = grad[++cnt]; y = grad[++cnt]; z = grad[++cnt];
 
269
      fprintf(fp_out,
 
270
              "%35.10lf%15.10f%15.10f\n",x,y,z);
 
271
    }
 
272
  }
 
273
  else if (print_flag == 10) {
 
274
    //if (new_geom_file) fprintf(fp_out,"%%%%\n");
 
275
    fprintf(fp_out,"%% %s\n",disp_label);
 
276
    fprintf(fp_out,"geometry%1d = (\n", disp_num);
 
277
    cnt = -1;
 
278
    for (i=0; i<natom; ++i) {
 
279
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
280
      fprintf(fp_out,"(%20.10f%20.10f%20.10f)\n",x,y,z);
 
281
    }
 
282
    fprintf(fp_out," )\n");
 
283
  }
 
284
  else if (print_flag == 11) {
 
285
    fprintf(fp_out,"%s\n",disp_label);
 
286
    fprintf(fp_out,"%5d%20.10lf\n",natom,energy);
 
287
    cnt = -1;
 
288
    for (i = 0; i < natom; ++i) {
 
289
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
290
      fprintf(fp_out,
 
291
              "%20.10lf%20.10lf%20.10lf%20.10lf\n",atomic_num[i],x,y,z);
 
292
    }
 
293
    cnt = -1;
 
294
    for (i = 0; i < natom; ++i) {
 
295
      x = grad[++cnt]; y = grad[++cnt]; z = grad[++cnt];
 
296
      fprintf(fp_out,"%40.10lf%20.10lf%20.10lf\n",x,y,z);
 
297
    }
 
298
  }
 
299
  else if (print_flag == 32) {
 
300
 
 
301
    fprintf(outfile,"\nGeometry written to chkpt\n");
 
302
 
 
303
    double **geom;
 
304
    geom = block_matrix(nallatom,3);
 
305
    for (i=0; i<nallatom; ++i)
 
306
      for (j=0; j<3; ++j)
 
307
        geom[optinfo.to_dummy[i]][j] = coord[3*i+j];
 
308
    chkpt_init(PSIO_OPEN_OLD);
 
309
    chkpt_wt_fgeom(geom);
 
310
    chkpt_close();
 
311
    free_block(geom);
 
312
 
 
313
    /*
 
314
    double **tmp2d;
 
315
    chkpt_init(PSIO_OPEN_OLD);
 
316
    tmp2d = chkpt_rd_fgeom();
 
317
    chkpt_close();
 
318
    for (i=0; i<nallatom; ++i)
 
319
      for (j=0; j<3; ++j)
 
320
        fprintf(outfile,"%15.10lf\n", tmp2d[i][j]);
 
321
    free_block(tmp2d);
 
322
    */
 
323
 
 
324
  }
 
325
  else if (print_flag == 12) { 
 
326
    for (i = 0; i < natom; ++i) {
 
327
      x = coord[++cnt]; y = coord[++cnt]; z = coord[++cnt];
 
328
      zval_to_symbol(atomic_num[i],sym);
 
329
      fprintf(fp_out,"  (%3s%15.10f%15.10f%15.10f )\n",sym,x,y,z);
 
330
    }
 
331
  }
 
332
  return;
 
333
}
 
334