diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..5b2aa47
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+*bin/*
+*/obj/*
+*mod.mod*
+irpf90.a
diff --git a/Makefile b/Makefile
index fba9d00..f34b0b9 100644
--- a/Makefile
+++ b/Makefile
@@ -1,23 +1,56 @@
 include ${SLEPC_DIR}/lib/slepc/conf/slepc_common
-CC=gcc
+#CC=gcc
+#FC = ifort
+CLINKER = mpicc -fPIC  #-wd1572 -O3 -axAVX,SSE4.2 -fno-alias -no-prec-div -no-prec-sqrt -ip
 MAKE = /usr/bin/make
-MKDIR_P = mkdir -p
-OBJ_DIR = obj
+MKDIR_P = /bin/mkdir -p
+OBJ_DIR := obj
+LIB_DIR := libs
+BIN_DIR := bin
+SRC_DIR := src
 
-irpf90.a:
-	cd src && irpf90 init && $(MAKE) && cp irpf90.a ../
+.PHONY: ex1
+
+ex1: ${BIN_DIR}/ex1
 
 ${OBJ_DIR}:
 	${MKDIR_P} ${OBJ_DIR}
 
-directories: ${OBJ_DIR}
+${LIB_DIR}:
+	${MKDIR_P} ${LIB_DIR}
 
-obj/read2.o: src/read2.c directories
-	${CC} -c -o $@ $<
+${BIN_DIR}:
+	${MKDIR_P} ${BIN_DIR}
 
-obj/ex1.o: src/ex1.c
-	-${CC} -c -o $@ $< ${SLEPC_EPS_LIB}
+directories: ${OBJ_DIR} ${LIB_DIR} ${BIN_DIR}
 
-ex1: obj/read2.o irpf90.a obj/ex1.o src/read2.h src/stimsyr.h chkopts
-	    -${CLINKER} -o ex1 obj/ex1.o obj/read2.o irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3
+${LIB_DIR}/irpf90.a: directories
+	cd ${SRC_DIR} && irpf90 init && $(MAKE) irpf90.a && cp irpf90.a ../${LIB_DIR}
+
+${OBJ_DIR}/get_ntot.o: ${SRC_DIR}/get_ntot.c directories chkopts
+	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
+
+${OBJ_DIR}/read2.o: ${SRC_DIR}/read2.c directories chkopts
+	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
+
+${OBJ_DIR}/get_s2.o: ${SRC_DIR}/get_s2.c directories chkopts
+	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
+
+${OBJ_DIR}/get_s2_cyclic.o: ${SRC_DIR}/get_s2_cyclic.c directories chkopts
+	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
+
+${OBJ_DIR}/get_s2_mov.o: ${SRC_DIR}/get_s2_mov.c directories chkopts
+	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
+
+${OBJ_DIR}/get_dmat.o: ${SRC_DIR}/get_dmat.c directories chkopts
+	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
+
+${OBJ_DIR}/get_val_iaa2.o: ${SRC_DIR}/get_val_iaa2.c directories chkopts
+	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
+
+${OBJ_DIR}/ex1.o: ${SRC_DIR}/ex1.c
+	-${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB}
+
+${BIN_DIR}/ex1: ${OBJ_DIR}/get_ntot.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${OBJ_DIR}/ex1.o ${SRC_DIR}/read2.h ${SRC_DIR}/get_ntot.h ${SRC_DIR}/stimsyr.h chkopts
+	    -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_ntot.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB}
 #    ${RM} ex1.o read2.o
diff --git a/README.md b/README.md
index 186955e..8e979d6 100644
--- a/README.md
+++ b/README.md
@@ -7,57 +7,63 @@ Double Exchange Hamiltonian: Complete Version
 
 (under GNU GENERAL PUBLIC LICENSE v2)
 
-1. Dependencies
+_Dependencies_
 ---------------
 
-  1. [PETSc](https://www.mcs.anl.gov/petsc/documentation/installation.html) and [SLEPc](http://slepc.upv.es/documentation/current/docs/instal.htm)
+  1. [PETSc](https://www.mcs.anl.gov/petsc/documentation/installation.html) and [SLEPc](http://slepc.upv.es/documentation/instal.htm)
 
   2. [IRPF90](https://github.com/scemama/irpf90)
 
-2. Compiling
+_Compiling_
 ------------
 
   1. Export environment variables for PETSc and SLEPc
 
-	```shell
-	export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION}
-	export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION}
-	```
+```shell
+export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION}
+export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION}
+export C_INCLUDE_PATH+=:$PETSC_DIR/include/:$SLEPC_DIR/include:$PETSC_DIR/arch-linux2-c-debug/include/:$SLEPC_DIR/arch-linux2-c-debug/include
+# The "arch-linux2-c-debug" directory can have different names depending on PETSC and SLEPC installation procedure.
+```
+
 
   2. Make the executable
 
-	```shell
-	make ex1
-	```
+```shell
+make ex1
+```
 
-3. Using DEHam
+_Using DEHam_
 ---------------
 
   1. The DEHam program requires an input file which 
    has the topology of the Hamiltonian and the various parameters
    as explained below in a sample inputfile:
 
-	```python
-	140					# The total number of determinants
-	7					# The largest number of non-zero elements per row
-	2					# The number of processors used in parallel
-	1					# The number of holes
-	0					# The isz (ms-1/2) value
-	1,2,3,1,2,3,4,5,6,7	# The topology of the system is specified here
-	2,3,4,8,7,6,5,6,7,8	# first and second line contain the two sites linked
-	1,1,1,2,2,2,2,3,3,3	# third line contains the type of link (1 for t, J 2 for K and 3 for none)
-	.1430,-0.20,0.0000	# The three types of links this line gives J, K
-	.1430,-0.20,0.0000	# 
-	-1.00,0.0,0.00		# This line gives t
-	```
+```python
+8     # The number of orbitals (total)
+140   # The largest number of non-zero elements per row (Multiple of Ndet)
+1     # The total number of processors used in parallel (Multiple of Ndet)
+1     # The number of holes
+0     # The isz (ms-1/2) value
+true     # Restrict the hole to the 1'st (i.e. half of natom) Family of states. *false* for no restrictions
+1,2,3,1,2,3,4,5,6,7	# The topology of the system is specified here
+2,3,4,8,7,6,5,6,7,8	# first and second line contain the two sites linked
+1,1,1,2,2,2,2,3,3,3	# third line contains the type of link (1 for t or J, 2 for K and 3 for none)
+.1430,-0.20,0.0000	# The three types of links this line gives J, K
+.1430,-0.20,0.0000	# 
+-1.00,0.0,0.00		# This line gives t
+1                   # Currently unused (Perhaps can be used for potential energy per site in the future.)
+1                   # The total number of roots
+```
 
   2. running DEHam
 
-	```shell
-	mpiexec -n [nprocs] ./ex1 inpfile 
-	```
+```shell
+mpiexec -n [nprocs] ./ex1 inpfile 
+```
 
-4. Publications using this code
+_Publications using this code_
 -------------------------------
 
   1. High-Spin Chains and Crowns from Double-Exchange Mechanism [doi:10.3390/cryst6040039](http://www.dx.doi.org/10.3390/cryst6040039)
diff --git a/examples/small.inp b/examples/small.inp
new file mode 100644
index 0000000..8e1a91e
--- /dev/null
+++ b/examples/small.inp
@@ -0,0 +1,14 @@
+8
+140
+1
+1
+0
+true
+1,2,3,1,2,3,4,5,6,7
+2,3,4,8,7,6,5,6,7,8
+1,1,1,2,2,2,2,3,3,3
+.1430,-0.20,0.0000
+.1430,-0.20,0.0000
+-1.00,0.0,0.00
+1
+2
diff --git a/examples/three.inp b/examples/three.inp
new file mode 100644
index 0000000..8706cf2
--- /dev/null
+++ b/examples/three.inp
@@ -0,0 +1,14 @@
+3
+1
+1
+1
+1
+true
+1,2,1,2,3,4,5
+2,3,6,5,4,5,6
+1,1,2,2,2,3,3
+.1430,-0.20,0.0000
+.1430,-0.20,0.0000
+-1.00,0.0,0.00
+1
+1
diff --git a/src/.gitignore b/src/.gitignore
new file mode 100644
index 0000000..7ac9fbf
--- /dev/null
+++ b/src/.gitignore
@@ -0,0 +1,5 @@
+IRPF90_temp/
+IRPF90_man/
+irpf90.make
+irpf90_entities
+tags
\ No newline at end of file
diff --git a/src/adr.irp.f b/src/adr.irp.f
index 76db98d..2f79c9d 100644
--- a/src/adr.irp.f
+++ b/src/adr.irp.f
@@ -1,4 +1,5 @@
 subroutine adr(ideter,add)
+    use iso_c_binding
     implicit none
     BEGIN_DOC
     ! this subroutine provides the address of a detrminant 
@@ -7,25 +8,25 @@ subroutine adr(ideter,add)
     ! matches the given determinant.
     END_DOC
     integer,INTENT(INOUT)::ideter(natomax)
-    integer(kind=selected_int_kind(16)),INTENT(INOUT)::add
-    integer(kind=selected_int_kind(16))::det,deth,addh,detnew
+    integer(C_SIZE_T),INTENT(INOUT)::add
+    integer(C_SIZE_T)::deti,dethi,addh,detnew
     integer::count,i,j
 
-    det=0
+    deti=0
     detnew=0
-    deth=0
+    dethi=0
     count=0
-    call conv(ideter,det,deth)
+    call conv(ideter,deti,dethi)
     Do i=0,natom-1
-        if(BTEST(deth,i))then
+        if(BTEST(dethi,i))then
             count=count+1
         endif
-        if(BTEST(det,i))then
+        if(BTEST(deti,i))then
             detnew=IBSET(detnew,i-count)
         endif
     enddo
-    det=detnew
-    call searchdet(det,add,deth,addh)
+    deti=detnew
+    call searchdet(deti,add,dethi,addh)
         add = add + (nt1-addh)*(nt2)
 
 
diff --git a/src/adrfull.irp.f b/src/adrfull_old
similarity index 100%
rename from src/adrfull.irp.f
rename to src/adrfull_old
diff --git a/src/adrnew.irp.f b/src/adrnew.irp.f
new file mode 100644
index 0000000..229cbc7
--- /dev/null
+++ b/src/adrnew.irp.f
@@ -0,0 +1,58 @@
+subroutine adrfull()
+    implicit none
+    use iso_c_binding
+    BEGIN_DOC
+    ! this subroutine provides the address of a detrminant 
+    ! given in old format.
+    ! It searches in a list of generated determinants and
+    ! matches the given determinant.
+    END_DOC
+    integer,dimension(natomax)::ideter
+    integer(C_SIZE_T)::add
+    integer(C_SIZE_T)::deti,dethi,addh,detnew
+    integer::count,i,j
+
+    deti=0
+    detnew=0
+    dethi=0
+    do j=1,detfound
+        detnew=0
+        count=0
+        ideter=foundet(:,j)
+        call conv(ideter,deti,dethi)
+        Do i=0,natom-1
+            if(BTEST(dethi,i))then
+                count=count+1
+            endif
+            if(BTEST(deti,i))then
+                detnew=IBSET(detnew,i-count)
+            endif
+        enddo
+        deti=detnew
+        foundadd(j,1)=deti
+        foundadd(j,3)=j
+        foundaddh(j,1)=dethi
+        foundaddh(j,3)=j
+        call searchdet(deti,add,dethi,addh)
+!   enddo
+!   call sort()
+!   call searchdetfull()
+!   call desort()
+!   do i=1,detfound
+!       add = foundadd(i,2)
+!       addh = foundaddh(i,2)
+        foundadd(j,2) = add
+        foundaddh(j,2)= addh
+        add = add + (nt1-addh)*(nt2)
+        foundetadr(j)=add
+    enddo
+
+
+10  FORMAT(B64,I8,F8.2)
+15  FORMAT(B64,I8,I8,I8)
+11  FORMAT(B64,I3,B64)
+12  FORMAT(I5,$)
+13  FORMAT(B64,B64)
+14  FORMAT(B64,I14)
+16  FORMAT(B64,I14,I14)
+end
diff --git a/src/analyse.irp.f b/src/analyse.irp.f
new file mode 100644
index 0000000..730a41d
--- /dev/null
+++ b/src/analyse.irp.f
@@ -0,0 +1,833 @@
+      SUBROUTINE ANALYSE(vect, dimvect, startvect, endvect, xymat2, norm2)
+!     INCLUDE "nbtots.prm"
+      use iso_c_binding
+      IMPLICIT NONE
+      INTEGER   dimvect, nbtots, startvect, endvect
+      REAL*8,dimension(dimvect)::vect
+      INTEGER (C_SIZE_T)::add,kvect
+      INTEGER (C_SIZE_T)::iaa2,i
+      INTEGER ,dimension(natomax)::ideter
+      INTEGER ,dimension(natomax)::ideter2
+      REAL*8,allocatable       ::xz(:)
+      REAL*8::xmat,xymat
+      REAL*8::xmat1,xymat1
+      REAL*8::xmat2
+      REAL*8,INTENT(INOUT)::xymat2
+      REAL*8::xmat3,xymat3
+      REAL*8::sym,nonsym,proj_trou
+      REAL*8,allocatable      ::xalpha1(:)
+!     REAL*8,allocatable      ::vect(:)
+      REAL*8,INTENT(INOUT)::norm2
+      REAL*8::norm,norm1,norm3,proj_2trou
+      REAL*8::t1,t2,XS,XS1,XS2,XS3
+      REAL*8::resta_mono,resta_one,resta_bi,delta
+      INTEGER ::kko,kok,kkio,j,eigen,nigen,count
+      INTEGER ::cntrou,countlvect,ndim,iaa
+      INTEGER ::ipt_1,ipt_2,iptemp_1,iptemp_2
+      INTEGER ::ipt_3,iptemp_3
+      INTEGER ::ibougetrou,jstart
+      REAL*8 ,allocatable::eigenvectors(:,:)
+      REAL*8 ,allocatable::eigenvalues(:)
+      REAL*8 ,allocatable::WORK(:)
+      REAL*8 ,allocatable::AP(:)
+      REAL*8 ,allocatable::densmat(:,:)
+      REAL*8 ,allocatable::densmat2(:,:)
+      REAL*8 proj_1,extradiag_dmat2,ionic,nonionic
+      REAL*8 proj_2,sum,conduction,prob,prob2
+      INTEGER INFO,nrow,ncol,mmo,mom,kk,k,omm,okk
+      CHARACTER*1 JOBZ,UPLO
+      INTEGER::RESTA=0
+  
+!     allocate(vect(nbtots))
+      allocate(xalpha1(natomax))
+      allocate(xz  (natom/2))
+
+!     OPEN (unit=59,file='FIL1',form='formatted',status='old')
+!     OPEN (unit=217,file='SBOX217',form='formatted',status='REPLACE')
+!     REWIND  59
+!     READ (59,*)
+  
+!     print *,' in analyse', startvect, endvect
+!     print *,'nalpha=',nalpha,'nbeta=',nbeta
+!     PRINT *,natom,ntrou,nbtots,nt1,nt2,isz
+!     ndim=3
+      ndim=(natom/2)*((natom/2)-1)/2
+      allocate(AP((ndim)*((ndim)+1)/2))
+      allocate(WORK(3*(ndim)))
+      allocate(eigenvectors((ndim)*(ndim),1))
+      allocate(eigenvalues((ndim)))
+      allocate(densmat(ndim,ndim))
+      allocate(densmat2(ndim,ndim))
+!     Touch isz maxdet maxial maxlien maxplac nalpha natom natomax nbeta nbtots nt1 nt2 ntrou
+!     PRINT *,(vect(j),j=1,30)
+  
+	IF(RESTA .eq. 1)THEN
+	do i=1,natom/2
+	if(mod(natom/2,2).eq.0)then
+		xz(i)=(((natom/2)/2)-0.5d0)-(i-1.0d0)
+		write(6  ,*)i,xz(i)
+	else
+		xz(i)=((natom/2)-1.0d0)/2.0d0-(i-1.0d0)
+		write(6  ,*)i,xz(i)
+	endif
+	enddo
+        ENDIF
+
+!!    PROVIDE det deth
+      nigen=1
+      DO eigen=1,nigen
+!     READ (59,10) (vect(j),j=1,nbtots)
+      IF (ntrou.eq.1) THEN
+        norm=0.d0
+        norm1=0.d0
+        norm2=0.d0
+        norm3=0.d0
+        count=0
+        cntrou=0
+        iaa=0
+        iaa2=0
+        countlvect=0
+        proj_2trou=0.d0
+        resta_bi=0.d0
+        resta_mono=0.d0
+        resta_one=0.d0
+        xymat = 0.0d0
+        xymat1 = 0.0d0
+        xymat2 = 0.d0
+        xymat3 = 0.d0
+        proj_1=0d0
+        proj_2=0d0
+        densmat=0d0
+        densmat2=0d0
+        extradiag_dmat2=0d0
+        conduction=0d0
+        nrow=0
+        ncol=0
+
+            jstart=1
+            ipt_1=0
+            ipt_2=0
+            ipt_3=0
+            iptemp_1=0
+            iptemp_2=0
+            iptemp_3=0
+            nonionic=0.d0
+            ionic=0.d0
+            prob=0.d0
+            prob2=0.d0
+            sym=0.d0
+            nonsym=0.d0
+            ibougetrou=0
+
+            DO kvect=1, endvect-startvect
+              CALL getdet(kvect+startvect,ideter)
+
+!!----------------------------------------
+!!  RESTA
+!!----------------------------------------
+!!! mono
+!                 proj_trou=vect(kvect)**2
+!                 DO i=1,natom/2
+!                   IF (ideter(i).eq.3) THEN
+!                     delta=0.0d0
+!                   ELSE
+!                     delta=1.0d0
+!                   END IF
+!                   resta_mono=resta_mono+delta*xz(i)*xz(i)*proj_trou
+!                   resta_one=resta_one+delta*xz(i)*proj_trou
+!                 END DO
+!!! bi
+!                 DO i=1,natom/2
+!                   DO j=1,natom/2
+!                     IF (ideter(i).eq.3.or.ideter(j).eq.3.or.i.eq.j)           &
+!                      THEN
+!                       delta=0.0d0
+!                     ELSE
+!                       delta=1.0d0
+!                     END IF
+!                     resta_bi=resta_bi+delta*xz(i)*xz(j)*proj_trou
+!                   END DO
+!                 END DO
+!!----------------------------------------
+
+
+!!----------------------------------------
+!! Prob ionic non-ionic
+!!----------------------------------------
+!             ipt_1=0
+!             ipt_2=0
+!             ipt_3=0
+!             DO kko=1,3
+!               IF(ideter(kko).eq.3)THEN
+!                 ipt_1=ipt_1+1
+!               ENDIF
+!             ENDDO
+
+!             IF(ipt_1.eq.1)THEN
+!               DO kko=4,6
+!                 IF(ideter(kko).eq.3)THEN
+!                   ipt_2=ipt_2+1
+!                 ENDIF
+!               ENDDO
+!             ENDIF
+
+!             IF(ipt_2.eq.1)THEN
+!               DO kko=7,9
+!                 IF(ideter(kko).eq.3)THEN
+!                   ipt_3=ipt_3+1
+!                 ENDIF
+!               ENDDO
+!             ENDIF
+
+!             IF(ipt_3 .eq. 1)THEN
+!               nonionic=nonionic+vect(kvect)**2
+!             ELSE
+!               ionic=ionic+vect(kvect)**2
+!             ENDIF
+!!----------------------------------------
+!! S_box
+!!----------------------------------------
+
+              xmat=0.0d0
+              xmat1=0.0d0
+              xmat2=0.0d0
+!             IF (.TRUE.)THEN
+!!              IF (ideter(6).eq.3 ) THEN
+!!                norm=norm+vect(kvect)**2
+!!                DO kko=5,7
+!!                  DO kok=kko,7
+!!                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+!!                      xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2)
+!!                    ELSE
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+!!                        xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+!!                        xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+!!                        xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=2
+!!                        ideter2(kok)=1
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat=xmat+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+!!                        xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=1
+!!                        ideter2(kok)=2
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat=xmat+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                    END IF
+!!                  END DO
+!!                END DO
+!!  
+!!                DO kko=16,18
+!!                  DO kok=kko,18
+!!                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+!!                      xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2)
+!!                    ELSE
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+!!                        xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+!!                        xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+!!                        xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=2
+!!                        ideter2(kok)=1
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat=xmat+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+!!                        xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=1
+!!                        ideter2(kok)=2
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat=xmat+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                    END IF
+!!                  END DO
+!!                END DO
+!!  
+!!                DO kko=5,7
+!!                  DO kok=16,18
+!!                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+!!                      xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2)
+!!                    ELSE
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+!!                        xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+!!                        xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+!!                        xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=2
+!!                        ideter2(kok)=1
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat=xmat+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+!!                        xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=1
+!!                        ideter2(kok)=2
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat=xmat+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                    END IF
+!!                  END DO
+!!                END DO
+!!              END IF
+!!!----------------------------------------
+!!              xymat=xymat+xmat
+!!
+!!              IF (ideter(7).eq.3 ) THEN
+!!                norm1=norm1+vect(kvect)**2
+!!                DO kko=5,9
+!!                  DO kok=kko,9
+!!                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+!!                      xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2)
+!!                    ELSE
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+!!                        xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+!!                        xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+!!                        xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=2
+!!                        ideter2(kok)=1
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat1=xmat1+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+!!                        xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=1
+!!                        ideter2(kok)=2
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat1=xmat1+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                    END IF
+!!                  END DO
+!!                END DO
+!!  
+!!                DO kko=14,18
+!!                  DO kok=kko,18
+!!                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+!!                      xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2)
+!!                    ELSE
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+!!                        xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+!!                        xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+!!                        xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=2
+!!                        ideter2(kok)=1
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat1=xmat1+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+!!                        xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=1
+!!                        ideter2(kok)=2
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat1=xmat1+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                    END IF
+!!                  END DO
+!!                END DO
+!!
+!!                DO kko=5,9
+!!                  DO kok=14,18
+!!                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+!!                      xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2)
+!!                    ELSE
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+!!                        xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+!!                        xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+!!                        xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=2
+!!                        ideter2(kok)=1
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat1=xmat1+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+!!                        xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
+!!                        DO kkio=1,natom
+!!                          ideter2(kkio)=ideter(kkio)
+!!                        END DO
+!!                        ideter2(kko)=1
+!!                        ideter2(kok)=2
+!!                        CALL adr(ideter2, iaa2)
+!!                        xmat1=xmat1+vect(kvect)*vect(iaa2)
+!!                      END IF
+!!                    END IF
+!!                  END DO
+!!                END DO
+!!              END IF
+!!!----------------------------------------
+!!              xymat1=xymat1+xmat1
+
+              IF (.TRUE.)THEN
+                norm2=norm2+vect(kvect)**2
+                DO kko=1,natom/2
+                  DO kok=kko,natom/2
+                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+                      xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2)
+                    ELSE
+                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+                        xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
+                      END IF
+                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+                        xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
+                      END IF
+                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+                        xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
+                        DO kkio=1,natom
+                          ideter2(kkio)=ideter(kkio)
+                        END DO
+                        ideter2(kko)=2
+                        ideter2(kok)=1
+                        CALL adr(ideter2, iaa2)
+                        xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
+                      END IF
+                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+                        xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
+                        DO kkio=1,natom
+                          ideter2(kkio)=ideter(kkio)
+                        END DO
+                        ideter2(kko)=1
+                        ideter2(kok)=2
+                        CALL adr(ideter2, iaa2)
+                        xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
+                      END IF
+                    END IF
+                  END DO
+                END DO
+
+                DO kko=(natom/2)+1,natom
+                  DO kok=kko,natom
+                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+                      xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2)
+                    ELSE
+                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+                        xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
+                      END IF
+                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+                        xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
+                      END IF
+                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+                        xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
+                        DO kkio=1,natom
+                          ideter2(kkio)=ideter(kkio)
+                        END DO
+                        ideter2(kko)=2
+                        ideter2(kok)=1
+                        CALL adr(ideter2, iaa2)
+                        xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
+                      END IF
+                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+                        xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
+                        DO kkio=1,natom
+                          ideter2(kkio)=ideter(kkio)
+                        END DO
+                        ideter2(kko)=1
+                        ideter2(kok)=2
+                        CALL adr(ideter2, iaa2)
+                        xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
+                      END IF
+                    END IF
+                  END DO
+                END DO
+  
+                DO kko=1,natom/2
+                  DO kok=(natom/2)+1,natom
+                    IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
+                      xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2)
+                    ELSE
+                      IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
+                        xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
+                      END IF
+                      IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
+                        xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
+                      END IF
+                      IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
+                        xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
+                        DO kkio=1,natom
+                          ideter2(kkio)=ideter(kkio)
+                        END DO
+                        ideter2(kko)=2
+                        ideter2(kok)=1
+                        CALL adr(ideter2, iaa2)
+                        xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
+                      END IF
+                      IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
+                        xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
+                        DO kkio=1,natom
+                          ideter2(kkio)=ideter(kkio)
+                        END DO
+                        ideter2(kko)=1
+                        ideter2(kok)=2
+                        CALL adr(ideter2, iaa2)
+                        xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
+                      END IF
+                    END IF
+                  END DO
+                END DO
+              END IF
+!----------------------------------------
+!             print *,"norm = ",norm2,"xmat2 = ",xmat2,"vect =",vect(kvect),"natom=",natom
+              xymat2=xymat2+xmat2
+
+!           XS=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat/norm)))
+!           XS1=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat1/norm1)))
+            XS2=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat2/norm2)))
+            END DO
+!           print *,eigen,xymat2,XS2
+
+!----------------------------------------
+! Resta and probabilities of dets
+!----------------------------------------
+!            DO kvect=1,nbtots
+!              ideter2=ideter
+!                   CALL getdet(kvect,ideter)
+!
+!              IF (jstart.eq.1) THEN
+!                jstart=0
+!                ideter2=ideter
+!                DO kko=1,natom
+!                  IF (ideter(kko).eq.3) THEN
+!                    ipt_1=kko
+!                    DO kok=kko+1,natom
+!                      IF (ideter(kok).eq.3) THEN
+!                        ipt_2=kok
+!                          DO okk=kok+1,natom
+!                            IF (ideter(okk).eq.3) THEN
+!                              ipt_3=okk
+!                              EXIT
+!                            END IF
+!                          END DO
+!                        EXIT
+!                      END IF
+!                    END DO
+!                    EXIT
+!                  END IF
+!                END DO
+!              END IF
+!  
+!
+!              DO kko=1,natom
+!                IF (ideter(kko).eq.3) THEN
+!                  iptemp_1=kko
+!                  DO kok=kko+1,natom
+!                    IF (ideter(kok).eq.3) THEN
+!                      iptemp_2=kok
+!                        DO okk=kok+1,natom
+!                          IF (ideter(okk).eq.3) THEN
+!                            iptemp_3=okk
+!                            EXIT
+!                          END IF
+!                        END DO
+!                      EXIT
+!                    END IF
+!                  END DO
+!                  EXIT
+!                END IF
+!              END DO
+!  
+!              IF (iptemp_1.ne.ipt_1.or.iptemp_2.ne.ipt_2                        &
+!                    .or. iptemp_3.ne.ipt_3) THEN
+!                ibougetrou=1
+!                ipt_1=iptemp_1
+!                ipt_2=iptemp_2
+!                ipt_3=iptemp_3
+!              ELSE
+!                proj_trou=proj_trou+vect(kvect)**2
+!                ibougetrou=0
+!              END IF
+!
+!              IF (iptemp_1.eq.(9-iptemp_3+1)                                    &
+!                        .and. iptemp_2 .eq. 5 .and. iptemp_1.ne.4) THEN
+!                sym=sym+vect(kvect)**2
+!              ELSEIF (ideter(4).eq.3 .and. ideter(5).eq.3                       &
+!                        .and. ideter(6).eq.3) THEN
+!                nonsym=nonsym+vect(kvect)**2
+!              ELSE
+!!             ELSEIF(ideter(2).eq.3 .and. ideter(3).eq.3                        &
+!!                 .and. ideter(6).eq.3)then
+!                nonsym=nonsym+vect(kvect)**2
+!              END IF
+!
+!              IF (ideter(2).eq.3 .and. ideter(4).eq.3                           &
+!                        .and. ideter(8).eq.3) THEN
+!                prob=prob+vect(kvect)**2
+!              END IF
+!  
+!!             IF (ideter(1).eq.3 .and. ideter(2).eq.3                           &
+!!                       .and. ideter(3).eq.3) THEN
+!!               prob2=prob2+vect(kvect)**2
+!!             END IF
+!
+!              IF (ibougetrou.eq.1.or.kvect.eq.nbtots) THEN
+!!----------------------------------------
+!! mono
+!                DO i=1,natom/2
+!                  IF (ideter2(i).eq.3) THEN
+!                    delta=1.0d0
+!                  ELSE
+!                    delta=0.0d0
+!                  END IF
+!                  resta_mono=resta_mono+delta*xz(i)*xz(i)*proj_trou
+!                  resta_one=resta_one+delta*xz(i)*proj_trou
+!                END DO
+!! bi
+!                DO i=1,natom/2
+!                  DO j=1,natom/2
+!                    IF (ideter2(i).ne.3.or.ideter2(j).ne.3.or.i.eq.j)           &
+!                     THEN
+!                      delta=0.0d0
+!                    ELSE
+!                      delta=1.0d0
+!                    END IF
+!                    resta_bi=resta_bi+delta*xz(i)*xz(j)*proj_trou
+!                  END DO
+!                END DO
+!!----------------------------------------
+!                proj_trou=0.d0
+!                proj_trou=proj_trou+vect(kvect)**2
+!              END IF
+!            END DO
+
+!----------------------------------------
+! One particle density matrix
+!----------------------------------------
+!            DO kko=1,3
+!              DO kok=1,3
+!
+!                DO kvect=1,nbtots
+
+!                  CALL getdet(kvect,ideter)
+!                  ideter2=ideter
+!                  IF (ideter(kko).ne.3) THEN
+!                    IF (ideter(kok).eq.3) THEN
+!                      ideter2(kok)=ideter(kko)
+!                      ideter2(kko)=3
+!                      CALL adr(ideter2, iaa2)
+!                      densmat(kko,kok)=densmat(kko,kok)+vect(kvect)*        &
+!                       vect(iaa2)
+!                    END IF
+!                  END IF
+!                  IF (kko.eq.kok.and.ideter(kko).ne.3) THEN
+!                    densmat(kko,kko)=densmat(kko,kko)+vect(kvect)**2
+!                  END IF
+
+!                END DO
+
+!              END DO
+!            END DO
+
+!----------------------------------------
+! two particle density matrix
+!----------------------------------------
+
+!           DO kko=1,(natom/2)-1
+!             DO kok=kko+1,natom/2
+
+!               nrow=nrow+1
+!               ncol=0
+!               DO mmo=1,(natom/2)-1
+!                 DO mom=mmo+1,natom/2
+
+!                   ncol=ncol+1
+
+!                   DO kvect=1,nbtots
+!                     CALL getdet(kvect,ideter)
+!                     ideter2=ideter
+!                     IF (ideter(kko).eq.3.and.ideter(kok)                      &
+!                      .eq.3.and.ideter(mmo).ne.3.and.ideter(mom).ne.3)         &
+!                      THEN
+!                      if(ideter(kok).ne.3 .and. ideter(mom).ne.3)then
+!                       ideter2(kko)=ideter(mmo)
+!                       ideter2(mmo)=3
+!                       ideter2(kok)=ideter(mom)
+!                       ideter2(mom)=3
+!                       CALL adr(ideter2, iaa2)
+!                       densmat2(nrow,ncol)=densmat2(nrow,ncol)+                &
+!                        vect(kvect)*vect(iaa2)
+!                           print *,nrow,ncol,kko,kok,mmo,mom,
+!    *                     densmat2(nrow,ncol)
+!                      endif
+!                     END IF
+
+
+!                     IF (nrow.eq.ncol.and.ideter(mmo)                          &
+!                      .ne.3.and.ideter(mom).ne.3) THEN
+!                       densmat2(nrow,ncol)=densmat2(nrow,ncol)+                &
+!                        vect(kvect)**2
+
+!                     END IF
+
+!                   END DO
+
+
+!                 END DO
+!               END DO
+
+!             END DO
+!           END DO
+
+!----------------------------------------
+
+!----------------------------------------
+! conduction
+!----------------------------------------
+
+!            count=0
+!            DO kko=1,(natom/2)-2
+!              DO kok=kko+1,(natom/2)-1
+!                DO okk=kok+1,natom/2
+!  
+!                nrow=nrow+1
+!                ncol=0
+!                DO mmo=kko,kko+1
+!                  DO mom=kok,kok+1
+!                    DO omm=okk,okk+1
+!  
+!                    ncol=ncol+1
+!                    DO kvect=1,nbtots
+!                      CALL getdet(kvect,ideter)
+!                      ideter2=ideter
+!                      IF (abs(kko-mmo).eq.1.or.abs(kok-mom).eq.1                        &
+!                            .or. abs(okk-omm).eq.1) THEN
+!                        IF (mmo.le.natom/2.and.mom.le.natom/2 .and.                     &
+!                            omm.le.natom/2) THEN
+!                          IF (mmo.ne.mom .and. mom.ne.omm) THEN
+!                            IF (ideter(kko).eq.3 .and. ideter(kok).eq.3                 &
+!                                .and. ideter(okk).eq.3) THEN
+!                              ideter2(okk)=ideter2(omm)
+!                              ideter2(omm)=3
+!                              ideter2(kok)=ideter2(mom)
+!                              ideter2(mom)=3
+!                              ideter2(kko)=ideter2(mmo)
+!                              ideter2(mmo)=3
+!                              CALL adr(ideter2, iaa2)
+!!                             count=0
+!!                             do i=1,natom/2
+!!                                   if(ideter2(i).eq.3)then
+!!                                     count+=1
+!!                                   endif
+!!                             enddo
+!!                             print *,kko,kok,okk,mmo,mom,omm,iaa2
+!                              conduction=conduction+dabs(vect(kvect)*vect(iaa2))
+!                            END IF
+!                          END IF
+!                        END IF
+!                      END IF
+!                    END DO
+!  
+!                    END DO
+!                  END DO
+!                END DO
+!  
+!                END DO
+!              END DO
+!            END DO
+  
+!----------------------------------------
+
+!            DO j=1,ndim
+!               write(217,1022)j,(densmat(j,kko),kko=1,ndim)
+!            END DO
+!----------------------------------------
+! diagonalisation de mat
+! affiche vecteur
+!           JOBZ='V'
+! matrice sup
+!           UPLO='U'
+
+! matrice en vecteur ligne ...
+!           extradiag_dmat2=0d0
+!           k=0
+!           DO j=1,ndim
+!             DO i=1,j-1
+!               if(i.ne.j)then
+!                   extradiag_dmat2 = extradiag_dmat2 + dabs(densmat2(i,j))
+!               endif
+!             END DO
+!           END DO
+
+! appel subroutine LAPACK de diagonalisation :: double précision !!
+!           INFO=0
+!           CALL DSPEV (JOBZ, UPLO, ndim, AP, eigenvalues, eigenvectors,            &
+!             ndim, WORK, INFO)
+
+!           IF (INFO.ne.0) THEN
+!             PRINT *,'SUBROUTINE MATRIX: Error at dspev',info
+!             CALL exit (1)
+!           END IF
+
+!           proj_2=0.d0
+!           sum=0d0
+!           DO j=1,ndim
+!             proj_2=proj_2-eigenvalues(j)*log(eigenvalues(j))
+!             sum+=eigenvalues(j)
+!             write(214,*)eigenvalues(j)
+!           END DO
+
+!           XS=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat/norm)))
+!           XS2=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat2/norm2)))
+!           XS3=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat3/norm3)))
+!           WRITE (217,*) eigen,XS,norm
+      END IF
+
+      END DO
+
+
+
+   10 FORMAT (E25.0)
+ 1022 FORMAT(3x,I3,6(2x,F12.4))
+      END
diff --git a/src/conv.irp.f b/src/conv.irp.f
index f296c60..416f1d7 100644
--- a/src/conv.irp.f
+++ b/src/conv.irp.f
@@ -1,20 +1,21 @@
-subroutine conv(ideter,det,deth)
+subroutine conv(ideter,deti,dethi)
+    use iso_c_binding
     implicit none
     BEGIN_DOC
     ! this routine converts a detrminant in the old
     ! format into the new one and returns the determinant.
     END_DOC
     integer,INTENT(INOUT)::ideter(natomax)
-    integer(kind=selected_int_kind(16)),INTENT(INOUT)::det
-    integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth
+    integer(C_SIZE_T),INTENT(INOUT)::deti
+    integer(C_SIZE_T),INTENT(INOUT)::dethi
     integer::i
-    det=0
-    deth=0
+    deti=0
+    dethi=0
     do i=1,natom
         if(ideter(natom-i+1).eq.2 .and. ideter(natom-i+1).ne.3)then
-            det=IBSET(det,i-1)
+            deti=IBSET(deti,i-1)
         elseif(ideter(natom-i+1).eq.3)then
-            deth=IBSET(deth,i-1)
+            dethi=IBSET(dethi,i-1)
         endif
     enddo
 end
diff --git a/src/desort.irp.f b/src/desort.irp.f
index c4d09b1..9fb936b 100644
--- a/src/desort.irp.f
+++ b/src/desort.irp.f
@@ -1,14 +1,15 @@
 subroutine desort()
+    use iso_c_binding
     implicit none
     integer::i,j,ord,ordh
-    integer(kind=selected_int_kind(16))::add,addh,det,deth,addt
+    integer(C_SIZE_T)::add,addh,deti,dethi,addt
 
     do i=1,detfound-1
         do j=i+1,detfound
             if(foundaddh(i,3).gt.foundaddh(j,3))then
-                deth = foundaddh(i,1)
+                dethi = foundaddh(i,1)
                 foundaddh(i,1) = foundaddh(j,1)
-                foundaddh(j,1) = deth
+                foundaddh(j,1) = dethi
                 addh = foundaddh(i,2)
                 foundaddh(i,2) = foundaddh(j,2)
                 foundaddh(j,2) = addh
@@ -17,9 +18,9 @@ subroutine desort()
                 foundaddh(j,3) = ordh
             endif
             if(foundadd(i,3).gt.foundadd(j,3))then
-                det = foundadd(i,1)
+                deti = foundadd(i,1)
                 foundadd(i,1) = foundadd(j,1)
-                foundadd(j,1) = det
+                foundadd(j,1) = deti
                 add = foundadd(i,2)
                 foundadd(i,2) = foundadd(j,2)
                 foundadd(j,2) = add
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 b29a2b5..c2a2659 100644
--- a/src/ex1.c
+++ b/src/ex1.c
@@ -1,56 +1,103 @@
-#include <slepceps.h>
 #include <petsctime.h>
-#include "stimsyr.h"
+#include <petscvec.h>
+
 #include "read2.h"
+#include "stimsyr.h"
+#include "get_s2.h"
+#include "get_ntot.h"
 
 #undef __FUNCT__
 #define __FUNCT__ "main"
 
+void solvequad(double *a, double *b, double *c, double *res){
+    *res = -*b/(2.0*(*a)) + sqrt((*b)*(*b) - 4.0*(*a)*(*c))/(2.0*(*a));
+}
+
 int main(int argc,char **argv)
 {
   Mat            A;           /* problem matrix */
   EPS            eps;         /* eigenproblem solver context */
   EPSType        type;
   PetscReal      error,tol,re,im;
