! Calculate, for each pair of XDS_ASCII.HKL files in XSCALE.INP, the nonisomorphism-angle
! given by cos(phi(i,j)))= cc(i,j)/(cc-star(i)*cc-star(j))
! This is sensible if the XDS_ASCII.HKL files each have significant multiplicity, because
! cc1/2 calculation requires at least 2 observations
!
! version 0.1 Kay Diederichs 18.5.2017
! compile: ifort -check all,noarg_temp_created -traceback -warn all,noexternal -fpe0 -static-intel -init=snan,arrays -sox xds_nonisomorphism.f90 hash.o resolution.o symops.o 
!       or gfortran -O2 -C -Wall -static-libgfortran -static-libgcc xds_nonisomorphism.f90 hash.o symops.o resolution.o 
! with previous compilation like: ifort -c hash.f90 symops.f90 resolution.f90 
!                             or: gfortran -O2 -C -Wall -c hash.f90 symops.f90 resolution.f90 
!
! 27.6.2018 check that unique reflections occur only once
! 28.6.2018 check that data sets are not duplicated

MODULE getresol_mod
      INTEGER, PARAMETER :: h=HUGE(h)
      INTEGER :: ispgr=0
      INTEGER :: ihmin=h,ikmin=h,ilmin=h,ihmax=-h,ikmax=-h,ilmax=-h
      REAL :: resmin=HUGE(resmin),resmax=-1.
      TYPE refdata
        REAL :: int,var=0.
      END TYPE refdata
      TYPE (refdata), ALLOCATABLE :: db(:,:,:,:)
contains
      SUBROUTINE get_nset_reso(unit,name)
      USE resolution_mod
      USE symops_mod
      IMPLICIT NONE
      CHARACTER line*80,name*(*)
      INTEGER :: ier,unit,ih,ik,il,ncheck,jh,jk,jl,ihkl,iop
      REAL :: s2,resol,fsq,sigma

 1000 FORMAT(a)
      
      ncheck=0
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,TRIM(name),' ier=',ier,' line=',line(:LEN_TRIM(line))
          STOP 'read error in xds_file'
        END IF
        IF (line(1:20)=='!SPACE_GROUP_NUMBER=') then
          READ(line(21:),*) ier
          IF (ispgr==0) THEN
            ispgr=ier
            CALL getops(ispgr)
            WRITE(*,'(a)') TRIM(line)
          ELSE IF (ispgr/=ier) THEN
            PRINT*,TRIM(name)
            STOP 'spacegroup mismatch between first and this XDS_ASCII.HKL'
          END IF
        ELSE IF (line(1:21).eq.'!UNIT_CELL_CONSTANTS=') then 
          IF (a==0.) THEN
            READ(line(22:),*) a,b,c,alf,bet,gam
            WRITE(*,'(a)') TRIM(line)
          END IF
        ELSE IF (line(1:9) == '!ITEM_H=1') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:9) == '!ITEM_K=2') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:9) == '!ITEM_L=3') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:12) == '!ITEM_IOBS=4') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:19) == '!ITEM_SIGMA(IOBS)=5') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO
      IF (ncheck <5) THEN
        STOP 'header items ITEM_H, ITEM_K, ITEM_L, ITEM_IOBS, ITEM_SIGMA(IOBS) not found'  
      END IF

! find min and max resolution
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,'ier=',ier
          PRINT*,'line=',TRIM(line)
          STOP 'read error in xds_file'
        END IF
        IF (line(1:12) == '!END_OF_DATA') EXIT

        READ(line,*)ih,ik,il,fsq,sigma
        IF (sigma<0.) CYCLE
        CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,size(op,3)/ncenop)
        ihmin=MIN(jh,ihmin)
        ikmin=MIN(jk,ikmin)
        ilmin=MIN(jl,ilmin)
        ihmax=MAX(jh,ihmax)
        ikmax=MAX(jk,ikmax)
        ilmax=MAX(jl,ilmax)
        CALL getres(ih,ik,il,s2,resol)
        IF (resol<resmin) resmin=resol
        IF (resol>resmax) resmax=resol      
      END DO      
      CLOSE(unit)
      END SUBROUTINE get_nset_reso
END MODULE getresol_mod

