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