-  PetscScalar    kr,ki,value[700];
+  PetscReal      norm=0.0;
+  PetscReal      norm2=0.0;
+  PetscReal      norm3=0.0;
+  PetscReal      norm4=0.0;
+  PetscReal      normfin=0.0;
+  PetscReal      normfin2=0.0;
+  PetscReal      normfin3=0.0;
+  PetscReal      normfin4=0.0;
+  const int            natomax=900;
+  PetscScalar    kr,ki,value[natomax];
   Vec            xr,xi;
-  PetscInt       i,Istart,Iend,col[700],maxit,its,nconv,countcol;
+  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];
   FILE* file = fopen(fileName, "r");
   Data getdata;
+  PetscInt		 nlocal;
   
+  /* gather the input data */
   Data_new(file, &getdata);
-//printf("n=%ld\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz);
+  getdata.n = get_ntot(getdata.FAM1, getdata.natom, getdata.isz, getdata.ntrou, getdata.fix_trou1, getdata.fix_trou2);
 
+  nlocal = getdata.n/getdata.npar;
 
-//Vec            Vr,Vi;
+  PetscScalar	 *valxr;
+  PetscInt	   	 indxr[nlocal];
   char           filename[PETSC_MAX_PATH_LEN]="FIL666";
   PetscViewer    viewer;
   PetscBool      ishermitian;
