← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/bin/cints/Fock/hf_fock.cc

  • 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/Fock/hf_fock.cc
 
1
/*! \file
 
2
    \ingroup CINTS
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include<cmath>
 
6
#include<cstdio>
 
7
#include<cstring>
 
8
#include<cstdlib>
 
9
#include<memory.h>
 
10
#include<pthread.h>
 
11
#include<libipv1/ip_lib.h>
 
12
#include<libciomr/libciomr.h>
 
13
#include<libpsio/psio.h>
 
14
#include<libint/libint.h>
 
15
 
 
16
#include"defines.h"
 
17
#define EXTERN
 
18
#include"global.h"
 
19
#include <stdexcept>
 
20
 
 
21
#include"quartet_data.h"      /* From Default_Ints */
 
22
#include"norm_quartet.h"
 
23
#include"hash.h"
 
24
#include"transmat.h"
 
25
#ifdef USE_TAYLOR_FM
 
26
  #include"taylor_fm_eval.h"
 
27
#else
 
28
  #include"int_fjt.h"
 
29
#endif
 
30
#include"schwartz.h"
 
31
#include"shell_block_matrix.h"
 
32
 
 
33
using namespace std;
 
34
 
 
35
extern pthread_mutex_t fock_mutex;
 
36
 
 
37
namespace psi { 
 
38
  namespace CINTS {
 
39
extern void *hf_fock_thread(void *);
 
40
 
 
41
/*--- To be accessed by all HF Fock threads ---*/
 
42
double ****Gskel, ****Gskel_o;       /* Shell-blocked skeleton G matrices */
 
43
 
 
44
void hf_fock()
 
45
{
 
46
  pthread_attr_t thread_attr;
 
47
  pthread_t *thread_id;
 
48
  
 
49
  int nstri;
 
50
  int count ;
 
51
  int dum;
 
52
  int g, i, j, k, l, m, ii, jj, kk, ll;
 
53
  int si, sj, ni, nj, li, lj, si_g, sj_g;
 
54
 
 
55
  double temp;
 
56
  double **tmpmat1;
 
57
  double ****Gfull, ****Gfull_o;  /* Shell-blocked G matrices in AO basis*/
 
58
  double ****Gsym, ****Gsym_o;    /* Shell-blocked symmetrized (Gskel + Gskel(transp.)) G matrices */
 
59
  double ***ao_type_transmat;
 
60
  double *Gtri, *Gtri_o;          /* Total G matrices in lower*/
 
61
                                  /* triagonal form*/
 
62
  /*---------------
 
63
    Initialization
 
64
   ---------------*/
 
65
#ifdef USE_TAYLOR_FM
 
66
/*  init_Taylor_Fm_Eval(BasisSet.max_am*4-4,UserOptions.cutoff);*/
 
67
  init_Taylor_Fm_Eval(BasisSet.max_am*4-4,1.0E-20);
 
68
#else
 
69
  init_fjt(BasisSet.max_am*4);
 
70
#endif
 
71
  init_libint_base();
 
72
 
 
73
  /*------------------------------------------
 
74
    Compute integrals for Schwartz inequality
 
75
   ------------------------------------------*/
 
76
  schwartz_eri();
 
77
 
 
78
  /*------------------------------------
 
79
    Allocate shell-blocked skeleton G's
 
80
   ------------------------------------*/
 
81
  Gskel = init_shell_block_matrix();
 
82
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf)
 
83
    Gskel_o = init_shell_block_matrix();
 
84
 
 
85
  thread_id = (pthread_t *) malloc(UserOptions.num_threads*sizeof(pthread_t));
 
86
  pthread_attr_init(&thread_attr);
 
87
  pthread_attr_setscope(&thread_attr,
 
88
                        PTHREAD_SCOPE_SYSTEM);
 
89
  pthread_mutex_init(&fock_mutex,NULL);
 
90
  for(long int i=0;i<UserOptions.num_threads-1;i++)
 
91
    pthread_create(&(thread_id[i]),&thread_attr,
 
92
                   hf_fock_thread,(void *)i);
 
93
  hf_fock_thread( (void *) (UserOptions.num_threads - 1) );
 
94
  for(i=0;i<UserOptions.num_threads-1;i++)
 
95
    pthread_join(thread_id[i], NULL);
 
96
  free(thread_id);
 
97
  pthread_mutex_destroy(&fock_mutex);
 
98
  
 
99
  /*-------------------------------
 
100
    Gskel = Gskel + Gskel(transp.)
 
101
   -------------------------------*/
 
102
  Gsym = init_shell_block_matrix();
 
103
  GplusGt(Gskel,Gsym);
 
104
  free_shell_block_matrix(Gskel);
 
105
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
106
    Gsym_o = init_shell_block_matrix();
 
107
    GplusGt(Gskel_o,Gsym_o);
 
108
    free_shell_block_matrix(Gskel_o);
 
109
  }
 
110
  
 
111
  /*-----------------
 
112
    Symmetrize Gskel
 
113
   -----------------*/
 
114
  if (Symmetry.nirreps > 1) {
 
115
    ao_type_transmat = build_transmat(Symmetry.sym_oper, Symmetry.nirreps, BasisSet.max_am);
 
116
    Gfull = init_shell_block_matrix();
 
117
    for(g=0;g<Symmetry.nirreps;g++) {
 
118
      for(si=0;si<BasisSet.num_shells;si++) {
 
119
        ni = ioff[BasisSet.shells[si].am];
 
120
        li = BasisSet.shells[si].am-1;
 
121
        si_g = BasisSet.shells[si].trans_vec[g] - 1;
 
122
        for(sj=0;sj<BasisSet.num_shells;sj++) {
 
123
          sj_g = BasisSet.shells[sj].trans_vec[g] - 1;
 
124
          nj = ioff[BasisSet.shells[sj].am];
 
125
          lj = BasisSet.shells[sj].am-1;
 
126
 
 
127
          for(i=0;i<ni;i++)
 
128
            for(j=0;j<nj;j++)
 
129
              Gfull[si_g][sj_g][i][j] += ao_type_transmat[li][g][i]*
 
130
                                           ao_type_transmat[lj][g][j]*
 
131
                                           Gsym[si][sj][i][j];
 
132
        }
 
133
      }
 
134
    }
 
135
  }
 
136
  else
 
137
    Gfull = Gsym;
 
138
 
 
139
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
140
    if (Symmetry.nirreps > 1) {
 
141
    Gfull_o = init_shell_block_matrix();
 
142
    for(g=0;g<Symmetry.nirreps;g++) {
 
143
      for(si=0;si<BasisSet.num_shells;si++) {
 
144
        ni = ioff[BasisSet.shells[si].am];
 
145
        li = BasisSet.shells[si].am-1;
 
146
        si_g = BasisSet.shells[si].trans_vec[g] - 1;
 
147
        for(sj=0;sj<BasisSet.num_shells;sj++) {
 
148
          sj_g = BasisSet.shells[sj].trans_vec[g] - 1;
 
149
          nj = ioff[BasisSet.shells[sj].am];
 
150
          lj = BasisSet.shells[sj].am-1;
 
151
 
 
152
          for(i=0;i<ni;i++)
 
153
            for(j=0;j<nj;j++)
 
154
              Gfull_o[si_g][sj_g][i][j] += ao_type_transmat[li][g][i]*
 
155
                                           ao_type_transmat[lj][g][j]*
 
156
                                           Gsym_o[si][sj][i][j];
 
157
        }
 
158
      }
 
159
    }
 
160
    }
 
161
    else
 
162
      Gfull_o = Gsym_o;
 
163
  }
 
