From 3278aabfeb7bf73cc8f2788d2d1dfe067976535b Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Sat, 27 Jan 2018 12:41:48 +0100 Subject: [PATCH] updating code to the current local version, it might not compile atm. Many new features added: 1. getting S2 values 2. possibility of setting position of hole 3. possibility of setting Sbox 4. three Sbox definitions at once 5. Doing only FAM1 or the full set of states 6. efficiency improvements --- README.md | 1 + src/adr.irp.f | 16 +- src/adrnew.irp.f | 57 +++ src/analyse.irp.f | 832 ++++++++++++++++++++++++++++++++++++++++ src/conv.irp.f | 14 +- src/desort.irp.f | 10 +- src/ex1.c | 448 ++++++++++------------ src/extra_diag.irp.f | 2 +- src/get_dmat.c | 145 +++++++ src/get_dmat.h | 9 + src/get_s2.c | 683 +++++++++++++++++++++++++++++++++ src/get_s2.h | 50 +++ src/get_s2_cyclic.c | 784 +++++++++++++++++++++++++++++++++++++ src/get_s2_mov.c | 683 +++++++++++++++++++++++++++++++++ src/getdet.irp.f | 57 +-- src/natom.irp.f | 4 +- src/providdet.irp.f | 85 ++++ src/provide_clust.irp.f | 1 + src/read2.c | 79 +++- src/read2.h | 21 +- src/searchdet.irp.f | 110 +++--- src/sort.irp.f | 10 +- src/unit_l1.irp.f | 4 + 23 files changed, 3721 insertions(+), 384 deletions(-) create mode 100644 src/adrnew.irp.f create mode 100644 src/analyse.irp.f create mode 100644 src/get_dmat.c create mode 100644 src/get_dmat.h create mode 100644 src/get_s2.c create mode 100644 src/get_s2.h create mode 100644 src/get_s2_cyclic.c create mode 100644 src/get_s2_mov.c create mode 100644 src/providdet.irp.f diff --git a/README.md b/README.md index 6170e0e..48e92e9 100644 --- a/README.md +++ b/README.md @@ -47,6 +47,7 @@ _Using DEHam_ 2 # 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, J 2 for K and 3 for none) diff --git a/src/adr.irp.f b/src/adr.irp.f index 76db98d..2a9cd5c 100644 --- a/src/adr.irp.f +++ b/src/adr.irp.f @@ -8,24 +8,24 @@ subroutine adr(ideter,add) 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(kind=selected_int_kind(16))::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/adrnew.irp.f b/src/adrnew.irp.f new file mode 100644 index 0000000..91ef8c0 --- /dev/null +++ b/src/adrnew.irp.f @@ -0,0 +1,57 @@ +subroutine adrfull() + implicit none + 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(kind=selected_int_kind(16))::add + integer(kind=selected_int_kind(16))::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..c621bc2 --- /dev/null +++ b/src/analyse.irp.f @@ -0,0 +1,832 @@ + SUBROUTINE ANALYSE(vect, dimvect, startvect, endvect, xymat2, norm2) +! INCLUDE "nbtots.prm" + IMPLICIT NONE + INTEGER dimvect, nbtots, startvect, endvect + REAL*8,dimension(dimvect)::vect + INTEGER (kind=selected_int_kind(16))::add,kvect + INTEGER (kind=selected_int_kind(16))::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..e01b839 100644 --- a/src/conv.irp.f +++ b/src/conv.irp.f @@ -1,20 +1,20 @@ -subroutine conv(ideter,det,deth) +subroutine conv(ideter,deti,dethi) 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(kind=selected_int_kind(16)),INTENT(INOUT)::deti + integer(kind=selected_int_kind(16)),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..643a7ad 100644 --- a/src/desort.irp.f +++ b/src/desort.irp.f @@ -1,14 +1,14 @@ subroutine desort() implicit none integer::i,j,ord,ordh - integer(kind=selected_int_kind(16))::add,addh,det,deth,addt + integer(kind=selected_int_kind(16))::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 +17,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/ex1.c b/src/ex1.c index cae2605..b37e99c 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -1,12 +1,17 @@ -#include #include #include -#include "stimsyr.h" + #include "read2.h" +#include "stimsyr.h" +#include "get_s2.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 */ @@ -14,7 +19,13 @@ int main(int argc,char **argv) EPSType type; PetscReal error,tol,re,im; 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; PetscScalar kr,ki,value[700]; Vec xr,xi; PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol; @@ -25,12 +36,7 @@ int main(int argc,char **argv) //PetscScalar eigr; //PetscScalar eigi; int mpiid; - int natomax=700, iaa2, iaa; - int ideter[natomax]; - int ideter2[natomax]; - PetscScalar densmat[4][4]={0.0}; - PetscScalar densmatfin[4][4]={0.0}; - PetscScalar trace; + int natomax=700; char const* const fileName = argv[1]; FILE* file = fopen(fileName, "r"); @@ -39,52 +45,83 @@ int main(int argc,char **argv) Data_new(file, &getdata); nlocal = getdata.n/getdata.npar; -//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); +//printf("n=%ld\t nsites=%d\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.natom, getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz); - PetscScalar valxr[nlocal]; - PetscInt indxr[nlocal]; + PetscScalar *valxr; + PetscInt indxr[nlocal]; //Vec Vr,Vi; 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 ii; long int tcountcol2,tcol[700],tcountcol[getdata.nnz]; double val[700]; - double xmat=0.0; 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; - int kko,kok,kkio; - Vec xiaa2; /* initial vector, destination vector */ - VecScatter scatter; /* scatter context */ - IS from, to; /* index sets that define the scatter */ - PetscScalar *values; - int idx_from[] = {0}, idx_to[] = {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)); + // fn = x^2+x^3 +//PetscInt range[] ={165094,310638,438886,551954,651820,740324,819168,889916,953994,1012690,1067154,1118398,1167296,1214584,1260860,1306584,1352078,1517172,1662716,1790964,1904032,2003898,2092402,2171246,2241994,2306072,2364768,2419232,2470476,2519374,2566662,2612938,2658662,2704156,2869250,3014794,3143042,3256110,3355976,3444480,3523324,3594072,3658150,3716846,3771310,3822554,3871452,3918740,3965016,4010740,4056234,4221328,4366872,4495120,4608188,4708054,4796558,4875402,4946150,5010228,5068924,5123388,5174632,5223530,5270818,5317094,5362818,5408312,5573406,5718950,5847198,5960266,6060132,6148636,6227480,6298228,6362306,6421002,6475466,6526710,6575608,6622896,6669172,6714896,6760390,6925484,7071028,7199276,7312344,7412210,7500714,7579558,7650306,7714384,7773080,7827544,7878788,7927686,7974974,8021250,8066974,8112468,8277562,8423106,8551354,8664422,8764288,8852792,8931636,9002384,9066462,9125158,9179622,9230866,9279764,9327052,9373328,9419052,9464546,9629640,9775184,9903432,10016500,10116366,10204870,10283714,10354462,10418540,10477236,10531700,10582944,10631842,10679130,10725406,10771130,10816624,10981718,11127262,11255510,11368578,11468444,11556948,11635792,11706540,11770618,11829314,11883778,11935022,11983920,12031208,12077484,12123208,12168702,12333796,12479340,12607588,12720656,12820522,12909026,12987870,13058618,13122696,13181392,13235856,13287100,13335998,13383286,13429562,13475286,13520780,13685874,13831418,13959666,14072734,14172600,14261104,14339948,14410696,14474774,14533470,14587934,14639178,14688076,14735364,14781640,14827364,14872858,15037952,15183496,15311744,15424812,15524678,15613182,15692026,15762774,15826852,15885548,15940012,15991256,16040154,16087442,16133718,16179442,16224936}; + 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 = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,2.0*getdata.natom,NULL,2.0*getdata.natom,NULL,&A);CHKERRQ(ierr); + ierr = MatMPIAIJSetPreallocation(A,getdata.natom,NULL,getdata.natom,NULL);CHKERRQ(ierr); //ierr = MatSetFromOptions(A);CHKERRQ(ierr); //ierr = MatSetUp(A);CHKERRQ(ierr); MPI_Comm_rank(MPI_COMM_WORLD,&mpiid); ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); ierr = PetscTime(&tt1);CHKERRQ(ierr); + ierr = PetscPrintf(PETSC_COMM_WORLD," start: %d end: %d\n",Istart, Iend);CHKERRQ(ierr); +// Iend = range[mpiid]; +// if(mpiid==0){ +// Istart = 0; +// } +// else{ +// Istart = range[mpiid-1]; +// } for (i=Istart; i0) { + if (0) { PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer); PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB); PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU); @@ -209,218 +239,125 @@ int main(int argc,char **argv) 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; - for (ii=Istart; ii +#include +#include +#include +#include +#include +#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=700; + 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=700; + 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..29539e5 --- /dev/null +++ b/src/get_dmat.h @@ -0,0 +1,9 @@ +#include +#include +#include +#include +#include +#include + +void get_1rdm(PetscScalar *, PetscInt *, PetscInt *, int *, PetscReal *); +void get_2rdm(PetscScalar *, PetscInt *, PetscInt *, int *, PetscReal *, double ****); diff --git a/src/get_s2.c b/src/get_s2.c new file mode 100644 index 0000000..4fb42f3 --- /dev/null +++ b/src/get_s2.c @@ -0,0 +1,683 @@ +#include +#include +#include +#include +#include +#include +#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=700; + long int iaa2, iaa; + long int iii; + long int ideter[natomax]; + long 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);CHKERRQ(ierr); +//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_s2.h b/src/get_s2.h new file mode 100644 index 0000000..2eb0415 --- /dev/null +++ b/src/get_s2.h @@ -0,0 +1,50 @@ +#include +#include +#include +#include +#include + +void get_s2(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 *); + +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 *); + +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 *); diff --git a/src/get_s2_cyclic.c b/src/get_s2_cyclic.c new file mode 100644 index 0000000..6204863 --- /dev/null +++ b/src/get_s2_cyclic.c @@ -0,0 +1,784 @@ +#include +#include +#include +#include +#include +#include +#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=700; + long int iaa2, iaa; + long int iii; + long int ideter[natomax]; + long 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);CHKERRQ(ierr); +//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..a7154fa --- /dev/null +++ b/src/get_s2_mov.c @@ -0,0 +1,683 @@ +#include +#include +#include +#include +#include +#include +#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=700; + long int iaa2, iaa; + long int iii; + long int ideter[natomax]; + long 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);CHKERRQ(ierr); +//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/getdet.irp.f b/src/getdet.irp.f index cd483af..a4e451c 100644 --- a/src/getdet.irp.f +++ b/src/getdet.irp.f @@ -7,7 +7,7 @@ subroutine getdet(add,ideter) 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::i,const,ia,ib, natom2 ib = MOD(add,nt2) if(MOD(add,nt2).eq.0)then @@ -20,33 +20,38 @@ 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) deta = ISHFT(deta,-(natom/2)) +! 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 + natom2 = natom/2 + 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/natom.irp.f b/src/natom.irp.f index 7b8091a..cd47a04 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,11 @@ 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 !--------------------------------------------- write(6,*)' ' write(6,*)' ' diff --git a/src/providdet.irp.f b/src/providdet.irp.f new file mode 100644 index 0000000..959940f --- /dev/null +++ b/src/providdet.irp.f @@ -0,0 +1,85 @@ +BEGIN_PROVIDER[integer(kind=selected_int_kind(16)),det,(nt2,2)] +&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)),deth,(nt1,2)] + BEGIN_DOC + ! provides det and deth array + END_DOC + implicit none +! integer(kind=selected_int_kind(16))::dethsh + integer(kind=selected_int_kind(16))::a + integer(kind=selected_int_kind(16))::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 + deth(count,1)=ISHFT(a,natom/2) + 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..a529887 100644 --- a/src/provide_clust.irp.f +++ b/src/provide_clust.irp.f @@ -6,6 +6,7 @@ BEGIN_PROVIDER[integer,l1, (maxlien)] &BEGIN_PROVIDER[real*8, xjjxy,(maxlien)] &BEGIN_PROVIDER [integer, ntrou] &BEGIN_PROVIDER [integer, isz] +&BEGIN_PROVIDER [logical*1, FAM1] implicit none ! integer::i ! open(unit=11,file="l1.dat",form="formatted") diff --git a/src/read2.c b/src/read2.c index ec92d41..f2308f0 100644 --- a/src/read2.c +++ b/src/read2.c @@ -1,7 +1,3 @@ -#include -#include -#include -#include #include "read2.h" void Data_new(FILE* file, Data* dat) { @@ -17,25 +13,31 @@ void Data_new(FILE* file, Data* dat) { /* note that fgets don't strip the terminating \n, checking its presence would allow to handle lines longer that sizeof(line) */ - if (count != 12){ + if (count != 26){ count++; switch(count){ case 1: dat->n=atol(line); break; case 2: - dat->nnz=atol(line); + dat->natom=atol(line); break; case 3: - dat->npar=atol(line); + dat->nnz=atol(line); break; case 4: - dat->ntrou=atol(line); + dat->npar=atol(line); break; case 5: - dat->isz=atol(line); + dat->ntrou=atol(line); break; case 6: + dat->isz=atol(line); + break; + case 7: + dat->FAM1 = to_bool(line); + break; + case 8: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -59,7 +61,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 7: + case 9: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -83,7 +85,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 8: + case 10: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -107,7 +109,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 9: + case 11: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -131,7 +133,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 10: + case 12: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -155,7 +157,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 11: + case 13: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -179,9 +181,48 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 12: + 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; } /* end of switch */ } /* end of the input file */ @@ -191,6 +232,14 @@ void Data_new(FILE* file, Data* dat) { //return dat; } +PetscBool to_bool(const char* str) { + PetscBool strflg; + PetscStrcmp("true\n",str, &strflg); + if(!strflg) PetscStrcmp("True\n",str, &strflg); + if(!strflg) PetscStrcmp("TRUE\n",str, &strflg); + return strflg; +} + /* int main(int argc, char* argv[]) { diff --git a/src/read2.h b/src/read2.h index db7021a..dedce9c 100644 --- a/src/read2.h +++ b/src/read2.h @@ -1,13 +1,18 @@ #include -#include #include #include #include +#include +#include + +PetscBool to_bool(const char* str); + typedef struct { PetscInt n; long int nnz,npar; long int ntrou,isz; + PetscBool FAM1; long int l1[700]; long int l2[700]; long int ktyp[700]; @@ -15,6 +20,20 @@ typedef struct { double xjjxy[700]; double xtt[700]; 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; } Data ; diff --git a/src/searchdet.irp.f b/src/searchdet.irp.f index d900984..f8ea9b2 100644 --- a/src/searchdet.irp.f +++ b/src/searchdet.irp.f @@ -1,80 +1,70 @@ -subroutine searchdet(det,add,deth,addh) +subroutine searchdet(deti,add,dethi,addh) 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)::deti 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)::dethi 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(kind=selected_int_kind(16))::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/sort.irp.f b/src/sort.irp.f index 1c8c130..14c7407 100644 --- a/src/sort.irp.f +++ b/src/sort.irp.f @@ -1,14 +1,14 @@ subroutine sort() implicit none integer::i,j,ord,ordh - integer(kind=selected_int_kind(16))::add,addh,det,deth,addt + integer(kind=selected_int_kind(16))::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 +17,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/unit_l1.irp.f b/src/unit_l1.irp.f index 5239a49..2bb5ecf 100644 --- a/src/unit_l1.irp.f +++ b/src/unit_l1.irp.f @@ -9,9 +9,11 @@ tcountcol, & tntrou, & tisz, & + tfam1, & tcol,tval) implicit none integer,INTENT(INOUT)::tistart, tnrows, tntrou, tisz + logical*1,INTENT(INOUT)::tfam1 integer::i real*8,INTENT(INOUT)::tval(maxlien) integer(kind=selected_int_kind(16)),INTENT(INOUT)::tcol(maxlien) @@ -31,6 +33,7 @@ enddo ntrou = tntrou isz = tisz + FAM1 = tfam1 tcol=0 tval=0d0 provide l1 l2 ktyp xtt xjjxy xjjz ntrou @@ -38,6 +41,7 @@ !print *,l1 !print *,"xjjz" !print *,xjjz +!print *,FAM1 call unit(tistart, tcountcol,tcol,tval) end