-  PetscInt       kk,ll,iii2;
+  PetscInt       kk,ll,mm,nn,iii2,iiii;
+  PetscInt       ii;
   long int       iii;
-  long int       tcountcol2,tcol[700],tcountcol[getdata.nnz];
-  double         val[700];
+  long int       tcountcol2,tcol[natomax],tcountcol[getdata.nnz];
+  double         val[natomax];
+  PetscReal      xymat=0.0;
+  PetscReal      xymat2=0.0;
+  PetscReal      xymat3=0.0;
+  PetscReal      xymat4=0.0;
+  PetscReal      xymatfin=0.0;
+  PetscReal      xymatfin2=0.0;
+  PetscReal      xymatfin3=0.0;
+  PetscReal      xymatfin4=0.0;
+  PetscReal      weight3fin = 0.0;
+  PetscReal      XS = 0.0;
+  PetscReal      XS2 = 0.0;
+  PetscReal      XS3 = 0.0;
+  PetscReal      XS4 = 0.0;
+  PetscReal      W3  = 0.0;
+  PetscReal      weight3 = 0.0;
+  PetscReal      trace1rdm=0.0;
+  PetscReal      trace1rdmfin=0.0;
+  PetscReal      trace2rdm=0.0;
+  PetscReal      trace2rdmfin=0.0;
+  IS                from, to; /* index sets that define the scatter */
+  PetscInt          idx_to[nlocal], idx_from[nlocal];
+  PetscScalar       *values;
+  int 		   	 ndim=(getdata.natom/2)*((getdata.natom/2)-1)/2;
+  double            a, b, c;
+  double            gamma_p = 0.0, gamma_m = 0.0;
+  double            gamma_pfin = 0.0, gamma_mfin = 0.0;
+  double            nel, s2dens;
+  double            nelfin, s2densfin;
+  double            densmat2[getdata.natom][getdata.natom][getdata.natom][getdata.natom];
+  memset(densmat2, 0, sizeof(densmat2));
 
   SlepcInitialize(&argc,&argv,(char*)0,NULL);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr);
   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
-  ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,getdata.nnz*getdata.npar,NULL,getdata.nnz*getdata.npar,NULL,&A);CHKERRQ(ierr);
-  ierr = MatMPIAIJSetPreallocation(A,getdata.nnz*getdata.npar,NULL,getdata.nnz*getdata.npar,NULL);CHKERRQ(ierr);
-//ierr = MatSetFromOptions(A);CHKERRQ(ierr);
-//ierr = MatSetUp(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);
 
   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 \n",Istart, Iend);CHKERRQ(ierr);
 
   for (i=Istart; i<Iend; i+=getdata.nnz) {
       tcountcol2=0;
@@ -58,53 +105,46 @@ int main(int argc,char **argv)
           tcountcol[kk]=0;
       }
       iii=i+1;
-      if(i%getdata.npar == 0 && mpiid==0){
-        ierr = PetscTime(&t1);CHKERRQ(ierr);
-      }
       unit_l1_(
             getdata.l1,
             getdata.l2,
             getdata.ktyp,
             &iii,
-			&getdata.nnz,
+            &getdata.nnz,
             getdata.xjjxy,
             getdata.xjjz ,
             getdata.xtt ,
+            getdata.E ,
             tcountcol,
             &getdata.ntrou,
             &getdata.isz,
+            &getdata.fix_trou1,
+            &getdata.fix_trou2,
+            &getdata.FAM1,
             tcol,
             val);
