C C $Id: pme_derecmc.F,v 1.2 1998/07/16 16:40:16 jjv5 Exp $ C C------------------------------------------------------------------------ subroutine pme_derecmc(ir, igr, n, ch, derec) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" C this routine calculates the intramolecular contribution to C the derivative of the reciprocal energy to the coordinates C or residue ir and subtracts this term from derec for all C residues in ir's group. 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 *) ir: residue number C *) igr: group number C *) n: offset for groupmembers C *) derec: the derivative of the energy to the coordinates (kcal/A) C *) ch: charges 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 C note that the convolution Conv(theta,Q) is stored in QPMEC_ C parameters: dimension derec(maxpar), ch(maxatm) C locals: logical ly1,ly2,lx1,lx2,lx,lyx1,lyx2 dimension xmn(maxspline2p1), ymn(maxspline2p1), zmn(maxspline2p1), & dxmn(maxsplinep1), dymn(maxsplinep1), dzmn(maxsplinep1), & gderec(maxatres3) #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 irr1 = irpnt(ir) nat = irpnt(ir+1)-irr1 nat3 = 3*nat do 5 i=1,nat3 gderec(i) = 0.0 5 continue f1 = 332.0704*dk1pme*recip1(1) f2 = 332.0704*dk2pme*recip2(2) f3 = 332.0704*dk3pme*recip3(3) do 210 i=1,nat iat = irr1+i-1 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 = ch(iat)*zmn(iz) dz0 = ch(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 gderec(i1) = gderec(i1) + QPMEC_(iw)*y*dxmn(ix) gderec(i2) = gderec(i2) + QPMEC_(iw)*dy*xmn(ix) gderec(i3) = gderec(i3) + QPMEC_(iw)*dz*xmn(ix) 40 continue do 50 kx=kxb2,k1pmem1 iw = 2*(iwy + kx) - 1 ix = ixpk1pmem1+kx gderec(i1) = gderec(i1) + QPMEC_(iw)*y*dxmn(ix) gderec(i2) = gderec(i2) + QPMEC_(iw)*dy*xmn(ix) gderec(i3) = gderec(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 = ch(iat)*zmn(iz) dz0 = ch(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 gderec(i1) = gderec(i1) + QPMEC_(iw)*y*dxmn(ix) gderec(i2) = gderec(i2) + QPMEC_(iw)*dy*xmn(ix) gderec(i3) = gderec(i3) + QPMEC_(iw)*dz*xmn(ix) 80 continue do 90 kx=kxb2,k1pmem1 iw = 2*(iwy + kx) - 1 ix = ixpk1pmem1+kx gderec(i1) = gderec(i1) + QPMEC_(iw)*y*dxmn(ix) gderec(i2) = gderec(i2) + QPMEC_(iw)*dy*xmn(ix) gderec(i3) = gderec(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 = ch(iat)*zmn(iz) dz0 = ch(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 gderec(i1) = gderec(i1) + QPMEC_(iw)*y*dxmn(ix) gderec(i2) = gderec(i2) + QPMEC_(iw)*dy*xmn(ix) gderec(i3) = gderec(i3) + QPMEC_(iw)*dz*xmn(ix) 120 continue do 130 kx=kxb2,k1pmem1 iw = 2*(iwy + kx) - 1 ix = ixpk1pmem1+kx gderec(i1) = gderec(i1) + QPMEC_(iw)*y*dxmn(ix) gderec(i2) = gderec(i2) + QPMEC_(iw)*dy*xmn(ix) gderec(i3) = gderec(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 = ch(iat)*zmn(iz) dz0 = ch(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 gderec(i1) = gderec(i1) + QPMEC_(iw)*y*dxmn(ix) gderec(i2) = gderec(i2) + QPMEC_(iw)*dy*xmn(ix) gderec(i3) = gderec(i3) + QPMEC_(iw)*dz*xmn(ix) 160 continue do 170 kx=kxb2,k1pmem1 iw = 2*(iwy + kx)-1 ix = ixpk1pmem1+kx gderec(i1) = gderec(i1) + QPMEC_(iw)*y*dxmn(ix) gderec(i2) = gderec(i2) + QPMEC_(iw)*dy*xmn(ix) gderec(i3) = gderec(i3) + QPMEC_(iw)*dz*xmn(ix) 170 continue 180 continue 190 continue endif gderec(i1) = f1*gderec(i1) gderec(i2) = f2*gderec(i2) gderec(i3) = f3*gderec(i3) 210 continue C now subtract this intramolecular contribution from the gradient do i=0,ingroupn(igr)-1 irr = ingroup(n+i) i1 = irpnt(irr) i2 = irpnt(irr+1)-1 i13 = 3*i1-2 i23 = 3*i2-2 k = 1 do j=i13,i23 derec(j) = derec(j) - gderec(k) k = k + 1 enddo enddo end #undef QPMEC_ #undef QPME_