SUBROUTINE read_reflection_file(unit,reslim1,reslim2,iset,name)
      USE symops_mod
      USE resolution_mod
      USE getresol_mod
      USE hash_mod
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) :: unit,iset
      REAL, INTENT(IN) :: reslim1,reslim2
      INTEGER, PARAMETER  :: maxref=100000   ! only needed for the hash function
      INTEGER :: ih,ik,il,oldhkl,jh,jk,jl,ihkl,iop,nref,nobs,na,ntot
      REAL :: s2,resol,fsq,sigma,xint,xvar,siga
      CHARACTER :: line*80,name*(*)
      LOGICAL :: last

 1000 FORMAT(A)
      last=.FALSE.
      nobs=0
      nref=0
      ntot=0
      na=0
      xint=0.
      xvar=0.     
      siga=0.
  ! the factor of 2 makes sure that the hash table is never more than half full:
      CALL hash_init(2*maxref)
      
! read header of XDS_ASCII.HKL
      DO
        READ(unit,1000) line
        IF (line(1:14) == '!END_OF_HEADER') EXIT
      END DO
 ! read reflections 
      DO
        READ(unit,1000) line
        IF (line(1:12) == '!END_OF_DATA') last=.TRUE. 
        IF (.NOT.last) THEN
          READ(line,*)ih,ik,il,fsq,sigma
          CALL getres(ih,ik,il,s2,resol)
          IF (resol>reslim1.OR.resol<reslim2) CYCLE
          IF (sigma<0.) CYCLE
          CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,size(op,3)/ncenop)
          nobs=nobs+1
          IF (nobs==1) oldhkl=ihkl  ! make sure oldhkl is sensibly initialized
        END IF
        IF (ihkl /= oldhkl.OR.last) THEN  ! process unique reflection
          nref=nref+1
          IF (nref==maxref+1) THEN
            PRINT*,'>100.000 unique reflections - no more checking for duplicates'
          ELSE IF (nref<=maxref) THEN
            IF (hash_exists(oldhkl)) THEN
              WRITE(*,*) 'error: duplicate unique reflection in ',TRIM(name)
              STOP ' file is not sorted by unique reflections - is it from XDS/XSCALE??'
            ELSE
              CALL hash_insert(oldhkl)
            END IF
          END IF          
! isomorphous signal: store this unique reflection
          IF (na>1) THEN   ! reflection is only stored if CC1/2 can be calculated
            xint=xint/siga
            xvar=xvar/siga
            db(jh,jk,jl,iset)%int=xint
            db(jh,jk,jl,iset)%var=MAX(0.,(xvar-xint**2)/(na-1))
            ntot=ntot+1
          END IF
          IF (last) EXIT  ! the big loop reading the XSCALE file    
          na=0
          xint=0.
          xvar=0.
          siga=0.
          oldhkl=ihkl
        END IF
! store this observation
        xint = xint + fsq/sigma**2
        xvar = xvar + fsq**2/sigma**2
        na = na + 1
        siga = siga + 1/sigma**2
      END DO   ! the big loop reading the XDS_ASCII.HKL file 
      WRITE(*,'(a,i4,3i8)')'iset,nobs,nunique,nunique w/ >1 observations=',iset,nobs,nref,ntot
    
END SUBROUTINE read_reflection_file

PROGRAM xds_nonisomorphism
      USE resolution_mod
      USE getresol_mod
      IMPLICIT NONE
      REAL :: reslim1=999,reslim2=1,tmp
      REAL, ALLOCATABLE :: reslim(:),rmsd(:)
      CHARACTER (LEN=132) :: arg='',name=''
      INTEGER :: unit1,unit2,unit3,ier,nbin=10,mode=2,i,narg,iset,ipass
      LOGICAL :: skiparg=.FALSE.
      INTEGER,DIMENSION(8) :: itime1
 
      narg=COMMAND_ARGUMENT_COUNT()
      WRITE(*,'(a)') 'xds_nonisomorphism KD 2020-07-23. Academic use only.'
      WRITE(*,'(a)') 'Pls cite Diederichs, K. (2017) Acta Cryst D73, 286-293. -h option shows options'  
      CALL DATE_AND_TIME(VALUES=itime1)
