1
ckbn ccsd_t.F modified to tce_mrcc_ccsdpt_uncoup_pt3.F
3
SUBROUTINE tce_mrcc_ccsdpt(d_t1,k_t1_offset,
8
ckbn + energy1,energy2,energy3,
14
#include "mafdecls.fh"
28
#include "tce_main.fh"
29
#include "tce_hetio.fh"
30
#include "tce_diis.fh"
31
#include "tce_prop.fh"
32
#include "tce_restart.fh"
33
#include "tce_mrcc.fh"
49
integer k_singles,l_singles
50
integer k_doubles,l_doubles
52
c integer k_q3fnt2,l_q3fnt2
58
integer offset_p4,offset_p5,offset_p6
59
integer offset_h1,offset_h2,offset_h3
60
integer range_p4,range_p5,range_p6
61
integer range_h1,range_h2,range_h3
62
double precision energy(3)
63
ckbn double precision energy1,energy2
65
double precision energy3
66
double precision factor,denom
67
double precision denom_p4,denom_p5,denom_p6
68
double precision denom_h1,denom_h2,denom_h3
72
integer d_f1,k_f1_offset
73
double precision dnorm1,dnorm2,dmaxt1,dmaxt2
74
double precision dmint1,dmint2
76
integer orbspin(6),orbindex(6),off(6),blck(6)
81
nodezero=(ga_nodeid().eq.0)
84
c if(nodezero) write(*,*)"Entering tce_mrcc_ccsdpt"
88
next = nxtask(nprocs,1)
98
do t_p4b = noab+1,noab+nvab
99
range_p4 = int_mb(k_range+t_p4b-1)
100
offset_p4 = k_evl_sorted+int_mb(k_offset+t_p4b-1)-1
101
do t_p5b = t_p4b,noab+nvab
102
range_p5 = int_mb(k_range+t_p5b-1)
103
offset_p5 = k_evl_sorted+int_mb(k_offset+t_p5b-1)-1
104
do t_p6b = t_p5b,noab+nvab
105
range_p6 = int_mb(k_range+t_p6b-1)
106
offset_p6 = k_evl_sorted+int_mb(k_offset+t_p6b-1)-1
108
range_h1 = int_mb(k_range+t_h1b-1)
109
offset_h1 = k_evl_sorted+int_mb(k_offset+t_h1b-1)-1
110
do t_h2b = t_h1b,noab
111
range_h2 = int_mb(k_range+t_h2b-1)
112
offset_h2 = k_evl_sorted+int_mb(k_offset+t_h2b-1)-1
113
do t_h3b = t_h2b,noab
114
range_h3 = int_mb(k_range+t_h3b-1)
115
offset_h3 = k_evl_sorted+int_mb(k_offset+t_h3b-1)-1
116
if (next.eq.count) then
117
if (int_mb(k_spin+t_p4b-1)
118
1 +int_mb(k_spin+t_p5b-1)
119
2 +int_mb(k_spin+t_p6b-1)
120
3 .eq.int_mb(k_spin+t_h1b-1)
121
4 +int_mb(k_spin+t_h2b-1)
122
5 +int_mb(k_spin+t_h3b-1)) then
123
if ((.not.restricted).or.
124
1 (int_mb(k_spin+t_p4b-1)
125
1 +int_mb(k_spin+t_p5b-1)
126
2 +int_mb(k_spin+t_p6b-1)
127
3 +int_mb(k_spin+t_h1b-1)
128
4 +int_mb(k_spin+t_h2b-1)
129
5 +int_mb(k_spin+t_h3b-1).le.8)) then
130
if (ieor(int_mb(k_sym+t_p4b-1),
131
1 ieor(int_mb(k_sym+t_p5b-1),
132
2 ieor(int_mb(k_sym+t_p6b-1),
133
3 ieor(int_mb(k_sym+t_h1b-1),
134
4 ieor(int_mb(k_sym+t_h2b-1),
135
5 int_mb(k_sym+t_h3b-1)))))).eq.0) then
137
IF (.not.((log_mb(k_active_tmpm(iref)+t_p4b-1).eqv..true.).and.
138
1 (log_mb(k_active_tmpm(iref)+t_p5b-1).eqv..true.).and.
139
1 (log_mb(k_active_tmpm(iref)+t_p6b-1).eqv..true.).and.
140
1 (log_mb(k_active_tmpm(iref)+t_h1b-1).eqv..true.).and.
141
1 (log_mb(k_active_tmpm(iref)+t_h2b-1).eqv..true.).and.
142
1 (log_mb(k_active_tmpm(iref)+t_h3b-1).eqv..true.))) then
144
size = range_p4 * range_p5 * range_p6
145
3 * range_h1 * range_h2 * range_h3
146
if (.not.MA_PUSH_GET(mt_dbl,size,'(T) singles',l_singles,
147
1 k_singles)) call errquit('ccsd_t: MA error',1,MA_ERR)
148
if (.not.MA_PUSH_GET(mt_dbl,size,'(T) doubles',l_doubles,
149
1 k_doubles)) call errquit('ccsd_t: MA error',2,MA_ERR)
151
dbl_mb(k_singles+i) = 0.0d0
152
dbl_mb(k_doubles+i) = 0.0d0
155
call tce_mrcc_q3vnt2(dbl_mb(k_doubles),d_t2,d_v2,k_t2_offset,
156
1 k_v2_offset,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2)
157
call tce_mrcc_q3vnt1(dbl_mb(k_singles),d_t1,d_v2,k_t1_offset,
158
1 k_v2_offset,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2)
159
call tce_mrcc_q3fnt2(dbl_mb(k_singles),d_f1,d_t2,k_f1_offset,
160
+ k_t2_offset,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2)
167
if ((t_p4b.eq.t_p5b).and.(t_p5b.eq.t_p6b)) then
168
factor = factor / 6.0d0
169
else if ((t_p4b.eq.t_p5b).or.(t_p5b.eq.t_p6b)) then
170
factor = factor / 2.0d0
172
if ((t_h1b.eq.t_h2b).and.(t_h2b.eq.t_h3b)) then
173
factor = factor / 6.0d0
174
else if ((t_h1b.eq.t_h2b).or.(t_h2b.eq.t_h3b)) then
175
factor = factor / 2.0d0
178
c factor = [ 1/36, 1/18, 1/12, 1/6, 1/4, 1/3, 1/2, 1, 2]
181
do t_p4 = 1, range_p4
182
denom_p4 = dbl_mb(offset_p4+t_p4)
183
do t_p5 = 1, range_p5
184
denom_p5 = dbl_mb(offset_p5+t_p5)
185
do t_p6 = 1, range_p6
186
denom_p6 = dbl_mb(offset_p6+t_p6)
187
do t_h1 = 1, range_h1
188
denom_h1 = dbl_mb(offset_h1+t_h1)
189
do t_h2 = 1, range_h2
190
denom_h2 = dbl_mb(offset_h2+t_h2)
191
do t_h3 = 1, range_h3
192
denom_h3 = dbl_mb(offset_h3+t_h3)
195
if(.not.lusesamefock_nonit) then
196
denom = 1.0d0 / ((denom_h1 + denom_h2 + denom_h3)
197
+ - (denom_p4 + denom_p5 + denom_p6))
201
orbspin(1) = int_mb(k_spin+t_p4b-1)-1
202
orbspin(2) = int_mb(k_spin+t_p5b-1)-1
203
orbspin(3) = int_mb(k_spin+t_p6b-1)-1
204
orbspin(4) = int_mb(k_spin+t_h1b-1)-1
205
orbspin(5) = int_mb(k_spin+t_h2b-1)-1
206
orbspin(6) = int_mb(k_spin+t_h3b-1)-1
208
orbindex(1) = (1-orbspin(1)+int_mb(k_mo_indexm(iref)+
209
1 int_mb(k_offset+t_p4b-1)+t_p4-1))/2
210
orbindex(2) = (1-orbspin(2)+int_mb(k_mo_indexm(iref)+
211
1 int_mb(k_offset+t_p5b-1)+t_p5-1))/2
212
orbindex(3) = (1-orbspin(3)+int_mb(k_mo_indexm(iref)+
213
1 int_mb(k_offset+t_p6b-1)+t_p6-1))/2
214
orbindex(4) = (1-orbspin(4)+int_mb(k_mo_indexm(iref)+
215
1 int_mb(k_offset+t_h1b-1)+t_h1-1))/2
216
orbindex(5) = (1-orbspin(5)+int_mb(k_mo_indexm(iref)+
217
1 int_mb(k_offset+t_h2b-1)+t_h2-1))/2
218
orbindex(6) = (1-orbspin(6)+int_mb(k_mo_indexm(iref)+
219
1 int_mb(k_offset+t_h3b-1)+t_h3-1))/2
222
blck(j) = orbinblck(orbindex(j),orbspin(j)+1,1)
223
off(j) = offsetinblck(orbindex(j),orbspin(j)+1,1)
227
+ -dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)
228
+ +blck(1)-1)+off(1))
229
+ -dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)
230
+ +blck(2)-1)+off(2))
231
+ -dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)
232
+ +blck(3)-1)+off(3))
233
+ +dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)
234
+ +blck(4)-1)+off(4))
235
+ +dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)
236
+ +blck(5)-1)+off(5))
237
+ +dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)
238
+ +blck(6)-1)+off(6)))
239
endif ! lusesamefock_nonit
241
if((abs(denom)).gt.1.0d2) then
244
+ 'Warning:Denominator is very low. 1/D=',abs(denom)
246
call util_flush(LuOut)
250
c energy1 = energy1 + factor*denom
251
c 1 * dbl_mb(k_doubles+i)*dbl_mb(k_doubles+i)
252
c energy2 = energy2 + factor*denom*dbl_mb(k_doubles+i)
253
c 1 * (dbl_mb(k_doubles+i)+dbl_mb(k_singles+i))
254
energy3 = energy3 + factor*denom*dbl_mb(k_doubles+i)
255
+ * (dbl_mb(k_doubles+i)+dbl_mb(k_singles+i))
265
c if (.not.MA_POP_STACK(l_q3fnt2))
266
c 1 call errquit('ccsd_t',2,MA_ERR)
267
if (.not.MA_POP_STACK(l_doubles))
268
1 call errquit('ccsd_t',3,MA_ERR)
269
if (.not.MA_POP_STACK(l_singles))
270
1 call errquit('ccsd_t',4,MA_ERR)
276
next = NXTASKsub(nprocs,1,mypgid)
278
next = nxtask(nprocs,1)
289
next = NXTASKsub(-nprocs,1,mypgid)
291
next = nxtask(-nprocs,1)
294
c energy(1) = energy1
295
c energy(2) = energy2
298
call ga_pgroup_sync(mypgid)
299
c call ga_pgroup_dgop(p_handle, type, buf, lenbuf, op)
300
call ga_pgroup_dgop(mypgid, mt_dbl,energy,3,'+')
302
call ga_dgop(mt_dbl,energy,3,'+')
304
c energy1 = energy(1)
305
c energy2 = energy(2)
309
c if(ga_pgroup_nodeid(mypgid).eq.0) then
310
cc write(*,'(A,F17.10,A,F17.10,A,F17.10)')"energy1 ",energy1,
311
cc + " energy2 ",energy2," energy3 ",energy3
312
c write(*,'(A,F17.10,A,F17.10,A,F17.10)')"energy3 ",energy3
316
cc write(*,'(A,F17.10,A,F17.10,A,F17.10)')"energy1 ",energy1,
317
cc + " energy2 ",energy2," energy3 ",energy3
318
c write(*,'(A,F17.10,A,F17.10,A,F17.10)')"energy3 ",energy3