LCPQ/DEHam
2
0
Fork 0
mirror of https://github.com/LCPQ/DEHam synced 2024-12-22 04:13:44 +01:00
Code Issues Releases Wiki Activity

Merge pull request #2 from v1j4y/master

Browse Source 
Merge Vijay
This commit is contained in:
Anthony Scemama 2020-03-05 15:42:32 +01:00 committed by GitHub
commit 27f031d20c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
42 changed files with 4084 additions and 299 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.gitignore vendored Normal file
View File

@ -0,0 +1,4 @@
*bin/*
*/obj/*
*mod.mod*
irpf90.a

57
Makefile
View File

@ -1,23 +1,56 @@
include ${SLEPC_DIR}/lib/slepc/conf/slepc_common
CC=gcc
#CC=gcc
#FC = ifort
CLINKER = mpicc -fPIC #-wd1572 -O3 -axAVX,SSE4.2 -fno-alias -no-prec-div -no-prec-sqrt -ip
MAKE = /usr/bin/make
MKDIR_P = mkdir -p
OBJ_DIR = obj
MKDIR_P = /bin/mkdir -p
OBJ_DIR := obj
LIB_DIR := libs
BIN_DIR := bin
SRC_DIR := src
irpf90.a:
cd src && irpf90 init && $(MAKE) && cp irpf90.a ../
.PHONY: ex1
ex1: ${BIN_DIR}/ex1
${OBJ_DIR}:
${MKDIR_P} ${OBJ_DIR}
directories: ${OBJ_DIR}
${LIB_DIR}:
${MKDIR_P} ${LIB_DIR}
obj/read2.o: src/read2.c directories
${CC} -c -o $@ $<
${BIN_DIR}:
${MKDIR_P} ${BIN_DIR}
obj/ex1.o: src/ex1.c
-${CC} -c -o $@ $< ${SLEPC_EPS_LIB}
directories: ${OBJ_DIR} ${LIB_DIR} ${BIN_DIR}
ex1: obj/read2.o irpf90.a obj/ex1.o src/read2.h src/stimsyr.h chkopts
-${CLINKER} -o ex1 obj/ex1.o obj/read2.o irpf90.a ${SLEPC_EPS_LIB}# -lifcore -lirc -lcomposerxe_gen_helpers_core_2.3
${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}
${OBJ_DIR}/get_s2.o: ${SRC_DIR}/get_s2.c directories chkopts
${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB}
${OBJ_DIR}/get_s2_cyclic.o: ${SRC_DIR}/get_s2_cyclic.c directories chkopts
${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB}
${OBJ_DIR}/get_s2_mov.o: ${SRC_DIR}/get_s2_mov.c directories chkopts
${CC} ${SLEPC_INCLUDE} ${PETSC_CC_INCLUDES} -c -o $@ $< ${SLEPC_EPS_LIB}
${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_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}
# ${RM} ex1.o read2.o

62
README.md
View File

@ -7,57 +7,63 @@ Double Exchange Hamiltonian: Complete Version
(under GNU GENERAL PUBLIC LICENSE v2)
1. Dependencies
_Dependencies_
---------------
1. [PETSc](https://www.mcs.anl.gov/petsc/documentation/installation.html) and [SLEPc](http://slepc.upv.es/documentation/current/docs/instal.htm)
1. [PETSc](https://www.mcs.anl.gov/petsc/documentation/installation.html) and [SLEPc](http://slepc.upv.es/documentation/instal.htm)
2. [IRPF90](https://github.com/scemama/irpf90)
2. Compiling
_Compiling_
------------
1. Export environment variables for PETSc and SLEPc
```shell
export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION}
export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION}
```
```shell
export PETSC_DIR=${PATH_TO_PETSC_INSTALLATION}
export SLEPC_DIR=${PATH_TO_SLEPC_INSTALLATION}
export C_INCLUDE_PATH+=:$PETSC_DIR/include/:$SLEPC_DIR/include:$PETSC_DIR/arch-linux2-c-debug/include/:$SLEPC_DIR/arch-linux2-c-debug/include
# The "arch-linux2-c-debug" directory can have different names depending on PETSC and SLEPC installation procedure.
```
2. Make the executable
```shell
make ex1
```
```shell
make ex1
```
3. Using DEHam
_Using DEHam_
---------------
1. The DEHam program requires an input file which
has the topology of the Hamiltonian and the various parameters
as explained below in a sample inputfile:
```python
140 # The total number of determinants
7 # The largest number of non-zero elements per row
2 # The number of processors used in parallel
1 # The number of holes
0 # The isz (ms-1/2) value
1,2,3,1,2,3,4,5,6,7 # The topology of the system is specified here
2,3,4,8,7,6,5,6,7,8 # first and second line contain the two sites linked
1,1,1,2,2,2,2,3,3,3 # third line contains the type of link (1 for t, J 2 for K and 3 for none)
.1430,-0.20,0.0000 # The three types of links this line gives J, K
.1430,-0.20,0.0000 #
-1.00,0.0,0.00 # This line gives t
```
```python
8 # The number of orbitals (total)
140 # The largest number of non-zero elements per row (Multiple of Ndet)
1 # The total number of processors used in parallel (Multiple of Ndet)
1 # The number of holes
0 # The isz (ms-1/2) value
true # Restrict the hole to the 1'st (i.e. half of natom) Family of states. *false* for no restrictions
1,2,3,1,2,3,4,5,6,7 # The topology of the system is specified here
2,3,4,8,7,6,5,6,7,8 # first and second line contain the two sites linked
1,1,1,2,2,2,2,3,3,3 # third line contains the type of link (1 for t or J, 2 for K and 3 for none)
.1430,-0.20,0.0000 # The three types of links this line gives J, K
.1430,-0.20,0.0000 #
-1.00,0.0,0.00 # This line gives t
1 # Currently unused (Perhaps can be used for potential energy per site in the future.)
1 # The total number of roots
```
2. running DEHam
```shell
mpiexec -n [nprocs] ./ex1 inpfile
```
```shell
mpiexec -n [nprocs] ./ex1 inpfile
```
4. Publications using this code
_Publications using this code_
-------------------------------
1. High-Spin Chains and Crowns from Double-Exchange Mechanism [doi:10.3390/cryst6040039](http://www.dx.doi.org/10.3390/cryst6040039)

14
examples/small.inp Normal file
View File

@ -0,0 +1,14 @@
8
140
1
1
0
true
1,2,3,1,2,3,4,5,6,7
2,3,4,8,7,6,5,6,7,8
1,1,1,2,2,2,2,3,3,3
.1430,-0.20,0.0000
.1430,-0.20,0.0000
-1.00,0.0,0.00
1
2

14
examples/three.inp Normal file
View File

@ -0,0 +1,14 @@
3
1
1
1
1
true
1,2,1,2,3,4,5
2,3,6,5,4,5,6
1,1,2,2,2,3,3
.1430,-0.20,0.0000
.1430,-0.20,0.0000
-1.00,0.0,0.00
1
1

5
src/.gitignore vendored Normal file
View File

@ -0,0 +1,5 @@
IRPF90_temp/
IRPF90_man/
irpf90.make
irpf90_entities
tags

19
src/adr.irp.f
View File

@ -1,4 +1,5 @@
subroutine adr(ideter,add)
use iso_c_binding
implicit none
BEGIN_DOC
! this subroutine provides the address of a detrminant
@ -7,25 +8,25 @@ subroutine adr(ideter,add)
! matches the given determinant.
END_DOC
integer,INTENT(INOUT)::ideter(natomax)
integer(kind=selected_int_kind(16)),INTENT(INOUT)::add
integer(kind=selected_int_kind(16))::det,deth,addh,detnew
integer(C_SIZE_T),INTENT(INOUT)::add
integer(C_SIZE_T)::deti,dethi,addh,detnew
integer::count,i,j
det=0
deti=0
detnew=0
deth=0
dethi=0
count=0
call conv(ideter,det,deth)
call conv(ideter,deti,dethi)
Do i=0,natom-1
if(BTEST(deth,i))then
if(BTEST(dethi,i))then
count=count+1
endif
if(BTEST(det,i))then
if(BTEST(deti,i))then
detnew=IBSET(detnew,i-count)
endif
enddo
det=detnew
call searchdet(det,add,deth,addh)
deti=detnew
call searchdet(deti,add,dethi,addh)
add = add + (nt1-addh)*(nt2)

0
src/adrfull.irp.f → src/adrfull_old
View File

58
src/adrnew.irp.f Normal file
View File

@ -0,0 +1,58 @@
subroutine adrfull()
implicit none
use iso_c_binding
BEGIN_DOC
! this subroutine provides the address of a detrminant
! given in old format.
! It searches in a list of generated determinants and
! matches the given determinant.
END_DOC
integer,dimension(natomax)::ideter
integer(C_SIZE_T)::add
integer(C_SIZE_T)::deti,dethi,addh,detnew
integer::count,i,j
deti=0
detnew=0
dethi=0
do j=1,detfound
detnew=0
count=0
ideter=foundet(:,j)
call conv(ideter,deti,dethi)
Do i=0,natom-1
if(BTEST(dethi,i))then
count=count+1
endif
if(BTEST(deti,i))then
detnew=IBSET(detnew,i-count)
endif
enddo
deti=detnew
foundadd(j,1)=deti
foundadd(j,3)=j
foundaddh(j,1)=dethi
foundaddh(j,3)=j
call searchdet(deti,add,dethi,addh)
! enddo
! call sort()
! call searchdetfull()
! call desort()
! do i=1,detfound
! add = foundadd(i,2)
! addh = foundaddh(i,2)
foundadd(j,2) = add
foundaddh(j,2)= addh
add = add + (nt1-addh)*(nt2)
foundetadr(j)=add
enddo
10 FORMAT(B64,I8,F8.2)
15 FORMAT(B64,I8,I8,I8)
11 FORMAT(B64,I3,B64)
12 FORMAT(I5,$)
13 FORMAT(B64,B64)
14 FORMAT(B64,I14)
16 FORMAT(B64,I14,I14)
end

833
src/analyse.irp.f Normal file
View File

@ -0,0 +1,833 @@
SUBROUTINE ANALYSE(vect, dimvect, startvect, endvect, xymat2, norm2)
! INCLUDE "nbtots.prm"
use iso_c_binding
IMPLICIT NONE
INTEGER dimvect, nbtots, startvect, endvect
REAL*8,dimension(dimvect)::vect
INTEGER (C_SIZE_T)::add,kvect
INTEGER (C_SIZE_T)::iaa2,i
INTEGER ,dimension(natomax)::ideter
INTEGER ,dimension(natomax)::ideter2
REAL*8,allocatable ::xz(:)
REAL*8::xmat,xymat
REAL*8::xmat1,xymat1
REAL*8::xmat2
REAL*8,INTENT(INOUT)::xymat2
REAL*8::xmat3,xymat3
REAL*8::sym,nonsym,proj_trou
REAL*8,allocatable ::xalpha1(:)
! REAL*8,allocatable ::vect(:)
REAL*8,INTENT(INOUT)::norm2
REAL*8::norm,norm1,norm3,proj_2trou
REAL*8::t1,t2,XS,XS1,XS2,XS3
REAL*8::resta_mono,resta_one,resta_bi,delta
INTEGER ::kko,kok,kkio,j,eigen,nigen,count
INTEGER ::cntrou,countlvect,ndim,iaa
INTEGER ::ipt_1,ipt_2,iptemp_1,iptemp_2
INTEGER ::ipt_3,iptemp_3
INTEGER ::ibougetrou,jstart
REAL*8 ,allocatable::eigenvectors(:,:)
REAL*8 ,allocatable::eigenvalues(:)
REAL*8 ,allocatable::WORK(:)
REAL*8 ,allocatable::AP(:)
REAL*8 ,allocatable::densmat(:,:)
REAL*8 ,allocatable::densmat2(:,:)
REAL*8 proj_1,extradiag_dmat2,ionic,nonionic
REAL*8 proj_2,sum,conduction,prob,prob2
INTEGER INFO,nrow,ncol,mmo,mom,kk,k,omm,okk
CHARACTER*1 JOBZ,UPLO
INTEGER::RESTA=0
! allocate(vect(nbtots))
allocate(xalpha1(natomax))
allocate(xz (natom/2))
! OPEN (unit=59,file='FIL1',form='formatted',status='old')
! OPEN (unit=217,file='SBOX217',form='formatted',status='REPLACE')
! REWIND 59
! READ (59,*)
! print *,' in analyse', startvect, endvect
! print *,'nalpha=',nalpha,'nbeta=',nbeta
! PRINT *,natom,ntrou,nbtots,nt1,nt2,isz
! ndim=3
ndim=(natom/2)*((natom/2)-1)/2
allocate(AP((ndim)*((ndim)+1)/2))
allocate(WORK(3*(ndim)))
allocate(eigenvectors((ndim)*(ndim),1))
allocate(eigenvalues((ndim)))
allocate(densmat(ndim,ndim))
allocate(densmat2(ndim,ndim))
! Touch isz maxdet maxial maxlien maxplac nalpha natom natomax nbeta nbtots nt1 nt2 ntrou
! PRINT *,(vect(j),j=1,30)
IF(RESTA .eq. 1)THEN
do i=1,natom/2
if(mod(natom/2,2).eq.0)then
xz(i)=(((natom/2)/2)-0.5d0)-(i-1.0d0)
write(6 ,*)i,xz(i)
else
xz(i)=((natom/2)-1.0d0)/2.0d0-(i-1.0d0)
write(6 ,*)i,xz(i)
endif
enddo
ENDIF
!! PROVIDE det deth
nigen=1
DO eigen=1,nigen
! READ (59,10) (vect(j),j=1,nbtots)
IF (ntrou.eq.1) THEN
norm=0.d0
norm1=0.d0
norm2=0.d0
norm3=0.d0
count=0
cntrou=0
iaa=0
iaa2=0
countlvect=0
proj_2trou=0.d0
resta_bi=0.d0
resta_mono=0.d0
resta_one=0.d0
xymat = 0.0d0
xymat1 = 0.0d0
xymat2 = 0.d0
xymat3 = 0.d0
proj_1=0d0
proj_2=0d0
densmat=0d0
densmat2=0d0
extradiag_dmat2=0d0
conduction=0d0
nrow=0
ncol=0
jstart=1
ipt_1=0
ipt_2=0
ipt_3=0
iptemp_1=0
iptemp_2=0
iptemp_3=0
nonionic=0.d0
ionic=0.d0
prob=0.d0
prob2=0.d0
sym=0.d0
nonsym=0.d0
ibougetrou=0
DO kvect=1, endvect-startvect
CALL getdet(kvect+startvect,ideter)
!!----------------------------------------
!! RESTA
!!----------------------------------------
!!! mono
! proj_trou=vect(kvect)**2
! DO i=1,natom/2
! IF (ideter(i).eq.3) THEN
! delta=0.0d0
! ELSE
! delta=1.0d0
! END IF
! resta_mono=resta_mono+delta*xz(i)*xz(i)*proj_trou
! resta_one=resta_one+delta*xz(i)*proj_trou
! END DO
!!! bi
! DO i=1,natom/2
! DO j=1,natom/2
! IF (ideter(i).eq.3.or.ideter(j).eq.3.or.i.eq.j) &
! THEN
! delta=0.0d0
! ELSE
! delta=1.0d0
! END IF
! resta_bi=resta_bi+delta*xz(i)*xz(j)*proj_trou
! END DO
! END DO
!!----------------------------------------
!!----------------------------------------
!! Prob ionic non-ionic
!!----------------------------------------
! ipt_1=0
! ipt_2=0
! ipt_3=0
! DO kko=1,3
! IF(ideter(kko).eq.3)THEN
! ipt_1=ipt_1+1
! ENDIF
! ENDDO
! IF(ipt_1.eq.1)THEN
! DO kko=4,6
! IF(ideter(kko).eq.3)THEN
! ipt_2=ipt_2+1
! ENDIF
! ENDDO
! ENDIF
! IF(ipt_2.eq.1)THEN
! DO kko=7,9
! IF(ideter(kko).eq.3)THEN
! ipt_3=ipt_3+1
! ENDIF
! ENDDO
! ENDIF
! IF(ipt_3 .eq. 1)THEN
! nonionic=nonionic+vect(kvect)**2
! ELSE
! ionic=ionic+vect(kvect)**2
! ENDIF
!!----------------------------------------
!! S_box
!!----------------------------------------
xmat=0.0d0
xmat1=0.0d0
xmat2=0.0d0
! IF (.TRUE.)THEN
!! IF (ideter(6).eq.3 ) THEN
!! norm=norm+vect(kvect)**2
!! DO kko=5,7
!! DO kok=kko,7
!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
!! xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2)
!! ELSE
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=2
!! ideter2(kok)=1
!! CALL adr(ideter2, iaa2)
!! xmat=xmat+vect(kvect)*vect(iaa2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=1
!! ideter2(kok)=2
!! CALL adr(ideter2, iaa2)
!! xmat=xmat+vect(kvect)*vect(iaa2)
!! END IF
!! END IF
!! END DO
!! END DO
!!
!! DO kko=16,18
!! DO kok=kko,18
!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
!! xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2)
!! ELSE
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=2
!! ideter2(kok)=1
!! CALL adr(ideter2, iaa2)
!! xmat=xmat+vect(kvect)*vect(iaa2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=1
!! ideter2(kok)=2
!! CALL adr(ideter2, iaa2)
!! xmat=xmat+vect(kvect)*vect(iaa2)
!! END IF
!! END IF
!! END DO
!! END DO
!!
!! DO kko=5,7
!! DO kok=16,18
!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
!! xmat=xmat+(3.d0/4.d0)*(vect(kvect)**2)
!! ELSE
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
!! xmat=xmat+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=2
!! ideter2(kok)=1
!! CALL adr(ideter2, iaa2)
!! xmat=xmat+vect(kvect)*vect(iaa2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
!! xmat=xmat-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=1
!! ideter2(kok)=2
!! CALL adr(ideter2, iaa2)
!! xmat=xmat+vect(kvect)*vect(iaa2)
!! END IF
!! END IF
!! END DO
!! END DO
!! END IF
!!!----------------------------------------
!! xymat=xymat+xmat
!!
!! IF (ideter(7).eq.3 ) THEN
!! norm1=norm1+vect(kvect)**2
!! DO kko=5,9
!! DO kok=kko,9
!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
!! xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2)
!! ELSE
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=2
!! ideter2(kok)=1
!! CALL adr(ideter2, iaa2)
!! xmat1=xmat1+vect(kvect)*vect(iaa2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=1
!! ideter2(kok)=2
!! CALL adr(ideter2, iaa2)
!! xmat1=xmat1+vect(kvect)*vect(iaa2)
!! END IF
!! END IF
!! END DO
!! END DO
!!
!! DO kko=14,18
!! DO kok=kko,18
!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
!! xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2)
!! ELSE
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=2
!! ideter2(kok)=1
!! CALL adr(ideter2, iaa2)
!! xmat1=xmat1+vect(kvect)*vect(iaa2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=1
!! ideter2(kok)=2
!! CALL adr(ideter2, iaa2)
!! xmat1=xmat1+vect(kvect)*vect(iaa2)
!! END IF
!! END IF
!! END DO
!! END DO
!!
!! DO kko=5,9
!! DO kok=14,18
!! IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
!! xmat1=xmat1+(3.d0/4.d0)*(vect(kvect)**2)
!! ELSE
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
!! xmat1=xmat1+(1.d0/2.d0)*(vect(kvect)**2)
!! END IF
!! IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=2
!! ideter2(kok)=1
!! CALL adr(ideter2, iaa2)
!! xmat1=xmat1+vect(kvect)*vect(iaa2)
!! END IF
!! IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
!! xmat1=xmat1-(1.d0/2.d0)*(vect(kvect)**2)
!! DO kkio=1,natom
!! ideter2(kkio)=ideter(kkio)
!! END DO
!! ideter2(kko)=1
!! ideter2(kok)=2
!! CALL adr(ideter2, iaa2)
!! xmat1=xmat1+vect(kvect)*vect(iaa2)
!! END IF
!! END IF
!! END DO
!! END DO
!! END IF
!!!----------------------------------------
!! xymat1=xymat1+xmat1
IF (.TRUE.)THEN
norm2=norm2+vect(kvect)**2
DO kko=1,natom/2
DO kok=kko,natom/2
IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2)
ELSE
IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
END IF
IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
END IF
IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
DO kkio=1,natom
ideter2(kkio)=ideter(kkio)
END DO
ideter2(kko)=2
ideter2(kok)=1
CALL adr(ideter2, iaa2)
xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
END IF
IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
DO kkio=1,natom
ideter2(kkio)=ideter(kkio)
END DO
ideter2(kko)=1
ideter2(kok)=2
CALL adr(ideter2, iaa2)
xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
END IF
END IF
END DO
END DO
DO kko=(natom/2)+1,natom
DO kok=kko,natom
IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2)
ELSE
IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
END IF
IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
END IF
IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
DO kkio=1,natom
ideter2(kkio)=ideter(kkio)
END DO
ideter2(kko)=2
ideter2(kok)=1
CALL adr(ideter2, iaa2)
xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
END IF
IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
DO kkio=1,natom
ideter2(kkio)=ideter(kkio)
END DO
ideter2(kko)=1
ideter2(kok)=2
CALL adr(ideter2, iaa2)
xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
END IF
END IF
END DO
END DO
DO kko=1,natom/2
DO kok=(natom/2)+1,natom
IF (kok.eq.kko.and.ideter(kok).ne.3) THEN
xmat2=xmat2+(3.d0/4.d0)*(vect(kvect)**2)
ELSE
IF (ideter(kko).eq.1.and.ideter(kok).eq.1) THEN
xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
END IF
IF (ideter(kko).eq.2.and.ideter(kok).eq.2) THEN
xmat2=xmat2+(1.d0/2.d0)*(vect(kvect)**2)
END IF
IF (ideter(kko).eq.1.and.ideter(kok).eq.2) THEN
xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
DO kkio=1,natom
ideter2(kkio)=ideter(kkio)
END DO
ideter2(kko)=2
ideter2(kok)=1
CALL adr(ideter2, iaa2)
xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
END IF
IF (ideter(kko).eq.2.and.ideter(kok).eq.1) THEN
xmat2=xmat2-(1.d0/2.d0)*(vect(kvect)**2)
DO kkio=1,natom
ideter2(kkio)=ideter(kkio)
END DO
ideter2(kko)=1
ideter2(kok)=2
CALL adr(ideter2, iaa2)
xmat2=xmat2+vect(kvect)*vect(iaa2-startvect)
END IF
END IF
END DO
END DO
END IF
!----------------------------------------
! print *,"norm = ",norm2,"xmat2 = ",xmat2,"vect =",vect(kvect),"natom=",natom
xymat2=xymat2+xmat2
! XS=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat/norm)))
! XS1=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat1/norm1)))
XS2=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat2/norm2)))
END DO
! print *,eigen,xymat2,XS2
!----------------------------------------
! Resta and probabilities of dets
!----------------------------------------
! DO kvect=1,nbtots
! ideter2=ideter
! CALL getdet(kvect,ideter)
!
! IF (jstart.eq.1) THEN
! jstart=0
! ideter2=ideter
! DO kko=1,natom
! IF (ideter(kko).eq.3) THEN
! ipt_1=kko
! DO kok=kko+1,natom
! IF (ideter(kok).eq.3) THEN
! ipt_2=kok
! DO okk=kok+1,natom
! IF (ideter(okk).eq.3) THEN
! ipt_3=okk
! EXIT
! END IF
! END DO
! EXIT
! END IF
! END DO
! EXIT
! END IF
! END DO
! END IF
!
!
! DO kko=1,natom
! IF (ideter(kko).eq.3) THEN
! iptemp_1=kko
! DO kok=kko+1,natom
! IF (ideter(kok).eq.3) THEN
! iptemp_2=kok
! DO okk=kok+1,natom
! IF (ideter(okk).eq.3) THEN
! iptemp_3=okk
! EXIT
! END IF
! END DO
! EXIT
! END IF
! END DO
! EXIT
! END IF
! END DO
!
! IF (iptemp_1.ne.ipt_1.or.iptemp_2.ne.ipt_2 &
! .or. iptemp_3.ne.ipt_3) THEN
! ibougetrou=1
! ipt_1=iptemp_1
! ipt_2=iptemp_2
! ipt_3=iptemp_3
! ELSE
! proj_trou=proj_trou+vect(kvect)**2
! ibougetrou=0
! END IF
!
! IF (iptemp_1.eq.(9-iptemp_3+1) &
! .and. iptemp_2 .eq. 5 .and. iptemp_1.ne.4) THEN
! sym=sym+vect(kvect)**2
! ELSEIF (ideter(4).eq.3 .and. ideter(5).eq.3 &
! .and. ideter(6).eq.3) THEN
! nonsym=nonsym+vect(kvect)**2
! ELSE
!! ELSEIF(ideter(2).eq.3 .and. ideter(3).eq.3 &
!! .and. ideter(6).eq.3)then
! nonsym=nonsym+vect(kvect)**2
! END IF
!
! IF (ideter(2).eq.3 .and. ideter(4).eq.3 &
! .and. ideter(8).eq.3) THEN
! prob=prob+vect(kvect)**2
! END IF
!
!! IF (ideter(1).eq.3 .and. ideter(2).eq.3 &
!! .and. ideter(3).eq.3) THEN
!! prob2=prob2+vect(kvect)**2
!! END IF
!
! IF (ibougetrou.eq.1.or.kvect.eq.nbtots) THEN
!!----------------------------------------
!! mono
! DO i=1,natom/2
! IF (ideter2(i).eq.3) THEN
! delta=1.0d0
! ELSE
! delta=0.0d0
! END IF
! resta_mono=resta_mono+delta*xz(i)*xz(i)*proj_trou
! resta_one=resta_one+delta*xz(i)*proj_trou
! END DO
!! bi
! DO i=1,natom/2
! DO j=1,natom/2
! IF (ideter2(i).ne.3.or.ideter2(j).ne.3.or.i.eq.j) &
! THEN
! delta=0.0d0
! ELSE
! delta=1.0d0
! END IF
! resta_bi=resta_bi+delta*xz(i)*xz(j)*proj_trou
! END DO
! END DO
!!----------------------------------------
! proj_trou=0.d0
! proj_trou=proj_trou+vect(kvect)**2
! END IF
! END DO
!----------------------------------------
! One particle density matrix
!----------------------------------------
! DO kko=1,3
! DO kok=1,3
!
! DO kvect=1,nbtots
! CALL getdet(kvect,ideter)
! ideter2=ideter
! IF (ideter(kko).ne.3) THEN
! IF (ideter(kok).eq.3) THEN
! ideter2(kok)=ideter(kko)
! ideter2(kko)=3
! CALL adr(ideter2, iaa2)
! densmat(kko,kok)=densmat(kko,kok)+vect(kvect)* &
! vect(iaa2)
! END IF
! END IF
! IF (kko.eq.kok.and.ideter(kko).ne.3) THEN
! densmat(kko,kko)=densmat(kko,kko)+vect(kvect)**2
! END IF
! END DO
! END DO
! END DO
!----------------------------------------
! two particle density matrix
!----------------------------------------
! DO kko=1,(natom/2)-1
! DO kok=kko+1,natom/2
! nrow=nrow+1
! ncol=0
! DO mmo=1,(natom/2)-1
! DO mom=mmo+1,natom/2
! ncol=ncol+1
! DO kvect=1,nbtots
! CALL getdet(kvect,ideter)
! ideter2=ideter
! IF (ideter(kko).eq.3.and.ideter(kok) &
! .eq.3.and.ideter(mmo).ne.3.and.ideter(mom).ne.3) &
! THEN
! if(ideter(kok).ne.3 .and. ideter(mom).ne.3)then
! ideter2(kko)=ideter(mmo)
! ideter2(mmo)=3
! ideter2(kok)=ideter(mom)
! ideter2(mom)=3
! CALL adr(ideter2, iaa2)
! densmat2(nrow,ncol)=densmat2(nrow,ncol)+ &
! vect(kvect)*vect(iaa2)
! print *,nrow,ncol,kko,kok,mmo,mom,
! * densmat2(nrow,ncol)
! endif
! END IF
! IF (nrow.eq.ncol.and.ideter(mmo) &
! .ne.3.and.ideter(mom).ne.3) THEN
! densmat2(nrow,ncol)=densmat2(nrow,ncol)+ &
! vect(kvect)**2
! END IF
! END DO
! END DO
! END DO
! END DO
! END DO
!----------------------------------------
!----------------------------------------
! conduction
!----------------------------------------
! count=0
! DO kko=1,(natom/2)-2
! DO kok=kko+1,(natom/2)-1
! DO okk=kok+1,natom/2
!
! nrow=nrow+1
! ncol=0
! DO mmo=kko,kko+1
! DO mom=kok,kok+1
! DO omm=okk,okk+1
!
! ncol=ncol+1
! DO kvect=1,nbtots
! CALL getdet(kvect,ideter)
! ideter2=ideter
! IF (abs(kko-mmo).eq.1.or.abs(kok-mom).eq.1 &
! .or. abs(okk-omm).eq.1) THEN
! IF (mmo.le.natom/2.and.mom.le.natom/2 .and. &
! omm.le.natom/2) THEN
! IF (mmo.ne.mom .and. mom.ne.omm) THEN
! IF (ideter(kko).eq.3 .and. ideter(kok).eq.3 &
! .and. ideter(okk).eq.3) THEN
! ideter2(okk)=ideter2(omm)
! ideter2(omm)=3
! ideter2(kok)=ideter2(mom)
! ideter2(mom)=3
! ideter2(kko)=ideter2(mmo)
! ideter2(mmo)=3
! CALL adr(ideter2, iaa2)
!! count=0
!! do i=1,natom/2
!! if(ideter2(i).eq.3)then
!! count+=1
!! endif
!! enddo
!! print *,kko,kok,okk,mmo,mom,omm,iaa2
! conduction=conduction+dabs(vect(kvect)*vect(iaa2))
! END IF
! END IF
! END IF
! END IF
! END DO
!
! END DO
! END DO
! END DO
!
! END DO
! END DO
! END DO
!----------------------------------------
! DO j=1,ndim
! write(217,1022)j,(densmat(j,kko),kko=1,ndim)
! END DO
!----------------------------------------
! diagonalisation de mat
! affiche vecteur
! JOBZ='V'
! matrice sup
! UPLO='U'
! matrice en vecteur ligne ...
! extradiag_dmat2=0d0
! k=0
! DO j=1,ndim
! DO i=1,j-1
! if(i.ne.j)then
! extradiag_dmat2 = extradiag_dmat2 + dabs(densmat2(i,j))
! endif
! END DO
! END DO
! appel subroutine LAPACK de diagonalisation :: double précision !!
! INFO=0
! CALL DSPEV (JOBZ, UPLO, ndim, AP, eigenvalues, eigenvectors, &
! ndim, WORK, INFO)
! IF (INFO.ne.0) THEN
! PRINT *,'SUBROUTINE MATRIX: Error at dspev',info
! CALL exit (1)
! END IF
! proj_2=0.d0
! sum=0d0
! DO j=1,ndim
! proj_2=proj_2-eigenvalues(j)*log(eigenvalues(j))
! sum+=eigenvalues(j)
! write(214,*)eigenvalues(j)
! END DO
! XS=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat/norm)))
! XS2=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat2/norm2)))
! XS3=(1.d0/2.d0)*(-1.d0+dsqrt(1.d0+(4.d0*xymat3/norm3)))
! WRITE (217,*) eigen,XS,norm
END IF
END DO
10 FORMAT (E25.0)
1022 FORMAT(3x,I3,6(2x,F12.4))
END

15
src/conv.irp.f
View File

@ -1,20 +1,21 @@
subroutine conv(ideter,det,deth)
subroutine conv(ideter,deti,dethi)
use iso_c_binding
implicit none
BEGIN_DOC
! this routine converts a detrminant in the old
! format into the new one and returns the determinant.
END_DOC
integer,INTENT(INOUT)::ideter(natomax)
integer(kind=selected_int_kind(16)),INTENT(INOUT)::det
integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth
integer(C_SIZE_T),INTENT(INOUT)::deti
integer(C_SIZE_T),INTENT(INOUT)::dethi
integer::i
det=0
deth=0
deti=0
dethi=0
do i=1,natom
if(ideter(natom-i+1).eq.2 .and. ideter(natom-i+1).ne.3)then
det=IBSET(det,i-1)
deti=IBSET(deti,i-1)
elseif(ideter(natom-i+1).eq.3)then
deth=IBSET(deth,i-1)
dethi=IBSET(dethi,i-1)
endif
enddo
end

11
src/desort.irp.f
View File

@ -1,14 +1,15 @@
subroutine desort()
use iso_c_binding
implicit none
integer::i,j,ord,ordh
integer(kind=selected_int_kind(16))::add,addh,det,deth,addt
integer(C_SIZE_T)::add,addh,deti,dethi,addt
do i=1,detfound-1
do j=i+1,detfound
if(foundaddh(i,3).gt.foundaddh(j,3))then
deth = foundaddh(i,1)
dethi = foundaddh(i,1)
foundaddh(i,1) = foundaddh(j,1)
foundaddh(j,1) = deth
foundaddh(j,1) = dethi
addh = foundaddh(i,2)
foundaddh(i,2) = foundaddh(j,2)
foundaddh(j,2) = addh
@ -17,9 +18,9 @@ subroutine desort()
foundaddh(j,3) = ordh
endif
if(foundadd(i,3).gt.foundadd(j,3))then
det = foundadd(i,1)
deti = foundadd(i,1)
foundadd(i,1) = foundadd(j,1)
foundadd(j,1) = det
foundadd(j,1) = deti
add = foundadd(i,2)
foundadd(i,2) = foundadd(j,2)
foundadd(j,2) = add

4
src/elem_diag.irp.f
View File

@ -28,6 +28,10 @@
! if(yw)write(6,*)iaa,'diag,v1'
! endif
enddo
do i=1, natom
if(deter(i).ne.3) xmatd = xmatd + E(i)
enddo
xmatd = xmatd - E(natom+1)
!-----stockage de l element diag

295
src/ex1.c
View File

@ -1,56 +1,103 @@
#include <slepceps.h>
#include <petsctime.h>
#include "stimsyr.h"
#include <petscvec.h>
#include "read2.h"
#include "stimsyr.h"
#include "get_s2.h"
#include "get_ntot.h"
#undef __FUNCT__
#define __FUNCT__ "main"
void solvequad(double *a, double *b, double *c, double *res){
*res = -*b/(2.0*(*a)) + sqrt((*b)*(*b) - 4.0*(*a)*(*c))/(2.0*(*a));
}
int main(int argc,char **argv)
{
Mat A; /* problem matrix */
EPS eps; /* eigenproblem solver context */
EPSType type;
PetscReal error,tol,re,im;
PetscScalar kr,ki,value[700];
PetscReal norm=0.0;
PetscReal norm2=0.0;
PetscReal norm3=0.0;
PetscReal norm4=0.0;
PetscReal normfin=0.0;
PetscReal normfin2=0.0;
PetscReal normfin3=0.0;
PetscReal normfin4=0.0;
const int natomax=900;
PetscScalar kr,ki,value[natomax];
Vec xr,xi;
PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol;
PetscInt i,Istart,Iend,col[natomax],maxit,its,nconv,countcol;
PetscInt nev, ncv, mpd;
PetscLogDouble t1,t2,tt1,tt2;
//PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
PetscErrorCode ierr;
//PetscScalar eigr;
//PetscScalar eigi;
int mpiid;
char const* const fileName = argv[1];
FILE* file = fopen(fileName, "r");
Data getdata;
PetscInt nlocal;
/* gather the input data */
Data_new(file, &getdata);
//printf("n=%ld\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz);
getdata.n = get_ntot(getdata.FAM1, getdata.natom, getdata.isz, getdata.ntrou, getdata.fix_trou1, getdata.fix_trou2);
nlocal = getdata.n/getdata.npar;
//Vec Vr,Vi;
PetscScalar *valxr;
PetscInt indxr[nlocal];
char filename[PETSC_MAX_PATH_LEN]="FIL666";
PetscViewer viewer;
PetscBool ishermitian;
PetscInt kk,ll,iii2;
PetscInt kk,ll,mm,nn,iii2,iiii;
PetscInt ii;
long int iii;
long int tcountcol2,tcol[700],tcountcol[getdata.nnz];
double val[700];
long int tcountcol2,tcol[natomax],tcountcol[getdata.nnz];
double val[natomax];
PetscReal xymat=0.0;
PetscReal xymat2=0.0;
PetscReal xymat3=0.0;
PetscReal xymat4=0.0;
PetscReal xymatfin=0.0;
PetscReal xymatfin2=0.0;
PetscReal xymatfin3=0.0;
PetscReal xymatfin4=0.0;
PetscReal weight3fin = 0.0;
PetscReal XS = 0.0;
PetscReal XS2 = 0.0;
PetscReal XS3 = 0.0;
PetscReal XS4 = 0.0;
PetscReal W3 = 0.0;
PetscReal weight3 = 0.0;
PetscReal trace1rdm=0.0;
PetscReal trace1rdmfin=0.0;
PetscReal trace2rdm=0.0;
PetscReal trace2rdmfin=0.0;
IS from, to; /* index sets that define the scatter */
PetscInt idx_to[nlocal], idx_from[nlocal];
PetscScalar *values;
int ndim=(getdata.natom/2)*((getdata.natom/2)-1)/2;
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;
double densmat2[getdata.natom][getdata.natom][getdata.natom][getdata.natom];
memset(densmat2, 0, sizeof(densmat2));
SlepcInitialize(&argc,&argv,(char*)0,NULL);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,getdata.nnz*getdata.npar,NULL,getdata.nnz*getdata.npar,NULL,&A);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(A,getdata.nnz*getdata.npar,NULL,getdata.nnz*getdata.npar,NULL);CHKERRQ(ierr);
//ierr = MatSetFromOptions(A);CHKERRQ(ierr);
//ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,10*getdata.natom,NULL,10*getdata.natom,NULL,&A);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(A,10*getdata.natom,NULL,10*getdata.natom,NULL);CHKERRQ(ierr);
MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
ierr = PetscTime(&tt1);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," start: %d end: %d \n",Istart, Iend);CHKERRQ(ierr);
for (i=Istart; i<Iend; i+=getdata.nnz) {
tcountcol2=0;
@ -58,53 +105,46 @@ int main(int argc,char **argv)
tcountcol[kk]=0;
}
iii=i+1;
if(i%getdata.npar == 0 && mpiid==0){
ierr = PetscTime(&t1);CHKERRQ(ierr);
}
unit_l1_(
getdata.l1,
getdata.l2,
getdata.ktyp,
&iii,
&getdata.nnz,
&getdata.nnz,
getdata.xjjxy,
getdata.xjjz ,
getdata.xtt ,
getdata.E ,
tcountcol,
&getdata.ntrou,
&getdata.isz,
&getdata.fix_trou1,
&getdata.fix_trou2,
&getdata.FAM1,
tcol,
val);
if(i%getdata.npar == 0 && mpiid==0){
ierr = PetscTime(&t2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," i: %d\n mpiid: %d\ntime: %f\n",i,mpiid,t2-t1);CHKERRQ(ierr);
}
// if(i%getdata.npar == 0 && mpiid==0){
// ierr = PetscPrintf(PETSC_COMM_WORLD," i: %d \n",i);CHKERRQ(ierr);
// }
for(ll=0;ll<getdata.nnz;ll++){
for(kk=0;kk<tcountcol[ll]+1;kk++){
value[kk] = val[kk+tcountcol2];
col[kk] = tcol[kk+tcountcol2]-1;
// PetscPrintf(PETSC_COMM_WORLD,"value = %f col = %d\n",value[kk],col[kk]);
}
for(kk=tcountcol2+tcountcol[ll]+1;kk<700;kk++){
value[kk] = 0.0;
col[kk] = 0;
}
tcountcol2=tcountcol2 + tcountcol[ll]+1;
countcol=tcountcol[ll]+1;
if(i%getdata.npar == 0 && mpiid==0){
ierr = PetscTime(&t1);CHKERRQ(ierr);
}
iii2=i+ll;
ierr = MatSetValues(A,1,&iii2,countcol,col,value,INSERT_VALUES);CHKERRQ(ierr);
if(i%getdata.npar == 0 && mpiid==0){
ierr = PetscTime(&t2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," processor \ntime: %f\n",t2-t1);CHKERRQ(ierr);
}
for(kk=0;kk<tcountcol[ll]+1;kk++){
value[kk] = val[kk+tcountcol2];
col[kk] = tcol[kk+tcountcol2]-1;
}
for(kk=tcountcol2+tcountcol[ll]+1;kk<natomax;kk++){
value[kk] = 0.0;
col[kk] = 0;
}
tcountcol2=tcountcol2 + tcountcol[ll]+1;
countcol=tcountcol[ll]+1;
iii2=i+ll;
ierr = MatSetValues(A,1,&iii2,countcol,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = PetscTime(&tt2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to build the matrix: %f\n",tt2-tt1);CHKERRQ(ierr);
printf("time = %f mpiid = %d \n",tt2-tt1, mpiid);
ierr = PetscTime(&tt1);CHKERRQ(ierr);
@ -112,8 +152,6 @@ int main(int argc,char **argv)
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = PetscTime(&tt2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to assemble the matrix: %f\n",tt2-tt1);CHKERRQ(ierr);
//ierr = MatGetVecs(A,NULL,&xr);CHKERRQ(ierr);
//ierr = MatGetVecs(A,NULL,&xi);CHKERRQ(ierr);
ierr = MatCreateVecs(A,NULL,&xr);CHKERRQ(ierr);
ierr = MatCreateVecs(A,NULL,&xi);CHKERRQ(ierr);
@ -121,15 +159,14 @@ int main(int argc,char **argv)
ierr = EPSSetOperators(eps,A,NULL);CHKERRQ(ierr);
ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr);
ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr);
ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
tol = 1.e-8;
tol = 1.e-9;
maxit = 10000000;
ierr = EPSSetTolerances(eps,tol,maxit);CHKERRQ(ierr);
nev = 4;
ncv = 10;
mpd = 10;
nev = getdata.nroots;
ierr = EPSSetDimensions(eps,nev,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
ierr = PetscTime(&t1);CHKERRQ(ierr);
@ -141,29 +178,145 @@ int main(int argc,char **argv)
ierr = EPSGetType(eps,&type);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRQ(ierr);
ierr = EPSGetDimensions(eps,&nev,NULL,NULL);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenvalues: %D\n",nev);CHKERRQ(ierr);
ierr = EPSGetTolerances(eps,&tol,&maxit);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);CHKERRQ(ierr);
ierr = EPSGetConverged(eps,&nconv);CHKERRQ(ierr);
//ierr = EPSPrintSolution(eps,NULL);CHKERRQ(ierr);
if (nconv>0) {
/*
Display eigenvalues and relative errors
*/
ierr = PetscPrintf(PETSC_COMM_WORLD,
" k ||Ax-kx||/||kx||\n"
" ----------------- ------------------\n");CHKERRQ(ierr);
/*
Save eigenvectors, if == ested
*/
EPSGetConverged(eps,&nconv);
if (getdata.print_wf) {
PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU);
EPSIsHermitian(eps,&ishermitian);
for (i=0;i<nev;i++) {
EPSGetEigenvector(eps,i,xr,xi);
VecView(xr,viewer);
#if !defined(PETSC_USE_COMPLEX)
if (!ishermitian) { VecView(xi,viewer); }
#endif
}
PetscViewerDestroy(&viewer);
}
/*
* now analyzing the eigenvector
*/
if (nconv>0) {
ierr = PetscPrintf(PETSC_COMM_WORLD,
" k ||Ax-kx||/||kx|| <S^2>\n"
" ----------------- ----------------- ------------------\n");CHKERRQ(ierr);
for(i=0;i<nev;i++){
for (i=0;i<nconv;i++) {
/*
Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
ki (imaginary part)
*/
Vec vec2;
VecScatter scatter; /* scatter context */
ierr = EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);CHKERRQ(ierr);
xymat = 0.0;
xymat2 = 0.0;
xymat3 = 0.0;
xymat4 = 0.0;
weight3 = 0.0;
norm = 0.0;
norm2 = 0.0;
norm3 = 0.0;
norm4 = 0.0;
// ierr = PetscTime(&tt1);CHKERRQ(ierr);
// ierr = VecGetArray(xr, &valxr);CHKERRQ(ierr);
VecScatterCreateToAll(xr,&scatter,&vec2);
VecScatterBegin(scatter,xr,vec2,INSERT_VALUES,SCATTER_FORWARD);
VecScatterEnd(scatter,xr,vec2,INSERT_VALUES,SCATTER_FORWARD);
ierr = VecGetArray(vec2,&values);CHKERRQ(ierr);
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);
// 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_2rdm(values, &Istart, &Iend, &getdata.natom, &trace2rdm, densmat2, natomax);
// analyse_(valxr, (Iend-Istart), &Istart, &Iend, &xymat, &norm);
VecRestoreArray(vec2,&values);
ierr = VecRestoreArray(xr, &valxr);CHKERRQ(ierr);
MPI_Reduce(&xymat, &xymatfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
MPI_Reduce(&xymat2, &xymatfin2, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
MPI_Reduce(&xymat3, &xymatfin3, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
MPI_Reduce(&xymat4, &xymatfin4, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
MPI_Reduce(&weight3, &weight3fin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
MPI_Reduce(&norm, &normfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
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);
// printf("done calc densmat\n");
// for(ll=0;ll<getdata.natom/2;ll++){
// for(kk=0;kk<getdata.natom/2;kk++){
// gamma_p = gamma_p + 0.5*(densmat2[ll][kk][kk][ll] + densmat2[ll][kk][ll][kk]);
// gamma_m = gamma_m + 0.5*(densmat2[ll][kk][kk][ll] - densmat2[ll][kk][ll][kk]);
// }
// }
// MPI_Reduce(&trace2rdm, &trace2rdmfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
// MPI_Reduce(&gamma_p, &gamma_pfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
// MPI_Reduce(&gamma_m, &gamma_mfin, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
// if(mpiid==0){
// for(kk=0;kk<getdata.natom;kk++){
// for(ll=0;ll<getdata.natom;ll++){
// for(mm=0;mm<getdata.natom;mm++){
// for(nn=0;nn<getdata.natom;nn++){
//// printf("%d\t%d\t%d\t%d\t%18f\n",kk,ll,mm,nn,densmat2[kk][ll][mm][nn]);
// }
// }
// }
// }
// /* calc nel */
// a=1.0;
// b=-1.0;
// c=-2.0*(gamma_mfin + gamma_pfin);
// printf("\n gp= %18f gm= %18f a=%18f b=%18f c=%18f\n", gamma_pfin, gamma_mfin, a, b, c);
// nel = -b/(2.0*(a)) + sqrt((b)*(b) - 4.0*(a)*(c))/(2.0*(a));
//// solvequad(&a, &b, &c, &nel);
//
// /* calc s^2 */
// a=1.0;
// b=1.0;
// c=-1.0*((gamma_mfin - gamma_pfin) - nel*(nel - 4.0)/4.0);
// s2dens = -b/(2.0*(a)) + sqrt((b)*(b) - 4.0*(a)*(c))/(2.0*(a));
//// solvequad(&a, &b, &c, &s2dens);
// printf("\n mpiid = %d # trace = %18f nel = %18f s2dens = %18f\n", mpiid, trace2rdmfin, nel, s2dens);
// }
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)));
XS4=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin4/normfin4)));
W3=weight3fin/normfin2;
// W3=weight3fin;
}
xymatfin = 0.0;
normfin = 0.0;
/*
Compute the relative error associated to each eigenpair
* Compute the relative error associated to each eigenpair
*/
ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRQ(ierr);
@ -177,32 +330,14 @@ 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\n",(double)re,(double)error);CHKERRQ(ierr);
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);
}
VecScatterDestroy(&scatter);
VecDestroy(&vec2);
}
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
}
/*
Save eigenvectors, if requested
*/
//PetscOptionsGetString(NULL,NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
EPSGetConverged(eps,&nconv);
if (nconv>0) {
PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU);
EPSIsHermitian(eps,&ishermitian);
for (i=0;i<nconv;i++) {
EPSGetEigenvector(eps,i,xr,xi);
VecView(xr,viewer);
#if !defined(PETSC_USE_COMPLEX)
if (!ishermitian) { VecView(xi,viewer); }
#endif
}
PetscViewerDestroy(&viewer);
}
ierr = EPSDestroy(&eps);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&xr);CHKERRQ(ierr);