!      IF (itime1(1)*10000 + itime1(2)*100 + itime1(3)>20200631) STOP 'expired'
      IF (narg==1) THEN
        CALL GET_COMMAND_ARGUMENT(1, arg)
        IF (arg(1:2)=='-h') THEN
        WRITE(*,'(a)') 'usage: xds_nonisomorphism -dmin <lowres> -dmax <highres> -nbin <nbin> -mode <1 or 2>'
        WRITE(*,'(a)') 'dmin (default 999A), dmax (default 1A) and nbin (default 10) have the usual meanings.'
        WRITE(*,'(a)') 'mode can be 1 (equal volumes of resolution shells) or 2 (increasing volumes; default).'
        WRITE(*,'(a)') 'The program reads XSCALE.INP in ./ to find the paths to the XDS_ASCII.HKL files.'
        WRITE(*,'(a)') 'Spacegroup and unit cell is taken from the _first_ XDS_ASCII.HKL file read.'
        STOP
        END IF
      END IF

      DO i=1,narg
        IF (skiparg) THEN
          skiparg=.FALSE.
          CYCLE
        END IF
        CALL GET_COMMAND_ARGUMENT(i, arg)
        IF      (arg(1:5)=='-DMIN'.OR.arg(1:5)=='-dmin') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) reslim1
          ELSE IF (i<narg) THEN
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) reslim1
            skiparg=.TRUE.
          ENDIF
        ELSE IF (arg(1:5)=='-DMAX'.OR.arg(1:5)=='-dmax') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) reslim1
          ELSE IF (i<narg) THEN
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) reslim2
            skiparg=.TRUE.
          END IF
        ELSE IF (arg(1:5)=='-NBIN'.OR.arg(1:5)=='-nbin') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) nbin
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) nbin 
            skiparg=.TRUE.
          END IF
          WRITE(*,'(a,i0)')' resolution shells: ',nbin
        ELSE IF (arg(1:5)=='-MODE'.OR.arg(1:5)=='-mode') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) mode
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) mode 
            skiparg=.TRUE.
          END IF
          WRITE(*,'(a,i2)')' mode for resolution shells is',mode
        END IF
      END DO 
      IF (reslim1<reslim2) THEN !swap
        tmp=reslim1
        reslim1=reslim2
        reslim2=tmp
      END IF
      WRITE(*,'(a,f8.3)')' low resolution limit :',reslim1
      WRITE(*,'(a,f8.3)')' high resolution limit:',reslim2

      WRITE(*,*)
! read XSCALE.INP and find XDS_ASCII.HKL files 
! - first pass:  get resol and array limits
! - second pass: read data and find nonisomorphism
! - third pass:  write XSCALE.INP.withweights
      OPEN(NEWUNIT=unit1,FILE='XSCALE.INP',STATUS='old',ACTION='READ')
      WRITE(*,*)'reading XSCALE.INP to find XDS_ASCII.HKL-type files'
      DO ipass=1,2
        iset=0
        DO 
          READ(unit1,'(a)',iostat=ier) arg
          IF (ier<0) EXIT  ! end-of-file
          IF (ier>0) STOP 'error reading XSCALE.INP'
          ier=INDEX(arg,'!') 
          IF (ier>0) arg=arg(:ier-1)
          ier=INDEX(arg,'INPUT_FILE=') 
          IF (ier==0) CYCLE
          name=arg(ier+11:)
          DO  ! remove leading blanks
            IF (name(1:1)/=' ')EXIT
            name=name(2:)
          END DO
          OPEN(NEWUNIT=unit2,FILE=name,STATUS='old',ACTION='READ',iostat=ier)
          IF (ier/=0) THEN
            WRITE(*,*) TRIM(name)
            STOP 'error when trying to open'
          END IF
          iset=iset+1
          IF (ipass==1) THEN
            resmin =HUGE(resmin)
            resmax =-1.
            CALL get_nset_reso(unit2,name) ! this sets resmax,resmin
            WRITE(*,'(a,i3,2f8.3,1x,a)') 'iset, dmin, dmax, name:',iset,resmax,resmin,TRIM(name)
            reslim1=MIN(reslim1,resmax)   ! overall max numerical resolution
            reslim2=MAX(reslim2,resmin)   ! overall min numerical resolution
          ELSE IF (ipass==2) THEN
            CALL read_reflection_file(unit2,reslim1,reslim2,iset,name)
          END IF
          CLOSE(unit2)  
        END DO
        REWIND(unit1)
        IF (ipass==1) THEN
          CALL setbin(nbin,reslim1,reslim2,mode)
          ALLOCATE(reslim(nbin))
          CALL getlim(reslim)
          ALLOCATE(db(ihmin:ihmax,ikmin:ikmax,ilmin:ilmax,iset),rmsd(iset))
        ELSE IF (ipass==2) THEN
          WRITE(*,*)
          WRITE(*,'(a,2f7.3)')'Lowest and highest resolution used:',reslim1,reslim2
          WRITE(*,'(i3,a)') nbin,' resolution shells:'
          WRITE(*,'(40f7.3)') reslim      
          CALL print_results(iset,nbin,rmsd)   
        END IF
      END DO ! ipass
