diff --git a/src/elem_diag.irp.f b/src/elem_diag.irp.f
index 679b0d2..85b244f 100644
--- a/src/elem_diag.irp.f
+++ b/src/elem_diag.irp.f
@@ -28,6 +28,10 @@
 !       if(yw)write(6,*)iaa,'diag,v1'
 !         endif
        enddo
+       do i=1, natom
+            if(deter(i).ne.3) xmatd = xmatd + E(i)
+       enddo
+       xmatd = xmatd - E(natom+1)
 
 !-----stockage de l element diag
 
diff --git a/src/ex1.c b/src/ex1.c
index 0652da0..e062938 100644
--- a/src/ex1.c
+++ b/src/ex1.c
@@ -33,10 +33,7 @@ int main(int argc,char **argv)
   PetscInt       i,Istart,Iend,col[natomax],maxit,its,nconv,countcol;
   PetscInt		 nev, ncv, mpd;
   PetscLogDouble t1,t2,tt1,tt2;
-//PetscBool      FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
   PetscErrorCode ierr;
-//PetscScalar    eigr;
-//PetscScalar    eigi;
   int            mpiid;
 
   char const* const fileName = argv[1];
@@ -44,16 +41,14 @@ int main(int argc,char **argv)
   Data getdata;
   PetscInt		 nlocal;
   
+  /* gather the input data */
   Data_new(file, &getdata);
   getdata.n = get_ntot(getdata.FAM1, getdata.natom, getdata.isz, getdata.ntrou, getdata.fix_trou1, getdata.fix_trou2);
 
   nlocal = getdata.n/getdata.npar;
-//printf("n=%ld\t nsites=%d\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.natom, getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz);
-
 
   PetscScalar	 *valxr;
   PetscInt	   	 indxr[nlocal];
-//Vec            Vr,Vi;
   char           filename[PETSC_MAX_PATH_LEN]="FIL666";
   PetscViewer    viewer;
   PetscBool      ishermitian;
@@ -98,30 +93,18 @@ int main(int argc,char **argv)
   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
   ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,10*getdata.natom,NULL,10*getdata.natom,NULL,&A);CHKERRQ(ierr);
   ierr = MatMPIAIJSetPreallocation(A,10*getdata.natom,NULL,10*getdata.natom,NULL);CHKERRQ(ierr);
-//ierr = MatSetFromOptions(A);CHKERRQ(ierr);
-//ierr = MatSetUp(A);CHKERRQ(ierr);
 
   MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
   ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
   ierr = PetscTime(&tt1);CHKERRQ(ierr);
-  ierr = PetscPrintf(PETSC_COMM_WORLD," start: %d end: %d | %d\n",Istart, Iend,getdata.natom);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD," start: %d end: %d \n",Istart, Iend);CHKERRQ(ierr);
 
