Skip to content

Commit 00d06f8

Browse files
cthree-40cthree-40
committed
Added starting files/routines to compute transition moments.
1 parent efba681 commit 00d06f8

13 files changed

+1664
-0
lines changed

‎mpi_source/atomic_orbitals.c

+313
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,313 @@
1+
// File: atomic_orbitals.c
2+
/*
3+
* atomic_orbitals: contains atomic orbital evaluation functions.
4+
* --------------------------------------------------------------
5+
*/
6+
7+
#include <stdio.h>
8+
#include <stdlib.h>
9+
#include <math.h>
10+
#include "errorlib.h"
11+
#include "buildao.h"
12+
#include "mathutil.h"
13+
#include "atomic_orbitals.h"
14+
15+
/*
16+
* compute_d_normconst: Compute the d orbital normalization constant.
17+
* (2048a^7/Pi^3)^(1/4).
18+
*/
19+
double compute_d_normconst(double alpha)
20+
{
21+
double val;
22+
val = 2048 * pow(alpha, 7) / M_PI_3;
23+
val = pow(val, 1.0/4.0);
24+
return val;
25+
}
26+
27+
/*
28+
* compute_f_normconst: Compute the d orbital normalization constant.
29+
* (32768a^9/Pi^3)^(1/4).
30+
*/
31+
double compute_f_normconst(double alpha)
32+
{
33+
double val;
34+
val = 32768 * pow(alpha, 9) / M_PI_3;
35+
val = pow(val, 1.0/4.0);
36+
return val;
37+
}
38+
39+
/*
40+
* compute_p_normconst: Compute the p orbital normalization constant.
41+
* (128a^5/Pi^3)^(1/4).
42+
*/
43+
double compute_p_normconst(double alpha)
44+
{
45+
double val;
46+
val = 128 * pow(alpha, 5) / M_PI_3;
47+
val = pow(val, 1.0/4.0);
48+
return val;
49+
}
50+
51+
/*
52+
* compute_s_normconst: Compute the s orbital normalization constant.
53+
* (8a^3/PI^3)^(1/4).
54+
*/
55+
double compute_s_normconst(double alpha)
56+
{
57+
double val;
58+
val = 8 * pow(alpha, 3) / M_PI_3;
59+
val = pow(val, 1.0/4.0);
60+
return val;
61+
}
62+
63+
/*
64+
* evaluate_orbital: evaluate an atomic orbital at some position (x, y, z).
65+
*
66+
* Input:
67+
* orb = orbital basis function
68+
* pos = (x, y, z)
69+
* scr = scratch array to hold |R - r| (3 * sizeof(double))
70+
*/
71+
double evaluate_orbital(struct ao_basisfunc orb, double *pos, double *scr)
72+
{
73+
double value = 0.0; /* Return value */
74+
75+
switch (orb.type) {
76+
case 1:
77+
value = evaluate_s(pos, orb.alpha, orb.ccoef, orb.ugaus,
78+
orb.gscale, orb.geom);
79+
break;
80+
case 2: case 3: case 4:
81+
value = evaluate_p(pos, orb.alpha, orb.ccoef, orb.ugaus,
82+
orb.gscale, orb.geom, orb.type, scr);
83+
break;
84+
case 5: case 6: case 7: case 8: case 9:
85+
value = evaluate_d(pos, orb.alpha, orb.ccoef, orb.ugaus,
86+
orb.gscale, orb.geom, orb.type, scr);
87+
break;
88+
case 10: case 11: case 12: case 13: case 14: case 15: case 16:
89+
value = evaluate_f(pos, orb.alpha, orb.ccoef, orb.ugaus,
90+
orb.gscale, orb.geom, orb.type, scr);
91+
break;
92+
default:
93+
error_message(0, "Unknown orbital type", "evaluate_orbital");
94+
break;
95+
}
96+
97+
return value;
98+
}
99+
100+
/*
101+
* evaluate_d: evaluate a contracted d orbital at the position (x, y, z).
102+
*
103+
* Input:
104+
* pos = x, y, z to evalute orbital at
105+
* alphac = alpha coefficients for contracted gaussians
106+
* ccoef = contraction coefficients
107+
* ng = number of gaussians in contraction
108+
* gscale = gaussian scaling coefficient
109+
* atompos = x, y, z of nuclear center
110+
* typ = 5=d2-(xy), 6=d1-(xz), 7=d0 (z2), 8=d1+(yz), 9=d2+(x2-y2)
111+
* scr = scratch array to hold |R - r|
112+
*/
113+
double evaluate_d(double *pos, double *alphac, double *ccoef, int ng,
114+
double gscale, double *atompos, int typ, double *scr)
115+
{
116+
double value = 0.0; /* Return value */
117+
double rval;
118+
double dcmpnt; /* d_{type} component {type} (xy, xz, z2, yz, x2-y2) */
119+
double nrmcnst;
120+
int i;
121+
/* Compute {X-x,Y-y,Z-z} and ||{X-x,Y-y,Z-z}||*/
122+
vector_difference(pos, atompos, 3, scr);
123+
rval = compute_vector_norm(scr, 3);
124+
switch (typ) {
125+
case 5:
126+
/* xy */
127+
dcmpnt = scr[0] * scr[1];
128+
break;
129+
case 6:
130+
/* xz */
131+
dcmpnt = scr[0] * scr[2];
132+
break;
133+
case 7:
134+
/* z2 = 2z^2 - x^2 - y^2*/
135+
dcmpnt = 2 * scr[2] * scr[2];
136+
dcmpnt -= scr[0] * scr[0];
137+
dcmpnt -= scr[1] * scr[1];
138+
break;
139+
case 8:
140+
/* yz */
141+
dcmpnt = scr[1] * scr[2];
142+
break;
143+
case 9:
144+
/* x2 - y2 */
145+
dcmpnt = scr[0] * scr[0];
146+
dcmpnt -= scr[1] * scr[1];
147+
break;
148+
default:
149+
error_message(0, "Unknown type!", "evaluate_d");
150+
printf(" type = %d\n", typ);
151+
return value;
152+
}
153+
for (i = 0; i < ng; i++) {
154+
nrmcnst = compute_d_normconst(alphac[i]);
155+
value = value + (gscale * ccoef[i] * nrmcnst *
156+
dcmpnt * gauss_fcn(alphac[i], rval));
157+
}
158+
return value;
159+
}
160+
161+
/*
162+
* evaluate_f: evaluate a contracted f orbital at the position (x, y, z).
163+
*
164+
* Input:
165+
* pos = x, y, z to evalute orbital at
166+
* alphac = alpha coefficients for contracted gaussians
167+
* ccoef = contraction coefficients
168+
* ng = number of gaussians in contraction
169+
* gscale = gaussian scaling coefficient
170+
* atompos = x, y, z of nuclear center
171+
* typ = 10=f3-, 11=f2-, 12=f1-, 13=f0,
172+
* 14=f1+, 15=f2+, 16=f3+
173+
* scr = scratch array to hold |R - r|
174+
*/
175+
double evaluate_f(double *pos, double *alphac, double *ccoef, int ng,
176+
double gscale, double *atompos, int typ, double *scr)
177+
{
178+
double value = 0.0; /* Return value */
179+
double rval;
180+
double fcmpnt; /* f_{type} component {type} () */
181+
double nrmcnst;
182+
int i;
183+
/* Compute {X-x,Y-y,Z-z} and ||{X-x,Y-y,Z-z}||*/
184+
vector_difference(pos, atompos, 3, scr);
185+
rval = compute_vector_norm(scr, 3);
186+
switch (typ) {
187+
case 10:
188+
/* m=-3, xxy - yyy */
189+
fcmpnt = pow(scr[0], 2) * scr[1];
190+
fcmpnt -= pow(scr[1], 3);
191+
break;
192+
case 11:
193+
/* m=-2, xyz */
194+
fcmpnt = scr[0] * scr[1] * scr[2];
195+
break;
196+
case 12:
197+
/* m=-1, 4*yzz - 3*yyy - xxy */
198+
fcmpnt = 4 * scr[1] * scr[2] * scr[2];
199+
fcmpnt -= 3 * pow(scr[1], 3);
200+
fcmpnt -= pow(scr[0], 2) * scr[1];
201+
break;
202+
case 13:
203+
/* m=0, xxz + yyz - 2*zzz */
204+
fcmpnt = pow(scr[0], 2);
205+
fcmpnt += pow(scr[1], 2) * scr[2];
206+
fcmpnt -= 2 * pow(scr[2], 3);
207+
break;
208+
case 14:
209+
/* m=+1, 4*xzz - 3*xxx - xyy */
210+
fcmpnt = 4 * scr[0] * pow(scr[2], 2);
211+
fcmpnt -= 3 * pow(scr[0], 3);
212+
fcmpnt -= scr[0] * pow(scr[1], 2);
213+
break;
214+
case 15:
215+
/* m=+2, xxz - yyz */
216+
fcmpnt = pow(scr[0], 2) * scr[2];
217+
fcmpnt -= pow(scr[1], 2) * scr[2];
218+
break;
219+
case 16:
220+
/* m=+3, xyy - xxx */
221+
fcmpnt = scr[0] * pow(scr[1], 2);
222+
fcmpnt -= pow(scr[0], 3);
223+
break;
224+
default:
225+
error_message(0, "Unknown type!", "evaluate_f");
226+
printf(" type = %d\n", typ);
227+
return value;
228+
}
229+
for (i = 0; i < ng; i++) {
230+
nrmcnst = compute_f_normconst(alphac[i]);
231+
value = value + (gscale * ccoef[i] * nrmcnst *
232+
fcmpnt * gauss_fcn(alphac[i], rval));
233+
}
234+
return value;
235+
}
236+
237+
/*
238+
* evaluate_p: evaluate a contracted p orbital at the position (x, y, z).
239+
*
240+
* Input:
241+
* pos = x, y, z to evalute orbital at
242+
* alphac = alpha coefficients for contracted gaussians
243+
* ccoef = contraction coefficients
244+
* ng = number of gaussians in contraction
245+
* gscale = gaussian scaling coefficient
246+
* atompos = x, y, z of nuclear center
247+
* typ = px=2, py=3, pz=4
248+
* scr = scratch array to hold |R - r|
249+
*/
250+
double evaluate_p(double *pos, double *alphac, double *ccoef, int ng,
251+
double gscale, double *atompos, int typ, double *scr)
252+
{
253+
double value = 0.0; /* Return value */
254+
double rval;
255+
double pcmpnt; /* p_{type} component {type} (x, y, or z) */
256+
double nrmcnst;
257+
int i;
258+
/* Compute {X-x,Y-y,Z-z} and ||{X-x,Y-y,Z-z}||*/
259+
vector_difference(pos, atompos, 3, scr);
260+
rval = compute_vector_norm(scr, 3);
261+
switch (typ) {
262+
case 2:
263+
pcmpnt = scr[0];
264+
break;
265+
case 3:
266+
pcmpnt = scr[1];
267+
break;
268+
case 4:
269+
pcmpnt = scr[2];
270+
break;
271+
default:
272+
error_message(0, "Unknown type!", "evaluate_p");
273+
printf(" type = %d\n", typ);
274+
return value;
275+
}
276+
for (i = 0; i < ng; i++) {
277+
nrmcnst = compute_p_normconst(alphac[i]);
278+
value = value + (gscale * ccoef[i] * nrmcnst *
279+
pcmpnt * gauss_fcn(alphac[i], rval));
280+
}
281+
return value;
282+
}
283+
284+
/*
285+
* evaluate_s: evaluate a contracted s orbital at the position (x, y, z).
286+
* Input:
287+
* pos = x, y, z to evalute orbital at
288+
* alphac = alpha coefficients for contracted gaussians
289+
* ccoef = contraction coefficients
290+
* ng = number of gaussians in contraction
291+
* gscale = gaussian scaling coefficient
292+
* atompos = x, y, z of nuclear center
293+
*/
294+
double evaluate_s(double *pos, double *alphac, double *ccoef, int ng,
295+
double gscale, double *atompos)
296+
{
297+
double value = 0.0; /* Return value */
298+
double rval; /* Sqrt(pos - atompos) */
299+
double nrmcnst;
300+
int i;
301+
rval = compute_3d_distance(atompos, pos);
302+
printf("%lf %lf %lf\n", pos[0], pos[1], pos[2]);
303+
printf("%lf %lf %lf\n", atompos[0], atompos[1], atompos[2]);
304+
printf("rval = %lf\n", rval);
305+
for (i = 0; i < ng; i++) {
306+
nrmcnst = compute_s_normconst(alphac[i]);
307+
value = value + (gscale * ccoef[i] *
308+
nrmcnst * gauss_fcn(alphac[i], rval));
309+
}
310+
return value;
311+
}
312+
313+

0 commit comments

Comments
 (0)