~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/intder/transform.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
#include "intco.h"
 
2
#include "transform.h"
 
3
#include "params.h"
 
4
#include "atom.h"
 
5
#include "bmat.h"
 
6
 
 
7
#define EXTERN
 
8
#include "globals.h"
 
9
#include "3dmatrix.h"
 
10
 
 
11
#include <stdio.h>
 
12
#include <stdlib.h>
 
13
#include <math.h>
 
14
 
 
15
extern "C" {
 
16
  #include <libciomr/libciomr.h>
 
17
  #include <libipv1/ip_lib.h>
 
18
  #include <physconst.h>
 
19
 
 
20
//  extern FILE *outfile;
 
21
};
 
22
 
 
23
extern InternalCoordinates gIntCo;
 
24
extern Params gParams;
 
25
extern BMat gBMat;
 
26
 
 
27
extern void open_PSIF();
 
28
extern void close_PSIF();
 
29
  
 
30
extern void secondDerivativeOut(int, double **);
 
31
extern void secondDerivativeIn(int, double **);
 
32
 
 
33
Transform::Transform()
 
34
{
 
35
/* Numerical testing stuff that is untested 
 
36
 
 
37
  if(gParams.derlvl + gParams.atEquilibrium > 1) {
 
38
    gParams.intco2ndDer = NULL;
 
39
    SR2Matrix = NULL;
 
40
    gParams.intco2ndDer = init_matrix(gParams.ncartesians, gParams.ncartesians);
 
41
    SR2Matrix = init_matrix(gParams.ncartesians, gParams.ncartesians);
 
42
  }
 
43
  if(gParams.derlvl + gParams.atEquilibrium > 2) {
 
44
    gParams.intco3rdDer = new C3DMatrix(gParams.ncartesians, gParams.ncartesians, gParams.ncartesians);
 
45
    gParams.intco3rdDer->array = init_array(gParams.intco3rdDer->Size(gParams.ncartesians,gParams.ncartesians,gParams.ncartesians));
 
46
    SR3Matrix = new C3DMatrix(gParams.ncartesians, gParams.ncartesians, gParams.ncartesians);
 
47
    SR3Matrix->array = init_array(SR3Matrix->Size(gParams.ncartesians,gParams.ncartesians,gParams.ncartesians));
 
48
  }
 
49
  if(gParams.derlvl + gParams.atEquilibrium > 3) {
 
50
    gParams.intco4thDer = new C4DMatrix(gParams.ncartesians, gParams.ncartesians, gParams.ncartesians, gParams.ncartesians);
 
51
    gParams.intco4thDer->array = init_array(gParams.intco4thDer->Size(gParams.ncartesians,gParams.ncartesians,gParams.ncartesians, gParams.ncartesians));
 
52
    gParams.intco4thDerBlock = new C3DMatrix(gParams.ncartesians, gParams.ncartesians, gParams.ncartesians);
 
53
    gParams.intco4thDerBlock->array = init_array(gParams.intco4thDerBlock->Size(gParams.ncartesians,gParams.ncartesians, gParams.ncartesians));
 
54
    SR4Matrix = new C4DMatrix(gParams.ncartesians, gParams.ncartesians, gParams.ncartesians, gParams.ncartesians);
 
55
    SR4Matrix->array = init_array(SR4Matrix->Size(gParams.ncartesians,gParams.ncartesians,gParams.ncartesians, gParams.ncartesians));
 
56
  }  */
 
57
 
 
58
}
 
59
 
 
60
Transform::~Transform()
 
61
{
 
62
/*  if(gParams.derlvl + gParams.atEquilibrium > 1) {
 
63
    //~SecondDerivative();
 
64
  }
 
65
  if(gParams.derlvl + gParams.atEquilibrium > 2) {
 
66
    //~ThirdDerivative();
 
67
  }
 
68
  if(gParams.derlvl + gParams.atEquilibrium > 3) {
 
69
    //FourthDerDelete();
 
70
  } 
 
71
*/
 
72
}
 
