Wick-Contraction/analysis.f90 at master · jjren/Wick-Contraction · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
Module analysis
	use para_mod
	use array_mod
	implicit none

	character(len=2),allocatable :: integralterm(:,:,:)

	contains
!===================================================
!===================================================

subroutine InterSingletFission
	implicit none
	integer :: i,igroup,j,k
	integer :: link1(noperators/2),link2(noperators/2)
	logical :: fullfill

	write(*,*) "orbital specific integral"
	allocate(integralterm(4,noutputterm,ngroups))
	integralterm="no"

	do igroup=1,ngroups,1
		do i=1,noutputterm,1
			if(outputterm(i).ntermoprs==0) then
				fullfill=.true.
				do j=1,outputterm(i).ndeltas,1
					do k=1,noperators,1
						! find the groupinfo terms corresponse the specific delta
						if(outputterm(i).deltapair(1,j)==groupinfo(1,k,igroup)) then
							link1(j)=k
						else if(outputterm(i).deltapair(2,j)==groupinfo(1,k,igroup)) then
							link2(j)=k
						end if
					end do
					if((outputterm(i).deltapair(3,j)/=groupinfo(2,link1(j),igroup) .and. groupinfo(2,link1(j),igroup)/='r') .or. &
						(outputterm(i).deltapair(3,j)/=groupinfo(2,link2(j),igroup) .and. groupinfo(2,link2(j),igroup)/='r')) then
						fullfill=.false.
						exit
					end if
				end do
				if(fullfill==.true.) then
					do j=1,outputterm(i).ndeltas,1
						if((groupinfo(3,link1(j),igroup)/=groupinfo(3,link2(j),igroup) .and. groupinfo(3,link2(j),igroup)/='r' .and. groupinfo(3,link1(j),igroup)/='r') .or. &
							(groupinfo(4,link1(j),igroup)/=groupinfo(4,link2(j),igroup))) then
							fullfill=.false.
							exit
						end if
					end do
				end if
				if(fullfill==.true.) then
					write(*,*) "======================="
					write(*,*) "igroup:",igroup
					write(*,*) "outputterm:",i,"sign:",outputterm(i).sign1
					do j=1,outputterm(i).ndeltas,1
						if(groupinfo(3,link1(j),igroup)=='r' .or. groupinfo(3,link2(j),igroup)=='r') then
							if(groupinfo(3,link1(j),igroup)=='r' .and. groupinfo(3,link2(j),igroup)=='r') then
								write(*,*) groupinfo(1,link1(j),igroup),outputterm(i).deltapair(3,j)
								call pqrssort(groupinfo(1,link1(j),igroup),outputterm(i).deltapair(3,j),i,igroup)
								write(*,*) groupinfo(1,link2(j),igroup),outputterm(i).deltapair(3,j)
								call pqrssort(groupinfo(1,link2(j),igroup),outputterm(i).deltapair(3,j),i,igroup)
							else if(groupinfo(3,link1(j),igroup)/='r') then
								write(*,*) groupinfo(1,link2(j),igroup),groupinfo(3,link1(j),igroup)
								call pqrssort(groupinfo(1,link2(j),igroup),groupinfo(3,link1(j),igroup),i,igroup)
							else if(groupinfo(3,link2(j),igroup)/='r') then
								write(*,*) groupinfo(1,link1(j),igroup),groupinfo(3,link2(j),igroup)
								call pqrssort(groupinfo(1,link1(j),igroup),groupinfo(3,link2(j),igroup),i,igroup)
							end if
						end if
					end do
				end if
			end if
		end do
	end do

	if(flag=='h' .or. flag=='g') then
		call sumintegral
	end if
	deallocate(integralterm)

return
end subroutine InterSingletFission

!===================================================
!===================================================