! pass 3:
      IF (iset>2) THEN   ! don't write XSCALE.INP.rename_me if just 2 data sets
        OPEN(NEWUNIT=unit3,FILE='XSCALE.INP.rename_me',STATUS='unknown',ACTION='WRITE')
        iset=0
        DO 
          READ(unit1,'(a)',iostat=ier) arg
          IF (ier<0) EXIT  ! end-of-file
! is it a line with INPUT_FILE= ?
          ier=INDEX(arg,'!') 
          name=arg
          IF (ier>0) name=arg(:ier-1)
          ier=INDEX(name,'INPUT_FILE=') 
          WRITE(unit3,'(a)') TRIM(arg)      ! copy line
          IF (ier>0) THEN
            iset=iset+1
            WRITE(unit3,'(a,f10.4)') '!RMSD=',rmsd(iset)
          END IF
        END DO
        STOP ' wrote XSCALE.INP.rename_me with additional RMSD= lines'
      END IF
END PROGRAM xds_nonisomorphism

SUBROUTINE print_results(nset,nbin,rmsd)
      USE getresol_mod
      USE resolution_mod
      IMPLICIT NONE
      REAL :: sumxint(nbin),sumyint(nbin),sumxintq(nbin),sumyintq(nbin), &
              sumxvar(nbin),sumyvar(nbin),sumxy(nbin),sumres(nbin), &
              cc12x,cc12y,t1,t2,ccxy,s2,resol,sigma,tau,diffatom,sigcc12x,sigcc12y,sigccxy
      INTEGER :: nset,nbin,nmatch(nbin),iset1,iset2,ih,ik,il,ibin,n,ibest,iworst,ipass,ipassmax
      REAL rmsd(nset),xbest,xworst
      REAL, ALLOCATABLE :: diffmat(:,:)
      INTEGER, ALLOCATABLE :: ndiffmat(:,:)
      
      ALLOCATE(diffmat(nset,nset),ndiffmat(nset,nset))
      ipassmax=1  ! 2
      DO ipass=1,ipassmax 
      diffmat=0
      ndiffmat=0
      IF (ipass==2) THEN
        PRINT*,'-------------------------------------------------------------------'
        PRINT*,'only reflection pairs with I/sigI>3 of each member of the pair:'
      END IF      
      DO iset1=1,nset-1
        DO iset2=iset1+1,nset
          PRINT*,''
          PRINT*,'iset1, iset2=',iset1,iset2
          WRITE(*,'(a)')'resol_shell nmatch   CC*_1     CC*_2     CC(1,2)  cos(phi) angle(deg) RMSD_coord' !  b'
          sumxint=0
          sumyint=0
          sumxintq=0
          sumyintq=0
          sumxvar=0
          sumyvar=0
          sumxy=0
          nmatch=0
          sumres=0
          DO il=ilmin,ilmax
          DO ik=ikmin,ikmax
          DO ih=ihmin,ihmax
            IF (db(ih,ik,il,iset1)%var * db(ih,ik,il,iset2)%var == 0) CYCLE
            IF (ipass==2) THEN
              IF (db(ih,ik,il,iset1)%int < 3*SQRT(db(ih,ik,il,iset1)%var)) CYCLE  ! weak
              IF (db(ih,ik,il,iset2)%int < 3*SQRT(db(ih,ik,il,iset2)%var)) CYCLE  ! weak
            END IF
! find resolution bin
            CALL getres(ih,ik,il,s2,resol)
            CALL getbin(s2,ibin)
! increment sums
            sumxint(ibin)  = sumxint(ibin)  + db(ih,ik,il,iset1)%int
            sumxintq(ibin) = sumxintq(ibin) + db(ih,ik,il,iset1)%int**2
            sumxvar(ibin)  = sumxvar(ibin)  + db(ih,ik,il,iset1)%var
            sumyint(ibin)  = sumyint(ibin)  + db(ih,ik,il,iset2)%int
            sumyintq(ibin) = sumyintq(ibin) + db(ih,ik,il,iset2)%int**2
            sumyvar(ibin)  = sumyvar(ibin)  + db(ih,ik,il,iset2)%var
            nmatch(ibin)   = nmatch(ibin)   + 1
            sumxy(ibin)=sumxy(ibin) + db(ih,ik,il,iset1)%int * db(ih,ik,il,iset2)%int
            sumres(ibin)   = sumres(ibin)   + resol
          END DO
          END DO
          END DO
          DO ibin=1,nbin
            IF (nmatch(ibin)<3) CYCLE