-//  Iend = range[mpiid];
-//  if(mpiid==0){
-//      Istart = 0;
-//  }
-//  else{
-//      Istart = range[mpiid-1];
-//  }
   for (i=Istart; i<Iend; i+=getdata.nnz) {
       tcountcol2=0;
       for(kk=0;kk<getdata.nnz;kk++){
           tcountcol[kk]=0;
       }
       iii=i+1;
-//    if(i%getdata.npar == 0 && mpiid==0){
-//      ierr = PetscTime(&t1);CHKERRQ(ierr);
-//    }
       unit_l1_(
             getdata.l1,
             getdata.l2,
@@ -131,6 +114,7 @@ int main(int argc,char **argv)
             getdata.xjjxy,
             getdata.xjjz ,
             getdata.xtt ,
+            getdata.E ,
             tcountcol,
             &getdata.ntrou,
             &getdata.isz,
@@ -147,7 +131,6 @@ int main(int argc,char **argv)
         for(kk=0;kk<tcountcol[ll]+1;kk++){
             value[kk] = val[kk+tcountcol2];
             col[kk] = tcol[kk+tcountcol2]-1;
-//          PetscPrintf(PETSC_COMM_WORLD,"value = %f col = %d\n",value[kk],col[kk]);
         }
         for(kk=tcountcol2+tcountcol[ll]+1;kk<natomax;kk++){
             value[kk] = 0.0;
@@ -169,8 +152,6 @@ int main(int argc,char **argv)
   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
   ierr = PetscTime(&tt2);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to assemble the matrix: %f\n",tt2-tt1);CHKERRQ(ierr);
-//ierr = MatGetVecs(A,NULL,&xr);CHKERRQ(ierr);
-//ierr = MatGetVecs(A,NULL,&xi);CHKERRQ(ierr);
   ierr = MatCreateVecs(A,NULL,&xr);CHKERRQ(ierr);
   ierr = MatCreateVecs(A,NULL,&xi);CHKERRQ(ierr);
 
@@ -208,9 +189,8 @@ int main(int argc,char **argv)
   /*
      Save eigenvectors, if  == ested
   */
-//PetscOptionsGetString(NULL,NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
   EPSGetConverged(eps,&nconv);
-  if (0) {
+  if (getdata.print_wf) {
     PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer);
   	PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
   	PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU);
@@ -255,10 +235,6 @@ int main(int argc,char **argv)
 	  norm4 = 0.0;
 
 
-//  	  for (ii=Istart; ii<Iend; ii+=1) {
-//		  indxr[ii-Istart] = ii;
-//	  }
-
 //        ierr = PetscTime(&tt1);CHKERRQ(ierr);
 //        ierr = VecGetArray(xr, &valxr);CHKERRQ(ierr);
           VecScatterCreateToAll(xr,&scatter,&vec2);
@@ -336,34 +312,12 @@ int main(int argc,char **argv)
             W3=weight3fin/normfin2;
 //          W3=weight3fin;
           }
-//        ierr = PetscTime(&tt2);CHKERRQ(ierr);
-//        ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to calc par S^2: %f\n",tt2-tt1);CHKERRQ(ierr);
 
-//PetscPrintf(PETSC_COMM_WORLD,"\n norm = %18f xymat =  %18f S^2 = %18f \n", norm4, xymat4, norm3);
   xymatfin = 0.0;
   normfin = 0.0;
 
-           /* sequential version of analyse */
-
-//           ierr = PetscTime(&tt1);CHKERRQ(ierr);
-//           VecScatterCreateToAll(xr,&scatter,&vec2);
-//           VecScatterBegin(scatter,xr,vec2,INSERT_VALUES,SCATTER_FORWARD);
-//           VecScatterEnd(scatter,xr,vec2,INSERT_VALUES,SCATTER_FORWARD);
-//           VecGetArray(vec2,&values);
-//           if(mpiid == 0){
-//             Istart = 0;
-//             Iend   = getdata.n;
-//             analyse_(values, (Iend-Istart), &Istart, &Iend, &xymat, &norm);
-//             XS=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin/normfin)));
-//             printf("\n norm = %18f xymat =  %18f S^2 = %18f \n", norm, xymat, XS);
-//           }
-//           VecRestoreArray(vec2,&values);
-//           ierr = PetscTime(&tt2);CHKERRQ(ierr);
-//           PetscPrintf(PETSC_COMM_WORLD,"seq time = %18f\n",tt2-tt1);
-
-
       /*
-         Compute the relative error associated to each eigenpair
+       * Compute the relative error associated to each eigenpair
       */
       ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRQ(ierr);
 
@@ -385,8 +339,6 @@ int main(int argc,char **argv)
     ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
   }
 
-//VecScatterDestroy(&scatter);
-//VecDestroy(&vec2); 
   ierr = EPSDestroy(&eps);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&xr);CHKERRQ(ierr);
diff --git a/src/get_ntot.c b/src/get_ntot.c
index 2ddea61..27b7e17 100644
--- a/src/get_ntot.c
+++ b/src/get_ntot.c
@@ -15,8 +15,8 @@ int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_t
       natom2 = natom;
   }
 
-  tnt1 =         (int)ceil(exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))));
-  printf("%10.5f | tnt1=%d\n",exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))),tnt1);
+  tnt1 =         (lrint)(exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))));
+//printf("%10.5f | tnt1=%d\n",exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))),tnt1);
   int            nalpha, nbeta;
 
     if((((natom-ntrou) + 2*isz) % 2) == 0){
@@ -34,7 +34,7 @@ int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_t
         }
     }
 