-     if(i%getdata.npar == 0 && mpiid==0){
-       ierr = PetscTime(&t2);CHKERRQ(ierr);
-       ierr = PetscPrintf(PETSC_COMM_WORLD," i: %d\n mpiid: %d\ntime: %f\n",i,mpiid,t2-t1);CHKERRQ(ierr);
-     }
+//    if(i%getdata.npar == 0 && mpiid==0){
+//      ierr = PetscPrintf(PETSC_COMM_WORLD," i: %d \n",i);CHKERRQ(ierr);
+//    }
       for(ll=0;ll<getdata.nnz;ll++){
 
-      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<700;kk++){
-          value[kk] = 0.0;
-          col[kk] = 0;
-      }
-      tcountcol2=tcountcol2 + tcountcol[ll]+1;
-      countcol=tcountcol[ll]+1;
-      if(i%getdata.npar == 0 && mpiid==0){
-        ierr = PetscTime(&t1);CHKERRQ(ierr);
-      }
-      iii2=i+ll;
-    ierr = MatSetValues(A,1,&iii2,countcol,col,value,INSERT_VALUES);CHKERRQ(ierr);
-      if(i%getdata.npar == 0 && mpiid==0){
-        ierr = PetscTime(&t2);CHKERRQ(ierr);
-        ierr = PetscPrintf(PETSC_COMM_WORLD," processor \ntime: %f\n",t2-t1);CHKERRQ(ierr);
-      }
+        for(kk=0;kk<tcountcol[ll]+1;kk++){
+            value[kk] = val[kk+tcountcol2];
+            col[kk] = tcol[kk+tcountcol2]-1;
+        }
+        for(kk=tcountcol2+tcountcol[ll]+1;kk<natomax;kk++){
+            value[kk] = 0.0;
+            col[kk] = 0;
+        }
+        tcountcol2=tcountcol2 + tcountcol[ll]+1;
+        countcol=tcountcol[ll]+1;
+        iii2=i+ll;
+        ierr = MatSetValues(A,1,&iii2,countcol,col,value,INSERT_VALUES);CHKERRQ(ierr);
       }
   }
   ierr = PetscTime(&tt2);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to build the matrix: %f\n",tt2-tt1);CHKERRQ(ierr);
+  printf("time = %f mpiid = %d \n",tt2-tt1, mpiid);
 
 
   ierr = PetscTime(&tt1);CHKERRQ(ierr);
@@ -112,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);
 
@@ -121,15 +159,14 @@ int main(int argc,char **argv)
   ierr = EPSSetOperators(eps,A,NULL);CHKERRQ(ierr);
   ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
   ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr);
-  ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr);
 
   ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
-  tol = 1.e-8;
+  tol = 1.e-9;
   maxit = 10000000;
   ierr = EPSSetTolerances(eps,tol,maxit);CHKERRQ(ierr);
-  nev  = 4;
   ncv  = 10;
   mpd  = 10;
+  nev  = getdata.nroots;
   ierr = EPSSetDimensions(eps,nev,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
 
   ierr = PetscTime(&t1);CHKERRQ(ierr);
@@ -141,29 +178,145 @@ int main(int argc,char **argv)
   ierr = EPSGetType(eps,&type);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRQ(ierr);
   ierr = EPSGetDimensions(eps,&nev,NULL,NULL);CHKERRQ(ierr);
-  ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenvalues: %D\n",nev);CHKERRQ(ierr);
   ierr = EPSGetTolerances(eps,&tol,&maxit);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);CHKERRQ(ierr);
 
   ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr);
 //ierr = EPSPrintSolution(eps,NULL);CHKERRQ(ierr);
 
-  if (nconv>0) {
-    /*
-       Display eigenvalues and relative errors
-    */
-    ierr = PetscPrintf(PETSC_COMM_WORLD,
-         "           k          ||Ax-kx||/||kx||\n"
-         "   ----------------- ------------------\n");CHKERRQ(ierr);
+  /*
+     Save eigenvectors, if  == ested
+  */
+  EPSGetConverged(eps,&nconv);
+  if (getdata.print_wf) {
+    PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer);
+  	PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
+  	PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU);
+    EPSIsHermitian(eps,&ishermitian);
+    for (i=0;i<nev;i++) {
+      EPSGetEigenvector(eps,i,xr,xi);
+      VecView(xr,viewer);
+#if !defined(PETSC_USE_COMPLEX)
+      if (!ishermitian) { VecView(xi,viewer); }
+#endif
+    }
+    PetscViewerDestroy(&viewer);
+  }
+
+  /*
+   * now analyzing the eigenvector
+   */
+
+  if (nconv>0) {
+
+    ierr = PetscPrintf(PETSC_COMM_WORLD,
+         "           k          ||Ax-kx||/||kx||        <S^2>\n"
+         "   ----------------- ----------------- ------------------\n");CHKERRQ(ierr);
+
+   for(i=0;i<nev;i++){
 
-    for (i=0;i<nconv;i++) {
       /*
         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
         ki (imaginary part)
       */
+      Vec               vec2;
+      VecScatter        scatter; /* scatter context */
       ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr);
+	  xymat = 0.0;
+	  xymat2 = 0.0;
+	  xymat3 = 0.0;
+	  xymat4 = 0.0;
+	  weight3 = 0.0;
+	  norm = 0.0;
+	  norm2 = 0.0;
+	  norm3 = 0.0;
+	  norm4 = 0.0;
+
+
+//        ierr = PetscTime(&tt1);CHKERRQ(ierr);
+//        ierr = VecGetArray(xr, &valxr);CHKERRQ(ierr);
+          VecScatterCreateToAll(xr,&scatter,&vec2);
+          VecScatterBegin(scatter,xr,vec2,INSERT_VALUES,SCATTER_FORWARD);
+          VecScatterEnd(scatter,xr,vec2,INSERT_VALUES,SCATTER_FORWARD);
+          ierr = VecGetArray(vec2,&values);CHKERRQ(ierr);
+          get_s2(xr, &Istart, &Iend, values, &getdata.natom, &norm, &norm2, &norm3, &norm4, &xymat, &xymat2, &xymat3, &xymat4, &weight3,
+                  &getdata.s21a1, &getdata.s21a2, &getdata.s21b1, &getdata.s21b2, &getdata.s22a1, &getdata.s22a2,
+                  &getdata.s22b1, &getdata.s22b2,  &getdata.s23a1, &getdata.s23a2,
+                  &getdata.s23b1, &getdata.s23b2, &getdata.postrou, natomax);
+//        get_s2_cyclic(xr, &Istart, &Iend, values, &getdata.natom, &norm, &norm2, &norm3, &norm4, &xymat, &xymat2, &xymat3, &xymat4,
+//                &getdata.s21a1, &getdata.s21a2, &getdata.s21b1, &getdata.s21b2, &getdata.s22a1, &getdata.s22a2,
+//                &getdata.s22b1, &getdata.s22b2,  &getdata.s23a1, &getdata.s23a2,
+//                &getdata.s23b1, &getdata.s23b2, &getdata.postrou, natomax);
+//        get_1rdm(values, &Istart, &Iend, &getdata.natom, &trace1rdm, natomax);
+//        get_2rdm(values, &Istart, &Iend, &getdata.natom, &trace2rdm, densmat2, natomax);
+//        analyse_(valxr, (Iend-Istart), &Istart, &Iend, &xymat, &norm);
+          VecRestoreArray(vec2,&values);
+          ierr = VecRestoreArray(xr, &valxr);CHKERRQ(ierr);
+          MPI_Reduce(&xymat, &xymatfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+          MPI_Reduce(&xymat2, &xymatfin2, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+          MPI_Reduce(&xymat3, &xymatfin3, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+          MPI_Reduce(&xymat4, &xymatfin4, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+          MPI_Reduce(&weight3, &weight3fin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+          MPI_Reduce(&norm, &normfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+          MPI_Reduce(&norm2, &normfin2, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+          MPI_Reduce(&norm3, &normfin3, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+          MPI_Reduce(&norm4, &normfin4, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+//        MPI_Reduce(&trace1rdm, &trace1rdmfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+//          printf("done calc densmat\n");
+//          for(ll=0;ll<getdata.natom/2;ll++){
+//              for(kk=0;kk<getdata.natom/2;kk++){
+//                  gamma_p = gamma_p + 0.5*(densmat2[ll][kk][kk][ll] + densmat2[ll][kk][ll][kk]);
+//                  gamma_m = gamma_m + 0.5*(densmat2[ll][kk][kk][ll] - densmat2[ll][kk][ll][kk]);
+//              }
+//          }
+//          MPI_Reduce(&trace2rdm, &trace2rdmfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+//          MPI_Reduce(&gamma_p, &gamma_pfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+//          MPI_Reduce(&gamma_m, &gamma_mfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
+//          if(mpiid==0){
+//            for(kk=0;kk<getdata.natom;kk++){
+//                for(ll=0;ll<getdata.natom;ll++){
+//                    for(mm=0;mm<getdata.natom;mm++){
+//                        for(nn=0;nn<getdata.natom;nn++){
+////                          printf("%d\t%d\t%d\t%d\t%18f\n",kk,ll,mm,nn,densmat2[kk][ll][mm][nn]);
+//                        }
+//                    }
+//                }
+//            }
+//            /* calc nel */
+//            a=1.0;
+//            b=-1.0;
+//            c=-2.0*(gamma_mfin + gamma_pfin);
+//            printf("\n  gp= %18f gm= %18f a=%18f b=%18f c=%18f\n", gamma_pfin, gamma_mfin, a, b, c);
+//            nel = -b/(2.0*(a)) + sqrt((b)*(b) - 4.0*(a)*(c))/(2.0*(a));
+////          solvequad(&a, &b, &c, &nel);
+//
+//            /* calc s^2 */
+//            a=1.0;
+//            b=1.0;
+//            c=-1.0*((gamma_mfin - gamma_pfin) - nel*(nel - 4.0)/4.0);
+//            s2dens = -b/(2.0*(a)) + sqrt((b)*(b) - 4.0*(a)*(c))/(2.0*(a));
+////          solvequad(&a, &b, &c, &s2dens);
+//            printf("\n mpiid = %d  # trace = %18f nel = %18f s2dens = %18f\n", mpiid, trace2rdmfin, nel, s2dens);
+//          }
+
+          if(!mpiid){
+            XS=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin/normfin)));
+//          XS2=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin2/normfin2)));
+//          XS3=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin3/normfin3)));
+            XS2=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin2)));
+            XS3=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin3)));
+            XS4=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin4/normfin4)));
+            XS4=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin4/normfin4)));
+            W3=weight3fin/normfin2;
+//          W3=weight3fin;
+          }
+
+  xymatfin = 0.0;
+  normfin = 0.0;
+
       /*
-         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);
 
@@ -177,32 +330,14 @@ int main(int argc,char **argv)
       if (im!=0.0) {
         ierr = PetscPrintf(PETSC_COMM_WORLD," %14f%+14fi %12g\n",(double)re,(double)im,(double)error);CHKERRQ(ierr);
       } else {
-        ierr = PetscPrintf(PETSC_COMM_WORLD,"   %18f       %12g\n",(double)re,(double)error);CHKERRQ(ierr);
+        ierr = PetscPrintf(PETSC_COMM_WORLD,"   %18f     %12g %18f %18f %18f %18f\n",(double)re,(double)error,(double)XS,(double)XS2,(double)XS3, (double)W3);CHKERRQ(ierr);
       }
+      VecScatterDestroy(&scatter);
+      VecDestroy(&vec2); 
     }
     ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
   }
 
-  /*
-     Save eigenvectors, if requested
-  */
-//PetscOptionsGetString(NULL,NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
-  EPSGetConverged(eps,&nconv);
-  if (nconv>0) {
-    PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer);
-  	PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
-  	PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU);
-    EPSIsHermitian(eps,&ishermitian);
-    for (i=0;i<nconv;i++) {
-      EPSGetEigenvector(eps,i,xr,xi);
-      VecView(xr,viewer);
-#if !defined(PETSC_USE_COMPLEX)
-      if (!ishermitian) { VecView(xi,viewer); }
-#endif
-    }
-    PetscViewerDestroy(&viewer);
-  }
-
   ierr = EPSDestroy(&eps);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&xr);CHKERRQ(ierr);
diff --git a/src/extra_diag.irp.f b/src/extra_diag.irp.f
index 3ee3ef2..0395347 100644
--- a/src/extra_diag.irp.f
+++ b/src/extra_diag.irp.f
@@ -1,8 +1,9 @@
         subroutine extra_diag(tistart)
+        use iso_c_binding
         implicit none
 
-        integer(kind=selected_int_kind(16)) :: iaa,iaa2,tistart,tistart2
-        integer(kind=selected_int_kind(16)) :: imat4,jmat4
+        integer(C_SIZE_T) :: iaa,iaa2,tistart,tistart2
+        integer(C_SIZE_T) :: imat4,jmat4
         integer :: i,ik,iik,j
         integer :: ik1,ik2,IC,k,ikmax,ikmin,count,count2,detfound2
         integer,allocatable :: ideter2(:)
@@ -99,7 +100,7 @@
             tistart=tistart+1
 
     enddo
-      Touch foundet foundetadr detfound foundadd foundaddh foundetdmat
+      Touch foundet foundetadr detfound foundadd foundaddh foundetdmat det deth
       call adrfull()
 
       do i=1,detfound
diff --git a/src/foundet.irp.f b/src/foundet.irp.f
index a7c73c7..e95c062 100644
--- a/src/foundet.irp.f
+++ b/src/foundet.irp.f
@@ -1,8 +1,9 @@
+use iso_c_binding
 BEGIN_PROVIDER[integer, foundet,(natomax,maxlien)]
-&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundetadr,(maxlien)]
+&BEGIN_PROVIDER[integer(C_SIZE_T), foundetadr,(maxlien)]
 &BEGIN_PROVIDER[real, foundetdmat,(maxlien)]
-&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundadd,(maxlien,3)]
-&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundaddh,(maxlien,3)]
+&BEGIN_PROVIDER[integer(C_SIZE_T), foundadd,(maxlien,3)]
+&BEGIN_PROVIDER[integer(C_SIZE_T), foundaddh,(maxlien,3)]
 &BEGIN_PROVIDER[integer, detfound]
     BEGIN_DOC
     ! provides all found determinants
