← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/tools/ga-5-2/gfutex/examples/scf/scf.c

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools/ga-5-2/gfutex/examples/scf/scf.c
 
1
/* C version of the GA SCF code */
 
2
 
 
3
#include <stdio.h>
 
4
#include <stdlib.h>
 
5
#include <math.h>
 
6
#include <assert.h>
 
7
 
 
8
#include <ga.h>
 
9
#include <macdecls.h>
 
10
 
 
11
#include "cscc.h"
 
12
#include "integ.h"
 
13
#include "input.h"
 
14
#include "output.h"
 
15
#include "oneel.h"
 
16
#include "twoel.h"
 
17
#include "scf.h"
 
18
#include "timer.h"
 
19
 
 
20
// #include "global.h"
 
21
// #include "tcgmsg.h"
 
22
 
 
23
#include "sum.h"
 
24
 
 
25
//function declaration
 
26
void makesz(double* schwmax);
 
27
void ininrm(void);
 
28
double s(int i,int j);
 
29
void makden(void);
 
30
int nxtask(int nproc);
 
31
int verify_zero_chunk(double chunk[][ichunk]);
 
32
void damp(double* fac);
 
33
void prnout(int iter, double energy, double deltad, double tester);
 
34
double dendif(void);
 
35
double testfock(void);
 
36
void shiftfock(double shift);
 
37
void prnfin(double energy);
 
38
void diagon(double* tester, int iter);
 
39
void makeob(void);
 
40
void denges(void);
 
41
void setarrays(void);
 
42
void closearrays(void);
 
43
void makoverlap(void);
 
44
void setoverlap(double *a,int* lo,int* hi, int ld1,int ld2);
 
45
void print_GA_block(int g_a);
 
46
void print_GA_block_ij(int g_a,int tlo);
 
47
void dump_chunk(double *a,int ld1,int ld2);
 
48
// double abs_double(double x); 
 
49
 
 
50
#ifdef CLOG
 
51
FILE *clog = NULL;
 
52
#endif
 
53
 
 
54
double enrep, q[maxatom], ax[maxatom], ay[maxatom], az[maxatom],
 
55
  x[maxnbfn], y[maxnbfn], z[maxnbfn], expnt[maxnbfn], rnorm[maxnbfn];
 
56
long long int iky[maxnbfn], nocc, nbfn, nnbfn;
 
57
long long int icut1,icut2,icut3,icut4; 
 
58
int natom; //long long int --> long 
 
59
 
 
60
double eigv[maxnbfn];
 
61
int g_counter, g_dens, g_fock, g_tfock, g_schwarz, g_work, g_ident, g_orbs;
 
62
int g_tmp, g_proc; //temporay global array for storage major transformation
 
63
 
 
64
int main(int argc, char **argv) 
 
