← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/tools/ga-5-2/tascel/examples/scf/twoelcpp.cc

  • 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/tascel/examples/scf/twoelcpp.cc
 
1
#include "ga.h"
 
2
#include "macdecls.h"
 
3
#include <stdio.h>
 
4
#include <math.h>
 
5
#include <unistd.h>
 
6
#include <assert.h>
 
7
 
 
8
// #define min(x,y) (((x)<=(y))?(x):(y))
 
9
 
 
10
#define maxatom  286
 
11
#define maxnbfn  (15*maxatom)
 
12
#define mxiter   30
 
13
#define maxnnbfn ((maxnbfn)*(maxnbfn+1)/2)
 
14
#define pi       3.141592653589793d0)
 
15
#define tol      (0.5e-3)
 
16
#define tol2e    (1.0e-6)
 
17
 
 
18
#define ichunk 20
 
19
 
 
20
#ifdef __cplusplus
 
21
extern "C" {
 
22
#endif
 
23
  extern int armci_me;
 
24
  void nga_get_(Integer *g_a,
 
25
                Integer *lo,
 
26
                Integer *hi,
 
27
                void    *buf,
 
28
                Integer *ld);
 
29
 
 
30
  void g_(double *gg, Integer *i, Integer *j, Integer *k, Integer *l);
 
31
  double contract_matrices_(Integer *g_a, Integer *g_b);
 
32
  void new_nga_acc_(Integer *g_a,
 
33
                    Integer *lo,
 
34
                    Integer *hi,
 
35
                    void    *buf,
 
36
                    Integer *ld,
 
37
                    void    *alpha);
 
38
 
 
39
  void nga_acc_(Integer *g_a,
 
40
                Integer *lo,
 
41
                Integer *hi,
 
42
                void    *buf,
 
43
                Integer *ld,
 
44
                void    *alpha);
 
45
 
 
46
  void ctwoel_(Integer *g_schwarz, Integer *g_dens, Integer *g_fock,
 
47
               double *schwmax, double *etwo,
 
48
               long long *pnbfn, long long *icut1, long long *icut2,
 
49
               long long *icut3, long long *icut4);
 
50
 
 
51
#ifdef __cplusplus
 
52
}
 
53
#endif
 
54
 
 
55
 
 
56
static void cpptwoel(Integer g_schwarz, Integer g_dens, Integer g_fock,
 
57
                     double schwmax, Integer nbfn,
 
58
                     long long *icut1, long long *icut2,
 
59
                     long long *icut3, long long *icut4);
 
60
 
 
61
void ctwoel_(Integer *g_schwarz, Integer *g_dens, Integer *g_fock,
 
62
             double *schwmax, double *etwo,
 
63
             long long *pnbfn, long long *icut1, long long *icut2,
 
64
             long long *icut3, long long *icut4) {
 
65
  const Integer nbfn = (Integer)*pnbfn;
 
66
  long long ijcnt,klcnt,ijklcnt;
 
67
 
 
68
  ijcnt = *icut1;
 
69
  klcnt = *icut2;
 
70
  ijklcnt = *icut3;
 
71
 
 
72
  cpptwoel(*g_schwarz, *g_dens, *g_fock, *schwmax, nbfn, 
 
73
           icut1, icut2, icut3, icut4);
 
74
 
 
75
  *etwo = 0.5*contract_matrices_(g_fock,g_dens);
 
76
 
 
77
  ijcnt = *icut1 - ijcnt;
 
78
  klcnt = *icut2 - klcnt;
 
79
  ijklcnt = (*icut3) - ijklcnt;
 
80
  *icut4 = *icut3;
 
81
  if (*icut3 <= 0L) {
 
82
    /*     no integrals may be calculated if there is no work for  */
 
83
    /*     this node (ichunk too big), or, something is wrong */
 
84
    printf("no two-electron integrals computed by node %d\n", GA_Nodeid());
 
85
    fflush(stdout);
 
86
  }
 
87
  return;
 
88
}
 
89
 
 
90
/*************************************************************************
 
91
 *
 
92
 * the code involving task pool stuff follows
 
93
 *
 
94
 *************************************************************************/
 
95
 
 
96
//#include <UniformTaskCollectionSplit.h>
 
97
#include <UniformTaskCollSplitData.h>
 
98
#include <DenseArray.h>
 
99
#include <algorithm>
 
100
 
 
101
using namespace std;
 
