1c
2c
3      Subroutine dim_eval_fnl(ptnl,nq,qwght,xyz,qxyz,npert,
4     $                        ipm, imag, muind)
5c
6      implicit none
7c
8#include "stdio.fh"
9#include "rtdb.fh"
10#include "mafdecls.fh"
11#include "dimqm_constants.fh"
12#include "dimqm.fh"
13#include "errquit.fh"
14#include "global.fh"
15#include "tcgmsg.fh"
16c
17c   Input Variables
18      integer nq                       ! Number of quadrature points
19      integer npert
20      double precision ptnl(nq, npert) ! DIM potential at each quadrature point for each perturbation
21      double precision qwght(nq)       ! Quadrature weights
22      double precision xyz(3, nDIM)    ! DIM atom coordinates
23      double precision qxyz(3, nq)     ! Quadrature point coodinates
24      integer ipm
25      integer imag
26      double precision muind(3, nDIM, npert)  ! Array that holds induced dipoles.  Structured as (component, atom, perturbation)
27c
28c   Local Variables
29      integer jj, n, i, idir, id
30      double precision r(3), dist, dist3
31      double precision screen_mu, screen
32      double precision util_erf
33c      character*(1) direction(3)
34c      character*(1) dpm(2)
35c      character*(1) dri(2)
36c      data direction /'x', 'y', 'z'/
37c      data dpm /'+', '-'/
38c      data dri /'r', 'i'/
39c
40      id = ga_nodeid()
41c      write(luout,*) "id:", npert, ipm, imag, nq
42c
43c   Zero potential array
44      call dfill(nq*npert, 0.d0, ptnl, 1)
45c      muind = ZERO
46      do idir = 1, npert ! Loop over perturbations
47c
48c        if(id .eq.0) then
49c        if(npert > 1) then ! only 1 perturbation for ground state, 3 for response
50c          if(ipm .gt. 0) then ! FD Response (real and imaginary, +/-)
51c            if(.not.rtdb_get(dimqm_rtdb,'dimqm:muind_'//direction(idir)
52c     $                        //'_'//dri(imag)//dpm(ipm), mt_dbl,
53c     $                                                  3*nDIM, muind))
54c     $      call errquit('get perturbed +/- dipoles failed',1,RTDB_ERR)
55c          else ! Static response
56c            if(.not.rtdb_get(dimqm_rtdb,'dimqm:muind_'//direction(idir),
57c     $                                        mt_dbl,3*nDIM, muind))
58c     $      call errquit('get perturbed xyz dipoles failed',1,RTDB_ERR)
59c          end if
60c        else ! Ground state
61c          if(.not.rtdb_get(dimqm_rtdb,'dimqm:muind',
62c     $                                        mt_dbl,3*nDIM, muind))
63c     $        call errquit('get perturbed dipoles failed',1,RTDB_ERR)
64c        endif
65c        end if
66c
67c        write(luout,*) "muind", id, muind
68        do n = 1, nq ! Loop over quadrature points
69          do i = 1, nDIM ! Loop over DIM atoms
70c
71c           Distance between two points
72            r(:) = xyz(:,i) - qxyz(:,n)
73            dist = SQRT(DOT_PRODUCT(r, r))
74c
75c           Screening
76            screen_mu = ONE
77            select case (scrnType)
78            case(ERFSCRN)
79              screen = util_erf(scrnFactor * dist)
80              screen_mu = screen * screen
81            case(EXPSCRN)
82              screen = 1.0d0 - EXP(-dist*scrnFactor)
83              screen_mu = screen * screen * screen
84            case default
85              if (dist < 1.0d-12) dist = 1.0d-12
86            end select
87            dist3 = dist * dist * dist
88c
89c           Calculate potential
90            ptnl(n,idir) = ptnl(n,idir)
91     $                           - screen_mu*muind(1,i,idir)*r(1)/dist3
92            ptnl(n,idir) = ptnl(n,idir)
93     $                           - screen_mu*muind(2,i,idir)*r(2)/dist3
94            ptnl(n,idir) = ptnl(n,idir)
95     $                           - screen_mu*muind(3,i,idir)*r(3)/dist3
96          end do ! End loop over DIM atoms
97        end do ! End loop over quadratrue points
98c
99c       Combine with quadrature weights
100        do jj = 1, nq
101          ptnl(jj,idir) = -ptnl(jj,idir)*qwght(jj)
102        enddo
103c
104      end do ! End loop over perturbations
105c      write(luout,*) "End eval_fnl"
106      return
107      end subroutine dim_eval_fnl
108