65
{
 
66
  long long int nints, maxint; 
 
67
 
 
68
  //     CAUTION: int precision requirements;
 
69
  //     nints, maxint, etc. are proportional to the number of basis functions;
 
70
  //     to the fourth power! 216^4 is greater than the largest number;
 
71
  //     that can be represented as a 32-bit signed interger, so 64-bit;
 
72
  //     arithmetic is needed to count integrals when calculating more than;
 
73
  //     14 Be atoms with 15 basis functions each. Since integrals are counted;
 
74
  //     over all iterations, 18 iterations with 7 atoms can result in precision;
 
75
  //     problems. Note that the wave function could be calculated correctly;
 
76
  //     for much larger basis sets without 64-bit ints because the required;
 
77
  //     indexing is usually proportional to nbfn^2, which is good to 46,340;
 
78
  //     basis functions, except that the task counter runs as (nbfn/ichunk)^4,;
 
79
  //     so with ichunk = 10, 32-bit ints yield correct wavefunctions out to;
 
80
  //     2145 basis functions (maxatom=143), or 4290 (maxatom=286) with ichunk =;
 
81
  //     20, ...;
 
82
  //;
 
83
  //     This warning applies to the Global Arrays implementation as well!;
 
84
  //     functions of special concern are GA_igop and NGA_Read_inc.;
 
85
  //;
 
86
#define USE_TRANSFORM 1
 
87
#define MPI 1
 
88
 
 
89
  int heap, stack;
 
90
  double tinit=0.0, tonel=0.0, ttwoel=0.0, tdiag=0.0, tdens=0.0, tprint=0.0;
 
91
  double eone=0.0, etwo=0.0, energy=0.0, deltad=0.0;
 
92
      
 
93
  //implicit variables in fortran
 
94
  int me;
 
95
  double frac;
 
96
  int iter;
 
97
  double scale;
 
98
  int nproc = 0;
 
99
  double rjunk;
 
100
  double totsec, tester;
 
101
  double schwmax;
 
102
#ifdef CLOG
 
103
  char fname[50];
 
104
#endif
 
105
 
 
106
  //     initalize the parallel message passing environment;
 
107
#ifdef MPI
 
108
  MPI_Init(&argc,&argv);
 
109
#else
 
110
  //pbeginf();
 
111
  tcg_pbeginf();
 
112
#endif
 
113
      
 
114
  GA_Initialize();
 
115
      
 
116
  //   Allocate memory;
 
117
  heap =  32000000;
 
118
  stack = 64000000;
 
119
  if (!MA_init(MT_DBL, stack, heap))
 
120
    GA_Error("ma_init failed",-1);
 
121
 
 
122
#if 1
 
123
  GFInitialize();
 
124
#endif
 
125
 
 
126
  me = GA_Nodeid();
 
127
  nproc = GA_Nnodes();
 
128
      
 
129
  //     initialize a bunch of stuff and initial density matrix;
 
130
  rjunk = timer();
 
131
 
 
132
#ifdef CLOG
 
133
  sprintf(fname, "clog.dat.%d", me);
 
134
  clog = fopen(fname, "w");
 
135
#endif
 
136
 
 
137
  //     get input from file be.inpt;
 
138
  input();
 
139
 
 
140
  //     create and allocate global arrays;
 
141
  setarrays();
 
142
 
 
143
  if (me == 0) {
 
144
    printf(" bytes of memory used by node 0: %lld\n\n", GA_Inquire_memory());
 
145
  }
 
146
 
 
147
  ininrm();
 
148
 
 
149
  //     create initial guess for density matrix by using single atom;
 
150
  //     densities;
 
151
  denges();
 
152
 
 
153
#if 1
 
154
  GA_Print_distribution(g_schwarz);
 
155
  GA_Print_distribution(g_dens);
 
156
#endif
 
157
 
 
158
#if USE_TRANSFORM
 
159
  //     make initial orthogonal orbital set for solution method using;
 
160
  //     similarity transform;
 
161
  makeob();
 
162
#endif
 
163
  //     make info for sparsity test;
 
164
  makesz(&schwmax);
 
165
 
 
166
  tinit = timer();
 
167
 
 
168
  //     print preliminary data before any long compute segments start;
 
169
  if (me == 0)
 
170
    printf("\n");
 
171
 
 
172
      //     ^* iterate ^*;
 
173
  for (iter = 0; iter < mxiter; iter++) {
 
174
    double s;
 
175
 
 
176
    //     make the one particle contribution to the fock matrix;
 
177
    //     and get the partial contribution to the energy;
 
178
    oneel(schwmax, &eone, nproc);
 
179
    tonel = tonel + timer();
 
180
 
 
181
#if 0
 
182
    s = sum(g_fock);
 
183
    GA_Dgop(&s, 1, "+");
 
184
 
 
185
    if (me == 0)
 
186
      printf("g_fock sum after oneel: %f\n", s);
 
187
 
 
188
    s = sum(g_orbs);
 
189
    GA_Dgop(&s, 1, "+");
 
190
 
 
191
    if (me == 0)
 
192
      printf("g_orbs sum after oneel: %f\n", s);
 
193
#endif
 
194
        
 
195
    //     compute the two particle contributions to the fock matrix and;
 
196
    //     get the total energy.;
 
197
    twoel(schwmax, &etwo, nproc);
 
198
    ttwoel = ttwoel + timer();
 
199
 
 
200
#if 0
 
201
    s = sum(g_fock);
 
202
    GA_Dgop(&s, 1, "+");
 
203
 
 
204
    if (me == 0)
 
205
      printf("g_fock sum after twoel: %f\n", s);
 
206
#endif
 
207
 
 
208
    //     Diagonalize the fock matrix. The diagonalizers used in this;
 
209
    //     are actually sequential, not parallel.;
 
210
    diagon(&tester, iter);
 
211
    tdiag = tdiag + timer();
 
212
 
 
213
#if 0
 
214
    s = sum(g_fock);
 
215
    GA_Dgop(&s, 1, "+");
 
216
 
 
217
    if (me == 0)
 
218
      printf("g_fock sum after diagon: %f\n", s);
 
219
 
 
220
    s = sum(g_orbs);
 
221
    GA_Dgop(&s, 1, "+");
 
222
 
 
223
    if (me == 0)
 
224
      printf("g_orbs sum after diagon: %f\n", s);
 
225
#endif
 
226
 
 
227
    //     make the new density matrix in g_work from orbitals in g_orbs,;
 
228
    //     compute the norm of the change in the density matrix and;
 
229
    //     then update the density matrix in g_dens with damping.;
 
230
    makden();
 
231
         
 
232
    deltad = dendif();
 
233
 
 
234
    if (iter == 0)
 
235
      scale = 0.00;
 
236
    else if (iter < 5) {
 
237
      if (nbfn > 60)
 
238
        scale = 0.50;
 
239
      else
 
240
        scale = 0.00;
 
241
    } 
 
242
    else
 
243
      scale = 0.00;
 
244
 
 
245
    damp(&scale);
 
246
    tdens = tdens + timer();
 
247
          
 
248
    //     add up energy and print out convergence information;
 
249
    if (me == 0) {
 
250
      energy = enrep + eone + etwo;
 
251
      prnout(iter, energy, deltad, tester);
 
252
      tprint = tprint + timer();
 
253
    }
 
254
 
 
255
    //     if converged exit iteration loop;
 
256
    if (deltad < tol)
 
257
      goto L20;
 
258
 
 
259
#if GA_VERSION_MAJOR >= 5
 
260
    GA_Llgop(&icut4, 1, "+");
 
261
#else
 
262
    GA_Igop(&icut4, 1, "+");
 
263
#endif
 
264
 
 
265
    if(icut4 == 0) {
 
266
      //     something has gone wrong--print what you know and quit.;
 
267
      printf("no two-electron integrals computed!\n");
 
268
      goto L20;
 
269
    }
 
270
 
 
271
  }
 
272
 
 
273
  iter = iter - 1;
 
274
  if(me == 0)
 
275
    printf("SCF failed to converge in %d iters\n", iter);
 
276
  //...v....1....v....2....v....3....v....4....v....5....v....6....v....7..;
 
277
  //;
 
278
  //     finished ... print out eigenvalues and occupied orbitals;
 
279
  //;
 
280
 
 
281
 L20:
 
282
 
 
283
#if GA_VERSION_MAJOR >= 5
 
284
  GA_Llgop(&icut1, 1, (char*) "+");
 
285
  GA_Llgop(&icut2, 1, (char*) "+");
 
286
  GA_Llgop(&icut3, 1, (char*) "+");
 
287
#else
 
288
  GA_Igop(&icut1, 1, (char*) "+");
 
289
  GA_Igop(&icut2, 1, (char*) "+");
 
290
  GA_Igop(&icut3, 1, (char*) "+");
 
291
#endif
 
292
 
 
293
  if (me == 0) {
 
294
    //     print out timing information;
 
295
    prnfin(energy);
 
296
    printf("      init       onel      twoel       diag       dens       print       ncpu    \n");
 
297
    printf("      ----       ----      -----       ----       ----       ----        ----    \n");
 
298
    totsec = tinit + tonel + ttwoel + tdiag + tdens + tprint;
 
299
    printf("%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f        %d",tinit, tonel, ttwoel, tdiag,
 
300
           tdens,tprint,nproc);
 
301
    printf("\n elapsed time in seconds %10.2f\n\n",totsec);
 
302
    //     print out information on # integrals evaluated each iteration;
 
303
    nints = icut1 + icut2 + icut3;
 
304
    frac = (double)icut3 / (double)nints;
 
305
    printf("No. of integrals screened or computed (all iters) \n\n");
 
306
    printf("-------------------------------------\n");
 
307
    printf("     failed #ij    test failed #kl    test #compute    #total    fraction\n");
 
308
    printf("     ----------    ----------------   -------------    ------    --------\n");
 
309
    printf("%15lld %15lld %15lld %15lld %9.6f\n", icut1, icut2, icut3, nints, frac);
 
310
         
 
311
    maxint = nbfn;
 
312
    maxint = pow(maxint, 4) * (iter + 1); 
 
313
    if(nints != maxint) {
 
314
      printf("Inconsistent number of integrals, should be %ld\n", maxint);
 
315
      printf("Note: largest 32-bit int is 2,147,483,647\n");
 
316
      printf("Probably due to insufficient int precision in GA.\n");
 
317
    }
 
318
#ifndef MPI
 
319
    tcg_stats();
 
320
#endif
 
321
  }
 
322
 
 
323
  closearrays();
 
324
 
 
325
#ifdef CLOG
 
326
  fprintf(stderr, "Before GFFinalize()\n");
 
327
  fflush(stderr);
 
328
#endif
 
329
 
 
330
#ifdef CLOG
 
331
  fclose(clog);
 
332
#endif
 
333
 
 
334
#if 1
 
335
  GFFinalize();
 
336
#endif
 
337
 
 
338
#ifdef CLOG
 
339
  fprintf(stderr, "After GFFinalize()\n");
 
340
  fflush(stderr);
 
341
#endif
 
342
      
 
343
  GA_Terminate();
 
344
 
 
345
#ifdef CLOG
 
346
  fprintf(stderr, "After GA_Terminate()\n");
 
347
  fflush(stderr);
 
348
#endif
 
349
 
 
350
#ifdef MPI
 
351
  MPI_Finalize();
 
352
#else
 
353
  tcg_pend();
 
354
#endif
 
355
 
 
356
}
 
357
 
 
358
void makesz(double *schwmax)
 
359
{
 
360
  double work[ichunk][ichunk];
 
361
  int lo[2],hi[2],i,j,iloc,jloc,ld;
 
362
  int dotask;
 
363
 
 
364
  int lo_c[2],hi_c[2]; //column-major inter-patch switch
 
365
 
 
366
  //implicit declared variable in fortran
 
367
  double gg;
 
368
    
 
369
  //schwarz(ij) = (ij|ij) for sparsity test;
 
370
  icut1 = 0;
 
371
  icut2 = 0;
 
372
  icut3 = 0;
 
373
 
 
374
  GA_Zero(g_schwarz); 
 
375
  GA_Zero(g_counter);
 
376
  *schwmax = 0.00;
 
377
  dotask = next_chunk(lo, hi);
 
378
  ld = ichunk;
 
379
 
 
380
  while (dotask) {
 
381
    for (i = lo[0]; i <= hi[0]; i++) {
 
382
      iloc = i - lo[0];
 
383
      for (j = lo[1]; j <= hi[1]; j++) {
 
384
        jloc = j - lo[1]; 
 
385
        g(&gg, i, j, i, j);
 
386
        work[iloc][jloc] = sqrt(gg);
 
387
        //work[jloc][iloc] = sqrt(gg);
 
388
        *schwmax = MAX((*schwmax), work[iloc][jloc]);
 
389
        // *schwmax = MAX((*schwmax), work[jloc][iloc]);
 
390
      }
 
391
    }
 
392
 
 
393
    // lo_c[0]=lo[1];
 
394
    // lo_c[1]=lo[0];
 
395
    // hi_c[0]=hi[1];
 
396
    // hi_c[1]=hi[0];
 
397
 
 
398
    NGA_Put(g_schwarz, lo, hi, work, &ld);
 
399
    // NGA_Put(g_schwarz,lo_c,hi_c,work,&ld); //column-major inter-patch switch
 
400
    dotask = next_chunk(lo, hi);
 
401
  }
 
402
  GA_Dgop(schwmax, 1, "max");
 
403
 
 
404
  return;
 
405
}
 
406
 
 
407
void ininrm(void)  
 