7
src/extra_diag.irp.f
View File

@ -1,8 +1,9 @@
subroutine extra_diag(tistart)
use iso_c_binding
implicit none
integer(kind=selected_int_kind(16)) :: iaa,iaa2,tistart,tistart2
integer(kind=selected_int_kind(16)) :: imat4,jmat4
integer(C_SIZE_T) :: iaa,iaa2,tistart,tistart2
integer(C_SIZE_T) :: imat4,jmat4
integer :: i,ik,iik,j
integer :: ik1,ik2,IC,k,ikmax,ikmin,count,count2,detfound2
integer,allocatable :: ideter2(:)
@ -99,7 +100,7 @@
tistart=tistart+1
enddo
Touch foundet foundetadr detfound foundadd foundaddh foundetdmat
Touch foundet foundetadr detfound foundadd foundaddh foundetdmat det deth
call adrfull()
do i=1,detfound

7
src/foundet.irp.f
View File

@ -1,8 +1,9 @@
use iso_c_binding
BEGIN_PROVIDER[integer, foundet,(natomax,maxlien)]
&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundetadr,(maxlien)]
&BEGIN_PROVIDER[integer(C_SIZE_T), foundetadr,(maxlien)]
&BEGIN_PROVIDER[real, foundetdmat,(maxlien)]
&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundadd,(maxlien,3)]
&BEGIN_PROVIDER[integer(kind=selected_int_kind(16)), foundaddh,(maxlien,3)]
&BEGIN_PROVIDER[integer(C_SIZE_T), foundadd,(maxlien,3)]
&BEGIN_PROVIDER[integer(C_SIZE_T), foundaddh,(maxlien,3)]
&BEGIN_PROVIDER[integer, detfound]
BEGIN_DOC
! provides all found determinants