diff --git a/src/get_dmat.c b/src/get_dmat.c
new file mode 100644
index 0000000..4fc5db5
--- /dev/null
+++ b/src/get_dmat.c
@@ -0,0 +1,143 @@
+#include <stdio.h>
+#include <petsctime.h>
+#include <slepceps.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include "get_dmat.h"
+
+
+/*
+ *---------------------------------------
+ * One particle density matrix
+ *----------------------------------------
+ *
+ * The One particle density matrix
+ * Input
+ * =====
+ * valxr    = The full vector
+ * Istart   = Local starting id
+ * Iend     = Local ending id
+ * natom    = number of sites
+ * Output
+ * =====
+ * trace    = trace
+ */
+void get_1rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace1rdm, const int natomax){
+
+  long int			 ideter[natomax];
+  long int			 ideter2[natomax];
+  int 		   	 kko,kok,kkio;
+  long int       ii;
+  PetscInt      iiii;
+  long int       iii;
+  long int       iaa2, iaa;
+  int 		   	 ndim=(*natom)*(*natom)/8-(*natom)/2;
+  double            densmat[ndim][ndim];
+  memset(densmat, 0, sizeof(densmat[0][0]) * ndim * ndim);
+
+            for(kko=0;kko<(*natom/2);kko++){
+              for(kok=0;kok<(*natom/2);kok++){
+
+          	  for(ii=*Istart;ii<*Iend;ii++) {
+                      iii = ii + 1;
+                      iiii = ii;
+                      getdet_(&iii, ideter);
+                      for(kkio=0;kkio<=*natom-1;kkio++){
+                        ideter2[kkio]=ideter[kkio];
+                      }
+                      if(ideter[kko] != 3){
+                        if(ideter[kok] == 3){
+                          ideter2[kok]=ideter[kko];
+                          ideter2[kko]=3;
+                          adr_(ideter2, &iaa2);
+                          densmat[kko][kok]=densmat[kko][kok]+valxr[iiii]*valxr[iaa2];
+                        }
+                      }
+                      if(kko == kok && ideter[kko] != 3){
+                        densmat[kko][kko]=densmat[kko][kko]+valxr[iiii]*valxr[iiii];
+                      }
+
+                  }
+                  if(kko == kok){
+                      *trace1rdm+=densmat[kko][kko];
+                  }
+
+              }
+            }
+} /** END **/
+
+/*
+ * 
+ *----------------------------------------
+ * two particle density matrix
+ *----------------------------------------
+ * Input
+ * =====
+ * valxr    = The full vector
+ * Istart   = Local starting id
+ * Iend     = Local ending id
+ * Output
+ * =====
+ * trace    = trace
+ */
+void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, double densmat2[*natom][*natom][*natom][*natom], const int natomax){
+
+  long int			 ideter[natomax];
+  long int			 ideter2[natomax];
+  int 		   	 kko,kok,kkio;
+  int 		   	 mmo,mom,mmio;
+  long int       ii;
+  PetscInt      iiii;
+  long int       iii;
+  long int       iaa2, iaa;
+  long int       nrow=-1, ncol=-1;
+//int 		   	 ndim=(*natom/2)*((*natom/2)-1)/2;
+//double            densmat2[ndim][ndim];
+//memset(densmat2, 0, sizeof(densmat2[0][0]) * ndim * ndim);
+
+            for(kko=0;kko<(*natom/2);kko++){
+              for(kok=0;kok<(*natom/2);kok++){
+
+               nrow=nrow+1;
+               ncol=-1;
+               for(mmo=0;mmo<(*natom/2);mmo++){
+                 for(mom=0;mom<(*natom/2);mom++){
+
+                   ncol=ncol+1;
+
+          	   for(ii=*Istart;ii<*Iend;ii++) {
+                       iii = ii + 1;
+                       iiii = ii;
+                       getdet_(&iii, ideter);
+                       for(kkio=0;kkio<=*natom-1;kkio++){
+                         ideter2[kkio]=ideter[kkio];
+                       }
+                       if(ideter[kko] == 3 && ideter[kok] == 3 && kko != kok && mmo != mom){
+                          ideter2[kko]=ideter[mmo];
+                          ideter2[mmo]=3;
+                          ideter2[kok]=ideter[mom];
+                          ideter2[mom]=3;
+                          adr_(ideter2, &iaa2);
+                          densmat2[kko][kok][mmo][mom]=densmat2[kko][kok][mmo][mom]+valxr[iiii]*valxr[iaa2];
+                       }
+
+
+                     if(kko == mmo && kok == mom && ideter[kko]==3 && ideter[kok]==3 && kko != kok){
+                       densmat2[kko][kok][mmo][mom]=densmat2[kko][kok][mmo][mom]+valxr[iiii]*valxr[iiii];
+                     }
+
+                   }
+                          printf("%d\t%d\t%d\t%d\t%18f\n",kko,kok,mmo,mom,densmat2[kko][kok][mmo][mom]);
+
+
+                    if(kko == mmo && kok == mom)*trace2rdm+=densmat2[kko][kok][mmo][mom];
+                 }
+               }
+
+
+              }
+            }
+
+} /** END **/
+
diff --git a/src/get_dmat.h b/src/get_dmat.h
new file mode 100644
index 0000000..990a4c0
--- /dev/null
+++ b/src/get_dmat.h
@@ -0,0 +1,9 @@
+#include <stdio.h>
+#include <petsctime.h>
+#include <slepceps.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+
+void get_1rdm(PetscScalar *, PetscInt *, PetscInt *, int *,  PetscReal *, const int natomax);
+void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, double densmat2[*natom][*natom][*natom][*natom], const int natomax);
diff --git a/src/get_ntot.c b/src/get_ntot.c
new file mode 100644
index 0000000..27b7e17
--- /dev/null
+++ b/src/get_ntot.c
@@ -0,0 +1,40 @@
+#include "get_ntot.h"
+
+int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2){
+  int            tnt1, tnt2;
+  int            natom2;
+  if(FAM1){
+      if(fix_trou1 == fix_trou2){
+          natom2 = natom/2;
+      }
+      else{
+          natom2 = fix_trou2 - fix_trou1;
+      }
+  }
+  else{
+      natom2 = natom;
+  }
+
+  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){
+        nalpha=(natom-ntrou+2*isz)/2;
+        nbeta=(natom -ntrou-2*isz)/2;
+        if(((natom-ntrou)/2) == isz){
+            nbeta=0;
+        }
+    }
+    else{
+        nalpha=(natom-ntrou+2*isz+1)/2;
+        nbeta=(natom -ntrou-2*isz-1)/2;
+        if(((natom-ntrou+1)/2) == isz){
+            nbeta=0;
+        }
+    }
+
+  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
new file mode 100644
index 0000000..67bca0e
--- /dev/null
+++ b/src/get_ntot.h
@@ -0,0 +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/get_s2.c b/src/get_s2.c
new file mode 100644
index 0000000..081213e
--- /dev/null
+++ b/src/get_s2.c
@@ -0,0 +1,685 @@
+#include <stdio.h>
+#include <petsctime.h>
+#include <slepceps.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include "get_s2.h"
+#include "get_val_iaa2.h"
+
+/* 
+ * This function simply calculates the S^2 value of the wavefunction
+ * Input
+ * =====
+ * Vr       = The full vector
+ * Istart   = Local starting id of the vector
+ * Iend     = Local vector ending id
+ * valxr    = Local vector values
+ * natom    = number of orbitals
+ * Output
+ * ======
+ * norm     = norm of the vector 
+ * xymat    = the S^2 value
+ */
+
+void get_s2(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *natom, 
+        PetscReal *norm, PetscReal *norm2, PetscReal *norm3, PetscReal *norm4, 
+        PetscReal *xymat, PetscReal *xymat2, PetscReal *xymat3, PetscReal *xymat4, PetscReal *weight3,
+            int *s21a1, int *s21a2, int *s21b1,  int *s21b2,  int *s22a1,  int *s22a2,  
+            int *s22b1, int *s22b2, int *s23a1,  int *s23a2,  int *s23b1, int *s23b2, int *postrou, const int natomax){
+  long int       iaa2, iaa;
+  long int       iii;
+  int			 ideter[natomax];
+  int			 ideter2[natomax];
+  int 		   	 kko,kok,kkio;
+  long int       ii;
+  double         xmat=0.0;
+  double         xmat2=0.0;
+  double         xmat3=0.0;
+  double         xmat4=0.0;
+  double         getvaliaa2;
+  PetscLogDouble t1,t2,tt1,tt2;
+  PetscErrorCode ierr;
+  PetscInt      iiii;
+  int            ntrouboit1=0;
+  int            ntrouboit2=0;
+  int            ntrouboit3=0;
+  int            okboit1=0;
+  int            okboit2=0;
+  int            okboit3=0;
+  int            mpiid;
+  int            pos1=0;
+  int            pos2=0;
+  int            pos3=0;
+  MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
+//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);}
+//ierr = PetscTime(&tt1);CHKERRQ(ierr);
+  	  for(ii=*Istart;ii<*Iend;ii++) {
+              iii = ii + 1;
+//            iiii = ii-*Istart;
+              iiii = ii;
+              xmat = 0.0;
+              xmat2 = 0.0;
+              xmat3 = 0.0;
+              xmat4 = 0.0;
+              ntrouboit1 = 0;
+              ntrouboit2 = 0;
+              ntrouboit3 = 0;
+              okboit1 = 0;
+              okboit2 = 0;
+              okboit3 = 0;
+              pos1 = 0;
+              pos2 = 0;
+              pos3 = 0;
+              getdet_(&iii, ideter);
+              *norm=*norm+valxr[iiii]*valxr[iiii];
+              for(kko=*s21a1;kko<=*s21a2;kko++){
+                  if(ideter[kko]==3){
+                      ntrouboit1++;
+                      pos1=kko;
+                  }
+              }
+              for(kko=*s22a1;kko<=*s22a2;kko++){
+                  if(ideter[kko]==3){
+                      ntrouboit2++;
+                      pos2=kko;
+                  }
+              }
+              for(kko=*s23a1;kko<=*s23a2;kko++){
+                  if(ideter[kko]==3){
+                      ntrouboit3++;
+                      pos3=kko;
+                  }
+              }
+              if(ntrouboit1==1 && pos1 == *postrou)okboit1=1;
+              if(ntrouboit2==1 && pos2 == *postrou)okboit2=1;
+              if(ntrouboit3==1 && pos3 == *postrou)okboit3=1;
+              if(okboit1){
+                *norm2=*norm2+valxr[iiii]*valxr[iiii];
+              }
+              if(okboit2){
+                *norm3=*norm3+valxr[iiii]*valxr[iiii];
+              }
+              if(okboit3){
+                *norm4=*norm4+valxr[iiii]*valxr[iiii];
+              }
+              /*
+               * calculate the weight of ms=5/2
+               *
+               * loop over the determinants to see if we have a S=5/2
+               */
+              int countw = 0;
+              for(kko=*s21a1;kko<=*s21a2;kko++){
+                  if(ideter[kko] == 2) countw=1;
+              }
+              for(kok=*s21b1;kok<=*s21b2;kok++){
+                  if(ideter[kok] == 2) countw=1;
+              }
+              if(countw==0 && okboit1){
+                  *weight3 += (valxr[iiii]*valxr[iiii]);
+              }
+              for(kko=0;kko<=(*natom/2)-1;kko++){
+                  for(kok=kko;kok<=(*natom/2)-1;kok++){
+                    if(kok == kko && ideter[kok] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21a1 && kok <=*s21a2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22a1 && kok <=*s22a2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23a1 && kok <=*s23a2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko] == 1 && ideter[kok] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko] == 1 && ideter[kok] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=2;
+                        ideter2[kok]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=1;
+                        ideter2[kok]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                       }
+                    }
+                }
+              }
+              for(kko=(*natom/2);kko<=*natom-1;kko++){
+                  for(kok=kko;kok<=*natom-1;kok++){
+                    if(kok == kko && ideter[kok] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko] == 1 && ideter[kok] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 1 && ideter[kok] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=2;
+                        ideter2[kok]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=1;
+                        ideter2[kok]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                       }
+                    }
+                }
+              }
+              for(kko=0;kko<=(*natom/2)-1;kko++){
+                  for(kok=(*natom/2);kok<=*natom-1;kok++){
+                    if(kok == kko && ideter[kok] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko] == 1 && ideter[kok] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 1 && ideter[kok] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=2;
+                        ideter2[kok]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=1;
+                        ideter2[kok]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+//                      if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);}
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                       }
+                    }
+                  }
+                }
+                *xymat=*xymat+xmat;
+                *xymat2=*xymat2+xmat2;
+                *xymat3=*xymat3+xmat3;
+                *xymat4=*xymat4+xmat4;
+//              if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3);
+          }
+
+  ierr = PetscTime(&tt2);
+//printf(" norm = %18f weight = %18f weight/N = %18f tmpwe = %18f\n", *norm2, *weight3, *weight3/(*norm2),tmpwe);
+//printf(" norm = %18f %18f xymat = %18f %18f | %d %d %d %d %d\n", *norm, *norm3, *xymat, *xymat3, *s22a1, *s22a2, *s22b1, *s22b2, *postrou);
+//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr);
+}
diff --git a/src/get_s2.h b/src/get_s2.h
new file mode 100644
index 0000000..facd7a2
--- /dev/null
+++ b/src/get_s2.h
@@ -0,0 +1,38 @@
+#include <stdio.h>
+#include <slepceps.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+
+void get_s2(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *natom, PetscReal *norm, PetscReal *norm2, PetscReal *norm3, PetscReal *norm4, PetscReal *xymat, PetscReal *xymat2, PetscReal *xymat3, PetscReal *xymat4, PetscReal *weight3,
+            int *s21a1, int *s21a2, int *s21b1,  int *s21b2,  int *s22a1,  int *s22a2,  int *s22b1, int *s22b2, int *s23a1,  int *s23a2,  int *s23b1, int *s23b2, int *postrou, const int natomax);
+
+void get_s2_mov(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *,  PetscReal *, PetscReal *, PetscReal *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *, const int natomax);
+
+void get_s2_cyclic(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *,  PetscReal *, PetscReal *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *,
+        int *, const int natomax);
diff --git a/src/get_s2_cyclic.c b/src/get_s2_cyclic.c
new file mode 100644
index 0000000..a02bf2d
--- /dev/null
+++ b/src/get_s2_cyclic.c
@@ -0,0 +1,783 @@
+#include <stdio.h>
+#include <petsctime.h>
+#include <slepceps.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include "get_s2.h"
+#include "get_val_iaa2.h"
+
+/* 
+ * This function simply calculates the S^2 value of the wavefunction
+ * Input
+ * =====
+ * Vr       = The full vector
+ * Istart   = Local starting id of the vector
+ * Iend     = Local vector ending id
+ * valxr    = Local vector values
+ * natom    = number of orbitals
+ * Output
+ * ======
+ * norm     = norm of the vector 
+ * xymat    = the S^2 value
+ */
+
+void get_s2_cyclic(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *natom, PetscReal *norm, PetscReal *norm2, PetscReal *norm3, PetscReal *norm4, PetscReal *xymat, PetscReal *xymat2, PetscReal *xymat3, PetscReal *xymat4,
+            int *s21a1, int *s21a2, int *s21b1,  int *s21b2,  int *s22a1,  int *s22a2,  int *s22b1, int *s22b2, int *s23a1,  int *s23a2,  int *s23b1, int *s23b2, int *postrou, const int natomax){
+  long int       iaa2, iaa;
+  long int       iii;
+  int			 ideter[natomax];
+  int			 ideter2[natomax];
+  int 		   	 kko,kok,kkio,kk;
+  int 		   	 kko2,kok2;
+  long int       ii;
+  double         xmat=0.0;
+  double         xmat2=0.0;
+  double         xmat3=0.0;
+  double         xmat4=0.0;
+  double         getvaliaa2;
+  PetscLogDouble t1,t2,tt1,tt2;
+  PetscErrorCode ierr;
+  PetscInt      iiii;
+  int            ntrouboit1=0;
+  int            ntrouboit2=0;
+  int            ntrouboit3=0;
+  int            okboit1=0;
+  int            okboit2=0;
+  int            okboit3=0;
+  int            mpiid;
+  int            pos1=0;
+  int            pos2=0;
+  int            pos3=0;
+  MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
+//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);}
+//ierr = PetscTime(&tt1);CHKERRQ(ierr);
+  	  for(ii=*Istart;ii<*Iend;ii++) {
+              iii = ii + 1;
+//            iiii = ii-*Istart;
+              iiii = ii;
+              xmat = 0.0;
+              xmat2 = 0.0;
+              xmat3 = 0.0;
+              xmat4 = 0.0;
+              ntrouboit1 = 0;
+              ntrouboit2 = 0;
+              ntrouboit3 = 0;
+              okboit1 = 0;
+              okboit2 = 0;
+              okboit3 = 0;
+              pos1 = 0;
+              pos2 = 0;
+              pos3 = 0;
+              getdet_(&iii, ideter);
+              *norm=*norm+valxr[iiii]*valxr[iiii];
+
+              for(kko=0;kko<*natom/2;kko++){
+                  if(ideter[kko]==3){
+                      kk=kko;
+                  }
+              }
+
+              *postrou  = kk;
+              *s21a1    = kk-1;
+//             if(*s21a1<0){
+//                 *s21a1 = (*natom/2) + (*s21a1)%(*natom/2);
+//             }
+//             else{
+//                 *s21a1 = (*s21a1)%(*natom/2);
+//             }
+               *s22a1    = kk-1;
+//             if(*s22a1<0){
+//                 *s22a1 = (*natom/2) + (*s22a1)%(*natom/2);
+//             }
+//             else{
+//                 *s22a1 = (*s22a1)%(*natom/2);
+//             }
+               *s23a1    = kk-2;
+//             if(*s23a1<0){
+//                 *s23a1 = (*natom/2) + (*s23a1)%(*natom/2);
+//             }
+//             else{
+//                 *s23a1 = (*s23a1)%(*natom/2);
+//             }
+//
+               *s21a2    = kk+1;
+//             if(*s21a2<0){
+//                 *s21a2 = (*natom/2) + (*s21a2)%(*natom/2);
+//             }
+//             else{
+//                 *s21a2 = (*s21a2)%(*natom/2);
+//             }
+               *s22a2    = kk+2;
+//             if(*s22a2<0){
+//                 *s22a2 = (*natom/2) + (*s22a2)%(*natom/2);
+//             }
+//             else{
+//                 *s22a2 = (*s22a2)%(*natom/2);
+//             }
+               *s23a2    = kk+2;
+//             if(*s23a2<0){
+//                 *s23a2 = (*natom/2) + (*s23a2)%(*natom/2);
+//             }
+//             else{
+//                 *s23a2 = (*s23a2)%(*natom/2);
+//             }
+//
+               *s21b1 =  *natom + *s21a1;
+               *s22b1 =  *natom + *s22a1;
+               *s23b1 =  *natom + *s23a1;
+
+               *s21b2 =  *natom + *s21a2;
+               *s22b2 =  *natom + *s22a2;
+               *s23b2 =  *natom + *s23a2;
+
+//            if(mpiid==0)printf("postrou = %d\n",*postrou);
+//            if(mpiid==0)printf("1a1 = %d, 1a2 = %d, 1b1 = %d, 1b2 = %d\n",*s21a1,*s21a2,*s21b1,*s21b2);
+//            if(mpiid==0)printf("2a1 = %d, 2a2 = %d, 2b1 = %d, 2b2 = %d\n",*s22a1,*s22a2,*s22b1,*s22b2);
+//            if(mpiid==0)printf("3a1 = %d, 3a2 = %d, 3b1 = %d, 3b2 = %d\n",*s23a1,*s23a2,*s23b1,*s23b2);
+
+//            for(kko=*s21a1;kko<=*s21a2;kko++){
+//                if(ideter[kko]==3){
+//                    ntrouboit1++;
+//                    pos1=kko;
+//                }
+//            }
+//            for(kko=*s22a1;kko<=*s22a2;kko++){
+//                if(ideter[kko]==3){
+//                    ntrouboit2++;
+//                    pos2=kko;
+//                }
+//            }
+//            for(kko=*s23a1;kko<=*s23a2;kko++){
+//                if(ideter[kko]==3){
+//                    ntrouboit3++;
+//                    pos3=kko;
+//                }
+//            }
+//            if(ntrouboit1==1 && pos1 == *postrou)okboit1=1;
+//            if(ntrouboit2==1 && pos2 == *postrou)okboit2=1;
+//            if(ntrouboit3==1 && pos3 == *postrou)okboit3=1;
+              okboit1 = 1;
+              okboit2 = 1;
+              okboit3 = 1;
+              if(okboit1){
+                *norm2=*norm2+valxr[iiii]*valxr[iiii];
+              }
+              if(okboit2){
+                *norm3=*norm3+valxr[iiii]*valxr[iiii];
+              }
+              if(okboit3){
+                *norm4=*norm4+valxr[iiii]*valxr[iiii];
+              }
+              for(kko=kk-2;kko<=kk+(*natom/2-3);kko++){
+                  for(kok=kko;kok<=kk+(*natom/2-3);kok++){
+                    kko2=kko;
+                    if(kko2<0){
+                        kko2 = (*natom/2) + (kko2)%(*natom/2);
+                    }
+                    else{
+                        kko2 = (kko2)%(*natom/2);
+                    }
+                    kok2=kok;
+                    if(kok2<0){
+                        kok2 = (*natom/2) + (kok2)%(*natom/2);
+                    }
+                    else{
+                        kok2 = (kok2)%(*natom/2);
+                    }
+
+                    if(kok == kko && ideter[kok2] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21a1 && kok <=*s21a2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22a1 && kok <=*s22a2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23a1 && kok <=*s23a2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko2] == 1 && ideter[kok2] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko2] == 2 && ideter[kok2] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko2] == 1 && ideter[kok2] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko2]=2;
+                        ideter2[kok2]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko2] == 2 && ideter[kok2] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko2]=1;
+                        ideter2[kok2]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                       }
+                    }
+                }
+              }
+              for(kko=*natom+kk-2;kko<=*natom+kk+(*natom/2-3);kko++){
+                  for(kok=kko;kok<=*natom+kk+(*natom/2-3);kok++){
+                    kko2=kko-*natom;
+                    if(kko2<0){
+                        kko2 = (*natom/2) + (kko2)%(*natom/2);
+                    }
+                    else{
+                        kko2 = (kko2)%(*natom/2);
+                    }
+                    kok2=kok-*natom;
+                    if(kok2<0){
+                        kok2 = (*natom/2) + (kok2)%(*natom/2);
+                    }
+                    else{
+                        kok2 = (kok2)%(*natom/2);
+                    }
+                    kko2 = (*natom) - 1 - kko2;
+                    kok2 = (*natom) - 1 - kok2;
+
+                    if(kok == kko && ideter[kok2] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko2] == 1 && ideter[kok2] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko2] == 2 && ideter[kok2] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko2] == 1 && ideter[kok2] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko2]=2;
+                        ideter2[kok2]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko2] == 2 && ideter[kok2] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko2]=1;
+                        ideter2[kok2]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                       }
+                    }
+                }
+              }
+              for(kko=kk-2;kko<=kk+(*natom/2-3);kko++){
+                  for(kok=*natom + kk-2;kok<=*natom + kk+(*natom/2-3);kok++){
+                    kko2=kko;
+                    if(kko2<0){
+                        kko2 = (*natom/2) + (kko2)%(*natom/2);
+                    }
+                    else{
+                        kko2 = (kko2)%(*natom/2);
+                    }
+                    kok2=kok-*natom;
+                    if(kok2<0){
+                        kok2 = (*natom/2) + (kok2)%(*natom/2);
+                    }
+                    else{
+                        kok2 = (kok2)%(*natom/2);
+                    }
+                    kok2 = (*natom) - 1 - kok2;
+
+                    if(kok == kko && ideter[kok2] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko2] == 1 && ideter[kok2] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko2] == 2 && ideter[kok2] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko2] == 1 && ideter[kok2] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko2]=2;
+                        ideter2[kok2]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko2] == 2 && ideter[kok2] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko2]=1;
+                        ideter2[kok2]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+//                      if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);}
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                       }
+                    }
+                  }
+                }
+                *xymat=*xymat+xmat;
+                *xymat2=*xymat2+xmat2;
+                *xymat3=*xymat3+xmat3;
+                *xymat4=*xymat4+xmat4;
+//              if(mpiid==0)printf(" ii = %d xmat3 = %18f xmat4 = %18f diff = %18f\n", ii, xmat3, xmat4, (xmat3-xmat4));
+          }
+
+  ierr = PetscTime(&tt2);
+//if(mpiid==0)printf(" norm3 = %18f norm4 = %18f xymat3= %18f xymat4= %18f\n", *norm3, *norm4, *xymat3, *xymat4);
+//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr);
+}
diff --git a/src/get_s2_mov.c b/src/get_s2_mov.c
new file mode 100644
index 0000000..5470745
--- /dev/null
+++ b/src/get_s2_mov.c
@@ -0,0 +1,682 @@
+#include <stdio.h>
+#include <petsctime.h>
+#include <slepceps.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include "get_s2.h"
+#include "get_val_iaa2.h"
+
+/* 
+ * This function simply calculates the S^2 value of the wavefunction
+ * Input
+ * =====
+ * Vr       = The full vector
+ * Istart   = Local starting id of the vector
+ * Iend     = Local vector ending id
+ * valxr    = Local vector values
+ * natom    = number of orbitals
+ * Output
+ * ======
+ * norm     = norm of the vector 
+ * xymat    = the S^2 value
+ */
+
+void get_s2_mov(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *natom, PetscReal *norm, PetscReal *norm2, PetscReal *norm3, PetscReal *norm4, PetscReal *xymat, PetscReal *xymat2, PetscReal *xymat3, PetscReal *xymat4, PetscReal *weight3,
+            int *s21a1, int *s21a2, int *s21b1,  int *s21b2,  int *s22a1,  int *s22a2,  int *s22b1, int *s22b2, int *s23a1,  int *s23a2,  int *s23b1, int *s23b2, int *postrou, const int natomax){
+  long int       iaa2, iaa;
+  long int       iii;
+  int			 ideter[natomax];
+  int			 ideter2[natomax];
+  int 		   	 kko,kok,kkio;
+  long int       ii;
+  double         xmat=0.0;
+  double         xmat2=0.0;
+  double         xmat3=0.0;
+  double         xmat4=0.0;
+  double         getvaliaa2;
+  PetscLogDouble t1,t2,tt1,tt2;
+  PetscErrorCode ierr;
+  PetscInt      iiii;
+  int            ntrouboit1=0;
+  int            ntrouboit2=0;
+  int            ntrouboit3=0;
+  int            okboit1=0;
+  int            okboit2=0;
+  int            okboit3=0;
+  int            mpiid;
+  int            pos1=0;
+  int            pos2=0;
+  int            pos3=0;
+  MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
+//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);}
+//ierr = PetscTime(&tt1);CHKERRQ(ierr);
+  	  for(ii=*Istart;ii<*Iend;ii++) {
+              iii = ii + 1;
+//            iiii = ii-*Istart;
+              iiii = ii;
+              xmat = 0.0;
+              xmat2 = 0.0;
+              xmat3 = 0.0;
+              xmat4 = 0.0;
+              ntrouboit1 = 0;
+              ntrouboit2 = 0;
+              ntrouboit3 = 0;
+              okboit1 = 0;
+              okboit2 = 0;
+              okboit3 = 0;
+              pos1 = 0;
+              pos2 = 0;
+              pos3 = 0;
+              getdet_(&iii, ideter);
+              *norm=*norm+valxr[iiii]*valxr[iiii];
+              for(kko=*s21a1;kko<=*s21a2;kko++){
+                  if(ideter[kko]==3){
+                      ntrouboit1++;
+                      pos1=kko;
+                  }
+              }
+              for(kko=*s22a1;kko<=*s22a2;kko++){
+                  if(ideter[kko]==3){
+                      ntrouboit2++;
+                      pos2=kko;
+                  }
+              }
+              for(kko=*s23a1;kko<=*s23a2;kko++){
+                  if(ideter[kko]==3){
+                      ntrouboit3++;
+                      pos3=kko;
+                  }
+              }
+              if(ntrouboit1==1 && *s21a1 <= pos1 && pos1 <= *s21a2)okboit1=1;
+              if(ntrouboit2==1 && pos2 == *postrou)okboit2=1;
+              if(ntrouboit3==1 && pos3 == *postrou)okboit3=1;
+              if(okboit1){
+                *norm2=*norm2+valxr[iiii]*valxr[iiii];
+              }
+              if(okboit2){
+                *norm3=*norm3+valxr[iiii]*valxr[iiii];
+              }
+              if(okboit3){
+                *norm4=*norm4+valxr[iiii]*valxr[iiii];
+              }
+              /*
+               * calculate the weight of ms=5/2
+               *
+               * loop over the determinants to see if we have a S=5/2
+               */
+              int countw = 0;
+              for(kko=*s21a1;kko<=*s21a2;kko++){
+                  if(ideter[kko] == 2) countw=1;
+              }
+              for(kok=*s21b1;kok<=*s21b2;kok++){
+                  if(ideter[kok] == 2) countw=1;
+              }
+              if(countw==0 && okboit1){
+                  *weight3 += (valxr[iiii]*valxr[iiii]);
+              }
+              for(kko=0;kko<=(*natom/2)-1;kko++){
+                  for(kok=kko;kok<=(*natom/2)-1;kok++){
+                    if(kok == kko && ideter[kok] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21a1 && kok <=*s21a2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22a1 && kok <=*s22a2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23a1 && kok <=*s23a2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko] == 1 && ideter[kok] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko] == 1 && ideter[kok] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=2;
+                        ideter2[kok]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                            }
+                          }
+                        }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=1;
+                        ideter2[kok]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                        if(okboit1){
+                          if( kko >=*s21a1 && kko <=*s21a2){
+                            if( kok >=*s21a1 && kok <=*s21a2){
+                                  xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit2){
+                          if( kko >=*s22a1 && kko <=*s22a2){
+                            if( kok >=*s22a1 && kok <=*s22a2){
+                                  xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                        if(okboit3){
+                          if( kko >=*s23a1 && kko <=*s23a2){
+                            if( kok >=*s23a1 && kok <=*s23a2){
+                                  xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                            }
+                          }
+                        }
+                       }
+                    }
+                }
+              }
+              for(kko=(*natom/2);kko<=*natom-1;kko++){
+                  for(kok=kko;kok<=*natom-1;kok++){
+                    if(kok == kko && ideter[kok] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko] == 1 && ideter[kok] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 1 && ideter[kok] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=2;
+                        ideter2[kok]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=1;
+                        ideter2[kok]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21b1 && kko <=*s21b2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22b1 && kko <=*s22b2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23b1 && kko <=*s23b2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                       }
+                    }
+                }
+              }
+              for(kko=0;kko<=(*natom/2)-1;kko++){
+                  for(kok=(*natom/2);kok<=*natom-1;kok++){
+                    if(kok == kko && ideter[kok] != 3){
+                      xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                    }
+                    else{
+                      if(ideter[kko] == 1 && ideter[kok] == 1){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 2){
+                        xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 1 && ideter[kok] == 2){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=2;
+                        ideter2[kok]=1;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      }
+                      if(ideter[kko] == 2 && ideter[kok] == 1){
+                        xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
+                          }
+                        }
+                      }
+                        for(kkio=0;kkio<=*natom-1;kkio++){
+                          ideter2[kkio]=ideter[kkio];
+                        }
+                        ideter2[kko]=1;
+                        ideter2[kok]=2;
+                        adr_(ideter2, &iaa2);
+                        iaa2 = iaa2 - 1;
+//                      if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);}
+                        xmat=xmat+valxr[iiii]*valxr[iaa2];
+                      if(okboit1){
+                        if( kko >=*s21a1 && kko <=*s21a2){
+                          if( kok >=*s21b1 && kok <=*s21b2){
+                                xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit2){
+                        if( kko >=*s22a1 && kko <=*s22a2){
+                          if( kok >=*s22b1 && kok <=*s22b2){
+                                xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                      if(okboit3){
+                        if( kko >=*s23a1 && kko <=*s23a2){
+                          if( kok >=*s23b1 && kok <=*s23b2){
+                                xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
+                          }
+                        }
+                      }
+                       }
+                    }
+                  }
+                }
+                *xymat=*xymat+xmat;
+                *xymat2=*xymat2+xmat2;
+                *xymat3=*xymat3+xmat3;
+                *xymat4=*xymat4+xmat4;
+//              if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3);
+          }
+
+  ierr = PetscTime(&tt2);
+//printf(" norm = %18f weight = %18f weight/N = %18f tmpwe = %18f\n", *norm2, *weight3, *weight3/(*norm2),tmpwe);
+//printf(" norm = %18f %18f xymat = %18f %18f\n", *norm2, *norm3, *xymat2, *xymat3);
+//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr);
+}
diff --git a/src/get_val_iaa2.c b/src/get_val_iaa2.c
new file mode 100644
index 0000000..e88a35f
--- /dev/null
+++ b/src/get_val_iaa2.c
@@ -0,0 +1,34 @@
+#include <stdio.h>
+#include <petsctime.h>
+#include <slepceps.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+#include "get_val_iaa2.h"
+
+/*
+ * This function gets vector from a different processor
+ * xr           the full vector
+ * iaa2         adresse to get value from
+ * getvaliaa2   the value
+ */
+
+void get_val_iaa2(Vec xr,long int *iaa2,double *getvaliaa2){
+    Vec x; /* initial vector, destination vector */
+    VecScatter scatter; /* scatter context */
+    IS from, to; /* index sets that define the scatter */
+    PetscScalar *values;
+    int idx_from[] = {*iaa2}, idx_to[] = {0}; 
+    VecCreateSeq(   PETSC_COMM_SELF,1,&x);
+    ISCreateGeneral(PETSC_COMM_SELF,1,idx_from,PETSC_COPY_VALUES,&from); 
+    ISCreateGeneral(PETSC_COMM_SELF,1,idx_to,  PETSC_COPY_VALUES,&to); 
+    printf("in get_val");
+    VecScatterCreate(xr,from,x,to,&scatter);
+    VecScatterBegin(scatter,xr,x,INSERT_VALUES,SCATTER_FORWARD); 
+    VecScatterEnd(scatter,xr,x,  INSERT_VALUES,SCATTER_FORWARD); 
+    VecGetArray(x,&values);
+    *getvaliaa2 = values[0];
+    ISDestroy(&from);
+    ISDestroy(&to);
+    VecScatterDestroy(&scatter);
+}
diff --git a/src/get_val_iaa2.h b/src/get_val_iaa2.h
new file mode 100644
index 0000000..f47bcaa
--- /dev/null
+++ b/src/get_val_iaa2.h
@@ -0,0 +1,8 @@
+#include <stdio.h>
+#include <petsctime.h>
+#include <slepceps.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <string.h>
+
+void get_val_iaa2(Vec, long int *, double *);
diff --git a/src/getdet.irp.f b/src/getdet.irp.f
index cd483af..acba73a 100644
--- a/src/getdet.irp.f
+++ b/src/getdet.irp.f
@@ -1,13 +1,14 @@
 subroutine getdet(add,ideter)
