From ee3afdd71e34e34a54676f20ec9f1292d6c8d316 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Tue, 20 Dec 2016 16:02:23 +0100 Subject: [PATCH 01/38] modified makefile to seperate into bin and lib dirs also now added neigenvalues option --- Makefile | 38 ++++++++++++++++++++++++++------------ src/ex1.c | 5 +++-- src/read2.c | 5 ++++- src/read2.h | 1 + 4 files changed, 34 insertions(+), 15 deletions(-) diff --git a/Makefile b/Makefile index fba9d00..99f2424 100644 --- a/Makefile +++ b/Makefile @@ -1,23 +1,37 @@ include ${SLEPC_DIR}/lib/slepc/conf/slepc_common -CC=gcc +#CC=gcc +#FC = ifort MAKE = /usr/bin/make -MKDIR_P = mkdir -p -OBJ_DIR = obj +MKDIR_P = /bin/mkdir -p +OBJ_DIR := obj +LIB_DIR := libs +BIN_DIR := bin +SRC_DIR := src -irpf90.a: - cd src && irpf90 init && $(MAKE) && cp irpf90.a ../ +.PHONY: ex1 + +ex1: ${BIN_DIR}/ex1 ${OBJ_DIR}: ${MKDIR_P} ${OBJ_DIR} -directories: ${OBJ_DIR} +${LIB_DIR}: + ${MKDIR_P} ${LIB_DIR} -obj/read2.o: src/read2.c directories - ${CC} -c -o $@ $< +${BIN_DIR}: + ${MKDIR_P} ${BIN_DIR} -obj/ex1.o: src/ex1.c - -${CC} -c -o $@ $< ${SLEPC_EPS_LIB} +directories: ${OBJ_DIR} ${LIB_DIR} ${BIN_DIR} -ex1: obj/read2.o irpf90.a obj/ex1.o src/read2.h src/stimsyr.h chkopts - -${CLINKER} -o ex1 obj/ex1.o obj/read2.o irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3 +${LIB_DIR}/irpf90.a: directories + cd ${SRC_DIR} && irpf90 init && $(MAKE) irpf90.a && cp irpf90.a ../${LIB_DIR} + +${OBJ_DIR}/read2.o: ${SRC_DIR}/read2.c directories chkopts + ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + +${OBJ_DIR}/ex1.o: ${SRC_DIR}/ex1.c + -${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + +${BIN_DIR}/ex1: ${OBJ_DIR}/read2.o ${LIB_DIR}/irpf90.a ${OBJ_DIR}/ex1.o ${SRC_DIR}/read2.h ${SRC_DIR}/stimsyr.h chkopts + -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3 # ${RM} ex1.o read2.o diff --git a/src/ex1.c b/src/ex1.c index b29a2b5..8d3c3a9 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -124,12 +124,13 @@ int main(int argc,char **argv) ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr); ierr = EPSSetFromOptions(eps);CHKERRQ(ierr); - tol = 1.e-8; + tol = 1.e-9; maxit = 10000000; ierr = EPSSetTolerances(eps,tol,maxit);CHKERRQ(ierr); nev = 4; ncv = 10; mpd = 10; + nev = getdata.nroots; ierr = EPSSetDimensions(eps,nev,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); ierr = PetscTime(&t1);CHKERRQ(ierr); @@ -156,7 +157,7 @@ int main(int argc,char **argv) " k ||Ax-kx||/||kx||\n" " ----------------- ------------------\n");CHKERRQ(ierr); - for (i=0;inroots=atol(line); + break; } /* end of switch */ } /* end of the input file */ diff --git a/src/read2.h b/src/read2.h index 6c98ab4..db7021a 100644 --- a/src/read2.h +++ b/src/read2.h @@ -14,6 +14,7 @@ typedef struct { double xjjz[700]; double xjjxy[700]; double xtt[700]; + long int nroots; } Data ; From 85b10c2f62241e7a48c0f59713bb925821bd645a Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Wed, 21 Dec 2016 11:41:15 +0100 Subject: [PATCH 02/38] added calculation of norm of resulting eigenvector --- src/ex1.c | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/ex1.c b/src/ex1.c index 8d3c3a9..af1e530 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -12,6 +12,7 @@ int main(int argc,char **argv) EPS eps; /* eigenproblem solver context */ EPSType type; PetscReal error,tol,re,im; + PetscReal norm; PetscScalar kr,ki,value[700]; Vec xr,xi; PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol; @@ -204,6 +205,26 @@ int main(int argc,char **argv) PetscViewerDestroy(&viewer); } + /* + * now analyzing the eigenvector + */ + + if (nconv>0) { + + for (i=0;i Date: Wed, 21 Dec 2016 11:53:13 +0100 Subject: [PATCH 03/38] trying to read array elements --- src/ex1.c | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index af1e530..36f0a77 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -211,17 +211,19 @@ int main(int argc,char **argv) if (nconv>0) { - for (i=0;i Date: Wed, 21 Dec 2016 11:54:14 +0100 Subject: [PATCH 04/38] add header to manipulate vectors --- src/ex1.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ex1.c b/src/ex1.c index 36f0a77..608933b 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -1,5 +1,6 @@ #include #include +#include #include "stimsyr.h" #include "read2.h" From 5185407a4c728896bce9e33dd1bc9b918a5f5c82 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Wed, 21 Dec 2016 11:55:04 +0100 Subject: [PATCH 05/38] values instead --- src/ex1.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ex1.c b/src/ex1.c index 608933b..9dfdf05 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -221,7 +221,7 @@ int main(int argc,char **argv) ierr = VecNorm(xr, NORM_2, &norm);CHKERRQ(ierr); PetscPrintf(PETSC_COMM_WORLD," Norm = %18f \n", (double)norm); for (i=Istart; i Date: Wed, 21 Dec 2016 11:56:39 +0100 Subject: [PATCH 06/38] add error checking --- src/ex1.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ex1.c b/src/ex1.c index 9dfdf05..88f7721 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -221,7 +221,7 @@ int main(int argc,char **argv) ierr = VecNorm(xr, NORM_2, &norm);CHKERRQ(ierr); PetscPrintf(PETSC_COMM_WORLD," Norm = %18f \n", (double)norm); for (i=Istart; i Date: Wed, 21 Dec 2016 12:17:40 +0100 Subject: [PATCH 07/38] now reads local values of eigenvector --- src/ex1.c | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index 88f7721..20b5a70 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -28,11 +28,15 @@ int main(int argc,char **argv) char const* const fileName = argv[1]; FILE* file = fopen(fileName, "r"); Data getdata; + PetscInt nlocal; 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); + PetscScalar valxr[nlocal]; + PetscInt indxr[nlocal]; //Vec Vr,Vi; char filename[PETSC_MAX_PATH_LEN]="FIL666"; PetscViewer viewer; @@ -220,9 +224,13 @@ int main(int argc,char **argv) ierr = VecNorm(xr, NORM_2, &norm);CHKERRQ(ierr); PetscPrintf(PETSC_COMM_WORLD," Norm = %18f \n", (double)norm); - for (i=Istart; i Date: Wed, 21 Dec 2016 13:08:15 +0100 Subject: [PATCH 08/38] now gets the determinant for a given adress --- src/ex1.c | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index 20b5a70..afb5e93 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -24,6 +24,8 @@ int main(int argc,char **argv) //PetscScalar eigr; //PetscScalar eigi; int mpiid; + int natomax=700; + int ideter[natomax]; char const* const fileName = argv[1]; FILE* file = fopen(fileName, "r"); @@ -43,6 +45,7 @@ int main(int argc,char **argv) PetscBool ishermitian; PetscInt kk,ll,iii2; long int iii; + long int ii; long int tcountcol2,tcol[700],tcountcol[getdata.nnz]; double val[700]; @@ -228,9 +231,14 @@ int main(int argc,char **argv) indxr[i-Istart] = i; } ierr = VecGetValues(xr, nlocal, indxr, valxr);CHKERRQ(ierr); - for (i=Istart; i Date: Wed, 21 Dec 2016 13:33:13 +0100 Subject: [PATCH 09/38] add error checking --- src/ex1.c | 69 +++++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 60 insertions(+), 9 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index afb5e93..5670b61 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -13,7 +13,7 @@ int main(int argc,char **argv) EPS eps; /* eigenproblem solver context */ EPSType type; PetscReal error,tol,re,im; - PetscReal norm; + PetscReal norm=0.0; PetscScalar kr,ki,value[700]; Vec xr,xi; PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol; @@ -24,8 +24,9 @@ int main(int argc,char **argv) //PetscScalar eigr; //PetscScalar eigi; int mpiid; - int natomax=700; + int natomax=700, iaa2, iaa; int ideter[natomax]; + int ideter2[natomax]; char const* const fileName = argv[1]; FILE* file = fopen(fileName, "r"); @@ -48,6 +49,9 @@ int main(int argc,char **argv) long int ii; long int tcountcol2,tcol[700],tcountcol[getdata.nnz]; double val[700]; + PetscReal xmat=0.0; + PetscReal xymat=0.0; + PetscInt kko,kok,kkio; SlepcInitialize(&argc,&argv,(char*)0,NULL); ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr); @@ -151,7 +155,7 @@ int main(int argc,char **argv) 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 = PetscPrintf(PETSC_COMM_WORLD," Number of == ested 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); @@ -160,7 +164,7 @@ int main(int argc,char **argv) if (nconv>0) { /* - Display eigenvalues and relative errors + Display eigenvalues && relative errors */ ierr = PetscPrintf(PETSC_COMM_WORLD, " k ||Ax-kx||/||kx||\n" @@ -194,7 +198,7 @@ int main(int argc,char **argv) } /* - Save eigenvectors, if requested + Save eigenvectors, if == ested */ //PetscOptionsGetString(NULL,NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs); EPSGetConverged(eps,&nconv); @@ -231,18 +235,65 @@ int main(int argc,char **argv) indxr[i-Istart] = i; } ierr = VecGetValues(xr, nlocal, indxr, valxr);CHKERRQ(ierr); + printf("istart = %d, iend = %d\n",Istart, Iend); for (ii=Istart; ii Date: Wed, 21 Dec 2016 13:59:01 +0100 Subject: [PATCH 10/38] now gets the determinant for a given adress --- src/ex1.c | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index 5670b61..d1f7283 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -14,6 +14,7 @@ int main(int argc,char **argv) EPSType type; PetscReal error,tol,re,im; PetscReal norm=0.0; + PetscReal normfin=0.0; PetscScalar kr,ki,value[700]; Vec xr,xi; PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol; @@ -49,9 +50,11 @@ int main(int argc,char **argv) long int ii; long int tcountcol2,tcol[700],tcountcol[getdata.nnz]; double val[700]; - PetscReal xmat=0.0; + double xmat=0.0; PetscReal xymat=0.0; - PetscInt kko,kok,kkio; + PetscReal xymatfin=0.0; + PetscReal XS = 0.0; + int kko,kok,kkio; SlepcInitialize(&argc,&argv,(char*)0,NULL); ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr); @@ -231,11 +234,11 @@ int main(int argc,char **argv) ierr = VecNorm(xr, NORM_2, &norm);CHKERRQ(ierr); PetscPrintf(PETSC_COMM_WORLD," Norm = %18f \n", (double)norm); + norm = 0.0; for (i=Istart; i Date: Wed, 21 Dec 2016 14:21:38 +0100 Subject: [PATCH 11/38] add S2 functionality but looks like theres a bug --- src/ex1.c | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index d1f7283..fa370f1 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -210,7 +210,7 @@ int main(int argc,char **argv) PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB); PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU); EPSIsHermitian(eps,&ishermitian); - for (i=0;i0) { + for(i=0;i Date: Wed, 21 Dec 2016 14:40:34 +0100 Subject: [PATCH 12/38] still working on the bug --- src/ex1.c | 89 +++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 80 insertions(+), 9 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index fa370f1..c0c7b23 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -234,9 +234,6 @@ int main(int argc,char **argv) */ ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr); xymat = 0.0; - - ierr = VecNorm(xr, NORM_2, &norm);CHKERRQ(ierr); - PetscPrintf(PETSC_COMM_WORLD," Norm = %18f \n", (double)norm); norm = 0.0; @@ -269,7 +266,7 @@ int main(int argc,char **argv) } if(ideter[kko] == 1 && ideter[kok] == 2){ xmat=xmat-(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); - for(kkio=1;kkio<=8;kkio++){ + for(kkio=0;kkio<=7;kkio++){ ideter2[kkio]=ideter[kkio]; } ideter2[kko]=2; @@ -280,7 +277,7 @@ int main(int argc,char **argv) } if(ideter[kko] == 2 && ideter[kok] == 1){ xmat=xmat-(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); - for(kkio=1;kkio<=8;kkio++){ + for(kkio=0;kkio<=7;kkio++){ ideter2[kkio]=ideter[kkio]; } ideter2[kko]=1; @@ -292,18 +289,92 @@ int main(int argc,char **argv) } } } +// for(kko=4;kko<=7;kko++){ +// for(kok=kko;kok<=7;kok++){ +// if(kok == kko && ideter[kok] != 3){ +// xmat=xmat+(3.0/4.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// } +// else{ +// if(ideter[kko] == 1 && ideter[kok] == 1){ +// xmat=xmat+(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// } +// if(ideter[kko] == 2 && ideter[kok] == 2){ +// xmat=xmat+(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// } +// if(ideter[kko] == 1 && ideter[kok] == 2){ +// xmat=xmat-(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// for(kkio=0;kkio<=7;kkio++){ +// ideter2[kkio]=ideter[kkio]; +// } +// ideter2[kko]=2; +// ideter2[kok]=1; +// adr_(ideter2, &iaa2); +// iaa2 = iaa2 - 1; +// xmat=xmat+valxr[ii-Istart]*valxr[iaa2]; +// } +// if(ideter[kko] == 2 && ideter[kok] == 1){ +// xmat=xmat-(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// for(kkio=0;kkio<=7;kkio++){ +// ideter2[kkio]=ideter[kkio]; +// } +// ideter2[kko]=1; +// ideter2[kok]=2; +// adr_(ideter2, &iaa2); +// iaa2 = iaa2 - 1; +// xmat=xmat+valxr[ii-Istart]*valxr[iaa2]; +// } +// } +// } +// } +// for(kko=0;kko<=3;kko++){ +// for(kok=4;kok<=7;kok++){ +// if(kok == kko && ideter[kok] != 3){ +// xmat=xmat+(3.0/4.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// } +// else{ +// if(ideter[kko] == 1 && ideter[kok] == 1){ +// xmat=xmat+(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// } +// if(ideter[kko] == 2 && ideter[kok] == 2){ +// xmat=xmat+(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// } +// if(ideter[kko] == 1 && ideter[kok] == 2){ +// xmat=xmat-(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// for(kkio=0;kkio<=7;kkio++){ +// ideter2[kkio]=ideter[kkio]; +// } +// ideter2[kko]=2; +// ideter2[kok]=1; +// adr_(ideter2, &iaa2); +// iaa2 = iaa2 - 1; +// xmat=xmat+valxr[ii-Istart]*valxr[iaa2]; +// } +// if(ideter[kko] == 2 && ideter[kok] == 1){ +// xmat=xmat-(1.0/2.0)*(valxr[ii-Istart]*valxr[ii-Istart]); +// for(kkio=0;kkio<=7;kkio++){ +// ideter2[kkio]=ideter[kkio]; +// } +// ideter2[kko]=1; +// ideter2[kok]=2; +// adr_(ideter2, &iaa2); +// iaa2 = iaa2 - 1; +// xmat=xmat+valxr[ii-Istart]*valxr[iaa2]; +// } +// } +// } +// } //--------------------------------------- xymat=xymat+xmat; } } - MPI_Reduce(&xymat, &xymatfin, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); - MPI_Reduce(&norm, &normfin, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&xymat, &xymatfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD); + MPI_Reduce(&norm, &normfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD); if(!mpiid){ XS=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin/normfin))); } -//PetscBarrier((PetscObject)&XS); -//PetscPrintf(PETSC_COMM_WORLD,"\n # root = %d norm = %18f xymat = %18f S^2 = %18f \n", i, normfin, xymatfin, XS); + + printf("\n mpiid = %d # root = %d norm = %18f xymat = %18f S^2 = %18f \n", mpiid, i, norm, xymat, XS); PetscPrintf(PETSC_COMM_WORLD,"\n norm = %18f xymat = %18f S^2 = %18f \n", normfin, xymatfin, XS); xymatfin = 0.0; normfin = 0.0; From 9a7a10c37eff2dd594de7d65a785b80a3a6df45f Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Wed, 21 Dec 2016 17:41:39 +0100 Subject: [PATCH 13/38] now S2 is working --- src/ex1.c | 193 +++++++++++++++++++++++++++++------------------------- 1 file changed, 104 insertions(+), 89 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index c0c7b23..c255bec 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -243,29 +243,39 @@ int main(int argc,char **argv) ierr = VecGetValues(xr, nlocal, indxr, valxr);CHKERRQ(ierr); for (ii=Istart; ii Date: Wed, 21 Dec 2016 17:42:56 +0100 Subject: [PATCH 14/38] removed redundant printing --- src/ex1.c | 20 ++------------------ 1 file changed, 2 insertions(+), 18 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index c255bec..de07d02 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -243,29 +243,16 @@ int main(int argc,char **argv) ierr = VecGetValues(xr, nlocal, indxr, valxr);CHKERRQ(ierr); for (ii=Istart; ii Date: Wed, 21 Dec 2016 17:49:39 +0100 Subject: [PATCH 15/38] printing energies and s2 together --- src/ex1.c | 63 +++++++++++++++++++++++-------------------------------- 1 file changed, 26 insertions(+), 37 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index de07d02..284ae91 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -165,41 +165,6 @@ int main(int argc,char **argv) ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr); //ierr = EPSPrintSolution(eps,NULL);CHKERRQ(ierr); - if (nconv>0) { - /* - Display eigenvalues && relative errors - */ - ierr = PetscPrintf(PETSC_COMM_WORLD, - " k ||Ax-kx||/||kx||\n" - " ----------------- ------------------\n");CHKERRQ(ierr); - - for (i=0;i0) { + ierr = PetscPrintf(PETSC_COMM_WORLD, + " k ||Ax-kx||/||kx|| \n" + " ----------------- ----------------- ------------------\n");CHKERRQ(ierr); + for(i=0;i Date: Wed, 21 Dec 2016 18:43:41 +0100 Subject: [PATCH 16/38] working on density matrix --- src/ex1.c | 65 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/src/ex1.c b/src/ex1.c index 284ae91..a94ab9a 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -28,6 +28,8 @@ int main(int argc,char **argv) int natomax=700, iaa2, iaa; int ideter[natomax]; int ideter2[natomax]; + PetscScalar densmat[4][4]={0.0}; + PetscScalar trace; char const* const fileName = argv[1]; FILE* file = fopen(fileName, "r"); @@ -347,6 +349,67 @@ int main(int argc,char **argv) xymatfin = 0.0; normfin = 0.0; + /* + * Calculating the one-particle density matrix + */ + + for(kko=0;kko<=3;kko++){ + for(kok=0;kok<=3;kok++){ + + for(ii=Istart;ii<=Iend;ii+=1){ + + iii = ii+1; + getdet_(&iii,ideter); + for(kkio=0;kkio<=7;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + if(ideter[kko] != 3){ + if(ideter[kok] == 3){ + ideter2[kok]=ideter[kko]; + ideter2[kko]=3; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + Vec xiaa2; /* initial vector, destination vector */ + VecScatter scatter; /* scatter context */ + IS from, to; /* index sets that define the scatter */ + PetscScalar *values; + int idx_from[] = {iaa2}, idx_to[] = {0}; + VecCreateSeq(PETSC_COMM_SELF,1,&xiaa2); + ISCreateGeneral(PETSC_COMM_SELF,1,idx_from,PETSC_COPY_VALUES,&from); + ISCreateGeneral(PETSC_COMM_SELF,1,idx_to, PETSC_COPY_VALUES,&to); + VecScatterCreate(xr,from,xiaa2,to,&scatter); + VecScatterBegin(scatter,xr,xiaa2,INSERT_VALUES,SCATTER_FORWARD); + VecScatterEnd(scatter,xr,xiaa2, INSERT_VALUES,SCATTER_FORWARD); + VecGetArray(xiaa2,&values); + ISDestroy(&from); + ISDestroy(&to); + VecScatterDestroy(&scatter); + densmat[kko][kok]=densmat[kko][kok]+valxr[ii]*values[0]; + if(1)printf("id = %d iaa2 = %d valxr = %18f\n", mpiid, iaa2, values[0]); + } + } + if(kko == kok && ideter[kko] != 3){ + densmat[kko][kko]=densmat[kko][kko]+valxr[ii]*valxr[ii]; + } + + } + + } + } + + trace = 0.0; + if(!mpiid){ + for(kko=0;kko<=3;kko++){ + trace = trace + densmat[kko][kko]; + } + for(kko=0;kko<=3;kko++){ + for(kok=0;kok<=3;kok++){ + densmat[kko][kok]=0.0; + } + } + } + + /* Compute the relative error associated to each eigenpair */ @@ -362,7 +425,7 @@ int main(int argc,char **argv) 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 %18f\n",(double)re,(double)error,(double)XS);CHKERRQ(ierr); + ierr = PetscPrintf(PETSC_COMM_WORLD," %18f %12g %18f %18f\n",(double)re,(double)error,(double)XS,(double)trace);CHKERRQ(ierr); } } ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr); From 934b7494e0aa893ff7a941b78dc2810ead5b7d40 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Wed, 21 Dec 2016 19:09:11 +0100 Subject: [PATCH 17/38] looks like its working with 1 proc --- src/ex1.c | 33 +++++++++++++++++++++++---------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index a94ab9a..cae2605 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -29,6 +29,7 @@ int main(int argc,char **argv) int ideter[natomax]; int ideter2[natomax]; PetscScalar densmat[4][4]={0.0}; + PetscScalar densmatfin[4][4]={0.0}; PetscScalar trace; char const* const fileName = argv[1]; @@ -57,6 +58,11 @@ int main(int argc,char **argv) PetscReal xymatfin=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}; SlepcInitialize(&argc,&argv,(char*)0,NULL); ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr); @@ -353,10 +359,11 @@ int main(int argc,char **argv) * Calculating the one-particle density matrix */ + for(ii=Istart;ii<=Iend;ii+=1){ + for(kko=0;kko<=3;kko++){ for(kok=0;kok<=3;kok++){ - for(ii=Istart;ii<=Iend;ii+=1){ iii = ii+1; getdet_(&iii,ideter); @@ -369,21 +376,23 @@ int main(int argc,char **argv) ideter2[kko]=3; adr_(ideter2, &iaa2); iaa2 = iaa2 - 1; - Vec xiaa2; /* initial vector, destination vector */ - VecScatter scatter; /* scatter context */ - IS from, to; /* index sets that define the scatter */ - PetscScalar *values; - int idx_from[] = {iaa2}, idx_to[] = {0}; + + /* get value from other processor */ + + idx_from[0] = iaa2; + IS from, to; /* index sets that define the scatter */ VecCreateSeq(PETSC_COMM_SELF,1,&xiaa2); ISCreateGeneral(PETSC_COMM_SELF,1,idx_from,PETSC_COPY_VALUES,&from); ISCreateGeneral(PETSC_COMM_SELF,1,idx_to, PETSC_COPY_VALUES,&to); VecScatterCreate(xr,from,xiaa2,to,&scatter); VecScatterBegin(scatter,xr,xiaa2,INSERT_VALUES,SCATTER_FORWARD); - VecScatterEnd(scatter,xr,xiaa2, INSERT_VALUES,SCATTER_FORWARD); + VecScatterEnd (scatter,xr,xiaa2,INSERT_VALUES,SCATTER_FORWARD); VecGetArray(xiaa2,&values); ISDestroy(&from); ISDestroy(&to); VecScatterDestroy(&scatter); + /* value stored in values */ + densmat[kko][kok]=densmat[kko][kok]+valxr[ii]*values[0]; if(1)printf("id = %d iaa2 = %d valxr = %18f\n", mpiid, iaa2, values[0]); } @@ -392,19 +401,23 @@ int main(int argc,char **argv) densmat[kko][kko]=densmat[kko][kko]+valxr[ii]*valxr[ii]; } - } - } } + printf("mpiid = %d ii = %d kko = %d kok = %d\n", mpiid, ii, kko, kok); + } + + printf("done----done"); trace = 0.0; + MPI_Reduce(densmat, densmatfin, 16, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD); if(!mpiid){ for(kko=0;kko<=3;kko++){ - trace = trace + densmat[kko][kko]; + trace = trace + densmatfin[kko][kko]; } for(kko=0;kko<=3;kko++){ for(kok=0;kok<=3;kok++){ densmat[kko][kok]=0.0; + densmatfin[kko][kok]=0.0; } } } From 4e5e4f36ecdf7067f335bb42f905229d29d81db4 Mon Sep 17 00:00:00 2001 From: Bouammali-MA <35452795+Bouammali-MA@users.noreply.github.com> Date: Mon, 15 Jan 2018 12:31:03 +0100 Subject: [PATCH 18/38] Makefile Update added in line 27 IRPF90_temp/ knowing that the requsted file should be in that directory. --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 99f2424..5982ec1 100644 --- a/Makefile +++ b/Makefile @@ -24,7 +24,7 @@ ${BIN_DIR}: directories: ${OBJ_DIR} ${LIB_DIR} ${BIN_DIR} ${LIB_DIR}/irpf90.a: directories - cd ${SRC_DIR} && irpf90 init && $(MAKE) irpf90.a && cp irpf90.a ../${LIB_DIR} + cd ${SRC_DIR} && irpf90 init && $(MAKE) irpf90.a && cp IRPF90_temp/irpf90.a ../${LIB_DIR} ${OBJ_DIR}/read2.o: ${SRC_DIR}/read2.c directories chkopts ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} From daa2efc25479e78b1e28cbe7fc7d1f168791749f Mon Sep 17 00:00:00 2001 From: Bouammali-MA <35452795+Bouammali-MA@users.noreply.github.com> Date: Mon, 15 Jan 2018 12:37:44 +0100 Subject: [PATCH 19/38] Update of Read.md file Addition of two lines 25 and 26. for solving C leading files path. and comment during export for process --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 186955e..389a955 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,8 @@ Double Exchange Hamiltonian: Complete Version ```shell export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION} export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION} + export C_INCLUDE_PATH+=:$PETSC_DIR/include/:$SLEPC_DIR/include:$PETSC_DIR/arch-linux2-c-debug/include/:$SLEPC_DIR/arch-linux2-c-debug/include + # The "arch-linux2-c-debug" directory can have different names depending on PETSC and SLEPC installation procedure. ``` 2. Make the executable From 31c7aa26946673afbfeb1ae648594d9d16998810 Mon Sep 17 00:00:00 2001 From: vijay Date: Mon, 15 Jan 2018 12:41:40 +0100 Subject: [PATCH 20/38] export path for c compiler --- README.md | 47 ++++++++++++++++++++++++----------------------- 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/README.md b/README.md index 186955e..38d4ff9 100644 --- a/README.md +++ b/README.md @@ -19,16 +19,17 @@ Double Exchange Hamiltonian: Complete Version 1. Export environment variables for PETSc and SLEPc - ```shell - export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION} - export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION} - ``` +```shell +export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION} +export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION} +export C_INCLUDE_PATH+=:$PETSC_DIR/include/:$SLEPC_DIR/include:$PETSC_DIR/arch-linux2-c-debug/include/:$SLEPC_DIR/arch-linux2-c-debug/include +``` 2. Make the executable - ```shell - make ex1 - ``` +```shell +make ex1 +``` 3. Using DEHam --------------- @@ -37,25 +38,25 @@ Double Exchange Hamiltonian: Complete Version has the topology of the Hamiltonian and the various parameters as explained below in a sample inputfile: - ```python - 140 # The total number of determinants - 7 # The largest number of non-zero elements per row - 2 # The number of processors used in parallel - 1 # The number of holes - 0 # The isz (ms-1/2) value - 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) - .1430,-0.20,0.0000 # The three types of links this line gives J, K - .1430,-0.20,0.0000 # - -1.00,0.0,0.00 # This line gives t - ``` +```python +140 # The total number of determinants +7 # The largest number of non-zero elements per row +2 # The number of processors used in parallel +1 # The number of holes +0 # The isz (ms-1/2) value +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) +.1430,-0.20,0.0000 # The three types of links this line gives J, K +.1430,-0.20,0.0000 # +-1.00,0.0,0.00 # This line gives t +``` 2. running DEHam - ```shell - mpiexec -n [nprocs] ./ex1 inpfile - ``` +```shell +mpiexec -n [nprocs] ./ex1 inpfile +``` 4. Publications using this code ------------------------------- From 5357361c1237dfabd9086ed2d6cb48075c99ec21 Mon Sep 17 00:00:00 2001 From: vijay Date: Mon, 15 Jan 2018 12:47:44 +0100 Subject: [PATCH 21/38] cosmetics --- README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 38d4ff9..9121fe5 100644 --- a/README.md +++ b/README.md @@ -7,14 +7,14 @@ Double Exchange Hamiltonian: Complete Version (under GNU GENERAL PUBLIC LICENSE v2) -1. Dependencies +_Dependencies_ --------------- 1. [PETSc](https://www.mcs.anl.gov/petsc/documentation/installation.html) and [SLEPc](http://slepc.upv.es/documentation/current/docs/instal.htm) 2. [IRPF90](https://github.com/scemama/irpf90) -2. Compiling +_Compiling_ ------------ 1. Export environment variables for PETSc and SLEPc @@ -31,7 +31,7 @@ export C_INCLUDE_PATH+=:$PETSC_DIR/include/:$SLEPC_DIR/include:$PETSC_DIR/arch-l make ex1 ``` -3. Using DEHam +_Using DEHam_ --------------- 1. The DEHam program requires an input file which @@ -58,7 +58,7 @@ make ex1 mpiexec -n [nprocs] ./ex1 inpfile ``` -4. Publications using this code +_Publications using this code_ ------------------------------- 1. High-Spin Chains and Crowns from Double-Exchange Mechanism [doi:10.3390/cryst6040039](http://www.dx.doi.org/10.3390/cryst6040039) From f13054d927a6b392e8ba30c51daf38fc585f42ab Mon Sep 17 00:00:00 2001 From: vijay Date: Mon, 15 Jan 2018 16:43:52 +0100 Subject: [PATCH 22/38] removed tabs after merge --- README.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index e18f659..21b069f 100644 --- a/README.md +++ b/README.md @@ -19,12 +19,12 @@ _Compiling_ 1. Export environment variables for PETSc and SLEPc - ```shell - export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION} - export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION} - export C_INCLUDE_PATH+=:$PETSC_DIR/include/:$SLEPC_DIR/include:$PETSC_DIR/arch-linux2-c-debug/include/:$SLEPC_DIR/arch-linux2-c-debug/include - # The "arch-linux2-c-debug" directory can have different names depending on PETSC and SLEPC installation procedure. - ``` +```shell +export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION} +export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION} +export C_INCLUDE_PATH+=:$PETSC_DIR/include/:$SLEPC_DIR/include:$PETSC_DIR/arch-linux2-c-debug/include/:$SLEPC_DIR/arch-linux2-c-debug/include +# The "arch-linux2-c-debug" directory can have different names depending on PETSC and SLEPC installation procedure. +``` 2. Make the executable From 0526613506f64f6acf6a12c4ee336ea7d82dd3df Mon Sep 17 00:00:00 2001 From: vijay Date: Mon, 15 Jan 2018 23:30:38 +0100 Subject: [PATCH 23/38] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 21b069f..7f0690d 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ Double Exchange Hamiltonian: Complete Version _Dependencies_ --------------- - 1. [PETSc](https://www.mcs.anl.gov/petsc/documentation/installation.html) and [SLEPc](http://slepc.upv.es/documentation/current/docs/instal.htm) + 1. [PETSc](https://www.mcs.anl.gov/petsc/documentation/installation.html) and [SLEPc](http://slepc.upv.es/documentation/instal.htm) 2. [IRPF90](https://github.com/scemama/irpf90) From 2e6f7d646ecdaf15f846f144c8f166f058020381 Mon Sep 17 00:00:00 2001 From: vijay Date: Fri, 26 Jan 2018 13:12:12 +0100 Subject: [PATCH 24/38] add new input file features --- README.md | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 7f0690d..be390be 100644 --- a/README.md +++ b/README.md @@ -41,9 +41,10 @@ _Using DEHam_ as explained below in a sample inputfile: ```python -140 # The total number of determinants -7 # The largest number of non-zero elements per row -2 # The number of processors used in parallel +140 # The total number of determinants (Ndet) +7 # The largest number of non-zero elements per row (Multiple of Ndet) +1 # The number of matrix elements to calculate at one (Multiple of Ndet) +2 # The total number of processors used in parallel (Multiple of Ndet) 1 # The number of holes 0 # The isz (ms-1/2) value 1,2,3,1,2,3,4,5,6,7 # The topology of the system is specified here From 198aa1bf4f4d2266c94ce111b4d63b68f5bd1a36 Mon Sep 17 00:00:00 2001 From: vijay Date: Fri, 26 Jan 2018 13:12:50 +0100 Subject: [PATCH 25/38] cosmetics --- README.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index be390be..c8afe94 100644 --- a/README.md +++ b/README.md @@ -42,11 +42,11 @@ _Using DEHam_ ```python 140 # The total number of determinants (Ndet) -7 # The largest number of non-zero elements per row (Multiple of Ndet) -1 # The number of matrix elements to calculate at one (Multiple of Ndet) -2 # The total number of processors used in parallel (Multiple of Ndet) -1 # The number of holes -0 # The isz (ms-1/2) value +7 # The largest number of non-zero elements per row (Multiple of Ndet) +1 # The number of matrix elements to calculate at onece (Multiple of Ndet) +2 # The total number of processors used in parallel (Multiple of Ndet) +1 # The number of holes +0 # The isz (ms-1/2) value 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) From 4fe71bb862694959a758b9d823403ef4c0dd9a15 Mon Sep 17 00:00:00 2001 From: vijay Date: Fri, 26 Jan 2018 13:15:45 +0100 Subject: [PATCH 26/38] fix a spelling mistake --- README.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index c8afe94..1d4e314 100644 --- a/README.md +++ b/README.md @@ -41,12 +41,12 @@ _Using DEHam_ as explained below in a sample inputfile: ```python -140 # The total number of determinants (Ndet) -7 # The largest number of non-zero elements per row (Multiple of Ndet) -1 # The number of matrix elements to calculate at onece (Multiple of Ndet) -2 # The total number of processors used in parallel (Multiple of Ndet) -1 # The number of holes -0 # The isz (ms-1/2) value +140 # The total number of determinants (Ndet) +7 # The largest number of non-zero elements per row (Multiple of Ndet) +1 # The number of matrix elements to calculate at once (Multiple of Ndet) +2 # The total number of processors used in parallel (Multiple of Ndet) +1 # The number of holes +0 # The isz (ms-1/2) value 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) From d1d6eae91d63f936a61c66ff0d7f112783addfef Mon Sep 17 00:00:00 2001 From: vijay Date: Fri, 26 Jan 2018 13:38:55 +0100 Subject: [PATCH 27/38] mistake in description of inputfile --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 1d4e314..6170e0e 100644 --- a/README.md +++ b/README.md @@ -42,8 +42,8 @@ _Using DEHam_ ```python 140 # The total number of determinants (Ndet) -7 # The largest number of non-zero elements per row (Multiple of Ndet) -1 # The number of matrix elements to calculate at once (Multiple of Ndet) +7 # The number of orbitals (total) +1 # The largest number of non-zero elements per row (Multiple of Ndet) 2 # The total number of processors used in parallel (Multiple of Ndet) 1 # The number of holes 0 # The isz (ms-1/2) value From 3278aabfeb7bf73cc8f2788d2d1dfe067976535b Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Sat, 27 Jan 2018 12:41:48 +0100 Subject: [PATCH 28/38] updating code to the current local version, it might not compile atm. Many new features added: 1. getting S2 values 2. possibility of setting position of hole 3. possibility of setting Sbox 4. three Sbox definitions at once 5. Doing only FAM1 or the full set of states 6. efficiency improvements --- README.md | 1 + src/adr.irp.f | 16 +- src/adrnew.irp.f | 57 +++ src/analyse.irp.f | 832 ++++++++++++++++++++++++++++++++++++++++ src/conv.irp.f | 14 +- src/desort.irp.f | 10 +- src/ex1.c | 448 ++++++++++------------ src/extra_diag.irp.f | 2 +- src/get_dmat.c | 145 +++++++ src/get_dmat.h | 9 + src/get_s2.c | 683 +++++++++++++++++++++++++++++++++ src/get_s2.h | 50 +++ src/get_s2_cyclic.c | 784 +++++++++++++++++++++++++++++++++++++ src/get_s2_mov.c | 683 +++++++++++++++++++++++++++++++++ src/getdet.irp.f | 57 +-- src/natom.irp.f | 4 +- src/providdet.irp.f | 85 ++++ src/provide_clust.irp.f | 1 + src/read2.c | 79 +++- src/read2.h | 21 +- src/searchdet.irp.f | 110 +++--- src/sort.irp.f | 10 +- src/unit_l1.irp.f | 4 + 23 files changed, 3721 insertions(+), 384 deletions(-) create mode 100644 src/adrnew.irp.f create mode 100644 src/analyse.irp.f create mode 100644 src/get_dmat.c create mode 100644 src/get_dmat.h create mode 100644 src/get_s2.c create mode 100644 src/get_s2.h create mode 100644 src/get_s2_cyclic.c create mode 100644 src/get_s2_mov.c create mode 100644 src/providdet.irp.f diff --git a/README.md b/README.md index 6170e0e..48e92e9 100644 --- a/README.md +++ b/README.md @@ -47,6 +47,7 @@ _Using DEHam_ 2 # The total number of processors used in parallel (Multiple of Ndet) 1 # The number of holes 0 # The isz (ms-1/2) value +true # Restrict the hole to the 1'st (i.e. half of natom) Family of states. *false* for no restrictions 1,2,3,1,2,3,4,5,6,7 # The topology of the system is specified here 2,3,4,8,7,6,5,6,7,8 # first and second line contain the two sites linked 1,1,1,2,2,2,2,3,3,3 # third line contains the type of link (1 for t, J 2 for K and 3 for none) diff --git a/src/adr.irp.f b/src/adr.irp.f index 76db98d..2a9cd5c 100644 --- a/src/adr.irp.f +++ b/src/adr.irp.f @@ -8,24 +8,24 @@ subroutine adr(ideter,add) END_DOC integer,INTENT(INOUT)::ideter(natomax) integer(kind=selected_int_kind(16)),INTENT(INOUT)::add - integer(kind=selected_int_kind(16))::det,deth,addh,detnew + integer(kind=selected_int_kind(16))::deti,dethi,addh,detnew integer::count,i,j - det=0 + deti=0 detnew=0 - deth=0 + dethi=0 count=0 - call conv(ideter,det,deth) + call conv(ideter,deti,dethi) Do i=0,natom-1 - if(BTEST(deth,i))then + if(BTEST(dethi,i))then count=count+1 endif - if(BTEST(det,i))then + if(BTEST(deti,i))then detnew=IBSET(detnew,i-count) endif enddo - det=detnew - call searchdet(det,add,deth,addh) + deti=detnew + call searchdet(deti,add,dethi,addh) add = add + (nt1-addh)*(nt2) diff --git a/src/adrnew.irp.f b/src/adrnew.irp.f new file mode 100644 index 0000000..91ef8c0 --- /dev/null +++ b/src/adrnew.irp.f @@ -0,0 +1,57 @@ +subroutine adrfull() + implicit none + BEGIN_DOC + ! this subroutine provides the address of a detrminant + ! given in old format. + ! It searches in a list of generated determinants and + ! matches the given determinant. + END_DOC + integer,dimension(natomax)::ideter + integer(kind=selected_int_kind(16))::add + integer(kind=selected_int_kind(16))::deti,dethi,addh,detnew + integer::count,i,j + + deti=0 + detnew=0 + dethi=0 + do j=1,detfound + detnew=0 + count=0 + ideter=foundet(:,j) + call conv(ideter,deti,dethi) + Do i=0,natom-1 + if(BTEST(dethi,i))then + count=count+1 + endif + if(BTEST(deti,i))then + detnew=IBSET(detnew,i-count) + endif + enddo + deti=detnew + foundadd(j,1)=deti + foundadd(j,3)=j + foundaddh(j,1)=dethi + foundaddh(j,3)=j + call searchdet(deti,add,dethi,addh) +! enddo +! call sort() +! call searchdetfull() +! call desort() +! do i=1,detfound +! add = foundadd(i,2) +! addh = foundaddh(i,2) + foundadd(j,2) = add + foundaddh(j,2)= addh + add = add + (nt1-addh)*(nt2) + foundetadr(j)=add + enddo + + +10 FORMAT(B64,I8,F8.2) +15 FORMAT(B64,I8,I8,I8) +11 FORMAT(B64,I3,B64) +12 FORMAT(I5,$) +13 FORMAT(B64,B64) +14 FORMAT(B64,I14) +16 FORMAT(B64,I14,I14) +end diff --git a/src/analyse.irp.f b/src/analyse.irp.f new file mode 100644 index 0000000..c621bc2 --- /dev/null +++ b/src/analyse.irp.f @@ -0,0 +1,832 @@ + SUBROUTINE ANALYSE(vect, dimvect, startvect, endvect, xymat2, norm2) +! INCLUDE "nbtots.prm" + IMPLICIT NONE + INTEGER dimvect, nbtots, startvect, endvect + REAL*8,dimension(dimvect)::vect + INTEGER (kind=selected_int_kind(16))::add,kvect + INTEGER (kind=selected_int_kind(16))::iaa2,i + INTEGER ,dimension(natomax)::ideter + INTEGER ,dimension(natomax)::ideter2 + REAL*8,allocatable ::xz(:) + REAL*8::xmat,xymat + REAL*8::xmat1,xymat1 + REAL*8::xmat2 + REAL*8,INTENT(INOUT)::xymat2 + REAL*8::xmat3,xymat3 + REAL*8::sym,nonsym,proj_trou + REAL*8,allocatable ::xalpha1(:) +! REAL*8,allocatable ::vect(:) + REAL*8,INTENT(INOUT)::norm2 + REAL*8::norm,norm1,norm3,proj_2trou + REAL*8::t1,t2,XS,XS1,XS2,XS3 + REAL*8::resta_mono,resta_one,resta_bi,delta + INTEGER ::kko,kok,kkio,j,eigen,nigen,count + INTEGER ::cntrou,countlvect,ndim,iaa + INTEGER ::ipt_1,ipt_2,iptemp_1,iptemp_2 + INTEGER ::ipt_3,iptemp_3 + INTEGER ::ibougetrou,jstart + REAL*8 ,allocatable::eigenvectors(:,:) + REAL*8 ,allocatable::eigenvalues(:) + REAL*8 ,allocatable::WORK(:) + REAL*8 ,allocatable::AP(:) + REAL*8 ,allocatable::densmat(:,:) + REAL*8 ,allocatable::densmat2(:,:) + REAL*8 proj_1,extradiag_dmat2,ionic,nonionic + REAL*8 proj_2,sum,conduction,prob,prob2 + INTEGER INFO,nrow,ncol,mmo,mom,kk,k,omm,okk + CHARACTER*1 JOBZ,UPLO + INTEGER::RESTA=0 + +! allocate(vect(nbtots)) + allocate(xalpha1(natomax)) + allocate(xz (natom/2)) + +! OPEN (unit=59,file='FIL1',form='formatted',status='old') +! OPEN (unit=217,file='SBOX217',form='formatted',status='REPLACE') +! REWIND 59 +! READ (59,*) + +! print *,' in analyse', startvect, endvect +! print *,'nalpha=',nalpha,'nbeta=',nbeta +! PRINT *,natom,ntrou,nbtots,nt1,nt2,isz +! ndim=3 + ndim=(natom/2)*((natom/2)-1)/2 + allocate(AP((ndim)*((ndim)+1)/2)) + allocate(WORK(3*(ndim))) + allocate(eigenvectors((ndim)*(ndim),1)) + allocate(eigenvalues((ndim))) + allocate(densmat(ndim,ndim)) + allocate(densmat2(ndim,ndim)) +! Touch isz maxdet maxial maxlien maxplac nalpha natom natomax nbeta nbtots nt1 nt2 ntrou +! PRINT *,(vect(j),j=1,30) + + IF(RESTA .eq. 1)THEN + do i=1,natom/2 + if(mod(natom/2,2).eq.0)then + xz(i)=(((natom/2)/2)-0.5d0)-(i-1.0d0) + write(6 ,*)i,xz(i) + else + xz(i)=((natom/2)-1.0d0)/2.0d0-(i-1.0d0) + write(6 ,*)i,xz(i) + endif + enddo + ENDIF + +!! PROVIDE det deth + nigen=1 + DO eigen=1,nigen +! READ (59,10) (vect(j),j=1,nbtots) + IF (ntrou.eq.1) THEN + norm=0.d0 + norm1=0.d0 + norm2=0.d0 + norm3=0.d0 + count=0 + cntrou=0 + iaa=0 + iaa2=0 + countlvect=0 + proj_2trou=0.d0 + resta_bi=0.d0 + resta_mono=0.d0 + resta_one=0.d0 + xymat = 0.0d0 + xymat1 = 0.0d0 + xymat2 = 0.d0 + xymat3 = 0.d0 + proj_1=0d0 + proj_2=0d0 + densmat=0d0 + densmat2=0d0 + extradiag_dmat2=0d0 + conduction=0d0 + nrow=0 + ncol=0 + + jstart=1 + ipt_1=0 + ipt_2=0 + ipt_3=0 + iptemp_1=0 + iptemp_2=0 + iptemp_3=0 + nonionic=0.d0 + ionic=0.d0 + prob=0.d0 + prob2=0.d0 + sym=0.d0 + nonsym=0.d0 + ibougetrou=0 + + DO kvect=1, endvect-startvect + CALL getdet(kvect+startvect,ideter) + +!!---------------------------------------- +!! RESTA +!!---------------------------------------- +!!! mono +! proj_trou=vect(kvect)**2 +! DO i=1,natom/2 +! IF (ideter(i).eq.3) THEN +! delta=0.0d0 +! ELSE +! delta=1.0d0 +! END IF +! resta_mono=resta_mono+delta*xz(i)*xz(i)*proj_trou +! resta_one=resta_one+delta*xz(i)*proj_trou +! END DO +!!! bi +! DO i=1,natom/2 +! DO j=1,natom/2 +! IF (ideter(i).eq.3.or.ideter(j).eq.3.or.i.eq.j) & +! THEN +! delta=0.0d0 +! ELSE +! delta=1.0d0 +! END IF +! resta_bi=resta_bi+delta*xz(i)*xz(j)*proj_trou +! END DO +! END DO +!!---------------------------------------- + + +!!---------------------------------------- +!! Prob ionic non-ionic +!!---------------------------------------- +! ipt_1=0 +! ipt_2=0 +! ipt_3=0 +! DO kko=1,3 +! IF(ideter(kko).eq.3)THEN +! ipt_1=ipt_1+1 +! ENDIF +! ENDDO + +! IF(ipt_1.eq.1)THEN +! DO kko=4,6 +! IF(ideter(kko).eq.3)THEN +! ipt_2=ipt_2+1 +! ENDIF +! ENDDO +! ENDIF + +! IF(ipt_2.eq.1)THEN +! DO kko=7,9 +! IF(ideter(kko).eq.3)THEN +! ipt_3=ipt_3+1 +! ENDIF +! ENDDO +! ENDIF + +! IF(ipt_3 .eq. 1)THEN +! nonionic=nonionic+vect(kvect)**2 +! ELSE +! ionic=ionic+vect(kvect)**2 +! ENDIF +!!---------------------------------------- +!! S_box +!!---------------------------------------- + + xmat=0.0d0 + xmat1=0.0d0 + xmat2=0.0d0 +! IF (.TRUE.)THEN +!! IF (ideter(6).eq.3 ) THEN +!! norm=norm+vect(kvect)**2 +!! DO kko=5,7 +!! DO kok=kko,7 +!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN +!! xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2) +!! ELSE +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN +!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN +!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN +!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=2 +!! ideter2(kok)=1 +!! CALL adr(ideter2, iaa2) +!! xmat=xmat+vect(kvect)*vect(iaa2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN +!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=1 +!! ideter2(kok)=2 +!! CALL adr(ideter2, iaa2) +!! xmat=xmat+vect(kvect)*vect(iaa2) +!! END IF +!! END IF +!! END DO +!! END DO +!! +!! DO kko=16,18 +!! DO kok=kko,18 +!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN +!! xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2) +!! ELSE +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN +!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN +!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN +!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=2 +!! ideter2(kok)=1 +!! CALL adr(ideter2, iaa2) +!! xmat=xmat+vect(kvect)*vect(iaa2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN +!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=1 +!! ideter2(kok)=2 +!! CALL adr(ideter2, iaa2) +!! xmat=xmat+vect(kvect)*vect(iaa2) +!! END IF +!! END IF +!! END DO +!! END DO +!! +!! DO kko=5,7 +!! DO kok=16,18 +!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN +!! xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2) +!! ELSE +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN +!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN +!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN +!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=2 +!! ideter2(kok)=1 +!! CALL adr(ideter2, iaa2) +!! xmat=xmat+vect(kvect)*vect(iaa2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN +!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=1 +!! ideter2(kok)=2 +!! CALL adr(ideter2, iaa2) +!! xmat=xmat+vect(kvect)*vect(iaa2) +!! END IF +!! END IF +!! END DO +!! END DO +!! END IF +!!!---------------------------------------- +!! xymat=xymat+xmat +!! +!! IF (ideter(7).eq.3 ) THEN +!! norm1=norm1+vect(kvect)**2 +!! DO kko=5,9 +!! DO kok=kko,9 +!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN +!! xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2) +!! ELSE +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN +!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN +!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN +!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=2 +!! ideter2(kok)=1 +!! CALL adr(ideter2, iaa2) +!! xmat1=xmat1+vect(kvect)*vect(iaa2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN +!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=1 +!! ideter2(kok)=2 +!! CALL adr(ideter2, iaa2) +!! xmat1=xmat1+vect(kvect)*vect(iaa2) +!! END IF +!! END IF +!! END DO +!! END DO +!! +!! DO kko=14,18 +!! DO kok=kko,18 +!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN +!! xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2) +!! ELSE +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN +!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN +!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN +!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=2 +!! ideter2(kok)=1 +!! CALL adr(ideter2, iaa2) +!! xmat1=xmat1+vect(kvect)*vect(iaa2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN +!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=1 +!! ideter2(kok)=2 +!! CALL adr(ideter2, iaa2) +!! xmat1=xmat1+vect(kvect)*vect(iaa2) +!! END IF +!! END IF +!! END DO +!! END DO +!! +!! DO kko=5,9 +!! DO kok=14,18 +!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN +!! xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2) +!! ELSE +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN +!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN +!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2) +!! END IF +!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN +!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=2 +!! ideter2(kok)=1 +!! CALL adr(ideter2, iaa2) +!! xmat1=xmat1+vect(kvect)*vect(iaa2) +!! END IF +!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN +!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2) +!! DO kkio=1,natom +!! ideter2(kkio)=ideter(kkio) +!! END DO +!! ideter2(kko)=1 +!! ideter2(kok)=2 +!! CALL adr(ideter2, iaa2) +!! xmat1=xmat1+vect(kvect)*vect(iaa2) +!! END IF +!! END IF +!! END DO +!! END DO +!! END IF +!!!---------------------------------------- +!! xymat1=xymat1+xmat1 + + IF (.TRUE.)THEN + norm2=norm2+vect(kvect)**2 + DO kko=1,natom/2 + DO kok=kko,natom/2 + IF (kok.eq.kko.and.ideter(kok).ne.3) THEN + xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2) + ELSE + IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN + xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2) + END IF + IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN + xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2) + END IF + IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN + xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2) + DO kkio=1,natom + ideter2(kkio)=ideter(kkio) + END DO + ideter2(kko)=2 + ideter2(kok)=1 + CALL adr(ideter2, iaa2) + xmat2=xmat2+vect(kvect)*vect(iaa2-startvect) + END IF + IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN + xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2) + DO kkio=1,natom + ideter2(kkio)=ideter(kkio) + END DO + ideter2(kko)=1 + ideter2(kok)=2 + CALL adr(ideter2, iaa2) + xmat2=xmat2+vect(kvect)*vect(iaa2-startvect) + END IF + END IF + END DO + END DO + + DO kko=(natom/2)+1,natom + DO kok=kko,natom + IF (kok.eq.kko.and.ideter(kok).ne.3) THEN + xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2) + ELSE + IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN + xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2) + END IF + IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN + xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2) + END IF + IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN + xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2) + DO kkio=1,natom + ideter2(kkio)=ideter(kkio) + END DO + ideter2(kko)=2 + ideter2(kok)=1 + CALL adr(ideter2, iaa2) + xmat2=xmat2+vect(kvect)*vect(iaa2-startvect) + END IF + IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN + xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2) + DO kkio=1,natom + ideter2(kkio)=ideter(kkio) + END DO + ideter2(kko)=1 + ideter2(kok)=2 + CALL adr(ideter2, iaa2) + xmat2=xmat2+vect(kvect)*vect(iaa2-startvect) + END IF + END IF + END DO + END DO + + DO kko=1,natom/2 + DO kok=(natom/2)+1,natom + IF (kok.eq.kko.and.ideter(kok).ne.3) THEN + xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2) + ELSE + IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN + xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2) + END IF + IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN + xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2) + END IF + IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN + xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2) + DO kkio=1,natom + ideter2(kkio)=ideter(kkio) + END DO + ideter2(kko)=2 + ideter2(kok)=1 + CALL adr(ideter2, iaa2) + xmat2=xmat2+vect(kvect)*vect(iaa2-startvect) + END IF + IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN + xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2) + DO kkio=1,natom + ideter2(kkio)=ideter(kkio) + END DO + ideter2(kko)=1 + ideter2(kok)=2 + CALL adr(ideter2, iaa2) + xmat2=xmat2+vect(kvect)*vect(iaa2-startvect) + END IF + END IF + END DO + END DO + END IF +!---------------------------------------- +! print *,"norm = ",norm2,"xmat2 = ",xmat2,"vect =",vect(kvect),"natom=",natom + xymat2=xymat2+xmat2 + +! XS=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat/norm))) +! XS1=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat1/norm1))) + XS2=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat2/norm2))) + END DO +! print *,eigen,xymat2,XS2 + +!---------------------------------------- +! Resta and probabilities of dets +!---------------------------------------- +! DO kvect=1,nbtots +! ideter2=ideter +! CALL getdet(kvect,ideter) +! +! IF (jstart.eq.1) THEN +! jstart=0 +! ideter2=ideter +! DO kko=1,natom +! IF (ideter(kko).eq.3) THEN +! ipt_1=kko +! DO kok=kko+1,natom +! IF (ideter(kok).eq.3) THEN +! ipt_2=kok +! DO okk=kok+1,natom +! IF (ideter(okk).eq.3) THEN +! ipt_3=okk +! EXIT +! END IF +! END DO +! EXIT +! END IF +! END DO +! EXIT +! END IF +! END DO +! END IF +! +! +! DO kko=1,natom +! IF (ideter(kko).eq.3) THEN +! iptemp_1=kko +! DO kok=kko+1,natom +! IF (ideter(kok).eq.3) THEN +! iptemp_2=kok +! DO okk=kok+1,natom +! IF (ideter(okk).eq.3) THEN +! iptemp_3=okk +! EXIT +! END IF +! END DO +! EXIT +! END IF +! END DO +! EXIT +! END IF +! END DO +! +! IF (iptemp_1.ne.ipt_1.or.iptemp_2.ne.ipt_2 & +! .or. iptemp_3.ne.ipt_3) THEN +! ibougetrou=1 +! ipt_1=iptemp_1 +! ipt_2=iptemp_2 +! ipt_3=iptemp_3 +! ELSE +! proj_trou=proj_trou+vect(kvect)**2 +! ibougetrou=0 +! END IF +! +! IF (iptemp_1.eq.(9-iptemp_3+1) & +! .and. iptemp_2 .eq. 5 .and. iptemp_1.ne.4) THEN +! sym=sym+vect(kvect)**2 +! ELSEIF (ideter(4).eq.3 .and. ideter(5).eq.3 & +! .and. ideter(6).eq.3) THEN +! nonsym=nonsym+vect(kvect)**2 +! ELSE +!! ELSEIF(ideter(2).eq.3 .and. ideter(3).eq.3 & +!! .and. ideter(6).eq.3)then +! nonsym=nonsym+vect(kvect)**2 +! END IF +! +! IF (ideter(2).eq.3 .and. ideter(4).eq.3 & +! .and. ideter(8).eq.3) THEN +! prob=prob+vect(kvect)**2 +! END IF +! +!! IF (ideter(1).eq.3 .and. ideter(2).eq.3 & +!! .and. ideter(3).eq.3) THEN +!! prob2=prob2+vect(kvect)**2 +!! END IF +! +! IF (ibougetrou.eq.1.or.kvect.eq.nbtots) THEN +!!---------------------------------------- +!! mono +! DO i=1,natom/2 +! IF (ideter2(i).eq.3) THEN +! delta=1.0d0 +! ELSE +! delta=0.0d0 +! END IF +! resta_mono=resta_mono+delta*xz(i)*xz(i)*proj_trou +! resta_one=resta_one+delta*xz(i)*proj_trou +! END DO +!! bi +! DO i=1,natom/2 +! DO j=1,natom/2 +! IF (ideter2(i).ne.3.or.ideter2(j).ne.3.or.i.eq.j) & +! THEN +! delta=0.0d0 +! ELSE +! delta=1.0d0 +! END IF +! resta_bi=resta_bi+delta*xz(i)*xz(j)*proj_trou +! END DO +! END DO +!!---------------------------------------- +! proj_trou=0.d0 +! proj_trou=proj_trou+vect(kvect)**2 +! END IF +! END DO + +!---------------------------------------- +! One particle density matrix +!---------------------------------------- +! DO kko=1,3 +! DO kok=1,3 +! +! DO kvect=1,nbtots + +! CALL getdet(kvect,ideter) +! ideter2=ideter +! IF (ideter(kko).ne.3) THEN +! IF (ideter(kok).eq.3) THEN +! ideter2(kok)=ideter(kko) +! ideter2(kko)=3 +! CALL adr(ideter2, iaa2) +! densmat(kko,kok)=densmat(kko,kok)+vect(kvect)* & +! vect(iaa2) +! END IF +! END IF +! IF (kko.eq.kok.and.ideter(kko).ne.3) THEN +! densmat(kko,kko)=densmat(kko,kko)+vect(kvect)**2 +! END IF + +! END DO + +! END DO +! END DO + +!---------------------------------------- +! two particle density matrix +!---------------------------------------- + +! DO kko=1,(natom/2)-1 +! DO kok=kko+1,natom/2 + +! nrow=nrow+1 +! ncol=0 +! DO mmo=1,(natom/2)-1 +! DO mom=mmo+1,natom/2 + +! ncol=ncol+1 + +! DO kvect=1,nbtots +! CALL getdet(kvect,ideter) +! ideter2=ideter +! IF (ideter(kko).eq.3.and.ideter(kok) & +! .eq.3.and.ideter(mmo).ne.3.and.ideter(mom).ne.3) & +! THEN +! if(ideter(kok).ne.3 .and. ideter(mom).ne.3)then +! ideter2(kko)=ideter(mmo) +! ideter2(mmo)=3 +! ideter2(kok)=ideter(mom) +! ideter2(mom)=3 +! CALL adr(ideter2, iaa2) +! densmat2(nrow,ncol)=densmat2(nrow,ncol)+ & +! vect(kvect)*vect(iaa2) +! print *,nrow,ncol,kko,kok,mmo,mom, +! * densmat2(nrow,ncol) +! endif +! END IF + + +! IF (nrow.eq.ncol.and.ideter(mmo) & +! .ne.3.and.ideter(mom).ne.3) THEN +! densmat2(nrow,ncol)=densmat2(nrow,ncol)+ & +! vect(kvect)**2 + +! END IF + +! END DO + + +! END DO +! END DO + +! END DO +! END DO + +!---------------------------------------- + +!---------------------------------------- +! conduction +!---------------------------------------- + +! count=0 +! DO kko=1,(natom/2)-2 +! DO kok=kko+1,(natom/2)-1 +! DO okk=kok+1,natom/2 +! +! nrow=nrow+1 +! ncol=0 +! DO mmo=kko,kko+1 +! DO mom=kok,kok+1 +! DO omm=okk,okk+1 +! +! ncol=ncol+1 +! DO kvect=1,nbtots +! CALL getdet(kvect,ideter) +! ideter2=ideter +! IF (abs(kko-mmo).eq.1.or.abs(kok-mom).eq.1 & +! .or. abs(okk-omm).eq.1) THEN +! IF (mmo.le.natom/2.and.mom.le.natom/2 .and. & +! omm.le.natom/2) THEN +! IF (mmo.ne.mom .and. mom.ne.omm) THEN +! IF (ideter(kko).eq.3 .and. ideter(kok).eq.3 & +! .and. ideter(okk).eq.3) THEN +! ideter2(okk)=ideter2(omm) +! ideter2(omm)=3 +! ideter2(kok)=ideter2(mom) +! ideter2(mom)=3 +! ideter2(kko)=ideter2(mmo) +! ideter2(mmo)=3 +! CALL adr(ideter2, iaa2) +!! count=0 +!! do i=1,natom/2 +!! if(ideter2(i).eq.3)then +!! count+=1 +!! endif +!! enddo +!! print *,kko,kok,okk,mmo,mom,omm,iaa2 +! conduction=conduction+dabs(vect(kvect)*vect(iaa2)) +! END IF +! END IF +! END IF +! END IF +! END DO +! +! END DO +! END DO +! END DO +! +! END DO +! END DO +! END DO + +!---------------------------------------- + +! DO j=1,ndim +! write(217,1022)j,(densmat(j,kko),kko=1,ndim) +! END DO +!---------------------------------------- +! diagonalisation de mat +! affiche vecteur +! JOBZ='V' +! matrice sup +! UPLO='U' + +! matrice en vecteur ligne ... +! extradiag_dmat2=0d0 +! k=0 +! DO j=1,ndim +! DO i=1,j-1 +! if(i.ne.j)then +! extradiag_dmat2 = extradiag_dmat2 + dabs(densmat2(i,j)) +! endif +! END DO +! END DO + +! appel subroutine LAPACK de diagonalisation :: double précision !! +! INFO=0 +! CALL DSPEV (JOBZ, UPLO, ndim, AP, eigenvalues, eigenvectors, & +! ndim, WORK, INFO) + +! IF (INFO.ne.0) THEN +! PRINT *,'SUBROUTINE MATRIX: Error at dspev',info +! CALL exit (1) +! END IF + +! proj_2=0.d0 +! sum=0d0 +! DO j=1,ndim +! proj_2=proj_2-eigenvalues(j)*log(eigenvalues(j)) +! sum+=eigenvalues(j) +! write(214,*)eigenvalues(j) +! END DO + +! XS=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat/norm))) +! XS2=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat2/norm2))) +! XS3=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat3/norm3))) +! WRITE (217,*) eigen,XS,norm + END IF + + END DO + + + + 10 FORMAT (E25.0) + 1022 FORMAT(3x,I3,6(2x,F12.4)) + END diff --git a/src/conv.irp.f b/src/conv.irp.f index f296c60..e01b839 100644 --- a/src/conv.irp.f +++ b/src/conv.irp.f @@ -1,20 +1,20 @@ -subroutine conv(ideter,det,deth) +subroutine conv(ideter,deti,dethi) implicit none BEGIN_DOC ! this routine converts a detrminant in the old ! format into the new one and returns the determinant. END_DOC integer,INTENT(INOUT)::ideter(natomax) - integer(kind=selected_int_kind(16)),INTENT(INOUT)::det - integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth + integer(kind=selected_int_kind(16)),INTENT(INOUT)::deti + integer(kind=selected_int_kind(16)),INTENT(INOUT)::dethi integer::i - det=0 - deth=0 + deti=0 + dethi=0 do i=1,natom if(ideter(natom-i+1).eq.2 .and. ideter(natom-i+1).ne.3)then - det=IBSET(det,i-1) + deti=IBSET(deti,i-1) elseif(ideter(natom-i+1).eq.3)then - deth=IBSET(deth,i-1) + dethi=IBSET(dethi,i-1) endif enddo end diff --git a/src/desort.irp.f b/src/desort.irp.f index c4d09b1..643a7ad 100644 --- a/src/desort.irp.f +++ b/src/desort.irp.f @@ -1,14 +1,14 @@ subroutine desort() implicit none integer::i,j,ord,ordh - integer(kind=selected_int_kind(16))::add,addh,det,deth,addt + integer(kind=selected_int_kind(16))::add,addh,deti,dethi,addt do i=1,detfound-1 do j=i+1,detfound if(foundaddh(i,3).gt.foundaddh(j,3))then - deth = foundaddh(i,1) + dethi = foundaddh(i,1) foundaddh(i,1) = foundaddh(j,1) - foundaddh(j,1) = deth + foundaddh(j,1) = dethi addh = foundaddh(i,2) foundaddh(i,2) = foundaddh(j,2) foundaddh(j,2) = addh @@ -17,9 +17,9 @@ subroutine desort() foundaddh(j,3) = ordh endif if(foundadd(i,3).gt.foundadd(j,3))then - det = foundadd(i,1) + deti = foundadd(i,1) foundadd(i,1) = foundadd(j,1) - foundadd(j,1) = det + foundadd(j,1) = deti add = foundadd(i,2) foundadd(i,2) = foundadd(j,2) foundadd(j,2) = add diff --git a/src/ex1.c b/src/ex1.c index cae2605..b37e99c 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -1,12 +1,17 @@ -#include #include #include -#include "stimsyr.h" + #include "read2.h" +#include "stimsyr.h" +#include "get_s2.h" #undef __FUNCT__ #define __FUNCT__ "main" +void solvequad(double *a, double *b, double *c, double *res){ + *res = -*b/(2.0*(*a)) + sqrt((*b)*(*b) - 4.0*(*a)*(*c))/(2.0*(*a)); +} + int main(int argc,char **argv) { Mat A; /* problem matrix */ @@ -14,7 +19,13 @@ int main(int argc,char **argv) EPSType type; PetscReal error,tol,re,im; PetscReal norm=0.0; + PetscReal norm2=0.0; + PetscReal norm3=0.0; + PetscReal norm4=0.0; PetscReal normfin=0.0; + PetscReal normfin2=0.0; + PetscReal normfin3=0.0; + PetscReal normfin4=0.0; PetscScalar kr,ki,value[700]; Vec xr,xi; PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol; @@ -25,12 +36,7 @@ int main(int argc,char **argv) //PetscScalar eigr; //PetscScalar eigi; int mpiid; - int natomax=700, iaa2, iaa; - int ideter[natomax]; - int ideter2[natomax]; - PetscScalar densmat[4][4]={0.0}; - PetscScalar densmatfin[4][4]={0.0}; - PetscScalar trace; + int natomax=700; char const* const fileName = argv[1]; FILE* file = fopen(fileName, "r"); @@ -39,52 +45,83 @@ int main(int argc,char **argv) Data_new(file, &getdata); nlocal = getdata.n/getdata.npar; -//printf("n=%ld\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz); +//printf("n=%ld\t nsites=%d\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.natom, getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz); - PetscScalar valxr[nlocal]; - PetscInt indxr[nlocal]; + PetscScalar *valxr; + PetscInt indxr[nlocal]; //Vec Vr,Vi; char filename[PETSC_MAX_PATH_LEN]="FIL666"; PetscViewer viewer; PetscBool ishermitian; - PetscInt kk,ll,iii2; + PetscInt kk,ll,mm,nn,iii2,iiii; + PetscInt ii; long int iii; - long int ii; long int tcountcol2,tcol[700],tcountcol[getdata.nnz]; double val[700]; - double xmat=0.0; PetscReal xymat=0.0; + PetscReal xymat2=0.0; + PetscReal xymat3=0.0; + PetscReal xymat4=0.0; PetscReal xymatfin=0.0; + PetscReal xymatfin2=0.0; + PetscReal xymatfin3=0.0; + PetscReal xymatfin4=0.0; + PetscReal weight3fin = 0.0; PetscReal XS = 0.0; - int kko,kok,kkio; - Vec xiaa2; /* initial vector, destination vector */ - VecScatter scatter; /* scatter context */ - IS from, to; /* index sets that define the scatter */ - PetscScalar *values; - int idx_from[] = {0}, idx_to[] = {0}; + PetscReal XS2 = 0.0; + PetscReal XS3 = 0.0; + PetscReal XS4 = 0.0; + PetscReal W3 = 0.0; + PetscReal weight3 = 0.0; + PetscReal trace1rdm=0.0; + PetscReal trace1rdmfin=0.0; + PetscReal trace2rdm=0.0; + PetscReal trace2rdmfin=0.0; + IS from, to; /* index sets that define the scatter */ + PetscInt idx_to[nlocal], idx_from[nlocal]; + PetscScalar *values; + int ndim=(getdata.natom/2)*((getdata.natom/2)-1)/2; + double a, b, c; + double gamma_p = 0.0, gamma_m = 0.0; + double gamma_pfin = 0.0, gamma_mfin = 0.0; + double nel, s2dens; + double nelfin, s2densfin; + double densmat2[getdata.natom][getdata.natom][getdata.natom][getdata.natom]; + memset(densmat2, 0, sizeof(densmat2)); + // fn = x^2+x^3 +//PetscInt range[] ={165094,310638,438886,551954,651820,740324,819168,889916,953994,1012690,1067154,1118398,1167296,1214584,1260860,1306584,1352078,1517172,1662716,1790964,1904032,2003898,2092402,2171246,2241994,2306072,2364768,2419232,2470476,2519374,2566662,2612938,2658662,2704156,2869250,3014794,3143042,3256110,3355976,3444480,3523324,3594072,3658150,3716846,3771310,3822554,3871452,3918740,3965016,4010740,4056234,4221328,4366872,4495120,4608188,4708054,4796558,4875402,4946150,5010228,5068924,5123388,5174632,5223530,5270818,5317094,5362818,5408312,5573406,5718950,5847198,5960266,6060132,6148636,6227480,6298228,6362306,6421002,6475466,6526710,6575608,6622896,6669172,6714896,6760390,6925484,7071028,7199276,7312344,7412210,7500714,7579558,7650306,7714384,7773080,7827544,7878788,7927686,7974974,8021250,8066974,8112468,8277562,8423106,8551354,8664422,8764288,8852792,8931636,9002384,9066462,9125158,9179622,9230866,9279764,9327052,9373328,9419052,9464546,9629640,9775184,9903432,10016500,10116366,10204870,10283714,10354462,10418540,10477236,10531700,10582944,10631842,10679130,10725406,10771130,10816624,10981718,11127262,11255510,11368578,11468444,11556948,11635792,11706540,11770618,11829314,11883778,11935022,11983920,12031208,12077484,12123208,12168702,12333796,12479340,12607588,12720656,12820522,12909026,12987870,13058618,13122696,13181392,13235856,13287100,13335998,13383286,13429562,13475286,13520780,13685874,13831418,13959666,14072734,14172600,14261104,14339948,14410696,14474774,14533470,14587934,14639178,14688076,14735364,14781640,14827364,14872858,15037952,15183496,15311744,15424812,15524678,15613182,15692026,15762774,15826852,15885548,15940012,15991256,16040154,16087442,16133718,16179442,16224936}; + SlepcInitialize(&argc,&argv,(char*)0,NULL); ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr); ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); - ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,getdata.nnz*getdata.npar,NULL,getdata.nnz*getdata.npar,NULL,&A);CHKERRQ(ierr); - ierr = MatMPIAIJSetPreallocation(A,getdata.nnz*getdata.npar,NULL,getdata.nnz*getdata.npar,NULL);CHKERRQ(ierr); + ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,2.0*getdata.natom,NULL,2.0*getdata.natom,NULL,&A);CHKERRQ(ierr); + ierr = MatMPIAIJSetPreallocation(A,getdata.natom,NULL,getdata.natom,NULL);CHKERRQ(ierr); //ierr = MatSetFromOptions(A);CHKERRQ(ierr); //ierr = MatSetUp(A);CHKERRQ(ierr); MPI_Comm_rank(MPI_COMM_WORLD,&mpiid); ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); ierr = PetscTime(&tt1);CHKERRQ(ierr); + ierr = PetscPrintf(PETSC_COMM_WORLD," start: %d end: %d\n",Istart, Iend);CHKERRQ(ierr); +// Iend = range[mpiid]; +// if(mpiid==0){ +// Istart = 0; +// } +// else{ +// Istart = range[mpiid-1]; +// } for (i=Istart; i0) { + if (0) { PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer); PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB); PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU); @@ -209,218 +239,125 @@ int main(int argc,char **argv) Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and ki (imaginary part) */ + Vec vec2; + VecScatter scatter; /* scatter context */ ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr); xymat = 0.0; + xymat2 = 0.0; + xymat3 = 0.0; + xymat4 = 0.0; + weight3 = 0.0; norm = 0.0; + norm2 = 0.0; + norm3 = 0.0; + norm4 = 0.0; - for (ii=Istart; ii +#include +#include +#include +#include +#include +#include "get_dmat.h" + + +/* + *--------------------------------------- + * One particle density matrix + *---------------------------------------- + * + * The One particle density matrix + * Input + * ===== + * valxr = The full vector + * Istart = Local starting id + * Iend = Local ending id + * natom = number of sites + * Output + * ===== + * trace = trace + */ +void get_1rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace1rdm){ + + const int natomax=700; + long int ideter[natomax]; + long int ideter2[natomax]; + int kko,kok,kkio; + long int ii; + PetscInt iiii; + long int iii; + long int iaa2, iaa; + int ndim=(*natom)*(*natom)/8-(*natom)/2; + double densmat[ndim][ndim]; + memset(densmat, 0, sizeof(densmat[0][0]) * ndim * ndim); + + for(kko=0;kko<(*natom/2);kko++){ + for(kok=0;kok<(*natom/2);kok++){ + + for(ii=*Istart;ii<*Iend;ii++) { + iii = ii + 1; + iiii = ii; + getdet_(&iii, ideter); + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + if(ideter[kko] != 3){ + if(ideter[kok] == 3){ + ideter2[kok]=ideter[kko]; + ideter2[kko]=3; + adr_(ideter2, &iaa2); + densmat[kko][kok]=densmat[kko][kok]+valxr[iiii]*valxr[iaa2]; + } + } + if(kko == kok && ideter[kko] != 3){ + densmat[kko][kko]=densmat[kko][kko]+valxr[iiii]*valxr[iiii]; + } + + } + if(kko == kok){ + *trace1rdm+=densmat[kko][kko]; + } + + } + } +} /** END **/ + +/* + * + *---------------------------------------- + * two particle density matrix + *---------------------------------------- + * Input + * ===== + * valxr = The full vector + * Istart = Local starting id + * Iend = Local ending id + * Output + * ===== + * trace = trace + */ +void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, double densmat2[*natom][*natom][*natom][*natom]){ + + const int natomax=700; + long int ideter[natomax]; + long int ideter2[natomax]; + int kko,kok,kkio; + int mmo,mom,mmio; + long int ii; + PetscInt iiii; + long int iii; + long int iaa2, iaa; + long int nrow=-1, ncol=-1; +//int ndim=(*natom/2)*((*natom/2)-1)/2; +//double densmat2[ndim][ndim]; +//memset(densmat2, 0, sizeof(densmat2[0][0]) * ndim * ndim); + + for(kko=0;kko<(*natom/2);kko++){ + for(kok=0;kok<(*natom/2);kok++){ + + nrow=nrow+1; + ncol=-1; + for(mmo=0;mmo<(*natom/2);mmo++){ + for(mom=0;mom<(*natom/2);mom++){ + + ncol=ncol+1; + + for(ii=*Istart;ii<*Iend;ii++) { + iii = ii + 1; + iiii = ii; + getdet_(&iii, ideter); + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + if(ideter[kko] == 3 && ideter[kok] == 3 && kko != kok && mmo != mom){ + ideter2[kko]=ideter[mmo]; + ideter2[mmo]=3; + ideter2[kok]=ideter[mom]; + ideter2[mom]=3; + adr_(ideter2, &iaa2); + densmat2[kko][kok][mmo][mom]=densmat2[kko][kok][mmo][mom]+valxr[iiii]*valxr[iaa2]; + } + + + if(kko == mmo && kok == mom && ideter[kko]==3 && ideter[kok]==3 && kko != kok){ + densmat2[kko][kok][mmo][mom]=densmat2[kko][kok][mmo][mom]+valxr[iiii]*valxr[iiii]; + } + + } + printf("%d\t%d\t%d\t%d\t%18f\n",kko,kok,mmo,mom,densmat2[kko][kok][mmo][mom]); + + + if(kko == mmo && kok == mom)*trace2rdm+=densmat2[kko][kok][mmo][mom]; + } + } + + + } + } + +} /** END **/ + diff --git a/src/get_dmat.h b/src/get_dmat.h new file mode 100644 index 0000000..29539e5 --- /dev/null +++ b/src/get_dmat.h @@ -0,0 +1,9 @@ +#include +#include +#include +#include +#include +#include + +void get_1rdm(PetscScalar *, PetscInt *, PetscInt *, int *, PetscReal *); +void get_2rdm(PetscScalar *, PetscInt *, PetscInt *, int *, PetscReal *, double ****); diff --git a/src/get_s2.c b/src/get_s2.c new file mode 100644 index 0000000..4fb42f3 --- /dev/null +++ b/src/get_s2.c @@ -0,0 +1,683 @@ +#include +#include +#include +#include +#include +#include +#include "get_s2.h" +#include "get_val_iaa2.h" + +/* + * This function simply calculates the S^2 value of the wavefunction + * Input + * ===== + * Vr = The full vector + * Istart = Local starting id of the vector + * Iend = Local vector ending id + * valxr = Local vector values + * natom = number of orbitals + * Output + * ====== + * norm = norm of the vector + * xymat = the S^2 value + */ + +void get_s2(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *natom, PetscReal *norm, PetscReal *norm2, PetscReal *norm3, PetscReal *norm4, PetscReal *xymat, PetscReal *xymat2, PetscReal *xymat3, PetscReal *xymat4, PetscReal *weight3, + int *s21a1, int *s21a2, int *s21b1, int *s21b2, int *s22a1, int *s22a2, int *s22b1, int *s22b2, int *s23a1, int *s23a2, int *s23b1, int *s23b2, int *postrou){ + const int natomax=700; + long int iaa2, iaa; + long int iii; + long int ideter[natomax]; + long int ideter2[natomax]; + int kko,kok,kkio; + long int ii; + double xmat=0.0; + double xmat2=0.0; + double xmat3=0.0; + double xmat4=0.0; + double getvaliaa2; + PetscLogDouble t1,t2,tt1,tt2; + PetscErrorCode ierr; + PetscInt iiii; + int ntrouboit1=0; + int ntrouboit2=0; + int ntrouboit3=0; + int okboit1=0; + int okboit2=0; + int okboit3=0; + int mpiid; + int pos1=0; + int pos2=0; + int pos3=0; + MPI_Comm_rank(MPI_COMM_WORLD,&mpiid); +//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);} +//ierr = PetscTime(&tt1);CHKERRQ(ierr); + for(ii=*Istart;ii<*Iend;ii++) { + iii = ii + 1; +// iiii = ii-*Istart; + iiii = ii; + xmat = 0.0; + xmat2 = 0.0; + xmat3 = 0.0; + xmat4 = 0.0; + ntrouboit1 = 0; + ntrouboit2 = 0; + ntrouboit3 = 0; + okboit1 = 0; + okboit2 = 0; + okboit3 = 0; + pos1 = 0; + pos2 = 0; + pos3 = 0; + getdet_(&iii, ideter); + *norm=*norm+valxr[iiii]*valxr[iiii]; + for(kko=*s21a1;kko<=*s21a2;kko++){ + if(ideter[kko]==3){ + ntrouboit1++; + pos1=kko; + } + } + for(kko=*s22a1;kko<=*s22a2;kko++){ + if(ideter[kko]==3){ + ntrouboit2++; + pos2=kko; + } + } + for(kko=*s23a1;kko<=*s23a2;kko++){ + if(ideter[kko]==3){ + ntrouboit3++; + pos3=kko; + } + } + if(ntrouboit1==1 && pos1 == *postrou)okboit1=1; + if(ntrouboit2==1 && pos2 == *postrou)okboit2=1; + if(ntrouboit3==1 && pos3 == *postrou)okboit3=1; + if(okboit1){ + *norm2=*norm2+valxr[iiii]*valxr[iiii]; + } + if(okboit2){ + *norm3=*norm3+valxr[iiii]*valxr[iiii]; + } + if(okboit3){ + *norm4=*norm4+valxr[iiii]*valxr[iiii]; + } + /* + * calculate the weight of ms=5/2 + * + * loop over the determinants to see if we have a S=5/2 + */ + int countw = 0; + for(kko=*s21a1;kko<=*s21a2;kko++){ + if(ideter[kko] == 2) countw=1; + } + for(kok=*s21b1;kok<=*s21b2;kok++){ + if(ideter[kok] == 2) countw=1; + } + if(countw==0 && okboit1){ + *weight3 += (valxr[iiii]*valxr[iiii]); + } + for(kko=0;kko<=(*natom/2)-1;kko++){ + for(kok=kko;kok<=(*natom/2)-1;kok++){ + if(kok == kko && ideter[kok] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko] == 1 && ideter[kok] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 1 && ideter[kok] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=2; + ideter2[kok]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=1; + ideter2[kok]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + for(kko=(*natom/2);kko<=*natom-1;kko++){ + for(kok=kko;kok<=*natom-1;kok++){ + if(kok == kko && ideter[kok] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko] == 1 && ideter[kok] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 1 && ideter[kok] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=2; + ideter2[kok]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=1; + ideter2[kok]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + for(kko=0;kko<=(*natom/2)-1;kko++){ + for(kok=(*natom/2);kok<=*natom-1;kok++){ + if(kok == kko && ideter[kok] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko] == 1 && ideter[kok] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 1 && ideter[kok] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=2; + ideter2[kok]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=1; + ideter2[kok]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; +// if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);} + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + *xymat=*xymat+xmat; + *xymat2=*xymat2+xmat2; + *xymat3=*xymat3+xmat3; + *xymat4=*xymat4+xmat4; +// if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3); + } + + ierr = PetscTime(&tt2);CHKERRQ(ierr); +//printf(" norm = %18f weight = %18f weight/N = %18f tmpwe = %18f\n", *norm2, *weight3, *weight3/(*norm2),tmpwe); +//printf(" norm = %18f %18f xymat = %18f %18f\n", *norm2, *norm3, *xymat2, *xymat3); +//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr); +} diff --git a/src/get_s2.h b/src/get_s2.h new file mode 100644 index 0000000..2eb0415 --- /dev/null +++ b/src/get_s2.h @@ -0,0 +1,50 @@ +#include +#include +#include +#include +#include + +void get_s2(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *); + +void get_s2_mov(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *); + +void get_s2_cyclic(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *, + int *); diff --git a/src/get_s2_cyclic.c b/src/get_s2_cyclic.c new file mode 100644 index 0000000..6204863 --- /dev/null +++ b/src/get_s2_cyclic.c @@ -0,0 +1,784 @@ +#include +#include +#include +#include +#include +#include +#include "get_s2.h" +#include "get_val_iaa2.h" + +/* + * This function simply calculates the S^2 value of the wavefunction + * Input + * ===== + * Vr = The full vector + * Istart = Local starting id of the vector + * Iend = Local vector ending id + * valxr = Local vector values + * natom = number of orbitals + * Output + * ====== + * norm = norm of the vector + * xymat = the S^2 value + */ + +void get_s2_cyclic(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *natom, PetscReal *norm, PetscReal *norm2, PetscReal *norm3, PetscReal *norm4, PetscReal *xymat, PetscReal *xymat2, PetscReal *xymat3, PetscReal *xymat4, + int *s21a1, int *s21a2, int *s21b1, int *s21b2, int *s22a1, int *s22a2, int *s22b1, int *s22b2, int *s23a1, int *s23a2, int *s23b1, int *s23b2, int *postrou){ + const int natomax=700; + long int iaa2, iaa; + long int iii; + long int ideter[natomax]; + long int ideter2[natomax]; + int kko,kok,kkio,kk; + int kko2,kok2; + long int ii; + double xmat=0.0; + double xmat2=0.0; + double xmat3=0.0; + double xmat4=0.0; + double getvaliaa2; + PetscLogDouble t1,t2,tt1,tt2; + PetscErrorCode ierr; + PetscInt iiii; + int ntrouboit1=0; + int ntrouboit2=0; + int ntrouboit3=0; + int okboit1=0; + int okboit2=0; + int okboit3=0; + int mpiid; + int pos1=0; + int pos2=0; + int pos3=0; + MPI_Comm_rank(MPI_COMM_WORLD,&mpiid); +//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);} +//ierr = PetscTime(&tt1);CHKERRQ(ierr); + for(ii=*Istart;ii<*Iend;ii++) { + iii = ii + 1; +// iiii = ii-*Istart; + iiii = ii; + xmat = 0.0; + xmat2 = 0.0; + xmat3 = 0.0; + xmat4 = 0.0; + ntrouboit1 = 0; + ntrouboit2 = 0; + ntrouboit3 = 0; + okboit1 = 0; + okboit2 = 0; + okboit3 = 0; + pos1 = 0; + pos2 = 0; + pos3 = 0; + getdet_(&iii, ideter); + *norm=*norm+valxr[iiii]*valxr[iiii]; + + for(kko=0;kko<*natom/2;kko++){ + if(ideter[kko]==3){ + kk=kko; + } + } + + *postrou = kk; + *s21a1 = kk-1; +// if(*s21a1<0){ +// *s21a1 = (*natom/2) + (*s21a1)%(*natom/2); +// } +// else{ +// *s21a1 = (*s21a1)%(*natom/2); +// } + *s22a1 = kk-1; +// if(*s22a1<0){ +// *s22a1 = (*natom/2) + (*s22a1)%(*natom/2); +// } +// else{ +// *s22a1 = (*s22a1)%(*natom/2); +// } + *s23a1 = kk-2; +// if(*s23a1<0){ +// *s23a1 = (*natom/2) + (*s23a1)%(*natom/2); +// } +// else{ +// *s23a1 = (*s23a1)%(*natom/2); +// } +// + *s21a2 = kk+1; +// if(*s21a2<0){ +// *s21a2 = (*natom/2) + (*s21a2)%(*natom/2); +// } +// else{ +// *s21a2 = (*s21a2)%(*natom/2); +// } + *s22a2 = kk+2; +// if(*s22a2<0){ +// *s22a2 = (*natom/2) + (*s22a2)%(*natom/2); +// } +// else{ +// *s22a2 = (*s22a2)%(*natom/2); +// } + *s23a2 = kk+2; +// if(*s23a2<0){ +// *s23a2 = (*natom/2) + (*s23a2)%(*natom/2); +// } +// else{ +// *s23a2 = (*s23a2)%(*natom/2); +// } +// + *s21b1 = *natom + *s21a1; + *s22b1 = *natom + *s22a1; + *s23b1 = *natom + *s23a1; + + *s21b2 = *natom + *s21a2; + *s22b2 = *natom + *s22a2; + *s23b2 = *natom + *s23a2; + +// if(mpiid==0)printf("postrou = %d\n",*postrou); +// if(mpiid==0)printf("1a1 = %d, 1a2 = %d, 1b1 = %d, 1b2 = %d\n",*s21a1,*s21a2,*s21b1,*s21b2); +// if(mpiid==0)printf("2a1 = %d, 2a2 = %d, 2b1 = %d, 2b2 = %d\n",*s22a1,*s22a2,*s22b1,*s22b2); +// if(mpiid==0)printf("3a1 = %d, 3a2 = %d, 3b1 = %d, 3b2 = %d\n",*s23a1,*s23a2,*s23b1,*s23b2); + +// for(kko=*s21a1;kko<=*s21a2;kko++){ +// if(ideter[kko]==3){ +// ntrouboit1++; +// pos1=kko; +// } +// } +// for(kko=*s22a1;kko<=*s22a2;kko++){ +// if(ideter[kko]==3){ +// ntrouboit2++; +// pos2=kko; +// } +// } +// for(kko=*s23a1;kko<=*s23a2;kko++){ +// if(ideter[kko]==3){ +// ntrouboit3++; +// pos3=kko; +// } +// } +// if(ntrouboit1==1 && pos1 == *postrou)okboit1=1; +// if(ntrouboit2==1 && pos2 == *postrou)okboit2=1; +// if(ntrouboit3==1 && pos3 == *postrou)okboit3=1; + okboit1 = 1; + okboit2 = 1; + okboit3 = 1; + if(okboit1){ + *norm2=*norm2+valxr[iiii]*valxr[iiii]; + } + if(okboit2){ + *norm3=*norm3+valxr[iiii]*valxr[iiii]; + } + if(okboit3){ + *norm4=*norm4+valxr[iiii]*valxr[iiii]; + } + for(kko=kk-2;kko<=kk+(*natom/2-3);kko++){ + for(kok=kko;kok<=kk+(*natom/2-3);kok++){ + kko2=kko; + if(kko2<0){ + kko2 = (*natom/2) + (kko2)%(*natom/2); + } + else{ + kko2 = (kko2)%(*natom/2); + } + kok2=kok; + if(kok2<0){ + kok2 = (*natom/2) + (kok2)%(*natom/2); + } + else{ + kok2 = (kok2)%(*natom/2); + } + + if(kok == kko && ideter[kok2] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko2] == 1 && ideter[kok2] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko2] == 2 && ideter[kok2] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko2] == 1 && ideter[kok2] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko2]=2; + ideter2[kok2]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko2] == 2 && ideter[kok2] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko2]=1; + ideter2[kok2]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + for(kko=*natom+kk-2;kko<=*natom+kk+(*natom/2-3);kko++){ + for(kok=kko;kok<=*natom+kk+(*natom/2-3);kok++){ + kko2=kko-*natom; + if(kko2<0){ + kko2 = (*natom/2) + (kko2)%(*natom/2); + } + else{ + kko2 = (kko2)%(*natom/2); + } + kok2=kok-*natom; + if(kok2<0){ + kok2 = (*natom/2) + (kok2)%(*natom/2); + } + else{ + kok2 = (kok2)%(*natom/2); + } + kko2 = (*natom) - 1 - kko2; + kok2 = (*natom) - 1 - kok2; + + if(kok == kko && ideter[kok2] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko2] == 1 && ideter[kok2] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko2] == 2 && ideter[kok2] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko2] == 1 && ideter[kok2] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko2]=2; + ideter2[kok2]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko2] == 2 && ideter[kok2] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko2]=1; + ideter2[kok2]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + for(kko=kk-2;kko<=kk+(*natom/2-3);kko++){ + for(kok=*natom + kk-2;kok<=*natom + kk+(*natom/2-3);kok++){ + kko2=kko; + if(kko2<0){ + kko2 = (*natom/2) + (kko2)%(*natom/2); + } + else{ + kko2 = (kko2)%(*natom/2); + } + kok2=kok-*natom; + if(kok2<0){ + kok2 = (*natom/2) + (kok2)%(*natom/2); + } + else{ + kok2 = (kok2)%(*natom/2); + } + kok2 = (*natom) - 1 - kok2; + + if(kok == kko && ideter[kok2] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko2] == 1 && ideter[kok2] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko2] == 2 && ideter[kok2] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko2] == 1 && ideter[kok2] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko2]=2; + ideter2[kok2]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko2] == 2 && ideter[kok2] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko2]=1; + ideter2[kok2]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; +// if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);} + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + *xymat=*xymat+xmat; + *xymat2=*xymat2+xmat2; + *xymat3=*xymat3+xmat3; + *xymat4=*xymat4+xmat4; +// if(mpiid==0)printf(" ii = %d xmat3 = %18f xmat4 = %18f diff = %18f\n", ii, xmat3, xmat4, (xmat3-xmat4)); + } + + ierr = PetscTime(&tt2);CHKERRQ(ierr); +//if(mpiid==0)printf(" norm3 = %18f norm4 = %18f xymat3= %18f xymat4= %18f\n", *norm3, *norm4, *xymat3, *xymat4); +//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr); +} diff --git a/src/get_s2_mov.c b/src/get_s2_mov.c new file mode 100644 index 0000000..a7154fa --- /dev/null +++ b/src/get_s2_mov.c @@ -0,0 +1,683 @@ +#include +#include +#include +#include +#include +#include +#include "get_s2.h" +#include "get_val_iaa2.h" + +/* + * This function simply calculates the S^2 value of the wavefunction + * Input + * ===== + * Vr = The full vector + * Istart = Local starting id of the vector + * Iend = Local vector ending id + * valxr = Local vector values + * natom = number of orbitals + * Output + * ====== + * norm = norm of the vector + * xymat = the S^2 value + */ + +void get_s2_mov(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *natom, PetscReal *norm, PetscReal *norm2, PetscReal *norm3, PetscReal *norm4, PetscReal *xymat, PetscReal *xymat2, PetscReal *xymat3, PetscReal *xymat4, PetscReal *weight3, + int *s21a1, int *s21a2, int *s21b1, int *s21b2, int *s22a1, int *s22a2, int *s22b1, int *s22b2, int *s23a1, int *s23a2, int *s23b1, int *s23b2, int *postrou){ + const int natomax=700; + long int iaa2, iaa; + long int iii; + long int ideter[natomax]; + long int ideter2[natomax]; + int kko,kok,kkio; + long int ii; + double xmat=0.0; + double xmat2=0.0; + double xmat3=0.0; + double xmat4=0.0; + double getvaliaa2; + PetscLogDouble t1,t2,tt1,tt2; + PetscErrorCode ierr; + PetscInt iiii; + int ntrouboit1=0; + int ntrouboit2=0; + int ntrouboit3=0; + int okboit1=0; + int okboit2=0; + int okboit3=0; + int mpiid; + int pos1=0; + int pos2=0; + int pos3=0; + MPI_Comm_rank(MPI_COMM_WORLD,&mpiid); +//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);} +//ierr = PetscTime(&tt1);CHKERRQ(ierr); + for(ii=*Istart;ii<*Iend;ii++) { + iii = ii + 1; +// iiii = ii-*Istart; + iiii = ii; + xmat = 0.0; + xmat2 = 0.0; + xmat3 = 0.0; + xmat4 = 0.0; + ntrouboit1 = 0; + ntrouboit2 = 0; + ntrouboit3 = 0; + okboit1 = 0; + okboit2 = 0; + okboit3 = 0; + pos1 = 0; + pos2 = 0; + pos3 = 0; + getdet_(&iii, ideter); + *norm=*norm+valxr[iiii]*valxr[iiii]; + for(kko=*s21a1;kko<=*s21a2;kko++){ + if(ideter[kko]==3){ + ntrouboit1++; + pos1=kko; + } + } + for(kko=*s22a1;kko<=*s22a2;kko++){ + if(ideter[kko]==3){ + ntrouboit2++; + pos2=kko; + } + } + for(kko=*s23a1;kko<=*s23a2;kko++){ + if(ideter[kko]==3){ + ntrouboit3++; + pos3=kko; + } + } + if(ntrouboit1==1 && *s21a1 <= pos1 && pos1 <= *s21a2)okboit1=1; + if(ntrouboit2==1 && pos2 == *postrou)okboit2=1; + if(ntrouboit3==1 && pos3 == *postrou)okboit3=1; + if(okboit1){ + *norm2=*norm2+valxr[iiii]*valxr[iiii]; + } + if(okboit2){ + *norm3=*norm3+valxr[iiii]*valxr[iiii]; + } + if(okboit3){ + *norm4=*norm4+valxr[iiii]*valxr[iiii]; + } + /* + * calculate the weight of ms=5/2 + * + * loop over the determinants to see if we have a S=5/2 + */ + int countw = 0; + for(kko=*s21a1;kko<=*s21a2;kko++){ + if(ideter[kko] == 2) countw=1; + } + for(kok=*s21b1;kok<=*s21b2;kok++){ + if(ideter[kok] == 2) countw=1; + } + if(countw==0 && okboit1){ + *weight3 += (valxr[iiii]*valxr[iiii]); + } + for(kko=0;kko<=(*natom/2)-1;kko++){ + for(kok=kko;kok<=(*natom/2)-1;kok++){ + if(kok == kko && ideter[kok] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko] == 1 && ideter[kok] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 1 && ideter[kok] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=2; + ideter2[kok]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=1; + ideter2[kok]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21a1 && kok <=*s21a2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22a1 && kok <=*s22a2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23a1 && kok <=*s23a2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + for(kko=(*natom/2);kko<=*natom-1;kko++){ + for(kok=kko;kok<=*natom-1;kok++){ + if(kok == kko && ideter[kok] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko] == 1 && ideter[kok] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 1 && ideter[kok] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=2; + ideter2[kok]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=1; + ideter2[kok]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21b1 && kko <=*s21b2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22b1 && kko <=*s22b2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23b1 && kko <=*s23b2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + for(kko=0;kko<=(*natom/2)-1;kko++){ + for(kok=(*natom/2);kok<=*natom-1;kok++){ + if(kok == kko && ideter[kok] != 3){ + xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + else{ + if(ideter[kko] == 1 && ideter[kok] == 1){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 2){ + xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + } + if(ideter[kko] == 1 && ideter[kok] == 2){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=2; + ideter2[kok]=1; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + if(ideter[kko] == 2 && ideter[kok] == 1){ + xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]); + } + } + } + for(kkio=0;kkio<=*natom-1;kkio++){ + ideter2[kkio]=ideter[kkio]; + } + ideter2[kko]=1; + ideter2[kok]=2; + adr_(ideter2, &iaa2); + iaa2 = iaa2 - 1; +// if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);} + xmat=xmat+valxr[iiii]*valxr[iaa2]; + if(okboit1){ + if( kko >=*s21a1 && kko <=*s21a2){ + if( kok >=*s21b1 && kok <=*s21b2){ + xmat2=xmat2+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit2){ + if( kko >=*s22a1 && kko <=*s22a2){ + if( kok >=*s22b1 && kok <=*s22b2){ + xmat3=xmat3+(valxr[iiii]*valxr[iaa2]); + } + } + } + if(okboit3){ + if( kko >=*s23a1 && kko <=*s23a2){ + if( kok >=*s23b1 && kok <=*s23b2){ + xmat4=xmat4+(valxr[iiii]*valxr[iaa2]); + } + } + } + } + } + } + } + *xymat=*xymat+xmat; + *xymat2=*xymat2+xmat2; + *xymat3=*xymat3+xmat3; + *xymat4=*xymat4+xmat4; +// if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3); + } + + ierr = PetscTime(&tt2);CHKERRQ(ierr); +//printf(" norm = %18f weight = %18f weight/N = %18f tmpwe = %18f\n", *norm2, *weight3, *weight3/(*norm2),tmpwe); +//printf(" norm = %18f %18f xymat = %18f %18f\n", *norm2, *norm3, *xymat2, *xymat3); +//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr); +} diff --git a/src/getdet.irp.f b/src/getdet.irp.f index cd483af..a4e451c 100644 --- a/src/getdet.irp.f +++ b/src/getdet.irp.f @@ -7,7 +7,7 @@ subroutine getdet(add,ideter) integer,INTENT(INOUT)::ideter(natomax) integer(kind=selected_int_kind(16)),INTENT(IN)::add integer(kind=selected_int_kind(16))::deta,detb - integer::i,const,ia,ib + integer::i,const,ia,ib, natom2 ib = MOD(add,nt2) if(MOD(add,nt2).eq.0)then @@ -20,33 +20,38 @@ subroutine getdet(add,ideter) detb=0 deta=0 i=1 - do while (i.le.(ib)) - const=1 - do while(popcnt(detb).ne.nbeta .or. const==1) - if(nbeta.eq.0)then - detb=0 - EXIT - endif - detb+=1 - const=0 - enddo - i+=1 -! write(6,14)detb,detb - enddo - i=1 - do while (i.le.(ia)) - const=1 - do while(popcnt(deta).ne.ntrou .or. const==1) - deta+=1 - const=0 - enddo - i+=1 -! write(6,14)deta,deta - enddo + detb = det(ib,1) + deta = deth(ia,1) + if(FAM1) deta = ISHFT(deta,-(natom/2)) +! do while (i.le.(ib)) +! const=1 +! do while(popcnt(detb).ne.nbeta .or. const==1) +! detb+=1 +! const=0 +! enddo +! i+=1 +! write(6,14)detb,detb +! enddo +! i=1 +! do while (i.le.(ia)) +! const=1 +! do while(popcnt(deta).ne.ntrou .or. const==1) +! deta+=1 +! const=0 +! enddo +! i+=1 +! write(6,14)deta,deta +! enddo const=0 - do i=0,(natom/2) - 1 + if(FAM1) then + natom2 = natom/2 + else + natom2 = natom + endif + + do i=0,(natom2) - 1 if(BTEST(deta,i))then - ideter((natom/2)-i)=3 + ideter((natom2)-i)=3 endif enddo do i=0,natom-1 diff --git a/src/natom.irp.f b/src/natom.irp.f index 7b8091a..cd47a04 100644 --- a/src/natom.irp.f +++ b/src/natom.irp.f @@ -2,7 +2,6 @@ BEGIN_PROVIDER [integer, natom] &BEGIN_PROVIDER [integer, natrest] &BEGIN_PROVIDER [integer, ial0] &BEGIN_PROVIDER [logical*1, yham] -&BEGIN_PROVIDER [logical*1, FAM1] &BEGIN_PROVIDER [integer, nlientot] &BEGIN_PROVIDER [real*8, xt,(maxlien)] &BEGIN_PROVIDER [real*8 , xjz,(maxlien)] @@ -93,10 +92,11 @@ BEGIN_PROVIDER [integer, natom] enddo !------------------Lecture Hamiltonien - FAM1=.TRUE. +! FAM1=.TRUE. yham=.TRUE. write(6,*)'HAMILTONIEN t-J' write(6,*)'Le nombre de trou est : ',ntrou + write(6,*)'Famille 1 : ',FAM1 !--------------------------------------------- write(6,*)' ' write(6,*)' ' diff --git a/src/providdet.irp.f b/src/providdet.irp.f new file mode 100644 index 0000000..959940f --- /dev/null +++ b/src/providdet.irp.f @@ -0,0 +1,85 @@ +BEGIN_PROVIDER[integer(kind=selected_int_kind(16)),det,(nt2,2)] +&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)),deth,(nt1,2)] + BEGIN_DOC + ! provides det and deth array + END_DOC + implicit none +! integer(kind=selected_int_kind(16))::dethsh + integer(kind=selected_int_kind(16))::a + integer(kind=selected_int_kind(16))::i,count + integer::const + i=1 + a=0 + const=0 + count=0 + + If(ntrou.ge.1)then + + const=0 +! dethsh = ISHFT(deth,-natom/2) +! i=nt1 + do while (i.le.(nt1)) +! if(a.eq.dethsh)then +! addh=i-1 +! EXIT +! endif + + i+=1 + a+=1 + do while(popcnt(a).ne.ntrou) + a+=1 + enddo + count+=1 + if(FAM1) then + deth(count,1)=ISHFT(a,natom/2) + else + deth(count, 1) = a + endif + deth(count,2)=i-1 +! write(6,16)ISHFT(a,natom/2),ISHFT(a,natom/2),i-1 + enddo +! if(a.eq.dethsh )then +! count+=1 +! deth(1,1)=ISHFT(a,natom/2) +! deth(1,2)=nt1 +! endif + + endif + + !C if det=0 then exit + a=0 + i=0 + count=0 + print *,'nt2=',nt2,'nbeta=',nbeta + do while (i.lt.(nt2)) + + i+=1 + a+=1 + do while(popcnt(a).ne.nbeta) + if(nbeta.eq.0)then + a=0 + count+=1 + det(count,1)=a + det(count,2)=i + EXIT + endif + a+=1 + enddo + + if(nbeta.ne.0)then + count+=1 + det(count,1)=a + det(count,2)=i + endif +! write(6,16)a,a,i + enddo + + +10 FORMAT(B64,I8,F8.2) +15 FORMAT(B64,I8,I8,I8) +11 FORMAT(B64,I3,B64) +12 FORMAT(I5,$) +13 FORMAT(B64,B64) +14 FORMAT(B64,I8) +16 FORMAT(B64,I8,I8) +END_PROVIDER diff --git a/src/provide_clust.irp.f b/src/provide_clust.irp.f index fe5e8ef..a529887 100644 --- a/src/provide_clust.irp.f +++ b/src/provide_clust.irp.f @@ -6,6 +6,7 @@ BEGIN_PROVIDER[integer,l1, (maxlien)] &BEGIN_PROVIDER[real*8, xjjxy,(maxlien)] &BEGIN_PROVIDER [integer, ntrou] &BEGIN_PROVIDER [integer, isz] +&BEGIN_PROVIDER [logical*1, FAM1] implicit none ! integer::i ! open(unit=11,file="l1.dat",form="formatted") diff --git a/src/read2.c b/src/read2.c index ec92d41..f2308f0 100644 --- a/src/read2.c +++ b/src/read2.c @@ -1,7 +1,3 @@ -#include -#include -#include -#include #include "read2.h" void Data_new(FILE* file, Data* dat) { @@ -17,25 +13,31 @@ void Data_new(FILE* file, Data* dat) { /* note that fgets don't strip the terminating \n, checking its presence would allow to handle lines longer that sizeof(line) */ - if (count != 12){ + if (count != 26){ count++; switch(count){ case 1: dat->n=atol(line); break; case 2: - dat->nnz=atol(line); + dat->natom=atol(line); break; case 3: - dat->npar=atol(line); + dat->nnz=atol(line); break; case 4: - dat->ntrou=atol(line); + dat->npar=atol(line); break; case 5: - dat->isz=atol(line); + dat->ntrou=atol(line); break; case 6: + dat->isz=atol(line); + break; + case 7: + dat->FAM1 = to_bool(line); + break; + case 8: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -59,7 +61,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 7: + case 9: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -83,7 +85,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 8: + case 10: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -107,7 +109,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 9: + case 11: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -131,7 +133,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 10: + case 12: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -155,7 +157,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 11: + case 13: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -179,9 +181,48 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 12: + case 14: dat->nroots=atol(line); break; + case 15: + dat->s21a1=atol(line); + break; + case 16: + dat->s21a2=atol(line); + break; + case 17: + dat->s21b1=atol(line); + break; + case 18: + dat->s21b2=atol(line); + break; + case 19: + dat->s22a1=atol(line); + break; + case 20: + dat->s22a2=atol(line); + break; + case 21: + dat->s22b1=atol(line); + break; + case 22: + dat->s22b2=atol(line); + break; + case 23: + dat->s23a1=atol(line); + break; + case 24: + dat->s23a2=atol(line); + break; + case 25: + dat->s23b1=atol(line); + break; + case 26: + dat->s23b2=atol(line); + break; + case 27: + dat->postrou=atol(line); + break; } /* end of switch */ } /* end of the input file */ @@ -191,6 +232,14 @@ void Data_new(FILE* file, Data* dat) { //return dat; } +PetscBool to_bool(const char* str) { + PetscBool strflg; + PetscStrcmp("true\n",str, &strflg); + if(!strflg) PetscStrcmp("True\n",str, &strflg); + if(!strflg) PetscStrcmp("TRUE\n",str, &strflg); + return strflg; +} + /* int main(int argc, char* argv[]) { diff --git a/src/read2.h b/src/read2.h index db7021a..dedce9c 100644 --- a/src/read2.h +++ b/src/read2.h @@ -1,13 +1,18 @@ #include -#include #include #include #include +#include +#include + +PetscBool to_bool(const char* str); + typedef struct { PetscInt n; long int nnz,npar; long int ntrou,isz; + PetscBool FAM1; long int l1[700]; long int l2[700]; long int ktyp[700]; @@ -15,6 +20,20 @@ typedef struct { double xjjxy[700]; double xtt[700]; long int nroots; + int natom; + int s21a1; + int s21a2; + int s21b1; + int s21b2; + int s22a1; + int s22a2; + int s22b1; + int s22b2; + int s23a1; + int s23a2; + int s23b1; + int s23b2; + int postrou; } Data ; diff --git a/src/searchdet.irp.f b/src/searchdet.irp.f index d900984..f8ea9b2 100644 --- a/src/searchdet.irp.f +++ b/src/searchdet.irp.f @@ -1,80 +1,70 @@ -subroutine searchdet(det,add,deth,addh) +subroutine searchdet(deti,add,dethi,addh) BEGIN_DOC ! this subroutine is at the heart of the idea ! it will generate all the determinants in a fixed order ! then find the posistion of the determinant given and ! return it's position in add. END_DOC - integer(kind=selected_int_kind(16)),INTENT(INOUT)::det + integer(kind=selected_int_kind(16)),INTENT(INOUT)::deti integer(kind=selected_int_kind(16)),INTENT(INOUT)::add - integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth + integer(kind=selected_int_kind(16)),INTENT(INOUT)::dethi integer(kind=selected_int_kind(16)),INTENT(INOUT)::addh integer(kind=selected_int_kind(16))::dethsh integer(kind=selected_int_kind(16))::a - integer(kind=selected_int_kind(16))::i - integer::const + integer(kind=selected_int_kind(16))::i,j + integer::count + logical::found + i=1 - a=0 - add=0 - const=0 - - If(ntrou.ge.1)then - - const=0 - dethsh = ISHFT(deth,-natom/2) - addh=0 -! i=nt1 - do while (i.le.(nt1)) - if(a.eq.dethsh)then - addh=i-1 - EXIT - endif - - i+=1 - a+=1 - do while(popcnt(a).ne.ntrou) - a+=1 - enddo - enddo - if(a.eq.dethsh .and. addh.eq.0)then - addh=nt1 + j=nt1 + found=.FALSE. + do while(.not.found) + if(deth((i+j)/2,1).eq.dethi)then + addh=deth((i+j)/2,2) + found=.TRUE. + EXIT + elseif (abs(i-j).eq.1)then + if(deth(i,1).eq.dethi)then + addh=deth(i,2) + elseif(deth(j,1).eq.dethi)then + addh=deth(j,2) endif - endif - - !C if det=0 then exit - a=0 - i=0 - count=0 - if(a.eq.det)then - add=1 - Return - endif - - do while (i.le.(nt2)) - if(a.eq.det)then - if(a.eq.1)then - add=i - count=-1 - EXIT + found=.TRUE. + EXIT + else + if(deth((i+j)/2,1).gt.dethi)then + j=(i+j)/2 else - add=i - count=-1 - EXIT + i=(i+j)/2 + endif + endif + enddo + + i=1 + j=nt2 + found=.FALSE. + do while(.not.found) + if(det((i+j)/2,1).eq.deti)then + add=det((i+j)/2,2) + found=.TRUE. + EXIT + elseif (abs(i-j).eq.1)then + if(det(i,1).eq.deti)then + add=det(i,2) + elseif(det(j,1).eq.deti)then + add=det(j,2) + endif + found=.TRUE. + EXIT + else + if(det((i+j)/2,1).gt.deti)then + j=(i+j)/2 + else + i=(i+j)/2 endif endif - - i+=1 - a+=1 -!C write(6,16)a,a,i-2 - do while(popcnt(a).ne.nbeta) - a+=1 - enddo enddo - if(a.eq.det .and. count.ne.-1)then - add=i-1 - endif - 10 FORMAT(B64,I8,F8.2) 15 FORMAT(B64,I8,I8,I8) diff --git a/src/sort.irp.f b/src/sort.irp.f index 1c8c130..14c7407 100644 --- a/src/sort.irp.f +++ b/src/sort.irp.f @@ -1,14 +1,14 @@ subroutine sort() implicit none integer::i,j,ord,ordh - integer(kind=selected_int_kind(16))::add,addh,det,deth,addt + integer(kind=selected_int_kind(16))::add,addh,deti,dethi,addt do i=1,detfound-1 do j=i+1,detfound if(foundaddh(i,1).gt.foundaddh(j,1))then - deth = foundaddh(i,1) + dethi = foundaddh(i,1) foundaddh(i,1) = foundaddh(j,1) - foundaddh(j,1) = deth + foundaddh(j,1) = dethi addh = foundaddh(i,2) foundaddh(i,2) = foundaddh(j,2) foundaddh(j,2) = addh @@ -17,9 +17,9 @@ subroutine sort() foundaddh(j,3) = ordh endif if(foundadd(i,1).gt.foundadd(j,1))then - det = foundadd(i,1) + deti = foundadd(i,1) foundadd(i,1) = foundadd(j,1) - foundadd(j,1) = det + foundadd(j,1) = deti add = foundadd(i,2) foundadd(i,2) = foundadd(j,2) foundadd(j,2) = add diff --git a/src/unit_l1.irp.f b/src/unit_l1.irp.f index 5239a49..2bb5ecf 100644 --- a/src/unit_l1.irp.f +++ b/src/unit_l1.irp.f @@ -9,9 +9,11 @@ tcountcol, & tntrou, & tisz, & + tfam1, & tcol,tval) implicit none integer,INTENT(INOUT)::tistart, tnrows, tntrou, tisz + logical*1,INTENT(INOUT)::tfam1 integer::i real*8,INTENT(INOUT)::tval(maxlien) integer(kind=selected_int_kind(16)),INTENT(INOUT)::tcol(maxlien) @@ -31,6 +33,7 @@ enddo ntrou = tntrou isz = tisz + FAM1 = tfam1 tcol=0 tval=0d0 provide l1 l2 ktyp xtt xjjxy xjjz ntrou @@ -38,6 +41,7 @@ !print *,l1 !print *,"xjjz" !print *,xjjz +!print *,FAM1 call unit(tistart, tcountcol,tcol,tval) end From 62dba9a3c59de70c45b07f3a25b908144ce999d5 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Sun, 28 Jan 2018 11:05:08 +0100 Subject: [PATCH 29/38] now limit the hole movement by giving range --- src/ex1.c | 4 +++- src/getdet.irp.f | 15 +++++++++++++-- src/natom.irp.f | 3 +++ src/nt1.irp.f | 8 +++++++- src/providdet.irp.f | 6 +++++- src/provide_clust.irp.f | 10 ++++++---- src/read2.c | 10 ++++++++-- src/read2.h | 8 +++++--- src/unit_l1.irp.f | 8 +++++++- 9 files changed, 57 insertions(+), 15 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index b37e99c..d27c41f 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -127,13 +127,15 @@ int main(int argc,char **argv) getdata.l2, getdata.ktyp, &iii, - &getdata.nnz, + &getdata.nnz, getdata.xjjxy, getdata.xjjz , getdata.xtt , tcountcol, &getdata.ntrou, &getdata.isz, + &getdata.fix_trou1, + &getdata.fix_trou2, &getdata.FAM1, tcol, val); diff --git a/src/getdet.irp.f b/src/getdet.irp.f index a4e451c..fa98b05 100644 --- a/src/getdet.irp.f +++ b/src/getdet.irp.f @@ -22,7 +22,14 @@ subroutine getdet(add,ideter) i=1 detb = det(ib,1) deta = deth(ia,1) - if(FAM1) deta = ISHFT(deta,-(natom/2)) + if(FAM1) then + if(fix_trou1 .eq. fix_trou2) then + deta = ISHFT(deta,-(natom/2)) + else + natom2 = natom - (fix_trou2 - fix_trou1) + deta = ISHFT(deta, -natom2) + endif + endif ! do while (i.le.(ib)) ! const=1 ! do while(popcnt(detb).ne.nbeta .or. const==1) @@ -44,7 +51,11 @@ subroutine getdet(add,ideter) ! enddo const=0 if(FAM1) then - natom2 = natom/2 + if(fix_trou1 .eq. fix_trou2) then + natom2 = natom/2 + else + natom2 = (fix_trou2 - fix_trou1) + endif else natom2 = natom endif diff --git a/src/natom.irp.f b/src/natom.irp.f index cd47a04..58bf864 100644 --- a/src/natom.irp.f +++ b/src/natom.irp.f @@ -97,6 +97,9 @@ BEGIN_PROVIDER [integer, natom] write(6,*)'HAMILTONIEN t-J' write(6,*)'Le nombre de trou est : ',ntrou write(6,*)'Famille 1 : ',FAM1 + if(FAM1) then + if(fix_trou1 .ne. fix_trou2) write(6,*)'Trou fixe entre :', fix_trou1, "et ", fix_trou2 + endif !--------------------------------------------- write(6,*)' ' write(6,*)' ' diff --git a/src/nt1.irp.f b/src/nt1.irp.f index 64c3e78..11dda96 100644 --- a/src/nt1.irp.f +++ b/src/nt1.irp.f @@ -7,7 +7,13 @@ BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt1] ! call combin(idet1(1,nt1+1),natom,ntrou,nt1,32,jrangmax) natom2=natom - if(FAM1)natom2=natom/2 + if(FAM1) then + if(fix_trou1 .eq. fix_trou2) then + natom2=natom/2 + else + natom2 = fix_trou2 - fix_trou1 + endif + endif 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 diff --git a/src/providdet.irp.f b/src/providdet.irp.f index 959940f..577890a 100644 --- a/src/providdet.irp.f +++ b/src/providdet.irp.f @@ -31,7 +31,11 @@ BEGIN_PROVIDER[integer(kind=selected_int_kind(16)),det,(nt2,2)] enddo count+=1 if(FAM1) then - deth(count,1)=ISHFT(a,natom/2) + if(fix_trou1 .eq. fix_trou2) then + deth(count,1)=ISHFT(a,natom/2) + else + deth(count,1)=ISHFT(a,natom - (fix_trou2-fix_trou1)) + endif else deth(count, 1) = a endif diff --git a/src/provide_clust.irp.f b/src/provide_clust.irp.f index a529887..9ed9988 100644 --- a/src/provide_clust.irp.f +++ b/src/provide_clust.irp.f @@ -1,12 +1,14 @@ 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, xtt ,(maxlien)] &BEGIN_PROVIDER[real*8, xjjz ,(maxlien)] &BEGIN_PROVIDER[real*8, xjjxy,(maxlien)] -&BEGIN_PROVIDER [integer, ntrou] -&BEGIN_PROVIDER [integer, isz] -&BEGIN_PROVIDER [logical*1, FAM1] +&BEGIN_PROVIDER[integer, ntrou] +&BEGIN_PROVIDER[integer, isz] +&BEGIN_PROVIDER[logical*1, FAM1] +&BEGIN_PROVIDER[integer, fix_trou1] +&BEGIN_PROVIDER[integer, fix_trou2] implicit none ! integer::i ! open(unit=11,file="l1.dat",form="formatted") diff --git a/src/read2.c b/src/read2.c index f2308f0..9fcd19e 100644 --- a/src/read2.c +++ b/src/read2.c @@ -13,7 +13,7 @@ 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 != 26){ + if (count != 29){ count++; switch(count){ case 1: @@ -223,6 +223,12 @@ void Data_new(FILE* file, Data* dat) { case 27: dat->postrou=atol(line); break; + case 28: + dat->fix_trou1=atol(line); + break; + case 29: + dat->fix_trou2=atol(line); + break; } /* end of switch */ } /* end of the input file */ @@ -232,7 +238,7 @@ void Data_new(FILE* file, Data* dat) { //return dat; } -PetscBool to_bool(const char* str) { +_Bool to_bool(const char* str) { PetscBool strflg; PetscStrcmp("true\n",str, &strflg); if(!strflg) PetscStrcmp("True\n",str, &strflg); diff --git a/src/read2.h b/src/read2.h index dedce9c..ce325db 100644 --- a/src/read2.h +++ b/src/read2.h @@ -6,13 +6,13 @@ #include #include -PetscBool to_bool(const char* str); +_Bool to_bool(const char* str); typedef struct { PetscInt n; long int nnz,npar; long int ntrou,isz; - PetscBool FAM1; + _Bool FAM1; long int l1[700]; long int l2[700]; long int ktyp[700]; @@ -33,7 +33,9 @@ typedef struct { int s23a2; int s23b1; int s23b2; - int postrou; + long int postrou; + long int fix_trou1; + long int fix_trou2; } Data ; diff --git a/src/unit_l1.irp.f b/src/unit_l1.irp.f index 2bb5ecf..cb61e46 100644 --- a/src/unit_l1.irp.f +++ b/src/unit_l1.irp.f @@ -9,10 +9,14 @@ tcountcol, & tntrou, & tisz, & + tfix_trou1, & + tfix_trou2, & tfam1, & tcol,tval) implicit none - integer,INTENT(INOUT)::tistart, tnrows, tntrou, tisz + integer,INTENT(INOUT)::tistart, tnrows + integer,INTENT(INOUT)::tntrou, tisz + integer,INTENT(INOUT)::tfix_trou1, tfix_trou2 logical*1,INTENT(INOUT)::tfam1 integer::i real*8,INTENT(INOUT)::tval(maxlien) @@ -34,6 +38,8 @@ ntrou = tntrou isz = tisz FAM1 = tfam1 + fix_trou1 = tfix_trou1 + fix_trou2 = tfix_trou2 tcol=0 tval=0d0 provide l1 l2 ktyp xtt xjjxy xjjz ntrou From 55e106f53db26d12a3c0f91c9ea579a5dd71321c Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Sun, 28 Jan 2018 11:16:22 +0100 Subject: [PATCH 30/38] now compiles but bug in the calculation of s2 --- .gitignore | 4 ++++ Makefile | 22 ++++++++++++++++--- src/.gitignore | 5 +++++ src/{adrfull.irp.f => adrfull_old} | 0 src/get_dmat.h | 2 +- src/get_s2.c | 2 +- src/get_s2_cyclic.c | 2 +- src/get_s2_mov.c | 2 +- src/get_val_iaa2.c | 34 ++++++++++++++++++++++++++++++ src/get_val_iaa2.h | 8 +++++++ src/stimsyr.h | 5 +++++ 11 files changed, 79 insertions(+), 7 deletions(-) create mode 100644 .gitignore create mode 100644 src/.gitignore rename src/{adrfull.irp.f => adrfull_old} (100%) create mode 100644 src/get_val_iaa2.c create mode 100644 src/get_val_iaa2.h diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5b2aa47 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*bin/* +*/obj/* +*mod.mod* +irpf90.a diff --git a/Makefile b/Makefile index 5982ec1..6624503 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,7 @@ include ${SLEPC_DIR}/lib/slepc/conf/slepc_common #CC=gcc #FC = ifort +CLINKER = mpicc -fPIC #-wd1572 -O3 -axAVX,SSE4.2 -fno-alias -no-prec-div -no-prec-sqrt -ip MAKE = /usr/bin/make MKDIR_P = /bin/mkdir -p OBJ_DIR := obj @@ -24,14 +25,29 @@ ${BIN_DIR}: directories: ${OBJ_DIR} ${LIB_DIR} ${BIN_DIR} ${LIB_DIR}/irpf90.a: directories - cd ${SRC_DIR} && irpf90 init && $(MAKE) irpf90.a && cp IRPF90_temp/irpf90.a ../${LIB_DIR} + cd ${SRC_DIR} && irpf90 init && $(MAKE) irpf90.a && cp irpf90.a ../${LIB_DIR} ${OBJ_DIR}/read2.o: ${SRC_DIR}/read2.c directories chkopts ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} +${OBJ_DIR}/get_s2.o: ${SRC_DIR}/get_s2.c directories chkopts + ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + +${OBJ_DIR}/get_s2_cyclic.o: ${SRC_DIR}/get_s2_cyclic.c directories chkopts + ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + +${OBJ_DIR}/get_s2_mov.o: ${SRC_DIR}/get_s2_mov.c directories chkopts + ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + +${OBJ_DIR}/get_dmat.o: ${SRC_DIR}/get_dmat.c directories chkopts + ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + +${OBJ_DIR}/get_val_iaa2.o: ${SRC_DIR}/get_val_iaa2.c directories chkopts + ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + ${OBJ_DIR}/ex1.o: ${SRC_DIR}/ex1.c -${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} -${BIN_DIR}/ex1: ${OBJ_DIR}/read2.o ${LIB_DIR}/irpf90.a ${OBJ_DIR}/ex1.o ${SRC_DIR}/read2.h ${SRC_DIR}/stimsyr.h chkopts - -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3 +${BIN_DIR}/ex1: ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${OBJ_DIR}/ex1.o ${SRC_DIR}/read2.h ${SRC_DIR}/stimsyr.h chkopts + -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3 # ${RM} ex1.o read2.o diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..7ac9fbf --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,5 @@ +IRPF90_temp/ +IRPF90_man/ +irpf90.make +irpf90_entities +tags \ No newline at end of file diff --git a/src/adrfull.irp.f b/src/adrfull_old similarity index 100% rename from src/adrfull.irp.f rename to src/adrfull_old diff --git a/src/get_dmat.h b/src/get_dmat.h index 29539e5..a60efa7 100644 --- a/src/get_dmat.h +++ b/src/get_dmat.h @@ -6,4 +6,4 @@ #include void get_1rdm(PetscScalar *, PetscInt *, PetscInt *, int *, PetscReal *); -void get_2rdm(PetscScalar *, PetscInt *, PetscInt *, int *, PetscReal *, double ****); +void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, double densmat2[*natom][*natom][*natom][*natom]); diff --git a/src/get_s2.c b/src/get_s2.c index 4fb42f3..087a578 100644 --- a/src/get_s2.c +++ b/src/get_s2.c @@ -676,7 +676,7 @@ void get_s2(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *n // if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3); } - ierr = PetscTime(&tt2);CHKERRQ(ierr); + ierr = PetscTime(&tt2); //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_cyclic.c b/src/get_s2_cyclic.c index 6204863..e2075c0 100644 --- a/src/get_s2_cyclic.c +++ b/src/get_s2_cyclic.c @@ -778,7 +778,7 @@ void get_s2_cyclic(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, // if(mpiid==0)printf(" ii = %d xmat3 = %18f xmat4 = %18f diff = %18f\n", ii, xmat3, xmat4, (xmat3-xmat4)); } - ierr = PetscTime(&tt2);CHKERRQ(ierr); + ierr = PetscTime(&tt2); //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 index a7154fa..a7f30ba 100644 --- a/src/get_s2_mov.c +++ b/src/get_s2_mov.c @@ -676,7 +676,7 @@ void get_s2_mov(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, in // if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3); } - ierr = PetscTime(&tt2);CHKERRQ(ierr); + ierr = PetscTime(&tt2); //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_val_iaa2.c b/src/get_val_iaa2.c new file mode 100644 index 0000000..e88a35f --- /dev/null +++ b/src/get_val_iaa2.c @@ -0,0 +1,34 @@ +#include +#include +#include +#include +#include +#include +#include "get_val_iaa2.h" + +/* + * This function gets vector from a different processor + * xr the full vector + * iaa2 adresse to get value from + * getvaliaa2 the value + */ + +void get_val_iaa2(Vec xr,long int *iaa2,double *getvaliaa2){ + Vec x; /* initial vector, destination vector */ + VecScatter scatter; /* scatter context */ + IS from, to; /* index sets that define the scatter */ + PetscScalar *values; + int idx_from[] = {*iaa2}, idx_to[] = {0}; + VecCreateSeq( PETSC_COMM_SELF,1,&x); + ISCreateGeneral(PETSC_COMM_SELF,1,idx_from,PETSC_COPY_VALUES,&from); + ISCreateGeneral(PETSC_COMM_SELF,1,idx_to, PETSC_COPY_VALUES,&to); + printf("in get_val"); + VecScatterCreate(xr,from,x,to,&scatter); + VecScatterBegin(scatter,xr,x,INSERT_VALUES,SCATTER_FORWARD); + VecScatterEnd(scatter,xr,x, INSERT_VALUES,SCATTER_FORWARD); + VecGetArray(x,&values); + *getvaliaa2 = values[0]; + ISDestroy(&from); + ISDestroy(&to); + VecScatterDestroy(&scatter); +} diff --git a/src/get_val_iaa2.h b/src/get_val_iaa2.h new file mode 100644 index 0000000..f47bcaa --- /dev/null +++ b/src/get_val_iaa2.h @@ -0,0 +1,8 @@ +#include +#include +#include +#include +#include +#include + +void get_val_iaa2(Vec, long int *, double *); diff --git a/src/stimsyr.h b/src/stimsyr.h index 84dc62f..2397381 100644 --- a/src/stimsyr.h +++ b/src/stimsyr.h @@ -1,3 +1,5 @@ +#include + void unit_l1_( long int *, long int *, @@ -11,5 +13,8 @@ void unit_l1_( long int *, long int *, long int *, + long int *, + _Bool *, + long int *, double *); From 9ebcefaa848db4246d8425b83c7aff480eac2642 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Sun, 28 Jan 2018 12:20:24 +0100 Subject: [PATCH 31/38] now s2 is correct and everything looks good --- src/ex1.c | 2 +- src/get_s2.c | 13 ++++++++----- src/get_s2.h | 16 ++-------------- src/get_s2_cyclic.c | 4 ++-- src/get_s2_mov.c | 4 ++-- src/read2.h | 2 +- 6 files changed, 16 insertions(+), 25 deletions(-) diff --git a/src/ex1.c b/src/ex1.c index d27c41f..88e86c6 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -104,7 +104,7 @@ int main(int argc,char **argv) 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); + ierr = PetscPrintf(PETSC_COMM_WORLD," start: %d end: %d | %d\n",Istart, Iend,getdata.natom);CHKERRQ(ierr); // Iend = range[mpiid]; // if(mpiid==0){ diff --git a/src/get_s2.c b/src/get_s2.c index 087a578..d882d00 100644 --- a/src/get_s2.c +++ b/src/get_s2.c @@ -22,13 +22,16 @@ * 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){ +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 ideter[natomax]; + int ideter2[natomax]; int kko,kok,kkio; long int ii; double xmat=0.0; @@ -678,6 +681,6 @@ void get_s2(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *n ierr = PetscTime(&tt2); //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); +//printf(" norm = %18f %18f xymat = %18f %18f | %d %d %d %d %d\n", *norm, *norm3, *xymat, *xymat3, *s22a1, *s22a2, *s22b1, *s22b2, *postrou); //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 index 2eb0415..5bf7e13 100644 --- a/src/get_s2.h +++ b/src/get_s2.h @@ -4,20 +4,8 @@ #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(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); void get_s2_mov(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, int *, diff --git a/src/get_s2_cyclic.c b/src/get_s2_cyclic.c index e2075c0..ae7a20f 100644 --- a/src/get_s2_cyclic.c +++ b/src/get_s2_cyclic.c @@ -27,8 +27,8 @@ void get_s2_cyclic(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, const int natomax=700; long int iaa2, iaa; long int iii; - long int ideter[natomax]; - long int ideter2[natomax]; + int ideter[natomax]; + int ideter2[natomax]; int kko,kok,kkio,kk; int kko2,kok2; long int ii; diff --git a/src/get_s2_mov.c b/src/get_s2_mov.c index a7f30ba..a829e77 100644 --- a/src/get_s2_mov.c +++ b/src/get_s2_mov.c @@ -27,8 +27,8 @@ void get_s2_mov(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, in const int natomax=700; long int iaa2, iaa; long int iii; - long int ideter[natomax]; - long int ideter2[natomax]; + int ideter[natomax]; + int ideter2[natomax]; int kko,kok,kkio; long int ii; double xmat=0.0; diff --git a/src/read2.h b/src/read2.h index ce325db..aa4aff3 100644 --- a/src/read2.h +++ b/src/read2.h @@ -33,7 +33,7 @@ typedef struct { int s23a2; int s23b1; int s23b2; - long int postrou; + int postrou; long int fix_trou1; long int fix_trou2; From 900765665f64797b6a4733676165478380aea30b Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Sun, 28 Jan 2018 18:41:28 +0100 Subject: [PATCH 32/38] now calculates ntot automatically --- Makefile | 7 +++++-- src/ex1.c | 6 +++--- src/get_ntot.c | 39 +++++++++++++++++++++++++++++++++++++++ src/get_ntot.h | 3 +++ 4 files changed, 50 insertions(+), 5 deletions(-) create mode 100644 src/get_ntot.c create mode 100644 src/get_ntot.h diff --git a/Makefile b/Makefile index 6624503..f07271e 100644 --- a/Makefile +++ b/Makefile @@ -27,6 +27,9 @@ directories: ${OBJ_DIR} ${LIB_DIR} ${BIN_DIR} ${LIB_DIR}/irpf90.a: directories cd ${SRC_DIR} && irpf90 init && $(MAKE) irpf90.a && cp irpf90.a ../${LIB_DIR} +${OBJ_DIR}/get_ntot.o: ${SRC_DIR}/get_ntot.c directories chkopts + ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + ${OBJ_DIR}/read2.o: ${SRC_DIR}/read2.c directories chkopts ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} @@ -48,6 +51,6 @@ ${OBJ_DIR}/get_val_iaa2.o: ${SRC_DIR}/get_val_iaa2.c directories chkopts ${OBJ_DIR}/ex1.o: ${SRC_DIR}/ex1.c -${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} -${BIN_DIR}/ex1: ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${OBJ_DIR}/ex1.o ${SRC_DIR}/read2.h ${SRC_DIR}/stimsyr.h chkopts - -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3 +${BIN_DIR}/ex1: ${OBJ_DIR}/get_ntot.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${OBJ_DIR}/ex1.o ${SRC_DIR}/read2.h ${SRC_DIR}/get_ntot.h ${SRC_DIR}/stimsyr.h chkopts + -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_ntot.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3 # ${RM} ex1.o read2.o diff --git a/src/ex1.c b/src/ex1.c index 88e86c6..4781021 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -4,6 +4,7 @@ #include "read2.h" #include "stimsyr.h" #include "get_s2.h" +#include "get_ntot.h" #undef __FUNCT__ #define __FUNCT__ "main" @@ -44,6 +45,8 @@ int main(int argc,char **argv) PetscInt nlocal; Data_new(file, &getdata); + getdata.n = get_ntot(getdata.FAM1, getdata.natom, getdata.isz, getdata.ntrou, getdata.fix_trou1, getdata.fix_trou2); + nlocal = getdata.n/getdata.npar; //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); @@ -89,9 +92,6 @@ int main(int argc,char **argv) 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); diff --git a/src/get_ntot.c b/src/get_ntot.c new file mode 100644 index 0000000..631cf65 --- /dev/null +++ b/src/get_ntot.c @@ -0,0 +1,39 @@ +#include "get_ntot.h" + +int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2){ + int tnt1, tnt2; + int natom2; + if(FAM1){ + if(fix_trou1 == fix_trou2){ + natom2 = natom/2; + } + else{ + natom2 = fix_trou2 - fix_trou1; + } + } + else{ + natom2 = natom; + } + + tnt1 = (int)exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))); + int nalpha, nbeta; + + if((((natom-ntrou) + 2*isz) % 2) == 0){ + nalpha=(natom-ntrou+2*isz)/2; + nbeta=(natom -ntrou-2*isz)/2; + if(((natom-ntrou)/2) == isz){ + nbeta=0; + } + } + else{ + nalpha=(natom-ntrou+2*isz+1)/2; + nbeta=(natom -ntrou-2*isz-1)/2; + if(((natom-ntrou+1)/2) == isz){ + nbeta=0; + } + } + + tnt2 = (int) exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1)))); +//printf("nalpha=%d nbeta=%d | | %d %d ntot=%d\n",nalpha, nbeta, tnt1, tnt2, tnt1*tnt2); + return tnt1*tnt2; +} diff --git a/src/get_ntot.h b/src/get_ntot.h new file mode 100644 index 0000000..8cf2069 --- /dev/null +++ b/src/get_ntot.h @@ -0,0 +1,3 @@ +#include + +int get_ntot(_Bool Fam1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2); From 2262d8dce6f18b5324b8d139283009627003f689 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Sun, 28 Jan 2018 18:44:56 +0100 Subject: [PATCH 33/38] removed ntot from read2 --- src/read2.c | 59 +++++++++++++++++++++++++---------------------------- 1 file changed, 28 insertions(+), 31 deletions(-) diff --git a/src/read2.c b/src/read2.c index 9fcd19e..b349799 100644 --- a/src/read2.c +++ b/src/read2.c @@ -13,31 +13,28 @@ 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 != 29){ + if (count != 28){ count++; switch(count){ case 1: - dat->n=atol(line); - break; - case 2: dat->natom=atol(line); break; - case 3: + case 2: dat->nnz=atol(line); break; - case 4: + case 3: dat->npar=atol(line); break; - case 5: + case 4: dat->ntrou=atol(line); break; - case 6: + case 5: dat->isz=atol(line); break; - case 7: + case 6: dat->FAM1 = to_bool(line); break; - case 8: + case 7: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -61,7 +58,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 9: + case 8: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -85,7 +82,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 10: + case 9: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -109,7 +106,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 11: + case 10: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -133,7 +130,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 12: + case 11: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -157,7 +154,7 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 13: + case 12: arrayIdx=0; for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim)) { @@ -181,52 +178,52 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 14: + case 13: dat->nroots=atol(line); break; - case 15: + case 14: dat->s21a1=atol(line); break; - case 16: + case 15: dat->s21a2=atol(line); break; - case 17: + case 16: dat->s21b1=atol(line); break; - case 18: + case 17: dat->s21b2=atol(line); break; - case 19: + case 18: dat->s22a1=atol(line); break; - case 20: + case 19: dat->s22a2=atol(line); break; - case 21: + case 20: dat->s22b1=atol(line); break; - case 22: + case 21: dat->s22b2=atol(line); break; - case 23: + case 22: dat->s23a1=atol(line); break; - case 24: + case 23: dat->s23a2=atol(line); break; - case 25: + case 24: dat->s23b1=atol(line); break; - case 26: + case 25: dat->s23b2=atol(line); break; - case 27: + case 26: dat->postrou=atol(line); break; - case 28: + case 27: dat->fix_trou1=atol(line); break; - case 29: + case 28: dat->fix_trou2=atol(line); break; } /* end of switch */ From aeabef5dedc7b7159668291f140178689fabb66d Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Mon, 29 Jan 2018 00:10:47 +0100 Subject: [PATCH 34/38] looks like now it works for natom > 16 --- src/adr.irp.f | 5 +++-- src/adrnew.irp.f | 5 +++-- src/analyse.irp.f | 5 +++-- src/conv.irp.f | 5 +++-- src/desort.irp.f | 3 ++- src/ex1.c | 24 ++++++++++++------------ src/extra_diag.irp.f | 5 +++-- src/foundet.irp.f | 7 ++++--- src/get_dmat.c | 6 ++---- src/get_dmat.h | 4 ++-- src/get_ntot.c | 7 ++++--- src/get_s2.c | 3 +-- src/get_s2.h | 6 +++--- src/get_s2_cyclic.c | 3 +-- src/get_s2_mov.c | 3 +-- src/getdet.irp.f | 5 +++-- src/main.irp.f | 9 +++++---- src/natomax.irp.f | 4 ++-- src/nt1.irp.f | 5 +++-- src/nt2.irp.f | 3 ++- src/providdet.irp.f | 10 ++++++---- src/searchdet.irp.f | 15 ++++++++------- src/searchdetfull.irp.f | 15 ++++++++------- src/sort.irp.f | 3 ++- src/unit_FIL44.irp.f | 12 ++++++------ src/unit_l1.irp.f | 7 ++++--- 26 files changed, 96 insertions(+), 83 deletions(-) diff --git a/src/adr.irp.f b/src/adr.irp.f index 2a9cd5c..2f79c9d 100644 --- a/src/adr.irp.f +++ b/src/adr.irp.f @@ -1,4 +1,5 @@ subroutine adr(ideter,add) + use iso_c_binding implicit none BEGIN_DOC ! this subroutine provides the address of a detrminant @@ -7,8 +8,8 @@ subroutine adr(ideter,add) ! 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))::deti,dethi,addh,detnew + integer(C_SIZE_T),INTENT(INOUT)::add + integer(C_SIZE_T)::deti,dethi,addh,detnew integer::count,i,j deti=0 diff --git a/src/adrnew.irp.f b/src/adrnew.irp.f index 91ef8c0..229cbc7 100644 --- a/src/adrnew.irp.f +++ b/src/adrnew.irp.f @@ -1,5 +1,6 @@ subroutine adrfull() implicit none + use iso_c_binding BEGIN_DOC ! this subroutine provides the address of a detrminant ! given in old format. @@ -7,8 +8,8 @@ subroutine adrfull() ! 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(C_SIZE_T)::add + integer(C_SIZE_T)::deti,dethi,addh,detnew integer::count,i,j deti=0 diff --git a/src/analyse.irp.f b/src/analyse.irp.f index c621bc2..730a41d 100644 --- a/src/analyse.irp.f +++ b/src/analyse.irp.f @@ -1,10 +1,11 @@ SUBROUTINE ANALYSE(vect, dimvect, startvect, endvect, xymat2, norm2) ! INCLUDE "nbtots.prm" + use iso_c_binding IMPLICIT NONE INTEGER dimvect, nbtots, startvect, endvect REAL*8,dimension(dimvect)::vect - INTEGER (kind=selected_int_kind(16))::add,kvect - INTEGER (kind=selected_int_kind(16))::iaa2,i + INTEGER (C_SIZE_T)::add,kvect + INTEGER (C_SIZE_T)::iaa2,i INTEGER ,dimension(natomax)::ideter INTEGER ,dimension(natomax)::ideter2 REAL*8,allocatable ::xz(:) diff --git a/src/conv.irp.f b/src/conv.irp.f index e01b839..416f1d7 100644 --- a/src/conv.irp.f +++ b/src/conv.irp.f @@ -1,12 +1,13 @@ subroutine conv(ideter,deti,dethi) + use iso_c_binding 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)::deti - integer(kind=selected_int_kind(16)),INTENT(INOUT)::dethi + integer(C_SIZE_T),INTENT(INOUT)::deti + integer(C_SIZE_T),INTENT(INOUT)::dethi integer::i deti=0 dethi=0 diff --git a/src/desort.irp.f b/src/desort.irp.f index 643a7ad..9fb936b 100644 --- a/src/desort.irp.f +++ b/src/desort.irp.f @@ -1,7 +1,8 @@ subroutine desort() + use iso_c_binding implicit none integer::i,j,ord,ordh - integer(kind=selected_int_kind(16))::add,addh,deti,dethi,addt + integer(C_SIZE_T)::add,addh,deti,dethi,addt do i=1,detfound-1 do j=i+1,detfound diff --git a/src/ex1.c b/src/ex1.c index 4781021..0652da0 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -27,9 +27,10 @@ int main(int argc,char **argv) PetscReal normfin2=0.0; PetscReal normfin3=0.0; PetscReal normfin4=0.0; - PetscScalar kr,ki,value[700]; + const int natomax=900; + PetscScalar kr,ki,value[natomax]; Vec xr,xi; - PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol; + PetscInt i,Istart,Iend,col[natomax],maxit,its,nconv,countcol; PetscInt nev, ncv, mpd; PetscLogDouble t1,t2,tt1,tt2; //PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE; @@ -37,7 +38,6 @@ int main(int argc,char **argv) //PetscScalar eigr; //PetscScalar eigi; int mpiid; - int natomax=700; char const* const fileName = argv[1]; FILE* file = fopen(fileName, "r"); @@ -60,8 +60,8 @@ int main(int argc,char **argv) PetscInt kk,ll,mm,nn,iii2,iiii; PetscInt ii; long int iii; - long int tcountcol2,tcol[700],tcountcol[getdata.nnz]; - double val[700]; + long int tcountcol2,tcol[natomax],tcountcol[getdata.nnz]; + double val[natomax]; PetscReal xymat=0.0; PetscReal xymat2=0.0; PetscReal xymat3=0.0; @@ -96,8 +96,8 @@ int main(int argc,char **argv) 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,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 = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,10*getdata.natom,NULL,10*getdata.natom,NULL,&A);CHKERRQ(ierr); + ierr = MatMPIAIJSetPreallocation(A,10*getdata.natom,NULL,10*getdata.natom,NULL);CHKERRQ(ierr); //ierr = MatSetFromOptions(A);CHKERRQ(ierr); //ierr = MatSetUp(A);CHKERRQ(ierr); @@ -149,7 +149,7 @@ int main(int argc,char **argv) 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++){ + for(kk=tcountcol2+tcountcol[ll]+1;kk #include -void get_1rdm(PetscScalar *, PetscInt *, PetscInt *, int *, PetscReal *); -void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, double densmat2[*natom][*natom][*natom][*natom]); +void get_1rdm(PetscScalar *, PetscInt *, PetscInt *, int *, PetscReal *, const int natomax); +void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, double densmat2[*natom][*natom][*natom][*natom], const int natomax); diff --git a/src/get_ntot.c b/src/get_ntot.c index 631cf65..2ddea61 100644 --- a/src/get_ntot.c +++ b/src/get_ntot.c @@ -15,7 +15,8 @@ int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_t natom2 = natom; } - tnt1 = (int)exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))); + tnt1 = (int)ceil(exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1))))); + printf("%10.5f | tnt1=%d\n",exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))),tnt1); int nalpha, nbeta; if((((natom-ntrou) + 2*isz) % 2) == 0){ @@ -33,7 +34,7 @@ int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_t } } - tnt2 = (int) exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1)))); -//printf("nalpha=%d nbeta=%d | | %d %d ntot=%d\n",nalpha, nbeta, tnt1, tnt2, tnt1*tnt2); + tnt2 = (int)ceil(exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1))))); + printf("natom2=%d fix_trou1=%d fix_trou2=%d nalpha=%d nbeta=%d | | %d %d ntot=%d\n",natom2, fix_trou1, fix_trou2, nalpha, nbeta, tnt1, tnt2, tnt1*tnt2); return tnt1*tnt2; } diff --git a/src/get_s2.c b/src/get_s2.c index d882d00..081213e 100644 --- a/src/get_s2.c +++ b/src/get_s2.c @@ -26,8 +26,7 @@ void get_s2(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *n 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; + int *s22b1, int *s22b2, int *s23a1, int *s23a2, int *s23b1, int *s23b2, int *postrou, const int natomax){ long int iaa2, iaa; long int iii; int ideter[natomax]; diff --git a/src/get_s2.h b/src/get_s2.h index 5bf7e13..facd7a2 100644 --- a/src/get_s2.h +++ b/src/get_s2.h @@ -5,7 +5,7 @@ #include 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); + 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); void get_s2_mov(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, int *, @@ -20,7 +20,7 @@ void get_s2_mov(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, int *, int *, int *, - int *); + int *, const int natomax); void get_s2_cyclic(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, int *, @@ -35,4 +35,4 @@ void get_s2_cyclic(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal int *, int *, int *, - int *); + int *, const int natomax); diff --git a/src/get_s2_cyclic.c b/src/get_s2_cyclic.c index ae7a20f..a02bf2d 100644 --- a/src/get_s2_cyclic.c +++ b/src/get_s2_cyclic.c @@ -23,8 +23,7 @@ */ 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; + 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){ long int iaa2, iaa; long int iii; int ideter[natomax]; diff --git a/src/get_s2_mov.c b/src/get_s2_mov.c index a829e77..5470745 100644 --- a/src/get_s2_mov.c +++ b/src/get_s2_mov.c @@ -23,8 +23,7 @@ */ 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; + 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){ long int iaa2, iaa; long int iii; int ideter[natomax]; diff --git a/src/getdet.irp.f b/src/getdet.irp.f index fa98b05..acba73a 100644 --- a/src/getdet.irp.f +++ b/src/getdet.irp.f @@ -1,12 +1,13 @@ subroutine getdet(add,ideter) + use iso_c_binding 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(C_SIZE_T),INTENT(IN)::add + integer(C_SIZE_T)::deta,detb integer::i,const,ia,ib, natom2 ib = MOD(add,nt2) diff --git a/src/main.irp.f b/src/main.irp.f index 9f26adb..307906a 100644 --- a/src/main.irp.f +++ b/src/main.irp.f @@ -1,13 +1,14 @@ program main + use iso_c_binding implicit none - integer(kind=selected_int_kind(16)),allocatable::tl1(:),tl2(:),tktyp(:) + integer(C_SIZE_T),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 + integer(C_SIZE_T),allocatable::tcol(:) + integer(C_SIZE_T),dimension(10)::tcountcol + integer(C_SIZE_T)::tistart allocate (tl1(maxlien)) allocate (tl2(maxlien)) allocate (tktyp(maxlien)) diff --git a/src/natomax.irp.f b/src/natomax.irp.f index c4c53ce..6f443a9 100644 --- a/src/natomax.irp.f +++ b/src/natomax.irp.f @@ -13,10 +13,10 @@ BEGIN_PROVIDER [integer, natomax] END_DOC implicit none - natomax=700 + natomax=900 jrangmax=10705432 maxial=20 - maxlien=700 + maxlien=900 maxplac=20 maxdet=10000 END_PROVIDER diff --git a/src/nt1.irp.f b/src/nt1.irp.f index 11dda96..55e7971 100644 --- a/src/nt1.irp.f +++ b/src/nt1.irp.f @@ -1,9 +1,10 @@ -BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt1] +use iso_c_binding +BEGIN_PROVIDER [integer(C_SIZE_T), nt1] BEGIN_DOC ! calculates the number of det the 3's moving END_DOC implicit none - integer(kind=selected_int_kind(16))::natom2 + integer(C_SIZE_T)::natom2 ! call combin(idet1(1,nt1+1),natom,ntrou,nt1,32,jrangmax) natom2=natom diff --git a/src/nt2.irp.f b/src/nt2.irp.f index 05e31c2..7becdac 100644 --- a/src/nt2.irp.f +++ b/src/nt2.irp.f @@ -1,4 +1,5 @@ -BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt2] +use iso_c_binding +BEGIN_PROVIDER [integer(C_SIZE_T), nt2] BEGIN_DOC ! calculates the number of det the 1's moving END_DOC diff --git a/src/providdet.irp.f b/src/providdet.irp.f index 577890a..d651244 100644 --- a/src/providdet.irp.f +++ b/src/providdet.irp.f @@ -1,12 +1,14 @@ -BEGIN_PROVIDER[integer(kind=selected_int_kind(16)),det,(nt2,2)] -&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)),deth,(nt1,2)] +use iso_c_binding +BEGIN_PROVIDER[integer(C_SIZE_T),det,(nt2,2)] +&BEGIN_PROVIDER[integer(C_SIZE_T),deth,(nt1,2)] BEGIN_DOC ! provides det and deth array END_DOC + use iso_c_binding 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(C_SIZE_t)::a + integer(C_SIZE_T)::i,count integer::const i=1 a=0 diff --git a/src/searchdet.irp.f b/src/searchdet.irp.f index f8ea9b2..66e9b6d 100644 --- a/src/searchdet.irp.f +++ b/src/searchdet.irp.f @@ -1,17 +1,18 @@ subroutine searchdet(deti,add,dethi,addh) + use iso_c_binding 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)::deti - integer(kind=selected_int_kind(16)),INTENT(INOUT)::add - 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,j + integer(C_SIZE_T),INTENT(INOUT)::deti + integer(C_SIZE_T),INTENT(INOUT)::add + integer(C_SIZE_T),INTENT(INOUT)::dethi + integer(C_SIZE_T),INTENT(INOUT)::addh + integer(C_SIZE_T)::dethsh + integer(C_SIZE_T)::a + integer(C_SIZE_T)::i,j integer::count logical::found diff --git a/src/searchdetfull.irp.f b/src/searchdetfull.irp.f index 362c527..0906e47 100644 --- a/src/searchdetfull.irp.f +++ b/src/searchdetfull.irp.f @@ -1,17 +1,18 @@ subroutine searchdetfull() + use iso_c_binding 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(C_SIZE_T),INTENT(INOUT)::foundetadr(maxlien,4) + integer(C_SIZE_T)::add +! integer(C_SIZE_T),INTENT(INOUT)::deth + integer(C_SIZE_T)::dethsh + integer(C_SIZE_T)::addh + integer(C_SIZE_T)::a + integer(C_SIZE_T)::i integer::const,count i=1 a=0 diff --git a/src/sort.irp.f b/src/sort.irp.f index 14c7407..6a1c1f2 100644 --- a/src/sort.irp.f +++ b/src/sort.irp.f @@ -1,7 +1,8 @@ subroutine sort() + use iso_c_binding implicit none integer::i,j,ord,ordh - integer(kind=selected_int_kind(16))::add,addh,deti,dethi,addt + integer(C_SIZE_T)::add,addh,deti,dethi,addt do i=1,detfound-1 do j=i+1,detfound diff --git a/src/unit_FIL44.irp.f b/src/unit_FIL44.irp.f index 647554b..96da0c8 100644 --- a/src/unit_FIL44.irp.f +++ b/src/unit_FIL44.irp.f @@ -3,19 +3,19 @@ BEGIN_DOC ! file units for writing END_DOC - + use iso_c_binding 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(C_SIZE_T)::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 + integer(C_SIZE_T),dimension(maxlien) ::tl1,tl2,tktyp + integer(C_SIZE_T),dimension(nrows)::tcountcol + integer(C_SIZE_T)::tistart real*8,dimension(maxlien)::tval - integer(kind=selected_int_kind(16)),dimension(maxlien)::tcol + integer(C_SIZE_T),dimension(maxlien)::tcol real*8 :: xmat integer :: ik,imat4,iaa2,iik integer :: ik1,ik2,jmat4,IC,ikmax,ikmin diff --git a/src/unit_l1.irp.f b/src/unit_l1.irp.f index cb61e46..b9d010b 100644 --- a/src/unit_l1.irp.f +++ b/src/unit_l1.irp.f @@ -13,6 +13,7 @@ tfix_trou2, & tfam1, & tcol,tval) + use iso_c_binding implicit none integer,INTENT(INOUT)::tistart, tnrows integer,INTENT(INOUT)::tntrou, tisz @@ -20,9 +21,9 @@ logical*1,INTENT(INOUT)::tfam1 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) + integer(C_SIZE_T),INTENT(INOUT)::tcol(maxlien) + integer(C_SIZE_T),INTENT(INOUT),dimension(tnrows)::tcountcol + integer(C_SIZE_T),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien) real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien) nrows = tnrows From c8e6e61603f046828ffc11b4329b27881bc71c8b Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Mon, 29 Jan 2018 16:21:56 +0100 Subject: [PATCH 35/38] many new features: 1. added wf printing support 2. fixed problem with factorial calculation --- src/elem_diag.irp.f | 4 +++ src/ex1.c | 58 +++---------------------------- src/get_ntot.c | 8 ++--- src/get_ntot.h | 1 + src/provide_clust.irp.f | 1 + src/read2.c | 77 +++++++++++++++++++++++++++++------------ src/read2.h | 16 +++++---- src/stimsyr.h | 1 + src/unit_l1.irp.f | 4 ++- 9 files changed, 83 insertions(+), 87 deletions(-) diff --git a/src/elem_diag.irp.f b/src/elem_diag.irp.f index 679b0d2..85b244f 100644 --- a/src/elem_diag.irp.f +++ b/src/elem_diag.irp.f @@ -28,6 +28,10 @@ ! if(yw)write(6,*)iaa,'diag,v1' ! endif enddo + do i=1, natom + if(deter(i).ne.3) xmatd = xmatd + E(i) + enddo + xmatd = xmatd - E(natom+1) !-----stockage de l element diag diff --git a/src/ex1.c b/src/ex1.c index 0652da0..e062938 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -33,10 +33,7 @@ int main(int argc,char **argv) PetscInt i,Istart,Iend,col[natomax],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]; @@ -44,16 +41,14 @@ int main(int argc,char **argv) Data getdata; PetscInt nlocal; + /* gather the input data */ Data_new(file, &getdata); getdata.n = get_ntot(getdata.FAM1, getdata.natom, getdata.isz, getdata.ntrou, getdata.fix_trou1, getdata.fix_trou2); nlocal = getdata.n/getdata.npar; -//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; PetscInt indxr[nlocal]; -//Vec Vr,Vi; char filename[PETSC_MAX_PATH_LEN]="FIL666"; PetscViewer viewer; PetscBool ishermitian; @@ -98,30 +93,18 @@ int main(int argc,char **argv) ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,10*getdata.natom,NULL,10*getdata.natom,NULL,&A);CHKERRQ(ierr); ierr = MatMPIAIJSetPreallocation(A,10*getdata.natom,NULL,10*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 | %d\n",Istart, Iend,getdata.natom);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; i +#include int get_ntot(_Bool Fam1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2); diff --git a/src/provide_clust.irp.f b/src/provide_clust.irp.f index 9ed9988..4f7e96e 100644 --- a/src/provide_clust.irp.f +++ b/src/provide_clust.irp.f @@ -4,6 +4,7 @@ BEGIN_PROVIDER[integer,l1, (maxlien)] &BEGIN_PROVIDER[real*8, xtt ,(maxlien)] &BEGIN_PROVIDER[real*8, xjjz ,(maxlien)] &BEGIN_PROVIDER[real*8, xjjxy,(maxlien)] +&BEGIN_PROVIDER[real*8, E,(maxlien)] &BEGIN_PROVIDER[integer, ntrou] &BEGIN_PROVIDER[integer, isz] &BEGIN_PROVIDER[logical*1, FAM1] diff --git a/src/read2.c b/src/read2.c index b349799..8f8b040 100644 --- a/src/read2.c +++ b/src/read2.c @@ -11,9 +11,9 @@ void Data_new(FILE* file, Data* dat) { 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 != 28){ + /* note that fgets doesn't strip the terminating \n, checking its + presence would allow to handle lines longer than sizeof(line) */ + if (count != 30){ count++; switch(count){ case 1: @@ -178,54 +178,84 @@ void Data_new(FILE* file, Data* dat) { } } break; - case 13: - dat->nroots=atol(line); + case 13: + 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->E[arrayIdx++] = val; + } + } break; case 14: - dat->s21a1=atol(line); + dat->nroots=atol(line); break; case 15: - dat->s21a2=atol(line); + dat->s21a1=atol(line); break; case 16: - dat->s21b1=atol(line); + dat->s21a2=atol(line); break; case 17: - dat->s21b2=atol(line); + dat->s21b1=atol(line); break; case 18: - dat->s22a1=atol(line); + dat->s21b2=atol(line); break; case 19: - dat->s22a2=atol(line); + dat->s22a1=atol(line); break; case 20: - dat->s22b1=atol(line); + dat->s22a2=atol(line); break; case 21: - dat->s22b2=atol(line); + dat->s22b1=atol(line); break; case 22: - dat->s23a1=atol(line); + dat->s22b2=atol(line); break; case 23: - dat->s23a2=atol(line); + dat->s23a1=atol(line); break; case 24: - dat->s23b1=atol(line); + dat->s23a2=atol(line); break; case 25: - dat->s23b2=atol(line); + dat->s23b1=atol(line); break; case 26: - dat->postrou=atol(line); + dat->s23b2=atol(line); break; case 27: - dat->fix_trou1=atol(line); + dat->postrou=atol(line); break; case 28: + dat->fix_trou1=atol(line); + break; + case 29: dat->fix_trou2=atol(line); break; + case 30: + dat->print_wf = atol(line); + break; + default: + printf("Done reading file\n"); + break; } /* end of switch */ } /* end of the input file */ @@ -236,10 +266,13 @@ void Data_new(FILE* file, Data* dat) { } _Bool to_bool(const char* str) { - PetscBool strflg; + PetscBool strflg=PETSC_FALSE; PetscStrcmp("true\n",str, &strflg); - if(!strflg) PetscStrcmp("True\n",str, &strflg); - if(!strflg) PetscStrcmp("TRUE\n",str, &strflg); + if(strflg != PETSC_TRUE) PetscStrcmp("True\n",str, &strflg); + if(strflg != PETSC_TRUE) PetscStrcmp("TRUE\n",str, &strflg); + if(strflg != PETSC_TRUE) PetscStrcmp("true",str, &strflg); + if(strflg != PETSC_TRUE) PetscStrcmp("True",str, &strflg); + if(strflg != PETSC_TRUE) PetscStrcmp("TRUE",str, &strflg); return strflg; } diff --git a/src/read2.h b/src/read2.h index aa4aff3..9a1a6ff 100644 --- a/src/read2.h +++ b/src/read2.h @@ -12,13 +12,14 @@ typedef struct { PetscInt n; long int nnz,npar; long int ntrou,isz; - _Bool FAM1; - long int l1[700]; - long int l2[700]; - long int ktyp[700]; - double xjjz[700]; - double xjjxy[700]; - double xtt[700]; + _Bool FAM1; + long int l1[900]; + long int l2[900]; + long int ktyp[900]; + double xjjz[900]; + double xjjxy[900]; + double xtt[900]; + double E[900]; long int nroots; int natom; int s21a1; @@ -36,6 +37,7 @@ typedef struct { int postrou; long int fix_trou1; long int fix_trou2; + long int print_wf; } Data ; diff --git a/src/stimsyr.h b/src/stimsyr.h index 2397381..5964517 100644 --- a/src/stimsyr.h +++ b/src/stimsyr.h @@ -9,6 +9,7 @@ void unit_l1_( double *, double *, double *, + double *, long int *, long int *, long int *, diff --git a/src/unit_l1.irp.f b/src/unit_l1.irp.f index b9d010b..fc0dc40 100644 --- a/src/unit_l1.irp.f +++ b/src/unit_l1.irp.f @@ -6,6 +6,7 @@ txjjxy, & txjjz , & txtt , & + tE , & tcountcol, & tntrou, & tisz, & @@ -24,7 +25,7 @@ integer(C_SIZE_T),INTENT(INOUT)::tcol(maxlien) integer(C_SIZE_T),INTENT(INOUT),dimension(tnrows)::tcountcol integer(C_SIZE_T),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien) - real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien) + real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien), tE(maxlien) nrows = tnrows provide nrows @@ -35,6 +36,7 @@ xtt(i) = txtt(i) xjjxy(i) = txjjxy(i) xjjz (i) = txjjz (i) + E (i) = tE (i) enddo ntrou = tntrou isz = tisz From ac24dd6dc1901ac5250eed11bc2201241f19e75e Mon Sep 17 00:00:00 2001 From: vijay Date: Fri, 21 Feb 2020 14:26:29 +0100 Subject: [PATCH 36/38] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 48e92e9..df2a7b8 100644 --- a/README.md +++ b/README.md @@ -50,7 +50,7 @@ _Using DEHam_ 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) +1,1,1,2,2,2,2,3,3,3 # third line contains the type of link (1 for t or J, 2 for K and 3 for none) .1430,-0.20,0.0000 # The three types of links this line gives J, K .1430,-0.20,0.0000 # -1.00,0.0,0.00 # This line gives t From 20390f5ad525caa00ca679a6c4ed373a2b40e40d Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Thu, 5 Mar 2020 15:25:26 +0100 Subject: [PATCH 37/38] Fixed issues with OSX and int to real conversion (nt1 and nt2) now compiles and runs on osx. Added two examples --- Makefile | 2 +- examples/small.inp | 14 ++++++++++++++ examples/three.inp | 14 ++++++++++++++ src/ex1.c | 7 +++---- src/nt1.irp.f | 2 +- src/nt2.irp.f | 2 +- 6 files changed, 34 insertions(+), 7 deletions(-) create mode 100644 examples/small.inp create mode 100644 examples/three.inp diff --git a/Makefile b/Makefile index f07271e..f34b0b9 100644 --- a/Makefile +++ b/Makefile @@ -52,5 +52,5 @@ ${OBJ_DIR}/ex1.o: ${SRC_DIR}/ex1.c -${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} ${BIN_DIR}/ex1: ${OBJ_DIR}/get_ntot.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${OBJ_DIR}/ex1.o ${SRC_DIR}/read2.h ${SRC_DIR}/get_ntot.h ${SRC_DIR}/stimsyr.h chkopts - -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_ntot.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3 + -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.o ${OBJ_DIR}/get_ntot.o ${OBJ_DIR}/get_s2.o ${OBJ_DIR}/get_s2_mov.o ${OBJ_DIR}/get_s2_cyclic.o ${OBJ_DIR}/get_dmat.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB} # ${RM} ex1.o read2.o diff --git a/examples/small.inp b/examples/small.inp new file mode 100644 index 0000000..4c0c4ae --- /dev/null +++ b/examples/small.inp @@ -0,0 +1,14 @@ +8 +140 +1 +1 +0 +true +1,2,3,1,2,3,4,5,6,7 +2,3,4,8,7,6,5,6,7,8 +1,1,1,2,2,2,2,3,3,3 +.1430,-0.20,0.0000 +.1430,-0.20,0.0000 +-1.00,0.0,0.00 +1 +1 diff --git a/examples/three.inp b/examples/three.inp new file mode 100644 index 0000000..8706cf2 --- /dev/null +++ b/examples/three.inp @@ -0,0 +1,14 @@ +3 +1 +1 +1 +1 +true +1,2,1,2,3,4,5 +2,3,6,5,4,5,6 +1,1,2,2,2,3,3 +.1430,-0.20,0.0000 +.1430,-0.20,0.0000 +-1.00,0.0,0.00 +1 +1 diff --git a/src/ex1.c b/src/ex1.c index e062938..c2a2659 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -123,9 +123,9 @@ int main(int argc,char **argv) &getdata.FAM1, tcol, val); - if(i%getdata.npar == 0 && mpiid==0){ - ierr = PetscPrintf(PETSC_COMM_WORLD," i: %d \n",i);CHKERRQ(ierr); - } +// if(i%getdata.npar == 0 && mpiid==0){ +// ierr = PetscPrintf(PETSC_COMM_WORLD," i: %d \n",i);CHKERRQ(ierr); +// } for(ll=0;ll