Fixed UPF support

src/upf.c

@@ -56,7 +56,6 @@ int pspio_upf_read(FILE *fp, pspio_pspdata_t *pspdata)
     SUCCEED_OR_RETURN( upf_read_nlcc(fp, np, pspdata) );
   }
   SUCCEED_OR_RETURN( upf_read_nonlocal(fp, np, pspdata) );
-  fprintf(stderr, "TODO: fix upf_read_pswfc\n");
   SUCCEED_OR_RETURN( upf_read_pswfc(fp, np, pspdata) );
   SUCCEED_OR_RETURN( upf_read_local(fp, np, pspdata) );
   SUCCEED_OR_RETURN( upf_read_rhoatom(fp, np, pspdata) );

src/upf_read.c

@@ -191,32 +191,12 @@ int upf_read_mesh(FILE *fp, int np, pspio_pspdata_t *pspdata)
      feasibility. */
   DEFER_FUNC_ERROR( upf_tag_init(fp, "PP_R", NO_GO_BACK) );
   if ( pspio_error_get_last(__func__) == PSPIO_SUCCESS ) {
-    for (i=0; i<np; i+=4) {
-      FULFILL_OR_BREAK( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL, PSPIO_EIO );
-      nargs = sscanf(line, "%lf %lf %lf %lf", &tmp[0], &tmp[1], &tmp[2], &tmp[3]);
-      FULFILL_OR_BREAK( nargs < 5 && nargs > 0, PSPIO_EFILE_CORRUPT );
-      for (j=0; j<nargs; j++) r[i+j] = tmp[j];
-    }
+    SKIP_FUNC_ON_ERROR( read_array_4by4(fp, r, np) );
     SKIP_FUNC_ON_ERROR( upf_tag_check_end(fp, "PP_R") );
     
     /* Read Rab */
     SKIP_FUNC_ON_ERROR( upf_tag_init(fp, "PP_RAB", NO_GO_BACK) );
-    if ( (np % 4) == 0 ) {
-      nsup = np - 4;
-    } else {
-      nsup = np - (np % 4);
-    }
-    for (i=0; i<nsup; i+=4) {
-      FULFILL_OR_BREAK( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL, PSPIO_EIO );
-      nargs = sscanf(line, "%lf %lf %lf %lf", &tmp[0], &tmp[1], &tmp[2], &tmp[3]);
-      fprintf(stderr, ">>> PP_RAB ( %d / %d / %d ) [%d] <<<\n", i, nsup, np, nargs);
-      FULFILL_OR_BREAK( nargs == 4, PSPIO_EFILE_CORRUPT );
-      for (j=0; j<nargs; j++) drdi[i+j] = tmp[j];
-    }
-    DEFER_TEST_ERROR( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL, PSPIO_EIO );
-    nargs = sscanf(line, "%lf %lf %lf %lf", &tmp[0], &tmp[1], &tmp[2], &tmp[3]);
-    fprintf(stderr, ">>> PP_RAB END ( %d / %d ) [%d] <<<\n", nsup, np, nargs);
-    for (j=0; j<nargs; j++) drdi[nsup+j] = tmp[j];
+    SKIP_FUNC_ON_ERROR( read_array_4by4(fp, drdi, np) );
     SKIP_FUNC_ON_ERROR( upf_tag_check_end(fp, "PP_RAB") );
 
     /* Store the mesh in the pspdata structure */
@@ -252,12 +232,7 @@ int upf_read_nlcc(FILE *fp, int np, pspio_pspdata_t *pspdata)
   FULFILL_OR_EXIT(rho != NULL, PSPIO_ENOMEM);
 
   /* Read core rho */
-  for (i=0; i<np; i+=4) {
-    FULFILL_OR_BREAK( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL, PSPIO_EIO );
-    nargs = sscanf(line,"%lf %lf %lf %lf",&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
-    FULFILL_OR_BREAK( nargs < 5 && nargs > 0, PSPIO_EFILE_CORRUPT );
-    for (j=0; j<nargs; j++) rho[i+j] = tmp[j];
-  }
+  SKIP_FUNC_ON_ERROR( read_array_4by4(fp, rho, np) );
 
   /* Store the non-linear core corrections in the pspdata structure */
   SKIP_FUNC_ON_ERROR( pspio_xc_set_nlcc_density(pspdata->xc, pspdata->mesh, rho, NULL, NULL) );
@@ -353,14 +328,8 @@ int upf_read_nonlocal(FILE *fp, int np, pspio_pspdata_t *pspdata)
     FULFILL_OR_BREAK( sscanf(line, "%d", &proj_np) == 1, PSPIO_EFILE_CORRUPT );
 
     /* Read the projector function */
-    for (j=0; j<proj_np; j+=4){
-      FULFILL_OR_BREAK( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL, PSPIO_EIO );
-      nargs = sscanf(line,"%lf %lf %lf %lf", &tmp[0], &tmp[1], &tmp[2],
-        &tmp[3]);
-      FULFILL_OR_BREAK( nargs < 5 && nargs > 0, PSPIO_EFILE_CORRUPT );
+    SKIP_FUNC_ON_ERROR( read_array_4by4(fp, projector_read, np) );
 
-      for (k=0; k<nargs; k++) projector_read[k+j] = tmp[k];
-    }
     BREAK_ON_DEFERRED_ERROR;
 
     /* Fill with zeros, if any left */
@@ -432,12 +401,7 @@ int upf_read_local(FILE *fp, int np, pspio_pspdata_t *pspdata)
   }
 
   /* Read local potential */
-  for (i=0; i<np; i+=4) {
-    FULFILL_OR_BREAK( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL,  PSPIO_EIO );
-    nargs = sscanf(line,"%lf %lf %lf %lf",&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
-    FULFILL_OR_BREAK( nargs < 5 && nargs > 0, PSPIO_EFILE_CORRUPT );
-    for (j=0; j<nargs; j++) vlocal[i+j] = tmp[j];
-  }
+  SKIP_FUNC_ON_ERROR( read_array_4by4(fp, vlocal, np) );
 
   /* Convert from Rydbergs */
   for (i=0; i<np; i++) vlocal[i] /= 2.0;
@@ -502,19 +466,22 @@ int upf_read_pswfc(FILE *fp, int np, pspio_pspdata_t *pspdata)
   for (is=0; is<pspdata->n_states; is++) {
     /* Read the quantum numbers and occupations */
     FULFILL_OR_BREAK( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL, PSPIO_EIO );
-    FULFILL_OR_BREAK( sscanf(line, "%1d%1c %d %lf", &n, &ll, &l, &occ) == 4, PSPIO_EFILE_CORRUPT );
+    // TODO: ask QE developers whether they still accept the old format
+    //       like in the He pseudo
+    //FULFILL_OR_BREAK( sscanf(line, "%1d%1c %d %lf", &n, &ll, &l, &occ) == 4, PSPIO_EFILE_CORRUPT );
+    nargs = sscanf(line, "%1d%1c %d %lf", &n, &ll, &l, &occ);
+    if ( nargs != 4 ) {
+        n = 1;
+        nargs = sscanf(line, "%1c %d %lf", &ll, &l, &occ);
+        FULFILL_OR_BREAK( nargs == 3, PSPIO_EFILE_CORRUPT );
+    }
     FULFILL_OR_BREAK( sprintf(label, "%1d%1c", n, ll) > 0, PSPIO_EIO );
     SUCCEED_OR_BREAK( pspio_qn_init(qn, n, l, j[is]) );
     lmax = max(l, lmax);
 
     /* Read wavefunction */
-    for (ir=0; ir<np; ir+=4) {
-      FULFILL_OR_BREAK( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL, PSPIO_EIO );
-      nargs = sscanf(line, "%lf %lf %lf %lf",
-        &tmp[0], &tmp[1], &tmp[2], &tmp[3]);
-      FULFILL_OR_BREAK( nargs < 5 && nargs > 0, PSPIO_EFILE_CORRUPT );
-      for (i=0; i<nargs; i++) wf[ir+i] = tmp[i];
-    }
+    SKIP_FUNC_ON_ERROR( read_array_4by4(fp, wf, np) );
+
     BREAK_ON_DEFERRED_ERROR;
 
     /* UPF uses the wavefunctions multiplied by r */
@@ -556,12 +523,7 @@ int upf_read_rhoatom(FILE *fp, int np, pspio_pspdata_t *pspdata)
   FULFILL_OR_EXIT(rho_read != NULL, PSPIO_ENOMEM);
 
   /* Read valence density */
-  for (i=0; i<np; i+=4){
-    FULFILL_OR_BREAK( fgets(line, PSPIO_STRLEN_LINE, fp) != NULL, PSPIO_EIO );
-    nargs = sscanf(line,"%lf %lf %lf %lf", &tmp[0], &tmp[1],&tmp[2], &tmp[3]);
-    FULFILL_OR_BREAK( nargs < 5 && nargs > 0, PSPIO_EFILE_CORRUPT );
-    for (j=0; j<nargs; j++) rho_read[i+j] = tmp[j];
-  }
+  SKIP_FUNC_ON_ERROR( read_array_4by4(fp, rho_read, np) );
 
   /* UPF uses the density multiplied by 4*Pi*r*r */
   for (i=0; i<np; i++) {

src/upf_tag.c

@@ -84,8 +84,6 @@ int upf_tag_check_end(FILE *fp, const char *tag)
     status = PSPIO_EFILE_CORRUPT;
   }
 
-  fprintf(stderr, "line: %s", read_string);
-  fprintf(stderr, "upf_tag_check_end status: %d\n\n", status);
   return status;
 }