1
/*! \file xc_grad_fock.cc
3
\brief Enter brief description of file here
9
#include<libipv1/ip_lib.h>
10
#include<libciomr/libciomr.h>
11
#include<libpsio/psio.h>
12
#include<libint/libint.h>
23
#include"calc_den_fast.h"
25
#include"calc_grad_fast.h"
26
#include"functional.h"
29
#include"free_grid_structs.h"
31
#include"init_close_shell_info.h"
32
#include"calc_close_basis.h"
34
namespace psi { namespace CINTS {
36
void xc_grad_fock(void){
37
int i,j,k,l,m,n,q,s,t,u;
45
int moff,noff,mtmp,ntmp;
46
int am2shell1,am2shell2;
47
int shell_type1,shell_type2;
59
double *Gtri; /* Total and open-shell
61
triagonal form in SO basis */
64
double rpoints_double;
73
double four_pi_div_by_rps;
89
struct coordinates geom;
90
struct den_info_s den_info;
91
struct xc_info_s xc_info;
93
struct atomic_grid_s *atm_grd;
94
struct leb_chunk_s *chnk;
99
num_ao = BasisSet.num_ao;
101
DFT_options.basis = (double *)malloc(sizeof(double)*num_ao);
102
DFT_options.gradx = (double *)malloc(sizeof(double)*num_ao);
103
DFT_options.grady = (double *)malloc(sizeof(double)*num_ao);
104
DFT_options.gradz = (double *)malloc(sizeof(double)*num_ao);
105
DFT_options.gamma_basis = (double *)malloc(sizeof(double)*num_ao);
106
omega_arr = (double *)malloc(sizeof(double)*num_ao);
108
G = init_matrix(num_ao,num_ao);
110
/* ----Initialize Close shell data structure ---- */
112
DFT_options.close_shell_info = init_close_shell_info();
119
/*-------------------------------------------------------
120
Loop over symmetry-unique atoms only since integration
121
domains around symm.-unique atoms are equivalent
122
We are NOT employing the symmetry of angular grids
123
about atoms in special positions (on symm. elements)
125
-------------------------------------------------------*/
126
for(ua=0;ua<Symmetry.num_unique_atoms;ua++){
127
atm_grd = &(DFT_options.grid.atomic_grid[ua]);
129
atom = atm_grd->atom_num;
131
/*--- Cheap trick to get degeneracies of each unique atom ---*/
133
ua_deg = atm_grd->atom_degen;
134
ua_deg_d = (double) ua_deg;
136
xa = atm_grd->atom_center.x;
137
ya = atm_grd->atom_center.y;
138
za = atm_grd->atom_center.z;
140
bragg = atm_grd->Bragg_radii;
142
for(j=0;j<atm_grd->chunk_num;j++){
143
chnk = &(atm_grd->leb_chunk[j]);
144
timer_on("close basis");
146
calc_close_basis(ua,j);
148
close_shells = DFT_options.close_shell_info.num_close_shells;
149
close_aos = DFT_options.close_shell_info.num_close_aos;
150
timer_off("close basis");
151
for(k=0;k<chnk->size;k++){
153
sphr = &(chnk->spheres[k]);
159
drdq = sphr->drdq*bragg*bragg*bragg;
161
for(l=0;l<sphr->n_ang_points;l++){
162
pnt = &(sphr->points[l]);
163
/* ----------------------------------
164
Calculate the cartesian points of the point
165
relative to the center of mass
166
-------------------------------------*/
167
geom.x = bragg*pnt->p_cart.x+xa;
168
geom.y = bragg*pnt->p_cart.y+ya;
169
geom.z = bragg*pnt->p_cart.z+za;
171
/*-----------------------------------
172
Calculate the weighting funtion
173
----------------------------------*/
174
Becke_weight = weight_calc(atom,geom,3);
176
if(Becke_weight> WEIGHT_CUTOFF){
177
/*-----------------------------------
178
Get the density information for this
180
----------------------------------*/
182
/*den_info = DFT_options.den_calc(geom,ua);*/
183
den_info = calc_grad_fast(geom,ua);
187
if(den_info.den > DEN_CUTOFF){
188
/*-------------------------------------
190
-----------------------------------*/
192
ang_quad = 4.0*_pi*pnt->ang_weight;
194
/*------------------------------------
195
Calculate the potential functional
196
and energy functional at this
198
-----------------------------------*/
199
/*fprintf(outfile,"\nua_deg = %10.10lf",ua_deg_d);*/
200
den_val += 2.0*ua_deg_d*drdq*ang_quad
201
*Becke_weight*den_info.den;
203
xc_info.exch_info = DFT_options.
204
exchange_func(den_info);
206
xc_info.corr_info = DFT_options.
207
correlation_func(den_info);
209
exch_pval = ua_deg_d*drdq*ang_quad
210
*Becke_weight*xc_info.exch_info.dpval;
211
exch_gval = ua_deg_d*drdq*ang_quad
212
*Becke_weight*xc_info.exch_info.dgval;
214
corr_pval = ua_deg_d*drdq*ang_quad
215
*Becke_weight*xc_info.corr_info.dpval;
216
corr_gval = ua_deg_d*drdq*ang_quad
217
*Becke_weight*xc_info.corr_info.dgval;
219
pval = exch_pval+corr_pval;
220
gval = exch_gval+corr_gval;
222
exch_eval += ua_deg_d*drdq*ang_quad
223
*Becke_weight*xc_info.exch_info.eval;
224
corr_eval += ua_deg_d*drdq*ang_quad
225
*Becke_weight*xc_info.corr_info.eval;
226
eval += ua_deg_d*drdq*ang_quad
227
*Becke_weight*(xc_info.exch_info.eval+
228
xc_info.corr_info.eval);
230
/*------------------------------------
232
-----------------------------------*/
236
/* Form omega array < p+dphiu > */
238
dpx = den_info.gradx;
239
dpy = den_info.grady;
240
dpz = den_info.gradz;
242
/* Omega is the df/dgamma * delp . delp(basu*basv) */
243
/* See Benny Johnson's Thesis page 109 */
245
/*for(m=0;m<close_aos;m++){
246
omega_arr[m] = 2.0*gval*(dpx*DFT_options.gradx[m]
247
+dpy*DFT_options.grady[m]
248
+dpz*DFT_options.gradz[m]);
251
for(m=0;m<close_aos;m++){
252
bas1 = pval*DFT_options.basis[m]
254
moff = DFT_options.close_shell_info.
255
aos_close_to_chunk[m];
256
for(n=m;n<close_aos;n++){
257
bas2 = bas1*DFT_options.basis[n]
258
+DFT_options.basis[m]*omega_arr[n];
259
noff = DFT_options.close_shell_info.
260
aos_close_to_chunk[n];
262
G[moff][noff] += bas2;
264
G[noff][moff] += bas2;
267
/*fprintf(outfile,"\ngval = %10.10lf",gval);*/
268
for(m=0;m<close_aos;m++){
269
for(n=0;n<close_aos;n++){
270
moff = DFT_options.close_shell_info.
271
aos_close_to_chunk[m];
272
noff = DFT_options.close_shell_info.
273
aos_close_to_chunk[n];
274
bas1 = pval*DFT_options.basis[m]
275
*DFT_options.basis[n];
276
bas2 = 2.0*gval*((dpx*DFT_options.gradx[m]+
277
dpy*DFT_options.grady[m]+
278
dpz*DFT_options.gradz[m])*
280
+(dpx*DFT_options.gradx[n]+
281
dpy*DFT_options.grady[n]+
282
dpz*DFT_options.gradz[n])*
283
DFT_options.basis[m]);
284
/*fprintf(outfile,"\nbas1 = %10.10lf bas2 = %10.10lf\n\n\n",bas1,bas2);*/
285
G[moff][noff] += pval*DFT_options.basis[m]
286
*DFT_options.basis[n]
287
+2.0*gval*((dpx*DFT_options.gradx[m]+
288
dpy*DFT_options.grady[m]+
289
dpz*DFT_options.gradz[m])*
291
+(dpx*DFT_options.gradx[n]+
292
dpy*DFT_options.grady[n]+
293
dpz*DFT_options.gradz[n])*
294
DFT_options.basis[m]);
306
free_close_shell_info(DFT_options.close_shell_info);
307
print_mat(G,num_ao,num_ao,outfile);
308
for(m=0;m<num_ao;m++)
309
for(n=m;n<num_ao;n++)
316
/*----------------------
317
Unsort the Fock matrix
318
back to shell ordering
319
---------------------*/
322
/*----------------------
323
Transform to SO basis
324
----------------------*/
325
if (Symmetry.nirreps > 1 || BasisSet.puream) {
326
tmpmat1 = block_matrix(Symmetry.num_so,BasisSet.num_ao);
327
mmult(Symmetry.usotao,0,G,0,tmpmat1,0,Symmetry.num_so,BasisSet.num_ao,BasisSet.num_ao,0);
328
mmult(tmpmat1,0,Symmetry.usotao,1,G,0,Symmetry.num_so,BasisSet.num_ao,Symmetry.num_so,0);
329
if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
330
mmult(Symmetry.usotao,0,Go,0,tmpmat1,0,Symmetry.num_so,BasisSet.num_ao,BasisSet.num_ao,0);
331
mmult(tmpmat1,0,Symmetry.usotao,1,Go,0,Symmetry.num_so,BasisSet.num_ao,Symmetry.num_so,0);
336
/*-------------------------
337
Write G-matrices to disk
338
-------------------------*/
340
nstri = ioff[Symmetry.num_so];
341
Gtri = init_array(nstri);
342
sq_to_tri(G,Gtri,Symmetry.num_so);
344
fprintf(outfile,"\nDFT_energy = %10.10lf",eval);
345
fprintf(outfile,"\nX-Energy = %10.10lf",exch_eval);
346
fprintf(outfile,"\nC-Energy = %10.10lf",corr_eval);
347
fprintf(outfile,"\ntrace of density = %10.10lf\n",den_val);
349
psio_open(IOUnits.itapDSCF, PSIO_OPEN_OLD);
350
psio_write_entry(IOUnits.itapDSCF,"DFT X-energy",
351
(char *) &(exch_eval), sizeof(double));
352
psio_write_entry(IOUnits.itapDSCF,"DFT C-energy",
353
(char *) &(corr_eval), sizeof(double));
354
psio_write_entry(IOUnits.itapDSCF,"DFT XC-energy",
355
(char *) &(eval), sizeof(double));
356
psio_write_entry(IOUnits.itapDSCF,"DFT Den",
357
(char *) &(den_val), sizeof(double));
359
psio_write_entry(IOUnits.itapDSCF, "Total XC G-matrix"
360
, (char *) Gtri, sizeof(double)*nstri);
363
/*-- Cleanup the DFT stuff and close files --*/
364
cleanup_grid_type(DFT_options.grid);
365
psio_close(IOUnits.itapDSCF, 1);
added
removed
src/bin/cints/DFT/xc_grad_fock.cc