73
 
 
74
/* linear_transform -- Performs a linear transformation from internal to cartesians *
 
75
 * or vice versa                                                                    *
 
76
 * Passed in variables are dependent on whether or not atomic masses are            *
 
77
 * included, i.e. if transformType is greater than or less than zero                *
 
78
 *                                                                                  *
 
79
 * if transformType > 0, coord1 = gParams.ncartesians, coord2 = gParams.nintco      *
 
80
 * transformMatrix = BMatrix                                                        *
 
81
 
 
82
 * if transformType < 0, coord1 = gParams.nintco, coord2 = gParams.ncartesians      *
 
83
 * transformMatrix = AMatrix                                                        *
 
84
 * transformType < 0 is not yet implemented!                                        */
 
85
 
 
86
void Transform::der_transform(int coord1, int coord2, double **transformMatrix)
 
87
{
 
88
  int i,j,m;
 
89
  double cf1, cf2, cf3, cf4; //conversion factors, which will be appropriate cf's or equal to 1
 
90
  double *temp_mass_array;
 
91
 
 
92
  temp_mass_array = init_array(gParams.nintco);
 
93
 
 
94
  if(abs(gParams.transformType != 3)) {
 
95
    cf1 = _hartree2J / (_bohr2angstroms * 1.0E-18);
 
96
    cf2 = cf1 / _bohr2angstroms;
 
97
    cf3 = cf2 / _bohr2angstroms;
 
98
    cf4 = cf3 / _bohr2angstroms;
 
99
  }
 
100
  else {
 
101
    cf1 = 1.0;
 
102
    cf2 = 1.0;
 
103
    cf3 = 1.0;
 
104
    cf4 = 1.0;
 
105
  }
 
106
  
 
107
  if(gParams.transformType >= 1) {
 
108
    cf1 = 1.0;
 
109
    cf2 = 1.0;
 
110
    cf3 = 1.0;
 
111
    cf4 = 1.0;
 
112
  }  //this seems dumb
 
113
  
 
114
  if(gParams.atEquilibrium) {
 
115
    if(gParams.transformType >= 1)
 
116
      gParams.FConst1 = init_array(coord1);
 
117
      }
 
118
  for(i = 0; i < coord2; i++) {
 
119
    temp_mass_array[i] = gParams.mass_array[i] * cf1;
 
120
    for(m = 0; m < coord1; m++)
 
121
      gParams.FConst1[m] += temp_mass_array[i] * gBMat.AMatrix[i][m];
 
122
  }
 
123
}
 