164
  /*--------------------
 
165
    Print out G for now
 
166
   --------------------*/
 
167
  G = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
168
  shell_block_to_block(Gfull,G);
 
169
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
170
    Go =  block_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
171
    shell_block_to_block(Gfull_o,Go);
 
172
  }
 
173
   /*----------------------
 
174
    Transform to SO basis
 
175
   ----------------------*/
 
176
  if (Symmetry.nirreps > 1 || BasisSet.puream) {
 
177
    tmpmat1 = block_matrix(Symmetry.num_so,BasisSet.num_ao);
 
178
    mmult(Symmetry.usotao,0,G,0,tmpmat1,0,Symmetry.num_so,BasisSet.num_ao,BasisSet.num_ao,0);
 
179
    mmult(tmpmat1,0,Symmetry.usotao,1,G,0,Symmetry.num_so,BasisSet.num_ao,Symmetry.num_so,0);
 
180
    if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
181
      mmult(Symmetry.usotao,0,Go,0,tmpmat1,0,Symmetry.num_so,BasisSet.num_ao,BasisSet.num_ao,0);
 
182
      mmult(tmpmat1,0,Symmetry.usotao,1,Go,0,Symmetry.num_so,BasisSet.num_ao,Symmetry.num_so,0);
 
183
    }
 