408
{
 
409
  long long int maxint;
 
410
  //implicit declared variables in fortran
 
411
  int i;
 
412
      
 
413
  //     write a little welcome message;
 
414
  maxint = nbfn;
 
415
  maxint = pow(maxint, 4);
 
416
 
 
417
  if (GA_Nodeid()==0) {
 
418
    printf(" Example Direct Self Consistent Field Program \n");
 
419
    printf(" -------------------------------------------- \n\n");
 
420
    printf(" no. of atoms .............. %5d\n", natom);
 
421
    printf(" no. of occupied orbitals .. %5ld\n", nocc);
 
422
    printf(" no. of basis functions .... %5ld\n", nbfn);
 
423
    printf(" basis functions^4 ......... %5ld\n", maxint);
 
424
    printf(" convergence threshold ..... %9.4f\n", tol);
 
425
    printf(" chunk size .................%5d\n", ichunk);
 
426
  }
 
427
 
 
428
  //     generate normalisation coefficients for the basis functions;
 
429
  //     and the index array iky;
 
430
  for (i = 0; i < nbfn; i++) {
 
431
    //iky[i] = i*(i-1)/2;
 
432
    iky[i] = (i + 1) * i / 2;
 
433
  }
 
434
 
 
435
  for (i = 0;i < nbfn;i++) {
 
436
    rnorm[i] = pow((expnt[i] * 2.00 / pi), 0.750);
 
437
  }
 
438
 
 
439
  //     initialize common for computing f0;
 
440
  setfm();
 
441
}
 
442
 
 
443
double h(int i,int j) 
 
444
{
 
445
//vd$r novector;
 
446
//vd$r noconcur;
 
447
//     generate the one particle hamiltonian matrix element;
 
448
//     over the normalized primitive 1s functions i and j;
 
449
      double ret;
 
450
      
 
451
      //implicit declared variables in fortran
 
452
      double f0val = 0.00;
 
453
      double sum = 0.00;
 
454
      double facij,expij,repij;
 
455
      int iat;
 
456
      //long iat;
 
457
      double xp,yp,zp,rpc2;
 
458
      double rab2;
 
459
     
 
460
      rab2 = (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) + (z[i]-z[j])*(z[i]-z[j]);
 
461
      facij = expnt[i]*expnt[j]/(expnt[i]+expnt[j]);
 
462
      expij = exprjh(-facij*rab2);
 
463
      repij = (2.00*pi/(expnt[i]+expnt[j])) * expij;
 
464
 
 
465
//     first do the nuclear attraction integrals;
 
466
 
 
467
      for (iat = 0;iat < natom;iat++) 
 
468
      {
 
469
         xp = (x[i]*expnt[i] + x[j]*expnt[j])/(expnt[i]+expnt[j]);
 
470
         yp = (y[i]*expnt[i] + y[j]*expnt[j])/(expnt[i]+expnt[j]);
 
471
         zp = (z[i]*expnt[i] + z[j]*expnt[j])/(expnt[i]+expnt[j]);
 
472
         rpc2 = (xp-ax[iat])*(xp-ax[iat]) + (yp-ay[iat])*(yp-ay[iat]) + (zp-az[iat])*(zp-az[iat]);
 
473
//;
 
474
         f0(&f0val, (expnt[i]+expnt[j])*rpc2);
 
475
         sum = sum - repij * q[iat] * f0val;
 
476
        }
 
477
//     add on the kinetic energy term;
 
478
 
 
479
      sum = sum + facij*(3.00-2.00*facij*rab2) *
 
480
            pow((pi/(expnt[i]+expnt[j])),1.50) * expij;
 
481
 
 
482
//     finally multiply by the normalization constants;
 
483
      ret = sum * rnorm[i] * rnorm[j];
 
484
      return ret;
 
485
}
 
486
 
 
487
double s(int i, int j) 
 
488
{
 
489
  //     generate the overlap matrix element between the normalized;
 
490
  //     primitve gaussian 1s functions i and j;
 
491
  double ret;
 
492
 
 
493
  //implicit declared variables in fortran
 
494
  double rab2, facij;
 
495
 
 
496
  rab2 = (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]) +
 
497
    (z[i] - z[j]) * (z[i] - z[j]);
 
498
  facij = expnt[i] * expnt[j] / (expnt[i] + expnt[j]);
 
499
  ret = pow((pi / (expnt[i] + expnt[j])), 1.50) * exprjh(-facij * rab2) * rnorm[i] * rnorm[j];
 
500
  return ret;
 
501
}
 
502
 
 
503
void makden(void) 
 
504
{
 
505
  double work[ichunk][ichunk], orbsi[maxnbfn][ichunk]; //maxnbfn
 
506
  double orbsj[maxnbfn][ichunk]; //maxnbfn
 
507
  int lo[2], hi[2], tlo[2], thi[2], i, j, iloc, jloc, ld;
 
508
  int dotask;
 
509
 
 
510
  //implicit declared variables in fortran
 
511
  double p;
 
512
  int k;
 
513
 
 
514
  //transform between column-major row-major
 
515
  // double a[ichunk][ichunk], b[ichunk][ichunk];
 
516
  // int lo_c[2],hi_c[2]; 
 
517
 
 
518
  //     generate density matrix from orbitals in g_orbs. the first;
 
519
  //     nocc orbitals are doubly occupied.;
 
520
  // GA_Zero(g_counter);
 
521
  // dotask = next_chunk(lo, hi);
 
522
  ld = ichunk;
 
523
 
 
524
  //transform g_orbs from column-major to row-major 
 
525
  /*   while (dotask) { */
 
526
  /*     NGA_Get(g_orbs, lo, hi, a, &ld); */
 
527
  /*     for (j = 0; j <= hi[1] - lo[1]; j++) { */
 
528
  /*       for (i = 0; i <= hi[0] - lo[0]; i++) { */
 
529
  /*    b[j][i] = a[i][j]; //intra-patch transfrom */
 
530
  /*       } */
 
531
  /*     } */
 
532
  /*     lo_c[0] = lo[1]; //inter-patch transform */
 
533
  /*     lo_c[1] = lo[0]; */
 
534
  /*     hi_c[0] = hi[1]; */
 
535
  /*     hi_c[1] = hi[0]; */
 
536
  /*     NGA_Put(g_tmp, lo_c, hi_c, b, &ld); */
 
537
  /*     dotask = next_chunk(lo, hi); */
 
538
  /*   } */
 
539
  /*   GA_Copy(g_tmp, g_orbs); */
 
540
 
 
541
  // GA_Transpose(g_orbs, g_tmp);
 
542
  // GA_Copy(g_tmp, g_orbs);
 
543
      
 
544
  GA_Zero(g_counter);
 
545
  dotask = next_chunk(lo, hi);
 
546
  while (dotask) {
 
547
    ld = ichunk;
 
548
 
 
549
    tlo[0] = 0;
 
550
    thi[0] = nocc - 1;
 
551
 
 
552
    tlo[1] = lo[0];
 
553
    thi[1] = hi[0];
 
554
    NGA_Get(g_orbs, tlo, thi, orbsi, &ld);
 
555
        
 
556
    tlo[1] = lo[1];
 
557
    thi[1] = hi[1];
 
558
    NGA_Get(g_orbs, tlo, thi, orbsj, &ld);
 
559
 
 
560
    for (i = lo[0]; i <= hi[0]; i++) {
 
561
      iloc = i - lo[0];// + 1;
 
562
      for (j = lo[1]; j <= hi[1]; j++) {
 
563
        jloc = j - lo[1];// + 1;
 
564
        p = 0.0;
 
565
        for (k = 0; k < nocc; k++) {
 
566
          p = p + orbsi[k][iloc] * orbsj[k][jloc];
 
567
        }
 
568
        work[iloc][jloc] = 2.0 * p;
 
569
      }
 
570
    }
 
571
    ld = ichunk;
 
572
    NGA_Put(g_work, lo, hi, work, &ld);
 
573
    dotask = next_chunk(lo, hi);
 
574
  }
 
575
 
 
576
  // GA_Transpose(g_orbs, g_tmp);
 
577
  // GA_Copy(g_tmp, g_orbs);
 
578
 
 
579
  // GA_Transpose(g_work, g_tmp);
 
580
  // GA_Copy(g_tmp, g_work);
 
581
 
 
582
  //transform g_orbs from row-major to column major
 
583
  /*   ld = ichunk; */
 
584
  /*   GA_Zero(g_counter); */
 
585
  /*   dotask = next_chunk(lo, hi); */
 
586
  /*   while (dotask) { */
 
587
  /*     NGA_Get(g_orbs, lo, hi, a, &ld); */
 
588
  /*     for (j = 0; j <= hi[1] - lo[1]; j++) { */
 
589
  /*       for (i = 0; i <= hi[0] - lo[0]; i++) { */
 
590
  /*    b[j][i] = a[i][j]; //intra-patch transfrom */
 
591
  /*       } */
 
592
  /*     } */
 
593
  /*     lo_c[0] = lo[1]; //inter-patch transform */
 
594
  /*     lo_c[1] = lo[0]; */
 
595
  /*     hi_c[0] = hi[1]; */
 
596
  /*     hi_c[1] = hi[0]; */
 
597
  /*     NGA_Put(g_tmp, lo_c, hi_c, b, &ld); */
 
598
  /*     dotask = next_chunk(lo, hi); */
 
599
  /*   } */
 
600
  /*   GA_Copy(g_tmp, g_orbs); */
 
601
  // transform g_work from row-major to column major
 
602
  /*   ld = ichunk; */
 
603
  /*   GA_Zero(g_counter); */
 
604
  /*   dotask = next_chunk(lo, hi); */
 
605
  /*   while (dotask) { */
 
606
  /*     NGA_Get(g_work, lo, hi, a, &ld); */
 
607
  /*     for (j = 0; j <= hi[1] - lo[1]; j++) { */
 
608
  /*       for (i = 0; i <= hi[0] - lo[0]; i++) { */
 
609
  /*    b[j][i] = a[i][j]; //intra-patch transfrom */
 
610
  /*       } */
 
611
  /*     } */
 
612
  /*     lo_c[0] = lo[1]; //inter-patch transform */
 
613
  /*     lo_c[1] = lo[0]; */
 
614
  /*     hi_c[0] = hi[1]; */
 
615
  /*     hi_c[1] = hi[0]; */
 
616
  /*     NGA_Put(g_tmp, lo_c, hi_c, b, &ld); */
 
617
  /*     dotask = next_chunk(lo, hi); */
 
618
  /*   } */
 
619
  /*   GA_Copy(g_tmp, g_work); */
 
620
  return;
 
621
}
 