124
 
 
125
/*
 
126
SecondDerivative::SecondDerivative()
 
127
{
 
128
  H11 = init_matrix(3,3);
 
129
  H21 = init_matrix(3,3);
 
130
  H22 = init_matrix(3,3);
 
131
  H31 = init_matrix(3,3);
 
132
  H32 = init_matrix(3,3);
 
133
  H33 = init_matrix(3,3);
 
134
  H41 = init_matrix(3,3);
 
135
  H42 = init_matrix(3,3);
 
136
  H43 = init_matrix(3,3);
 
137
  H44 = init_matrix(3,3);
 
138
}
 
139
 
 
140
SecondDerivative::~SecondDerivative()
 
141
{
 
142
  free_matrix(H11,3);
 
143
  free_matrix(H21,3);
 
144
  free_matrix(H22,3);
 
145
  free_matrix(H31,3);
 
146
  free_matrix(H32,3);
 
147
  free_matrix(H33,3);
 
148
  free_matrix(H41,3);
 
149
  free_matrix(H42,3);
 
150
  free_matrix(H43,3);
 
151
  free_matrix(H44,3);
 
152
 
 
153
  free_matrix(SR2Matrix, gParams.ncartesians);
 
154
}
 
155
 
 
156
ThirdDerivative::ThirdDerivative()
 
157
{
 
158
  H111 = new C3DMatrix(3,3,3);
 
159
  H112 = new C3DMatrix(3,3,3); 
 
160
  H221 = new C3DMatrix(3,3,3); 
 
161
  H222 = new C3DMatrix(3,3,3); 
 
162
  H113 = new C3DMatrix(3,3,3); 
 
163
  H123 = new C3DMatrix(3,3,3); 
 
164
  H223 = new C3DMatrix(3,3,3); 
 
165
  H331 = new C3DMatrix(3,3,3); 
 
166
  H332 = new C3DMatrix(3,3,3); 
 
167
  H333 = new C3DMatrix(3,3,3); 
 
168
  H411 = new C3DMatrix(3,3,3); 
 
169
  H421 = new C3DMatrix(3,3,3); 
 
170
  H422 = new C3DMatrix(3,3,3);
 
171
  H431 = new C3DMatrix(3,3,3); 
 
172
  H432 = new C3DMatrix(3,3,3); 
 
173
  H433 = new C3DMatrix(3,3,3); 
 
174
  H441 = new C3DMatrix(3,3,3); 
 
175
  H442 = new C3DMatrix(3,3,3); 
 
176
  H443 = new C3DMatrix(3,3,3); 
 
177
  H444 = new C3DMatrix(3,3,3); 
 
178
 
 
179
}
 
180
 
 
181
ThirdDerivative::~ThirdDerivative()
 
182
{
 
183
  delete H111;
 
184
  delete H112;
 
185
  delete H221;
 
186
  delete H222;
 
187
  delete H113;
 
188
  delete H123;
 
189
  delete H223;
 
190
  delete H331;
 
191
  delete H332;
 
192
  delete H333;
 
193
  delete H411;
 
194
  delete H421;
 
195
  delete H431;
 
196
  delete H432;
 
197
  delete H433;
 
198
  delete H441;
 
199
  delete H442;
 
200
  delete H443;
 
201
  delete H444;
 
202
 
 
203
  delete SR3Matrix;
 
204
}
 
205
 
 
206
FourthDerivative::FourthDerivative()
 
207
{
 
208
 
 
209
}
 
210
 
 
211
FourthDerivative::~FourthDerivative()
 
212
{
 
213
 
 
214
}
 
215
 
 
216
//This will be MACHX
 
217
void SecondDerivative::i_to_c()
 
218
{
 
219
  int i,j,r;
 
220
  double s = 0.0;
 
221
  
 
222
  XRow(SR2Matrix, 0, r); 
 
223
  
 
224
  //Does this command ever get called in intder2000.f?!
 
225
  if(gParams.intcoInclude[r]) {
 
226
    secondDerivativeIn(-r, SR2Matrix);
 
227
    secondDerivativeIn(r, gParams.intco2ndDer);
 
228
    fprintf(outfile, "Numerical SR(I,J) and X(M,N) matrices used for simple internal coordinate %i\n",r);
 
229
  } 
 
230
  else
 
231
    fprintf(outfile, "SR(I,J) and X(M,N) Matrices set to zero for simple internal coordinate %i\n",r);
 
232
  
 
233
  secondDerivativeOut(-r, SR2Matrix);
 
234
  secondDerivativeOut(r, gParams.intco2ndDer);
 
235
  //Again, matrix dimensions are not right. Also need to add lines dealing with symmetric internals
 
236
}
 
237
 
 
238
//MACHY
 
239
void ThirdDerivative::i_to_c()
 
240
{
 
241
  int i,j,r;
 
242
  
 
243
  YRow(SR3Matrix, 0, r);
 
244
  
 
245
  if(gParams.intcoInclude[r]) {
 
246
    thirdDerivativeIn(-r, SR3Matrix);
 
247
    thirdDerivativeIn(r, gParams.intco3rdDer);
 
248
    fprintf(outfile, "Numerical SR(I,J,K) and Y(M,N,P) matrices used for simple internal coordinate %i\n",r);
 
249
  } 
 
250
  else
 
251
    fprintf(outfile, "SR(I,J,K) and Y(M,N,P) Matrices set to zero for simple internal coordinate %i\n",r);
 
252
  
 
253
  fill3a(gParams.ncartesians, gParams.ncartesians, gParams.intco3rdDer);
 
254
  thirdDerivativeOut(-r, SR3Matrix);
 
255
  thirdDerivativeOut(r, gParams.intco3rdDer);
 
256
}
 
257
 
 
258
void FourthDerivative::i_to_c()
 
259
{
 
260
 
 
261
 
 
262
} */
 
263