New implementation for evaluating dyson orbitals. Needs debugging.

Author: cthree-40

mpi_source/binarystr.c

@@ -324,19 +324,19 @@ void print_determinant(struct det d, int aelec, int belec)
  * print_occstring: print occupation string
  *  byte1, virtx, nvrtx
  */
-void print_occstring(struct occstr ostr, int nelec, int ndocc, int nactv)
+void print_occstring(struct occstr *ostr, int nelec, int ndocc, int nactv)
 {
         int *string = NULL;
         int i = 0, j = 0;
         int ecnt = 0;
         string = malloc(sizeof(int) * nelec);
         init_int_array_0(string, nelec);
-        nonzerobits(ostr.byte1, (ndocc + nactv), string);
+        nonzerobits(ostr->byte1, (ndocc + nactv), string);
         for (i = 0; i < nelec; i++) {
                 if (string[i] > 0) ecnt++;
         }
         for (i = ecnt; i < nelec; i++) {
-                string[i] = ostr.virtx[j];
+                string[i] = ostr->virtx[j];
                 j++;
         }
 

mpi_source/dysoncomp.c

@@ -95,89 +95,121 @@ int comparestrings_dyson(struct occstr str0, struct occstr str1, int ninto)
 }
 
 /*
- * compute_dyson_orbital: compute the dyson orbital between electronic
- * states of N+1 and N electron wavefuntions by comparing alpha/beta strings.
+ * compute_dyson_orbital_b: compute the dyson orbital between electronic
+ * states of N+1 and N electron wavefuntions by comparing beta strings.
  * Input:
  *
  */