622
 
 
623
int nxtask(int nproc) //not used
 
624
{
 
625
      const int ichunk_local = 10; //in the header file, ichunk is defined as 20, so I rename it as ichunk_local
 
626
      static double nleft=0, icount=0;
 
627
      
 
628
      int ret = 0;
 
629
 
 
630
      //implicit declared variables in fortran
 
631
      int junk;
 
632
 
 
633
//     wrapper round nxtval() to increase granularity;
 
634
//     and thus reduce no. of requests to shared counter;
 
635
#ifndef MPI
 
636
      if(nproc>0) 
 
637
      {
 
638
         if(nleft==0) 
 
639
         {
 
640
            icount = tcg_nxtval(nproc) * ichunk_local;
 
641
            nleft = ichunk_local;
 
642
         }
 
643
         ret = icount;
 
644
         icount = icount + 1;
 
645
         nleft = nleft -1;
 
646
      } 
 
647
      else 
 
648
      {
 
649
          nleft = 0;
 
650
          ret = 0;
 
651
          junk = tcg_nxtval(nproc);
 
652
      }
 
653
#endif
 
654
//     following does dumb static load balancing;
 
655
//;
 
656
//$$$      if(nproc>0) {
 
657
//$$$         if (nleft == 0) {
 
658
//$$$            icount = GA_Nodeid();
 
659
//$$$            nleft = 1;
 
660
//$$$         }
 
661
//$$$         nxtask = icount;
 
662
//$$$         icount = icount + GA_Nnodes();
 
663
//$$$      } else {
 
664
//$$$          nleft = 0;
 
665
//$$$          nxtask = 0;
 
666
//$$$      }
 
667
 
 
668
     return ret;
 
669
}
 
670
 
 
671
int next_chunk(int* lo,int* hi) 
 
672
{
 
673
  //const int one=1;
 
674
  int one = 0;
 
675
  int imax, ilo, jlo;
 
676
 
 
677
  //implicit declared variables in fortran
 
678
  int itask;   
 
679
      
 
680
  int ret;   
 
681
 
 
682
  itask = NGA_Read_inc(g_counter, &one, 1);
 
683
  imax = nbfn / ichunk;
 
684
  if (nbfn - ichunk * imax > 0)
 
685
    imax = imax + 1;
 
686
  if (itask < imax * imax) {
 
687
    ilo = itask % imax;
 
688
    jlo = (itask - ilo) / imax;
 
689
    lo[0] = ilo * ichunk; 
 
690
    lo[1] = jlo * ichunk; 
 
691
    hi[0] = MIN((ilo + 1) * ichunk - 1, nbfn); 
 
692
    hi[1] = MIN((jlo + 1) * ichunk - 1, nbfn); 
 
693
    ret = 1;
 
694
  } 
 
695
  else 
 
696
    ret = 0;
 
697
 
 
698
  return ret;
 
699
}
 
700
 
 
701
long int acquire_tasks(int numTasks)
 
702
{
 
703
  int one = 0;
 
704
  long int itask;
 
705
 
 
706
  itask = NGA_Read_inc(g_counter, &one, numTasks);
 
707
 
 
708
  return itask;
 
709
} // acquire_tasks
 
710
 
 
711
int translate_task(long int itask, int *lo, int *hi, int *ilo, int *jlo, int *klo, int *llo)
 
712
{
 
713
  long int imax;
 
714
  int itmp;
 
715
  int ret;
 
716
 
 
717
  imax = nbfn / ichunk;
 
718
  if (nbfn - ichunk * imax > 0)
 
719
    imax = imax + 1;
 
720
  if (itask < 0) {
 
721
    printf("translate_task: itask negative: %ld imax: %ld nbfn: %ld ichunk: %d\n", itask, imax,
 
722
           nbfn, ichunk);
 
723
    printf("probable GA int precision problem if imax^4 > 2^31\n");
 
724
    printf("\n"); 
 
725
    printf("translate_task\n");
 
726
    exit(0);
 
727
  }
 
728
  if (itask < pow(imax, 4)) {
 
729
    *ilo = itask % imax;
 
730
    itmp = (itask - (*ilo)) / imax;
 
731
    *jlo = itmp % imax;
 
732
    itmp = (itmp - (*jlo)) / imax;
 
733
    *klo = itmp % imax;
 
734
    *llo = (itmp - (*klo)) / imax;
 
735
    lo[0] = (*ilo) * ichunk;
 
736
    lo[1] = (*jlo) * ichunk;
 
737
    lo[2] = (*klo) * ichunk;
 
738
    lo[3] = (*llo) * ichunk;
 
739
    hi[0] = MIN(((*ilo) + 1) * ichunk - 1, nbfn);
 
740
    hi[1] = MIN(((*jlo) + 1) * ichunk - 1, nbfn);
 
741
    hi[2] = MIN(((*klo) + 1) * ichunk  -1, nbfn);
 
742
    hi[3] = MIN(((*llo) + 1) * ichunk - 1, nbfn);
 
743
    ret = 1;
 
744
  }
 
745
  else
 
746
    ret = 0;
 
747
 
 
748
  return ret;
 
749
} // translate_task
 
750
 
 
751
int next_4chunk(int *lo,int *hi,int* ilo,int* jlo,int* klo,int* llo, long int *pitask)
 
752
{
 
753
  int one = 0;
 
754
  long int imax;
 
755
  long int itask;
 
756
  int itmp;
 
757
 
 
758
  int ret;
 
759
 
 
760
  itask = NGA_Read_inc(g_counter, &one, 1);
 
761
  *pitask = itask;
 
762
 
 
763
  imax = nbfn / ichunk;
 
764
  if (nbfn - ichunk * imax > 0)
 
765
    imax = imax + 1;
 
766
  if (itask < 0) {
 
767
    printf("next_4chunk: itask negative: %ld imax: %ld nbfn: %ld ichunk: %d\n", itask, imax,
 
768
           nbfn, ichunk);
 
769
    printf("probable GA int precision problem if imax^4 > 2^31\n");
 
770
    printf("\n"); 
 
771
    printf("next_4chunk\n");
 
772
    exit(0);
 
773
  }
 
774
  if (itask < pow(imax, 4)) {
 
775
    *ilo = itask % imax;
 
776
    itmp = (itask - (*ilo)) / imax;
 
777
    *jlo = itmp % imax;
 
778
    itmp = (itmp - (*jlo)) / imax;
 
779
    *klo = itmp % imax;
 
780
    *llo = (itmp - (*klo)) / imax;
 
781
    lo[0] = (*ilo) * ichunk;
 
782
    lo[1] = (*jlo) * ichunk;
 
783
    lo[2] = (*klo) * ichunk;
 
784
    lo[3] = (*llo) * ichunk;
 
785
    hi[0] = MIN(((*ilo) + 1) * ichunk - 1, nbfn);
 
786
    hi[1] = MIN(((*jlo) + 1) * ichunk - 1, nbfn);
 
787
    hi[2] = MIN(((*klo) + 1) * ichunk  -1, nbfn);
 
788
    hi[3] = MIN(((*llo) + 1) * ichunk - 1, nbfn);
 
789
    ret = 1;
 
790
  }
 
791
  else
 
792
    ret = 0;
 
793
 
 
794
  return ret;
 
795
}
 
