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