← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/cints/DFT/xc_fock.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cints/DFT/xc_fock.c
 
1
#include<math.h>
 
2
#include<stdio.h>
 
3
#include<string.h>
 
4
#include<stdlib.h>
 
5
#include<libipv1/ip_lib.h>
 
6
#include<libciomr/libciomr.h>
 
7
#include<libpsio/psio.h>
 
8
#include<libint/libint.h>
 
9
#include<pthread.h>
 
10
#include<libqt/qt.h>
 
11
 
 
12
#include"defines.h"
 
13
#define EXTERN
 
14
#include"global.h"
 
15
 
 
16
#include"dft_init.h"
 
17
#include"weighting.h"
 
18
#include"calc_den_fast.h"
 
19
#include"calc_den.h"
 
20
#include"functional.h"
 
21
#include"physconst.h"
 
22
#include"grid_init.h"
 
23
#include"free_grid_structs.h"
 
24
#include"dcr.h"
 
25
#include"init_close_shell_info.h"
 
26
#include"calc_close_basis.h"
 
27
 
 
28
void xc_fock(void){
 
29
  int i,j,k,l,m,n,q,s,t,u;
 
30
  int ua, atom, ua_deg;
 
31
  int rpoints;
 
32
  int ang_points;
 
33
  int num_ao;
 
34
  int point_count=0;
 
35
  int dum;
 
36
  int num;
 
37
  int moff,noff,mtmp,ntmp;
 
38
  int am2shell1,am2shell2;
 
39
  int shell_type1,shell_type2;
 
40
  int ioff1,ioff2;
 
41
  int chek = 1;
 
42
  int close_shells;
 
43
  int close_aos;
 
44
  int tmp;
 
45
  
 
46
  int nstri;
 
47
  double temp;
 
48
  double *temp_arr;
 
49
  double **tmpmat1;
 
50
  double *Gtri;  /* Total and open-shell 
 
51
                             G matrices and lower 
 
52
                             triagonal form in SO basis */  
 
53
  double r;
 
54
  double rind;
 
55
  double rpoints_double;
 
56
  double ua_deg_d;
 
57
  double bragg;
 
58
  double qr;
 
59
  double drdq;
 
60
  double jacobian;
 
61
  double xa,ya,za;
 
62
  double Becke_weight;
 
63
  double ang_quad;
 
64
  double four_pi_div_by_rps;
 
65
  double den_val=0.0;
 
66
  double exch_vval=0.0;
 
67
  double corr_vval=0.0;
 
68
  double exch_eval=0.0;
 
69
  double corr_eval=0.0;
 
70
  double vval = 0.0;
 
71
  double eval = 0.0;
 
72
  double bas1 = 0.0;
 
73
  double bas2 = 0.0;
 
74
  double vvalbas = 0.0;
 
75
  
 
76
  struct coordinates geom;
 
77
  struct den_info_s den_info;
 
78
  struct xc_info_s xc_info;
 
79
  
 
80
  struct atomic_grid_s *atm_grd;
 
81
  struct leb_chunk_s *chnk;
 
82
  leb_sphere_t *sphr;
 
83
  leb_point_t *pnt;
 
84
  
 
85
  
 
86
  num_ao = BasisSet.num_ao;
 
87
  DFT_options.basis = (double *)malloc(sizeof(double)*num_ao);
 
88
  
 
89
  if(UserOptions.reftype == rhf)
 
90
      G = init_matrix(num_ao,num_ao);
 
91
  
 
92
  /* ----Initialize Close shell data structure ---- */
 
93
  
 
94
  DFT_options.close_shell_info = init_close_shell_info();
 
95
 
 
96
  timer_init();
 
97
  timer_on("DFT");
 
98
  
 
99
  grid_init();
 
100
  
 
101
/*-------------------------------------------------------
 
102
    Loop over symmetry-unique atoms only since integration
 
103
    domains around symm.-unique atoms are equivalent
 
104
    We are NOT employing the symmetry of angular grids
 
105
    about atoms in special positions (on symm. elements)
 
106
    like Handy & co. do
 
107
   -------------------------------------------------------*/
 
108
  for(ua=0;ua<Symmetry.num_unique_atoms;ua++){
 
109
      atm_grd = &(DFT_options.grid.atomic_grid[ua]);
 
110
      
 
111
      atom = atm_grd->atom_num;
 
112
      
 
113
/*--- Cheap trick to get degeneracies of each unique atom ---*/
 
114
      
 
115
      ua_deg = atm_grd->atom_degen;
 
116
      ua_deg_d = (double) ua_deg;
 
117
 
 
118
      xa = atm_grd->atom_center.x;
 
119
      ya = atm_grd->atom_center.y;
 
120
      za = atm_grd->atom_center.z;
 
121
      
 
122
      bragg = atm_grd->Bragg_radii;
 
123
      
 
124
      for(j=0;j<atm_grd->chunk_num;j++){
 
125
          chnk = &(atm_grd->leb_chunk[j]);
 
126
          timer_on("close basis");
 
127
          
 
128
          calc_close_basis(ua,j);
 
129
          
 
130
          close_shells = DFT_options.close_shell_info.num_close_shells;
 
131
          close_aos = DFT_options.close_shell_info.num_close_aos;
 
132
          timer_off("close basis");
 
133
          for(k=0;k<chnk->size;k++){
 
134
              
 
135
              sphr = &(chnk->spheres[k]);
 
136
              
 
137
              r = sphr->r*bragg;
 
138
              /*r = 0.87895;*/
 
139
              
 
140
              
 
141
              drdq = sphr->drdq*bragg*bragg*bragg;
 
142
              fprintf(outfile,"\nr = %e drdq = %e ang = %d",r,drdq,sphr->n_ang_points);
 
143
              for(l=0;l<sphr->n_ang_points;l++){
 
144
                  pnt = &(sphr->points[l]);
 
145
                  /* ----------------------------------
 
146
                     Calculate the cartesian points of the point
 
147
                     relative to the center of mass
 
148
                     -------------------------------------*/
 
149
                  geom.x = bragg*pnt->p_cart.x+xa;
 
150
                  geom.y = bragg*pnt->p_cart.y+ya;
 
151
                  geom.z = bragg*pnt->p_cart.z+za;
 
152
                  
 
153
                  /*-----------------------------------
 
154
                    Calculate the weighting funtion 
 
155
                    ----------------------------------*/
 
156
                  Becke_weight = weight_calc(atom,geom,3);
 
157
                  
 
158
                  if(Becke_weight> WEIGHT_CUTOFF){
 
159
                      /*-----------------------------------
 
160
                        Get the density information for this 
 
161
                        point
 
162
                        ----------------------------------*/
 
163
                      timer_on("DEN1");
 
164
                      den_info = DFT_options.den_calc(geom,ua);
 
165
                      
 
166
                      timer_off("DEN1");
 
167
                      
 
168
                      if(den_info.den > DEN_CUTOFF){
 
169
                          /*-------------------------------------
 
170
                            Weight from Lebedev
 
171
                            -----------------------------------*/
 
172
                      
 
173
                          ang_quad = 4.0*_pi*pnt->ang_weight;
 
174
                      
 
175
                          /*------------------------------------
 
176
                            Calculate the potential functional
 
177
                            and energy functional at this
 
178
                            point
 
179
                            -----------------------------------*/
 
180
                          /*fprintf(outfile,"\nua_deg = %10.10lf",ua_deg_d);*/
 
181
                          den_val += 2.0*ua_deg_d*drdq*ang_quad
 
182
                              *Becke_weight*den_info.den;
 
183
                          
 
184
                          xc_info.exch_info = DFT_options.
 
185
                              exchange_func(den_info);
 
186
                          
 
187
                          xc_info.corr_info = DFT_options.
 
188
                              correlation_func(den_info);
 
189
                          
 
190
                          exch_vval = ua_deg_d*drdq*ang_quad
 
191
                              *Becke_weight*xc_info.exch_info.dval;
 
192
                          corr_vval = ua_deg_d*drdq*ang_quad
 
193
                              *Becke_weight*xc_info.corr_info.dval;
 
194
                          
 
195
                          vval = exch_vval+corr_vval;
 
196
                          
 
197
                          exch_eval += ua_deg_d*drdq*ang_quad
 
198
                              *Becke_weight*xc_info.exch_info.eval;
 
199
                          corr_eval += ua_deg_d*drdq*ang_quad
 
200
                              *Becke_weight*xc_info.corr_info.eval;
 
201
                          eval += ua_deg_d*drdq*ang_quad
 
202
                              *Becke_weight*(xc_info.exch_info.eval+
 
203
                                             xc_info.corr_info.eval);
 
204
                          
 
205
                          /*------------------------------------
 
206
                            Update the G matrix
 
207
                            -----------------------------------*/
 
208
                          timer_on("FOCK");
 
209
                          t=0;
 
210
                          
 
211
                          for(m=0;m<close_aos;m++){
 
212
                              bas1 = vval*DFT_options.basis[m];
 
213
                              moff = DFT_options.close_shell_info.
 
214
                                  aos_close_to_chunk[m];
 
215
                              for(n=m;n<close_aos;n++){
 
216
                                  bas2 = bas1*DFT_options.basis[n];
 
217
                                  noff = DFT_options.close_shell_info.
 
218
                                      aos_close_to_chunk[n];
 
219
                                  if(noff > moff)
 
220
                                      G[moff][noff] += bas2;
 
221
                                  else
 
222
                                      G[noff][moff] += bas2;
 
223
                              }
 
224
                          }
 
225
                          timer_off("FOCK");
 
226
                      }
 
227
                      
 
228
                  }
 
229
              }
 
230
          }
 
231
      }
 
232
      
 
233
  }
 
234
 
 
235
  
 
236
  free_close_shell_info(DFT_options.close_shell_info);
 
237
  /*print_mat(G,num_ao,num_ao,outfile);*/
 
238
   for(m=0;m<num_ao;m++)
 
239
      for(n=m;n<num_ao;n++)
 
240
          G[n][m]=G[m][n];
 
241
  
 
242
  timer_off("DFT");
 
243
  timer_done();
 
244
  
 
245
  
 
246
  /*----------------------
 
247
    Unsort the Fock matrix
 
248
    back to shell ordering
 
249
    ---------------------*/
 
250
  
 
251
  
 
252
/*----------------------
 
253
  Transform to SO basis
 
254
  ----------------------*/
 
255
  if (Symmetry.nirreps > 1 || BasisSet.puream) {
 
256
      tmpmat1 = block_matrix(Symmetry.num_so,BasisSet.num_ao);
 
257
      mmult(Symmetry.usotao,0,G,0,tmpmat1,0,Symmetry.num_so,BasisSet.num_ao,BasisSet.num_ao,0);
 
258
      mmult(tmpmat1,0,Symmetry.usotao,1,G,0,Symmetry.num_so,BasisSet.num_ao,Symmetry.num_so,0);
 
259
      if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
260
          mmult(Symmetry.usotao,0,Go,0,tmpmat1,0,Symmetry.num_so,BasisSet.num_ao,BasisSet.num_ao,0);
 
261
          mmult(tmpmat1,0,Symmetry.usotao,1,Go,0,Symmetry.num_so,BasisSet.num_ao,Symmetry.num_so,0);
 
262
      }
 
263
      free_block(tmpmat1);
 
264
  }
 