796
 
 
797
void clean_chunk(double chunk[][ichunk])
 
798
{
 
799
  int i, j;
 
800
 
 
801
  for (i = 0; i < ichunk; i++)
 
802
    for (j = 0; j < ichunk; j++)
 
803
      chunk[i][j] = 0.0;
 
804
} // clean_chunk
 
805
 
 
806
int verify_zero_chunk(double chunk[][ichunk])
 
807
{
 
808
  int i, j, flag;
 
809
 
 
810
  flag = 1;
 
811
 
 
812
  for (i = 0; (i < ichunk) && flag; i++)
 
813
    for (j = 0; (j < ichunk) && flag; j++)
 
814
      if (chunk[i][j] != 0.0) {
 
815
        flag = 0;
 
816
        fprintf(stderr, "chunk[%d][%d] != 0.0, %.16f\n", i, j, chunk[i][j]);
 
817
        break;
 
818
      }
 
819
 
 
820
  return flag;
 
821
}
 
822
 
 
823
void damp(double *fac)
 
824
{
 
825
  //    create damped density matrix as a linear combination of;
 
826
  //    old density matrix and density matrix formed from new orbitals;
 
827
      
 
828
  //implicit declared variables in fortran
 
829
  double ofac;
 
830
 
 
831
  ofac = 1.00 - (*fac);
 
832
  GA_Add(fac, g_dens, &ofac, g_work, g_dens);
 
833
  return;
 
834
}
 
835
 
 
836
void prnout(int iter, double energy, double deltad, double tester) 
 
837
{
 
838
  //     printout results of each iteration;
 
839
  if (GA_Nodeid() != 0)
 
840
    return;
 
841
 
 
842
  printf(" iter= %3d, energy=%15.8f, deltad= %9.7f, deltaf=%9.7f\n",
 
843
         iter, energy, deltad, tester);
 
844
  return;
 
845
}
 
846
 
 
847
double dendif(void)
 
848
{
 
849
  double xdiff;
 
850
  double dens_c[ichunk][ichunk],work_c[ichunk][ichunk];
 
851
  int lo[2], hi[2], i, j, ld;
 
852
  int dotask;
 
853
 
 
854
  // int lo_c[2],hi_c[2]; //column-major inter-patch switch
 
855
 
 
856
  double ret;
 
857
 
 
858
  //implicit declared variables in fortran
 
859
  double denmax;
 
860
  //     compute largest change in density matrix elements;
 
861
  denmax = 0.00;
 
862
  GA_Zero(g_counter);
 
863
  dotask = next_chunk(lo,hi);
 
864
  ld = ichunk;
 
865
  while(dotask) {
 
866
    //column-major inter-path switch
 
867
    // lo_c[0] = lo[1];
 
868
    // lo_c[1] = lo[0];
 
869
    // hi_c[0] = hi[1];
 
870
    // hi_c[1] = hi[0];
 
871
        
 
872
    NGA_Get(g_dens, lo, hi, dens_c, &ld);
 
873
    NGA_Get(g_work, lo, hi, work_c, &ld);
 
874
 
 
875
    //column-major inter-patch switch        
 
876
    // NGA_Get(g_dens,lo_c,hi_c,dens_c,&ld);
 
877
    // NGA_Get(g_work,lo_c,hi_c,work_c,&ld);
 
878
        
 
879
    for (j = 0; j <= hi[1] - lo[1]; j++) {
 
880
      for (i = 0; i <= hi[0] - lo[0]; i++) {
 
881
        //xdiff = abs(dens_c[i][j]-work_c[i][j]); 
 
882
        //xdiff = abs_double(dens_c[j][i]-work_c[j][i]); //column-major intra-patch switch 
 
883
        xdiff = fabs(dens_c[i][j] - work_c[i][j]); //column-major intra-patch switch 
 
884
        if (xdiff > denmax)
 
885
          denmax = xdiff;
 
886
      }
 
887
    }
 
888
    dotask = next_chunk(lo, hi);
 
889
  }
 
890
  GA_Dgop(&denmax, 1, "max");
 
891
  ret = denmax;
 
892
  return ret;
 
893
}
 
894
 
 
895
double testfock(void) 
 
896
{
 
897
  double xmax, xtmp;
 
898
  double work[ichunk][ichunk];
 
899
  int lo[2], hi[2], i, j, iloc, jloc, ld;
 
900
  int dotask;
 
901
 
 
902
  // int lo_c[2],hi_c[2]; //column-major inter-patch switch
 
903
 
 
904
  double ret;
 
905
  //     compute largest change in density matrix elements;
 
906
  xmax = 0.00;
 
907
  GA_Zero(g_counter);
 
908
  dotask = next_chunk(lo, hi);
 
909
  ld = ichunk;
 
910
  while(dotask) {
 
911
    //column-major inter-patch switch
 
912
    // lo_c[0]=lo[1];
 
913
    // lo_c[1]=lo[0];
 
914
    // hi_c[0]=hi[1];
 
915
    // hi_c[1]=hi[0];
 
916
    
 
917
    NGA_Get(g_fock, lo, hi, work, &ld); 
 
918
    // NGA_Get(g_fock, lo_c, hi_c, work, &ld);//column-major inter-patch switch
 
919
 
 
920
    for (j = lo[1]; j <= hi[1]; j++) {
 
921
      jloc = j - lo[1];
 
922
      for (i = lo[0]; i <= hi[0]; i++) {
 
923
        iloc = i - lo[0];
 
924
        if (i != j) {
 
925
          //xtmp = abs(work[iloc][jloc]);
 
926
          // xtmp = abs_double(work[jloc][iloc]); //column-major intra-patch switch
 
927
          xtmp = fabs(work[iloc][jloc]); //column-major intra-patch switch
 
928
          if (xtmp > xmax)
 
929
            xmax = xtmp;
 
930
        }
 
931
      }
 
932
    }
 
933
    dotask = next_chunk(lo, hi);
 
934
  }
 
935
  GA_Dgop(&xmax, 1, "max");
 
936
  ret = xmax;
 
937
 
 
938
  return ret;
 
939
}
 
940
 
 
941
void shiftfock(double shift)
 
942
{
 
943
  double work[ichunk][ichunk];
 
944
  int lo[2], hi[2], i, j, iloc, jloc, ld, icnt;
 
945
  int dotask;
 
946
 
 
947
  // int lo_c[2],hi_c[2]; //column-major inter-patch switch
 
948
 
 
949
  //     compute largest change in density matrix elements;
 
950
  GA_Zero(g_counter);
 
951
  dotask = next_chunk(lo, hi);
 
952
  ld = ichunk;
 
953
  while (dotask) {
 
954
    //column-major inter-path switch
 
955
    // lo_c[0]=lo[1];
 
956
    // lo_c[1]=lo[0];
 
957
    // hi_c[0]=hi[1];
 
958
    // hi_c[1]=hi[0];
 
959
 
 
960
    // NGA_Get(g_fock, lo_c, hi_c, work, &ld); //column-major inter-patch switch
 
961
    NGA_Get(g_fock, lo, hi, work, &ld);
 
962
    icnt = 0;
 
963
    for (j = lo[1]; j <= hi[1]; j++) {
 
964
      jloc = j - lo[1];
 
965
      for (i = lo[0]; i <= hi[0]; i++) {
 
966
        iloc = i - lo[0];
 
967
        if ((i == j) && (i > (nocc - 1))) {
 
968
          work[iloc][jloc] = work[iloc][jloc] + shift;
 
969
          // work[jloc][iloc] = work[jloc][iloc] + shift; //column-major inter-patch switch
 
970
          icnt = icnt + 1;
 
971
        }
 
972
      }
 
973
    }
 
974
    if (icnt > 0)
 
975
      //NGA_Put(g_fock, lo_c, hi_c, work, &ld);//column-major inter-patch switch
 
976
      NGA_Put(g_fock, lo, hi, work, &ld);
 
977
    dotask = next_chunk(lo, hi);
 
978
  }
 
979
  return;
 
980
}
 
