2
0
mirror of https://github.com/LCPQ/DEHam synced 2024-07-07 03:45:59 +02:00

initial commit

This commit is contained in:
vijay gopal chilkuri 2016-12-19 23:26:16 +01:00
commit 79bdcfa221
32 changed files with 2056 additions and 0 deletions

39
adr.irp.f Normal file
View File

@ -0,0 +1,39 @@
subroutine adr(ideter,add)
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,INTENT(INOUT)::ideter(natomax)
integer(kind=selected_int_kind(16)),INTENT(INOUT)::add
integer(kind=selected_int_kind(16))::det,deth,addh,detnew
integer::count,i,j
det=0
detnew=0
deth=0
count=0
call conv(ideter,det,deth)
Do i=0,natom-1
if(BTEST(deth,i))then
count=count+1
endif
if(BTEST(det,i))then
detnew=IBSET(detnew,i-count)
endif
enddo
det=detnew
call searchdet(det,add,deth,addh)
add = add + (nt1-addh)*(nt2)
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

54
adrfull.irp.f Normal file
View File

@ -0,0 +1,54 @@
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))::det,deth,addh,detnew
integer::count,i,j
det=0
detnew=0
deth=0
do j=1,detfound
detnew=0
count=0
ideter=foundet(:,j)
call conv(ideter,det,deth)
Do i=0,natom-1
if(BTEST(deth,i))then
count=count+1
endif
if(BTEST(det,i))then
detnew=IBSET(detnew,i-count)
endif
enddo
det=detnew
foundadd(j,1)=det
foundadd(j,3)=j
foundaddh(j,1)=deth
foundaddh(j,3)=j
enddo
call sort()
call searchdetfull()
call desort()
do i=1,detfound
add = foundadd(i,2)
addh = foundaddh(i,2)
add = add + (nt1-addh)*(nt2)
foundetadr(i)=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

28
combin.irp.f Normal file
View File

@ -0,0 +1,28 @@
SUBROUTINE COMBIN(NS1,N,M,NT,MD,ND)
INTEGER NS1(MD,ND)
INTEGER :: i,j,MMJ1,J1
I=1
DO J=1,M
NS1(j,i)=J
enddo
J1=0
do while (J1.LT.M)
J1=0
DO J=1,M
IF(NS1(j,i).EQ.N-M+J) J1=J1+1
enddo
IF(J1.EQ.M) EXIT
I=I+1
DO J=1,M
NS1(j,i)=NS1(j,I-1)
enddo
NS1(m-j1,I)=NS1(m-j1,I-1)+1
IF(J1.EQ.0) CYCLE
MMJ1=M-J1+1
DO J=MMJ1,M
NS1(j,I)=NS1(j-1,I)+1
enddo
enddo
NT=I
RETURN
END

20
conv.irp.f Normal file
View File

@ -0,0 +1,20 @@
subroutine conv(ideter,det,deth)
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::i
det=0
deth=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)
elseif(ideter(natom-i+1).eq.3)then
deth=IBSET(deth,i-1)
endif
enddo
end

32
desort.irp.f Normal file
View File

@ -0,0 +1,32 @@
subroutine desort()
implicit none
integer::i,j,ord,ordh
integer(kind=selected_int_kind(16))::add,addh,det,deth,addt
do i=1,detfound-1
do j=i+1,detfound
if(foundaddh(i,3).gt.foundaddh(j,3))then
deth = foundaddh(i,1)
foundaddh(i,1) = foundaddh(j,1)
foundaddh(j,1) = deth
addh = foundaddh(i,2)
foundaddh(i,2) = foundaddh(j,2)
foundaddh(j,2) = addh
ordh = foundaddh(i,3)
foundaddh(i,3) = foundaddh(j,3)
foundaddh(j,3) = ordh
endif
if(foundadd(i,3).gt.foundadd(j,3))then
det = foundadd(i,1)
foundadd(i,1) = foundadd(j,1)
foundadd(j,1) = det
add = foundadd(i,2)
foundadd(i,2) = foundadd(j,2)
foundadd(j,2) = add
ord = foundadd(i,3)
foundadd(i,3) = foundadd(j,3)
foundadd(j,3) = ord
endif
enddo
enddo
end

9
deter.irp.f Normal file
View File

@ -0,0 +1,9 @@
BEGIN_PROVIDER [integer, deter, (natomax)]
implicit none
BEGIN_DOC
! provides ideter and iaa
END_DOC
integer ::ideter(natomax)
deter=0
END_PROVIDER

43
elem_diag.irp.f Normal file
View File

@ -0,0 +1,43 @@
subroutine elem_diag(xmatd)
implicit none
integer :: i
real*8 :: xmatd
logical :: yw
! write(6,*)'in elem_diag'
!---
! if(yw)write(6,*)'ideter',(det(iaa,i),i=1,natom)
yw=.FALSE.
xmatd=0.d0
! testing code
! xmatd=1000.d0
! if(.not.redo)write(6,*)'vijayyves'
do i=1,nlientot
if(yalt(i))then
! xmatd=-(xj1+xeneJ(i)*xbJ)+xmatd
xmatd= -(xjz(i))+xmatd
if(yw)write(6,*)xmatd,'xmatd'
if(yw)write(6,*)'xjz',xjz(i)
endif
! if(yrep1(i))then
! xmat=xv1+xmat
! if(yw)write(6,*)iaa,'diag,v1'
! endif
enddo
!-----stockage de l element diag
! imat=iaa
! imat3(isto3+1)=imat
! jmat3(isto3+1)=imat
! dmat3(isto3+1)=xmat/2
return
end

212
ex1.c Normal file
View File

