← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/detcas/gradient.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/detcas/gradient.c
1
 
#include <stdlib.h>
2
 
#include <stdio.h>
3
 
#include <math.h>
4
 
#include <libipv1/ip_lib.h>
5
 
#include <libciomr/libciomr.h>
6
 
#include <libqt/qt.h>
7
 
#include "globaldefs.h"
8
 
#include "globals.h"
9
 
 
10
 
 
11
 
/*
12
 
** calc_grad_1
13
 
**
14
 
** Calculates the orbital gradient from the MO Lagrangian read in
15
 
** from the CLAG program
16
 
**
17
 
** An alternative method is available using calc_grad_2()
18
 
**/
19
 
void calc_grad_1(int npair, int *ppair, int *qpair, double **lag, double *grad)
20
 
{
21
 
  int pair, p, q;
22
 
  double value;
23
 
 
24
 
  /* this def of the gradient is consistent with orbital rotations defined
25
 
   * in the same sense as our VBD paper...also, what we would call the
26
 
   * Lagrangian in the VBD paper is actually twice the Lagrangian computed
27
 
   * by CLAG, so double the contribution.
28
 
   */
29
 
  for (pair=0; pair<npair; pair++) {
30
 
    p = *ppair++;
31
 
    q = *qpair++;
32
 
    value = 2.0 * (lag[q][p] - lag[p][q]);
33
 
    grad[pair] = value;
34
 
  }
35
 
 
36
 
}
37
 
 
38
 
 
39
 
 
40
 
 
41
 
/*
42
 
** calc_grad_2
43
 
**
44
 
** Calculates the orbital gradient from the F_core and F_act quantities,
45
 
** not from CLAG
46
 
**
47
 
** I am assuming that the pairs (p,q) are always given such that p>=q
48
 
**
49
 
*/
50
 
void calc_grad_2(int npairs, int *ppair, int *qpair, double *F_core, 
51
 
                 double *tei, double **opdm, double *tpdm, double *F_act, 
52
 
                 int firstact, int lastact, double *grad)
53
 
{
54
 
 
55
 
  int pair, p, q, pq, pp, qq;
56
 
  int i,ii,a,aa,t,tt,ti,u,tu,au,iu,v,w,vw,tuvw,auvw,iuvw;
57
 
  double value;
58
 
 
59
 
  /* this def of the gradient is consistent with orbital rotations defined
60
 
   * in the same sense as our VBD paper...we have to reverse the sign
61
 
   * on all their terms
62
 
   */
63
 
 
64
 
  /* loop over the independent pairs */
65
 
  for (pair=0; pair<npairs; pair++) {
66
 
    p = ppair[pair];
67
 
    q = qpair[pair];
68
 
    pq = ioff[p] + q;
69
 
    pp = ioff[p] + p;
70
 
    qq = ioff[q] + q;
71
 
  
72
 
    /* g_{ai}, i.e., inactive virt/inactive occ */
73
 
    if (p >= lastact && q < firstact) {
74
 
      grad[pair] = -4.0 * (F_core[pq] + F_act[pq]);
75
 
    }
76
 
 
77
 
    /* g_{at}, i.e., inactive virt with active orb */
78
 
    else if (p >= lastact && q >= firstact) {
79
 
      a = p;  t = q;
80
 
      aa = ioff[a] + a;
81
 
 
82
 
      value = 0.0;
83
 
      for (u=firstact; u<lastact; u++) {
84
 
        au = INDEX(a,u);
85
 
        value += opdm[t][u] * F_core[au];
86
 
      }
87
 
      grad[pair] = -2.0 * value;
88
 
 
89
 
      value = 0.0;
90
 
      for (u=firstact; u<lastact; u++) {
91
 
        tu = INDEX(t,u);
92
 
        au = INDEX(a,u); 
93
 
        /* loop over active v,w ... later restrict to symmetry allowed */
94
 
        for (v=firstact; v<lastact; v++) {
95
 
          for (w=firstact; w<lastact; w++) {
96
 
            vw = INDEX(v,w);
97
 
            tuvw = INDEX(tu,vw); 
98
 
            auvw = INDEX(au,vw);
99
 
            value += tpdm[tuvw] * tei[auvw];
100
 
          }
101
 
        }
102
 
      }
103
 
      grad[pair] -= 2.0 * value; 
104
 
       
105
 
    } 
106
 
 
107
 
    /* g_{ti}, i.e., active orb with inactive occ */
108
 
    else if (p >= firstact && q < firstact) {
109
 
      t = p;  i = q;
110
 
      tt = ioff[t] + t;
111
 
      ii = ioff[i] + i; 
112
 
      ti = ioff[t] + i;
113
 
      
114
 
      grad[pair] = -4.0 * (F_core[ti] + F_act[ti]);
115
 
 
116
 
      value = 0.0;
117
 
      for (u=firstact; u<lastact; u++) {
118
 
        tu = INDEX(t,u);
119
 
        iu = INDEX(i,u);
120
 
        value += opdm[t][u] * F_core[iu];
121
 
      }
122
 
      grad[pair] += 2.0 * value;
123
 
 
124
 
      value = 0.0;
125
 
      for (u=firstact; u<lastact; u++) {
126
 
        tu = INDEX(t,u);
127
 
        iu = INDEX(i,u);
128
 
 
129
 
        /* loop over active v,w ... later restrict to symmetry allowed */
130
 
        for (v=firstact; v<lastact; v++) {
131
 
          for (w=firstact; w<lastact; w++) {
132
 
            vw = INDEX(v,w);
133
 
            tuvw = INDEX(tu,vw); 
134
 
            iuvw = INDEX(iu,vw);
135
 
            value += tpdm[tuvw] * tei[iuvw];
136
 
          }
137
 
        }
138
 
         
139
 
      }
140
 
      grad[pair] += 2.0 * value; 
141
 
    }
142
 
 
143
 
    else {
144
 
      fprintf(outfile, 
145
 
             "(calc_grad_2): Error, unrecognized class of indep pair\n");
146
 
    }
147
 
 
148
 
  } /* end loop over pairs */
149
 
 
150
 
}
151
 
 
152