← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccresponse/linresp.c

  • 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
http://bugs.debian.org/464867
* 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:
src/bin/ccresponse/linresp.c
 
1
#include <stdio.h>
 
2
#include <string.h>
 
3
#include <stdlib.h>
 
4
#include <libciomr/libciomr.h>
 
5
#include <libpsio/psio.h>
 
6
#define EXTERN
 
7
#include "globals.h"
 
8
 
 
9
double LCX(char *pert_c, char *cart_c, int irrep_c, 
 
10
           char *pert_x, char *cart_x, int irrep_x, double omega);
 
11
double HXY(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
 
12
           char *pert_y, char *cart_y, int irrep_y, double omega_y);
 
13
double LHX1Y1(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
 
14
              char *pert_y, char *cart_y, int irrep_y, double omega_y);
 
15
double LHX2Y2(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
 
16
              char *pert_y, char *cart_y, int irrep_y, double omega_y);
 
17
double LHX1Y2(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
 
18
              char *pert_y, char *cart_y, int irrep_y, double omega_y);
 
19
double cc2_LHX1Y1(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
 
20
                  char *pert_y, char *cart_y, int irrep_y, double omega_y);
 
21
double cc2_LHX1Y2(char *pert_x, char *cart_x, int irrep_x, double omega_x, 
 
22
                  char *pert_y, char *cart_y, int irrep_y, double omega_y);
 
23
 
 
24
void linresp(double **tensor, double A, double B,
 
25
             char *pert_x, int *x_irreps, double omega_x,
 
26
             char *pert_y, int *y_irreps, double omega_y)
 