+    use iso_c_binding
     implicit none
     BEGIN_DOC
     ! this routing gives the determinant in
     ! the traditional form given its address
     END_DOC
     integer,INTENT(INOUT)::ideter(natomax)
-    integer(kind=selected_int_kind(16)),INTENT(IN)::add
-    integer(kind=selected_int_kind(16))::deta,detb
-    integer::i,const,ia,ib
+    integer(C_SIZE_T),INTENT(IN)::add
+    integer(C_SIZE_T)::deta,detb
+    integer::i,const,ia,ib, natom2
 
     ib = MOD(add,nt2)
     if(MOD(add,nt2).eq.0)then
@@ -20,33 +21,49 @@ subroutine getdet(add,ideter)
     detb=0
     deta=0
     i=1
-    do while (i.le.(ib))
-        const=1
-        do while(popcnt(detb).ne.nbeta .or. const==1)
-            if(nbeta.eq.0)then
-                detb=0
-                EXIT
-            endif
-            detb+=1
-            const=0
-        enddo
-        i+=1
-!       write(6,14)detb,detb
-    enddo
-    i=1
-    do while (i.le.(ia))
-        const=1
-        do while(popcnt(deta).ne.ntrou .or. const==1)
-            deta+=1
-            const=0
-        enddo
-        i+=1
-!       write(6,14)deta,deta
-    enddo
+    detb = det(ib,1)
+    deta = deth(ia,1)
+    if(FAM1) then
+        if(fix_trou1 .eq. fix_trou2) then
+            deta = ISHFT(deta,-(natom/2))
+        else
+            natom2 = natom - (fix_trou2 - fix_trou1)
+            deta = ISHFT(deta, -natom2)
+        endif
+    endif
+!    do while (i.le.(ib))
+!        const=1
+!        do while(popcnt(detb).ne.nbeta .or. const==1)
+!            detb+=1
+!            const=0
+!        enddo
+!        i+=1
+!        write(6,14)detb,detb
+!    enddo
+!    i=1
+!    do while (i.le.(ia))
+!        const=1
+!        do while(popcnt(deta).ne.ntrou .or. const==1)
+!            deta+=1
+!            const=0
+!        enddo
+!        i+=1
+!        write(6,14)deta,deta
+!    enddo
     const=0