-void compute_dyson_orbital(int v0_hndl, int v1_hndl, int w0_hndl, int w1_hndl,
-			   struct occstr *pstr0, struct eospace *peosp0, int npe0,
-			   struct occstr *qstr0, struct eospace *qeosp0, int nqe0,
-			   struct occstr *pstr1, struct eospace *peosp1, int npe1,
-			   struct occstr *qstr1, struct eospace *qeosp1, int nqe1,
-			   int norbs, int ndocc, int nactv, int ndyst0,
-			   int *dysnst0, int ndyst1, int *dysnst1, int ndets0,
-			   int ndets1, int sp1, int sp2)
+void compute_dyson_orbital_b(int v0_hndl, int v1_hndl, int w0_hndl, int w1_hndl,
+                             struct occstr *pstr0, struct eospace *peosp0, int npe0,
+                             struct occstr *qstr0, struct eospace *qeosp0, int nqe0,
+                             struct occstr *pstr1, struct eospace *peosp1, int npe1,
+                             struct occstr *qstr1, struct eospace *qeosp1, int nqe1,
+                             int **pq1, int npq1,
+                             int norbs, int ndocc, int nactv, int ndyst0,
+                             int *dysnst0, int ndyst1, int *dysnst1, int ndets0,
+                             int ndets1, int **strcont, int nbelec0,
+                             double **dyorb)
 {
 #define MAXBUFFER 10000
-    /* Local V and W buffers */
-    double **v0 = NULL, *v0data = NULL;
-    int v0cols = 0, v0rows = 0;
-    double **v1 = NULL, *v1data = NULL;
-    int v1cols = 0, v1rows = 0;
-    int **w0 = NULL, *w0data = NULL;
-    int **w1 = NULL, *w1data = NULL;
-    /* GA dimensions for buffers */
-    int v0_lo[2] = {0, 0}, v0_hi[2] = {0, 0}, v0_ld[1] = {0};
-    int v1_lo[2] = {0, 0}, v1_hi[2] = {0, 0}, v1_ld[1] = {0};
-    int w0_lo[2] = {0, 0}, w0_hi[2] = {0, 0}, w0_ld[1] = {0};
-    int w1_lo[2] = {0, 0}, w1_hi[2] = {0, 0}, w1_ld[1] = {0};
-
-    int buflen;
-    int i, imax, j, jmax;
-    int nvirt;
-
-    nvirt = norbs - ndocc - nactv;
+    /* (N+1) and (N) Vector arrays */
+    double **v0 = NULL,  *v01d = NULL; 
+    double **v1 = NULL,  *v11d = NULL; 
+    int v0_lo[2], v0_hi[2], v0_ld[1];
+    int v1_lo[2], v1_hi[2], v1_ld[1];
+    int v0_rows, v0_cols;
+    int v1_rows, v1_cols;
+    /* (N+1) and (N) wavefunction arrays */
+    int **w0 = NULL, *w01d = NULL;
+    int **w1 = NULL, *w11d = NULL;
+    int w0_lo[2], w0_hi[2], w0_ld[1];
+    int w1_lo[2], w1_hi[2], w1_ld[1];
+    /* Local variables */
+    int p0, q0; /* w0 values */
+    int p1, q1; /* w1 values */
+    int *j0indx = NULL;  /* Determinant index of |p0,q0> */
+    int *j1indx = NULL;  /* Determinant index of |p1,q1> */
+    int *o1indx = NULL;  /* Orbital difference of <p0,q0|p1,q1> */
+    int nj1 = 0;
+    int **vindx = NULL, *vindx1d = NULL;
+    int i, j, k;
+    int cnt;
 
-    /* Compute dimensions of V0. Allocate array. Get buffer. */
+    /* Get local distribution of V0. Allocate array and read in values. */
     NGA_Distribution(v0_hndl, mpi_proc_rank, v0_lo, v0_hi);
-    v0rows = v0_hi[1] - v0_lo[1] + 1;
-    v0cols = v0_hi[0] - v0_lo[0] + 1;
-    v0_ld[0]= v0rows;
-    v0data = allocate_mem_double_cont(&v0, v0rows, v0cols);
-    NGA_Get(v0_hndl, v0_lo, v0_hi, v0data, v0_ld);
-    /* Get W0 wavefunction information */
+    v0_rows = v0_hi[1] - v0_lo[1] + 1;
+    v0_cols = v0_hi[0] - v0_lo[0] + 1;
+    v0_ld[0]= v0_rows;
+    v01d = allocate_mem_double_cont(&v0, v0_rows, v0_cols);
+    NGA_Get(v0_hndl, v0_lo, v0_hi, v01d, v0_ld);
+    /* Get corresponding wavefunction information */
+    w01d = allocate_mem_int_cont(&w0, 3, v0_rows);
     w0_lo[0] = v0_lo[1];
     w0_lo[1] = 0;
     w0_hi[0] = v0_hi[1];
     w0_hi[1] = 2;
-    w0data = allocate_mem_int_cont(&w0, 3, v0rows);
-    NGA_Get(w0_hndl, w0_lo, w0_hi, w0data, w0_ld);
+    w0_ld[0] = 3;
+    NGA_Get(w0_hndl, w0_lo, w0_hi, w01d, w0_ld); 
 
-    buflen = MAXBUFFER; /* Set buffer length */
+    /* Allocate J index array and j-index array to grab V1 values from
+     * global array. */
+    nj1 = v0_rows * nbelec0;
+    j0indx = malloc(sizeof(int) * nj1);
+    j1indx = malloc(sizeof(int) * nj1);
+    o1indx = malloc(sizeof(int) * nj1);
+    vindx1d = allocate_mem_int_cont(&vindx, 2, (ndyst1 * nj1));
 
-    /* Get state number for v1 by getting number of columns
-     * on this process. (Array is distributed by row.) */
-    NGA_Distribution(v1_hndl, mpi_proc_rank, v1_lo, v1_hi);
-    v1cols = v1_hi[0] - v1_lo[0] + 1;
-    v1rows = buflen;
-    v1_ld[0] = v1rows;
-    v1_lo[0] = 0;
-    v1_hi[0] = 0;
-    v1data = allocate_mem_double_cont(&v1, v1rows, v1cols);
-    w1_lo[1] = 0;
-    w1_hi[1] = 2;
-    w1data = allocate_mem_int_cont(&w1, 3, v1rows);
-    
-    GA_Sync();
+    /* Allocate local buffer of V1 */
+    v11d = allocate_mem_double_cont(&v1, nj1, ndyst1);
 
-    for (i = 0; i < ndets1; i += buflen) {
-	imax = int_min((i + buflen - 1), (ndets1 - 1));
-	v1_lo[1] = i;
-	v1_hi[1] = imax;
-	v1rows = v1_hi[1] - v1_lo[1] + 1;
-	NGA_Get(v1_hndl, v1_lo, v1_hi, v1data, v1_ld);
-	w1_lo[0] = i;
-	w1_hi[0] = imax;
-	NGA_Get(w1_hndl, w1_lo, w1_hi, w1data, w1_ld);
-
+    /* Loop over w0 = (p,q), finding q'/p' in w1 */
+    cnt = 0;
+    for (i = 0; i < v0_rows; i++) {
+        p0 = w0[i][0];
+        q0 = w0[i][1];
+        /* Loop over q' */
+        for (j = 0; j < nbelec0; j++) {
+            p1 = p0;
+            q1 = strcont[q0][j]; /* New string */
+            j1indx[cnt] = string_info_to_determinant(p1, q1, peosp1, npe1,
+                                                     qeosp1, nqe1, pq1, npq1);
+            o1indx[cnt] = qstr0[q0].istr[j] - 1;
+            j0indx[cnt] = i;
+            cnt++;
+#ifdef DEBUGGING
+            //print_occstring(&(qstr0[q0]), nbelec0, ndocc, nactv);
+            //print_occstring(&(qstr1[q1]), (nbelec0 - 1), ndocc, nactv);
+            //printf(" | %d >, orbital = %d\n", j1indx[cnt - 1], o1indx[cnt - 1]);
+#endif
+        }
     }
 
-    GA_Sync();
-    deallocate_mem_cont_int(&w0, w0data);
-    deallocate_mem_cont_int(&w1, w1data);
-    deallocate_mem_cont(&v0, v0data);
-    deallocate_mem_cont(&v1, v1data);
+    /* Get indexes for all columns of J that we need */
+    set_ga_det_indexes_spec(j1indx, cnt, ndyst1, dysnst1, vindx);
+
+    /* Get V1 elements */
+    NGA_Gather(v1_hndl, v11d, vindx, (cnt * ndyst1));
+
+    /* Loop over contributions */
+    for (i = 0; i < cnt; i++) {
+
+        for (j = 0; j < ndyst0; j++) {
+            for (k = 0; k < ndyst1; k++) {
+                dyorb[j * ndyst1 + k][o1indx[i]] += v0[j][j0indx[i]] *
+                    v1[k][j1indx[i]];
+            }
+        }
+    }
+    
+               
+    /* Deallocate arrays */
+    deallocate_mem_cont(&v0, v01d);
+    deallocate_mem_cont(&v1, v11d);
+    deallocate_mem_cont_int(&w0, w01d);
+    deallocate_mem_cont_int(&vindx, vindx1d);
+    free(j0indx);
+    free(j1indx);
+    free(o1indx);
     return;
 }
 
@@ -299,10 +331,10 @@ void generate_strcontlist(struct occstr *str, int nstr, struct eospace *eosp0,
 	    newstr.virtx[1] = str[i].virtx[1];
 
 #ifdef DEBUGGING
-	    printf(" ORIG: ");
-	    print_occstring(str[i], nelec0, ndocc, nactv);
-	    printf(" NEW:  ");
-	    print_occstring(newstr, nelec1, ndocc, nactv);
+	    //printf(" ORIG: ");
+	    //print_occstring(&(str[i]), nelec0, ndocc, nactv);
+	    //printf(" NEW:  ");
+	    //print_occstring(&(newstr), nelec1, ndocc, nactv);
 #endif
 	    neospx = get_string_eospace(&newstr, ndocc, nactv, eosp1, ne1);
 	    strcont[i][cnt] = occstr2address(&newstr, eosp1[neospx], ndocc,
@@ -311,8 +343,111 @@ void generate_strcontlist(struct occstr *str, int nstr, struct eospace *eosp0,
 	}
 	/* Loop over external orbital electrons, removing them */
 	for (j = 0; j < str[i].nvrtx; j++) {
-	    
+            newstr.byte1 = str[i].byte1;
+            newstr.nvrtx = str[i].nvrtx - 1;
+            newstr.virtx[0] = str[i].virtx[0];
+            newstr.virtx[1] = str[i].virtx[1];
+            newstr.virtx[j] = 0;
+            if (j == 0 && str[i].nvrtx == 2) {
+                newstr.virtx[0] = newstr.virtx[1];
+                newstr.virtx[1] = 0;
+            }
+#ifdef DEBUGGING
+	    //printf(" ORIG: ");
+	    //print_occstring(&(str[i]), nelec0, ndocc, nactv);
+	    //printf(" NEW:  ");
+	    //print_occstring(&(newstr), nelec1, ndocc, nactv);
+#endif
+	    neospx = get_string_eospace(&newstr, ndocc, nactv, eosp1, ne1);
+	    strcont[i][cnt] = occstr2address(&newstr, eosp1[neospx], ndocc,
+					     nactv, nvirt, nelec1, elecx);
+	    cnt++;
+            
 	}
+        //printf("Count: %d\n", cnt);
+    }
+    return;
+}
+
+/*
+ * string_info_to_determinant: compute the determinant index given
+ * the p and q string information.
+ * Input:
+ *  pval    = p string
+ *  qval    = q string
+ *  peosp   = alpha electron orbital spaces
+ *  pegrps  = number of alpha electron orbital spaces
+ *  qeosp   = beta  electron orbital spaces
+ *  qegrps  = number of beta  electron orbital spaces
+ *  pq      = (p,q)-space pairings
+ *  npq     = number of (p,q)-space pairings
+ * Output:
+ *  detindx = determinant index in expansion.
+ */
+int string_info_to_determinant(int pval, int qval, struct eospace *peosp,
+                               int pegrps, struct eospace *qeosp, int qegrps,
+                               int **pq, int npq)
+{
+        int i, j, k;
+        int jstart, jmax, kstart, kmax;
+        int iloc, jloc;
+        int cntr = 0; /* Counter */
+        /* Find space combo for p,q */
+        for (i = 0; i < pegrps - 1; i++) {
+                if (pval >= peosp[i].start && pval < peosp[i + 1].start) break;
+        }
+        iloc = i;
+        for (j = 0; j < qegrps - 1; j++) {
+                if (qval >= qeosp[j].start && qval < qeosp[j + 1].start) break;
+        }
+        jloc = j;
+        for (i = 0; i < npq; i++) {
+                if (pq[i][0] == iloc && pq[i][1] == jloc) break;
+        }
+        iloc = i;
+        /* p,q pairing is found: i */
+        for (i = 0; i < iloc; i++) {
+                cntr = cntr + peosp[(pq[i][0])].nstr * qeosp[(pq[i][1])].nstr;
+        }
+        /* Search for index */
+        jstart = peosp[(pq[iloc][0])].start;
+        kstart = qeosp[(pq[iloc][1])].start;
+        jmax   = jstart + peosp[(pq[iloc][0])].nstr;
+        kmax   = kstart + qeosp[(pq[iloc][1])].nstr;
+        for (j = jstart; j < jmax; j++) {
+                for (k = kstart; k < kmax; k++) {
+                        if (j == pval && k == qval) {
+                                k = kmax;
+                                j = jmax;
+                                break;
+                        }
+                        cntr++;
+                }
+        }
+        return cntr;
+}
+
+/*
+ * set_ga_det_indexes_spec: set the array of indices to gather from global array.
+ * Input:
+ *  jindx   = list of row numbers in vector V
+ *  num     = number of row numbers
+ *  cols    = number of columns of V
+ *  colindx = column indices
+ * Output:
+ *  vindx   = indices of global array V to gather
+ */
+void set_ga_det_indexes_spec(int *jindx, int num, int cols, int *colindx,
+                             int **vindx)
+{
+    int cntr = 0; /* Index counter */
+    /* Arrays in global arrays are stored like: [column][row] */
+    for (int i = 0; i < cols; i++) {
+        for (int j = 0; j < num; j++) {
+            vindx[cntr][0] = colindx[i];
+            vindx[cntr][1] = jindx[j];
+            cntr++;
+        }
     }
     return;
 }

mpi_source/include/binarystr.h

@@ -84,7 +84,7 @@ void print_determinant(struct det d, int aelec, int belec);
  * print_occstring: print occupation string
  *  byte1, virtx, nvrtx
  */
-void print_occstring(struct occstr ostr, int nelec, int ndocc, int nactv);
+void print_occstring(struct occstr *ostr, int nelec, int ndocc, int nactv);
 
 /*
  * init_detlist: initialize determinant list