! CC1/2 and CC-star of iset1:
            sigma = (sumxintq(ibin) - sumxint(ibin)**2/nmatch(ibin)) / (nmatch(ibin)-1)
            IF (sigma<=0.) CYCLE   ! this also takes care of t1 (below)
            tau   = sumxvar(ibin) / nmatch(ibin)
            cc12x = (sigma-tau)/(sigma+tau)
            IF (cc12x<=0) CYCLE    
! sigcc12x = upper limit of cc12x
            sigcc12x = MIN(1.,cc12x + 2*SQRT((1.-cc12x**2)/(nmatch(ibin)-2)))
            sigcc12x = SQRT(2*sigcc12x/(1+sigcc12x))

            cc12x = SQRT(2*cc12x/(1+cc12x))
! CC1/2 and CC-star of iset2:
            sigma = (sumyintq(ibin) - sumyint(ibin)**2/nmatch(ibin)) / (nmatch(ibin)-1)
            IF (sigma<=0.) CYCLE   ! this also takes care of t2 (below)
            tau   = sumyvar(ibin) / nmatch(ibin)
            cc12y = (sigma-tau)/(sigma+tau)
            IF (cc12y<=0) CYCLE
! sigcc12y = upper limit of cc12y
            sigcc12y = MIN(1.,cc12y + 2*SQRT((1.-cc12y**2)/(nmatch(ibin)-2)))
            sigcc12y = SQRT(2*sigcc12y/(1+sigcc12y))
            
            cc12y = SQRT(2*cc12y/(1+cc12y))
! CC of iset1 vs iset2:
            t1 = sumxintq(ibin) - sumxint(ibin)**2/nmatch(ibin)
            t2 = sumyintq(ibin) - sumyint(ibin)**2/nmatch(ibin)
            ccxy = (sumxy(ibin) - sumxint(ibin)*sumyint(ibin)/nmatch(ibin))/SQRT(t1*t2)
            t1=ccxy/(cc12x*cc12y)   ! cos(deltaphi)   or can be interpreted as random_noise-corrected correlation
            IF (ABS(t1)>1.) THEN
!standard error of cc is SQRT(1.-cc**2)/(n-2))
!if (cc12x**2>1 .OR. cc12y**2>1 .OR. ccxy**2>1) print*,'cc12x,cc12y,ccxy,nmatch(ibin)=',cc12x,cc12y,ccxy,nmatch(ibin)
     
                sigccxy =SQRT(MAX(0.,1.-ccxy**2) /(nmatch(ibin)-2))    ! sometimes (often if n=3), 1-ccxy is -0 !
                IF (sigcc12x*sigcc12y < ccxy) THEN
                PRINT '(a,2i5,a)','CC(1,2) is significantly higher than CC*_1 * CC*_2 . Are',iset1,iset2,' duplicates ??'
              END IF
              t1=t1/ABS(t1)
            END IF
            sumres(ibin)=sumres(ibin)/nmatch(ibin)
            t2=ACOS(t1)*180/3.14159
!            diffatom=t2/180*sumres(ibin)      !  interpretiert t1 als cos(phasedifference) und interpretiert letzteres als Koordinaten-RMS
! 
! was ist das Verhaeltnis von rms-coord difference zur correlation?
! https://www.phenix-online.org/documentation/reference/phaser.html   gives defaults Bsol=300 and Fsol=0.95 for eqn 1a
! http://www.pnas.org/content/114/14/3637.full.pdf equation 1
! actually a coordinate deviation has a similar effect as a B-factor
            diffatom=sumres(ibin)*SQRT(-3./(2.*3.14159**2)*log(MAX(0.01,t1)))  ! eqn 1b with fp=1
            IF (t1<0.01) diffatom=10.
            PRINT '(i10,i8,7f10.4)',ibin,nmatch(ibin),cc12x,cc12y,ccxy,t1,t2,diffatom   !,1-t1**2           
            IF (ipass==ipassmax.AND.diffatom<9.9) THEN           
              diffmat(iset1,iset2)=diffmat(iset1,iset2)+diffatom
              ndiffmat(iset1,iset2)=ndiffmat(iset1,iset2)+1
            END IF
          END DO
        END DO
      END DO
      PRINT*,' '
      END DO ! ipass