981
 
 
982
void prnfin(double energy) 
 
983
{
 
984
  double orbs[maxnbfn][maxnbfn];
 
985
  int lo[2],hi[2],ld;
 
986
 
 
987
  //     printout final results;
 
988
  if (GA_Nodeid() != 0)
 
989
    return;
 
990
 
 
991
  printf("\n\nfinal energy = %18.11f\n",energy);
 
992
  printf("\neigenvalues\n\n");
 
993
  output(eigv, 0, MIN(nbfn,nocc+5), 0, 1, nbfn, 1, 1);
 
994
  //output(eigv, 1, MIN(nbfn,nocc+5), 1, 1, nbfn, 1, 1);
 
995
 
 
996
  return;
 
997
}
 
998
 
 
999
void g(double *value, int i, int j, int k, int l)
 
1000
{
 
1001
  //implicit declared variables in fortran
 
1002
  double f0val, rab2, rcd2, facij,fackl,exijkl,denom, fac, xp,yp,zp,xq,yq,zq,rpq2;
 
1003
 
 
1004
  //     compute the two electon integral (ij|kl) over normalized;
 
1005
  //     primitive 1s gaussians;
 
1006
  f0val = 0.00;
 
1007
  rab2 = (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]) +
 
1008
    (z[i] - z[j]) * (z[i] - z[j]);
 
1009
  rcd2 = (x[k] - x[l]) * (x[k] - x[l]) + (y[k] - y[l]) * (y[k] - y[l]) +
 
1010
    (z[k] - z[l]) * (z[k] - z[l]);
 
1011
  facij = expnt[i] * expnt[j] / (expnt[i] + expnt[j]);
 
1012
  fackl = expnt[k] * expnt[l] / (expnt[k] + expnt[l]);
 
1013
  exijkl = exprjh(-facij * rab2 - fackl * rcd2);
 
1014
  denom = (expnt[i] + expnt[j]) * (expnt[k] + expnt[l]) *
 
1015
    sqrt(expnt[i] + expnt[j] + expnt[k] + expnt[l]);
 
1016
  fac = (expnt[i] + expnt[j]) * (expnt[k] + expnt[l]) /
 
1017
    (expnt[i] + expnt[j] + expnt[k] + expnt[l]);
 
1018
 
 
1019
  xp = (x[i] * expnt[i] + x[j] * expnt[j]) / (expnt[i] + expnt[j]);
 
1020
  yp = (y[i] * expnt[i] + y[j] * expnt[j]) / (expnt[i] + expnt[j]);
 
1021
  zp = (z[i] * expnt[i] + z[j] * expnt[j]) / (expnt[i] + expnt[j]);
 
1022
  xq = (x[k] * expnt[k] + x[l] * expnt[l]) / (expnt[k] + expnt[l]);
 
1023
  yq = (y[k] * expnt[k] + y[l] * expnt[l]) / (expnt[k] + expnt[l]);
 
1024
  zq = (z[k] * expnt[k] + z[l] * expnt[l]) / (expnt[k] + expnt[l]);
 
1025
  rpq2 = (xp - xq) * (xp - xq) + (yp - yq) * (yp - yq) + (zp - zq) * (zp - zq);
 
1026
 
 
1027
  f0(&f0val, fac * rpq2);
 
1028
  *value = (2.00 * pow(pi, 2.50) / denom) * exijkl * f0val *
 
1029
    rnorm[i] * rnorm[j] * rnorm[k] * rnorm[l];
 
1030
  return;
 
1031
}
 
1032
 
 
1033
void diagon(double *tester, int iter)
 
1034
{
 
1035
  //      diagon(fock, orbs, evals, work, tester, iter);
 
1036
  double r_zero, r_one, shift;
 
1037
 
 
1038
  //implicit declared variables in fortran
 
1039
  int i, me;
 
1040
 
 
1041
  double test = -1.0, s;
 
1042
 
 
1043
  me = GA_Nodeid();
 
1044
 
 
1045
#if USE_TRANSFORM
 
1046
  //     use similarity transform to solve standard eigenvalue problem;
 
1047
  //     (overlap matrix has been transformed out of the problem);
 
1048
  r_one = 1.00;
 
1049
  r_zero = 0.00;
 
1050
  GA_Dgemm('n', 'n', nbfn, nbfn, nbfn, r_one, g_fock, g_orbs, r_zero, g_tfock);
 
1051
 
 
1052
#if 0
 
1053
  s = sum(g_tfock);
 
1054
  GA_Dgop(&s, 1, "+");
 
1055
 
 
1056
  if (me == 0)
 
1057
    printf("g_tfock sum after first dgemm(): %f\n", s);
 
1058
#endif
 
1059
 
 
1060
  GA_Dgemm('t', 'n', nbfn, nbfn, nbfn, r_one, g_orbs, g_tfock, r_zero, g_fock);
 
1061
 
 
1062
#if 0
 
1063
  s = sum(g_fock);
 
1064
  GA_Dgop(&s, 1, "+");
 
1065
 
 
1066
  if (me == 0)
 
1067
    printf("g_fock sum after second dgemm(): %f\n", s);
 
1068
#endif
 
1069
 
 
1070
  *tester = testfock();
 
1071
  shift = 0.00;
 
1072
  if ((*tester) > 0.30) {
 
1073
    shift = 0.30;
 
1074
  }
 
1075
#if 1
 
1076
  else {
 
1077
    if (nbfn > 60)
 
1078
      shift = 0.10;
 
1079
    else
 
1080
      shift = 0.00;
 
1081
  }
 
1082
#endif
 
1083
      
 
1084
  //if (iter>=2 && shift!=0.00) 
 
1085
  if (iter >= 1 && shift != 0.00) { //iter 2 in Fotran is iter 1 in C
 
1086
    shiftfock(shift);
 
1087
  }
 
1088
  GA_Copy(g_orbs, g_tfock);
 
1089
  GA_Diag_std_seq(g_fock, g_work, eigv);
 
1090
 
 
1091
  //     Back transform eigenvectors;
 
1092
  GA_Dgemm('n', 'n', nbfn, nbfn, nbfn, r_one, g_tfock, g_work, r_zero, g_orbs);
 
1093
 
 
1094
  if (iter>= 1 && shift != 0.00) { //>=2 --> >=1
 
1095
    for (i = nocc; i < nbfn; i++) {
 
1096
      eigv[i] = eigv[i] - shift;
 
1097
    }
 
1098
  }
 
1099
#else
 
1100
  //;
 
1101
  //     Keep remaking overlap matrix since GA_Diag_seq does not;
 
1102
  //     guarantee that g_ident is preserved.;
 
1103
  //;
 
1104
  makoverlap();
 
1105
  GA_Diag_seq(g_fock, g_ident, g_orbs, eigv); 
 
1106
  *tester = 0.00;
 
1107
#endif
 
1108
  return;
 
1109
}
 
1110
 
 
1111
void makeob(void) 
 
1112
{
 
1113
  double work[ichunk][ichunk],orbs[ichunk][ichunk];
 
1114
  double eval[maxnbfn];
 
1115
  int lo[2],hi[2],ld,me,i,j,iloc,jloc;
 
1116
  int dotask;
 
1117
 
 
1118
  int lo_c[2],hi_c[2]; //for column-major switch 
 
1119
  //     generate set of orthonormal vectors by creating a random;
 
1120
  //     symmetric matrix and solving associated generalized eigenvalue;
 
1121
  //     problem using the correct overlap matrix.;
 
1122
 
 
1123
  me = GA_Nodeid();
 
1124
  GA_Zero(g_counter);
 
1125
  dotask = next_chunk(lo, hi);
 
1126
  ld = ichunk;
 
1127
 
 
1128
  while (dotask) {
 
1129
    for (i = lo[0]; i <= hi[0]; i++) {
 
1130
      iloc = i - lo[0];
 
1131
      for (j = lo[1]; j <= hi[1]; j++) {
 
1132
        jloc = j - lo[1];
 
1133
        work[iloc][jloc] = s(i, j);
 
1134
        //work[jloc][iloc] = s(i,j); //column-major intra-patch switch
 
1135
        orbs[iloc][jloc] = 0.5;//drand48();
 
1136
        //orbs[jloc][iloc] = 0.5; //column-major intra-patch switch
 
1137
      }
 
1138
    }
 
1139
 
 
1140
    NGA_Put(g_ident, lo, hi, work, &ld);
 
1141
    NGA_Put(g_fock, lo, hi, orbs, &ld);
 
1142
    // lo_c[0]=lo[1];
 
1143
    //lo_c[1]=lo[0];
 
1144
    // hi_c[0]=hi[1];
 
1145
    // hi_c[1]=hi[0];
 
1146
    // NGA_Put(g_ident,lo_c,hi_c,work,&ld);
 
1147
    // NGA_Put(g_fock,lo_c,hi_c,orbs,&ld);
 
1148
 
 
1149
    dotask = next_chunk(lo, hi);
 
1150
  }
 
1151
  GA_Symmetrize(g_fock);
 
1152
  GA_Diag_seq(g_fock, g_ident, g_orbs, eval);
 
1153
  return;
 
1154
}
 
