C C $Id: pme_derec.F,v 1.2 1998/07/16 16:40:15 jjv5 Exp $ C C------------------------------------------------------------------------ subroutine pme_derec(derec) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" C this routine calculates the derivative of the reciprocal energy C to the coordinates. C C in order to save memory, dQ(k1,k2,k3)/dr is not calculated in C pme_calcq, but calculated here where Conv(theta,Q) is known. C C C Written by Arjan van der Vaart, Dec. '97 C C parameters: C C *) derec, a wrapped-around real array, will be filled with C . the derivative dErec/dr 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 derec needs to be recalculated whenever the coordinates or C charges change. C C note that the convolution Conv(theta,Q) is stored in QPMEC_ dimension derec(maxpar) logical ly1,ly2,lx1,lx2,lx,lyx1,lyx2 dimension xmn(maxspline2p1), ymn(maxspline2p1), zmn(maxspline2p1), & dxmn(maxsplinep1), dymn(maxsplinep1), dzmn(maxsplinep1) #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 f1 = 332.0704*dk1pme*recip1(1) f2 = 332.0704*dk2pme*recip2(2) f3 = 332.0704*dk3pme*recip3(3) do 210 iat=1,natoms i3 = 3*iat i2 = i3-1 i1 = i2-1 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) and derivatives. C . Only one call is needed to fill these in the x,y and z-direction C . for all the values of k for which Mn(upme-k) is NOT zero. call mbspline(upmex,upmey,upmez,xmn,ymn,zmn,dxmn,dymn,dzmn) 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 derivatives dQ(kx,ky,kz)/di, i=x,y,z C . are 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) dz0 = chgpme(iat)*dzmn(iz) iwz = kz*k1pmek2pme+1 do 60 ky=kyb1,kye1 iy = iyp+1+ky y = z*ymn(iy) dy = z*dymn(iy) dz = dz0*ymn(iy) iwy = iwz + ky*k1pme do 40 kx=kxb1,kxe1 iw = 2*(iwy + kx) - 1 ix = ixp1+kx derec(i1) = derec(i1) + QPMEC_(iw)*y*dxmn(ix) derec(i2) = derec(i2) + QPMEC_(iw)*dy*xmn(ix) derec(i3) = derec(i3) + QPMEC_(iw)*dz*xmn(ix) 40 continue do 50 kx=kxb2,k1pmem1 iw = 2*(iwy + kx) - 1 ix = ixpk1pmem1+kx derec(i1) = derec(i1) + QPMEC_(iw)*y*dxmn(ix) derec(i2) = derec(i2) + QPMEC_(iw)*dy*xmn(ix) derec(i3) = derec(i3) + QPMEC_(iw)*dz*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) dz0 = chgpme(iat)*dzmn(iz) iwz = kz*k1pmek2pme+1 do 100 ky=kyb2,k2pmem1 iy = iyp+ky-k2pmem1 y = z*ymn(iy) dy = z*dymn(iy) dz = dz0*ymn(iy) iwy = iwz + ky*k1pme do 80 kx=kxb1,kxe1 iw = 2*(iwy + kx) - 1 ix = ixp1+kx derec(i1) = derec(i1) + QPMEC_(iw)*y*dxmn(ix) derec(i2) = derec(i2) + QPMEC_(iw)*dy*xmn(ix) derec(i3) = derec(i3) + QPMEC_(iw)*dz*xmn(ix) 80 continue do 90 kx=kxb2,k1pmem1 iw = 2*(iwy + kx) - 1 ix = ixpk1pmem1+kx derec(i1) = derec(i1) + QPMEC_(iw)*y*dxmn(ix) derec(i2) = derec(i2) + QPMEC_(iw)*dy*xmn(ix) derec(i3) = derec(i3) + QPMEC_(iw)*dz*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) dz0 = chgpme(iat)*dzmn(iz) iwz = kz*k1pmek2pme+1 do 140 ky=kyb1,kye1 iy = iyp+1+ky y = z*ymn(iy) dy = z*dymn(iy) dz = dz0*ymn(iy) iwy = iwz + ky*k1pme do 120 kx=kxb1,kxe1 iw = 2*(iwy + kx) - 1 ix = ixp1+kx derec(i1) = derec(i1) + QPMEC_(iw)*y*dxmn(ix) derec(i2) = derec(i2) + QPMEC_(iw)*dy*xmn(ix) derec(i3) = derec(i3) + QPMEC_(iw)*dz*xmn(ix) 120 continue do 130 kx=kxb2,k1pmem1 iw = 2*(iwy + kx) - 1 ix = ixpk1pmem1+kx derec(i1) = derec(i1) + QPMEC_(iw)*y*dxmn(ix) derec(i2) = derec(i2) + QPMEC_(iw)*dy*xmn(ix) derec(i3) = derec(i3) + QPMEC_(iw)*dz*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) dz0 = chgpme(iat)*dzmn(iz) iwz = kz*k1pmek2pme+1 do 180 ky=kyb2,k2pmem1 iy = iyp+ky-k2pmem1 y = z*ymn(iy) dy = z*dymn(iy) dz = dz0*ymn(iy) iwy = iwz + ky*k1pme do 160 kx=kxb1,kxe1 iw = 2*(iwy + kx) - 1 ix = ixp1+kx derec(i1) = derec(i1) + QPMEC_(iw)*y*dxmn(ix) derec(i2) = derec(i2) + QPMEC_(iw)*dy*xmn(ix) derec(i3) = derec(i3) + QPMEC_(iw)*dz*xmn(ix) 160 continue do 170 kx=kxb2,k1pmem1 iw = 2*(iwy + kx)-1 ix = ixpk1pmem1+kx derec(i1) = derec(i1) + QPMEC_(iw)*y*dxmn(ix) derec(i2) = derec(i2) + QPMEC_(iw)*dy*xmn(ix) derec(i3) = derec(i3) + QPMEC_(iw)*dz*xmn(ix) 170 continue 180 continue 190 continue endif derec(i1) = f1*derec(i1) derec(i2) = f2*derec(i2) derec(i3) = f3*derec(i3) 210 continue end #undef QPMEC_ #undef QPME_