!
      DO iset1=1,nset-1
        DO iset2=iset1+1,nset
          IF (ndiffmat(iset1,iset2)>0) diffmat(iset1,iset2)=diffmat(iset1,iset2)/ndiffmat(iset1,iset2)
          diffmat(iset2,iset1)=diffmat(iset1,iset2)     ! diffmat is symmetric
          ndiffmat(iset2,iset1)=ndiffmat(iset1,iset2)
        END DO
      END DO
      PRINT*,''
      PRINT*,'using average RMSD values (excluding unreasonable table entries):'
      ibest=0
      iworst=0
      xbest=HUGE(xbest)
      xworst=-HUGE(xworst)
      DO iset1=1,nset
        n=0
        rmsd(iset1)=0
        DO iset2=1,nset
          IF (ndiffmat(iset2,iset1)==0) CYCLE
          n=n+1
          rmsd(iset1)=rmsd(iset1)+diffmat(iset2,iset1)
        END DO
        IF (n>0) rmsd(iset1)=rmsd(iset1)/n
        PRINT*,'dataset #, mean RMSD to all other datasets:',iset1,rmsd(iset1)
        IF (rmsd(iset1)<xbest) THEN
          xbest=rmsd(iset1)
          ibest=iset1
        END IF
        IF (rmsd(iset1)>xworst) THEN
          xworst=rmsd(iset1)
          iworst=iset1
        END IF
      END DO
      IF (nset>2) THEN
        PRINT*,'central dataset (most isomorphous) is number',ibest
        PRINT*,'most distant dataset (least isom.) is number',iworst
        PRINT*,'RMSD= lines in XSCALE.INP.rename_me will be specified w.r.t. to central dataset'
      END IF
      DO iset1=1,nset
        diffmat(iset1,iset1)=0        ! diagonal of diffmat has zeroes
        ndiffmat(iset1,iset1)=1 
        IF (ndiffmat(iset1,ibest)>0) rmsd(iset1)=diffmat(iset1,ibest)
      END DO
! at this point, sets that have no relative RMSD to ibest will just have their average RMSD to all others
      IF (ANY(ndiffmat==0)) THEN
        PRINT*,'no calculation of 3D coordinates because some relative RMSDs are undefined'
        RETURN
      END IF
      IF (nset<3) RETURN
BLOCK   
REAL, ALLOCATABLE :: result(:,:)
      ALLOCATE (result(3,nset))
      CALL distance2coord(diffmat,nset,result,3)
      PRINT*,'coordinates in 3D that best reproduce the angles as distances:'
      DO iset1=1,nset
        PRINT*,result(:,iset1)
      END DO
      CALL wpdb(100*result,nset,3,'noniso.pdb','A')
      PRINT*,'noniso.pdb=representation of data set arrangement in 3D (coords*100)'
END BLOCK

END SUBROUTINE print_results

SUBROUTINE wpdb(xyz,n,idim,fn,segid)
      implicit none
      integer i,idim,n,unit
      real xyz(idim,n)
      character*(*) fn
      character segid*1
! 102   format('ATOM',i7,'  CA  SER ',a1,i4,t31,3f8.3,t60,f6.2,t55,f6.2)
 102  format('HETATM',i5,'  O   HOH ',a1,i4,t31,3f8.3,t60,f6.2,t55,f6.2)
 
      open(newunit=unit,file=fn)
      write(unit,'(a)') 'CRYST1  100.000  100.000  100.000  90.00  90.00  90.00 P 1'
      do i=1,n
        IF (idim==1) THEN
          write(unit,102) i,segid,i,xyz(1,i),0.,0.,100.,1.
        ELSE IF (idim==2) THEN
          write(unit,102) i,segid,i,xyz(1,i),xyz(2,i),0.,100.,1.
        ELSE IF (idim==3) THEN
          write(unit,102) i,segid,i,xyz(1,i),xyz(2,i),xyz(3,i),100.,1.
        ELSE IF (idim==4) THEN
          write(unit,102) i,segid,i,xyz(1,i),xyz(2,i),xyz(3,i),100+xyz(4,i),1.
        ELSE IF (idim==5) THEN
          write(unit,102) i,segid,i,xyz(1,i),xyz(2,i),xyz(3,i),xyz(4,i),(100+xyz(5,i))/100.
        END IF
      end do
      print*,'wrote ',fn(:LEN_TRIM(fn))
      close(unit)