102
using namespace tascel;
 
103
 
 
104
static int next_4chunk(int g_counter, Integer nbfn, 
 
105
                       Integer lo[4], Integer hi[4], 
 
106
                       Integer *ilo, Integer *jlo,
 
107
                       Integer *klo, Integer *llo) {
 
108
  const long one=1;
 
109
  const int zero=0;
 
110
  long imax, itask;
 
111
  Integer itmp;
 
112
  long long total_ntsks;
 
113
  imax = nbfn/ichunk;
 
114
  if (nbfn - ichunk*imax > 0) {
 
115
    imax = imax + 1;
 
116
  }
 
117
  total_ntsks = imax*imax*imax*imax;
 
118
 
 
119
  itask = NGA_Read_inc(g_counter, (int *)&zero, one);
 
120
/*   printf("%d: itask=%d total_ntsks=%lf\n", armci_me, (int)itask, (double)total_ntsks); */
 
121
  if (itask < 0) {
 
122
    printf("next_4chunk: itask negative: %d imax:%d nbfn:%d ichunk:%d\n",
 
123
           (int) itask, (int)imax, (int)nbfn, (int)ichunk);
 
124
    printf("probable GA integer precision problem if %lf > 2^31\n", 
 
125
           pow(imax*1.0, 4));
 
126
    fflush(stdout);
 
127
    assert(0);
 
128
  }
 
129
  if (itask < total_ntsks) {
 
130
    itmp = itask;
 
131
    *ilo = itmp % imax;
 
132
    itmp /= imax;
 
133
    *jlo = itmp %imax;
 
134
    itmp /= imax;
 
135
    *klo = itmp % imax;
 
136
    itmp /= imax;
 
137
    assert(itmp < imax);
 
138
    *llo = itmp;
 
139
 
 
140
    lo[0] = *ilo*ichunk;
 
141
    lo[1] = *jlo*ichunk;
 
142
    lo[2] = *klo*ichunk;
 
143
    lo[3] = *llo*ichunk;
 
144
    hi[0] = min((*ilo+1)*ichunk,nbfn) - 1;
 
145
    hi[1] = min((*jlo+1)*ichunk,nbfn) - 1;
 
146
    hi[2] = min((*klo+1)*ichunk,nbfn) - 1;
 
147
    hi[3] = min((*llo+1)*ichunk,nbfn) - 1;
 
148
 
 
149
    return 1;
 
150
  }
 
151
  else {
 
152
    return 0;
 
153
  }
 
154
}
 
155
 
 
156
static void clean_chunk(double *chunk) {
 
157
  int i;
 
158
  for(i=0; i<ichunk*ichunk; i++) {
 
159
    chunk[i] = 0.0;
 
160
  }
 
161
}
 
162
 
 
163
typedef struct {
 
164
  int id;
 
165
} task_dscr_t;
 
166
 
 
167
typedef struct {
 
168
  Integer nbfn;
 
169
  double schwmax;
 
170
  int cnt1, cnt2;
 
171
  long long *icut1, *icut2, *icut3, *icut4;
 
172
  Integer g_schwarz, g_dens, g_fock;
 
173
} task_plo_t;
 
174
 
 
175
static void compute_lo_hi(long itask, long nbfn, long lo[4], long hi[4], 
 
176
                         long *ilo, long *jlo, long *klo, long *llo) {
 
177
  long imax;
 
178
  Integer itmp;
 
179
  long long total_ntsks;
 
180
  imax = nbfn/ichunk;
 
181
  if (nbfn - ichunk*imax > 0) {
 
182
    imax = imax + 1;
 
183
  }
 
184
  total_ntsks = imax*imax*imax*imax;
 
185
 
 
186
  if (itask < 0) {
 
187
    printf("next_4chunk: itask negative: %d imax:%d nbfn:%d ichunk:%d\n",
 
188
           (int) itask, (int)imax, (int)nbfn, (int)ichunk);
 
189
    printf("probable GA integer precision problem if %lf > 2^31\n", 
 
190
           pow(imax*1.0, 4));
 
191
    fflush(stdout);
 
192
    assert(0);
 
193
  }
 
194
  assert(itask < total_ntsks);
 
195
 
 
196
  itmp = itask;
 
197
  *ilo = itmp % imax;
 
198
  itmp /= imax;
 
199
  *jlo = itmp %imax;
 
200
  itmp /= imax;
 
201
  *klo = itmp % imax;
 
202
  itmp /= imax;
 
203
  assert(itmp < imax);
 
204
  *llo = itmp;
 
205
  
 
206
  lo[0] = *ilo*ichunk;
 
207
  lo[1] = *jlo*ichunk;
 
208
  lo[2] = *klo*ichunk;
 
209
  lo[3] = *llo*ichunk;
 
210
  hi[0] = min((*ilo+1)*ichunk,nbfn) - 1;
 
211
  hi[1] = min((*jlo+1)*ichunk,nbfn) - 1;
 
212
  hi[2] = min((*klo+1)*ichunk,nbfn) - 1;
 
213
  hi[3] = min((*llo+1)*ichunk,nbfn) - 1;  
 
214
}
 