1155
 
 
1156
void denges(void)
 
1157
{
 
1158
  //     Form guess density from superposition of atomic densities in the AO;
 
1159
  //     basis set ... instead of doing the atomic SCF hardwire for this;
 
1160
  //     small basis set for the Be atom.;
 
1161
  int one, itask, lo[2], hi[2], ld;
 
1162
 
 
1163
  double atdens[15][15] = {
 
1164
    {0.000002,0.000027,0.000129,0.000428,0.000950,0.001180,
 
1165
     0.000457,-0.000270,-0.000271,0.000004,0.000004,0.000004,
 
1166
     0.000004,0.000004,0.000004},
 
1167
    {0.000027,0.000102,0.000987,
 
1168
     0.003269,0.007254,0.009007,0.003492,-0.002099,-0.002108,
 
1169
     0.000035,0.000035,0.000035,0.000035,0.000035,0.000035},
 
1170
    {0.000129,0.000987,0.002381,0.015766,0.034988,0.043433,
 
1171
     0.016835,-0.010038,-0.010082,0.000166,0.000166,0.000166,
 
1172
     0.000166,0.000166,0.000166},
 
1173
    {0.000428,0.003269,0.015766,
 
1174
     0.026100,0.115858,0.144064,0.055967,-0.035878,-0.035990,
 
1175
     0.000584,0.000584,0.000584,0.000584,0.000584,0.000584},
 
1176
    {0.000950,0.007254,0.034988,0.115858,0.128586,0.320120,
 
1177
     0.124539,-0.083334,-0.083536,0.001346,0.001346,0.001346,
 
1178
     0.001346,0.001346,0.001346},
 
1179
    {0.001180,0.009007,0.043433,
 
1180
     0.144064,0.320120,0.201952,0.159935,-0.162762,-0.162267,
 
1181
     0.002471,0.002471,0.002471,0.002471,0.002471,0.002471},
 
1182
    {0.000457,0.003492,0.016835,0.055967,0.124539,0.159935,
 
1183
     0.032378,-0.093780,-0.093202,0.001372,0.001372,0.001372,
 
1184
     0.001372,0.001372,0.001372},
 
1185
    {-0.000270,-0.002099,-0.010038,
 
1186
     -0.035878,-0.083334,-0.162762,-0.093780,0.334488,0.660918,
 
1187
     -0.009090,-0.009090,-0.009090,-0.009090,-0.009090,-0.009090},
 
1188
    {-0.000271,-0.002108,-0.010082,-0.035990,-0.083536,-0.162267,
 
1189
     -0.093202,0.660918,0.326482,-0.008982,-0.008982,-0.008981,
 
1190
     -0.008981,-0.008981,-0.008982},
 
1191
    {0.000004,0.000035,0.000166,
 
1192
     0.000584,0.001346,0.002471,0.001372,-0.009090,-0.008982,
 
1193
     0.000062,0.000124,0.000124,0.000124,0.000124,0.000124},
 
1194
    {0.000004,0.000035,0.000166,0.000584,0.001346,0.002471,
 
1195
     0.001372,-0.009090,-0.008982,0.000124,0.000062,0.000124,
 
1196
     0.000124,0.000124,0.000124},
 
1197
    {0.000004,0.000035,0.000166,
 
1198
     0.000584,0.001346,0.002471,0.001372,-0.009090,-0.008981,
 
1199
     0.000124,0.000124,0.000062,0.000124,0.000124,0.000124},
 
1200
    {0.000004,0.000035,0.000166,0.000584,0.001346,0.002471,
 
1201
     0.001372,-0.009090,-0.008981,0.000124,0.000124,0.000124,
 
1202
     0.000062,0.000124,0.000124},
 
1203
    {0.000004,0.000035,0.000166,
 
1204
     0.000584,0.001346,0.002471,0.001372,-0.009090,-0.008981,
 
1205
     0.000124,0.000124,0.000124,0.000124,0.000062,0.000124},
 
1206
    {0.000004,0.000035,0.000166,0.000584,0.001346,0.002471,
 
1207
     0.001372,-0.009090,-0.008982,0.000124,0.000124,0.000124,
 
1208
     0.000124,0.000124,0.000062}};
 
1209
      
 
1210
  //implicit declared variables in fortran
 
1211
  int ioff,i,j;
 
1212
      
 
1213
  double atdens_f[15][15]; 
 
1214
 
 
1215
  //   Create initial guess for density matrix in global array;
 
1216
  GA_Zero(g_dens);
 
1217
  GA_Zero(g_counter);
 
1218
  //one = 1;
 
1219
  one = 0;
 
1220
  ld = 15;
 
1221
 
 
1222
  //   Correct for a factor of two along the diagonal;
 
1223
  for (i = 0; i < ld; i++) {
 
1224
    atdens[i][i] = 2.00 * atdens[i][i];
 
1225
  }
 
1226
 
 
1227
  //to get a fortran based array by switching the indices
 
1228
  for (i = 0; i < 15; i++) {
 
1229
    for (j = 0; j < 15; j++) {
 
1230
      atdens_f[j][i] = atdens[i][j];
 
1231
    }
 
1232
  }
 
1233
 
 
1234
  itask = NGA_Read_inc(g_counter, &one, 1);
 
1235
      
 
1236
  while(itask < natom) {
 
1237
    ioff = itask * 15;
 
1238
    lo[0] = ioff + 0; //1--->0
 
1239
    lo[1] = ioff + 0; //1--->0
 
1240
    hi[0] = ioff + 14;//15--->14
 
1241
    hi[1] = ioff + 14;//15--->14 
 
1242
    //NGA_Put(g_dens,lo,hi,atdens,&ld);
 
1243
    NGA_Put(g_dens, lo, hi, atdens_f, &ld);
 
1244
    itask = NGA_Read_inc(g_counter, &one, 1);
 
1245
  }
 
1246
  GA_Sync();
 
1247
  return;
 
1248
}
 
1249
 
 
1250
void setarrays(void)
 
1251
{
 
1252
  int one, two, dims[2];
 
1253
  int status;
 
1254
  int nproc;
 
1255
 
 
1256
  one = 1;
 
1257
  two = 2;
 
1258
 
 
1259
  nproc = GA_Nnodes();
 
1260
 
 
1261
  g_counter = GA_Create_handle();
 
1262
  GA_Set_data(g_counter,one,&one,MT_INT);
 
1263
  status = GA_Allocate(g_counter);
 
1264
  GA_Zero(g_counter);
 
1265
  
 
1266
  g_proc = GA_Create_handle();
 
1267
  GA_Set_data(g_proc, one, &nproc, MT_INT);
 
1268
  status = GA_Allocate(g_proc);
 
1269
  GA_Zero(g_proc);
 
1270
      
 
1271
  dims[0] = nbfn;
 
1272
  dims[1] = nbfn;
 
1273
  g_dens = GA_Create_handle();
 
1274
  GA_Set_data(g_dens, two, dims, MT_DBL);
 
1275
  status = GA_Allocate(g_dens);
 
1276
  GA_Zero(g_dens);
 
1277
 
 
1278
  g_schwarz = GA_Create_handle();
 
1279
  GA_Set_data(g_schwarz, two, dims, MT_DBL);
 
1280
  status = GA_Allocate(g_schwarz);
 
1281
  GA_Zero(g_schwarz);
 
1282
 
 
1283
  g_fock = GA_Create_handle();
 
1284
  GA_Set_data(g_fock, two, dims, MT_DBL);
 
1285
  status = GA_Allocate(g_fock);
 
1286
  GA_Zero(g_fock);
 
1287
 
 
1288
  g_tfock = GA_Create_handle();
 
1289
  GA_Set_data(g_tfock, two, dims, MT_DBL);
 
1290
  status = GA_Allocate(g_tfock);
 
1291
  GA_Zero(g_tfock);
 
1292
 
 
1293
  g_work = GA_Create_handle();
 
1294
  GA_Set_data(g_work, two, dims, MT_DBL);
 
1295
  status = GA_Allocate(g_work);
 
1296
  GA_Zero(g_work);
 
1297
 
 
1298
  g_ident = GA_Create_handle();
 
1299
  GA_Set_data(g_ident, two, dims, MT_DBL);
 
1300
  status = GA_Allocate(g_ident);
 
1301
  GA_Zero(g_ident);
 
1302
 
 
1303
  g_orbs = GA_Create_handle();
 
1304
  GA_Set_data(g_orbs, two, dims, MT_DBL);
 
1305
  status = GA_Allocate(g_orbs);
 
1306
  GA_Zero(g_orbs);
 
1307
 
 
1308
  //temporay global array for storage major transformation
 
1309
  g_tmp = GA_Create_handle();
 
1310
  GA_Set_data(g_tmp, two, dims, MT_DBL);
 
1311
  status = GA_Allocate(g_tmp);
 
1312
  GA_Zero(g_tmp);
 
1313
 
 
1314
  return;
 
1315
}
 
