Hotfix for Qfep analysis tool (#1)

* Hot fix for Qfep analysis tool

Was missing fix from old version to use of different number of
off-diagonal elements. Added now

* Fix last issues in qfep

Fixed the remaining issues from previous fix.

* I'm obviously blind

see title

src/qfep.f90

@@ -34,7 +34,7 @@ program qfep
 	integer		:: mxpts=200,mxbin=10,mxstates=4
 	character(80)      ::filnam, line
 	integer           ::i,j,ifile,ipt,istate,ibin,nfiles,nstates,ERR, &
-	                      nskip,nbins,nmin,idum,noffd,nnoffd,offel
+	                      nskip,nbins,nmin,idum,noffd,nnoffd,offel,num_offd
 
 	type(OFFDIAG_SAVE),allocatable		:: offd(:)
 
@@ -195,7 +195,7 @@ program qfep
 		read (*,*) alfa(istate)
 		write (*,6) istate,alfa(istate)
 	end do
-6	format('# Alpha for state ',i2,'              =',f6.2)
+6	format('# Alpha for state ',i2,'              =',f9.2)
 7	format('--> Give alpha for state ',i2,':')
 
 	scale_Hij=0.0
@@ -204,6 +204,8 @@ program qfep
 		read (*,*) scale_Hij
 		write (*,8) scale_Hij
 8	format('# Scale factor for Hij            =',f6.2)
+! variable to store the actual offdiagonals used later
+                num_offd = noffd
 	else
 		call prompt ('--> Number of off-diagonal elements:')
 		read (*,*) nnoffd
@@ -214,6 +216,8 @@ program qfep
 			read (*,*) i, j, A(i,j), mu(i,j), eta(i,j), rxy0(i,j)
 			write(*,12) i, j, A(i,j), mu(i,j), eta(i,j), rxy0(i,j)
 		end do
+! variable to store the actual offdiagonals used later
+                num_offd = nnoffd
 9	format('# Number of off-diagonal elements =',i6)
 11	format('#   i   j   A(i,j)  mu(i,j) eta(i,j) rxy0(i,j)')
 12	format('#',2i4,1x,f9.2,1x,f9.2,1x,f9.2,1x,f9.2)
@@ -344,7 +348,7 @@ program qfep
 		write(*,'(a,a)') '# Energy files from Qdyn, version: ',fileheader%version
 		!read 1st record to get lambdas
 		if(is_old_file) rewind(f)
-		idum = get_ene(f, EQ(:), offd, nstates,nnoffd,fileheader%arrays)
+		idum = get_ene(f, EQ(:), offd, nstates,num_offd,fileheader%arrays)
 		if(idum /= 0) then
 			!we have a problem
 			write(*,'(a,a)') 'ERROR: Unexpected end-of-file in first record of ', &
@@ -395,7 +399,7 @@ program qfep
                                 fileheader%gcnum(1:fileheader%arrays),fileheader%version
 
 		end if
-		do while(get_ene(f, EQ(:), offd, nstates, nnoffd,fileheader%arrays) == 0) !keep reading till EOF
+		do while(get_ene(f, EQ(:), offd, nstates, num_offd,fileheader%arrays) == 0) !keep reading till EOF
 			ipt = ipt + 1
 			if (ipt .eq. (mxpts-1)) then
 !ops, we need larger arrays to keep track of the temporary points
@@ -446,7 +450,7 @@ program qfep
 				do j=1,fileheader%arrays
 				EQ(i)%total(j)=EQ(i)%total(j)+alfa(i)
 				FEPtmp%v(j,i, ipt) = EQ(i)%total(j)
-				if (nnoffd .ne. 0) then
+				if (num_offd .ne. 0) then
 					FEPtmp%r(j,i,ipt) = offd(i)%rkl
 			        end if
 				end do
@@ -488,10 +492,10 @@ program qfep
 			do j=1,fileheader%arrays
 				if (nstates==2) then
 					FEPtmp%vg(j,ipt)=0.5_prec*(EQ(1)%total(j)+EQ(2)%total(j))-  &
-						0.5_prec*q_sqrt( (EQ(1)%total(j)-EQ(2)%total(j))**2 + 4.*Hij(1,2,j)**2 )
-					if(nnoffd > 0) then
-						FEPtmp%c1(j,ipt)=1./(1.+((FEPtmp%vg(j,ipt)-EQ(1)%total(j))/Hij(1,2,j))**2)
-						FEPtmp%c2(j,ipt)=1-FEPtmp%c1(j,ipt)
+						0.5_prec*q_sqrt( (EQ(1)%total(j)-EQ(2)%total(j))**2 + 4.0_prec*Hij(1,2,j)**2 )
+					if(num_offd > 0) then
+						FEPtmp%c1(j,ipt)=one/(one+((FEPtmp%vg(j,ipt)-EQ(1)%total(j))/Hij(1,2,j))**2)
+						FEPtmp%c2(j,ipt)=one-FEPtmp%c1(j,ipt)
 					end if
 				else
 					call tred2(Hij(:,:,j),nstates,d,e)
@@ -528,7 +532,7 @@ program qfep
 		FEP(ifile)%gap(:,:) = FEPtmp%gap(1:fileheader%arrays,1:FEP(ifile)%npts)
 		FEP(ifile)%c1(:,:) = FEPtmp%c1(1:fileheader%arrays,1:FEP(ifile)%npts)
 		FEP(ifile)%c2(:,:) = FEPtmp%c2(1:fileheader%arrays,1:FEP(ifile)%npts)
-		if(nnoffd > 0) then
+		if(num_offd > 0) then
 			allocate(FEP(ifile)%r(fileheader%arrays,nstates, FEP(ifile)%npts), STAT=ERR)
 			if(ERR /= 0) then
 				stop 'Qfep5 terminated abnormally: Out of memory when allocating arrays. 1'
@@ -741,7 +745,7 @@ program qfep
 	         		avc1(ibin)=avc1(ibin)+FEP(ifile)%c1(j,ipt)
                 		avc2(ibin)=avc2(ibin)+FEP(ifile)%c2(j,ipt)
 				!Only gives first r_xy distance
-				if(nnoffd > 0)	avr(ibin)=avr(ibin)+FEP(ifile)%r(j,1,ipt)
+				if(num_offd > 0) avr(ibin)=avr(ibin)+FEP(ifile)%r(j,1,ipt)
 				nbinpts(ibin)=nbinpts(ibin)+1
 			end do          !ipt
 			do ibin=1,nbins
@@ -827,7 +831,7 @@ program qfep
 	do ifile=1,nfiles
 		deallocate(FEP(ifile)%v, FEP(ifile)%vg, FEP(ifile)%gap, &
 		FEP(ifile)%c1, FEP(ifile)%c2,FEP(ifile)%lambda)
-		if(nnoffd > 0) deallocate(FEP(ifile)%r)
+		if(num_offd > 0) deallocate(FEP(ifile)%r)
 	end do
 	deallocate(FEP)
 	deallocate(FEPtmp%v,FEPtmp%vg,FEPtmp%gap,FEPtmp%c1,FEPtmp%c2,FEPtmp%lambda,FEPtmp%r)