1      subroutine matvecReal(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   ==================================================================
8      implicit none
9#include "errquit.fh"
10#include "inp.fh"
11#include "rtdb.fh"
12#include "stdio.fh"
13#include "nwc_const.fh"
14#include "mafdecls.fh"
15#include "global.fh"
16#include "testutil.fh"
17#include "dimqm_constants.fh"
18#include "dimqm.fh"
19c
20c     Input Variables
21      integer rtdb
22      integer nlen
23      integer g_vec
24      integer g_rhs
25      logical ldbg
26      double precision xyz(3,nDIM)
27c
28c     Variables from dimqm common blocks
29c      integer nDIM, nDIMTypes, nOrder, g_DIM_diag
30c
31c     Local Variables
32      integer ntype, mtype
33      integer m, n, ms, me, ns, ne, id
34      integer mLo, mUp, mOnNode, nodes
35      double precision time
36      double precision sPol(nDIMTypes)
37      double precision vec(nOrder)
38      integer tx(nDIM) ! Type Index
39      double precision diagS(nOrder)
40      double precision fld(3)
41      logical stat
42      integer lo, hi
43      call ga_sync()
44c
45c      Determine node ID
46      id = ga_nodeid()
47      if(ldbg .and. id.eq.0) write(LuOut,*)
48     $   "Start MatvecReal routine"
49      time = util_timer()
50c
51c      Zero right-hand side
52      call ga_zero(g_rhs)
53c
54c      Determine chunk to be done on this node
55      call pphilo(id, nDIM, mLo, mUp, mOnNode)
56c
57c      Pull values from RTDB and GAs
58      stat = rtdb_get(rtdb, 'dimpar:sPol', mt_dbl, nDIMTypes, sPol)
59      stat = rtdb_get(rtdb, 'dimpar:typeindex', mt_int, nDIM, tx)
60      call ga_get(g_DIM_diag, 1, nOrder, 1, 1, diagS, 1)
61      call ga_get(g_vec, 1, nOrder, 1, 1, vec, 1)
62      call ga_init_fence()
63c
64c      Loop over atoms on this node
65      do m = mLo, mUp
66        mtype = tx(m)
67        fld = ZERO
68        ms = 3 * (m - 1) + 1
69        me = 3 * (m - 1) + 3
70        do n = 1, nDIM
71          ntype = tx(n)
72          ns = 3 * (n - 1) + 1
73          ne = 3 * (n - 1) + 3
74          fld(1:3) = fld(1:3) + atomPairField(m, n,
75     $                                     xyz(:,m),     ! Coordinates to atom m
76     $                                     xyz(:,n),     ! Coordinates to atom n
77     $                                     sPol(mtype),  ! Static polarizability of atom m
78     $                                     sPol(ntype),  ! Static polarizability of atom n
79     $                                     diagS(ns:ne), ! Diagonal elements
80     $                                     vec(ns:ne),   ! Vector x in A*x = b
81     $                                     ldbg)         ! debug
82        end do
83        call ga_acc(g_rhs, ms, me, 1, 1, fld, 1, ONE)
84      end do
85      call ga_fence()
86      time = util_timer() - time
87      if(ldbg .and. id.eq.0) write(LuOut,*)
88     $   "End MatvecReal Routine"
89
90      contains
91      function atomPairField(m, n, mxyz, nxyz, mpol, npol,
92     $                        diag, vec, debug) result (fld)
93      implicit none
94#include "stdio.fh"
95#include "dimqm_constants.fh"
96      integer, intent(in) :: m, n
97      double precision, intent(in) :: mxyz(3), nxyz(3) ! Coordinates of the atoms
98      double precision, intent(in) :: mpol, npol ! Static polarizability of the atoms
99      double precision, intent(in) :: diag(3) ! Diagonal elements
100      double precision, intent(in) :: vec(3)  ! Vector being multiplied
101      logical, intent(in) :: debug
102c
103c     Output
104      double precision fld(3)
105c
106c     Local variables
107      double precision r(3)
108      double precision t2(3,3)
109      double precision dist, invdist
110      interface
111        function t2r(r, dist, invdist, a1, a2, pol1, pol2)
112          integer, intent(in) :: a1, a2
113          double precision, intent(in) :: dist, invdist
114          double precision, intent(in) :: pol1, pol2
115          double precision t2r(3,3)
116          double precision, intent(in) :: r(3)
117        end function t2r
118      end interface
119c
120c   Diagonal element
121      if(m .eq. n) then
122        fld(1:3) = vec(1:3) * diag(1:3)
123      else
124C   Off diagonal
125        r(1:3) = nxyz(1:3) - mxyz(1:3)
126        dist = SQRT(DOT_PRODUCT(r, r))
127        invdist = ONE / dist
128        t2 = t2r(r, dist, invdist, m, n, mpol, npol)
129        fld(1:3) = - MATMUL(t2, vec)
130      end if
131c      write(LuOut,*) "Fld:", fld
132      end function atomPairField
133      end subroutine matvecReal
134
135