@ -0,0 +1,212 @@
#include <slepceps.h>
#include <petsctime.h>
#include "stimsyr.h"
#include "read2.h"
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
Mat A; /* problem matrix */
EPS eps; /* eigenproblem solver context */
EPSType type;
PetscReal error,tol,re,im;
PetscScalar kr,ki,value[700];
Vec xr,xi;
PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol;
PetscInt nev, ncv, mpd;
PetscLogDouble t1,t2,tt1,tt2;
//PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
PetscErrorCode ierr;
//PetscScalar eigr;
//PetscScalar eigi;
int mpiid;
char const* const fileName = argv[1];
FILE* file = fopen(fileName, "r");
Data getdata;
Data_new(file, &getdata);
//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);
//Vec Vr,Vi;
char filename[PETSC_MAX_PATH_LEN]="FIL666";
PetscViewer viewer;
PetscBool ishermitian;
PetscInt kk,ll,iii2;
long int iii;
long int tcountcol2,tcol[700],tcountcol[getdata.nnz];
double val[700];
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 = 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);
for (i=Istart; i<Iend; i+=getdata.nnz) {
tcountcol2=0;
for(kk=0;kk<getdata.nnz;kk++){
tcountcol[kk]=0;
}
iii=i+1;
if(i%getdata.npar == 0 && mpiid==0){
ierr = PetscTime(&t1);CHKERRQ(ierr);
}
unit_l1_(
getdata.l1,
getdata.l2,
getdata.ktyp,
&iii,
&getdata.nnz,
getdata.xjjxy,
getdata.xjjz ,
getdata.xtt ,
tcountcol,
&getdata.ntrou,
&getdata.isz,
tcol,
val);
if(i%getdata.npar == 0 && mpiid==0){
ierr = PetscTime(&t2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," i: %d\n mpiid: %d\ntime: %f\n",i,mpiid,t2-t1);CHKERRQ(ierr);
}
for(ll=0;ll<getdata.nnz;ll++){
for(kk=0;kk<tcountcol[ll]+1;kk++){
value[kk] = val[kk+tcountcol2];
col[kk] = tcol[kk+tcountcol2]-1;
// PetscPrintf(PETSC_COMM_WORLD,"value = %f col = %d\n",value[kk],col[kk]);
}
for(kk=tcountcol2+tcountcol[ll]+1;kk<700;kk++){
value[kk] = 0.0;
col[kk] = 0;
}
tcountcol2=tcountcol2 + tcountcol[ll]+1;
countcol=tcountcol[ll]+1;
if(i%getdata.npar == 0 && mpiid==0){
ierr = PetscTime(&t1);CHKERRQ(ierr);
}
iii2=i+ll;
ierr = MatSetValues(A,1,&iii2,countcol,col,value,INSERT_VALUES);CHKERRQ(ierr);
if(i%getdata.npar == 0 && mpiid==0){
ierr = PetscTime(&t2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," processor \ntime: %f\n",t2-t1);CHKERRQ(ierr);
}
}
}
ierr = PetscTime(&tt2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to build the matrix: %f\n",tt2-tt1);CHKERRQ(ierr);
ierr = PetscTime(&tt1);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = PetscTime(&tt2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to assemble the matrix: %f\n",tt2-tt1);CHKERRQ(ierr);
//ierr = MatGetVecs(A,NULL,&xr);CHKERRQ(ierr);
//ierr = MatGetVecs(A,NULL,&xi);CHKERRQ(ierr);
ierr = MatCreateVecs(A,NULL,&xr);CHKERRQ(ierr);
ierr = MatCreateVecs(A,NULL,&xi);CHKERRQ(ierr);
ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
ierr = EPSSetOperators(eps,A,NULL);CHKERRQ(ierr);
ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr);
ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr);
ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
tol = 1.e-8;
maxit = 10000000;
ierr = EPSSetTolerances(eps,tol,maxit);CHKERRQ(ierr);
nev = 4;
ncv = 10;
mpd = 10;
ierr = EPSSetDimensions(eps,nev,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
ierr = PetscTime(&t1);CHKERRQ(ierr);
ierr = EPSSolve(eps);CHKERRQ(ierr);
ierr = PetscTime(&t2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Time used: %f\n",t2-t1);CHKERRQ(ierr);
ierr = EPSGetIterationNumber(eps,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);CHKERRQ(ierr);
ierr = EPSGetType(eps,&type);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRQ(ierr);
ierr = EPSGetDimensions(eps,&nev,NULL,NULL);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRQ(ierr);
ierr = EPSGetTolerances(eps,&tol,&maxit);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);CHKERRQ(ierr);
ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr);
//ierr = EPSPrintSolution(eps,NULL);CHKERRQ(ierr);
if (nconv>0) {
/*
Display eigenvalues and relative errors
*/
ierr = PetscPrintf(PETSC_COMM_WORLD,
" k ||Ax-kx||/||kx||\n"
" ----------------- ------------------\n");CHKERRQ(ierr);
for (i=0;i<nconv;i++) {
/*
Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
ki (imaginary part)
*/
ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr);
/*
Compute the relative error associated to each eigenpair
*/
ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
re = PetscRealPart(kr);
im = PetscImaginaryPart(kr);
#else
re = kr;
im = ki;
#endif
if (im!=0.0) {
ierr = PetscPrintf(PETSC_COMM_WORLD," %14f%+14fi %12g\n",(double)re,(double)im,(double)error);CHKERRQ(ierr);
} else {
ierr = PetscPrintf(PETSC_COMM_WORLD," %18f %12g\n",(double)re,(double)error);CHKERRQ(ierr);
}
}
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
}
/*
Save eigenvectors, if requested
*/
//PetscOptionsGetString(NULL,NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
EPSGetConverged(eps,&nconv);
if (nconv>0) {
PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU);
EPSIsHermitian(eps,&ishermitian);
for (i=0;i<nconv;i++) {
EPSGetEigenvector(eps,i,xr,xi);
VecView(xr,viewer);
#if !defined(PETSC_USE_COMPLEX)
if (!ishermitian) { VecView(xi,viewer); }
#endif
}
PetscViewerDestroy(&viewer);
}
ierr = EPSDestroy(&eps);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&xr);CHKERRQ(ierr);
ierr = VecDestroy(&xi);CHKERRQ(ierr);
ierr = SlepcFinalize();
return 0;
}

135
extra_diag.irp.f Normal file
View File

@ -0,0 +1,135 @@
subroutine extra_diag(tistart)
implicit none
integer(kind=selected_int_kind(16)) :: iaa,iaa2,tistart,tistart2
integer(kind=selected_int_kind(16)) :: imat4,jmat4
integer :: i,ik,iik,j
integer :: ik1,ik2,IC,k,ikmax,ikmin,count,count2,detfound2
integer,allocatable :: ideter2(:)
real*8 :: dmat4,xmat
logical :: yw
!---------------------------------------------------------------------
! Calcul des elements extradiagonaux
!---------------------------------------------------------------------
!----------boucle premier voisin
allocate (ideter2(natomax))
foundet=0
foundetadr=0
foundetdmat=0d0
detfound=0
detfound2=0
foundadd=0
foundaddh=0
count=0
count2=1
tistart2=tistart
do j=1,nrows
call getdet(tistart,ideter2)
deter=ideter2
ideter2=0
Touch deter
! call adr(deter,iaa)
! call elem_diag(xmat)
! countcol+=1
! col(countcol)=iaa
! val(countcol)=xmat*1.0d0
count=0
yw=.FALSE.
do ik=1,nlientot
ik1=iliatom1(ik)
ik2=iliatom2(ik)
do k=1,natom
ideter2(k)=deter(k)
enddo
if(yalt(ik)) then
if (deter(ik1).eq.2) then
ideter2(ik1)=1
ideter2(ik2)=2
else
ideter2(ik1)=2
ideter2(ik2)=1
endif
dmat4=xjz(ik)
if(dmat4.ne.0d0)then
count+=1
foundet(:,detfound+count)=ideter2
foundetdmat(detfound+count)=dmat4
endif
endif
if(ytrou(ik)) then
if(deter(ik2).eq.3)then
if (deter(ik1).eq.1) then
ideter2(ik1)=3
ideter2(ik2)=1
else
ideter2(ik1)=3
ideter2(ik2)=2
endif
else
if (deter(ik2).eq.1) then
ideter2(ik1)=1
ideter2(ik2)=3
else
ideter2(ik1)=2
ideter2(ik2)=3
endif
endif
ikmin=min(ik1,ik2)
ikmax=max(ik1,ik2)
IC=0
do iik=ikmin+1,ikmax-1
if(deter(iik).ne.3)IC=IC+1
enddo
dmat4=(xt(ik))*(-1)**(IC)
if(dmat4.ne.0d0)then
count+=1
foundet(:,detfound+count)=ideter2
foundetdmat(detfound+count)=dmat4
endif
endif
enddo
detfound+=count
countcolfull(j)=count
tistart=tistart+1
enddo
Touch foundet foundetadr detfound foundadd foundaddh foundetdmat
call adrfull()
do i=1,detfound
if(i.eq.1 .or. i-1.eq.detfound2)then
call getdet(tistart2,ideter2)
deter=ideter2
ideter2=0
Touch deter
call adr(deter,iaa)
call elem_diag(xmat)
countcol+=1
col(countcol)=iaa
val(countcol)=xmat*1.0d0
tistart2+=1
detfound2+=countcolfull(count2)
! if(i.ne.1)then
count2+=1
! endif
endif
imat4=iaa+1
jmat4=foundetadr(i)
dmat4=foundetdmat(i)
if(jmat4.le.(nt1*nt2) .and. dmat4 .ne. 0d0)then
countcol+=1
col(countcol)=jmat4
val(countcol)=dmat4
endif
enddo
return
end

