now calculates ntot automatically

Makefile

@@ -27,6 +27,9 @@ directories: ${OBJ_DIR} ${LIB_DIR} ${BIN_DIR}
 ${LIB_DIR}/irpf90.a: directories
 	cd ${SRC_DIR} && irpf90 init && $(MAKE) irpf90.a && cp irpf90.a ../${LIB_DIR}
 
+${OBJ_DIR}/get_ntot.o: ${SRC_DIR}/get_ntot.c directories chkopts
+	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
+
 ${OBJ_DIR}/read2.o: ${SRC_DIR}/read2.c directories chkopts
 	${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o  $@ $< ${SLEPC_EPS_LIB} 
 
@@ -48,6 +51,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}/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}/stimsyr.h chkopts
-	    -${CLINKER} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -o ${BIN_DIR}/ex1 ${OBJ_DIR}/ex1.o ${OBJ_DIR}/read2.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}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3
+${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}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3
 #    ${RM} ex1.o read2.o

src/ex1.c

@@ -4,6 +4,7 @@
 #include "read2.h"
 #include "stimsyr.h"
 #include "get_s2.h"
+#include "get_ntot.h"
 
 #undef __FUNCT__
 #define __FUNCT__ "main"
@@ -44,6 +45,8 @@ int main(int argc,char **argv)
   PetscInt		 nlocal;
   
   Data_new(file, &getdata);
+  getdata.n = get_ntot(getdata.FAM1, getdata.natom, getdata.isz, getdata.ntrou, getdata.fix_trou1, getdata.fix_trou2);
+
   nlocal = getdata.n/getdata.npar;
 //printf("n=%ld\t nsites=%d\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.natom, getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz);
 
@@ -89,9 +92,6 @@ int main(int argc,char **argv)
   double            nelfin, s2densfin;
   double            densmat2[getdata.natom][getdata.natom][getdata.natom][getdata.natom];
   memset(densmat2, 0, sizeof(densmat2));
-  // fn = x^2+x^3
-//PetscInt          range[] ={165094,310638,438886,551954,651820,740324,819168,889916,953994,1012690,1067154,1118398,1167296,1214584,1260860,1306584,1352078,1517172,1662716,1790964,1904032,2003898,2092402,2171246,2241994,2306072,2364768,2419232,2470476,2519374,2566662,2612938,2658662,2704156,2869250,3014794,3143042,3256110,3355976,3444480,3523324,3594072,3658150,3716846,3771310,3822554,3871452,3918740,3965016,4010740,4056234,4221328,4366872,4495120,4608188,4708054,4796558,4875402,4946150,5010228,5068924,5123388,5174632,5223530,5270818,5317094,5362818,5408312,5573406,5718950,5847198,5960266,6060132,6148636,6227480,6298228,6362306,6421002,6475466,6526710,6575608,6622896,6669172,6714896,6760390,6925484,7071028,7199276,7312344,7412210,7500714,7579558,7650306,7714384,7773080,7827544,7878788,7927686,7974974,8021250,8066974,8112468,8277562,8423106,8551354,8664422,8764288,8852792,8931636,9002384,9066462,9125158,9179622,9230866,9279764,9327052,9373328,9419052,9464546,9629640,9775184,9903432,10016500,10116366,10204870,10283714,10354462,10418540,10477236,10531700,10582944,10631842,10679130,10725406,10771130,10816624,10981718,11127262,11255510,11368578,11468444,11556948,11635792,11706540,11770618,11829314,11883778,11935022,11983920,12031208,12077484,12123208,12168702,12333796,12479340,12607588,12720656,12820522,12909026,12987870,13058618,13122696,13181392,13235856,13287100,13335998,13383286,13429562,13475286,13520780,13685874,13831418,13959666,14072734,14172600,14261104,14339948,14410696,14474774,14533470,14587934,14639178,14688076,14735364,14781640,14827364,14872858,15037952,15183496,15311744,15424812,15524678,15613182,15692026,15762774,15826852,15885548,15940012,15991256,16040154,16087442,16133718,16179442,16224936};
-
 
   SlepcInitialize(&argc,&argv,(char*)0,NULL);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr);

src/get_ntot.c

@@ -0,0 +1,39 @@
+#include "get_ntot.h"
+
+int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2){
+  int            tnt1, tnt2;
+  int            natom2;
+  if(FAM1){
+      if(fix_trou1 == fix_trou2){
+          natom2 = natom/2;
+      }
+      else{
+          natom2 = fix_trou2 - fix_trou1;
+      }
+  }
+  else{
+      natom2 = natom;
+  }
+
+  tnt1 =         (int)exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1))));
+  int            nalpha, nbeta;
+
+    if((((natom-ntrou) + 2*isz) % 2) == 0){
+        nalpha=(natom-ntrou+2*isz)/2;
+        nbeta=(natom -ntrou-2*isz)/2;
+        if(((natom-ntrou)/2) == isz){
+            nbeta=0;
+        }
+    }
+    else{
+        nalpha=(natom-ntrou+2*isz+1)/2;
+        nbeta=(natom -ntrou-2*isz-1)/2;
+        if(((natom-ntrou+1)/2) == isz){
+            nbeta=0;
+        }
+    }
+
+  tnt2 =         (int) exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1))));
+//printf("nalpha=%d nbeta=%d | | %d %d ntot=%d\n",nalpha, nbeta, tnt1, tnt2, tnt1*tnt2);
+  return tnt1*tnt2;
+}

src/get_ntot.h

@@ -0,0 +1,3 @@
+#include <tgmath.h>
+
+int get_ntot(_Bool Fam1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2);