184
    free_block(tmpmat1);
 
185
  }
 
186
  /*
 
187
  fprintf(outfile,"  Closed-shell Fock matrix in SO basis:\n");
 
188
  print_mat(G,Symmetry.num_so,Symmetry.num_so,outfile);
 
189
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
190
      fprintf(outfile,"  Open-shell Fock matrix in SO basis:\n");
 
191
      print_mat(Go,Symmetry.num_so,Symmetry.num_so,outfile);
 
192
      } */
 
193
  
 
194
  /*-------------------------
 
195
    Write G-matrices to disk
 
196
   -------------------------*/
 
197
  nstri = ioff[Symmetry.num_so];
 
198
  Gtri = init_array(nstri);
 
199
  sq_to_tri(G,Gtri,Symmetry.num_so);
 
200
  free_block(G);
 
201
  psio_open(IOUnits.itapDSCF, PSIO_OPEN_OLD);
 
202
  switch (UserOptions.reftype) {
 
203
  case rohf:
 
204
      Gtri_o = init_array(nstri);
 
205
      sq_to_tri(Go,Gtri_o,Symmetry.num_so);
 
206
      free_block(Go);
 
207
      psio_write_entry(IOUnits.itapDSCF, "Open-shell JX G-matrix", (char *) Gtri_o, sizeof(double)*nstri);
 
208
      free(Gtri_o);
 
209
  case rhf:
 
210
      psio_write_entry(IOUnits.itapDSCF, "Total JX G-matrix", (char *) Gtri, sizeof(double)*nstri);
 
211
      free(Gtri);
 
212
      break;
 
213
 
 
214
  case uhf:
 
215
      Gtri_o = init_array(nstri);
 
216
      sq_to_tri(Go,Gtri_o,Symmetry.num_so);
 
217
      free_block(Go);
 
218
      /*--- Form alpha and beta Fock matrices first and then write them out ---*/
 
219
      for(i=0;i<nstri;i++) {
 
220
        temp = Gtri[i] + Gtri_o[i];
 
221
        Gtri[i] = Gtri[i] - Gtri_o[i];
 
222
        Gtri_o[i] = temp;
 
223
      }
 
224
      
 
225
      psio_write_entry(IOUnits.itapDSCF, "Alpha JX G-matrix", (char *) Gtri, sizeof(double)*nstri);
 
226
      psio_write_entry(IOUnits.itapDSCF, "Beta JX G-matrix", (char *) Gtri_o, sizeof(double)*nstri);
 
227
      free(Gtri);
 
228
      free(Gtri_o);
 
229
      break;
 
230
  }
 
231
  psio_close(IOUnits.itapDSCF, 1);
 
232
 
 
233
  /*---------
 
234
    Clean-up
 
235
   ---------*/
 
236
  free_shell_block_matrix(Gsym);
 
237
  if (Symmetry.nirreps > 1) {
 
238
    free_shell_block_matrix(Gfull);
 
239
  }
 
240
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
241
    free_shell_block_matrix(Gsym_o);
 
242
    if (Symmetry.nirreps > 1)
 
243
      free_shell_block_matrix(Gfull_o);
 
244
  }
 
245
#ifdef USE_TAYLOR_FM
 
246
  free_Taylor_Fm_Eval();
 
247
#else
 
248
  free_fjt();
 
249
#endif
 
250
 
 
251
  return;
 
252
}
 
253
 
 
254
};};