← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/Tools/moinfo.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/cints/Tools/moinfo.c
1
 
#include<stdio.h>
2
 
#include<stdlib.h>
3
 
#include<math.h>
4
 
#include<libciomr/libciomr.h>
5
 
#include<libchkpt/chkpt.h>
6
 
#include<libint/libint.h>
7
 
 
8
 
#include"defines.h"
9
 
#define EXTERN
10
 
#include"global.h"
11
 
#include"small_fns.h"
12
 
 
13
 
/*-------------------------------
14
 
  Explicit function declarations
15
 
 -------------------------------*/
16
 
static void print_ccoefs(void);
17
 
 
18
 
void init_moinfo()
19
 
{
20
 
  int i, j;
21
 
  int irrep, iopen;
22
 
  int openirrs;
23
 
  double **ccvecs, *alpha, *beta;
24
 
  double **scf_evec_so;
25
 
 
26
 
  MOInfo.Escf = chkpt_rd_escf();
27
 
  MOInfo.Eref = 0.0;
28
 
  MOInfo.Ecorr = 0.0;
29
 
 
30
 
  /* CDS: I revised this stuff about correlation and SCF energies */
31
 
 
32
 
  /* If SCF, can say Eref = Escf (I guess...) */
33
 
  if (strcmp(UserOptions.wfn,"SCF")==0) {
34
 
    MOInfo.Eref = MOInfo.Escf;
35
 
  }
36
 
  else {
37
 
    MOInfo.Eref = chkpt_rd_eref();
38
 
  }
39
 
 
40
 
  /* Note: this init_moinfo() routine is not always called!  We'll need
41
 
     to re-grab some of the above energies in other subroutines to be
42
 
     positive we have them on-hand */
43
 
 
44
 
  MOInfo.num_mo = chkpt_rd_nmo();
45
 
  MOInfo.orbspi = chkpt_rd_orbspi();
46
 
  MOInfo.clsdpi = chkpt_rd_clsdpi();
47
 
  MOInfo.openpi = chkpt_rd_openpi();
48
 
  MOInfo.virtpi = init_int_array(Symmetry.nirreps);
49
 
  for(irrep=0;irrep<Symmetry.nirreps;irrep++)
50
 
    MOInfo.virtpi[irrep] = MOInfo.orbspi[irrep] - MOInfo.clsdpi[irrep] - MOInfo.openpi[irrep];
51
 
  MOInfo.scf_evals[0] = chkpt_rd_evals();
52
 
  scf_evec_so = chkpt_rd_scf();
53
 
  MOInfo.scf_evec[0] = block_matrix(MOInfo.num_mo,BasisSet.num_ao);
54
 
  mmult(Symmetry.usotao,1,scf_evec_so,0,MOInfo.scf_evec[0],1,BasisSet.num_ao,Symmetry.num_so,MOInfo.num_mo,0);
55
 
  free_block(scf_evec_so);
56
 
  if (UserOptions.reftype == uhf) {
57
 
    MOInfo.scf_evals[1] = chkpt_rd_beta_evals();
58
 
    scf_evec_so = chkpt_rd_beta_scf();
59
 
    MOInfo.scf_evec[1] = block_matrix(MOInfo.num_mo,BasisSet.num_ao);
60
 
    mmult(Symmetry.usotao,1,scf_evec_so,0,MOInfo.scf_evec[1],1,BasisSet.num_ao,Symmetry.num_so,MOInfo.num_mo,0);
61
 
    free_block(scf_evec_so);
62
 
  }
63
 
 
64
 
  /*--- Check the validity of the checkpoint file ---*/
65
 
  iopen = chkpt_rd_iopen();
66
 
  switch (UserOptions.reftype) {
67
 
  case rhf:     if (iopen != 0) punt("Content of checkpoint file inconsistent with REFERENCE\n"); break;
68
 
  case uhf:     if (iopen != 0) punt("Content of checkpoint file inconsistent with REFERENCE\n"); break;
69
 
  case rohf:    if (iopen <= 0) punt("Content of checkpoint file inconsistent with REFERENCE\n"); break;
70
 
  case twocon:  if (iopen >= 0) punt("Content of checkpoint file inconsistent with REFERENCE\n"); break;
71
 
  }
72
 
 
73
 
  /*--- Number of d.-o. MOs and s.-o. MOs ---*/
74
 
  MOInfo.ndocc = 0;
75
 
  openirrs = 0;
76
 
  MOInfo.nsocc = 0;
77
 
  for (i=0;i<Symmetry.nirreps;i++) {
78
 
    MOInfo.ndocc += MOInfo.clsdpi[i];
79
 
    if (MOInfo.openpi[i] != 0)
80
 
      openirrs++;
81
 
    MOInfo.nsocc += MOInfo.openpi[i];
82
 
  }
83
 
  MOInfo.nuocc = MOInfo.num_mo - MOInfo.ndocc - MOInfo.nsocc;
84
 
  
85
 
  /*--- Number of closed and open shells (no virtuals) ---*/
86
 
  MOInfo.num_moshells = MOInfo.num_openmoshells = ( -1 + (int)sqrt(1.0+8.0*abs(iopen)) )/2;    /* number of open shells */
87
 
  if (MOInfo.ndocc > 0)  /* add closed shells as well */
88
 
    MOInfo.num_moshells++;
89
 
 
90
 
  
91
 
  /*--- Read in open-shell coupling coeffcients ---*/
92
 
  ccvecs = chkpt_rd_ccvecs();
93
 
  if (iopen != 0) {    /*--- NOTE! These are Pitzer's coupling constants (a and b).
94
 
                         To get Yamaguchi's constants (alpha and beta) use this:
95
 
                         alpha = (1-a)/2  beta = (b-1)/4
96
 
                        ---*/
97
 
    alpha = ccvecs[0];
98
 
    beta = ccvecs[1];
99
 
  }
100
 
 
101
 
  if (UserOptions.reftype == twocon) {
102
 
    MOInfo.tcscf_occ[0] = 2.0/(1.0-alpha[0]);
103
 
    MOInfo.tcscf_occ[1] = 2.0/(1.0-alpha[2]);
104
 
  }
105
 
 
106
 
  if (UserOptions.reftype == rohf || UserOptions.reftype == twocon) {
107
 
    /*--- Form square matrices of coupling coeffcients ---*/
108
 
    MOInfo.Alpha = block_matrix(MOInfo.num_moshells,MOInfo.num_moshells);
109
 
    MOInfo.Beta  = block_matrix(MOInfo.num_moshells,MOInfo.num_moshells);
110
 
    /* Put alpha's and beta's in Alpha and Beta */
111
 
    if (MOInfo.ndocc > 0) { /* There are closed shells */
112
 
      /* Closed-Closed CCs */
113
 
      MOInfo.Alpha[0][0] = 2.0;
114
 
      MOInfo.Beta[0][0] = -1.0;
115
 
      if (UserOptions.reftype == rohf) {           /*--- Highspin and opeh-shell singlet cases ---*/
116
 
        /* Closed-Open CCs */
117
 
        for(i=1;i<MOInfo.num_moshells;i++) {
118
 
          MOInfo.Alpha[0][i] = MOInfo.Alpha[i][0] = 1.0;
119
 
          MOInfo.Beta[0][i] = MOInfo.Beta[i][0] = -0.5;
120
 
        }
121
 
        /* Open-Open blocks */
122
 
        for(i=0;i<MOInfo.num_openmoshells;i++)
123
 
          for(j=0;j<=i;j++) {
124
 
            MOInfo.Alpha[i+1][j+1] = MOInfo.Alpha[j+1][i+1] = (1.0 - alpha[ioff[i]+j]) * 0.5;
125
 
            MOInfo.Beta[i+1][j+1] = MOInfo.Beta[j+1][i+1] = (beta[ioff[i]+j] + 1.0) * -0.25;
126
 
          }
127
 
      }
128
 
      else if (UserOptions.reftype == twocon) {      /*--- TCSCF ---*/
129
 
        MOInfo.Alpha[0][1] = MOInfo.Alpha[1][0] = MOInfo.tcscf_occ[0];
130
 
        MOInfo.Beta[0][1] = MOInfo.Beta[1][0] = (-0.5) * MOInfo.tcscf_occ[0];
131
 
        MOInfo.Alpha[0][2] = MOInfo.Alpha[2][0] = MOInfo.tcscf_occ[1];
132
 
        MOInfo.Beta[0][2] = MOInfo.Beta[2][0] = (-0.5) * MOInfo.tcscf_occ[1];
133
 
        MOInfo.Alpha[1][1] = MOInfo.tcscf_occ[0] * 0.5;
134
 
        MOInfo.Alpha[2][2] = MOInfo.tcscf_occ[1] * 0.5;
135
 
        MOInfo.Beta[1][2] = MOInfo.Beta[2][1] = sqrt(MOInfo.tcscf_occ[0]*MOInfo.tcscf_occ[1]) * (-0.5);
136
 
      }
137
 
    }
138
 
    else { /* There are no closed shells */
139
 
      if (UserOptions.reftype == rohf) {           /*--- Highspin and opeh-shell singlet cases ---*/
140
 
        /* Open-Open blocks */
141
 
        for(i=0;i<MOInfo.num_openmoshells;i++)
142
 
          for(j=0;j<=i;j++) {
143
 
            MOInfo.Alpha[i][j] = MOInfo.Alpha[j][i] = (1.0 - alpha[ioff[i]+j]) * 0.5;
144
 
            MOInfo.Beta[i][j] = MOInfo.Beta[j][i] = (beta[ioff[i]+j] - 1.0) * 0.25;
145
 
          }
146
 
      }
147
 
      else if (UserOptions.reftype == twocon) {      /*--- TCSCF ---*/
148
 
        MOInfo.Alpha[0][0] = MOInfo.tcscf_occ[0] * 0.5;
149
 
        MOInfo.Alpha[1][1] = MOInfo.tcscf_occ[1] * 0.5;
150
 
        MOInfo.Beta[0][1] = MOInfo.Beta[1][0] = sqrt(MOInfo.tcscf_occ[0]*MOInfo.tcscf_occ[1]) * (-0.5);
151
 
      }
152
 
    }
153
 
    print_ccoefs();
154
 
  }
155
 
 
156
 
  free_block(ccvecs);
157
 
 
158
 
  return;
159
 
}
160
 
 
161
 
 
162
 
