diff --git a/src/get_proj.c b/src/get_proj.c new file mode 100644 index 0000000..7985c95 --- /dev/null +++ b/src/get_proj.c @@ -0,0 +1,193 @@ +#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(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, int iroot, double *projvec, const int natomax){ + + int ideter[natomax]; + int ideter2[natomax]; + int kko,kok,kkio; + int mmo,mom,mmio; + int p, q; + long int ii; + PetscInt iiii; + long int iii; + long int iaa2, iaa; + long int nrow=-1, ncol=-1; + double ms=0.0, rest=0.0, norm=0.0; + double sq2; + double normvec[6]; + int sumMs=0; + int ntrouGauch = 0; + int ntrouDroit = 0; + double projvecmat[3*3*6]; + for(ii=0;ii<3*3*6;++ii) projvecmat[ii] = 0.0; + sq2 = sqrt(2.0); + + for(ii=*Istart;ii<*Iend;ii++) { + iii = ii + 1; + iiii = ii; + getdet_(&iii, ideter); + + ms = 0.0; + p = -1; + q = -1; + // Find if Ms=5/2 + sumMs = 0; + ntrouGauch = 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){ + ntrouDroit += 1; + q = kko - 3; + } + } + if(ntrouGauch == 1 && ntrouDroit == 1){ + sumMs = 1; + } + else{ + sumMs = 0; + } + for(kko=0;kko<=2;++kko){ + if(sumMs==1){ + if(ideter[kko]==1){ + ms = ms + 0.5; + } + else if(ideter[kko]==2){ + ms = ms - 0.5; + } + } + } + for(kok= 11;kok>=9;--kok){ + if(sumMs==1){ + if(ideter[kok]==1){ + ms = ms + 0.5; + } + else if(ideter[kok]==2){ + ms = ms - 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.h b/src/get_proj.h new file mode 100644 index 0000000..8b226b5 --- /dev/null +++ b/src/get_proj.h @@ -0,0 +1,8 @@ +#include +#include +#include +#include +#include +#include + +void get_proj(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, int iroot, double *projvec, const int natomax);