143
src/get_dmat.c Normal file
View File

@ -0,0 +1,143 @@
#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 **/

9
src/get_dmat.h Normal file
View File

@ -0,0 +1,9 @@
#include <stdio.h>
#include <petsctime.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#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);

40
src/get_ntot.c Normal file
View File

@ -0,0 +1,40 @@
#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 = (lrint)(exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))));
//printf("%10.5f | tnt1=%d\n",exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))),tnt1);
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 = (lrint)(exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1)))));
//printf("natom2=%d fix_trou1=%d fix_trou2=%d nalpha=%d nbeta=%d | | %d %d ntot=%d\n",natom2, fix_trou1, fix_trou2, nalpha, nbeta, tnt1, tnt2, tnt1*tnt2);
return tnt1*tnt2;
}

4
src/get_ntot.h Normal file
View File

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

685
src/get_s2.c Normal file
View File

@ -0,0 +1,685 @@
#include <stdio.h>
#include <petsctime.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include "get_s2.h"
#include "get_val_iaa2.h"
/*
* This function simply calculates the S^2 value of the wavefunction
* Input
* =====
* Vr = The full vector
* Istart = Local starting id of the vector
* Iend = Local vector ending id
* valxr = Local vector values
* natom = number of orbitals
* Output
* ======
* norm = norm of the vector
* xymat = the S^2 value
*/
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){
long int iaa2, iaa;
long int iii;
int ideter[natomax];
int ideter2[natomax];
int kko,kok,kkio;
long int ii;
double xmat=0.0;
double xmat2=0.0;
double xmat3=0.0;
double xmat4=0.0;
double getvaliaa2;
PetscLogDouble t1,t2,tt1,tt2;
PetscErrorCode ierr;
PetscInt iiii;
int ntrouboit1=0;
int ntrouboit2=0;
int ntrouboit3=0;
int okboit1=0;
int okboit2=0;
int okboit3=0;
int mpiid;
int pos1=0;
int pos2=0;
int pos3=0;
MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);}
//ierr = PetscTime(&tt1);CHKERRQ(ierr);
for(ii=*Istart;ii<*Iend;ii++) {
iii = ii + 1;
// iiii = ii-*Istart;
iiii = ii;
xmat = 0.0;
xmat2 = 0.0;
xmat3 = 0.0;
xmat4 = 0.0;
ntrouboit1 = 0;
ntrouboit2 = 0;
ntrouboit3 = 0;
okboit1 = 0;
okboit2 = 0;
okboit3 = 0;
pos1 = 0;
pos2 = 0;
pos3 = 0;
getdet_(&iii, ideter);
*norm=*norm+valxr[iiii]*valxr[iiii];
for(kko=*s21a1;kko<=*s21a2;kko++){
if(ideter[kko]==3){
ntrouboit1++;
pos1=kko;
}
}
for(kko=*s22a1;kko<=*s22a2;kko++){
if(ideter[kko]==3){
ntrouboit2++;
pos2=kko;
}
}
for(kko=*s23a1;kko<=*s23a2;kko++){
if(ideter[kko]==3){
ntrouboit3++;
pos3=kko;
}
}
if(ntrouboit1==1 && pos1 == *postrou)okboit1=1;
if(ntrouboit2==1 && pos2 == *postrou)okboit2=1;
if(ntrouboit3==1 && pos3 == *postrou)okboit3=1;
if(okboit1){
*norm2=*norm2+valxr[iiii]*valxr[iiii];
}
if(okboit2){
*norm3=*norm3+valxr[iiii]*valxr[iiii];
}
if(okboit3){
*norm4=*norm4+valxr[iiii]*valxr[iiii];
}
/*
* calculate the weight of ms=5/2
*
* loop over the determinants to see if we have a S=5/2
*/
int countw = 0;
for(kko=*s21a1;kko<=*s21a2;kko++){
if(ideter[kko] == 2) countw=1;
}
for(kok=*s21b1;kok<=*s21b2;kok++){
if(ideter[kok] == 2) countw=1;
}
if(countw==0 && okboit1){
*weight3 += (valxr[iiii]*valxr[iiii]);
}
for(kko=0;kko<=(*natom/2)-1;kko++){
for(kok=kko;kok<=(*natom/2)-1;kok++){
if(kok == kko && ideter[kok] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko] == 1 && ideter[kok] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 1 && ideter[kok] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=2;
ideter2[kok]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=1;
ideter2[kok]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
for(kko=(*natom/2);kko<=*natom-1;kko++){
for(kok=kko;kok<=*natom-1;kok++){
if(kok == kko && ideter[kok] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko] == 1 && ideter[kok] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 1 && ideter[kok] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=2;
ideter2[kok]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=1;
ideter2[kok]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
for(kko=0;kko<=(*natom/2)-1;kko++){
for(kok=(*natom/2);kok<=*natom-1;kok++){
if(kok == kko && ideter[kok] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko] == 1 && ideter[kok] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 1 && ideter[kok] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=2;
ideter2[kok]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=1;
ideter2[kok]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
// if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);}
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
*xymat=*xymat+xmat;
*xymat2=*xymat2+xmat2;
*xymat3=*xymat3+xmat3;
*xymat4=*xymat4+xmat4;
// if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3);
}
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);
//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr);
}

38
src/get_s2.h Normal file
View File

@ -0,0 +1,38 @@
#include <stdio.h>
#include <slepceps.h>
#include <stdlib.h>
#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_mov(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *, const int natomax);
void get_s2_cyclic(Vec, PetscInt *, PetscInt *, PetscScalar *, int *, PetscReal *, PetscReal *,PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *,
int *, const int natomax);

783
src/get_s2_cyclic.c Normal file
View File

@ -0,0 +1,783 @@
#include <stdio.h>
#include <petsctime.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include "get_s2.h"
#include "get_val_iaa2.h"
/*
* This function simply calculates the S^2 value of the wavefunction
* Input
* =====
* Vr = The full vector
* Istart = Local starting id of the vector
* Iend = Local vector ending id
* valxr = Local vector values
* natom = number of orbitals
* Output
* ======
* norm = norm of the vector
* xymat = the S^2 value
*/
void get_s2_cyclic(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,
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){
long int iaa2, iaa;
long int iii;
int ideter[natomax];
int ideter2[natomax];
int kko,kok,kkio,kk;
int kko2,kok2;
long int ii;
double xmat=0.0;
double xmat2=0.0;
double xmat3=0.0;
double xmat4=0.0;
double getvaliaa2;
PetscLogDouble t1,t2,tt1,tt2;
PetscErrorCode ierr;
PetscInt iiii;
int ntrouboit1=0;
int ntrouboit2=0;
int ntrouboit3=0;
int okboit1=0;
int okboit2=0;
int okboit3=0;
int mpiid;
int pos1=0;
int pos2=0;
int pos3=0;
MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);}
//ierr = PetscTime(&tt1);CHKERRQ(ierr);
for(ii=*Istart;ii<*Iend;ii++) {
iii = ii + 1;
// iiii = ii-*Istart;
iiii = ii;
xmat = 0.0;
xmat2 = 0.0;
xmat3 = 0.0;
xmat4 = 0.0;
ntrouboit1 = 0;
ntrouboit2 = 0;
ntrouboit3 = 0;
okboit1 = 0;
okboit2 = 0;
okboit3 = 0;
pos1 = 0;
pos2 = 0;
pos3 = 0;
getdet_(&iii, ideter);
*norm=*norm+valxr[iiii]*valxr[iiii];
for(kko=0;kko<*natom/2;kko++){
if(ideter[kko]==3){
kk=kko;
}
}
*postrou = kk;
*s21a1 = kk-1;
// if(*s21a1<0){
// *s21a1 = (*natom/2) + (*s21a1)%(*natom/2);
// }
// else{
// *s21a1 = (*s21a1)%(*natom/2);
// }
*s22a1 = kk-1;
// if(*s22a1<0){
// *s22a1 = (*natom/2) + (*s22a1)%(*natom/2);
// }
// else{
// *s22a1 = (*s22a1)%(*natom/2);
// }
*s23a1 = kk-2;
// if(*s23a1<0){
// *s23a1 = (*natom/2) + (*s23a1)%(*natom/2);
// }
// else{
// *s23a1 = (*s23a1)%(*natom/2);
// }
//
*s21a2 = kk+1;
// if(*s21a2<0){
// *s21a2 = (*natom/2) + (*s21a2)%(*natom/2);
// }
// else{
// *s21a2 = (*s21a2)%(*natom/2);
// }
*s22a2 = kk+2;
// if(*s22a2<0){
// *s22a2 = (*natom/2) + (*s22a2)%(*natom/2);
// }
// else{
// *s22a2 = (*s22a2)%(*natom/2);
// }
*s23a2 = kk+2;
// if(*s23a2<0){
// *s23a2 = (*natom/2) + (*s23a2)%(*natom/2);
// }
// else{
// *s23a2 = (*s23a2)%(*natom/2);
// }
//
*s21b1 = *natom + *s21a1;
*s22b1 = *natom + *s22a1;
*s23b1 = *natom + *s23a1;
*s21b2 = *natom + *s21a2;
*s22b2 = *natom + *s22a2;
*s23b2 = *natom + *s23a2;
// if(mpiid==0)printf("postrou = %d\n",*postrou);
// if(mpiid==0)printf("1a1 = %d, 1a2 = %d, 1b1 = %d, 1b2 = %d\n",*s21a1,*s21a2,*s21b1,*s21b2);
// if(mpiid==0)printf("2a1 = %d, 2a2 = %d, 2b1 = %d, 2b2 = %d\n",*s22a1,*s22a2,*s22b1,*s22b2);
// if(mpiid==0)printf("3a1 = %d, 3a2 = %d, 3b1 = %d, 3b2 = %d\n",*s23a1,*s23a2,*s23b1,*s23b2);
// for(kko=*s21a1;kko<=*s21a2;kko++){
// if(ideter[kko]==3){
// ntrouboit1++;
// pos1=kko;
// }
// }
// for(kko=*s22a1;kko<=*s22a2;kko++){
// if(ideter[kko]==3){
// ntrouboit2++;
// pos2=kko;
// }
// }
// for(kko=*s23a1;kko<=*s23a2;kko++){
// if(ideter[kko]==3){
// ntrouboit3++;
// pos3=kko;
// }
// }
// if(ntrouboit1==1 && pos1 == *postrou)okboit1=1;
// if(ntrouboit2==1 && pos2 == *postrou)okboit2=1;
// if(ntrouboit3==1 && pos3 == *postrou)okboit3=1;
okboit1 = 1;
okboit2 = 1;
okboit3 = 1;
if(okboit1){
*norm2=*norm2+valxr[iiii]*valxr[iiii];
}
if(okboit2){
*norm3=*norm3+valxr[iiii]*valxr[iiii];
}
if(okboit3){
*norm4=*norm4+valxr[iiii]*valxr[iiii];
}
for(kko=kk-2;kko<=kk+(*natom/2-3);kko++){
for(kok=kko;kok<=kk+(*natom/2-3);kok++){
kko2=kko;
if(kko2<0){
kko2 = (*natom/2) + (kko2)%(*natom/2);
}
else{
kko2 = (kko2)%(*natom/2);
}
kok2=kok;
if(kok2<0){
kok2 = (*natom/2) + (kok2)%(*natom/2);
}
else{
kok2 = (kok2)%(*natom/2);
}
if(kok == kko && ideter[kok2] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko2] == 1 && ideter[kok2] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko2] == 2 && ideter[kok2] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko2] == 1 && ideter[kok2] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko2]=2;
ideter2[kok2]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko2] == 2 && ideter[kok2] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko2]=1;
ideter2[kok2]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
for(kko=*natom+kk-2;kko<=*natom+kk+(*natom/2-3);kko++){
for(kok=kko;kok<=*natom+kk+(*natom/2-3);kok++){
kko2=kko-*natom;
if(kko2<0){
kko2 = (*natom/2) + (kko2)%(*natom/2);
}
else{
kko2 = (kko2)%(*natom/2);
}
kok2=kok-*natom;
if(kok2<0){
kok2 = (*natom/2) + (kok2)%(*natom/2);
}
else{
kok2 = (kok2)%(*natom/2);
}
kko2 = (*natom) - 1 - kko2;
kok2 = (*natom) - 1 - kok2;
if(kok == kko && ideter[kok2] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko2] == 1 && ideter[kok2] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko2] == 2 && ideter[kok2] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko2] == 1 && ideter[kok2] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko2]=2;
ideter2[kok2]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko2] == 2 && ideter[kok2] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko2]=1;
ideter2[kok2]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
for(kko=kk-2;kko<=kk+(*natom/2-3);kko++){
for(kok=*natom + kk-2;kok<=*natom + kk+(*natom/2-3);kok++){
kko2=kko;
if(kko2<0){
kko2 = (*natom/2) + (kko2)%(*natom/2);
}
else{
kko2 = (kko2)%(*natom/2);
}
kok2=kok-*natom;
if(kok2<0){
kok2 = (*natom/2) + (kok2)%(*natom/2);
}
else{
kok2 = (kok2)%(*natom/2);
}
kok2 = (*natom) - 1 - kok2;
if(kok == kko && ideter[kok2] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko2] == 1 && ideter[kok2] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko2] == 2 && ideter[kok2] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko2] == 1 && ideter[kok2] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko2]=2;
ideter2[kok2]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko2] == 2 && ideter[kok2] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko2]=1;
ideter2[kok2]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
// if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);}
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
*xymat=*xymat+xmat;
*xymat2=*xymat2+xmat2;
*xymat3=*xymat3+xmat3;
*xymat4=*xymat4+xmat4;
// if(mpiid==0)printf(" ii = %d xmat3 = %18f xmat4 = %18f diff = %18f\n", ii, xmat3, xmat4, (xmat3-xmat4));
}
ierr = PetscTime(&tt2);
//if(mpiid==0)printf(" norm3 = %18f norm4 = %18f xymat3= %18f xymat4= %18f\n", *norm3, *norm4, *xymat3, *xymat4);
//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr);
}