-    do i=0,(natom/2) - 1
+    if(FAM1) then
+        if(fix_trou1 .eq. fix_trou2) then
+            natom2 = natom/2
+        else
+            natom2 = (fix_trou2 - fix_trou1)
+        endif
+    else
+        natom2 = natom
+    endif
+
+    do i=0,(natom2) - 1
         if(BTEST(deta,i))then
-            ideter((natom/2)-i)=3
+            ideter((natom2)-i)=3
         endif
     enddo
     do i=0,natom-1
diff --git a/src/main.irp.f b/src/main.irp.f
index 9f26adb..307906a 100644
--- a/src/main.irp.f
+++ b/src/main.irp.f
@@ -1,13 +1,14 @@
         program main
+        use iso_c_binding
         implicit none
-        integer(kind=selected_int_kind(16)),allocatable::tl1(:),tl2(:),tktyp(:)
+        integer(C_SIZE_T),allocatable::tl1(:),tl2(:),tktyp(:)
         real*8,allocatable::txtt(:),txjjxy(:),txjjz(:)
         integer::i, tnrows, tntrou,tisz
         real*4::t1,t2
         real*8,allocatable::tval(:)
-        integer(kind=selected_int_kind(16)),allocatable::tcol(:)
-        integer(kind=selected_int_kind(16)),dimension(10)::tcountcol
-        integer(kind=selected_int_kind(16))::tistart
+        integer(C_SIZE_T),allocatable::tcol(:)
+        integer(C_SIZE_T),dimension(10)::tcountcol
+        integer(C_SIZE_T)::tistart
         allocate (tl1(maxlien))
         allocate (tl2(maxlien))
         allocate (tktyp(maxlien))
diff --git a/src/natom.irp.f b/src/natom.irp.f
index 7b8091a..58bf864 100644
--- a/src/natom.irp.f
+++ b/src/natom.irp.f
@@ -2,7 +2,6 @@ BEGIN_PROVIDER [integer, natom]
 &BEGIN_PROVIDER [integer, natrest]
 &BEGIN_PROVIDER [integer, ial0]
 &BEGIN_PROVIDER [logical*1, yham]
-&BEGIN_PROVIDER [logical*1, FAM1]
 &BEGIN_PROVIDER [integer, nlientot]
 &BEGIN_PROVIDER [real*8, xt,(maxlien)]
 &BEGIN_PROVIDER [real*8 , xjz,(maxlien)]
@@ -93,10 +92,14 @@ BEGIN_PROVIDER [integer, natom]
       enddo
 !------------------Lecture Hamiltonien
 
-       FAM1=.TRUE.
+!      FAM1=.TRUE. 
        yham=.TRUE.
        write(6,*)'HAMILTONIEN t-J'
        write(6,*)'Le nombre de trou est : ',ntrou
+       write(6,*)'Famille 1 : ',FAM1
+       if(FAM1) then
+        if(fix_trou1 .ne. fix_trou2) write(6,*)'Trou fixe entre :', fix_trou1, "et ", fix_trou2
+       endif
 !---------------------------------------------
       write(6,*)' '
       write(6,*)' '
diff --git a/src/natomax.irp.f b/src/natomax.irp.f
index c4c53ce..6f443a9 100644
--- a/src/natomax.irp.f
+++ b/src/natomax.irp.f
@@ -13,10 +13,10 @@ BEGIN_PROVIDER [integer, natomax]
     END_DOC
 
     implicit none
-    natomax=700
+    natomax=900
     jrangmax=10705432
     maxial=20
-    maxlien=700
+    maxlien=900
     maxplac=20
     maxdet=10000
 END_PROVIDER
diff --git a/src/nt1.irp.f b/src/nt1.irp.f
index 64c3e78..8bd0124 100644
--- a/src/nt1.irp.f
+++ b/src/nt1.irp.f
@@ -1,13 +1,20 @@
-BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt1]
+use iso_c_binding
+BEGIN_PROVIDER [integer(C_SIZE_T), nt1]
     BEGIN_DOC
     ! calculates the number of det the 3's moving
     END_DOC
     implicit none
-    integer(kind=selected_int_kind(16))::natom2
+    integer(C_SIZE_T)::natom2
 
 !   call combin(idet1(1,nt1+1),natom,ntrou,nt1,32,jrangmax)
     natom2=natom
-    if(FAM1)natom2=natom/2
-    nt1=   nint(gamma(real(natom2+1,16))/(gamma(real(natom2-ntrou+1,16))*gamma(real(ntrou+1,16))),selected_int_kind(16))
+    if(FAM1) then
+        if(fix_trou1 .eq. fix_trou2) then
+            natom2=natom/2
+        else
+            natom2 = fix_trou2 - fix_trou1
+        endif
+    endif
+    nt1=   nint(gamma(1.0*(natom2+1))/(gamma(1.0*(natom2-ntrou+1))*gamma(1.0*(ntrou+1))),selected_int_kind(16))
     write(6,*)'nt1',nt1
 END_PROVIDER
diff --git a/src/nt2.irp.f b/src/nt2.irp.f
index 05e31c2..0ff4c5d 100644
--- a/src/nt2.irp.f
+++ b/src/nt2.irp.f
@@ -1,10 +1,11 @@
-BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt2]
+use iso_c_binding
+BEGIN_PROVIDER [integer(C_SIZE_T), nt2]
     BEGIN_DOC
     ! calculates the number of det the 1's moving
     END_DOC
 
 
 !   call combin(idet2(1,nt2+1),natrest,ial0,nt2,32,jrangmax)
-    nt2=   nint(gamma(real(natom-ntrou+1,16))/((gamma(real(nalpha+1,16))*gamma(real(nbeta+1,16)))),selected_int_kind(16))
+    nt2=   nint(gamma(1.0*(natom-ntrou+1))/((gamma(1.0*(nalpha+1))*gamma(1.0*(nbeta+1)))),selected_int_kind(16))
     print *,"nt2=",nt2
 END_PROVIDER
diff --git a/src/providdet.irp.f b/src/providdet.irp.f
new file mode 100644
index 0000000..d651244
--- /dev/null
+++ b/src/providdet.irp.f
@@ -0,0 +1,91 @@
+use iso_c_binding
+BEGIN_PROVIDER[integer(C_SIZE_T),det,(nt2,2)]
+&BEGIN_PROVIDER[integer(C_SIZE_T),deth,(nt1,2)]
+    BEGIN_DOC
+    ! provides det and deth array
+    END_DOC
+    use iso_c_binding
+    implicit none
+!   integer(kind=selected_int_kind(16))::dethsh
+    integer(C_SIZE_t)::a
+    integer(C_SIZE_T)::i,count
+    integer::const
+    i=1
+    a=0
+    const=0
+    count=0
+
+    If(ntrou.ge.1)then
+
+            const=0
+!           dethsh = ISHFT(deth,-natom/2)
+!           i=nt1
+            do while (i.le.(nt1))
+!               if(a.eq.dethsh)then
+!                   addh=i-1
+!                   EXIT
+!               endif
+
+                i+=1
+                a+=1
+                do while(popcnt(a).ne.ntrou)
+                    a+=1
+                enddo
+                count+=1
+                if(FAM1) then
+                    if(fix_trou1 .eq. fix_trou2) then
+                        deth(count,1)=ISHFT(a,natom/2)
+                    else
+                        deth(count,1)=ISHFT(a,natom - (fix_trou2-fix_trou1))
+                    endif
+                else 
+                    deth(count, 1) = a
+                endif
+                deth(count,2)=i-1
+!       write(6,16)ISHFT(a,natom/2),ISHFT(a,natom/2),i-1
+            enddo
+!           if(a.eq.dethsh )then
+!           count+=1
+!           deth(1,1)=ISHFT(a,natom/2)
+!           deth(1,2)=nt1
+!           endif
+
+    endif
+
+    !C if det=0 then exit
+    a=0
+    i=0
+    count=0
+    print *,'nt2=',nt2,'nbeta=',nbeta
+    do while (i.lt.(nt2))
+
+        i+=1
+        a+=1
+        do while(popcnt(a).ne.nbeta)
+            if(nbeta.eq.0)then
+                a=0
+                count+=1
+                det(count,1)=a
+                det(count,2)=i
+                EXIT
+            endif
+            a+=1
+        enddo
+        
+        if(nbeta.ne.0)then
+        count+=1
+        det(count,1)=a
+        det(count,2)=i
+        endif
+!       write(6,16)a,a,i
+    enddo
+    
+
+10  FORMAT(B64,I8,F8.2)
+15  FORMAT(B64,I8,I8,I8)
+11  FORMAT(B64,I3,B64)
+12  FORMAT(I5,$)
+13  FORMAT(B64,B64)
+14  FORMAT(B64,I8)
+16  FORMAT(B64,I8,I8)
+END_PROVIDER
diff --git a/src/provide_clust.irp.f b/src/provide_clust.irp.f
index fe5e8ef..4f7e96e 100644
--- a/src/provide_clust.irp.f
+++ b/src/provide_clust.irp.f
@@ -1,11 +1,15 @@
 BEGIN_PROVIDER[integer,l1,   (maxlien)]
 &BEGIN_PROVIDER[integer,l2,  (maxlien)]
 &BEGIN_PROVIDER[integer,ktyp,(maxlien)]
-&BEGIN_PROVIDER [real*8, xtt  ,(maxlien)]
+&BEGIN_PROVIDER[real*8, xtt  ,(maxlien)]
 &BEGIN_PROVIDER[real*8, xjjz ,(maxlien)]
 &BEGIN_PROVIDER[real*8, xjjxy,(maxlien)]
-&BEGIN_PROVIDER [integer, ntrou]
-&BEGIN_PROVIDER [integer, isz]
+&BEGIN_PROVIDER[real*8, E,(maxlien)]
+&BEGIN_PROVIDER[integer, ntrou]
+&BEGIN_PROVIDER[integer, isz]
+&BEGIN_PROVIDER[logical*1, FAM1]
+&BEGIN_PROVIDER[integer, fix_trou1]
+&BEGIN_PROVIDER[integer, fix_trou2]
       implicit none
 !     integer::i
 !       open(unit=11,file="l1.dat",form="formatted")
