16
#include <libciomr/libciomr.h>
17
#include <libipv1/ip_lib.h>
18
#include <physconst.h>
20
// extern FILE *outfile;
23
extern InternalCoordinates gIntCo;
24
extern Params gParams;
27
extern void open_PSIF();
28
extern void close_PSIF();
30
extern void secondDerivativeOut(int, double **);
31
extern void secondDerivativeIn(int, double **);
33
Transform::Transform()
35
/* Numerical testing stuff that is untested
37
if(gParams.derlvl + gParams.atEquilibrium > 1) {
38
gParams.intco2ndDer = NULL;
40
gParams.intco2ndDer = init_matrix(gParams.ncartesians, gParams.ncartesians);
41
SR2Matrix = init_matrix(gParams.ncartesians, gParams.ncartesians);
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));
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));
60
Transform::~Transform()
62
/* if(gParams.derlvl + gParams.atEquilibrium > 1) {
63
//~SecondDerivative();
65
if(gParams.derlvl + gParams.atEquilibrium > 2) {
68
if(gParams.derlvl + gParams.atEquilibrium > 3) {
74
/* linear_transform -- Performs a linear transformation from internal to cartesians *
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 *
79
* if transformType > 0, coord1 = gParams.ncartesians, coord2 = gParams.nintco *
80
* transformMatrix = BMatrix *
82
* if transformType < 0, coord1 = gParams.nintco, coord2 = gParams.ncartesians *
83
* transformMatrix = AMatrix *
84
* transformType < 0 is not yet implemented! */
86
void Transform::der_transform(int coord1, int coord2, double **transformMatrix)
89
double cf1, cf2, cf3, cf4; //conversion factors, which will be appropriate cf's or equal to 1
90
double *temp_mass_array;
92
temp_mass_array = init_array(gParams.nintco);
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;
107
if(gParams.transformType >= 1) {
114
if(gParams.atEquilibrium) {
115
if(gParams.transformType >= 1)
116
gParams.FConst1 = init_array(coord1);
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];
126
SecondDerivative::SecondDerivative()
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);
140
SecondDerivative::~SecondDerivative()
153
free_matrix(SR2Matrix, gParams.ncartesians);
156
ThirdDerivative::ThirdDerivative()
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);
181
ThirdDerivative::~ThirdDerivative()
206
FourthDerivative::FourthDerivative()
211
FourthDerivative::~FourthDerivative()
217
void SecondDerivative::i_to_c()
222
XRow(SR2Matrix, 0, r);
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);
231
fprintf(outfile, "SR(I,J) and X(M,N) Matrices set to zero for simple internal coordinate %i\n",r);
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
239
void ThirdDerivative::i_to_c()
243
YRow(SR3Matrix, 0, r);
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);
251
fprintf(outfile, "SR(I,J,K) and Y(M,N,P) Matrices set to zero for simple internal coordinate %i\n",r);
253
fill3a(gParams.ncartesians, gParams.ncartesians, gParams.intco3rdDer);
254
thirdDerivativeOut(-r, SR3Matrix);
255
thirdDerivativeOut(r, gParams.intco3rdDer);
258
void FourthDerivative::i_to_c()