1316
 
 
1317
void closearrays(void)
 
1318
{
 
1319
  GA_Destroy(g_counter);
 
1320
  GA_Destroy(g_proc);
 
1321
  GA_Destroy(g_dens);
 
1322
  GA_Destroy(g_schwarz);
 
1323
  GA_Destroy(g_fock);
 
1324
  GA_Destroy(g_tfock);
 
1325
  GA_Destroy(g_work);
 
1326
  GA_Destroy(g_ident);
 
1327
  GA_Destroy(g_orbs);
 
1328
  //temporay global array for storage major transformation
 
1329
  GA_Destroy(g_tmp);
 
1330
 
 
1331
  return;
 
1332
}
 
1333
 
 
1334
void makoverlap(void) //not used due to USE_TRANSFORM 
 
1335
{
 
1336
      int me, lo[2], hi[2], ld[2];
 
1337
      int ld1, ld2;
 
1338
      double *ptr; //int--->double
 
1339
 
 
1340
      me = GA_Nodeid();
 
1341
      NGA_Distribution(g_ident, me, lo, hi);
 
1342
      NGA_Access(g_ident, lo, hi, &ptr, ld);
 
1343
      ld1 = hi[0] - lo[0] + 1;
 
1344
      ld2 = hi[1] - lo[1] + 1;
 
1345
      
 
1346
      setoverlap(ptr,lo,hi,ld1,ld2); 
 
1347
      NGA_Release(g_ident,lo,hi);
 
1348
      return;
 
1349
}
 
1350
 
 
1351
void setoverlap(double *a,int* lo,int* hi, int ld1,int ld2) //not used due to USE_TRANSFORM 
 
1352
{
 
1353
      int i,j,ii, jj;
 
1354
      
 
1355
      for (i = 0;i < ld1;i++) 
 
1356
      {
 
1357
        ii = i + lo[0]; 
 
1358
        for (j = 0;j < ld2;j++) 
 
1359
        {
 
1360
          jj = j + lo[1];
 
1361
#if USE_TRANSFORM
 
1362
          if (ii==jj)
 
1363
            //a[i][j] = 1.0;
 
1364
            a[ld2*i+j]=1.0;
 
1365
          else
 
1366
            //a[i][j] = 0.0;
 
1367
            a[ld2*i+j] = 0.0;
 
1368
#else
 
1369
          //a[i][j] = s[ii][jj];
 
1370
          a[ld1*i+j] = s[ii][jj];
 
1371
#endif
 
1372
        }
 
1373
      }
 
1374
      return;
 
1375
}
 
1376
 
 
1377
void print_GA_block(int g_a) //not used 
 
1378
{
 
1379
      int lo[2], hi[2], ld1, ld2, ld;
 
1380
      double *ptr; //int--->double
 
1381
 
 
1382
      //implicit declared variables in fortran
 
1383
      int me;
 
1384
 
 
1385
      me = GA_Nodeid();
 
1386
      NGA_Distribution(g_a, me, lo, hi);
 
1387
      ld1 = hi[0] - lo[0] + 1;
 
1388
      ld2 = hi[1] - lo[1] + 1;
 
1389
      NGA_Access(g_a, lo, hi, &ptr, &ld);
 
1390
      dump_chunk(ptr,ld1,ld2);
 
1391
      NGA_Release(g_a,lo,hi);
 
1392
      return;
 
1393
}
 
1394
 
 
1395
void print_GA_block_ij(int g_a,int tlo) //not used
 
1396
{
 
1397
      int lo[2], hi[2], ld1, ld2, ld;
 
1398
      double *ptr; //int--->double
 
1399
 
 
1400
      //implicit declared variables in fortran
 
1401
      int me;
 
1402
 
 
1403
      me = GA_Nodeid();
 
1404
      NGA_Distribution(g_a, me, lo, hi);
 
1405
      ld1 = hi[0] - lo[0] + 1;
 
1406
      ld2 = hi[1] - lo[1] + 1;
 
1407
      NGA_Access(g_a, &tlo, hi, &ptr, &ld);
 
1408
      dump_chunk(ptr,ld1,ld2);
 
1409
      NGA_Release(g_a,lo,hi);
 
1410
      return;
 
1411
}
 
1412
 
 
1413
void dump_chunk(double *a,int ld1,int ld2) //not used 
 
1414
{
 
1415
      //implicit declared variables in fortran
 
1416
      int i,j;
 
1417
      double trace;
 
1418
 
 
1419
      for (i = 0; i < MIN(10,ld1);i++) 
 
1420
      {
 
1421
         for (j = 0; j < MIN(10,ld2); j++)
 
1422
         {
 
1423
            //printf("%10.4f\n",a[i][j]);
 
1424
            printf("%10.4f\n",a[MIN(10,ld2)*i+j]);
 
1425
         }
 
1426
      }
 
1427
      trace = 0.00;
 
1428
      for (i=0;i<ld2;i++) 
 
1429
      {
 
1430
         //trace = trace +a[i][i];
 
1431
         trace = trace + a[MIN(10,ld2)*i+j];
 
1432
      }
 
1433
      printf("trace= %10.4f\n",trace);
 
1434
      return;
 
1435
}
 
1436
 
 
1437
double contract_matrices(int g_a, int g_b)
 
1438
{
 
1439
  int lo[2], hi[2], ptr_a, ptr_b, ld, ld1, ld2;
 
1440
  double a[ichunk][ichunk],b[ichunk][ichunk];
 
1441
  double value;
 
1442
  int dotask;
 
1443
 
 
1444
  //implicit declared variables in fortran
 
1445
  int i,j;
 
1446
 
 
1447
  double ret;
 
1448
 
 
1449
  //   evalute sum_ij a_ij*b_ij;
 
1450
  value = 0.00;
 
1451
  GA_Zero(g_counter);
 
1452
  dotask = next_chunk(lo, hi);
 
1453
  ld = ichunk;
 
1454
  while (dotask) {
 
1455
    NGA_Get(g_a, lo, hi, a, &ld);
 
1456
    NGA_Get(g_b, lo, hi, b, &ld);
 
1457
 
 
1458
    for (i = 0; i <= hi[0] - lo[0]; i++) {
 
1459
      for (j = 0; j <= hi[1] - lo[1]; j++) {
 
1460
        value += a[i][j] * b[i][j];   
 
1461
      }
 
1462
    }
 
1463
 
 
1464
    /*     for (j = 0; j <= hi[0] - lo[0]; j++) //column-major switch */
 
1465
    /*       { */
 
1466
    /*  for (i = 0;i <= hi[1]-lo[1];i++)  */
 
1467
    /*           { */
 
1468
    /*             value = value + a[j][i]*b[j][i];    */
 
1469
    /*           } */
 
1470
    /*       } */
 
1471
 
 
1472
    dotask = next_chunk(lo, hi);
 
1473
  }
 
1474
  GA_Dgop(&value, 1, "+");
 
1475
  ret = value;
 
1476
  return ret;
 
1477
}
 
1478
 
 
1479
#if 0
 
1480
double abs_double(double x)
 
1481
{
 
1482
   double ret;
 
1483
   if (x<0.0)
 
1484
      ret=-x;
 
1485
   else
 
1486
      ret=x;
 
1487
 
 
1488
   return ret;   
 
1489
  
 
1490
}
 
1491
#endif