diff --git a/src/read2.c b/src/read2.c
index 70836dc..8f8b040 100644
--- a/src/read2.c
+++ b/src/read2.c
@@ -1,7 +1,3 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <ctype.h>
-#include <string.h>
 #include "read2.h"
 
 void Data_new(FILE* file, Data* dat) { 
@@ -15,13 +11,13 @@ 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 != 11){
+	    /* 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:
-					dat->n=atol(line);
+					dat->natom=atol(line);
 					break;
 				case 2:
 					dat->nnz=atol(line);
@@ -36,6 +32,9 @@ void Data_new(FILE* file, Data* dat) {
 					dat->isz=atol(line);
 					break;
 				case 6:
+					dat->FAM1 = to_bool(line);
+					break;
+				case 7:
 					arrayIdx=0;
 					for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
 					{
@@ -59,7 +58,7 @@ void Data_new(FILE* file, Data* dat) {
 					  }
 					}
 					break;
-				case 7:
+				case 8:
 					arrayIdx=0;
 					for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
 					{
@@ -83,7 +82,7 @@ void Data_new(FILE* file, Data* dat) {
 					  }
 					}
 					break;
-				case 8:
+				case 9:
 					arrayIdx=0;
 					for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
 					{
@@ -107,7 +106,7 @@ void Data_new(FILE* file, Data* dat) {
 					  }
 					}
 					break;
-				case 9:
+				case 10:
 					arrayIdx=0;
 					for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
 					{
@@ -131,7 +130,7 @@ void Data_new(FILE* file, Data* dat) {
 					  }
 					}
 					break;
-				case 10:
+				case 11:
 					arrayIdx=0;
 					for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
 					{
@@ -155,7 +154,7 @@ void Data_new(FILE* file, Data* dat) {
 					  }
 					}
 					break;
-				case 11:
+				case 12:
 					arrayIdx=0;
 					for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
 					{
@@ -179,6 +178,84 @@ void Data_new(FILE* file, Data* dat) {
 					  }
 					}
 					break;
+				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->nroots=atol(line);
+					break;
+                case 15:
+					dat->s21a1=atol(line);
+					break;
+                case 16:
+					dat->s21a2=atol(line);
+					break;
+                case 17:
+					dat->s21b1=atol(line);
+					break;
+                case 18:
+					dat->s21b2=atol(line);
+					break;
+                case 19:
+					dat->s22a1=atol(line);
+					break;
+                case 20:
+					dat->s22a2=atol(line);
+					break;
+                case 21:
+					dat->s22b1=atol(line);
+					break;
+                case 22:
+					dat->s22b2=atol(line);
+					break;
+                case 23:
+					dat->s23a1=atol(line);
+					break;
+                case 24:
+					dat->s23a2=atol(line);
+					break;
+                case 25:
+					dat->s23b1=atol(line);
+					break;
+                case 26:
+					dat->s23b2=atol(line);
+					break;
+                case 27:
+					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 */
@@ -188,6 +265,17 @@ void Data_new(FILE* file, Data* dat) {
 //return dat;
 }
 
+_Bool to_bool(const char* str) {
+    PetscBool strflg=PETSC_FALSE;
+    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;
+}
+
 /*
 int main(int argc, char* argv[])
 {
diff --git a/src/read2.h b/src/read2.h
index 6c98ab4..9a1a6ff 100644
--- a/src/read2.h
+++ b/src/read2.h
@@ -1,19 +1,43 @@
 #include <stdio.h>
-#include <slepceps.h>
 #include <stdlib.h>
 #include <ctype.h>
 #include <string.h>
 
+#include <petscsys.h>
+#include <slepceps.h>
+
+_Bool to_bool(const char* str);
+
 typedef struct {
 	PetscInt n;
 	long int nnz,npar;
 	long int ntrou,isz;
-	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;
+	int s21a2;
+	int s21b1;
+	int s21b2;
+	int s22a1;
+	int s22a2;
+	int s22b1;
+	int s22b2;
+	int s23a1;
+	int s23a2;
+	int s23b1;
+	int s23b2;
+	int postrou;
+	long int fix_trou1;
+	long int fix_trou2;
+    long int print_wf;
 
 } Data ;
 
diff --git a/src/searchdet.irp.f b/src/searchdet.irp.f
index d900984..66e9b6d 100644
--- a/src/searchdet.irp.f
+++ b/src/searchdet.irp.f
@@ -1,80 +1,71 @@
-subroutine searchdet(det,add,deth,addh)
+subroutine searchdet(deti,add,dethi,addh)
+    use iso_c_binding
     BEGIN_DOC
     ! this subroutine is at the heart of the idea
     ! it will generate all the determinants in a fixed order
     ! then find the posistion of the determinant given and
     ! return it's position in add.
     END_DOC
-    integer(kind=selected_int_kind(16)),INTENT(INOUT)::det
-    integer(kind=selected_int_kind(16)),INTENT(INOUT)::add
-    integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth
-    integer(kind=selected_int_kind(16)),INTENT(INOUT)::addh
-    integer(kind=selected_int_kind(16))::dethsh
-    integer(kind=selected_int_kind(16))::a
-    integer(kind=selected_int_kind(16))::i
-    integer::const
+    integer(C_SIZE_T),INTENT(INOUT)::deti
+    integer(C_SIZE_T),INTENT(INOUT)::add
+    integer(C_SIZE_T),INTENT(INOUT)::dethi
+    integer(C_SIZE_T),INTENT(INOUT)::addh
+    integer(C_SIZE_T)::dethsh
+    integer(C_SIZE_T)::a
+    integer(C_SIZE_T)::i,j
+    integer::count
+    logical::found
+    
     i=1
-    a=0
-    add=0
-    const=0
-
-    If(ntrou.ge.1)then
-
-            const=0
-            dethsh = ISHFT(deth,-natom/2)
-            addh=0
-!           i=nt1
-            do while (i.le.(nt1))
-                if(a.eq.dethsh)then
-                    addh=i-1
-                    EXIT
-                endif
-
-                i+=1
-                a+=1
-                do while(popcnt(a).ne.ntrou)
-                    a+=1
-                enddo
-            enddo
-            if(a.eq.dethsh .and. addh.eq.0)then
-            addh=nt1
+    j=nt1
+    found=.FALSE.
+    do while(.not.found)
+        if(deth((i+j)/2,1).eq.dethi)then
+            addh=deth((i+j)/2,2)
+            found=.TRUE.
+            EXIT
+        elseif (abs(i-j).eq.1)then
+            if(deth(i,1).eq.dethi)then
+            addh=deth(i,2)
+            elseif(deth(j,1).eq.dethi)then
+            addh=deth(j,2)
             endif
 
-    endif
-
-    !C if det=0 then exit
-    a=0
-    i=0
-    count=0
-    if(a.eq.det)then
-    add=1
-    Return
-    endif
-
-    do while (i.le.(nt2))
-        if(a.eq.det)then
-            if(a.eq.1)then
-               add=i
-               count=-1
-               EXIT
+            found=.TRUE.
+            EXIT
+        else
+            if(deth((i+j)/2,1).gt.dethi)then
+                j=(i+j)/2
             else
-               add=i
-               count=-1
-               EXIT
+                i=(i+j)/2
+            endif
+        endif
+    enddo
+
+    i=1
+    j=nt2
+    found=.FALSE.
+    do while(.not.found)
+        if(det((i+j)/2,1).eq.deti)then
+            add=det((i+j)/2,2)
+            found=.TRUE.
+            EXIT
+        elseif (abs(i-j).eq.1)then
+            if(det(i,1).eq.deti)then
+            add=det(i,2)
+            elseif(det(j,1).eq.deti)then
+            add=det(j,2)
+            endif
+            found=.TRUE.
+            EXIT
+        else
+            if(det((i+j)/2,1).gt.deti)then
+                j=(i+j)/2
+            else
+                i=(i+j)/2
             endif
         endif
-
-        i+=1
-        a+=1
-!C      write(6,16)a,a,i-2
-        do while(popcnt(a).ne.nbeta)
-            a+=1
-        enddo
     enddo
-    if(a.eq.det .and. count.ne.-1)then
-    add=i-1
-    endif
-    
 
 10  FORMAT(B64,I8,F8.2)
 15  FORMAT(B64,I8,I8,I8)
diff --git a/src/searchdetfull.irp.f b/src/searchdetfull.irp.f
index 362c527..0906e47 100644
--- a/src/searchdetfull.irp.f
+++ b/src/searchdetfull.irp.f
@@ -1,17 +1,18 @@
 subroutine searchdetfull()
+    use iso_c_binding
     BEGIN_DOC
     ! this subroutine is at the heart of the idea
     ! it will generate all the determinants in a fixed order
     ! then find the posistion of the determinant given and
     ! return it's position in add.
     END_DOC
-!   integer(kind=selected_int_kind(16)),INTENT(INOUT)::foundetadr(maxlien,4)
-    integer(kind=selected_int_kind(16))::add
-!   integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth
-    integer(kind=selected_int_kind(16))::dethsh
-    integer(kind=selected_int_kind(16))::addh
-    integer(kind=selected_int_kind(16))::a
-    integer(kind=selected_int_kind(16))::i
+!   integer(C_SIZE_T),INTENT(INOUT)::foundetadr(maxlien,4)
+    integer(C_SIZE_T)::add
+!   integer(C_SIZE_T),INTENT(INOUT)::deth
+    integer(C_SIZE_T)::dethsh
+    integer(C_SIZE_T)::addh
+    integer(C_SIZE_T)::a
+    integer(C_SIZE_T)::i
     integer::const,count
     i=1
     a=0
diff --git a/src/sort.irp.f b/src/sort.irp.f
index 1c8c130..6a1c1f2 100644
--- a/src/sort.irp.f
+++ b/src/sort.irp.f
@@ -1,14 +1,15 @@
 subroutine sort()
+    use iso_c_binding
     implicit none
     integer::i,j,ord,ordh
-    integer(kind=selected_int_kind(16))::add,addh,det,deth,addt
+    integer(C_SIZE_T)::add,addh,deti,dethi,addt
     
     do i=1,detfound-1
         do j=i+1,detfound
             if(foundaddh(i,1).gt.foundaddh(j,1))then
-                deth = foundaddh(i,1)
+                dethi = foundaddh(i,1)
                 foundaddh(i,1) = foundaddh(j,1)
-                foundaddh(j,1) = deth
+                foundaddh(j,1) = dethi
                 addh = foundaddh(i,2)
                 foundaddh(i,2) = foundaddh(j,2)
                 foundaddh(j,2) = addh
@@ -17,9 +18,9 @@ subroutine sort()
                 foundaddh(j,3) = ordh
             endif
             if(foundadd(i,1).gt.foundadd(j,1))then
-                det = foundadd(i,1)
+                deti = foundadd(i,1)
                 foundadd(i,1) = foundadd(j,1)
-                foundadd(j,1) = det
+                foundadd(j,1) = deti
                 add = foundadd(i,2)
                 foundadd(i,2) = foundadd(j,2)
                 foundadd(j,2) = add
diff --git a/src/stimsyr.h b/src/stimsyr.h
index 84dc62f..5964517 100644
--- a/src/stimsyr.h
+++ b/src/stimsyr.h
@@ -1,3 +1,5 @@
+#include <petscsys.h>
+
 void unit_l1_(
             long int *,
             long int *,
@@ -7,9 +9,13 @@ void unit_l1_(
             double *,
             double *,
             double *,
+            double *,
             long int *,
             long int *,
             long int *,
             long int *,
+            long int *,
+            _Bool *,
+            long int *,
             double *);
 
diff --git a/src/unit_FIL44.irp.f b/src/unit_FIL44.irp.f
index 647554b..96da0c8 100644
--- a/src/unit_FIL44.irp.f
+++ b/src/unit_FIL44.irp.f
@@ -3,19 +3,19 @@
     BEGIN_DOC
     ! file units for writing
     END_DOC
-
+    use iso_c_binding
     implicit none
     integer :: i,j,k,ia1,ia2,l,m,chcind,chcval,ii,tistart2
     integer :: count,unit_44,unit_33
     integer :: iat,nbtots
-    integer(kind=selected_int_kind(16))::iaa
+    integer(C_SIZE_T)::iaa
     integer :: kkio,kkiok,n,nz,cdiag,cexdiag
     integer,allocatable ::ideter1(:),ideter2(:),deti(:),detj(:)
-    integer(kind=selected_int_kind(16)),dimension(maxlien) ::tl1,tl2,tktyp
-    integer(kind=selected_int_kind(16)),dimension(nrows)::tcountcol
-    integer(kind=selected_int_kind(16))::tistart
+    integer(C_SIZE_T),dimension(maxlien) ::tl1,tl2,tktyp
+    integer(C_SIZE_T),dimension(nrows)::tcountcol
+    integer(C_SIZE_T)::tistart
     real*8,dimension(maxlien)::tval
-    integer(kind=selected_int_kind(16)),dimension(maxlien)::tcol
+    integer(C_SIZE_T),dimension(maxlien)::tcol
     real*8 :: xmat
         integer :: ik,imat4,iaa2,iik
         integer :: ik1,ik2,jmat4,IC,ikmax,ikmin
diff --git a/src/unit_l1.irp.f b/src/unit_l1.irp.f
index 5239a49..fc0dc40 100644
--- a/src/unit_l1.irp.f
+++ b/src/unit_l1.irp.f
@@ -6,18 +6,26 @@
         txjjxy,                &
         txjjz ,                &
         txtt  ,                &
+        tE    ,                &
         tcountcol,             &
         tntrou,                &
         tisz,                  &
+        tfix_trou1,            &
+        tfix_trou2,            &
+        tfam1,                 &
         tcol,tval)
+        use iso_c_binding
         implicit none
-        integer,INTENT(INOUT)::tistart, tnrows, tntrou, tisz
+        integer,INTENT(INOUT)::tistart, tnrows
+        integer,INTENT(INOUT)::tntrou, tisz
+        integer,INTENT(INOUT)::tfix_trou1, tfix_trou2
+        logical*1,INTENT(INOUT)::tfam1
         integer::i
         real*8,INTENT(INOUT)::tval(maxlien)
-        integer(kind=selected_int_kind(16)),INTENT(INOUT)::tcol(maxlien)
-        integer(kind=selected_int_kind(16)),INTENT(INOUT),dimension(tnrows)::tcountcol
-        integer(kind=selected_int_kind(16)),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien)
-        real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien)
+        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), tE(maxlien)
    
 		nrows = tnrows
 		provide nrows
@@ -28,9 +36,13 @@
         xtt(i)     = txtt(i)
         xjjxy(i)   = txjjxy(i)
         xjjz (i)   = txjjz (i)
+        E (i)   = tE (i)
         enddo
         ntrou = tntrou
         isz = tisz
+        FAM1 = tfam1
+        fix_trou1 = tfix_trou1
+        fix_trou2 = tfix_trou2
         tcol=0
         tval=0d0
         provide l1 l2 ktyp xtt xjjxy xjjz ntrou
@@ -38,6 +50,7 @@
 !print *,l1
 !print *,"xjjz"
 !print *,xjjz
+!print *,FAM1
         call unit(tistart, tcountcol,tcol,tval)
 
         end