1      subroutine matvecComplex(rtdb, nlen, g_vec, g_rhs, xyz, ldbg)
2c   =====================================================================
3c   Performs the operation A*x = b, where x and b are real.
4c
5c   In  : rtdb      - Handle to the RTDB
6c         nlen      - Length of the work array
7c         g_vec     - Global array handle for the vector to be multiplied
8c         xyz       - Cartesian coordinates to DIM atoms
9c         ldebug    - Debug flag
10c   Out : g_rhs     - Right-hand side of matrix-vector multiply (solution)
11c   =====================================================================
12      implicit none
13#include "errquit.fh"
14#include "inp.fh"
15#include "rtdb.fh"
16#include "stdio.fh"
17#include "nwc_const.fh"
18#include "mafdecls.fh"
19#include "global.fh"
20#include "testutil.fh"
21#include "dimqm_constants.fh"
22#include "dimqm.fh"
23c
24c     Input Variables
25      integer rtdb
26      integer nlen
27      integer g_vec
28      integer g_rhs
29      logical ldbg
30      double precision xyz(3,nDIM)
31c
32c     Local Variables
33      integer ntype, mtype
34      integer m, n, ms, me, ns, ne, id
35      integer mLo, mUp, mOnNode, nodes
36      double precision time
37      double precision sPol(nDIMTypes)
38      double complex vec(nOrder)
39      integer tx(nDIM) ! Type Index
40      double complex diagFD(nOrder)
41      double complex fld(3)
42      logical stat
43      integer lo, hi
44c
45c     Common variables used from dimqm.fh
46c
47c     integer nDIM, nDIMTypes, nOrder, d_DIM_diag
48c
49c      Determine node ID
50      id = ga_nodeid()
51      if(id.eq.0 .and. ldbg) write(LuOut,*)
52     $   "Start MatvecComplex routine"
53      time = util_timer()
54c
55c      Initialize right hand side
56      call ga_zero(g_rhs)
57c
58c      Determine chunk to be done on this node
59      call pphilo(id, nDIM, mLo, mUp, mOnNode)
60c
61c      Fill arrays
62      stat = rtdb_get(rtdb, 'dimpar:sPol', mt_dbl, nDIMTypes, sPol)
63      stat = rtdb_get(rtdb, 'dimpar:typeindex', mt_int, nDIM, tx)
64      call ga_get(g_DIM_diag, 1, nOrder, 1, 1, diagFD, 1)
65      call ga_get(g_vec, 1, nOrder, 1, 1, vec, 1)
66      call ga_sync()
67
68      call ga_init_fence()
69c
70c      Loop over atoms on this node
71      do m = mLo, mUp
72        mtype = tx(m)
73        fld = ZERO_C
74        ms = 3 * (m - 1) + 1
75        me = 3 * (m - 1) + 3
76        do n = 1, nDIM
77          ntype = tx(n)
78          ns = 3 * (n - 1) + 1
79          ne = 3 * (n - 1) + 3
80          fld(1:3) = fld(1:3) + atomPairField(m, n,
81     $                                     xyz(:,m),      ! xyz1
82     $                                     xyz(:,n),      ! xyz2
83     $                                     sPol(mtype),   ! static pol1 (for screening)
84     $                                     sPol(ntype),   ! static pol2 (for screening)
85     $                                     diagFD(ns:ne), ! ainv
86     $                                     vec(ns:ne),    ! vec
87     $                                     ldbg)          ! debug
88        end do
89        call ga_acc(g_rhs, ms, me, 1, 1, fld, 1, ONE_C)
90      end do
91      call ga_fence()
92c
93c      Clean up arrays
94      time = util_timer() - time
95      if(id.eq.0.and.ldbg) write(LuOut,*)
96     $   "End MatvecComplex Routine"
97
98      contains
99      function atomPairField(m, n, mxyz, nxyz, mpol, npol,
100     $                        diag, vec, debug) result (fld)
101c
102c   ============================================
103c   Calculates the field at atom m due to atom n
104c   ============================================
105c
106
107      implicit none
108#include "stdio.fh"
109#include "dimqm_constants.fh"
110      integer, intent(in) :: m, n
111      double precision, intent(in) :: mxyz(3), nxyz(3)
112      double precision, intent(in) :: mpol, npol
113      double complex,   intent(in) :: diag(3)
114      double complex,   intent(in) :: vec(3)
115      logical,          intent(in) :: debug
116c
117c     Output
118      double complex fld(3)
119c
120c     Local variables
121      double precision r(3)
122      double precision t2(3,3)
123      double precision dist, invdist
124      interface
125        function t2r(r, dist, invdist, a1, a2, pol1, pol2)
126          integer, intent(in) :: a1, a2
127          double precision, intent(in) :: dist, invdist
128          double precision, intent(in) :: pol1, pol2
129          double precision t2r(3,3)
130          double precision, intent(in) :: r(3)
131        end function t2r
132      end interface
133c
134c   Diagonal element
135      if(m .eq. n) then
136        fld(1:3) = vec(1:3) * diag(1:3)
137      else
138c   Off diagonal
139        r(1:3) = nxyz(1:3) - mxyz(1:3)
140        dist = SQRT(DOT_PRODUCT(r, r))
141        invdist = ONE / dist
142        t2 = t2r(r, dist, invdist, m, n, mpol, npol)
143        fld(1:3) = - MATMUL(t2, vec)
144      end if
145      end function atomPairField
146      end subroutine matvecComplex
147
148