C C $Id: pme_calcq.F,v 1.7 1998/07/16 16:40:13 jjv5 Exp $ C C------------------------------------------------------------------------ subroutine pme_calcq implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" C this routine calculates Q as defined in Eq.4.6 of C J.Chem.Phys. 103 ('95) 8577. C C C compiler directives: C *) MEMORY_OVERLAP if defined, the pme variables will C . share memory with the eigenvectors / C . eigenvalues C C QPME_, a wrapped-around real array, will be filled C . with the values of Q. C C Written by Arjan van der Vaart, Dec. '97 C C C Note that the performance of this algorithm is very high C by demanding that nspline <= K_i, i=1,2,3 and by C shifting the scaled fractional coordinates upme such that C 0 <= upme_i <= K_i, i=1,2,3. C C Q needs to be recalculated whenever the coordinates or C charges change. C C logical ly1,ly2,lx1,lx2,lx,lyx1,lyx2 dimension xmn(maxspline2p1), & ymn(maxspline2p1), & zmn(maxspline2p1) #ifdef MEMORY_OVERLAP C now the following arrays share memory: C array1 common | array2 common | shared memory array C ----------------------------------------------------------------------- C qpmec(maxkpme32) /pme/ | ff(msorb2) /work/ | ff(maxovlp1) C qpme(maxkpme3) /pme/ | ww(msorb2) /work/ | ww(maxovlp2) C #define QPMEC_ ff #define QPME_ ww #else C no sharing of memory of the pme arrays with the /work/ arrays #define QPMEC_ qpmec #define QPME_ qpme #endif CC----------------------------------------------------------CC C initialize do 20 i=1,k123pme QPME_(i) = 0.0 20 continue do 200 iat=1,natoms C . calculate scaled fractional coordinates x = 0.0 y = 0.0 z = 0.0 do 30 j=1,3 x = x + recip1(j)*(xyz(j,iat)-xyzmin(j)) y = y + recip2(j)*(xyz(j,iat)-xyzmin(j)) z = z + recip3(j)*(xyz(j,iat)-xyzmin(j)) 30 continue upmex = dk1pme*x iupmex = upmex-nspline upmey = dk2pme*y iupmey = upmey-nspline upmez = dk3pme*z iupmez = upmez-nspline C . calculate the cardinal B-Splines Mn(upme-k). Only one call C . is needed to fill these in the x,y and z-direction for C . all the values of k for which Mn(upme-k) is NOT zero. call mbspline1(upmex,upmey,upmez,xmn,ymn,zmn) C . calculate the k values for which Mn(u-k-nK) is NOT zero. C . n = 0 for k=kxb1,kxe1 (same in y and z direction) C . n = -1 for k=kxb2,kxe2 (same in y and z direction) C . These are the only possible values of n (since 0 <= upme <= K C . and nspline <= K). kxe1 = upmex kxb2 = kxe1+k1pmenspl kxe1 = min(kxe1,k1pmem1) kxb1 = max(izero,iupmex) ixp = nspline-upmex ixp1 = ixp+1 ixpk1pmem1 = ixp-k1pmem1 kye1 = upmey kyb2 = kye1+k2pmenspl kye1 = min(kye1,k2pmem1) kyb1 = max(izero,iupmey) iyp = nspline-upmey kze1 = upmez kzb2 = kze1+k3pmenspl kze1 = min(kze1,k3pmem1) kzb1 = max(izero,iupmez) izp = nspline-upmez C . computation time can be safed if the program doesn't C . have to go in triple loops of which the inner loops C . are empty lx1 = kxb1.le.kxe1 lx2 = kxb2.le.k1pmem1 lx = lx1.or.lx2 ly1 = kyb1.le.kye1 ly2 = kyb2.le.k2pmem1 lyx1 = ly1.and.lx lyx2 = ly2.and.lx C . in the following the summation Q = q_i*Mn_x*Mn_y*Mn_z C . is performed for Mn_i .ne. 0 (i=x,y,z) if (lyx1) then do 70 kz=kzb1,kze1 iz = izp+1+kz z = chgpme(iat)*zmn(iz) iwz = kz*k1pmek2pme+1 do 60 ky=kyb1,kye1 iy = iyp+1+ky y = z*ymn(iy) iwy = iwz + ky*k1pme do 40 kx=kxb1,kxe1 iw = iwy + kx ix = ixp1+kx QPME_(iw) = QPME_(iw) + y*xmn(ix) 40 continue do 50 kx=kxb2,k1pmem1 iw = iwy + kx ix = ixpk1pmem1+kx QPME_(iw) = QPME_(iw) + y*xmn(ix) 50 continue 60 continue 70 continue endif if (lyx2) then do 110 kz=kzb1,kze1 iz = izp+1+kz z = chgpme(iat)*zmn(iz) iwz = kz*k1pmek2pme+1 do 100 ky=kyb2,k2pmem1 iy = iyp+ky-k2pmem1 y = z*ymn(iy) iwy = iwz + ky*k1pme do 80 kx=kxb1,kxe1 iw = iwy + kx ix = ixp1+kx QPME_(iw) = QPME_(iw) + y*xmn(ix) 80 continue do 90 kx=kxb2,k1pmem1 iw = iwy + kx ix = ixpk1pmem1+kx QPME_(iw) = QPME_(iw) + y*xmn(ix) 90 continue 100 continue 110 continue endif if (lyx1) then do 150 kz=kzb2,k3pmem1 iz = izp+kz-k3pmem1 z = chgpme(iat)*zmn(iz) iwz = kz*k1pmek2pme+1 do 140 ky=kyb1,kye1 iy = iyp+1+ky y = z*ymn(iy) iwy = iwz + ky*k1pme do 120 kx=kxb1,kxe1 iw = iwy + kx ix = ixp1+kx QPME_(iw) = QPME_(iw) + y*xmn(ix) 120 continue do 130 kx=kxb2,k1pmem1 iw = iwy + kx ix = ixpk1pmem1+kx QPME_(iw) = QPME_(iw) + y*xmn(ix) 130 continue 140 continue 150 continue endif if (lyx2) then do 190 kz=kzb2,k3pmem1 iz = izp+kz-k3pmem1 z = chgpme(iat)*zmn(iz) iwz = kz*k1pmek2pme+1 do 180 ky=kyb2,k2pmem1 iy = iyp+ky-k2pmem1 y = z*ymn(iy) iwy = iwz + ky*k1pme do 160 kx=kxb1,kxe1 iw = iwy + kx ix = ixp1+kx QPME_(iw) = QPME_(iw) + y*xmn(ix) 160 continue do 170 kx=kxb2,k1pmem1 iw = iwy + kx ix = ixpk1pmem1+kx QPME_(iw) = QPME_(iw) + y*xmn(ix) 170 continue 180 continue 190 continue endif 200 continue end CC----------------------------------------------------------------CC CC----------------------------------------------------------------CC #undef QPMEC_ #undef QPME_