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