!----------------------------------------------------------------------- ! ! SUBROUTINE force(p,nterms,iterms,pos_com,dx,dy,dz,drdotdr,pmass,pquad,acc_g,option) ! ! !----------------------------------------------------------------------- ! ! ! Subroutine to compute the force on the p particle ! !----------------------------------------------------------------------- USE fly_h implicit none ! Declaration of local variables. ! ------------------------------- INTEGER(KIND=4) :: p,i, nterms,nqterms,j,n,iqin INTEGER(KIND=4), DIMENSION (maxnterm) :: qindex, smindex INTEGER(KIND=4), DIMENSION (ndim) ::i_mesh,ip,csign REAL(KIND=8) :: drdeldrg,phsm,drsm,accsm REAL(KIND=8) :: sdrdotdr, epsp,epsi,rpinveff REAL(KIND=8) :: xx,yy,zz,wt1,wt2,wt3,wt4,wt5,wt6,wt7,wt8 REAL(KIND=8) :: wi1,wj1,wk1,wi2,wj2,wk2 REAL(KIND=8), DIMENSION (maxnterm) ::r3inveff,rinveff,qr5inv,phiquad REAL(KIND=8), DIMENSION (maxnterm) :: r2inveff,acci REAL(KIND=8), DIMENSION (1:2*ndim-1) ::quad_l REAL(KIND=8), DIMENSION (ndim) ::pos_l,r_abs,f,acc_local,pos_s REAL(KIND=8), DIMENSION (ndim) ::acc_local_ew, acc_local_dir,acc_local_quad CHARACTER(LEN=4) :: option INTEGER :: il_sh REAL(KIND=8), INTENT(inout) :: pos_comm(:),dx(:),dy(:),dz(:),drdotdr(:) REAL(KIND=8), INTENT(inout) :: pmass(:) INTEGER(KIND=4), INTENT(inout) :: iterms(:) REAL(KIND=8), INTENT(inout) ::pquad(:,:) REAL(KIND=8), INTENT(inout) :: acc_g(:) !======================================================================= acc_local=0. acc_local_ew=0. acc_local_dir=0. acc_local_quad=0. !CACHE IF(group_access.EQ.0) THEN !came from grouping pos_l(1:ndim)=pos(1:ndim,p) ELSE !NOT came from grouping pos_l=pos_comm ENDIF !came from grouping IF(nterms.GT.maxnterm) THEN CALL error('force error1: overflow') ENDIF !----------------------------------------------------------------------- ! Compute monopole contribution !----------------------------------------------------------------------- rpinveff=1./(0.+tiny) epsp=eps !----------------------------------------------------------------------- ! Loop over interaction list. !----------------------------------------------------------------------- DO 30 i=1,nterms IF(group_access.EQ.0) THEN !came from grouping il_sh=iterms(i) IF(il_sh.gt.nbodsmax) THEN il_sh=il_sh-nbodsmax pos_s(1:ndim)=pos_cell(1:ndim,il_sh) ELSE !il_sh.gt.nbodsmax pos_s(1:ndim)=pos(1:ndim,il_sh) ENDIF !il_sh.gt.nbodsmax dx(i)=pos_l(1)-pos_s(1) dy(i)=pos_l(2)-pos_s(2) dz(i)=pos_l(3)-pos_s(3) drdotdr(i)=dx(i)**2+dy(i)**2+dz(i)**2 ENDIF ! come from grouping sdrdotdr=SQRT(drdotdr(i)) rinveff(i)=1./(sdrdotdr+tiny) r3inveff(i)=rinveff(i)/(drdotdr(i)+tiny) epsi=eps drdeldrg=sdrdotdr*ninterp/(epsp+epsi) smindex(i)=drdeldrg smindex(i)=MIN(ninterp,smindex(i)) drsm=MIN(one,drdeldrg-smindex(i)) phsm=(1.-drsm)*phsmooth(smindex(i))+ & drsm*phsmooth(1+smindex(i)) accsm=(1.-drsm)*acsmooth(smindex(i))+ & drsm*acsmooth(1+smindex(i)) rinveff(i)=phsm*rinveff(i) r3inveff(i)=accsm*r3inveff(i) acci(i)=pmass(i)*r3inveff(i) !----------------------------------------------------------------------- ! Ewald boundary periodic conditions section !----------------------------------------------------------------------- xx=dx(i) yy=dy(i) zz=dz(i) r_abs(1)=abs(xx) r_abs(2)=abs(yy) r_abs(3)=abs(zz) if(r_abs(1).eq.0.0) then csign(1)=0 else csign(1)=xx/r_abs(1) endif if(r_abs(2).eq.0.0) then csign(2)=0. else csign(2)=yy/r_abs(2) endif if(r_abs(3).eq.0.0) then csign(3)=0. else csign(3)=zz/r_abs(3) endif i_mesh=r_abs*linv ip=i_mesh+1 wi1=1.-linv*(r_abs(1)-rk(i_mesh(1),1)) wi2=1.-linv*(rk(ip(1),1)-r_abs(1)) wj1=1.-linv*(r_abs(2)-rk(i_mesh(2),2)) wj2=1.-linv*(rk(ip(2),2)-r_abs(2)) wk1=1.-linv*(r_abs(3)-rk(i_mesh(3),3)) wk2=1.-linv*(rk(ip(3),3)-r_abs(3)) wt1=wi1*wj1*wk1 wt2=wi2*wj1*wk1 wt3=wi1*wj2*wk1 wt4=wi1*wj1*wk2 wt5=wi2*wj2*wk1 wt6=wi2*wj1*wk2 wt7=wi1*wj2*wk2 wt8=wi2*wj2*wk2 f(1:3)=wt1*fc(i_mesh(1),i_mesh(2),i_mesh(3),1:3) f(1:3)=f(1:3)+wt2*fc(ip(1),i_mesh(2),i_mesh(3),1:3) f(1:3)=f(1:3)+wt3*fc(i_mesh(1),ip(2),i_mesh(3),1:3) f(1:3)=f(1:3)+wt4*fc(i_mesh(1),i_mesh(2),ip(3),1:3) f(1:3)=f(1:3)+wt5*fc(ip(1),ip(2),i_mesh(3),1:3) f(1:3)=f(1:3)+wt6*fc(ip(1),i_mesh(2),ip(3),1:3) f(1:3)=f(1:3)+wt7*fc(i_mesh(1),ip(2),ip(3),1:3) f(1:3)=f(1:3)+wt8*fc(ip(1),ip(2),ip(3),1:3) f=csign*f*pmass(i) acc_local_ew=acc_local_ew+f !----------------------------------------------------------------------- ! END Ewald section !----------------------------------------------------------------------- 30 CONTINUE IF(option.NE.'pot ') THEN DO 50 i=1,nterms acc_local_dir(1)=acc_local_dir(1)-dx(i)*acci(i) acc_local_dir(2)=acc_local_dir(2)-dy(i)*acci(i) acc_local_dir(3)=acc_local_dir(3)-dz(i)*acci(i) 50 CONTINUE ENDIF !----------------------------------------------------------------------- ! If required, compute quadrupole contribution. !----------------------------------------------------------------------- IF(usequad) THEN !----------------------------------------------------------------------- ! Filter out bodies. ! Compute quadrupole interaction from cells. !----------------------------------------------------------------------- nqterms=0 DO i=1,nterms IF(iterms(i).GT.nbodsmax) THEN nqterms=nqterms+1 qindex(nqterms)=i ENDIF ENDDO IF(option.NE.'pot ') THEN DO 80 i=1,nqterms iqin=qindex(i) quad_l(1:2*ndim-1)=pquad(1:2*ndim-1,qindex(i)) r2inveff(i)=rinveff(iqin)*rinveff(iqin) qr5inv(i)=r3inveff(iqin)*r2inveff(i) phiquad(i)=((-.5*((dx(iqin)**2-dz(iqin)**2)* & quad_l(1)+(dy(iqin)**2- & dz(iqin)**2)*quad_l(4))- & (dx(iqin)*dy(iqin)*quad_l(2)+ & dx(iqin)*dz(iqin)*quad_l(3)+ & dy(iqin)*dz(iqin)*quad_l(5)))* & qr5inv(i))*5.*r2inveff(i) acc_local_quad(1)=acc_local_quad(1)+dx(iqin)*phiquad(i)+ & (dx(iqin)*quad_l(1)+ & dy(iqin)*quad_l(2)+ & dz(iqin)*quad_l(3))*qr5inv(i) acc_local_quad(2)=acc_local_quad(2)+dy(iqin)*phiquad(i)+ & (dy(iqin)*quad_l(4)+ & dx(iqin)*quad_l(2)+ & dz(iqin)*quad_l(5))*qr5inv(i) acc_local_quad(3)=acc_local_quad(3)+dz(iqin)*phiquad(i)+ & (dz(iqin)*(-quad_l(1)- & quad_l(4))+dx(iqin)* & quad_l(3)+dy(iqin)* & quad_l(5))*qr5inv(i) 80 CONTINUE ENDIF ENDIF acc_local=acc_local_ew+acc_local_dir+acc_local_quad IF(group_access.EQ.0) THEN !access from grouping section acc_local(1)=acc_local(1)+acc_g(1) acc_local(2)=acc_local(2)+acc_g(2) acc_local(3)=acc_local(3)+acc_g(3) ENDIF IF(rmt_acc) THEN acc_rmt(1:ndim,p)=acc_rmt(1:ndim,p)+acc_local(1:ndim) ELSE acc(1:ndim,p)=acc_local(1:ndim) ENDIF RETURN END