-  tnt2 =         (int)ceil(exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1)))));
-  printf("natom2=%d fix_trou1=%d fix_trou2=%d nalpha=%d nbeta=%d | | %d %d ntot=%d\n",natom2, fix_trou1, fix_trou2, nalpha, nbeta, tnt1, tnt2, tnt1*tnt2);
+  tnt2 =         (lrint)(exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1)))));
+//printf("natom2=%d fix_trou1=%d fix_trou2=%d nalpha=%d nbeta=%d | | %d %d ntot=%d\n",natom2, fix_trou1, fix_trou2, nalpha, nbeta, tnt1, tnt2, tnt1*tnt2);
   return tnt1*tnt2;
 }
diff --git a/src/get_ntot.h b/src/get_ntot.h
index 8cf2069..67bca0e 100644
--- a/src/get_ntot.h
+++ b/src/get_ntot.h
@@ -1,3 +1,4 @@
 #include <tgmath.h>
+#include <math.h>
 
 int get_ntot(_Bool Fam1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2);
diff --git a/src/provide_clust.irp.f b/src/provide_clust.irp.f
index 9ed9988..4f7e96e 100644
--- a/src/provide_clust.irp.f
+++ b/src/provide_clust.irp.f
@@ -4,6 +4,7 @@ BEGIN_PROVIDER[integer,l1,   (maxlien)]
 &BEGIN_PROVIDER[real*8, xtt  ,(maxlien)]
 &BEGIN_PROVIDER[real*8, xjjz ,(maxlien)]
 &BEGIN_PROVIDER[real*8, xjjxy,(maxlien)]
+&BEGIN_PROVIDER[real*8, E,(maxlien)]
 &BEGIN_PROVIDER[integer, ntrou]
 &BEGIN_PROVIDER[integer, isz]
 &BEGIN_PROVIDER[logical*1, FAM1]
