1      subroutine eomccsd_energy(d_rx1,d_rx2,size_x1,size_x2,
2     1           k_x1_offset,k_x2_offset,
3     1           d_f1,d_v2,d_t1,d_t2,k_f1_offset,k_v2_offset,
4     2           k_t1_offset,k_t2_offset,k_omegax)
5      implicit none
6#include "global.fh"
7#include "inp.fh"
8#include "mafdecls.fh"
9#include "util.fh"
10#include "errquit.fh"
11#include "stdio.fh"
12#include "tce.fh"
13#include "tce_main.fh"
14#include "tce_diis.fh"
15c
16      integer size_x1,size_x2
17      integer k_x1_offset,k_x2_offset
18      integer d_f1,d_v2
19      integer d_t1,d_t2
20      integer k_f1_offset,k_v2_offset
21      integer k_t1_offset,k_t2_offset
22c
23c      integer iter,ivec
24      integer ivec
25      logical nodezero
26      double precision cpu, wall
27      logical needt1,needt2,needt3,needt4
28      integer size_x3,size_x4
29      integer k_x3_offset,k_x4_offset
30      integer d_rx1,d_rx2,d_rx3,d_rx4
31      integer k_residual,l_residual
32      integer k_omegax, l_omegax
33      double precision au2ev   ! Conversion factor from a.u. to eV
34      parameter (au2ev=27.2113961d0)
35c
36      character*255 filename
37      character*255 modelname
38      logical converged
39c
40      nodezero=(ga_nodeid().eq.0)
41      needt1=.true.
42      needt2=.true.
43      needt3=.false.
44      needt4=.false.
45c
46      call tce_eom_xguess(needt1,needt2,needt3,needt4,
47     1     size_x1,size_x2,size_x3,size_x4,
48     2     k_x1_offset,k_x2_offset,k_x3_offset,k_x4_offset)
49      if (nxtrials.eq.0) goto 200
50      modelname = "EOM-CCSD right-hand side"
51      if (nodezero) write(LuOut,9220)
52     1   modelname(1:inp_strlen(modelname))
53      do iter=1,maxiter
54         if (nodezero.and.util_print('eom',print_default))
55     1      write(LuOut,9210) iter,nxtrials
56         do ivec = 1,nxtrials
57            if (.not.xp1_exist(ivec)) then
58               call tce_filenameindexed(ivec,'xp1',filename)
59               call createfile(filename,xp1(ivec),size_x1)
60               xp1_exist(ivec) = .true.
61               call dratoga(x1(ivec))
62               call dratoga(x2(ivec))
63               call eomccsd_x1_grad(d_f1,xp1(ivec),d_t1,d_t2,d_v2,
64     1            x1(ivec),x2(ivec),k_f1_offset,k_x1_offset,
65     2            k_t1_offset,k_t2_offset,k_v2_offset,
66     3            k_x1_offset,k_x2_offset)
67               call reconcilefile(xp1(ivec),size_x1)
68               call gatodra(x2(ivec))
69               call gatodra(x1(ivec))
70               call gatodra(xp1(ivec))
71            endif
72            if (.not.xp2_exist(ivec)) then
73               call tce_filenameindexed(ivec,'xp2',filename)
74               call createfile(filename,xp2(ivec),size_x2)
75               xp2_exist(ivec) = .true.
76               call dratoga(x1(ivec))
77               call dratoga(x2(ivec))
78               call eomccsd_x2_grad(d_f1,xp2(ivec),d_t1,d_t2,d_v2,
79     1              x1(ivec),x2(ivec),k_f1_offset,k_x2_offset,
80     2              k_t1_offset,k_t2_offset,k_v2_offset,
81     3              k_x1_offset,k_x2_offset)
82               call reconcilefile(xp2(ivec),size_x2)
83               call gatodra(x2(ivec))
84               call gatodra(x1(ivec))
85               call gatodra(xp2(ivec))
86            endif
87         enddo
88         if (.not.ma_push_get(mt_dbl,nxtrials,'residual',
89     1       l_residual,k_residual))
90     2       call errquit('eomccsd_energy: ma problem',10,ma_err)
91         call tce_eom_xdiagon_grad(needt1,needt2,needt3,needt4,
92     1        size_x1,size_x2,size_x3,size_x4,
93     2        k_x1_offset,k_x2_offset,k_x3_offset,k_x4_offset,
94     3        d_rx1,d_rx2,d_rx3,d_rx4,
95     4        dbl_mb(k_omegax),dbl_mb(k_residual))
96         cpu=cpu+util_cpusec()
97         wall=wall+util_wallsec()
98         converged = .true.
99         do ivec = 1,nroots_reduced
100            if (nodezero.and.(ivec.ne.nroots_reduced))
101     1         write(LuOut,9230) dbl_mb(k_residual+ivec-1),
102     2         dbl_mb(k_omegax+ivec-1),
103     3         dbl_mb(k_omegax+ivec-1)*au2ev
104            if (nodezero.and.(ivec.eq.nroots_reduced))
105     1         write(LuOut,9230) dbl_mb(k_residual+ivec-1),
106     2         dbl_mb(k_omegax+ivec-1),
107     3         dbl_mb(k_omegax+ivec-1)*au2ev,cpu,wall
108            if (nodezero) call util_flush(LuOut)
109            if (dbl_mb(k_residual+ivec-1).gt.thresh)
110     1         converged = .false.
111         enddo
112         cpu=-util_cpusec()
113         wall=-util_wallsec()
114         if (.not.ma_pop_stack(l_residual))
115     1      call errquit("eomccsd_energy: ma problem",20,ma_err)
116         if (converged) then
117            call tce_eom_xtidy
118            if (nodezero) then
119               write(LuOut,9240)
120               call util_flush(LuOut)
121            endif
122            goto 200
123         endif
124      enddo
125      call errquit('tce_energy: maxiter exceeded',iter,CALC_ERR)
126  200 continue
127 9210 format(/,1x,'Iteration ',i3,' using ',i4,' trial vectors')
128 9220 format(/,1x,A,' iterations',/,1x,
129     1'--------------------------------------------------------------'
130     2,/,1x,
131     3'     Residuum       Omega / hartree  Omega / eV    Cpu    Wall'
132     4,/,1x,
133     5'--------------------------------------------------------------')
134 9230 format(1x,f17.13,f18.13,f11.5,2f8.1)
135 9240 format(1x,
136     1'--------------------------------------------------------------'
137     2,/,1x,'Iterations converged')
138      end
139c $Id$
140