682
src/get_s2_mov.c Normal file
View File

@ -0,0 +1,682 @@
#include <stdio.h>
#include <petsctime.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include "get_s2.h"
#include "get_val_iaa2.h"
/*
* This function simply calculates the S^2 value of the wavefunction
* Input
* =====
* Vr = The full vector
* Istart = Local starting id of the vector
* Iend = Local vector ending id
* valxr = Local vector values
* natom = number of orbitals
* Output
* ======
* norm = norm of the vector
* xymat = the S^2 value
*/
void get_s2_mov(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){
long int iaa2, iaa;
long int iii;
int ideter[natomax];
int ideter2[natomax];
int kko,kok,kkio;
long int ii;
double xmat=0.0;
double xmat2=0.0;
double xmat3=0.0;
double xmat4=0.0;
double getvaliaa2;
PetscLogDouble t1,t2,tt1,tt2;
PetscErrorCode ierr;
PetscInt iiii;
int ntrouboit1=0;
int ntrouboit2=0;
int ntrouboit3=0;
int okboit1=0;
int okboit2=0;
int okboit3=0;
int mpiid;
int pos1=0;
int pos2=0;
int pos3=0;
MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
//if(!mpiid){printf("istart= %d ind = %d\n",*Istart,*Iend);}
//ierr = PetscTime(&tt1);CHKERRQ(ierr);
for(ii=*Istart;ii<*Iend;ii++) {
iii = ii + 1;
// iiii = ii-*Istart;
iiii = ii;
xmat = 0.0;
xmat2 = 0.0;
xmat3 = 0.0;
xmat4 = 0.0;
ntrouboit1 = 0;
ntrouboit2 = 0;
ntrouboit3 = 0;
okboit1 = 0;
okboit2 = 0;
okboit3 = 0;
pos1 = 0;
pos2 = 0;
pos3 = 0;
getdet_(&iii, ideter);
*norm=*norm+valxr[iiii]*valxr[iiii];
for(kko=*s21a1;kko<=*s21a2;kko++){
if(ideter[kko]==3){
ntrouboit1++;
pos1=kko;
}
}
for(kko=*s22a1;kko<=*s22a2;kko++){
if(ideter[kko]==3){
ntrouboit2++;
pos2=kko;
}
}
for(kko=*s23a1;kko<=*s23a2;kko++){
if(ideter[kko]==3){
ntrouboit3++;
pos3=kko;
}
}
if(ntrouboit1==1 && *s21a1 <= pos1 && pos1 <= *s21a2)okboit1=1;
if(ntrouboit2==1 && pos2 == *postrou)okboit2=1;
if(ntrouboit3==1 && pos3 == *postrou)okboit3=1;
if(okboit1){
*norm2=*norm2+valxr[iiii]*valxr[iiii];
}
if(okboit2){
*norm3=*norm3+valxr[iiii]*valxr[iiii];
}
if(okboit3){
*norm4=*norm4+valxr[iiii]*valxr[iiii];
}
/*
* calculate the weight of ms=5/2
*
* loop over the determinants to see if we have a S=5/2
*/
int countw = 0;
for(kko=*s21a1;kko<=*s21a2;kko++){
if(ideter[kko] == 2) countw=1;
}
for(kok=*s21b1;kok<=*s21b2;kok++){
if(ideter[kok] == 2) countw=1;
}
if(countw==0 && okboit1){
*weight3 += (valxr[iiii]*valxr[iiii]);
}
for(kko=0;kko<=(*natom/2)-1;kko++){
for(kok=kko;kok<=(*natom/2)-1;kok++){
if(kok == kko && ideter[kok] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko] == 1 && ideter[kok] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 1 && ideter[kok] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=2;
ideter2[kok]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=1;
ideter2[kok]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21a1 && kok <=*s21a2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22a1 && kok <=*s22a2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23a1 && kok <=*s23a2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
for(kko=(*natom/2);kko<=*natom-1;kko++){
for(kok=kko;kok<=*natom-1;kok++){
if(kok == kko && ideter[kok] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko] == 1 && ideter[kok] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 1 && ideter[kok] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=2;
ideter2[kok]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=1;
ideter2[kok]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21b1 && kko <=*s21b2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22b1 && kko <=*s22b2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23b1 && kko <=*s23b2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
for(kko=0;kko<=(*natom/2)-1;kko++){
for(kok=(*natom/2);kok<=*natom-1;kok++){
if(kok == kko && ideter[kok] != 3){
xmat=xmat+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(3.0/4.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
else{
if(ideter[kko] == 1 && ideter[kok] == 1){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 2){
xmat=xmat+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
}
if(ideter[kko] == 1 && ideter[kok] == 2){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=2;
ideter2[kok]=1;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
if(ideter[kko] == 2 && ideter[kok] == 1){
xmat=xmat-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4-(1.0/2.0)*(valxr[iiii]*valxr[iiii]);
}
}
}
for(kkio=0;kkio<=*natom-1;kkio++){
ideter2[kkio]=ideter[kkio];
}
ideter2[kko]=1;
ideter2[kok]=2;
adr_(ideter2, &iaa2);
iaa2 = iaa2 - 1;
// if(!mpiid){if(iaa2 > *Iend || iaa2 < *Istart)printf("out iaa2 = %d\n",iaa2);}
xmat=xmat+valxr[iiii]*valxr[iaa2];
if(okboit1){
if( kko >=*s21a1 && kko <=*s21a2){
if( kok >=*s21b1 && kok <=*s21b2){
xmat2=xmat2+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit2){
if( kko >=*s22a1 && kko <=*s22a2){
if( kok >=*s22b1 && kok <=*s22b2){
xmat3=xmat3+(valxr[iiii]*valxr[iaa2]);
}
}
}
if(okboit3){
if( kko >=*s23a1 && kko <=*s23a2){
if( kok >=*s23b1 && kok <=*s23b2){
xmat4=xmat4+(valxr[iiii]*valxr[iaa2]);
}
}
}
}
}
}
}
*xymat=*xymat+xmat;
*xymat2=*xymat2+xmat2;
*xymat3=*xymat3+xmat3;
*xymat4=*xymat4+xmat4;
// if(mpiid==3)printf(" ii = %d norm = %18f %18f 3 = %18f 4 = %18f\n", ii, *norm2, *norm3, *xymat2, *xymat3);
}
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\n", *norm2, *norm3, *xymat2, *xymat3);
//ierr = PetscPrintf(PETSC_COMM_WORLD," Time used for the s2 loop: %f\n",tt2-tt1);CHKERRQ(ierr);
}

34
src/get_val_iaa2.c Normal file
View File

@ -0,0 +1,34 @@
#include <stdio.h>
#include <petsctime.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include "get_val_iaa2.h"
/*
* This function gets vector from a different processor
* xr the full vector
* iaa2 adresse to get value from
* getvaliaa2 the value
*/
void get_val_iaa2(Vec xr,long int *iaa2,double *getvaliaa2){
Vec x; /* initial vector, destination vector */
VecScatter scatter; /* scatter context */
IS from, to; /* index sets that define the scatter */
PetscScalar *values;
int idx_from[] = {*iaa2}, idx_to[] = {0};
VecCreateSeq( PETSC_COMM_SELF,1,&x);
ISCreateGeneral(PETSC_COMM_SELF,1,idx_from,PETSC_COPY_VALUES,&from);
ISCreateGeneral(PETSC_COMM_SELF,1,idx_to, PETSC_COPY_VALUES,&to);
printf("in get_val");
VecScatterCreate(xr,from,x,to,&scatter);
VecScatterBegin(scatter,xr,x,INSERT_VALUES,SCATTER_FORWARD);
VecScatterEnd(scatter,xr,x, INSERT_VALUES,SCATTER_FORWARD);
VecGetArray(x,&values);
*getvaliaa2 = values[0];
ISDestroy(&from);
ISDestroy(&to);
VecScatterDestroy(&scatter);
}

8
src/get_val_iaa2.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_val_iaa2(Vec, long int *, double *);

73
src/getdet.irp.f
View File

@ -1,13 +1,14 @@
subroutine getdet(add,ideter)
use iso_c_binding
implicit none
BEGIN_DOC
! this routing gives the determinant in
! the traditional form given its address
END_DOC
integer,INTENT(INOUT)::ideter(natomax)
integer(kind=selected_int_kind(16)),INTENT(IN)::add
integer(kind=selected_int_kind(16))::deta,detb
integer::i,const,ia,ib
integer(C_SIZE_T),INTENT(IN)::add
integer(C_SIZE_T)::deta,detb
integer::i,const,ia,ib, natom2
ib = MOD(add,nt2)
if(MOD(add,nt2).eq.0)then
@ -20,33 +21,49 @@ subroutine getdet(add,ideter)
detb=0
deta=0
i=1
do while (i.le.(ib))
const=1
do while(popcnt(detb).ne.nbeta .or. const==1)
if(nbeta.eq.0)then
detb=0
EXIT
endif
detb+=1
const=0
enddo
i+=1
! write(6,14)detb,detb
enddo
i=1
do while (i.le.(ia))
const=1
do while(popcnt(deta).ne.ntrou .or. const==1)
deta+=1
const=0
enddo
i+=1
! write(6,14)deta,deta
enddo
detb = det(ib,1)
deta = deth(ia,1)
if(FAM1) then
if(fix_trou1 .eq. fix_trou2) then
deta = ISHFT(deta,-(natom/2))
else
natom2 = natom - (fix_trou2 - fix_trou1)
deta = ISHFT(deta, -natom2)
endif
endif
! do while (i.le.(ib))
! const=1
! do while(popcnt(detb).ne.nbeta .or. const==1)
! detb+=1
! const=0
! enddo
! i+=1
! write(6,14)detb,detb
! enddo
! i=1
! do while (i.le.(ia))
! const=1
! do while(popcnt(deta).ne.ntrou .or. const==1)
! deta+=1
! const=0
! enddo
! i+=1
! write(6,14)deta,deta
! enddo
const=0
do i=0,(natom/2) - 1
if(FAM1) then
if(fix_trou1 .eq. fix_trou2) then
natom2 = natom/2
else
natom2 = (fix_trou2 - fix_trou1)
endif
else
natom2 = natom
endif
do i=0,(natom2) - 1
if(BTEST(deta,i))then
ideter((natom/2)-i)=3
ideter((natom2)-i)=3
endif
enddo
do i=0,natom-1

9
src/main.irp.f
View File

@ -1,13 +1,14 @@
program main
use iso_c_binding
implicit none
integer(kind=selected_int_kind(16)),allocatable::tl1(:),tl2(:),tktyp(:)
integer(C_SIZE_T),allocatable::tl1(:),tl2(:),tktyp(:)
real*8,allocatable::txtt(:),txjjxy(:),txjjz(:)
integer::i, tnrows, tntrou,tisz
real*4::t1,t2
real*8,allocatable::tval(:)
integer(kind=selected_int_kind(16)),allocatable::tcol(:)
integer(kind=selected_int_kind(16)),dimension(10)::tcountcol
integer(kind=selected_int_kind(16))::tistart
integer(C_SIZE_T),allocatable::tcol(:)
integer(C_SIZE_T),dimension(10)::tcountcol
integer(C_SIZE_T)::tistart
allocate (tl1(maxlien))
allocate (tl2(maxlien))
allocate (tktyp(maxlien))

7
src/natom.irp.f
View File

@ -2,7 +2,6 @@ BEGIN_PROVIDER [integer, natom]
&BEGIN_PROVIDER [integer, natrest]
&BEGIN_PROVIDER [integer, ial0]
&BEGIN_PROVIDER [logical*1, yham]
&BEGIN_PROVIDER [logical*1, FAM1]
&BEGIN_PROVIDER [integer, nlientot]
&BEGIN_PROVIDER [real*8, xt,(maxlien)]
&BEGIN_PROVIDER [real*8 , xjz,(maxlien)]
@ -93,10 +92,14 @@ BEGIN_PROVIDER [integer, natom]
enddo
!------------------Lecture Hamiltonien
FAM1=.TRUE.
! FAM1=.TRUE.
yham=.TRUE.
write(6,*)'HAMILTONIEN t-J'
write(6,*)'Le nombre de trou est : ',ntrou
write(6,*)'Famille 1 : ',FAM1
if(FAM1) then
if(fix_trou1 .ne. fix_trou2) write(6,*)'Trou fixe entre :', fix_trou1, "et ", fix_trou2
endif
!---------------------------------------------
write(6,*)' '
write(6,*)' '

4
src/natomax.irp.f
View File

@ -13,10 +13,10 @@ BEGIN_PROVIDER [integer, natomax]
END_DOC
implicit none
natomax=700
natomax=900
jrangmax=10705432
maxial=20
maxlien=700
maxlien=900
maxplac=20
maxdet=10000
END_PROVIDER

15
src/nt1.irp.f
View File

@ -1,13 +1,20 @@
BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt1]
use iso_c_binding
BEGIN_PROVIDER [integer(C_SIZE_T), nt1]
BEGIN_DOC
! calculates the number of det the 3's moving
END_DOC
implicit none
integer(kind=selected_int_kind(16))::natom2
integer(C_SIZE_T)::natom2
! call combin(idet1(1,nt1+1),natom,ntrou,nt1,32,jrangmax)
natom2=natom
if(FAM1)natom2=natom/2
nt1= nint(gamma(real(natom2+1,16))/(gamma(real(natom2-ntrou+1,16))*gamma(real(ntrou+1,16))),selected_int_kind(16))
if(FAM1) then
if(fix_trou1 .eq. fix_trou2) then
natom2=natom/2
else
natom2 = fix_trou2 - fix_trou1
endif
endif
nt1= nint(gamma(1.0*(natom2+1))/(gamma(1.0*(natom2-ntrou+1))*gamma(1.0*(ntrou+1))),selected_int_kind(16))
write(6,*)'nt1',nt1
END_PROVIDER

5
src/nt2.irp.f
View File

@ -1,10 +1,11 @@
BEGIN_PROVIDER [integer(kind=selected_int_kind(16)), nt2]
use iso_c_binding
BEGIN_PROVIDER [integer(C_SIZE_T), nt2]
BEGIN_DOC
! calculates the number of det the 1's moving
END_DOC
! call combin(idet2(1,nt2+1),natrest,ial0,nt2,32,jrangmax)
nt2= nint(gamma(real(natom-ntrou+1,16))/((gamma(real(nalpha+1,16))*gamma(real(nbeta+1,16)))),selected_int_kind(16))
nt2= nint(gamma(1.0*(natom-ntrou+1))/((gamma(1.0*(nalpha+1))*gamma(1.0*(nbeta+1)))),selected_int_kind(16))
print *,"nt2=",nt2
END_PROVIDER

91
src/providdet.irp.f Normal file
View File

@ -0,0 +1,91 @@
use iso_c_binding
BEGIN_PROVIDER[integer(C_SIZE_T),det,(nt2,2)]
&BEGIN_PROVIDER[integer(C_SIZE_T),deth,(nt1,2)]
BEGIN_DOC
! provides det and deth array
END_DOC
use iso_c_binding
implicit none
! integer(kind=selected_int_kind(16))::dethsh
integer(C_SIZE_t)::a
integer(C_SIZE_T)::i,count
integer::const
i=1
a=0
const=0
count=0
If(ntrou.ge.1)then
const=0
! dethsh = ISHFT(deth,-natom/2)
! i=nt1
do while (i.le.(nt1))
! if(a.eq.dethsh)then
! addh=i-1
! EXIT
! endif
i+=1
a+=1
do while(popcnt(a).ne.ntrou)
a+=1
enddo
count+=1
if(FAM1) then
if(fix_trou1 .eq. fix_trou2) then
deth(count,1)=ISHFT(a,natom/2)
else
deth(count,1)=ISHFT(a,natom - (fix_trou2-fix_trou1))
endif
else
deth(count, 1) = a
endif
deth(count,2)=i-1
! write(6,16)ISHFT(a,natom/2),ISHFT(a,natom/2),i-1
enddo
! if(a.eq.dethsh )then
! count+=1
! deth(1,1)=ISHFT(a,natom/2)
! deth(1,2)=nt1
! endif
endif
!C if det=0 then exit
a=0
i=0
count=0
print *,'nt2=',nt2,'nbeta=',nbeta
do while (i.lt.(nt2))
i+=1
a+=1
do while(popcnt(a).ne.nbeta)
if(nbeta.eq.0)then
a=0
count+=1
det(count,1)=a
det(count,2)=i
EXIT
endif
a+=1
enddo
if(nbeta.ne.0)then
count+=1
det(count,1)=a
det(count,2)=i
endif
! write(6,16)a,a,i
enddo
10 FORMAT(B64,I8,F8.2)
15 FORMAT(B64,I8,I8,I8)
11 FORMAT(B64,I3,B64)
12 FORMAT(I5,$)
13 FORMAT(B64,B64)
14 FORMAT(B64,I8)
16 FORMAT(B64,I8,I8)
END_PROVIDER

10
src/provide_clust.irp.f
View File

@ -1,11 +1,15 @@
BEGIN_PROVIDER[integer,l1, (maxlien)]
&BEGIN_PROVIDER[integer,l2, (maxlien)]
&BEGIN_PROVIDER[integer,ktyp,(maxlien)]
&BEGIN_PROVIDER [real*8, xtt ,(maxlien)]
&BEGIN_PROVIDER[real*8, xtt ,(maxlien)]
&BEGIN_PROVIDER[real*8, xjjz ,(maxlien)]
&BEGIN_PROVIDER[real*8, xjjxy,(maxlien)]
&BEGIN_PROVIDER [integer, ntrou]
&BEGIN_PROVIDER [integer, isz]
&BEGIN_PROVIDER[real*8, E,(maxlien)]
&BEGIN_PROVIDER[integer, ntrou]
&BEGIN_PROVIDER[integer, isz]
&BEGIN_PROVIDER[logical*1, FAM1]
&BEGIN_PROVIDER[integer, fix_trou1]
&BEGIN_PROVIDER[integer, fix_trou2]
implicit none
! integer::i
! open(unit=11,file="l1.dat",form="formatted")

114
src/read2.c
View File

@ -1,7 +1,3 @@
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include "read2.h"
void Data_new(FILE* file, Data* dat) {
@ -15,13 +11,13 @@ void Data_new(FILE* file, Data* dat) {
while (fgets(line, sizeof(line), file)) {
/* note that fgets don't strip the terminating \n, checking its
presence would allow to handle lines longer that sizeof(line) */
if (count != 11){
/* note that fgets doesn't strip the terminating \n, checking its
presence would allow to handle lines longer than sizeof(line) */
if (count != 30){
count++;
switch(count){
case 1:
dat->n=atol(line);
dat->natom=atol(line);
break;
case 2:
dat->nnz=atol(line);
@ -36,6 +32,9 @@ void Data_new(FILE* file, Data* dat) {
dat->isz=atol(line);
break;
case 6:
dat->FAM1 = to_bool(line);
break;
case 7:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
@ -59,7 +58,7 @@ void Data_new(FILE* file, Data* dat) {
}
}
break;
case 7:
case 8:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
@ -83,7 +82,7 @@ void Data_new(FILE* file, Data* dat) {
}
}
break;
case 8:
case 9:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
@ -107,7 +106,7 @@ void Data_new(FILE* file, Data* dat) {
}
}
break;
case 9:
case 10:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
@ -131,7 +130,7 @@ void Data_new(FILE* file, Data* dat) {
}
}
break;
case 10:
case 11:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
@ -155,7 +154,7 @@ void Data_new(FILE* file, Data* dat) {
}
}
break;
case 11:
case 12:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
@ -179,6 +178,84 @@ void Data_new(FILE* file, Data* dat) {
}
}
break;
case 13:
arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
double val;
char *unconverted;
/**
* Convert the next token to a float value
*/
val = strtof(token, &unconverted);
if (!isspace(*unconverted) && *unconverted != 0)
{
/**
* Bad input string. Again, we just bail.
*/
fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
break;
}
else
{
dat->E[arrayIdx++] = val;
}
}
break;
case 14:
dat->nroots=atol(line);
break;
case 15:
dat->s21a1=atol(line);
break;
case 16:
dat->s21a2=atol(line);
break;
case 17:
dat->s21b1=atol(line);
break;
case 18:
dat->s21b2=atol(line);
break;
case 19:
dat->s22a1=atol(line);
break;
case 20:
dat->s22a2=atol(line);
break;
case 21:
dat->s22b1=atol(line);
break;
case 22:
dat->s22b2=atol(line);
break;
case 23:
dat->s23a1=atol(line);
break;
case 24:
dat->s23a2=atol(line);
break;
case 25:
dat->s23b1=atol(line);
break;
case 26:
dat->s23b2=atol(line);
break;
case 27:
dat->postrou=atol(line);
break;
case 28:
dat->fix_trou1=atol(line);
break;
case 29:
dat->fix_trou2=atol(line);
break;
case 30:
dat->print_wf = atol(line);
break;
default:
printf("Done reading file\n");
break;
} /* end of switch */
} /* end of the input file */
@ -188,6 +265,17 @@ void Data_new(FILE* file, Data* dat) {
//return dat;
}
_Bool to_bool(const char* str) {
PetscBool strflg=PETSC_FALSE;
PetscStrcmp("true\n",str, &strflg);
if(strflg != PETSC_TRUE) PetscStrcmp("True\n",str, &strflg);
if(strflg != PETSC_TRUE) PetscStrcmp("TRUE\n",str, &strflg);
if(strflg != PETSC_TRUE) PetscStrcmp("true",str, &strflg);
if(strflg != PETSC_TRUE) PetscStrcmp("True",str, &strflg);
if(strflg != PETSC_TRUE) PetscStrcmp("TRUE",str, &strflg);
return strflg;
}
/*
int main(int argc, char* argv[])
{

38
src/read2.h
View File

@ -1,19 +1,43 @@
#include <stdio.h>
#include <slepceps.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <petscsys.h>
#include <slepceps.h>
_Bool to_bool(const char* str);
typedef struct {
PetscInt n;
long int nnz,npar;
long int ntrou,isz;
long int l1[700];
long int l2[700];
long int ktyp[700];
double xjjz[700];
double xjjxy[700];
double xtt[700];
_Bool FAM1;
long int l1[900];
long int l2[900];
long int ktyp[900];
double xjjz[900];
double xjjxy[900];
double xtt[900];
double E[900];
long int nroots;
int natom;
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;
long int fix_trou1;
long int fix_trou2;
long int print_wf;
} Data ;

117
src/searchdet.irp.f
View File

@ -1,80 +1,71 @@
subroutine searchdet(det,add,deth,addh)
subroutine searchdet(deti,add,dethi,addh)
use iso_c_binding
BEGIN_DOC
! this subroutine is at the heart of the idea
! it will generate all the determinants in a fixed order
! then find the posistion of the determinant given and
! return it's position in add.
END_DOC
integer(kind=selected_int_kind(16)),INTENT(INOUT)::det
integer(kind=selected_int_kind(16)),INTENT(INOUT)::add
integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth
integer(kind=selected_int_kind(16)),INTENT(INOUT)::addh
integer(kind=selected_int_kind(16))::dethsh
integer(kind=selected_int_kind(16))::a
integer(kind=selected_int_kind(16))::i
integer::const
integer(C_SIZE_T),INTENT(INOUT)::deti
integer(C_SIZE_T),INTENT(INOUT)::add
integer(C_SIZE_T),INTENT(INOUT)::dethi
integer(C_SIZE_T),INTENT(INOUT)::addh
integer(C_SIZE_T)::dethsh
integer(C_SIZE_T)::a
integer(C_SIZE_T)::i,j
integer::count
logical::found
i=1
a=0
add=0
const=0
If(ntrou.ge.1)then
const=0
dethsh = ISHFT(deth,-natom/2)
addh=0
! i=nt1
do while (i.le.(nt1))
if(a.eq.dethsh)then
addh=i-1
EXIT
endif
i+=1
a+=1
do while(popcnt(a).ne.ntrou)
a+=1
enddo
enddo
if(a.eq.dethsh .and. addh.eq.0)then
addh=nt1
j=nt1
found=.FALSE.
do while(.not.found)
if(deth((i+j)/2,1).eq.dethi)then
addh=deth((i+j)/2,2)
found=.TRUE.
EXIT
elseif (abs(i-j).eq.1)then
if(deth(i,1).eq.dethi)then
addh=deth(i,2)
elseif(deth(j,1).eq.dethi)then
addh=deth(j,2)
endif
endif
!C if det=0 then exit
a=0
i=0
count=0
if(a.eq.det)then
add=1
Return
endif
do while (i.le.(nt2))
if(a.eq.det)then
if(a.eq.1)then
add=i
count=-1
EXIT
found=.TRUE.
EXIT
else
if(deth((i+j)/2,1).gt.dethi)then
j=(i+j)/2
else
add=i
count=-1
EXIT
i=(i+j)/2
endif
endif
i+=1
a+=1
!C write(6,16)a,a,i-2
do while(popcnt(a).ne.nbeta)
a+=1
enddo
enddo
if(a.eq.det .and. count.ne.-1)then
add=i-1
endif
i=1
j=nt2
found=.FALSE.
do while(.not.found)
if(det((i+j)/2,1).eq.deti)then
add=det((i+j)/2,2)
found=.TRUE.
EXIT
elseif (abs(i-j).eq.1)then
if(det(i,1).eq.deti)then
add=det(i,2)
elseif(det(j,1).eq.deti)then
add=det(j,2)
endif
found=.TRUE.
EXIT
else
if(det((i+j)/2,1).gt.deti)then
j=(i+j)/2
else
i=(i+j)/2
endif
endif
enddo
10 FORMAT(B64,I8,F8.2)
15 FORMAT(B64,I8,I8,I8)

15
src/searchdetfull.irp.f
View File

@ -1,17 +1,18 @@
subroutine searchdetfull()
use iso_c_binding
BEGIN_DOC
! this subroutine is at the heart of the idea
! it will generate all the determinants in a fixed order
! then find the posistion of the determinant given and
! return it's position in add.
END_DOC
! integer(kind=selected_int_kind(16)),INTENT(INOUT)::foundetadr(maxlien,4)
integer(kind=selected_int_kind(16))::add
! integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth
integer(kind=selected_int_kind(16))::dethsh
integer(kind=selected_int_kind(16))::addh
integer(kind=selected_int_kind(16))::a
integer(kind=selected_int_kind(16))::i
! integer(C_SIZE_T),INTENT(INOUT)::foundetadr(maxlien,4)
integer(C_SIZE_T)::add
! integer(C_SIZE_T),INTENT(INOUT)::deth
integer(C_SIZE_T)::dethsh
integer(C_SIZE_T)::addh
integer(C_SIZE_T)::a
integer(C_SIZE_T)::i
integer::const,count
i=1
a=0

11
src/sort.irp.f
View File

@ -1,14 +1,15 @@
subroutine sort()
use iso_c_binding
implicit none
integer::i,j,ord,ordh
integer(kind=selected_int_kind(16))::add,addh,det,deth,addt
integer(C_SIZE_T)::add,addh,deti,dethi,addt
do i=1,detfound-1
do j=i+1,detfound
if(foundaddh(i,1).gt.foundaddh(j,1))then
deth = foundaddh(i,1)
dethi = foundaddh(i,1)
foundaddh(i,1) = foundaddh(j,1)
foundaddh(j,1) = deth
foundaddh(j,1) = dethi
addh = foundaddh(i,2)
foundaddh(i,2) = foundaddh(j,2)
foundaddh(j,2) = addh
@ -17,9 +18,9 @@ subroutine sort()
foundaddh(j,3) = ordh
endif
if(foundadd(i,1).gt.foundadd(j,1))then
det = foundadd(i,1)
deti = foundadd(i,1)
foundadd(i,1) = foundadd(j,1)
foundadd(j,1) = det
foundadd(j,1) = deti
add = foundadd(i,2)
foundadd(i,2) = foundadd(j,2)
foundadd(j,2) = add

6
src/stimsyr.h
View File

@ -1,3 +1,5 @@
#include <petscsys.h>
void unit_l1_(
long int *,
long int *,
@ -7,9 +9,13 @@ void unit_l1_(
double *,
double *,
double *,
double *,
long int *,
long int *,
long int *,
long int *,
long int *,
_Bool *,
long int *,
double *);

12
src/unit_FIL44.irp.f
View File

@ -3,19 +3,19 @@
BEGIN_DOC
! file units for writing
END_DOC
use iso_c_binding
implicit none
integer :: i,j,k,ia1,ia2,l,m,chcind,chcval,ii,tistart2
integer :: count,unit_44,unit_33
integer :: iat,nbtots
integer(kind=selected_int_kind(16))::iaa
integer(C_SIZE_T)::iaa
integer :: kkio,kkiok,n,nz,cdiag,cexdiag
integer,allocatable ::ideter1(:),ideter2(:),deti(:),detj(:)
integer(kind=selected_int_kind(16)),dimension(maxlien) ::tl1,tl2,tktyp
integer(kind=selected_int_kind(16)),dimension(nrows)::tcountcol
integer(kind=selected_int_kind(16))::tistart
integer(C_SIZE_T),dimension(maxlien) ::tl1,tl2,tktyp
integer(C_SIZE_T),dimension(nrows)::tcountcol
integer(C_SIZE_T)::tistart
real*8,dimension(maxlien)::tval
integer(kind=selected_int_kind(16)),dimension(maxlien)::tcol
integer(C_SIZE_T),dimension(maxlien)::tcol
real*8 :: xmat
integer :: ik,imat4,iaa2,iik
integer :: ik1,ik2,jmat4,IC,ikmax,ikmin

23
src/unit_l1.irp.f
View File

@ -6,18 +6,26 @@
txjjxy, &
txjjz , &
txtt , &
tE , &
tcountcol, &
tntrou, &
tisz, &
tfix_trou1, &
tfix_trou2, &
tfam1, &
tcol,tval)
use iso_c_binding
implicit none
integer,INTENT(INOUT)::tistart, tnrows, tntrou, tisz
integer,INTENT(INOUT)::tistart, tnrows
integer,INTENT(INOUT)::tntrou, tisz
integer,INTENT(INOUT)::tfix_trou1, tfix_trou2
logical*1,INTENT(INOUT)::tfam1
integer::i
real*8,INTENT(INOUT)::tval(maxlien)
integer(kind=selected_int_kind(16)),INTENT(INOUT)::tcol(maxlien)
integer(kind=selected_int_kind(16)),INTENT(INOUT),dimension(tnrows)::tcountcol
integer(kind=selected_int_kind(16)),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien)
real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien)
integer(C_SIZE_T),INTENT(INOUT)::tcol(maxlien)
integer(C_SIZE_T),INTENT(INOUT),dimension(tnrows)::tcountcol
integer(C_SIZE_T),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien)
real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien), tE(maxlien)
nrows = tnrows
provide nrows
@ -28,9 +36,13 @@
xtt(i) = txtt(i)
xjjxy(i) = txjjxy(i)
xjjz (i) = txjjz (i)
E (i) = tE (i)
enddo
ntrou = tntrou
isz = tisz
FAM1 = tfam1
fix_trou1 = tfix_trou1
fix_trou2 = tfix_trou2
tcol=0
tval=0d0
provide l1 l2 ktyp xtt xjjxy xjjz ntrou
@ -38,6 +50,7 @@
!print *,l1
!print *,"xjjz"
!print *,xjjz
!print *,FAM1
call unit(tistart, tcountcol,tcol,tval)
end