27
{
 
28
  int alpha, beta, j;
 
29
  double polar, polar_LCX, polar_HXY, polar_LHX1Y1, polar_LHX1Y2, polar_LHX2Y2;
 
30
  char **cartcomp;
 
31
 
 
32
  cartcomp = (char **) malloc(3 * sizeof(char *));
 
33
  cartcomp[0] = strdup("X");
 
34
  cartcomp[1] = strdup("Y");
 
35
  cartcomp[2] = strdup("Z");
 
36
 
 
37
  for(alpha=0; alpha < 3; alpha++) {
 
38
    for(beta=0; beta < 3; beta++) {
 
39
 
 
40
      /* clear out scratch space */
 
41
      for(j=CC_TMP; j <= CC_TMP11; j++) {
 
42
         psio_close(j,0); psio_open(j,0);
 
43
      }
 
44
 
 
45
      polar_LCX = 0.0;
 
46
      polar_HXY = 0.0;
 
47
      polar_LHX1Y1 = 0.0;
 
48
      polar_LHX2Y2 = 0.0;
 
49
      polar_LHX1Y2 = 0.0;
 
50
 
 
51
      if(!(x_irreps[alpha]^y_irreps[beta])) {
 
52
 
 
53
        if(omega_y != 0.0) {  /* we assume omega_x = -omega_y */
 
54
          polar_LCX = LCX(pert_x, cartcomp[alpha], x_irreps[alpha], 
 
55
                          pert_y, cartcomp[beta], y_irreps[beta], omega_y);
 
56
          polar_LCX += LCX(pert_y, cartcomp[beta], y_irreps[beta], 
 
57
                           pert_x, cartcomp[alpha], x_irreps[alpha], omega_x);
 
58
 
 
59
          if(!params.sekino) {
 
60
            if (!strcmp(params.wfn,"CC2")) {
 
61
              polar_HXY = HXY(pert_x, cartcomp[alpha], x_irreps[alpha], omega_x,
 
62
                              pert_y, cartcomp[beta], y_irreps[beta], omega_y);
 
63
              polar_LHX1Y1 = cc2_LHX1Y1(pert_x, cartcomp[alpha], x_irreps[alpha], omega_x,
 
64
                                        pert_y, cartcomp[beta], y_irreps[beta], omega_y);
 
65
              polar_LHX1Y2 = cc2_LHX1Y2(pert_x, cartcomp[alpha], x_irreps[alpha], omega_x,
 
66
                                        pert_y, cartcomp[beta], y_irreps[beta], omega_y);
 
67
              polar_LHX1Y2 += cc2_LHX1Y2(pert_y, cartcomp[beta], y_irreps[beta], omega_y,
 
68
                                         pert_x, cartcomp[alpha], x_irreps[alpha], omega_x);
 
69
            }
 
70
            else {
 
71
              polar_LHX1Y1 = LHX1Y1(pert_x, cartcomp[alpha], x_irreps[alpha], omega_x,
 
72
                                    pert_y, cartcomp[beta], y_irreps[beta], omega_y);
 
73
              polar_LHX2Y2 = LHX2Y2(pert_x, cartcomp[alpha], x_irreps[alpha], omega_x,
 
74
                                    pert_y, cartcomp[beta], y_irreps[beta], omega_y);
 
75
              polar_LHX1Y2 = LHX1Y2(pert_x, cartcomp[alpha], x_irreps[alpha], omega_x,
 
76
                                    pert_y, cartcomp[beta], y_irreps[beta], omega_y);
 
77
              polar_LHX1Y2 += LHX1Y2(pert_y, cartcomp[beta], y_irreps[beta], omega_y,
 
78
                                     pert_x, cartcomp[alpha], x_irreps[alpha], omega_x);
 
79
            }
 
80
          }
 
81
        }
 
82
        else {
 
83
          polar_LCX = LCX(pert_x, cartcomp[alpha], x_irreps[alpha], pert_y, cartcomp[beta],
 
84
                          y_irreps[beta], 0.0);
 
85
          polar_LCX += LCX(pert_y, cartcomp[beta], y_irreps[beta], pert_x, cartcomp[alpha],
 
86
                           x_irreps[alpha], 0.0);
 
87
          if(!params.sekino) {
 
88
            if (!strcmp(params.wfn,"CC2")) {
 
89
              polar_HXY = HXY(pert_x, cartcomp[alpha], x_irreps[alpha], 0.0,
 
90
                              pert_y, cartcomp[beta], y_irreps[beta], 0.0);
 
91
              polar_LHX1Y1 = cc2_LHX1Y1(pert_x, cartcomp[alpha], x_irreps[alpha], 0.0,
 
92
                                        pert_y, cartcomp[beta], y_irreps[beta], 0.0);
 
93
              polar_LHX1Y2 = cc2_LHX1Y2(pert_x, cartcomp[alpha], x_irreps[alpha], 0.0,
 
94
                                        pert_y, cartcomp[beta], y_irreps[beta], 0.0);
 
95
              polar_LHX1Y2 += cc2_LHX1Y2(pert_y, cartcomp[beta], y_irreps[beta], 0.0,
 
96
                                         pert_x, cartcomp[alpha], x_irreps[alpha], 0.0);
 
97
            }
 
98
            else {
 
99
              polar_LHX1Y1 = LHX1Y1(pert_x, cartcomp[alpha], x_irreps[alpha], 0.0,
 
100
                                    pert_y, cartcomp[beta], y_irreps[beta], 0.0);
 
101
              polar_LHX2Y2 = LHX2Y2(pert_x, cartcomp[alpha], x_irreps[alpha], 0.0,
 
102
                                    pert_y, cartcomp[beta], y_irreps[beta], 0.0);
 
103
              polar_LHX1Y2 = LHX1Y2(pert_x, cartcomp[alpha], x_irreps[alpha], 0.0,
 
104
                                    pert_y, cartcomp[beta], y_irreps[beta], 0.0);
 
105
              polar_LHX1Y2 += LHX1Y2(pert_y, cartcomp[beta], y_irreps[beta], 0.0,
 
106
                                     pert_x, cartcomp[alpha], x_irreps[alpha], 0.0);
 
107
            }
 
108
          }
 
109
        }
 
110
 
 
111
        if(params.sekino) /* only linear term needed in Sekino-Bartlett model III */
 
112
          polar = polar_LCX;
 
113
        else 
 
114
          polar = polar_LCX + polar_HXY + polar_LHX1Y1 + polar_LHX2Y2 + polar_LHX1Y2;
 
115
 
 
116
        if(params.print & 2) {
 
117
          fprintf(outfile, "\n\tLinresp tensor[%s][%s]\n", cartcomp[alpha], cartcomp[beta]);
 
118
          fprintf(outfile, "\tpolar_LCX    = %20.12f\n", polar_LCX);
 
119
          if(!strcmp(params.wfn,"CC2"))
 
120
            fprintf(outfile, "\tpolar_HXY    = %20.12f\n", polar_HXY);
 
121
          fprintf(outfile, "\tpolar_LHX1Y1 = %20.12f\n", polar_LHX1Y1);
 
122
          fprintf(outfile, "\tpolar_LHX1Y2 = %20.12f\n", polar_LHX1Y2);
 
123
          fprintf(outfile, "\tpolar_LHX2Y2 = %20.12f\n", polar_LHX2Y2);
 
124
          fflush(outfile);
 
125
        }
 
126
 
 
127
        tensor[alpha][beta] = A * polar + B * tensor[alpha][beta];
 
128
      }
 
129
    }
 
130
  }
 
131
 
 
132
  free(cartcomp[0]);
 
133
  free(cartcomp[1]);
 
134
  free(cartcomp[2]);
 
135
  free(cartcomp);
 
136
}