2
0
mirror of https://github.com/LCPQ/DEHam synced 2024-07-30 09:04:19 +02:00

Projection working for 6_2h.

This commit is contained in:
v1j4y 2023-02-07 22:00:00 +01:00
parent 51d2fb62ec
commit 0398124bfe
8 changed files with 86 additions and 37 deletions

View File

@ -45,12 +45,15 @@ ${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
${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}/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}
${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}
# ${RM} ex1.o read2.o

View File

@ -5,6 +5,7 @@
#include "stimsyr.h"
#include "get_s2.h"
#include "get_ntot.h"
#include "get_proj.h"
#undef __FUNCT__
#define __FUNCT__ "main"
@ -84,13 +85,18 @@ int main(int argc,char **argv)
//PetscInt idx_to[nlocal], idx_from[nlocal];
PetscScalar *values;
int ndim=(getdata.natom/2)*((getdata.natom/2)-1)/2;
int ndimdmat2=(getdata.natom)*(getdata.natom)*(getdata.natom)*(getdata.natom);
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;
int nstates = getdata.ntrou*3;
double projvec[getdata.nroots*nstates], weightproj=0.0;
//double densmat2[getdata.natom][getdata.natom][getdata.natom][getdata.natom];
//memset(densmat2, 0, sizeof(densmat2));
////double *densmat2;
////densmat2 = (double *)malloc(sizeof(double)*ndimdmat2);
//memset(densmat2, 0, sizeof(double)*ndimdmat2);
if(mpiid==0)printf("Initializing Slepc vars\n");
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr);
@ -250,13 +256,17 @@ int main(int argc,char **argv)
get_s2(xr, &Istart, &Iend, values, &getdata.natom, &norm, &norm2, &norm3, &norm4, &xymat, &xymat2, &xymat3, &xymat4, &weight3,
&getdata.s21a1, &getdata.s21a2, &getdata.s21b1, &getdata.s21b2, &getdata.s22a1, &getdata.s22a2,
&getdata.s22b1, &getdata.s22b2, &getdata.s23a1, &getdata.s23a2,
&getdata.s23b1, &getdata.s23b2, &getdata.postrou, natomax);
&getdata.s23b1, &getdata.s23b2, &getdata.postrou1, &getdata.postrou2, &getdata.postrou3, natomax);
// get_s2_cyclic(xr, &Istart, &Iend, values, &getdata.natom, &norm, &norm2, &norm3, &norm4, &xymat, &xymat2, &xymat3, &xymat4,
// &getdata.s21a1, &getdata.s21a2, &getdata.s21b1, &getdata.s21b2, &getdata.s22a1, &getdata.s22a2,
// &getdata.s22b1, &getdata.s22b2, &getdata.s23a1, &getdata.s23a2,
// &getdata.s23b1, &getdata.s23b2, &getdata.postrou, natomax);
// get_1rdm(values, &Istart, &Iend, &getdata.natom, &trace1rdm, natomax);
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);
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);
VecRestoreArray(vec2,&values);
ierr = VecRestoreArray(xr, &valxr);CHKERRQ(ierr);
@ -269,7 +279,7 @@ int main(int argc,char **argv)
MPI_Reduce(&norm2, &normfin2, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
MPI_Reduce(&norm3, &normfin3, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
MPI_Reduce(&norm4, &normfin4, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
// MPI_Reduce(&trace1rdm, &trace1rdmfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
MPI_Reduce(&trace1rdm, &trace1rdmfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
// printf("done calc densmat\n");
// for(ll=0;ll<getdata.natom/2;ll++){
// for(kk=0;kk<getdata.natom/2;kk++){
@ -309,11 +319,11 @@ int main(int argc,char **argv)
if(!mpiid){
XS=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin/normfin)));
// XS2=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin2/normfin2)));
// XS3=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin3/normfin3)));
XS2=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin2)));
XS3=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin3)));
XS4=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin4/normfin4)));
XS2=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin2/normfin2)));
XS3=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin3/normfin3)));
// XS2=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin2)));
// XS3=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin3)));
// XS4=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin4)));
XS4=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin4/normfin4)));
W3=weight3fin/normfin2;
// W3=weight3fin;
@ -337,7 +347,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 %18f %18f %18f\n",(double)re,(double)error,(double)XS,(double)XS2,(double)XS3, (double)W3);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," %18f %12g %18f %18f %18f %18f %18f %18f %18f %18f %18f %18f %18f\n",(double)re,(double)error,(double)XS,(double)XS2,(double)XS3, (double)XS4, (double)weightproj, (double)projvec[i*6 + 0], (double)projvec[i*6 + 1], (double)projvec[i*6 + 2], (double)projvec[i*6 + 3], (double)projvec[i*6 + 4], (double)projvec[i*6 + 5]);CHKERRQ(ierr);
}
VecScatterDestroy(&scatter);
VecDestroy(&vec2);
@ -345,6 +355,7 @@ int main(int argc,char **argv)
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
}
//free(densmat2);
ierr = EPSDestroy(&eps);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&xr);CHKERRQ(ierr);