16
foundet.irp.f Normal file
View File

@ -0,0 +1,16 @@
BEGIN_PROVIDER[integer, foundet,(natomax,maxlien)]
&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundetadr,(maxlien)]
&BEGIN_PROVIDER[real, foundetdmat,(maxlien)]
&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundadd,(maxlien,3)]
&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundaddh,(maxlien,3)]
&BEGIN_PROVIDER[integer, detfound]
BEGIN_DOC
! provides all found determinants
END_DOC
detfound=0
founddet=0
foundetdmat=0d0
founddetadr=0
foundadd=0
foundaddh=0
END_PROVIDER

100
gamma.irp.f Normal file
View File

@ -0,0 +1,100 @@
DOUBLE PRECISION FUNCTION DGAMMA(X)
INTEGER I,N
LOGICAL PARITY
DIMENSION C(7),P(8),Q(8)
DOUBLE PRECISION &
C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, &
TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO ! gamma.irp.f: 95
DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,&
SQRTPI/0.9189385332046727417803297D0/, &
PI/3.1415926535897932384626434D0/ ! gamma.irp.f: 105
DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,&
XINF/1.79D308/ ! gamma.irp.f: 113
DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,&
-3.79804256470945635097577D+2,6.29331155312818442661052D+2, &
8.66966202790413211295064D+2,-3.14512729688483675254357D+4, &
-3.61444134186911729807069D+4,6.64561438202405440627855D+4/ ! gamma.irp.f: 127
DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,&
-1.01515636749021914166146D+3,-3.10777167157231109440444D+3, &
2.25381184209801510330112D+4,4.75584627752788110767815D+3, &
-1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ ! gamma.irp.f: 131
DATA C/-1.910444077728D-03,8.4171387781295D-04,&
-5.952379913043012D-04,7.93650793500350248D-04, &
-2.777777777777681622553D-03,8.333333333333333331554247D-02, &
5.7083835261D-03/ ! gamma.irp.f: 142
I = DBLE(I) ! gamma.irp.f: 150
PARITY = .FALSE. ! gamma.irp.f: 151
FACT = ONE ! gamma.irp.f: 152
N = 0 ! gamma.irp.f: 153
Y = X ! gamma.irp.f: 154
IF (Y .LE. ZERO) THEN ! gamma.irp.f: 155
Y = -X ! gamma.irp.f: 159
Y1 = AINT(Y) ! gamma.irp.f: 160
RES = Y - Y1 ! gamma.irp.f: 161
IF (RES .NE. ZERO) THEN ! gamma.irp.f: 162
IF (Y1 .NE. AINT(Y1*HALF)*TWO) then ! gamma.irp.f: 163
PARITY = .TRUE. ! gamma.irp.f: 163
endif ! gamma.irp.f: 163
FACT = -PI / SIN(PI*RES) ! gamma.irp.f: 164
Y = Y + ONE ! gamma.irp.f: 165
ELSE ! gamma.irp.f: 166
RES = XINF ! gamma.irp.f: 167
GO TO 900 ! gamma.irp.f: 168
END IF ! gamma.irp.f: 169
END IF ! gamma.irp.f: 170
IF (Y .LT. EPS) THEN ! gamma.irp.f: 174
IF (Y .GE. XMININ) THEN ! gamma.irp.f: 178
RES = ONE / Y ! gamma.irp.f: 179
ELSE ! gamma.irp.f: 180
RES = XINF ! gamma.irp.f: 181
GO TO 900 ! gamma.irp.f: 182
END IF ! gamma.irp.f: 183
ELSE IF (Y .LT. TWELVE) THEN ! gamma.irp.f: 184
Y1 = Y ! gamma.irp.f: 185
IF (Y .LT. ONE) THEN ! gamma.irp.f: 186
Z = Y ! gamma.irp.f: 190
Y = Y + ONE ! gamma.irp.f: 191
ELSE ! gamma.irp.f: 192
N = INT(Y) - 1 ! gamma.irp.f: 196
Y = Y - N ! gamma.irp.f: 197
Z = Y - ONE ! gamma.irp.f: 198
END IF ! gamma.irp.f: 199
XNUM = ZERO ! gamma.irp.f: 203
XDEN = ONE ! gamma.irp.f: 204
DO I = 1, 8 ! gamma.irp.f: 205
XNUM = (XNUM + P(I)) * Z ! gamma.irp.f: 206
XDEN = XDEN * Z + Q(I) ! gamma.irp.f: 207
enddo ! gamma.irp.f: 208
RES = XNUM / XDEN + ONE ! gamma.irp.f: 209
IF (Y1 .LT. Y) THEN ! gamma.irp.f: 210
RES = RES / Y1 ! gamma.irp.f: 214
ELSE IF (Y1 .GT. Y) THEN ! gamma.irp.f: 215
DO I = 1, N ! gamma.irp.f: 219
RES = RES * Y ! gamma.irp.f: 220
Y = Y + ONE ! gamma.irp.f: 221
enddo ! gamma.irp.f: 222
END IF ! gamma.irp.f: 223
ELSE ! gamma.irp.f: 224
IF (Y .LE. XBIG) THEN ! gamma.irp.f: 228
YSQ = Y * Y ! gamma.irp.f: 229
SUM = C(7) ! gamma.irp.f: 230
DO I = 1, 6 ! gamma.irp.f: 231
SUM = SUM / YSQ + C(I) ! gamma.irp.f: 232
enddo ! gamma.irp.f: 233
SUM = SUM/Y - Y + SQRTPI ! gamma.irp.f: 234
SUM = SUM + (Y-HALF)*LOG(Y) ! gamma.irp.f: 235
RES = EXP(SUM) ! gamma.irp.f: 236
ELSE ! gamma.irp.f: 237
RES = XINF ! gamma.irp.f: 238
GO TO 900 ! gamma.irp.f: 239
END IF ! gamma.irp.f: 240
END IF ! gamma.irp.f: 241
IF (PARITY) then ! gamma.irp.f: 245
RES = -RES ! gamma.irp.f: 245
endif ! gamma.irp.f: 245
IF (FACT .NE. ONE) then ! gamma.irp.f: 246
RES = FACT / RES ! gamma.irp.f: 246
endif ! gamma.irp.f: 246
900 DGAMMA = RES ! gamma.irp.f: 248
RETURN ! gamma.irp.f: 249
end ! gamma.irp.f: 251

69
getdet.irp.f Normal file
View File

@ -0,0 +1,69 @@
subroutine getdet(add,ideter)
implicit none
BEGIN_DOC
! this routing gives the determinant in
! the traditional form given its address
END_DOC
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
ib = MOD(add,nt2)
if(MOD(add,nt2).eq.0)then
ib=nt2
endif
ia = (add-ib)/nt2
ia = nt1 - ia
ideter=1
const=0
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
const=0
do i=0,(natom/2) - 1
if(BTEST(deta,i))then
ideter((natom/2)-i)=3
endif
enddo
do i=0,natom-1
if(ideter(natom-i).eq.1)then
if(BTEST(detb,const))then
ideter(natom-i)=2
endif
const=const+1
endif
enddo
return
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 subroutine getdet

