From 6df166e7191632cdc295b4e64eeacc4e4b0f76f4 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Thu, 23 Feb 2023 18:17:30 +0100 Subject: [PATCH] Working on 9_3h. --- Makefile | 9 +- src/ex1.c | 5 +- src/get_proj_9_3h.c | 224 ++++++++++++++++++++++++++++++++++++++++++++ src/get_proj_9_3h.h | 8 ++ 4 files changed, 241 insertions(+), 5 deletions(-) create mode 100644 src/get_proj_9_3h.c create mode 100644 src/get_proj_9_3h.h diff --git a/Makefile b/Makefile index 49fb7c4..7ad10c8 100644 --- a/Makefile +++ b/Makefile @@ -45,7 +45,10 @@ ${OBJ_DIR}/get_s2_mov.o: ${SRC_DIR}/get_s2_mov.c directories chkopts ${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_proj.o: ${SRC_DIR}/get_proj.c directories chkopts +${OBJ_DIR}/get_proj_9_3h.o: ${SRC_DIR}/get_proj_9_3h.c directories chkopts + ${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB} + +${OBJ_DIR}/get_proj_9_3h.o: ${SRC_DIR}/get_proj_9_3h.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 @@ -54,6 +57,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}/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_proj.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_proj.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${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_proj_9_3h.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_proj_9_3h.o ${OBJ_DIR}/get_val_iaa2.o ${LIB_DIR}/irpf90.a ${SLEPC_EPS_LIB} # ${RM} ex1.o read2.o diff --git a/src/ex1.c b/src/ex1.c index ac16b83..6eaa86f 100644 --- a/src/ex1.c +++ b/src/ex1.c @@ -5,7 +5,8 @@ #include "stimsyr.h" #include "get_s2.h" #include "get_ntot.h" -#include "get_proj.h" +//#include "get_proj.h" +#include "get_proj_9_3h.h" #undef __FUNCT__ #define __FUNCT__ "main" @@ -264,7 +265,7 @@ int main(int argc,char **argv) get_1rdm(values, &Istart, &Iend, &getdata.natom, &trace1rdm, natomax); //get_2rdm(values, &Istart, &Iend, &getdata.natom, &trace2rdm, densmat2, natomax); //get_2rdm(values, &Istart, &Iend, &getdata.natom, &trace2rdm, natomax); - get_proj(values, &Istart, &Iend, &getdata.natom, i, projvec, natomax); + get_proj_9_3h(values, &Istart, &Iend, &getdata.natom, i, projvec, natomax); weightproj = 0.0; for(ii=0;ii<6;++ii) weightproj += projvec[(i)*6 + ii]*projvec[(i)*6+ii]; // analyse_(valxr, (Iend-Istart), &Istart, &Iend, &xymat, &norm); diff --git a/src/get_proj_9_3h.c b/src/get_proj_9_3h.c new file mode 100644 index 0000000..89e571c --- /dev/null +++ b/src/get_proj_9_3h.c @@ -0,0 +1,224 @@ +#include +#include +#include +#include +#include +#include +#include "get_dmat.h" + + +/* + * + *---------------------------------------- + * Calculate the projection for 6_2h + *---------------------------------------- + * Input + * ===== + * valxr = The full vector + * Istart = Local starting id + * Iend = Local ending id + * Output + * ===== + * projvec + */ +void get_proj_9_3h(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, int iroot, double *projvec, const int natomax){ + + int ideter[natomax]; + int ideter1[natomax]; + int kko,kok,kkio; + int mmo,mom,mmio; + int p, q, r; + long int ii; + PetscInt iiii; + long int iii; + long int iaa2, iaa; + long int nrow=-1, ncol=-1; + double ms1=0.0, ms2=0.0 ,rest=0.0, norm=0.0; + double sq2; + double normvec[21]; + int sumMs=0; + int ntrouGauch = 0; + int ntrouMilieu = 0; + int ntrouDroit = 0; + double projvecmat[3*3*3*21]; + for(ii=0;ii<3*3*3*21;++ii) projvecmat[ii] = 0.0; + sq2 = sqrt(2.0); + + for(ii=*Istart;ii<*Iend;ii++) { + iii = ii + 1; + iiii = ii; + getdet_(&iii, ideter); + + ms1 = 0.0; + ms2 = 0.0; + p = -1; + q = -1; + // Find if Ms=5/2 + sumMs = 0; + ntrouGauch = 0; + ntrouMilieu = 0; + ntrouDroit = 0; + for(kko=0;kko<3;++kko){ + if(ideter[kko]==3){ + ntrouGauch += 1; + p = kko; + } + } + for(kko=3;kko<6;++kko){ + if(ideter[kko]==3){ + ntrouMilieu += 1; + q = kko - 3; + } + } + for(kko=6;kko<9;++kko){ + if(ideter[kko]==3){ + ntrouDroit += 1; + r = kko - 3; + } + } + if(ntrouGauch == 1 && ntrouMilieu == 1 && ntrouDroit == 1){ + sumMs = 1; + } + else{ + sumMs = 0; + } + // First box ms + for(kko=0;kko<=2;++kko){ + if(sumMs==1){ + if(ideter[kko]==1){ + ms1 = ms1 + 0.5; + } + else if(ideter[kko]==2){ + ms1 = ms1 - 0.5; + } + } + } + for(kok= 11;kok>=9;--kok){ + if(sumMs==1){ + if(ideter[kok]==1){ + ms1 = ms1 + 0.5; + } + else if(ideter[kok]==2){ + ms1 = ms1 - 0.5; + } + } + } + // Second box ms + for(kko=0;kko<=2;++kko){ + if(sumMs==1){ + if(ideter[kko]==1){ + ms2 = ms2 + 0.5; + } + else if(ideter[kko]==2){ + ms2 = ms2 - 0.5; + } + } + } + for(kok= 11;kok>=9;--kok){ + if(sumMs==1){ + if(ideter[kok]==1){ + ms2 = ms2 + 0.5; + } + else if(ideter[kok]==2){ + ms2 = ms2 - 0.5; + } + } + } + if(ideter[1] ==3 && ideter[4] == 3){ + norm += valxr[iiii]*valxr[iiii]; + } + + if(fabs(ms-2.5) < 1e-10){ + projvecmat[6*3*p + 6*q + 0] = valxr[iiii]; + normvec[0] += 1.0; + } + else if(fabs(ms+2.5) < 1e-10){ + projvecmat[6*3*p + 6*q + 1] = valxr[iiii]; + normvec[1] += 1.0; + } + else if(fabs(ms-1.5) < 1e-10){ + projvecmat[6*3*p + 6*q + 2] += valxr[iiii]; + normvec[2] += 1.0; + } + else if(fabs(ms+1.5) < 1e-10){ + projvecmat[6*3*p + 6*q + 3] += valxr[iiii]; + normvec[3] += 1.0; + } + else if(fabs(ms-0.5) < 1e-10){ + projvecmat[6*3*p + 6*q + 4] += valxr[iiii]; + normvec[4] += 1.0; + } + else if(fabs(ms+0.5) < 1e-10){ + projvecmat[6*3*p + 6*q + 5] += valxr[iiii]; + normvec[5] += 1.0; + } + else{ + rest+=valxr[iiii]; + } + + ms = 0.0; + + } + + for(p=0;p<3;++p){ + for(q=0;q<3;++q) { + projvecmat[p*3*6 + q*6 + 0] = projvecmat[p*3*6 + q*6 + 0]/sqrt(1.0); + projvecmat[p*3*6 + q*6 + 1] = projvecmat[p*3*6 + q*6 + 1]/sqrt(1.0); + projvecmat[p*3*6 + q*6 + 2] = projvecmat[p*3*6 + q*6 + 2]/sqrt(25.0); + projvecmat[p*3*6 + q*6 + 3] = projvecmat[p*3*6 + q*6 + 3]/sqrt(25.0); + projvecmat[p*3*6 + q*6 + 4] = projvecmat[p*3*6 + q*6 + 4]/sqrt(100.0); + projvecmat[p*3*6 + q*6 + 5] = projvecmat[p*3*6 + q*6 + 5]/sqrt(100.0); + } + } + //projvec[(iroot)*6 + 0] = projvec[(iroot)*6 + 0]/sqrt(1.0); + //projvec[(iroot)*6 + 1] = projvec[(iroot)*6 + 1]/sqrt(1.0); + //projvec[(iroot)*6 + 2] = projvec[(iroot)*6 + 2]/sqrt(25.0); + //projvec[(iroot)*6 + 3] = projvec[(iroot)*6 + 3]/sqrt(25.0); + //projvec[(iroot)*6 + 4] = projvec[(iroot)*6 + 4]/sqrt(100.0); + //projvec[(iroot)*6 + 5] = projvec[(iroot)*6 + 5]/sqrt(100.0); + //for(ii=0;ii<6;++ii){ + // projvec[(iroot)*6 + ii] = \ + // projvecmat[1*3*6 + 1*6 + ii]; + //} + for(ii=0;ii<6;++ii){ + projvecmat[0*3*6 + 0*6 + ii] *= 1.0/4.0; // (1,1) = 1.0/4.0 + projvecmat[0*3*6 + 1*6 + ii] *= sq2/4.0; // (1,2) = 1.4/4.0 + projvecmat[0*3*6 + 2*6 + ii] *= 1.0/4.0; // (1,3) = 1.0/4.0 + + projvecmat[1*3*6 + 0*6 + ii] *= sq2/4.0; // (2,1) = 1.4/4.0 + projvecmat[1*3*6 + 1*6 + ii] *= 2.0/4.0; // (2,2) = 2.0/4.0 + projvecmat[1*3*6 + 2*6 + ii] *= sq2/4.0; // (2,3) = 1.4/4.0 + + projvecmat[2*3*6 + 0*6 + ii] *= 1.0/4.0; // (3,1) = 1.0/4.0 + projvecmat[2*3*6 + 1*6 + ii] *= sq2/4.0; // (3,2) = 1.4/4.0 + projvecmat[2*3*6 + 2*6 + ii] *= 1.0/4.0; // (3,3) = 1.0/4.0 + } + + for(ii=0;ii<6;++ii){ + projvec[(iroot)*6 + ii] = \ + projvecmat[0*3*6 + 0*6 + ii] \ + + projvecmat[0*3*6 + 1*6 + ii] \ + + projvecmat[0*3*6 + 2*6 + ii] \ + + projvecmat[1*3*6 + 0*6 + ii] \ + + projvecmat[1*3*6 + 1*6 + ii] \ + + projvecmat[1*3*6 + 2*6 + ii] \ + + projvecmat[2*3*6 + 0*6 + ii] \ + + projvecmat[2*3*6 + 1*6 + ii] \ + + projvecmat[2*3*6 + 2*6 + ii]; + } + //projvec[(iroot)*6 + 0] = projvec[(iroot)*6 *6+ 0]*6/sqrt(9.0); + //projvec[(iroot)*6 + 1] = projvec[(iroot)*6 + 1]/sqrt(9.0); + //projvec[(iroot)*6 + 2] = projvec[(iroot)*6 + 2]/sqrt(25.0*9); + //projvec[(iroot)*6 + 3] = projvec[(iroot)*6 + 3]/sqrt(25.0*9); + //projvec[(iroot)*6 + 4] = projvec[(iroot)*6 + 4]/sqrt(100.0*9); + //projvec[(iroot)*6 + 5] = projvec[(iroot)*6 + 5]/sqrt(100.0*9); + //printf(" norm = %4.4f projnorm = %4.4f\n",norm,projvec[5]*projvec[5] +projvec[4]*projvec[4] +projvec[3]*projvec[3] +projvec[2]*projvec[2] +projvec[1]*projvec[1] +projvec[0]*projvec[0]); + //printf(" rest=%4.4f norm1=%2.2f norm2=%2.2f norm3=%2.2f norm4=%2.2f norm5=%2.2f norm6=%2.2f\n",rest,normvec[0],normvec[1],normvec[2],normvec[3],normvec[4],normvec[5]); + //printf(" proj ( 5/2,-5/2) : 1,4 %2.3f 1,5 %2.3f 1,6 %2.3f 2,4 %2.3f 2,5 %2.3f 2,6 %2.3f 3,4 %2.3f 3,5 %2.3f 3,6 %2.3f \n",projvecseparate[0], projvecseparate[1], projvecseparate[2], projvecseparate[3], projvecseparate[4], projvecseparate[5], projvecseparate[6], projvecseparate[7], projvecseparate[8]); + //printf(" proj (-5/2, 5/2) : 1,4 %2.3f 1,5 %2.3f 1,6 %2.3f 2,4 %2.3f 2,5 %2.3f 2,6 %2.3f 3,4 %2.3f 3,5 %2.3f 3,6 %2.3f \n",projvecseparate[0+9], projvecseparate[1+9], projvecseparate[2+9], projvecseparate[3+9], projvecseparate[4+9], projvecseparate[5+9], projvecseparate[6+9], projvecseparate[7+9], projvecseparate[8+9]); + //printf(" proj ( 3/2,-3/2) : 1,4 %2.3f 1,5 %2.3f 1,6 %2.3f 2,4 %2.3f 2,5 %2.3f 2,6 %2.3f 3,4 %2.3f 3,5 %2.3f 3,6 %2.3f \n",projvecseparate[0+18], projvecseparate[1+18], projvecseparate[2+18], projvecseparate[3+18], projvecseparate[4+18], projvecseparate[5+18], projvecseparate[6+18], projvecseparate[7+18], projvecseparate[8+18]); + //printf(" proj (-3/2, 3/2) : 1,4 %2.3f 1,5 %2.3f 1,6 %2.3f 2,4 %2.3f 2,5 %2.3f 2,6 %2.3f 3,4 %2.3f 3,5 %2.3f 3,6 %2.3f \n",projvecseparate[0+27], projvecseparate[1+27], projvecseparate[2+27], projvecseparate[3+27], projvecseparate[4+27], projvecseparate[5+27], projvecseparate[6+27], projvecseparate[7+27], projvecseparate[8+27]); + //printf(" proj ( 1/2,-1/2) : 1,4 %2.3f 1,5 %2.3f 1,6 %2.3f 2,4 %2.3f 2,5 %2.3f 2,6 %2.3f 3,4 %2.3f 3,5 %2.3f 3,6 %2.3f \n",projvecseparate[0+36], projvecseparate[1+36], projvecseparate[2+36], projvecseparate[3+36], projvecseparate[4+36], projvecseparate[5+36], projvecseparate[6+36], projvecseparate[7+36], projvecseparate[8+36]); + //printf(" proj (-1/2, 1/2) : 1,4 %2.3f 1,5 %2.3f 1,6 %2.3f 2,4 %2.3f 2,5 %2.3f 2,6 %2.3f 3,4 %2.3f 3,5 %2.3f 3,6 %2.3f \n",projvecseparate[0+45], projvecseparate[1+45], projvecseparate[2+45], projvecseparate[3+45], projvecseparate[4+45], projvecseparate[5+45], projvecseparate[6+45], projvecseparate[7+45], projvecseparate[8+45]); +} /** END **/ + diff --git a/src/get_proj_9_3h.h b/src/get_proj_9_3h.h new file mode 100644 index 0000000..2612138 --- /dev/null +++ b/src/get_proj_9_3h.h @@ -0,0 +1,8 @@ +#include +#include +#include +#include +#include +#include + +void get_proj_9_3h(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, int iroot, double *projvec, const int natomax);