265
  
 
266
  /*-------------------------
 
267
    Write G-matrices to disk
 
268
    -------------------------*/
 
269
  
 
270
  nstri = ioff[Symmetry.num_so];
 
271
  Gtri = init_array(nstri);
 
272
  sq_to_tri(G,Gtri,Symmetry.num_so);
 
273
  free_block(G);
 
274
  fprintf(outfile,"\nDFT_energy = %10.10lf",eval);
 
275
  fprintf(outfile,"\nX-Energy = %10.10lf",exch_eval);
 
276
  fprintf(outfile,"\nC-Energy = %10.10lf",corr_eval);
 
277
  fprintf(outfile,"\ntrace of density = %10.10lf\n",den_val);
 
278
  
 
279
  psio_open(IOUnits.itapDSCF, PSIO_OPEN_OLD);
 
280
  psio_write_entry(IOUnits.itapDSCF,"DFT X-energy",
 
281
                   (char *) &(exch_eval), sizeof(double));
 
282
  psio_write_entry(IOUnits.itapDSCF,"DFT C-energy",
 
283
                   (char *) &(corr_eval), sizeof(double));
 
284
  psio_write_entry(IOUnits.itapDSCF,"DFT XC-energy",
 
285
                   (char *) &(eval), sizeof(double));
 
286
  psio_write_entry(IOUnits.itapDSCF,"DFT Den",
 
287
                   (char *) &(den_val), sizeof(double));
 
288
  
 
289
  psio_write_entry(IOUnits.itapDSCF, "Total XC G-matrix"
 
290
                   , (char *) Gtri, sizeof(double)*nstri);
 
291
  free(Gtri);
 
292
 
 
293
  /*-- Cleanup the DFT stuff and close files --*/
 
294
  cleanup_grid_type(DFT_options.grid);
 
295
  psio_close(IOUnits.itapDSCF, 1);
 
296
  return;
 
297
}
 
298
 
 
299
 
 
300
 
 
301