215
 
 
216
static int compute_owner(long itask, const vector<DataColl*> &colls, int arrid, long nbfn)
 
217
{
 
218
  const int size = 2*sizeof(int);
 
219
  char idx[size];
 
220
  long int ilo, jlo, klo, llo, lo[4], hi[4], tmp;
 
221
  int *iidx = (int *)idx;
 
222
 
 
223
  compute_lo_hi(itask, nbfn, lo, hi, &ilo, &jlo, &klo, &llo);
 
224
 
 
225
  switch(arrid) {
 
226
  case 0:
 
227
    iidx[0] = lo[1]/ichunk;
 
228
    iidx[1] = lo[0]/ichunk;
 
229
    break;
 
230
  case 1:
 
231
    iidx[0] = lo[3]/ichunk;
 
232
    iidx[1] = lo[2]/ichunk;
 
233
    break;
 
234
  case 2:
 
235
    iidx[0] = lo[3]/ichunk;
 
236
    iidx[1] = lo[2]/ichunk;
 
237
    break;
 
238
  case 3:
 
239
    iidx[0] = lo[3]/ichunk;
 
240
    iidx[1] = lo[1]/ichunk;
 
241
    break;
 
242
  case 4:
 
243
    iidx[0] = lo[1]/ichunk;
 
244
    iidx[1] = lo[0]/ichunk;
 
245
    break;
 
246
  case 5:
 
247
    iidx[0] = lo[2]/ichunk;
 
248
    iidx[1] = lo[0]/ichunk;
 
249
    break;
 
250
  default:
 
251
    assert(0);
 
252
  } 
 
253
 
 
254
  return colls.at(arrid)->getProc(idx,size);
 
255
}
 
256
 
 
257
static void compute_index(void *dscr, int dscr_len, 
 
258
                         void *pldata, int pldata_len, int arrid, void *idx, int idxlen) {
 
259
  assert(dscr_len == sizeof(task_dscr_t));
 
260
  assert(arrid >=0 && arrid<6);
 
261
  assert(idxlen ==  2*sizeof(int));
 
262
  const long itask = ((task_dscr_t*)dscr)->id;
 
263
  int *iidx = (int *)idx;
 
264
  task_plo_t *ptplo = (task_plo_t*)pldata;
 
265
  assert(pldata_len == sizeof(task_plo_t));
 
266
  const long int nbfn = ptplo->nbfn;
 
267
  long int ilo, jlo, klo, llo, lo[4], hi[4], tmp;
 
268
 
 
269
  compute_lo_hi(itask, nbfn, lo, hi, &ilo, &jlo, &klo, &llo);
 
270
 
 
271
  switch(arrid) {
 
272
  case 0:
 
273
    iidx[0] = lo[1]/ichunk;
 
274
    iidx[1] = lo[0]/ichunk;
 
275
    break;
 
276
  case 1:
 
277
    iidx[0] = lo[3]/ichunk;
 
278
    iidx[1] = lo[2]/ichunk;
 
279
    break;
 
280
  case 2:
 
281
    iidx[0] = lo[3]/ichunk;
 
282
    iidx[1] = lo[2]/ichunk;
 
283
    break;
 
284
  case 3:
 
285
    iidx[0] = lo[3]/ichunk;
 
286
    iidx[1] = lo[1]/ichunk;
 
287
    break;
 
288
  case 4:
 
289
    iidx[0] = lo[1]/ichunk;
 
290
    iidx[1] = lo[0]/ichunk;
 
291
    break;
 
292
  case 5:
 
293
    iidx[0] = lo[2]/ichunk;
 
294
    iidx[1] = lo[0]/ichunk;
 
295
    break;
 
296
  default:
 
297
    assert(0);
 
298
  }  
 
299
}
 