END SUBROUTINE wpdb


SUBROUTINE distance2coord(z,n,result,idim)
! determine result(idim,n) such that Euclidean distances between 
! result(:,i) and result(:,j) matches z(i,j)
! This implements the algebraic solution in article
! http://de.wikipedia.org/wiki/Multidimensionale_Skalierung#Algorithmus
! but it (deviating from the article) does not normalize the eigenvectors.
! z is assumed to be symmetric
      INTEGER, INTENT(in) :: n,idim
      REAL :: z(n,n)
      REAL, INTENT(out) :: result(idim,n)
! local variables
      REAL, ALLOCATABLE :: a(:,:),w(:)
      REAL v
      INTEGER i,j,it_num,num_rot
      
      IF (idim>n) STOP 'dim>n'    !  das war >= ! warum? KD 16.10.17
      ALLOCATE(a(n,n),w(n))
! step 1 in Wikipedia article:
      DO i=1,n
        DO j=1,n
          z(j,i)=-0.5d0*z(j,i)**2
        END DO
      END DO
! step 2:
      w=SUM(z,dim=2)/n               ! a(j,.)          in the article  
      DO i=1,n
        v=SUM(z(:,i))/n -SUM(w)/n    ! a(.,i) + a(.,.) in the article
        DO j=1,n
          a(j,i)=z(j,i)-w(j)-v
        END DO
      END DO
!      CALL arpack_syev(n,a,w,z,idim,info)! w=eigenvalues z=eigenvectors nev=# eigenvalues
CALL jacobi_eigenvalue ( n, a, 100, z, w, it_num, num_rot )
print*,'Jacobi it_num,num_rot:',it_num,num_rot
print*,'Eigenvalues:',w
!      IF (info/=0) PRINT*,'ERROR: info=',info

      DO j=1,idim   ! highest eigenvalue is x, second highest y, ...
        DO i=1,n
!          result(j,i)= SQRT(w(idim+1-j))*z(i,idim+1-j) ! scale Eigenvectors by sqrt(Eigenvalue) ; code for syev
          result(j,i)= SQRT(MAX(0.,w(n+1-j)))*z(i,n+1-j) ! scale Eigenvectors by sqrt(Eigenvalue)
        END DO
      END DO

END SUBROUTINE distance2coord

subroutine jacobi_eigenvalue ( n, a, it_max, v, d, it_num, rot_num )
! https://people.sc.fsu.edu/~jburkardt/f_src/jacobi_eigenvalue/jacobi_eigenvalue.f90

!*****************************************************************************80
!
!! JACOBI_EIGENVALUE carries out the Jacobi eigenvalue iteration.
!
!  Discussion:
!
!    This function computes the eigenvalues and eigenvectors of a
!    real symmetric matrix, using Rutishauser's modfications of the classical
!    Jacobi rotation method with threshold pivoting. 
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    17 September 2013
!
!  Author:
!
!    FORTRAN90 version by John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the order of the matrix.
!
!    Input, real ( kind = 8 ) A(N,N), the matrix, which must be square, real,
!    and symmetric.
!
!    Input, integer ( kind = 4 ) IT_MAX, the maximum number of iterations.
!
!    Output, real ( kind = 8 ) V(N,N), the matrix of eigenvectors.
!
!    Output, real ( kind = 8 ) D(N), the eigenvalues, in descending order.
!
!    Output, integer ( kind = 4 ) IT_NUM, the total number of iterations.
!
!    Output, integer ( kind = 4 ) ROT_NUM, the total number of rotations.
!
  implicit none

  integer ( kind = 4 ),INTENT(IN) :: n

  real ( kind = 4 ),INTENT(INOUT) :: a(n,n)
  real ( kind = 8 ) bw(n)
  real ( kind = 8 ) c
  real ( kind = 4 ),INTENT(OUT) :: d(n)
  real ( kind = 8 ) g
  real ( kind = 8 ) gapq
  real ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ),INTENT(IN) :: it_max
  integer ( kind = 4 ),INTENT(OUT) :: it_num
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) p
  integer ( kind = 4 ) q
  integer ( kind = 4 ),INTENT(OUT) :: rot_num
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) tau
  real ( kind = 8 ) term
  real ( kind = 8 ) termp
  real ( kind = 8 ) termq
  real ( kind = 8 ) theta
  real ( kind = 8 ) thresh
  real ( kind = 4 ),INTENT(OUT) :: v(n,n)
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) zw(n)

  do j = 1, n
    do i = 1, n
      v(i,j) = 0.0D+00
    end do
    v(j,j) = 1.0D+00
  end do

  do i = 1, n
    d(i) = a(i,i)
  end do

  bw(1:n) = d(1:n)
  zw(1:n) = 0.0D+00
  it_num = 0
  rot_num = 0

  do while ( it_num < it_max )

    it_num = it_num + 1