void cleanup_moinfo()
163
 
{
164
 
  free(MOInfo.scf_evals[0]);
165
 
  free_block(MOInfo.scf_evec[0]);
166
 
  if (UserOptions.reftype == uhf) {
167
 
    free(MOInfo.scf_evals[1]);
168
 
    free_block(MOInfo.scf_evec[1]);
169
 
  }
170
 
  free(MOInfo.orbspi);
171
 
  free(MOInfo.clsdpi);
172
 
  free(MOInfo.openpi);
173
 
  free(MOInfo.virtpi);
174
 
  if (UserOptions.reftype == rohf || UserOptions.reftype == twocon) {
175
 
    free_block(MOInfo.Alpha);
176
 
    free_block(MOInfo.Beta);
177
 
  }
178
 
 
179
 
  return;
180
 
}
181
 
 
182
 
 
183
 
/*----------------------
184
 
  Print coupling coeffs
185
 
 ----------------------*/
186
 
void print_ccoefs()
187
 
{
188
 
  int i,j;
189
 
 
190
 
  if (UserOptions.print_lvl >= PRINT_CCOEFF) {
191
 
    fprintf(outfile,"  -Yamaguchi's coupling coefficients :\n\n");
192
 
    fprintf(outfile,"    i  j    alpha     beta \n");
193
 
    fprintf(outfile,"   -------------------------\n");
194
 
    for (i=0;i<MOInfo.num_moshells;i++)
195
 
      for (j=0;j<=i;j++)
196
 
        fprintf(outfile,"   %2d %2d    %5.3f    %5.3f\n",
197
 
                i+1,j+1,MOInfo.Alpha[i][j],MOInfo.Beta[i][j]);
198
 
    fprintf(outfile,"\n\n");
199
 
  }
200
 
 
201
 
  return;
202
 
}