83
main.irp.f Normal file
View File

@ -0,0 +1,83 @@
program main
implicit none
integer(kind=selected_int_kind(16)),allocatable::tl1(:),tl2(:),tktyp(:)
real*8,allocatable::txtt(:),txjjxy(:),txjjz(:)
integer::i, tnrows, tntrou,tisz
real*4::t1,t2
real*8,allocatable::tval(:)
integer(kind=selected_int_kind(16)),allocatable::tcol(:)
integer(kind=selected_int_kind(16)),dimension(10)::tcountcol
integer(kind=selected_int_kind(16))::tistart
allocate (tl1(maxlien))
allocate (tl2(maxlien))
allocate (tktyp(maxlien))
allocate (tcol(maxlien))
allocate (tval(maxlien))
allocate (txjjxy(maxlien))
allocate (txjjz (maxlien))
allocate (txtt (maxlien))
tl1=0
tl2=0
tktyp=0
txjjxy =0
txjjz =0
txtt =0
tcountcol=0
tcol=0
tval=0d0
tntrou=1
tnrows=10
tisz=0
do i=1,30,tnrows
istart = i
tistart = istart
! tl1=(/1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,12, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,25,0,0,0,0,0/)
! tl2=(/2, 3, 4, 5, 6, 7, 8, 9,10, 11, 12,13,26,25,24,23,22,21,20,19,18,17, 16, 15, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,26,0,0,0,0,0/)
! tktyp=(/1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,0,0,0,0,0/)
! tl1=(/1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,0,0,0,0,0,0,0,0/)
! tl2=(/2, 3, 4, 5, 6, 7, 8, 9,10, 11, 12,24,23,22,21,20,19,18,17,16,15, 14, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,0,0,0,0,0,0,0,0/)
! tktyp=(/1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,0,0,0,0,0,0,0,0/)
! tl1=(/1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tl2=(/2, 3, 4, 5, 6, 7, 8, 9,18,17,16,15,14,13,12,11,10, 11, 12, 13, 14, 15, 16, 17, 18,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tktyp=(/1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tl1=(/1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,18,19,20,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tl2=(/2, 3, 4, 5, 6, 7, 8, 9,10, 11,22,21,20,19,18,17,16,15,14, 13, 12, 13, 14, 15, 16, 17, 18,19,20,21,22,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tktyp=(/1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,3,3,3,3,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
tl1= (/1, 1, 2, 6, 2, 5, 3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
tl2= (/2, 6, 5, 5, 3, 4, 4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
tktyp=(/1, 2, 2, 3, 1, 3, 2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tl1= (/1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5,6,7,14,13,12,11,10, 9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tl2= (/2, 3, 4, 5, 6, 7,14,13,12,11,10,9,8,13,12,11,10, 9, 8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tktyp=(/1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,2,2,3,3,3,3,3,3,3,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
txjjz= (/.1000d0,-0.8d0,0.000d0/)
txjjxy=(/.1000d0,-0.8d0,0.000d0/)
txtt= (/-1.0000d0,0.d0,0.0d0/)
! tl1= (/1, 2, 3, 4, 1, 2, 3,4,5,6,7,8,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tl2= (/2, 3, 4, 5,10, 9, 8,7,6,7,8,9,10,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
! tktyp=(/1, 1, 1, 1, 2, 2, 2,2,2,3,3,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
t1=secnds(0.0)
call unit_l1(tl1, &
tl2, &
tktyp, &
tistart, &
tnrows, &
txjjxy, &
txjjz , &
txtt , &
tcountcol, &
tntrou, &
tisz, &
tcol,tval)
t2=secnds(t1)
print *,'time=',t2
enddo
! deallocate (tl1)
! deallocate (tl2)
! deallocate (tktyp)
! deallocate (tcol)
! deallocate (tval)
! deallocate (txjjxy)
! deallocate (txjjz )
! deallocate (txtt )
end

367
natom.irp.f Normal file
View File

@ -0,0 +1,367 @@
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)]
&BEGIN_PROVIDER [real*8 , xv1]
&BEGIN_PROVIDER [integer, iliatom1,(maxlien)]
&BEGIN_PROVIDER [integer, iliatom2,(maxlien)]
&BEGIN_PROVIDER [integer, ivoisl1,(6,maxlien)]
&BEGIN_PROVIDER [real*8, Emin]
&BEGIN_PROVIDER [real*8, Emax]
&BEGIN_PROVIDER [integer, M0]
&BEGIN_PROVIDER [integer, nalpha]
&BEGIN_PROVIDER [integer, nbeta]
BEGIN_DOC
! read data
END_DOC
implicit none
logical*1 ysuiv
logical*1 yec
logical*1 yw
logical y(natomax)
logical ykl(maxlien)
logical yplac2(maxlien)
logical yplac(maxlien)
! integer,allocatable :: l1(:),l2(:),ktyp(:)
integer :: ityp
!C real*8,allocatable :: xjjz(:),xjjxy(:)
!C real*8,allocatable :: xtt(:)
integer :: iprec,maxda2,niter,ngao,nvec,numero,nes4
integer :: ilien(natomax,natomax)
integer :: nlien,ilien2(natomax,natomax)
real*8 xjxy(maxlien)
integer :: nvois1(natomax),ikx(maxlien)
integer :: iech(natomax,natomax)
integer :: ltyp(maxlien),nvois(natomax),ivois1(6,natomax)
integer :: nvoisl1(maxlien)
integer :: i,ik,ikk,ikl,ik1,ik2,iiki,j,il,il2,il3,il4
integer :: kko,kkio1,kkio2,kkio1ou,klo1,l,nlientot2
integer :: klo2,kat1,kat2,ii,jclust,kl,jk,nplac
integer :: jplac(2,maxlien),iplac(4,maxplac)
real*8 xj1
real*8 xj2
real*8 xt1
real*8 xt2
real*8 xv2
! real*8 xv1
real*8 xv3
real*8 xbJ
real*8 xbT
real*8 xspar
real*8 xsperp
real*8 xseuil
real*8 xeneparJ
real*8 xeneparT
real*8 xeneperpJ
real*8 xeneperpT
real*8 xenediagJ
real*8 xenediagT
!C allocate (xtt(maxlien))
!C allocate (xjjz(maxlien))
!C allocate (xjjxy(maxlien))
!C NAMELIST/hamilton/yham,FAM1
!C NAMELIST/clust/l1,l2,ktyp
!C NAMELIST/integ/ityp,xjjz,xjjxy,xtt
!C NAMELIST/param/xj1,xj2,xt1,xt2,xv1,xv2,xv3,xbJ,xbT,xeneperpJ, &
!C xeneparJ,xenediagJ,xeneperpT,xeneparT,xenediagT,xspar,xsperp
!C NAMELIST/infdia/isz,iprec,maxda2,niter,ngao,nvec,numero,nes4, &
!C xseuil,ysuiv,yec
! NAMELIST/scsrev1/Emin,Emax,M0
ilien2=0
yw=.FALSE.
!-----mise a zero
jclust=0
natom=0
nlientot=0
do ikk=1,maxlien
iliatom1(ikk)=0
iliatom2(ikk)=0
enddo
do ik=1,natomax
do j=1,natomax
ilien(j,ik)=0
iech(j,ik)=0
enddo
enddo
!------------------Lecture Hamiltonien
FAM1=.TRUE.
yham=.TRUE.
write(6,*)'HAMILTONIEN t-J'
write(6,*)'Le nombre de trou est : ',ntrou
!---------------------------------------------
write(6,*)' '
write(6,*)' '
write(6,*)'LECTURE DES ATOMES, DES LIAISONS, DES INTEGRALES'
write(6,*)' '
write(6,*)' '
!-----------Lecture 1ier voisin
! read(5,clust)
jclust=jclust+1
if(yw)write(6,*)' '
write(6,*)'================ CLUSTER',jclust,'=================='
write(6,*)' '
!------------------------
do i=1,maxlien
if(l1(i).eq.0)EXIT
nlien=i
enddo
if(yw)write(6,*)(l1(i),i=1,maxlien)
write(6,*)'Liaisons entre les atomes',nlien
do i=1,natomax
y(i)=.false.
enddo
do i=1,nlien
ilien(l1(i),l2(i))=nlientot+i
ilien(l2(i),l1(i))=nlientot+i
if(l1(i).lt.l2(i))then
iliatom1(nlientot+i)=l1(i)
iliatom2(nlientot+i)=l2(i)
else
iliatom1(nlientot+i)=l2(i)
iliatom2(nlientot+i)=l1(i)
endif
iech(l1(i),l2(i))=(-1)**(l2(i)-l1(i)+1)
iech(l2(i),l1(i))=(-1)**(l2(i)-l1(i)+1)
y(l1(i))=.true.
y(l2(i))=.true.
ltyp(nlientot+i)=ktyp(i)
write(6,*)'Les atomes ',l1(i),' et ',l2(i),' forment la liaison', &
ilien(l1(i),l2(i)),'qui est de type',ltyp(i)
if(yw)write(6,*)'iliatom1',iliatom1(i)
enddo
nlientot=nlientot+nlien
do i=1,natomax
if(y(i))natom=natom+1
enddo
write(6,*)'============================================='
write(6,*)'Le nombre total d atomes est ',natom
write(6,*)'============================================='
!----------------------------------------------------------------
!----------------Voisins---------------------
! nvois(ik)=nbre de voisins de ik
! (ivois(ii,ik),ii=1,nvois(ik))=No des voisins
!--------------------------------------------
do ik=1,natom
nvois1(ik)=0
do ii=1,6
ivois1(ii,ik)=0
enddo
enddo
do ik=1,natom
do ii=1,natom
if (ilien(ik,ii).ne.0) then
nvois1(ik)=nvois1(ik)+1
ivois1(nvois1(ik),ik)=ii
endif
enddo
enddo
!----------------------------------------------------------------
! calcul des liaisons voisines de type 1
!-------------------------------------------------------
do i=1,nlientot
nvoisl1(i)=0
do j=1,6
ivoisl1(j,i)=0
enddo
enddo
do i=1,natom
do j=i+1,natom
if(ilien(i,j).ne.0)then
l=ilien(i,j)
else
EXIT
endif
do ik=1,nvois1(i)
if(ivois1(ik,i).ne.j)then
nvoisl1(l)=nvoisl1(l)+1
ivoisl1(nvoisl1(l),l)=ilien(ivois1(ik,i),i)
endif
enddo
do jk=1,nvois1(j)
if(ivois1(jk,j).ne.i)then
nvoisl1(l)=nvoisl1(l)+1
ivoisl1(nvoisl1(l),l)=ilien(ivois1(jk,j),j)
endif
enddo
enddo
enddo
ityp=3
!C xjjz=(/.1000d0,-0.8d0,0.000d0/)
!C xjjxy=(/.1000d0,-0.8d0,0.000d0/)
!C xtt=(/-1.0000d0,0.d0,0.0d0/)
write(6,*)'Nombre de J differents',ityp
do ikl=1,nlientot
Ykl(ikl)=.false.
enddo
do il=1,nlientot
write(6,*)'type de liaison',il,ltyp(il)
enddo
do iiki=1,ityp
write(6,*)'type de J',xjjz(iiki)
enddo
do il=1,nlientot
if(.not.ykl(il))then
xjz(il)=xjjz(ltyp(il))
write(6,*)'Parametres : Jz',il,'=',xjz(il)
xjxy(il)=xjjxy(ltyp(il))
write(6,*)'Parametres : Jxy',il,'=',xjxy(il)
xt(il)=xtt(ltyp(il))
write(6,*)'Parametre : t',il,'=',xt(il)
ykl(il)=.true.
endif
enddo
xbJ=0.00d0
xbT=0.0d0
xv1=0.0d0
xv2=0.0d0
xv3=0.0d0
xt1=-0.20d0
xt2=0.0d0
xj1=0.01d0
xj2=0.d0
xeneparJ=0.000d0
xeneparT=0.000d0
xeneperpJ=0.000d0
xeneperpT=0.000d0
xenediagJ=0.000d0
xenediagT=0.000d0
xspar=-0.00d0
xsperp=-0.00d0
write(6,*)'coucoudslect3'
write(6,*)'coucou'
write(6,*)'Parametres pour le t-J'
write(6,*)'xj1 = ',xj1
write(6,*)'xj2 = ',xj2
write(6,*)'xt1 = ',xt1
write(6,*)'xt2 = ',xt2
write(6,*)'xv1 = ',xv1
write(6,*)'xv2 = ',xv2
write(6,*)'xv3 = ',xv3
write(6,*)'xbj = ',xbj
write(6,*)'xbt = ',xbt
write(6,*)'xeneparJ = ',xeneparJ
write(6,*)'xeneperpJ = ',xeneperpJ
write(6,*)'xeneparT = ',xeneparT
write(6,*)'xeneperpT = ',xeneperpT
write(6,*)'xenediagJ = ',xenediagJ
write(6,*)'xenediagT = ',xenediagT
write(6,*)'xspar = ',xspar
write(6,*)'xsperp = ',xsperp
!===================================================================
!====================================================================
! calcul des plaquettes pour un 2D t-J
!===================================================================
do ikl=1,nlientot
yplac(ikl)=.false.
IKX(ikl)=0
enddo
nplac=0
do ikl=1,nlientot
ik1=iliatom1(ikl)
ik2=iliatom2(ikl)
do kkio1=1,nvois1(ik1)
do kkio2=1,nvois1(ik2)
kat1=ivois1(kkio1,ik1)
kat2=ivois1(kkio2,ik2)
if(kat1.ne.ik2.and.kat2.ne.ik1)then
if(ilien(kat1,kat2).ne.0)then
klo1=ilien2(kat1,ik2)
if (klo1 == 0) cycle
klo2=ilien2(kat2,ik1)
if (klo2 == 0) cycle
if(.not.yplac2(klo1))then
yplac2(klo1)=.true.
yplac2(klo2)=.true.
nplac=nplac+1
il2=ilien(ik1,kat1)
il3=ilien(ik2,kat2)
il4=ilien(kat1,kat2)
iplac(1,nplac)=ik1
iplac(2,nplac)=ik2
iplac(3,nplac)=kat1
iplac(4,nplac)=kat2
IKX(ikl)=IKX(ikl)+1
IKX(il2)=IKX(il2)+1
IKX(il3)=IKX(il3)+1
IKX(il4)=IKX(il4)+1
jplac(IKX(ikl),ikl)=nplac
jplac(IKX(il2),il2)=nplac
jplac(IKX(il3),il3)=nplac
jplac(IKX(il4),il4)=nplac
endif
endif
endif
enddo
enddo
enddo
write(6,*)'Le systeme comporte ',nplac,' plaquettes.'
do kko=1,nplac
write(6,*)'La plaquette ',kko,' est contituee des atomes',&
iplac(1,kko),' ',iplac(2,kko),' ',iplac(3,kko),' et ',iplac(4,kko)
enddo
!===================================================================
! isz=0
IPREC=8
maxda2=20
NITER=280
ngao=100
nvec=8
numero=1
NES4=0
xseuil=1.0E-008
ysuiv=.FALSE.
yec=.TRUE.
write(6,*)'Spin total',isz
write(6,*)'Nombre de vecteurs demande',nvec
write(6,*)'Nombre maximal d iterations de Davidson',niter
write(6,*)'Vecteur calcule le plus bas',numero
write(6,*)'Variable Nes4 (vecteurs d essai)',nes4
write(6,*)'Nombre de determinants en donnees',ngao
write(6,*)'Variable Ysuiv (suivre le vecteur initial)',ysuiv
write(6,*)'Seuil au dela duquel seront ecrits les vecteurs',xseuil
write(6,*)'Option d ecriture des determinants sur FIL2',yec
! write(6,*)Emin
! write(6,*)Emax
! write(6,*)M0
if(yham)then
IAL0=(natom-ntrou)/2+mod(natom-ntrou,2)+isz
else
IAL0=(natom)/2+mod(natom,2)+isz
endif
write(6,*)'=======nombre de centres de spin alpha=====',ial0
natrest=natom-ntrou
!C calculating nalpha and nbeta
if(mod(natom-ntrou+2*isz,2).eq.0)then
nalpha=(natom-ntrou+2*isz)/2
nbeta=(natom -ntrou-2*isz)/2
if(((natom-ntrou)/2).eq.isz)then
nbeta=0
endif
else
nalpha=(natom-ntrou+2*isz+1)/2
nbeta=(natom -ntrou-2*isz-1)/2
if(((natom-ntrou+1)/2).eq.isz)then
nbeta=0
endif
endif
END_PROVIDER

22
natomax.irp.f Normal file
View File

@ -0,0 +1,22 @@
BEGIN_PROVIDER [integer, natomax]
&BEGIN_PROVIDER [integer, jrangmax]
&BEGIN_PROVIDER [integer, maxial]
&BEGIN_PROVIDER [integer, maxlien]
&BEGIN_PROVIDER [integer, maxplac]
&BEGIN_PROVIDER [integer, maxdet]
&BEGIN_PROVIDER [integer, nrows]
BEGIN_DOC
! provides natomax=22
! provides jrangmax=705432
! provides maxial=20
! provides maxlien=42
END_DOC
implicit none
natomax=700
jrangmax=10705432
maxial=20
maxlien=700
maxplac=20
maxdet=10000
END_PROVIDER

9
nnk.irp.f Normal file
View File

@ -0,0 +1,9 @@
BEGIN_PROVIDER [integer, nnk]
implicit none
BEGIN_DOC
! provides nnk = total number of non zero elements in H
END_DOC
nnk=0
END_PROVIDER

13
nt1.irp.f Normal file
View File

@ -0,0 +1,13 @@
BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt1]
BEGIN_DOC
! calculates the number of det the 3's moving
END_DOC
implicit none
integer(kind=selected_int_kind(16))::natom2
! call combin(idet1(1,nt1+1),natom,ntrou,nt1,32,jrangmax)
natom2=natom
if(FAM1)natom2=natom/2
nt1= nint(gamma(real(natom2+1,16))/(gamma(real(natom2-ntrou+1,16))*gamma(real(ntrou+1,16))),selected_int_kind(16))
write(6,*)'nt1',nt1
END_PROVIDER

10
nt2.irp.f Normal file
View File

@ -0,0 +1,10 @@
BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt2]
BEGIN_DOC
! calculates the number of det the 1's moving
END_DOC
! call combin(idet2(1,nt2+1),natrest,ial0,nt2,32,jrangmax)
nt2= nint(gamma(real(natom-ntrou+1,16))/((gamma(real(nalpha+1,16))*gamma(real(nbeta+1,16)))),selected_int_kind(16))
print *,"nt2=",nt2
END_PROVIDER

8
prodcol.irp.f Normal file
View File

@ -0,0 +1,8 @@
BEGIN_PROVIDER[integer, col,(natomax)]
&BEGIN_PROVIDER[integer, countcol]
&BEGIN_PROVIDER[integer, countcolfull,(nrows)]
implicit none
col=0d0
countcol=0
countcolfull=0
END_PROVIDER

4
prodval.irp.f Normal file
View File

@ -0,0 +1,4 @@
BEGIN_PROVIDER[real*8, val,(natomax)]
implicit none
val=0d0
END_PROVIDER

26
provide_clust.irp.f Normal file
View File

@ -0,0 +1,26 @@
BEGIN_PROVIDER[integer,l1, (maxlien)]
&BEGIN_PROVIDER[integer,l2, (maxlien)]
&BEGIN_PROVIDER[integer,ktyp,(maxlien)]
&BEGIN_PROVIDER [real*8, xtt ,(maxlien)]
&BEGIN_PROVIDER[real*8, xjjz ,(maxlien)]
&BEGIN_PROVIDER[real*8, xjjxy,(maxlien)]
&BEGIN_PROVIDER [integer, ntrou]
&BEGIN_PROVIDER [integer, isz]
implicit none
! integer::i
! open(unit=11,file="l1.dat",form="formatted")
! open(unit=12,file="l2.dat",form="formatted")
! open(unit=13,file="ktyp.dat",form="formatted")
! rewind(11)
! rewind(12)
! rewind(13)
! do i=1,42
! read(11,*)l1(i)
! print *,l1(i)
! read(12,*)l2(i)
! read(13,*)ktyp(i)
! enddo
! close(11)
! close(12)
! close(13)
END_PROVIDER

6
pstartend.irp.f Normal file
View File

@ -0,0 +1,6 @@
BEGIN_PROVIDER[integer, istart]
implicit none
istart=0
END_PROVIDER

15
rank.irp.f Normal file
View File

@ -0,0 +1,15 @@
BEGIN_PROVIDER [integer, rank]
&BEGIN_PROVIDER [integer, rank_16]
BEGIN_DOC
! calculates the rank of matrix
END_DOC
implicit none
rank=nt1*nt2
if(MOD(rank,16).eq.0)then
rank_16=rank
else
rank_16=rank+16-MOD(rank,16)
endif
END_PROVIDER

211
read2.c Normal file
View File

@ -0,0 +1,211 @@
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include "read2.h"
void Data_new(FILE* file, Data* dat) {
//ata* dat = (Data*)malloc(sizeof(Data));
char line[256];
char *token;
const char *delim=",";
size_t arrayIdx=0;
int count=0;
while (fgets(line, sizeof(line), file)) {
/* note that fgets don't strip the terminating \n, checking its
presence would allow to handle lines longer that sizeof(line) */
if (count != 11){
count++;
switch(count){
case 1:
dat->n=atol(line);
break;
case 2:
dat->nnz=atol(line);
break;
case 3:
dat->npar=atol(line);
break;
case 4:
dat->ntrou=atol(line);
break;
case 5:
dat->isz=atol(line);
break;
case 6:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
double val;
char *unconverted;
/**
* Convert the next token to a float value
*/
val = strtol(token, &unconverted, 10);
if (!isspace(*unconverted) && *unconverted != 0)
{
/**
* Bad input string. Again, we just bail.
*/
fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
break;
}
else
{
dat->l1[arrayIdx++] = val;
}
}
break;
case 7:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
double val;
char *unconverted;
/**
* Convert the next token to a float value
*/
val = strtol(token, &unconverted, 10);
if (!isspace(*unconverted) && *unconverted != 0)
{
/**
* Bad input string. Again, we just bail.
*/
fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
break;
}
else
{
dat->l2[arrayIdx++] = val;
}
}
break;
case 8:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
double val;
char *unconverted;
/**
* Convert the next token to a float value
*/
val = strtol(token, &unconverted, 10);
if (!isspace(*unconverted) && *unconverted != 0)
{
/**
* Bad input string. Again, we just bail.
*/
fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
break;
}
else
{
dat->ktyp[arrayIdx++] = val;
}
}
break;
case 9:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
double val;
char *unconverted;
/**
* Convert the next token to a float value
*/
val = strtof(token, &unconverted);
if (!isspace(*unconverted) && *unconverted != 0)
{
/**
* Bad input string. Again, we just bail.
*/
fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
break;
}
else
{
dat->xjjz[arrayIdx++] = val;
}
}
break;
case 10:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
double val;
char *unconverted;
/**
* Convert the next token to a float value
*/
val = strtof(token, &unconverted);
if (!isspace(*unconverted) && *unconverted != 0)
{
/**
* Bad input string. Again, we just bail.
*/
fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
break;
}
else
{
dat->xjjxy[arrayIdx++] = val;
}
}
break;
case 11:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
double val;
char *unconverted;
/**
* Convert the next token to a float value
*/
val = strtof(token, &unconverted);
if (!isspace(*unconverted) && *unconverted != 0)
{
/**
* Bad input string. Again, we just bail.
*/
fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
break;
}
else
{
dat->xtt[arrayIdx++] = val;
}
}
break;
} /* end of switch */
} /* end of the input file */
} /* end of loop */
//return dat;
}
/*
int main(int argc, char* argv[])
{
char const* const fileName = argv[1];
FILE* file = fopen(fileName, "r");
Data *getdata=Data_new(file);
fclose(file);
printf("n=%d,nnz=%d,npar=%d\n",getdata->n,getdata->nnz,getdata->npar);
for(int i=0;i<7;i++){
printf("l1(%d)=%d l2(%d)=%d l3(%d)=%d\n",i,getdata->l1[i],i,getdata->l2[i],i,getdata->l3[i]);
}
for(int i=0;i<7;i++){
printf("xjjz(%d)=%f xjjxy(%d)=%f xtt(%d)=%f\n",i,getdata->xjjz[i],i,getdata->xjjxy[i],i,getdata->xtt[i]);
}
return 0;
}
*/

20
read2.h Normal file
View File

@ -0,0 +1,20 @@
#include <stdio.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
typedef struct {
PetscInt n;
long int nnz,npar;
long int ntrou,isz;
long int l1[700];
long int l2[700];
long int ktyp[700];
double xjjz[700];
double xjjxy[700];
double xtt[700];
} Data ;
void Data_new(FILE* , Data* );

86
searchdet.irp.f Normal file
View File

@ -0,0 +1,86 @@
subroutine searchdet(det,add,deth,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)::add
integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth
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
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
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
else
add=i
count=-1
EXIT
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)
11 FORMAT(B64,I3,B64)
12 FORMAT(I5,$)
13 FORMAT(B64,B64)
14 FORMAT(B64,I8)
16 FORMAT(B64,I8,I8)
end

140
searchdetfull.irp.f Normal file
View File

@ -0,0 +1,140 @@
subroutine searchdetfull()
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)::foundetadr(maxlien,4)
integer(kind=selected_int_kind(16))::add
! integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth
integer(kind=selected_int_kind(16))::dethsh
integer(kind=selected_int_kind(16))::addh
integer(kind=selected_int_kind(16))::a
integer(kind=selected_int_kind(16))::i
integer::const,count
i=1
a=0
add=0
const=0
count=1
If(ntrou.ge.1)then
const=0
a=0
addh=0
i=1
dethsh = ISHFT(foundaddh(count,1),-natom/2)
do while (i.le.(2*nt1))
if(a.eq.dethsh)then
addh=i-1
foundaddh(count,2)=addh
count+=1
dethsh = ISHFT(foundaddh(count,1),-natom/2)
do while(count .le. detfound .and. a .eq. dethsh)
foundaddh(count,2)=addh
count+=1
dethsh = ISHFT(foundaddh(count,1),-natom/2)
enddo
if(count.gt.detfound)then
EXIT
endif
endif
i+=1
! const=1
a+=1
! do while(popcnt(a).ne.ntrou .or. const==1)
do while(popcnt(a).ne.ntrou)
a+=1
! const=0
enddo
enddo
if(a.eq.foundaddh(count,1))then
addh=i-1
foundaddh(count,2)=addh
count+=1
dethsh = ISHFT(foundaddh(count,1),-natom/2)
do while(count .le. detfound .and. a .eq. foundaddh(count,1))
foundaddh(count,2)=addh
count+=1
dethsh = ISHFT(foundaddh(count,1),-natom/2)
enddo
endif
endif
!C if det=0 then exit
a=0
i=1
count=1
const=0
if(a.eq.foundadd(count,1))then
add=1
foundadd(count,2)=add
count+=1
do while(foundadd(count,1).eq.a .and. count.le.detfound)
foundadd(count,2)=add
count+=1
enddo
endif
do while (i.le.(nt2) .and. count .le.detfound)
if(a.eq.foundadd(count,1))then
if(a.eq.1)then
add=i-1
foundadd(count,2)=add
count+=1
do while(foundadd(count,1).eq.a .and. count.le.detfound)
foundadd(count,2)=add
count+=1
enddo
if(count.eq.detfound)then
const=-1
EXIT
endif
else
add=i-1
foundadd(count,2)=add
count+=1
do while(foundadd(count,1).eq.a .and. count.le.detfound)
foundadd(count,2)=add
count+=1
enddo
if(count.gt.detfound)then
const=-1
EXIT
endif
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.foundadd(count,1) .and. foundadd(count,1).eq.0)then
if(a.eq.foundadd(count,1) .and. foundadd(count,2) .eq. 0)then
foundadd(count,2) = nt2
count+=1
do while(foundadd(count,1).eq.a .and. count.le.detfound)
foundadd(count,2)=nt2
count+=1
enddo
endif
! do i=1,detfound
! write(6,16)foundadd(i,1) ,foundadd(i,2),foundadd(i,3)
! write(6,16)foundaddh(i,1),foundaddh(i,2),foundaddh(i,3)
! 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

32
sort.irp.f Normal file
View File

@ -0,0 +1,32 @@
subroutine sort()
implicit none
integer::i,j,ord,ordh
integer(kind=selected_int_kind(16))::add,addh,det,deth,addt
do i=1,detfound-1
do j=i+1,detfound
if(foundaddh(i,1).gt.foundaddh(j,1))then
deth = foundaddh(i,1)
foundaddh(i,1) = foundaddh(j,1)
foundaddh(j,1) = deth
addh = foundaddh(i,2)
foundaddh(i,2) = foundaddh(j,2)
foundaddh(j,2) = addh
ordh = foundaddh(i,3)
foundaddh(i,3) = foundaddh(j,3)
foundaddh(j,3) = ordh
endif
if(foundadd(i,1).gt.foundadd(j,1))then
det = foundadd(i,1)
foundadd(i,1) = foundadd(j,1)
foundadd(j,1) = det
add = foundadd(i,2)
foundadd(i,2) = foundadd(j,2)
foundadd(j,2) = add
ord = foundadd(i,3)
foundadd(i,3) = foundadd(j,3)
foundadd(j,3) = ord
endif
enddo
enddo
end

15
stimsyr.h Normal file
View File

@ -0,0 +1,15 @@
void unit_l1_(
long int *,
long int *,
long int *,
long int *,
long int *,
double *,
double *,
double *,
long int *,
long int *,
long int *,
long int *,
double *);

93
unit_FIL44.irp.f Normal file
View File

@ -0,0 +1,93 @@
subroutine unit(tistart,tcountcol,tcol,tval)
BEGIN_DOC
! file units for writing
END_DOC
implicit none
integer :: i,j,k,ia1,ia2,l,m,chcind,chcval,ii,tistart2
integer :: count,unit_44,unit_33
integer :: iat,nbtots
integer(kind=selected_int_kind(16))::iaa
integer :: kkio,kkiok,n,nz,cdiag,cexdiag
integer,allocatable ::ideter1(:),ideter2(:),deti(:),detj(:)
integer(kind=selected_int_kind(16)),dimension(maxlien) ::tl1,tl2,tktyp
integer(kind=selected_int_kind(16)),dimension(nrows)::tcountcol
integer(kind=selected_int_kind(16))::tistart
real*8,dimension(maxlien)::tval
integer(kind=selected_int_kind(16)),dimension(maxlien)::tcol
real*8 :: xmat
integer :: ik,imat4,iaa2,iik
integer :: ik1,ik2,jmat4,IC,ikmax,ikmin
real*8 :: dmat4
logical :: yw
! BEGIN_DOC
! provides unit of FIL33 & FIL44
! END_DOC
allocate (ideter2(natomax))
! allocate (tcol(natomax))
! allocate (tval(natomax))
countcol=0
unit_44=44
unit_33=33
nnk=0
xmat=0d0
count=0
cdiag=1
cexdiag=0
tistart2=tistart
do i=1,natomax
col(i)=0
val(i)=0d0
tval(i)=0d0
tcol(i)=0
enddo
tcountcol=0
countcol=0
xmat=0d0
count=0
! tistart=tistart
i=1+tistart/nt2
k=1+mod(tistart , nt2)
! call getdet(tistart,ideter2)
! deter=ideter2
! Touch deter
! call adr(deter,iaa)
! call elem_diag(xmat)
! countcol+=1
! col(countcol)=iaa
! val(countcol)=xmat*1.0d0
call extra_diag(tistart)
tcountcol=countcolfull
do i=1,maxlien
if(col(i).ne.0)then
if(val(i) .ne. 0 .or. col(i).eq.tistart2)then
tcol(i)=col(i)
tval(i)=val(i)
endif
if(col(i).eq.tistart2)then
cexdiag+=1
elseif(cexdiag .eq. countcolfull(cdiag))then
cexdiag=0
if(cdiag.lt.nrows)then
cdiag+=1
endif
tistart2+=1
else
cexdiag+=1
endif
endif
enddo
! print *,"tistart =", tistart,"countcol =", countcol,"\n",(tcountcol(i),i=1,nrows)
! print *,""
! print *,(tcol(i),i=1,maxlien)
deallocate(ideter2)
end

45
unit_l1.irp.f Normal file
View File

@ -0,0 +1,45 @@
subroutine unit_l1(tl1,&
tl2, &
tktyp, &
tistart, &
tnrows, &
txjjxy, &
txjjz , &
txtt , &
tcountcol, &
tntrou, &
tisz, &
tcol,tval)
implicit none
integer,INTENT(INOUT)::tistart, tnrows, tntrou, tisz
integer::i
real*8,INTENT(INOUT)::tval(maxlien)
integer(kind=selected_int_kind(16)),INTENT(INOUT)::tcol(maxlien)
integer(kind=selected_int_kind(16)),INTENT(INOUT),dimension(tnrows)::tcountcol
integer(kind=selected_int_kind(16)),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien)
real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien)
nrows = tnrows
provide nrows
do i=1,maxlien
l1(i)=tl1(i)
l2(i)=tl2(i)
ktyp(i)=tktyp(i)
xtt(i) = txtt(i)
xjjxy(i) = txjjxy(i)
xjjz (i) = txjjz (i)
enddo
ntrou = tntrou
isz = tisz
tcol=0
tval=0d0
provide l1 l2 ktyp xtt xjjxy xjjz ntrou
!print *,"l1"
!print *,l1
!print *,"xjjz"
!print *,xjjz
call unit(tistart, tcountcol,tcol,tval)
end