300
 
 
301
 
 
302
static void twoel_task(UniformTaskCollection *utc, void *bigd, int bigd_len, 
 
303
                       void *pldata, int pldata_len, vector<void *> data_bufs);
 
304
 
 
305
static void cpptwoel(Integer g_schwarz, Integer g_dens, Integer g_fock,
 
306
                     double schwmax, Integer nbfn,
 
307
                     long long *icut1, long long *icut2,
 
308
                     long long *icut3, long long *icut4) {
 
309
  int dotask;
 
310
  int g_counter, ione = 1;
 
311
  task_dscr_t tdscr;
 
312
  task_plo_t tplo;
 
313
 
 
314
  {
 
315
    tplo.nbfn = nbfn;
 
316
    tplo.schwmax = schwmax;
 
317
    tplo.icut1 = icut1;
 
318
    tplo.icut2 = icut2;
 
319
    tplo.icut3 = icut3;
 
320
    tplo.icut4 = icut4;
 
321
    tplo.g_schwarz = g_schwarz;
 
322
    tplo.g_dens = g_dens;
 
323
    tplo.g_fock = g_fock;
 
324
 
 
325
    TslFuncRegTbl frt;
 
326
    TslFunc tf = frt.add(twoel_task);
 
327
    vector<DataColl*> colls;
 
328
    int block[2] = {ichunk, ichunk};
 
329
    DenseArray coll_schwarz(g_schwarz, block, 2);
 
330
    DenseArray coll_dens(g_dens, block, 2);
 
331
    DenseArray coll_fock(g_fock, block, 2);
 
332
    /**FIXME: modes are associated with arrays and not specific data
 
333
       movement operation*/
 
334
    colls.push_back(&coll_schwarz);
 
335
    colls.push_back(&coll_schwarz);
 
336
    colls.push_back(&coll_dens);
 
337
    colls.push_back(&coll_dens);
 
338
    colls.push_back(&coll_fock);
 
339
    colls.push_back(&coll_fock);
 
340
 
 
341
    vector<AccessMode> modes(4, MODE_RONLY);
 
342
    vector<int> idxlens(6, 2*sizeof(int));
 
343
 
 
344
    modes.push_back(MODE_ACC);
 
345
    modes.push_back(MODE_ACC);
 
346
 
 
347
    const long me = GA_Nodeid();
 
348
    const long nproc = GA_Nnodes();
 
349
    const long imax = nbfn/ichunk + ((nbfn%ichunk) ? 1 : 0);
 
350
    const long total_ntasks = imax*imax*imax*imax;
 
351
    const long ntasks_per_proc = (long)ceil(1.0*total_ntasks/nproc);
 
352
    const long tasklo = ntasks_per_proc * me;
 
353
    const long taskhi = min(tasklo+ntasks_per_proc,total_ntasks)-1;
 
354
 
 
355
    TaskCollProps props;
 
356
    props.functions(tf,frt).taskSize(sizeof(task_dscr_t))
 
357
#define EVEN_DISTRIBUTION 0
 
358
#if EVEN_DISTRIBUTION
 
359
        .maxTasks(ntasks_per_proc)
 
360
#else
 
361
        .maxTasks(total_ntasks)
 
362
#endif
 
363
        .localData(&tplo,sizeof(tplo));
 
364
    UniformTaskCollSplitData utc(props, colls, modes, idxlens, compute_index);
 
365
 
 
366
#if EVEN_DISTRIBUTION
 
367
    for(long i=tasklo; i<=taskhi; i++)
 
368
#else
 
369
    for (long i=0; i<total_ntasks; ++i)
 
370
      if (me == compute_owner(i, colls, 0, nbfn))
 
371
#endif
 
372
      {
 
373
        tdscr.id = i;
 
374
        utc.addTask(&tdscr, sizeof(tdscr));
 
375
      }
 
376
    GA_Sync();
 
377
    utc.process();
 
378
  }
 
379
}
 
380
 
 
381
typedef struct {
 
382
  double buf[ichunk][ichunk];
 
383
} tbuf_t;
 