!
!  The convergence threshold is based on the size of the elements in
!  the strict upper triangle of the matrix.
!
    thresh = 0.0D+00
    do j = 1, n
      do i = 1, j - 1
        thresh = thresh + a(i,j) ** 2
      end do
    end do

    thresh = sqrt ( thresh ) / real ( 4 * n, kind = 8 )

    if ( thresh == 0.0D+00 ) then
      exit 
    end if

    do p = 1, n
      do q = p + 1, n

        gapq = 10.0D+00 * abs ( a(p,q) )
        termp = gapq + abs ( d(p) )
        termq = gapq + abs ( d(q) )
!
!  Annihilate tiny offdiagonal elements.
!
        if ( 4 < it_num .and. &
             termp == abs ( d(p) ) .and. &
             termq == abs ( d(q) ) ) then

          a(p,q) = 0.0D+00
!
!  Otherwise, apply a rotation.
!
        else if ( thresh <= abs ( a(p,q) ) ) then

          h = d(q) - d(p)
          term = abs ( h ) + gapq

          if ( term == abs ( h ) ) then
            t = a(p,q) / h
          else
            theta = 0.5D+00 * h / a(p,q)
            t = 1.0D+00 / ( abs ( theta ) + sqrt ( 1.0D+00 + theta * theta ) )
            if ( theta < 0.0D+00 ) then 
              t = - t
            end if
          end if

          c = 1.0D+00 / sqrt ( 1.0D+00 + t * t )
          s = t * c
          tau = s / ( 1.0D+00 + c )
          h = t * a(p,q)
!
!  Accumulate corrections to diagonal elements.
!
          zw(p) = zw(p) - h                  
          zw(q) = zw(q) + h
          d(p) = d(p) - h
          d(q) = d(q) + h

          a(p,q) = 0.0D+00
!
!  Rotate, using information from the upper triangle of A only.
!
          do j = 1, p - 1
            g = a(j,p)
            h = a(j,q)
            a(j,p) = g - s * ( h + g * tau )
            a(j,q) = h + s * ( g - h * tau )
          end do

          do j = p + 1, q - 1
            g = a(p,j)
            h = a(j,q)
            a(p,j) = g - s * ( h + g * tau )
            a(j,q) = h + s * ( g - h * tau )
          end do

!DIR$ SIMD
          do j = q + 1, n
            g = a(p,j)
            h = a(q,j)
            a(p,j) = g - s * ( h + g * tau )
            a(q,j) = h + s * ( g - h * tau )
          end do
!
!  Accumulate information in the eigenvector matrix.
!
          do j = 1, n
            g = v(j,p)
            h = v(j,q)
            v(j,p) = g - s * ( h + g * tau )
            v(j,q) = h + s * ( g - h * tau )
          end do

          rot_num = rot_num + 1

        end if

      end do
    end do

    bw(1:n) = bw(1:n) + zw(1:n)
    d(1:n) = bw(1:n)
    zw(1:n) = 0.0D+00

  end do
!
!  Restore upper triangle of input matrix.
!
  do j = 1, n
    do i = 1, j - 1
      a(i,j) = a(j,i)
    end do
  end do
!
!  Ascending sort the eigenvalues and eigenvectors.
!
  do k = 1, n - 1

    m = k

    do l = k + 1, n
      if ( d(l) < d(m) ) then
        m = l
      end if
    end do

    if ( m /= k ) then

      t    = d(m)
      d(m) = d(k)
      d(k) = t

      w(1:n)   = v(1:n,m)
      v(1:n,m) = v(1:n,k)
      v(1:n,k) = w(1:n)

    end if

  end do

  return
end