94
ylogic.irp.f Normal file
View File

@ -0,0 +1,94 @@
! SUBROUTINE ylogic(ideter,yalt,ytrou,yrep1)
BEGIN_PROVIDER [logical, yalt,(maxlien)]
&BEGIN_PROVIDER [logical, ytrou,(maxlien)]
&BEGIN_PROVIDER [logical, yrep1,(maxlien)]
implicit none
BEGIN_DOC
! provides the important logical units
! for the presence of holes (ytrou) and
! an alternance (yalt) or a rep (yrep)
END_DOC
integer :: i,il
integer :: ik1,ik2
integer :: ntr,nalt,naltmax
integer :: itypl(maxlien)
logical :: yw
! write(6,*)'in elem_diag'
!---mise a zero
nalt=0
ntr=0
yw=.FALSE.
do il=1,nlientot
itypl(il)=0
yalt(il)=.false.
ytrou(il)=.false.
yrep1(il)=.false.
enddo
!---
if(yw)write(6,*)'ideter',(deter(il),il=1,natom)
do il=1,nlientot
ik1=iliatom1(il)
ik2=iliatom2(il)
if(yw)write(6,*)'ik1',ik1,'ik2',ik2
if(deter(ik1).eq.deter(ik2).and. &
deter(ik1).eq.3)then
yrep1(il)=.true.
endif
if(deter(ik1).ne.deter(ik2))then
if(deter(ik1).eq.1.and.deter(ik2).eq.2)then
yalt(il)=.true.
nalt=nalt+1
itypl(il)=1
! goto 32
CYCLE
endif
if(deter(ik1).eq.2.and.deter(ik2).eq.1)then
yalt(il)=.true.
nalt=nalt+1
itypl(il)=2
goto 32
CYCLE
endif
32 continue
if(deter(ik1).eq.1.and.deter(ik2).eq.3)then
ytrou(il)=.true.
ntr=ntr+1
itypl(il)=3
! goto 33
CYCLE
endif
if(deter(ik1).eq.3.and.deter(ik2).eq.1)then
ytrou(il)=.true.
ntr=ntr+1
itypl(il)=4
! goto 33
CYCLE
endif
if(deter(ik1).eq.2.and.deter(ik2).eq.3)then
ytrou(il)=.true.
ntr=ntr+1
itypl(il)=5
! goto 33
CYCLE
endif
if(deter(ik1).eq.3.and.deter(ik2).eq.2)then
ytrou(il)=.true.
ntr=ntr+1
itypl(il)=6
! goto 33
CYCLE
endif
endif
!33 continue
enddo
END_PROVIDER