View File

@ -25,16 +25,19 @@
*/
void get_1rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace1rdm, const int natomax){
long int ideter[natomax];
long int ideter2[natomax];
int ideter[natomax];
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;
//int ndim=(*natom)*(*natom)/8-(*natom)/2;
int ndim=(*natom)/2;
double densmat[ndim][ndim];
memset(densmat, 0, sizeof(densmat[0][0]) * ndim * ndim);
memset(ideter , 0, sizeof(ideter[0] ) * natomax);
memset(ideter2, 0, sizeof(ideter2[0]) * natomax);
for(kko=0;kko<(*natom/2);kko++){
for(kok=0;kok<(*natom/2);kok++){
@ -54,17 +57,23 @@ void get_1rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom,
densmat[kko][kok]=densmat[kko][kok]+valxr[iiii]*valxr[iaa2];
}
}
if(kko == kok && ideter[kko] != 3){
if(kko == kok && ideter[kko] != 3l){
densmat[kko][kko]=densmat[kko][kko]+valxr[iiii]*valxr[iiii];
}
}
if(kko == kok){
*trace1rdm+=densmat[kko][kko];
}
}
}
//printf("\nDensmat 1body\n");
//for(kko=0;kko<(*natom/2);kko++){
// printf(" %4.4f ", densmat[kko][kko]);
//}
//printf("\n");
} /** END **/
/*
@ -81,10 +90,11 @@ void get_1rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom,
* =====
* trace = trace
*/
void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, double densmat2[*natom][*natom][*natom][*natom], 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){
void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, const int natomax){
long int ideter[natomax];
long int ideter2[natomax];
int ideter[natomax];
int ideter2[natomax];
int kko,kok,kkio;
int mmo,mom,mmio;
long int ii;
@ -92,6 +102,12 @@ void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom,
long int iii;
long int iaa2, iaa;
long int nrow=-1, ncol=-1;
int ndimdmat2=(*natom)*(*natom)*(*natom)*(*natom)/16;
double dmat2=0.0;
double densmat2[*natom/2][*natom/2][*natom/2][*natom/2];
//double *densmat2;
//densmat2 = (double *)malloc(sizeof(double)*ndimdmat2);
memset(densmat2, 0, sizeof(double)*ndimdmat2);
//int ndim=(*natom/2)*((*natom/2)-1)/2;
//double densmat2[ndim][ndim];
//memset(densmat2, 0, sizeof(densmat2[0][0]) * ndim * ndim);
@ -119,21 +135,28 @@ void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom,
ideter2[kok]=ideter[mom];
ideter2[mom]=3;
adr_(ideter2, &iaa2);
densmat2[kko][kok][mmo][mom]=densmat2[kko][kok][mmo][mom]+valxr[iiii]*valxr[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];
//densmat2[kko][kok][mmo][mom]=densmat2[kko][kok][mmo][mom]+valxr[iiii]*valxr[iiii];
}
if(kko == mmo && kok == mom && ideter[kko]==3 && ideter[kok]==3){
dmat2 += 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];
if(kko == mmo && kok == mom){
printf("%d\t%d\t%d\t%d\t%18f\n",kko,mmo,kok,mom,dmat2);
//printf("%d\t%d\t%d\t%d\t%18f\n",kko,mmo,kok,mom,densmat2[kko][kok][mmo][mom]);
}
}
}
dmat2=0.0;
}

