1 /*----------------------------------------------------------------------
2   erl.h
3 
4   16 Feb. 2009, M. Toyoda
5 ----------------------------------------------------------------------*/
6 #ifndef LIBERI_ERI_H_INCLUDED
7 #define LIBERI_ERI_H_INCLUDED
8 
9 #include <stdlib.h>
10 #include <float.h>
11 
12 
13 #define ERI_SBT_LOG    1 /* Siegman-Talman method */
14 #define ERI_SBT_LINEAR 2 /* Toyoda-Ozaki method */
15 #define ERI_NOSCREEN -1.0
16 
17 #define ERI_SH_REAL    1  /* Real-valued SH is used */
18 #define ERI_SH_COMPLEX 2  /* Complex-valued SH is used */
19 
20 
21 #if defined(__cplusplus)
22 extern "C" {
23 #endif
24 
25 
26 /*----------------------------------------------------------------------
27   Initializer and Finalizer
28 ----------------------------------------------------------------------*/
29 typedef struct ERI_Struct ERI_t;
30 
31 const char* ERI_Version(void);
32 
33 
34 typedef struct {
35   int    sbttype;
36   double rmax;
37   double rho0;
38   double qmin;
39   double qmax;
40   int    nq;
41 } ERI_Init_Misc;
42 
43 
44 /*----------------------------------------------------------------------
45   ERI_Init
46 
47   Initializes LIBERI.
48 
49   IN:
50     lmax    : maximum angular momentum
51     lmax_gl : maximum angular momenrum for the final summation
52     ngrid   : number of radial mesh points
53     ngl     : number of abscissas for GL quadrature
54     sbttyp  : type of SBT routine (ERI_SBT_LOG or ERI_SBT_LINEAR)
55     rmax    : range of r-mesh
56     info    : additional information
57               (this is optional. if this is NULL, the default settings
58                are used).
59     rcsh    : real- or complex-valued SH is used.
60   RETURN:
61     pointer to an ERI_Struct object.
62 ----------------------------------------------------------------------*/
63 ERI_t* ERI_Init(
64   int lmax,
65   int lmax_gl,
66   int ngrid,
67   int ngl,
68   int rcsh,
69   const ERI_Init_Misc* info
70 );
71 
72 
73 /*----------------------------------------------------------------------
74   ERI_Required_Size
75 
76   Estimates the byte size which is used internally by LIBERI.
77 
78   See ERI_Init for the explantation of the arguments.
79 ----------------------------------------------------------------------*/
80 size_t ERI_Required_Size(
81   int lmax,
82   int lmax_gl,
83   int ngrid,
84   int ngl,
85   int rcsh,
86   const ERI_Init_Misc* info
87 );
88 
89 
90 /*----------------------------------------------------------------------
91   ERI_Free
92 
93   Finalizes LIBERI.
94 
95   IN:
96     ptr : pointer returned by ERI_Init.
97 ----------------------------------------------------------------------*/
98 void ERI_Free(ERI_t *ptr);
99 
100 
101 /*----------------------------------------------------------------------
102   ERI_lmax, ERI_lmax_gl, ERI_ngird, ERI_ngl, ERI_Misc
103 
104   Returns the initialization parameters.
105 ----------------------------------------------------------------------*/
106 int                  ERI_lmax    (const ERI_t *ptr);
107 int                  ERI_lmax_gl (const ERI_t *ptr);
108 int                  ERI_ngl     (const ERI_t *ptr);
109 const ERI_Init_Misc* ERI_Misc    (const ERI_t *ptr);
110 int                  ERI_ngrid   (const ERI_t *ptr);
111 int                  ERI_rcsh    (const ERI_t *ptr);
112 
113 /*----------------------------------------------------------------------
114   ERI_Mesh_r, ERI_Mesh_k, ERI_Mesh_dr, ERI_Mesh_dk
115 
116   Returns the radial mesh points or the intervals of the mesh.
117 ----------------------------------------------------------------------*/
118 double        ERI_Mesh_r       (const ERI_t *ptr, int i);
119 double        ERI_Mesh_k       (const ERI_t *ptr, int i);
120 double        ERI_Mesh_dr      (const ERI_t *ptr, int i);
121 double        ERI_Mesh_dk      (const ERI_t *ptr, int i);
122 
123 
124 /*----------------------------------------------------------------------
125   ERI_Mesh_Array_r, ERI_Mesh_Array_k
126 
127   Returns the all radial mesh points.
128 ----------------------------------------------------------------------*/
129 const double* ERI_Mesh_Array_r (const ERI_t *ptr);
130 const double* ERI_Mesh_Array_k (const ERI_t *ptr);
131 
132 
133 /*----------------------------------------------------------------------
134   ERI_Mesh_Array_glx, ERI_Mesh_Array_glw
135 
136   Returns the abscissas and the weight factors for the GL quadrature.
137 ----------------------------------------------------------------------*/
138 const double* ERI_Mesh_Array_glx(const ERI_t *ptr);
139 const double* ERI_Mesh_Array_glw(const ERI_t *ptr);
140 
141 
142 
143 
144 /*----------------------------------------------------------------------
145   High-level `black box' functions
146 ----------------------------------------------------------------------*/
147 
148 typedef struct {
149   double       *fr;   /* radial function */
150   double       *xr;   /* radial grid */
151   unsigned int  ngrid; /* number of radial grid */
152   int           l;     /* angular momentum */
153   int           m;     /* angular momentum */
154   double        c[3];  /* position */
155 } ERI_Orbital_t;
156 
157 
158 
159 /*----------------------------------------------------------------------
160   ERI_Overlap
161 
162   Calculated the transformed overlap function for given orbital
163   functions.
164 
165   IN:
166     ptr  : pointer returned by ERI_Init
167     orb1 : an orbital function
168     orb2 : another orbital fuction
169     cx   : parameter for expansion center
170 
171   OUT:
172     F    : claculated overlap function
173     dF   : derivatives of F
174            (this is optional. if you do not need the derivatives,
175             specify NULL here.)
176 ----------------------------------------------------------------------*/
177 void ERI_Overlap(
178   ERI_t  *ptr,
179   double *F,
180   double *dF[3],
181   const ERI_Orbital_t *orb1,
182   const ERI_Orbital_t *orb2,
183   double cx
184 );
185 
186 
187 
188 /*----------------------------------------------------------------------
189   ERI_Integral
190 
191   Calculates ERI for given four orbital functions.
192 
193   IN:
194     ptr  : pointer returned by ERI_Init.
195     orb1 : orbital function #1
196     orb2 : orbital function #2
197     orb3 : orbital function #3
198     orb4 : orbital function #4
199 
200   OUT:
201     I4   : calculated ERI (complex number)
202     dI4  : derivatives of ERI
203            (this is optional. if you do not need the derivatives,
204             specify NULL here.)
205 ----------------------------------------------------------------------*/
206 void ERI_Integral(
207   ERI_t  *ptr,
208   double I4[2],
209   double *dI4,
210   const ERI_Orbital_t *orb1,
211   const ERI_Orbital_t *orb2,
212   const ERI_Orbital_t *orb3,
213   const ERI_Orbital_t *orb4,
214   double               scr
215 );
216 
217 
218 void ERI_Coordinate_Transform(
219   double out_R[3],  /* (OUT) displacement in spherical coord. */
220   double out_a1[3], /* (OUT) translated center 1 in spherical coord. */
221   double out_a2[3], /* (OUT) translated center 2 in spherical coord. */
222   double out_a3[3], /* (OUT) translated center 3 in spherical coord. */
223   double out_a4[3], /* (OUT) translated center 4 in spherical coord. */
224   const double in_a1[3],  /* (IN) center 1 in Cartesian coord. */
225   const double in_a2[3],  /* (IN) center 2 in Cartesian coord. */
226   const double in_a3[3],  /* (IN) center 3 in Cartesian coord. */
227   const double in_a4[3],  /* (IN) center 4 in Cartesian coord. */
228   double x12,   /* (IN) ratio */
229   double x34    /* (IN) ratio */
230 );
231 
232 
233 
234 /*----------------------------------------------------------------------
235   ERI_Center_r2
236 
237   IN:
238     orb1 : orbital #1
239     orb2 : orbital #2
240 
241   RETURN:
242     ratio for the expansion center.
243 ----------------------------------------------------------------------*/
244 double ERI_Center_r2(
245   const ERI_Orbital_t *orb1,
246   const ERI_Orbital_t *orb2
247 );
248 
249 
250 /*----------------------------------------------------------------------
251   ERI_Center_DK
252 
253 
254   IN:
255     ptr  : pointer returned by ERI_Init
256     orb1 : orbital #1
257     orb1 : orbital #2
258 
259   RETURN:
260     ratio for the expansion center.
261 ----------------------------------------------------------------------*/
262 double ERI_Center_DK(
263   ERI_t        *ptr,
264   const ERI_Orbital_t *orb1,
265   const ERI_Orbital_t *orb2
266 );
267 
268 
269 
270 
271 /*----------------------------------------------------------------------
272   Low-level functions
273 ----------------------------------------------------------------------*/
274 
275 /*----------------------------------------------------------------------
276   ERI_Size_of_*
277 
278   Returns required byte size for the matrices.
279 ----------------------------------------------------------------------*/
280 size_t ERI_Size_of_Orbital (const ERI_t *ptr);
281 size_t ERI_Size_of_Gamma   (const ERI_t *ptr);
282 size_t ERI_Size_of_Alpha   (const ERI_t *ptr);
283 size_t ERI_Size_of_Overlap (const ERI_t *ptr);
284 size_t ERI_Size_of_GLF     (const ERI_t *ptr);
285 
286 
287 /*----------------------------------------------------------------------
288   ERI_Transform_Orbital
289 
290   Given input function, this performs the SBT.
291 
292   Note that the radial mesh for input function (fr) is NOT the r-mesh
293   defined by LIBERI. It is the user-defined mesh (xr). The user-defined
294   radial mesh (xr) can be of any type, as long as it is in ascenfing
295   order.
296 
297   The k-mesh for the output function (fk) is defined by LIBERI.
298 
299   IN:
300     ptr : pointer returned by ERI_Init
301     fr  : values of input function at each
302     xr  : points of the radial mesh for fr.
303     n   : number of the radial mesh for fr.
304     l   : angular momenrum
305 
306   OUT:
307     fk  : transformed orbital
308           (The size of the array should be no smaller then the size
309            returned by ERI_Size_of_Orbital.)
310 ----------------------------------------------------------------------*/
311 void ERI_Transform_Orbital(
312   ERI_t        *ptr,
313   double       *fk,  /* (OUT) transformed orbital */
314   const double *fr,  /* (IN) radial function */
315   const double *xr,  /* (IN) radial mesh */
316   int           n,   /* (IN) number of mesh poitns */
317   int           l    /* (IN) angular momentum */
318 );
319 
320 
321 #if 0
322 /*----------------------------------------------------------------------
323   ERI_LL_Gamma
324 
325   Calculates the gamma-functions.
326 
327   IN:
328     ptr : pointer returned by ERI_Init
329     fk  : transformed function (returned by ERI_Transform_Orbital)
330     rd  : displacement from the center
331 
332   OUT:
333     g   : gamma-functions.
334           (The size of the array should be no smaller than the size
335            returned by ERI_Size_of_Gamma.)
336 ----------------------------------------------------------------------*/
337 void ERI_LL_Gamma(
338   ERI_t        *ptr,
339   double       *g,   /* (OUT) Gamma term */
340   const double *fk,  /* (IN)  Transformed radial function */
341   double        rd   /* (IN)  Displacemenent from the center */
342 );
343 #endif
344 
345 
346 /*----------------------------------------------------------------------
347   ERI_LL_Alpha
348 
349   Calculates the alpha-functions.
350 
351   IN:
352     ptr        : pointer returned by ERI_Init
353     g          : gamma-functions (returned by ERI_LL_Gamma)
354     theta, phi : angle displacement from the center
355     m          : magnetic angular momentum of the orbital
356 
357   OIT:
358     alp        : alpha-functions
359                  (The size of the array shuold be no smaller than the
360                   size returned by ERI_Size_of_Alpha.)
361 ----------------------------------------------------------------------*/
362 void ERI_LL_Alpha(
363   ERI_t        *ptr,
364   double       *alp,   /* (OUT) Alpha terms */
365   const double *g,     /* (IN)  Gamma terms */
366   double        theta,
367   double        phi,
368   int           l,     /* (IN)  Angular momentum of the orbital */
369   int           m
370 );
371 
372 
373 /*----------------------------------------------------------------------
374   ERI_LL_Overlap
375 
376   Calculatets the overlap function.
377 
378   IN:
379     ptr : pointer returned by ERI_Init
380     a1  : alpha-function of orbital 1
381     a2  : alpha-function of orbital 2
382 
383   OUT:
384     p   : overlap function
385           (The size of the array should be no smaller than the size
386            returned by ERI_Size_of_Overlap.)
387 ----------------------------------------------------------------------*/
388 void ERI_LL_Overlap(
389   ERI_t        *ptr,
390   double       *p,
391   const double *a1,
392   const double *a2
393 );
394 
395 
396 /*----------------------------------------------------------------------
397   ERI_Transform_Overlap
398 
399   Performes SBT of the overlap functions.
400 
401   If you do not need the derivatives, specify NULL for dp and df.
402 
403   IN:
404     ptr : pointer returned by ERI_Init
405     p   : overlap function
406     dp  : derivatives of overlap function (optional)
407 
408   OUT:
409     f   : transformed overlap functions
410     df  : derivatives wrt nucleus position in x, y and z direction,
411           respectively (optional).
412           (The size of each array should be no smaller than the size
413            returned by ERI_Size_of_Overlap.)
414 ----------------------------------------------------------------------*/
415 void ERI_Transform_Overlap(
416   ERI_t        *ptr,
417   double       *f,
418   const double *p
419 );
420 #if 0
421 void ERI_Transform_Overlap(
422   ERI_t        *ptr,
423   double       *f,
424   double       *df[3],
425   const double *p,
426   const double *dp[3]
427 );
428 #endif
429 
430 /*----------------------------------------------------------------------
431   ERI_GL_Interpolate
432 
433   Interpolates the transformed overlap function from SBT-mesh to
434   GL-mesh.
435 
436   IN:
437     ptr : poitner returned by ERI_Init
438     f   : overlap function on SBT-mesh
439 
440   OUT:
441     glf : overlap function on GL-mesh
442 ----------------------------------------------------------------------*/
443 void ERI_GL_Interpolate(
444   ERI_t        *ptr,
445   double       *glf,
446   const double *f
447 );
448 
449 
450 /*----------------------------------------------------------------------
451   ERI_LL_Gamma
452 
453   Calculates the gamma functions and the derivatives.
454 
455   IN:
456     ptr : poitner returned by ERI_Init
457     fk  : transformed function returned by ERI_Transform_Orbital.
458     rd  : displacement from the center
459 
460   OUT:
461     g   : gamma-functions
462     dg  : derivatives of the gamma-functions
463           (The size of each array should be no smaller than the size
464            returned by ERI_Size_of_Gamma)
465 ----------------------------------------------------------------------*/
466 void ERI_LL_Gamma(
467   ERI_t        *ptr,
468   double       *g,   /* (OUT) Gamma terms */
469   double       *dg,  /* (OUT) Derivatives of the Gamma terms */
470   const double *fk,  /* (IN) Transformed radidal function */
471   double        rd   /* (IN) Displacement from the center */
472 );
473 
474 
475 /*----------------------------------------------------------------------
476   ERI_LL_Alpha_d
477 
478   Calculates the alpha-funtions and the derivatives.
479 
480   IN:
481     ptr            : pointer returned by ERI_Init
482     g              : gamma-functions
483     dg             : derivatives of gamma-functions
484     rd, theta, phi : displacement from the center in spherical coord.
485     l, m           : angular momentum or the orbital.
486 
487   OUT:
488     alp  : alpha-functions
489     dalp : derivatives of alpha-functions wrt nuclear position
490            in x, y and z directions, respectively.
491            (The size of each array should be no smaller than the size
492             returned by ERI_Size_of_Alpha.)
493 ----------------------------------------------------------------------*/
494 void ERI_LL_Alpha_d(
495   ERI_t        *ptr,
496   double       *alp,
497   double       *dalp[3],
498   const double *g,
499   const double *dg,
500   double        rd,
501   double        theta,
502   double        phi,
503   int           l,
504   int           m
505 );
506 
507 
508 /*----------------------------------------------------------------------
509   ERI_LL_Overlap_d
510 
511   Calculates the overlap functions and the derivatives.
512 
513   IN:
514     ptr : pointer returned by ERI_Init
515     a1  : alpha-function of orbtial 1
516     da1 : derivatives of alpha-function of orbital 1
517     a2  : alpha-function of orbtial 2
518     da2 : derivatives of alpha-function of orbital 2
519     cx  : ratio for expansion center
520 
521   OUT:
522     p   : overlap function
523     dp  : derivatives of overlap function
524           (The size of each array should be no smaller than the size
525            returned by ERI_Size_of_Overlap.)
526 ----------------------------------------------------------------------*/
527 void ERI_LL_Overlap_d(
528   ERI_t        *ptr,
529   double       *p,
530   double       *dp[3],
531   const double *a1,
532   const double *da1[3],
533   const double *a2,
534   const double *da2[3],
535   double        x
536 );
537 
538 
539 
540 
541 /*----------------------------------------------------------------------
542   ERI_Integral_GL
543 
544   Calculates ERI for given overlap functions.
545 
546   IN:
547     ptr : pointer returned by ERI_Init
548     p12 : transformed overlap function for orbiatal 1 and 2
549     p34 : transformed overlap function for orbiatal 3 and 4
550     rd, theta, phi : displacement of expansion centers in sph. coord.
551     omega          : screening factor
552                      (speficy ERI_NOSCREEN if you do not need it)
553     lmax           : cut-off of the angular momentum summation
554 
555   OUT:
556     I4  : ERI (complex number)
557 ----------------------------------------------------------------------*/
558 void ERI_Integral_GL(
559   ERI_t        *ptr,
560   double        I4[2],
561   const double *p12,
562   const double *p34,
563   double        rd,
564   double        theta,
565   double        phi,
566   double        omega,
567   int           lmax1
568 );
569 
570 
571 void ERI_Integral_GL_d(
572   ERI_t        *ptr,
573   double        I4[2],  /* (OUT) */
574   double        dI4[4][3][2],
575   const double *F1,     /* (IN) Overlap matrix */
576   const double *F2,     /* (IN) */
577   const double *dF1[3],    /* (IN) Overlap matrix */
578   const double *dF2[3],     /* (IN) */
579   double        R,      /* (IN) Displacement of two expansion centers */
580   double        theta,
581   double        phi,
582   double        cx12,
583   double        cx34,
584   double        delta,
585   double        omega,  /* (IN) screening parameter */
586   int           lmax1
587 );
588 
589 
590 
591 #if 0
592 void ERI_Integral_GL_PrejY(
593   ERI_t        *solver,
594   const double *F1,    /* (IN) Overlap matrix */
595   const double *F2,    /* (IN) */
596   double       *R,     /* (IN) Displacement of two expansion centers */
597   double       *theta,
598   double       *phi,
599   int           numR,
600   double        omega,  /* (IN) screening parameter */
601   int          *mul_j2, /* [jmax1*jmax1*jmax1] */
602   double       *mul_gc, /* [jmax1*jmax1*jmax1] */
603   int          *mul_n,  /* [jmax1*jmax1] */
604   double       *prej,   /* [lmax*ngl*numR] */
605   double       *preY,   /* [numR*jmax1] */
606   int          *num_minimalR,
607   int          *normalR,      /* [numR*numR] */
608   int          *num_normalR    /* [numR] */
609 );
610 
611 
612 
613 
614 void ERI_Integral_GL_Post(
615   ERI_t        *solver,
616   double       *I4,    /* (OUT) [numR] */
617   const double *F1,    /* (IN) Overlap matrix */
618   const double *F2,    /* (IN) */
619   int           numR,
620   int          *mul_j2, /* [jmax1*jmax1*jmax1] */
621   double       *mul_gc, /* [jmax1*jmax1*jmax1] */
622   int          *mul_n,  /* [jmax1*jmax1] */
623   double       *prej,   /* [lmax*ngl*numR] */
624   double       *preY,    /* [numR*jmax1] */
625   int           num_minimalR,
626   int          *normalR,      /* [numR*numR] */
627   int          *num_normalR    /* [numR] */
628 );
629 #endif
630 
631 
632 void ERI_Integral_GL_PrejY(
633   ERI_t        *solver,
634   double       *R,     /* (IN) Displacement of two expansion centers */
635   double       *theta,
636   double       *phi,
637   int           numR,
638   double        omega,  /* (IN) screening parameter */
639   double       *prej,   /* [lmax*ngl*numR] */
640   double       *preY,   /* [numR*jmax1] */
641   int          *mul_j2, /* [jmax1*jmax1*jmax1] */
642   double       *mul_gc, /* [jmax1*jmax1*jmax1] */
643   int          *mul_n,  /* [jmax1*jmax1] */
644   int          *minimalR,     /* [numR] */
645   int          *num_minimalR
646 );
647 
648 
649 void ERI_Integral_GL_Post(
650   ERI_t        *solver,
651   double       *I4,    /* (OUT) [numR] */
652   const double *F1,    /* (IN) Overlap matrix */
653   const double *F2,    /* (IN) */
654   /*const double *Q1,*/    /* (IN) Overlap matrix */
655   /*const double *Q2,*/    /* (IN) */
656   int           numR,
657   double       *prej,   /* [lmax*ngl*numR] */
658   double       *preY,    /* [numR*jmax1] */
659   int          *mul_j2, /* [jmax1*jmax1*jmax1] */
660   double       *mul_gc, /* [jmax1*jmax1*jmax1] */
661   int          *mul_n,  /* [jmax1*jmax1] */
662   int          *minimalR,     /* [numR] */
663   int           num_minimalR
664 );
665 
666 void ERI_Integral_GL_Post2(
667   ERI_t        *solver,
668   double       *I4,    /* (OUT) [numR] */
669   const double *F1,    /* (IN) Overlap matrix */
670   const double *F2,    /* (IN) */
671   /*const double *Q1,*/    /* (IN) Overlap matrix */
672   /*const double *Q2,*/    /* (IN) */
673   int           numR,
674   double       *prej,   /* [lmax*ngl*numR] */
675   double       *preY,    /* [numR*jmax1] */
676   int          *mul_j2, /* [jmax1*jmax1*jmax1] */
677   double       *mul_gc, /* [jmax1*jmax1*jmax1] */
678   int          *mul_n,  /* [jmax1*jmax1] */
679   int          *minimalR,     /* [numR] */
680   int           num_minimalR
681 );
682 
683 
684 
685 void ERI_Integral_GL_X(
686   ERI_t        *solver,
687   double       *X4,    /* (OUT) [numR] */
688   const double *F1,    /* (IN) Overlap matrix */
689   const double *F2,    /* (IN) */
690   int           numR,
691   double       *prej,   /* [lmax*ngl*numR] */
692   int          *mul_j2, /* [jmax1*jmax1*jmax1] */
693   double       *mul_gc, /* [jmax1*jmax1*jmax1] */
694   int          *mul_n,  /* [jmax1*jmax1] */
695   int          *minimalR,    /* [numR] */
696   int           num_minimalR
697 );
698 
699 void ERI_Integral_GL_X_Post(
700   ERI_t        *solver,
701   double       *I4,       /* (OUT) [numR] */
702   const double *X,        /* (IN)  [numR] */
703   int           numR,
704   const double *preY,     /* (IN)  [numR*jmax1] */
705   const int    *minimalR  /* (IN)  [numR] */
706 );
707 
708 
709 #if defined(__cplusplus)
710 }
711 #endif
712 
713 #endif /* LIBERI_ERI_H_INCLUDED */
714 /* EOF */
715