subroutine pqrssort(oindex,orbital,termindex,igroup)
	implicit none
	character(len=2) :: oindex,orbital
	integer :: termindex,igroup

	if(oindex=='p') then
		integralterm(1,termindex,igroup)=orbital
	else if(oindex=='q') then
		integralterm(2,termindex,igroup)=orbital
	else if(oindex=='r') then
		integralterm(3,termindex,igroup)=orbital
	else if(oindex=='s') then
		integralterm(4,termindex,igroup)=orbital
	end if
return
end subroutine pqrssort

!===================================================
!===================================================

subroutine sumintegral
	implicit none
	character(len=2) :: sumintegralterm(4,ngroups*noutputterm)
	integer :: realterms
	real(kind=8) :: sumintegralcoeff(ngroups*noutputterm)
	logical :: ifexist
	integer :: igroup,i,j,k

	realterms=0
	sumintegralcoeff=0.0D0

	do igroup=1,ngroups,1
	do j=1,noutputterm,1
		if(integralterm(1,j,igroup)/='no') then
			ifexist=.false.
			do k=1,realterms,1
				if(flag=='h') then
					if((sumintegralterm(1,k)==integralterm(1,j,igroup) &
						.and. sumintegralterm(2,k)==integralterm(2,j,igroup)) .or. &
						(sumintegralterm(1,k)==integralterm(2,j,igroup) &
						.and. sumintegralterm(2,k)==integralterm(1,j,igroup))) then
						ifexist=.true.
						exit
					end if
				else if(flag=='g') then
					if (  &
						! combine the same two electron integral terms
						(((sumintegralterm(1,k)==integralterm(1,j,igroup) &
						.and. sumintegralterm(4,k)==integralterm(4,j,igroup)) .or. &
						(sumintegralterm(1,k)==integralterm(4,j,igroup) &
						.and. sumintegralterm(4,k)==integralterm(1,j,igroup))) .and. &
						((sumintegralterm(3,k)==integralterm(3,j,igroup) &
						.and. sumintegralterm(2,k)==integralterm(2,j,igroup)) .or. &
						(sumintegralterm(3,k)==integralterm(2,j,igroup) &
						.and. sumintegralterm(2,k)==integralterm(3,j,igroup)))) &
						 				.or. &
						(((sumintegralterm(1,k)==integralterm(3,j,igroup) &
						.and. sumintegralterm(4,k)==integralterm(2,j,igroup)) .or. &
						(sumintegralterm(4,k)==integralterm(3,j,igroup) &
						.and. sumintegralterm(1,k)==integralterm(2,j,igroup))) .and. &
						((sumintegralterm(2,k)==integralterm(1,j,igroup) &
						.and. sumintegralterm(3,k)==integralterm(4,j,igroup)) .or. &
						(sumintegralterm(3,k)==integralterm(4,j,igroup) &
						.and. sumintegralterm(2,k)==integralterm(1,j,igroup)))) &
										) then
						ifexist=.true.
						exit
					end if
				end if
			end do
			if(ifexist==.true.) then
				sumintegralcoeff(k)=sumintegralcoeff(k)+groupcoeff(igroup)*DBLE(outputterm(j).sign1)
			else
				realterms=realterms+1
				sumintegralterm(:,realterms)=integralterm(:,j,igroup)
				sumintegralcoeff(realterms)=groupcoeff(igroup)*DBLE(outputterm(j).sign1)
			end if
		end if
	end do
	end do

	write(*,*) "============================================================"
	write(*,*) "in the two electron term do not forget about the 1/2 factor!"
	write(*,*) "============================================================"
	do i=1,realterms,1
		write(*,*) "==========================="
		write(*,*) "integral term",i
		write(*,*) sumintegralcoeff(i)
		if(flag=='g') then
			write(*,*) "(",sumintegralterm(1,i),sumintegralterm(4,i),"|",sumintegralterm(2,i),sumintegralterm(3,i),")"
		else if(flag=='h') then
			write(*,*) "(",sumintegralterm(1,i),"|",sumintegralterm(2,i),")"
		end if
	end do

return
end subroutine sumintegral
!===================================================
!===================================================
end module