View File

@ -6,4 +6,5 @@
#include <string.h>
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);
//void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, double densmat2[*natom][*natom][*natom][*natom], const int natomax);
void get_2rdm(PetscScalar *valxr, PetscInt *Istart, PetscInt *Iend, int *natom, PetscReal *trace2rdm, const int natomax);

View File

@ -26,7 +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){
int *s22b1, int *s22b2, int *s23a1, int *s23a2, int *s23b1, int *s23b2, int *postrou1, int *postrou2, int *postrou3, const int natomax){
long int iaa2, iaa;
long int iii;
int ideter[natomax];
@ -91,9 +91,9 @@ void get_s2(Vec xr, PetscInt *Istart, PetscInt *Iend, PetscScalar *valxr, int *n
pos3=kko;
}
}
if(ntrouboit1==1 && pos1 == *postrou)okboit1=1;
if(ntrouboit2==1 && pos2 == *postrou)okboit2=1;
if(ntrouboit3==1 && pos3 == *postrou)okboit3=1;
if(ntrouboit1==1 && pos1 == *postrou1)okboit1=1;
if(ntrouboit2==1 && pos2 == *postrou2)okboit2=1;
if(ntrouboit3==1 && pos3 == *postrou3)okboit3=1;
if(okboit1){
*norm2=*norm2+valxr[iiii]*valxr[iiii];
}
@ -681,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 | %d %d %d %d %d\n", *norm, *norm3, *xymat, *xymat3, *s22a1, *s22a2, *s22b1, *s22b2, *postrou);
//printf(" norm = %18f %18f %18f %18f xymat = %18f %18f %18f %18f | %d %d %d %d %d\n", *norm, *norm2, *norm3, *norm4, *xymat, *xymat2, *xymat3, *xymat4, *s21a1, *s21a2, *s21b1, *s21b2, *postrou);
//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr);
}

View File

@ -4,8 +4,11 @@
#include <ctype.h>
#include <string.h>
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);
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 *postrou1, int *postrou2, int *postrou3, const int natomax);
void get_s2_mov(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *,
int *,

View File

@ -34,7 +34,7 @@ void Data_new(FILE* file, Data* dat) {
/* note that fgets doesn't strip the terminating \n, checking its
presence would allow to handle lines longer than sizeof(line) */
if (count != 30){
if (count != 32){
count++;
switch(count){
case 1:
@ -266,15 +266,21 @@ void Data_new(FILE* file, Data* dat) {
dat->s23b2=atol(line);
break;
case 27:
dat->postrou=atol(line);
dat->postrou1=atol(line);
break;
case 28:
dat->fix_trou1=atol(line);
dat->postrou2=atol(line);
break;
case 29:
dat->fix_trou2=atol(line);
dat->postrou3=atol(line);
break;
case 30:
dat->fix_trou1=atol(line);
break;
case 31:
dat->fix_trou2=atol(line);
break;
case 32:
dat->print_wf = atol(line);
break;
default:

View File

@ -34,7 +34,9 @@ typedef struct {
int s23a2;
int s23b1;
int s23b2;
int postrou;
int postrou1;
int postrou2;
int postrou3;
long int fix_trou1;
long int fix_trou2;
long int print_wf;