diff --git a/src/read2.c b/src/read2.c
index b349799..8f8b040 100644
--- a/src/read2.c
+++ b/src/read2.c
@@ -11,9 +11,9 @@ void Data_new(FILE* file, Data* dat) {
 
 	while (fgets(line, sizeof(line), file)) {
 	
-	    /* note that fgets don't strip the terminating \n, checking its
-	       presence would allow to handle lines longer that sizeof(line) */
-	  if (count != 28){
+	    /* note that fgets doesn't strip the terminating \n, checking its
+	       presence would allow to handle lines longer than sizeof(line) */
+	  if (count != 30){
 			count++;
 			switch(count){
 				case 1:
@@ -178,54 +178,84 @@ void Data_new(FILE* file, Data* dat) {
 					  }
 					}
 					break;
-                case 13:
-					dat->nroots=atol(line);
+				case 13:
+					arrayIdx=0;
+					for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
+					{
+					  double val;
+					  char *unconverted;
+					  /**
+					   * Convert the next token to a float value
+					   */
+					  val = strtof(token, &unconverted);
+					  if (!isspace(*unconverted) && *unconverted != 0)
+					  {
+					    /**
+					     * Bad input string.  Again, we just bail.
+					     */
+					    fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
+					    break;
+					  }
+					  else
+					  {
+					    dat->E[arrayIdx++] = val;
+					  }
+					}
 					break;
                 case 14:
-					dat->s21a1=atol(line);
+					dat->nroots=atol(line);
 					break;
                 case 15:
-					dat->s21a2=atol(line);
+					dat->s21a1=atol(line);
 					break;
                 case 16:
-					dat->s21b1=atol(line);
+					dat->s21a2=atol(line);
 					break;
                 case 17:
-					dat->s21b2=atol(line);
+					dat->s21b1=atol(line);
 					break;
                 case 18:
-					dat->s22a1=atol(line);
+					dat->s21b2=atol(line);
 					break;
                 case 19:
-					dat->s22a2=atol(line);
+					dat->s22a1=atol(line);
 					break;
                 case 20:
-					dat->s22b1=atol(line);
+					dat->s22a2=atol(line);
 					break;
                 case 21:
-					dat->s22b2=atol(line);
+					dat->s22b1=atol(line);
 					break;
                 case 22:
-					dat->s23a1=atol(line);
+					dat->s22b2=atol(line);
 					break;
                 case 23:
-					dat->s23a2=atol(line);
+					dat->s23a1=atol(line);
 					break;
                 case 24:
-					dat->s23b1=atol(line);
+					dat->s23a2=atol(line);
 					break;
                 case 25:
-					dat->s23b2=atol(line);
+					dat->s23b1=atol(line);
 					break;
                 case 26:
-					dat->postrou=atol(line);
+					dat->s23b2=atol(line);
 					break;
                 case 27:
-					dat->fix_trou1=atol(line);
+					dat->postrou=atol(line);
 					break;
                 case 28:
+					dat->fix_trou1=atol(line);
+					break;
+                case 29:
 					dat->fix_trou2=atol(line);
 					break;
+				case 30:
+					dat->print_wf = atol(line);
+					break;
+                default:
+                    printf("Done reading file\n");
+                    break;
 		} /* end of switch */
 
 	  } /* end of the input file */
@@ -236,10 +266,13 @@ void Data_new(FILE* file, Data* dat) {
 }
 
 _Bool to_bool(const char* str) {
-    PetscBool strflg;
+    PetscBool strflg=PETSC_FALSE;
     PetscStrcmp("true\n",str, &strflg);
-    if(!strflg) PetscStrcmp("True\n",str, &strflg);
-    if(!strflg) PetscStrcmp("TRUE\n",str, &strflg);
+    if(strflg != PETSC_TRUE) PetscStrcmp("True\n",str, &strflg);
+    if(strflg != PETSC_TRUE) PetscStrcmp("TRUE\n",str, &strflg);
+    if(strflg != PETSC_TRUE) PetscStrcmp("true",str, &strflg);
+    if(strflg != PETSC_TRUE) PetscStrcmp("True",str, &strflg);
+    if(strflg != PETSC_TRUE) PetscStrcmp("TRUE",str, &strflg);
     return strflg;
 }
 
diff --git a/src/read2.h b/src/read2.h
index aa4aff3..9a1a6ff 100644
--- a/src/read2.h
+++ b/src/read2.h
@@ -12,13 +12,14 @@ typedef struct {
 	PetscInt n;
 	long int nnz,npar;
 	long int ntrou,isz;
-        _Bool FAM1;
-	long int l1[700];
-	long int l2[700];
-	long int ktyp[700];
-	double	xjjz[700];
-	double	xjjxy[700]; 
-	double	xtt[700]; 
+    _Bool FAM1;
+	long int l1[900];
+	long int l2[900];
+	long int ktyp[900];
+	double	xjjz[900];
+	double	xjjxy[900]; 
+	double	xtt[900]; 
+	double  E[900]; 
 	long int nroots;
 	int natom;
 	int s21a1;
@@ -36,6 +37,7 @@ typedef struct {
 	int postrou;
 	long int fix_trou1;
 	long int fix_trou2;
+    long int print_wf;
 
 } Data ;
 
diff --git a/src/stimsyr.h b/src/stimsyr.h
index 2397381..5964517 100644
--- a/src/stimsyr.h
+++ b/src/stimsyr.h
@@ -9,6 +9,7 @@ void unit_l1_(
             double *,
             double *,
             double *,
+            double *,
             long int *,
             long int *,
             long int *,
diff --git a/src/unit_l1.irp.f b/src/unit_l1.irp.f
index b9d010b..fc0dc40 100644
--- a/src/unit_l1.irp.f
+++ b/src/unit_l1.irp.f
@@ -6,6 +6,7 @@
         txjjxy,                &
         txjjz ,                &
         txtt  ,                &
+        tE    ,                &
         tcountcol,             &
         tntrou,                &
         tisz,                  &
@@ -24,7 +25,7 @@
         integer(C_SIZE_T),INTENT(INOUT)::tcol(maxlien)
         integer(C_SIZE_T),INTENT(INOUT),dimension(tnrows)::tcountcol
         integer(C_SIZE_T),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien)
-        real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien)
+        real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien), tE(maxlien)
    
 		nrows = tnrows
 		provide nrows
@@ -35,6 +36,7 @@
         xtt(i)     = txtt(i)
         xjjxy(i)   = txjjxy(i)
         xjjz (i)   = txjjz (i)
+        E (i)   = tE (i)
         enddo
         ntrou = tntrou
         isz = tisz