LCPQ/DEHam
2
0
Fork 0
mirror of https://github.com/LCPQ/DEHam synced 2024-07-30 09:04:19 +02:00
Code Issues Releases Wiki Activity

Working on 9_3h.

Browse Source
This commit is contained in:
v1j4y 2023-02-23 18:17:30 +01:00
parent 0398124bfe
commit 6df166e719
4 changed files with 241 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Makefile
View File

@ -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

5
src/ex1.c
View File

@ -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);

224
src/get_proj_9_3h.c Normal file
View File

@ -0,0 +1,224 @@
#include <stdio.h>
#include <petsctime.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#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 **/

8
src/get_proj_9_3h.h Normal file
View File

@ -0,0 +1,8 @@
#include <stdio.h>
#include <petsctime.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
void get_proj_9_3h(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, int iroot, double *projvec, const int natomax);