384
 
 
385
static void twoel_task(UniformTaskCollection *utc, void *_bigd, int bigd_len, 
 
386
                       void *pldata, int pldata_len, vector<void *> data_bufs) {
 
387
  assert(_bigd!=NULL);
 
388
  assert(bigd_len == sizeof(task_dscr_t));
 
389
  assert(pldata!=NULL);
 
390
  assert(pldata_len == sizeof(task_plo_t));
 
391
  assert(data_bufs.size()==6);
 
392
 
 
393
  task_dscr_t *ptdscr = (task_dscr_t*)_bigd;
 
394
  task_plo_t *ptplo = (task_plo_t*)pldata;
 
395
 
 
396
  tbuf_t *s_ij = (tbuf_t*)data_bufs[0];
 
397
  tbuf_t *s_kl = (tbuf_t*)data_bufs[1];
 
398
  tbuf_t *d_kl = (tbuf_t*)data_bufs[2];
 
399
  tbuf_t *d_jl = (tbuf_t*)data_bufs[3];
 
400
  tbuf_t *f_ij = (tbuf_t*)data_bufs[4];
 
401
  tbuf_t *f_ik = (tbuf_t*)data_bufs[5];
 
402
 
 
403
  int ich;
 
404
  Integer lch;
 
405
  Integer lo_ik[2],hi_ik[2],lo_jl[2],hi_jl[2];
 
406
  Integer i,j,k,l,iloc,jloc,kloc,lloc,ld;
 
407
  double gg;
 
408
  Integer ione = 1;
 
409
  double one = 1.0;
 
410
  Integer nbfn;
 
411
  long lo[4], hi[4];
 
412
  long it, jt, kt, lt;
 
413
  long long *icut1, *icut2, *icut3, *icut4;
 
414
  Integer g_schwarz, g_dens, g_fock;
 
415
  double schwmax;
 
416
 
 
417
  long itask = ptdscr->id;
 
418
  nbfn = ptplo->nbfn;
 
419
  schwmax = ptplo->schwmax;
 
420
  icut1 = ptplo->icut1;
 
421
  icut2 = ptplo->icut2;
 
422
  icut3 = ptplo->icut3;
 
423
  icut4 = ptplo->icut4;
 
424
  g_schwarz = ptplo->g_schwarz;
 
425
  g_dens = ptplo->g_dens;
 
426
  g_fock = ptplo->g_fock;  
 
427
 
 
428
  ld = maxnbfn;
 
429
  ich = lch = ichunk;
 
430
 
 
431
  compute_lo_hi(itask, nbfn, lo, hi, &it, &jt, &kt, &lt);
 
432
 
 
433
  assert(ich == hi[0]-lo[0]+1);
 
434
  assert(ich == hi[2]-lo[2]+1);
 
435
  clean_chunk((double*)f_ij->buf);
 
436
  clean_chunk((double *)f_ik->buf);
 
437
  for(i = lo[0]; i<=hi[0]; i++) {
 
438
    iloc = i-lo[0];
 
439
    for(j = lo[1]; j<= hi[1]; j++) {
 
440
      jloc = j-lo[1];
 
441
      if (s_ij->buf[jloc][iloc]*(schwmax) < tol2e) {
 
442
        *icut1 = *icut1 + (hi[2]-lo[2]+1)*(hi[3]-lo[3]+1);
 
443
      }
 
444
      else {
 
445
        /* *cnt1+=1;*/
 
446
        for(k = lo[2]; k<=hi[2]; k++) {
 
447
          kloc = k-lo[2];
 
448
          for(l = lo[3]; l<=hi[3]; l++) {
 
449
            lloc = l-lo[3];
 
450
            if (s_ij->buf[jloc][iloc] * s_kl->buf[lloc][kloc] < tol2e) {
 
451
              *icut2 = *icut2 + 1;
 
452
            }
 
453
            else  {
 
454
              Integer _i=i+1, _j=j+1, _k=k+1, _l=l+1;
 
455
              g_(&gg, &_i, &_j, &_k, &_l);
 
456
              /* *cnt2 += 1;*/
 
457
              f_ij->buf[jloc][iloc] = f_ij->buf[jloc][iloc] + gg*d_kl->buf[lloc][kloc];
 
458
              f_ik->buf[kloc][iloc] = f_ik->buf[kloc][iloc] - 0.5*gg*d_jl->buf[lloc][jloc];
 
459
              *icut3 = *icut3 + 1;
 
460
            }
 
461
          }
 
462
        }
 
463
      }
 
464
    }
 
465
  }
 
466
}
 
467