mirror of https://github.com/LCPQ/DEHam
144 lines
4.4 KiB
C
144 lines
4.4 KiB
C
#include <stdio.h>
|
|
#include <petsctime.h>
|
|
#include <slepceps.h>
|
|
#include <stdlib.h>
|
|
#include <ctype.h>
|
|
#include <string.h>
|
|
#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){
|
|
|
|
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){
|
|
|
|
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 **/
|
|
|