! #include "copyright.h" #include "dprec.h" #include "assert.h" !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine do_pmesh_kspace here] subroutine do_pmesh_kspace( numatoms,crd,charge, & frc,prefac1,prefac2,prefac3,fftable,stack,istack & #ifdef PSANDER ,known_atm & #endif ) implicit none # include "flocntrl.h" # include "ew_pme_recip.h" # include "ew_unitcell.h" # include "ew_frc.h" # include "def_time.h" # include "box.h" #ifdef MPI # include "parallel.h" # include "ew_parallel.h" #ifdef MPI_DOUBLE_PRECISION #undef MPI_DOUBLE_PRECISION #endif # include "mpif.h" #ifdef CRAY_PVP #define MPI_DOUBLE_PRECISION MPI_REAL8 #endif #endif ! INPUT ! numatoms: number of atoms ! crd atomic coords ! charge atomic charges integer numatoms,num_ks_trial _REAL_ crd(3,numatoms),charge(numatoms) integer nmine,nderiv #ifdef PSANDER integer known_atm #endif ! OUTPUT ! eer: ewald reciprocal or k-space energy ! frc forces incremented by k-space sum _REAL_ frc(3,numatoms) ! HEAP STORAGE: These arrays need to be preserved throughout simulation _REAL_ prefac1(*),prefac2(*),prefac3(*),fftable(*) ! STACK STORAGE: These arrays can be tossed after leaving this routine _REAL_ d2th1,d2th2,d2th3 integer nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork,sfft,sffw _REAL_ stack(*) integer istack(*) integer l_fr1,l_fr2,l_fr3 integer l_th1,l_th2,l_th3,l_dth1,l_dth2,l_dth3 integer l_fftw,l_q integer imy_cg integer num_ks save num_ks data num_ks/0/ integer l_tmpy,l_alpha,l_beta nderiv = 1 ! FLOW CONTROL FLAG (debug, possible vacuum sims and future respa) if ( do_rec == 0 )return ! get some integer array dimensions call get_fftdims(nfft1,nfft2,nfft3, & nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork,sfft,sffw) # ifdef MPI # ifdef LES num_ks = numatoms # else if(num_ks == 0) & num_ks = min((((nxyslab(0)+order-1)*numatoms*4)/ & (3*nfft3)),numatoms) # endif # else num_ks = numatoms # endif call timer_start(TIME_BSPL) num_ks_trial = 0 100 call get_stack(l_fftw,sizffwrk) call get_stack(l_q,siz_q) call get_stack(l_fr1,num_ks) call get_stack(l_fr2,num_ks) call get_stack(l_fr3,num_ks) call get_stack(l_th1,num_ks*order) call get_stack(l_th2,num_ks*order) call get_stack(l_th3,num_ks*order) if ( nderiv == 1 )then call get_stack(l_dth1,num_ks*order) call get_stack(l_dth2,num_ks*order) call get_stack(l_dth3,num_ks*order) end if call get_istack(imy_cg,num_ks) call get_grid_weights( & numatoms,crd,recip,nfft1,nfft2,nfft3, & stack(l_fr1),stack(l_fr2),stack(l_fr3),order, & stack(l_th1),stack(l_th2),stack(l_th3), & stack(l_dth1),stack(l_dth2),stack(l_dth3), & d2th1,d2th2,d2th3, & #ifdef PSANDER known_atm, & #endif istack(imy_cg),nmine,nderiv,num_ks) if(nmine > num_ks)then #ifdef MPI write(6,'("***** Processor ",i6)') mytaskid #endif write(6,'("***** System must be very inhomogeneous.")') write(6,'("***** Readjusting recip sizes.")') write(6,'(A,i9,A,i9/)') & " In this slab, Atoms found: ",nmine, & " Allocated: ",num_ks if( num_ks_trial >= 2 ) call mexit( 6,1 ) if ( nderiv == 1 )then call free_stack(l_dth3) call free_stack(l_dth2) call free_stack(l_dth1) end if call free_stack(l_th3) call free_stack(l_th2) call free_stack(l_th1) call free_stack(l_fr3) call free_stack(l_fr2) call free_stack(l_fr1) call free_stack(l_q) call free_stack(l_fftw) call free_istack(imy_cg) num_ks= nmine*4/3 num_ks_trial = num_ks_trial + 1 goto 100 end if call timer_stop_start(TIME_BSPL,TIME_FILLG) !........Fill Charge Grid ! charges are approximated on an even grid call fill_charge_grid(numatoms,charge, & stack(l_th1),stack(l_th2),stack(l_th3), & stack(l_fr1),stack(l_fr2),stack(l_fr3), & order,nfft1,nfft2,nfft3,nfftdim1,nfftdim2 & # ifdef MPI ,mxyslabs & # else ,nfftdim3 & # endif ,stack(l_q),istack(imy_cg),nmine) call timer_stop_start(TIME_FILLG,TIME_FFT) call get_stack(l_tmpy,2*nfftdim1) call get_stack(l_alpha,nfft1) call get_stack(l_beta,nfft1) call fft_backrc( & stack(l_q),fftable,stack(l_fftw), & nfft1,nfft2,nfft3, & nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork & ,stack(l_tmpy), & stack(l_alpha),stack(l_beta) & ) call timer_stop_start(TIME_FFT,TIME_SCSUM) ! -------------SCALAR SUM------------------ if( ifbox == 1 ) then call scalar_sumrc_orthog( & stack(l_q), & ew_coeff,volume,recip,prefac1,prefac2,prefac3, & nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,eer,rec_vir) else call scalar_sumrc( & stack(l_q), & ew_coeff,volume,recip,prefac1,prefac2,prefac3, & nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,eer,rec_vir) end if call timer_stop_start(TIME_SCSUM,TIME_FFT) ! -----------FFT FORWARD-------------------- call fft_forwardrc( & stack(l_q),fftable,stack(l_fftw), & nfft1,nfft2,nfft3, & nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork & ,stack(l_tmpy), & stack(l_alpha),stack(l_beta) & ) call free_stack(l_beta) call free_stack(l_alpha) call free_stack(l_tmpy) call timer_stop_start(TIME_FFT,TIME_GRADS) ! -----------GRAD SUM-------------------- call grad_sumrc( & numatoms,charge,recip, & stack(l_th1),stack(l_th2),stack(l_th3), & stack(l_dth1),stack(l_dth2),stack(l_dth3), & frc,stack(l_fr1),stack(l_fr2),stack(l_fr3), & frcx,frcy,frcz, & order,nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3, & stack(l_q), & istack(imy_cg),nmine & ) call timer_stop(TIME_GRADS) if ( nderiv == 1 )then call free_stack(l_dth3) call free_stack(l_dth2) call free_stack(l_dth1) end if call free_stack(l_th3) call free_stack(l_th2) call free_stack(l_th1) call free_stack(l_fr3) call free_stack(l_fr2) call free_stack(l_fr1) call free_stack(l_q) call free_stack(l_fftw) call free_istack(imy_cg) #ifdef MPI call mpi_barrier(recip_comm,ierr) #endif return end subroutine do_pmesh_kspace !------------------------------------------------------------------- ! ----- ! --- FILL_CHARGE_GRID -- RC------- ! ----- !------------------------------------------------------------------- !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine fill_charge_grid here] subroutine fill_charge_grid(numatoms,charge, & theta1,theta2,theta3,fr1,fr2,fr3, & order,nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,q,my_cg,nmine) !--------------------------------------------------------------------- ! INPUT: ! numatoms: number of atoms ! charge: the array of atomic charges ! theta1,theta2,theta3: the spline coeff arrays ! fr1,fr2,fr3 the scaled and shifted fractional coords ! nfft1,nfft2,nfft3: the charge grid dimensions ! nfftdim1,nfftdim2,nfftdim3: physical charge grid dims ! order: the order of spline interpolation ! OUTPUT: ! Q the charge grid !--------------------------------------------------------------------- implicit none integer numatoms,order,nfft1,nfft2,nfft3 integer nfftdim1,nfftdim2,nfftdim3 _REAL_ fr1(numatoms),fr2(numatoms),fr3(numatoms) _REAL_ theta1(order,numatoms),theta2(order,numatoms), & theta3(order,numatoms),charge(numatoms) _REAL_ q(*) integer my_cg(*),nmine #ifdef MPI # include "parallel.h" # include "ew_parallel.h" # include "ew_bspline.h" integer kbot0 #endif integer n,ith1,ith2,ith3,i0,j0,k0,i,j,k,i00,j00,iqk,iqj _REAL_ prod integer im,kq,ntot !........Zero the Charge grids ntot = 2*nfftdim1*nfftdim2*nfftdim3 call zero_array(q,ntot) #ifdef MPI kbot0 = mxystart(mytaskid) kbot = kbot0 + 1 ktop = kbot0 + mxyslabs #endif !ocl novrec do im = 1,nmine #ifdef MPI n = my_cg(im) #else n = im #endif k0 = int(fr3(im)) - order j00 = int(fr2(im)) - order i00 = int(fr1(im)) - order !ocl novrec do ith3 = 1,order k0 = k0 + 1 if (k0 >= 0) then k = k0 + 1 else k = k0 + 1 + nfft3 end if #ifdef MPI if ( k >= kbot .and. k <= ktop ) then kq = k - kbot0 #else kq = k #endif iqk = (kq-1)*2*nfftdim1*nfftdim2 j0 = j00 !ocl novrec do ith2 = 1,order j0 = j0 + 1 iqj = iqk + j0*2*nfftdim1 if( j0 < 0 ) iqj = iqj + 2*nfft2*nfftdim1 prod = theta2(ith2,im)*theta3(ith3,im)*charge(n) i0 = i00 + 1 !ocl novrec do ith1 = 1,order i0 = i0 + 1 if (i0 >= 1) then q(i0+iqj) = q(i0+iqj) + theta1(ith1,im)*prod else q(i0+nfft1+iqj) = q(i0+nfft1+iqj) + theta1(ith1,im)*prod end if end do end do #ifdef MPI end if #endif end do ! ith3 = 1,order end do ! im = 1,nmine return end subroutine fill_charge_grid !------------------------------------------------------------------- ! GRAD_SUM RC **** REAL not complex **** !----------------------------------------------------------------- !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine grad_sumrc here] subroutine grad_sumrc( & numatoms,charge,recip,theta1,theta2,theta3, & dtheta1,dtheta2,dtheta3,frc,fr1,fr2,fr3, & frcx,frcy,frcz, & order,nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3, & q,my_cg,nmine & ) implicit none integer numatoms,order,nfft1,nfft2,nfft3 integer nfftdim1,nfftdim2,nfftdim3 integer kq #include "box.h" #ifdef MPI # include "ew_parallel.h" # include "parallel.h" # ifdef MPI_DOUBLE_PRECISION # undef MPI_DOUBLE_PRECISION # endif # include "mpif.h" # ifdef CRAY_PVP # define MPI_DOUBLE_PRECISION MPI_REAL8 # endif # include "ew_bspline.h" ! _REAL_ Q(nfftdim1*2,nfftdim2,mxyslabs) _REAL_ q(*) #else ! _REAL_ Q(nfftdim1*2,nfftdim2,nfft3) _REAL_ q(*) ! If MPI only some atoms are touched by the local process. Hence loop is over ! im=1,nmine, arrays are filled in order of being seen by process ! and n is obtained by a pointer. #endif _REAL_ recip(3,3),recip11,recip22,recip33 _REAL_ fr1(numatoms),fr2(numatoms),fr3(numatoms) _REAL_ frc(3,numatoms) _REAL_ theta1(order,numatoms),theta2(order,numatoms), & theta3(order,numatoms),charge(numatoms) _REAL_ dtheta1(order,numatoms),dtheta2(order,numatoms), & dtheta3(order,numatoms) integer my_cg(*),nmine integer iqk,iqj,j00,i00 integer n,ith1,ith2,ith3,i0,j0,k0,i,j,k,im _REAL_ f1,f2,f3,term,chargen _REAL_ frcx,frcy,frcz _REAL_ dfx,dfy,dfz,dnfft1,dnfft2,dnfft3 _REAL_ f1fac,f2fac,f3fac dnfft1 = nfft1 dnfft2 = nfft2 dnfft3 = nfft3 recip11 = recip(1,1)*dnfft1 recip22 = recip(2,2)*dnfft2 recip33 = recip(3,3)*dnfft3 do im = 1,nmine #ifdef MPI n = my_cg(im) #else n = im #endif f1 = 0.d0 f2 = 0.d0 f3 = 0.d0 i00 = int(fr1(im)) - order j00 = int(fr2(im)) - order k0 = int(fr3(im)) - order chargen = charge(n) do ith3 = 1,order k0 = k0 + 1 if (k0 >= 0) then k = k0 + 1 else k = k0 + 1 + nfft3 end if #ifdef MPI if ( k >= kbot .and. k <= ktop )then kq = k - mxystart(mytaskid) #else kq = k #endif iqk = (kq-1)*2*nfftdim1*nfftdim2 j0 = j00 do ith2 = 1,order j0 = j0 + 1 iqj = iqk + j0*2*nfftdim1 if( j0 < 0 ) iqj = iqj + 2*nfft2*nfftdim1 i0 = i00 + 1 f1fac = theta2(ith2,im) * theta3(ith3,im) * chargen f2fac = dtheta2(ith2,im) * theta3(ith3,im) * chargen f3fac = theta2(ith2,im) * dtheta3(ith3,im) * chargen do ith1 = 1,order i0 = i0 + 1 if (i0 >= 1) then term = q(i0 + iqj) else term = q(i0 + nfft1 + iqj) end if ! ---force is negative of grad f1 = f1 - term * dtheta1(ith1,im) * f1fac f2 = f2 - term * theta1(ith1,im) * f2fac f3 = f3 - term * theta1(ith1,im) * f3fac end do end do #ifdef MPI end if ! ( k >= kbot .and. k <= ktop ) #endif end do ! ith3 = 1,order if ( ifbox == 1 ) then ! orthogonal unit cell dfx = recip11 * f1 dfy = recip22 * f2 dfz = recip33 * f3 else f1 = dnfft1*f1 f2 = dnfft2*f2 f3 = dnfft3*f3 dfx=recip(1,1)*f1+recip(1,2)*f2+recip(1,3)*f3 dfy=recip(2,1)*f1+recip(2,2)*f2+recip(2,3)*f3 dfz=recip(3,1)*f1+recip(3,2)*f2+recip(3,3)*f3 end if frc(1,n) = frc(1,n) + dfx frc(2,n) = frc(2,n) + dfy frc(3,n) = frc(3,n) + dfz frcx=frcx+dfx frcy=frcy+dfy frcz=frcz+dfz end do ! im = 1,nmine return end subroutine grad_sumrc !================================================================= ! --- SCALAR_SUM RC--- !------------------------------------------------------------------- !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine scalar_sumrc here] subroutine scalar_sumrc( & q, & ewaldcof,volume,recip,prefac1,prefac2,prefac3, & nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,eer,vir) implicit none # include "extra.h" # include "box.h" # include "md.h" #ifdef MPI # include "ew_parallel.h" # include "parallel.h" integer ier,i #endif integer nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3 _REAL_ prefac1(nfft1),prefac2(nfft2), & prefac3(nfft3),ewaldcof,volume _REAL_ eer,vir(3,3) _REAL_ recip(3,3) integer k2q _REAL_ q(2,nfft3,nfftdim1,nfft2) _REAL_ pi,fac,denom,eterm,vterm,energy integer k1,k2,k3,m1,m2,m3,nff,indtop integer nf1,nf2,nf3 integer k10 _REAL_ mhat1,mhat2,mhat3,msq,struc2,msq_inv,piv_inv integer k1s,k2s,k3s,m1s,m2s,m3s _REAL_ mhat1s,mhat2s,mhat3s,msqs _REAL_ struc2s,eterms,vterms,denoms,tmp1,tmp2 indtop = nfft1*nfft2*nfft3 pi = 3.14159265358979323846d0 piv_inv = 1.d0/(pi*volume) fac = (pi*pi)/(ewaldcof*ewaldcof) nff = nfft1*nfft2 nf1 = nfft1/2 if ( 2*nf1 < nfft1 )nf1 = nf1+1 nf2 = nfft2/2 if ( 2*nf2 < nfft2 )nf2 = nf2+1 nf3 = nfft3/2 if ( 2*nf3 < nfft3 )nf3 = nf3+1 energy = 0.d0 k10 = 1 do m2 = 1,3 do m1 = 1,3 vir(m1,m2) = 0.d0 end do end do !........Insist that Q(1,1,1,1) is set to 0 (true already for neutral) if(master)then q(1,1,1,1) = 0.d0 q(2,1,1,1) = 0.d0 end if !====================================================================== ! BIG LOOP !====================================================================== #ifdef MPI do k2q = 1, mxzslabs if(master)then k2=k2q k2s=mod(nfft2-k2+1,nfft2)+1 else k2 = k2q + mxzstart(mytaskid) k2s=mod(nfft2-k2+1,nfft2)+1 end if #else do k2q = 1, nfft2 k2s=mod(nfft2-k2q+1,nfft2)+1 k2=k2q #endif m2 = k2 - 1 if ( k2 > nf2 )m2 = k2 - 1 - nfft2 do k3 = 1,nfft3 k3s=mod(nfft3-k3+1,nfft3)+1 m3 = k3 - 1 if ( k3 > nf3 )m3 = k3 - 1 - nfft3 k10 = 1 if(master)then if(k3+k2 == 2) k10 = 2 end if do k1 = k10, nf1+1 k1s=nfft1-k1+2 m1 = k1 - 1 if ( k1 > nf1 ) m1 = k1 - 1 - nfft1 mhat1 = recip(1,1)*m1+recip(1,2)*m2+recip(1,3)*m3 mhat2 = recip(2,1)*m1+recip(2,2)*m2+recip(2,3)*m3 mhat3 = recip(3,1)*m1+recip(3,2)*m2+recip(3,3)*m3 msq = mhat1*mhat1+mhat2*mhat2+mhat3*mhat3 msq_inv = 1.d0/msq eterm = exp(-fac*msq)*prefac1(k1)*prefac2(k2)*prefac3(k3) & *piv_inv*msq_inv struc2 = q(1,k3,k1,k2q)*q(1,k3,k1,k2q) + & q(2,k3,k1,k2q)*q(2,k3,k1,k2q) tmp1 = eterm*struc2 energy = energy + tmp1 #ifndef noVIRIAL vterm = 2.d0*(fac*msq + 1.d0)*msq_inv tmp2 = tmp1*vterm vir(1,1) = vir(1,1) + tmp1 * (vterm*mhat1*mhat1 - 1.d0) vir(1,2) = vir(1,2) + tmp2*mhat1*mhat2 vir(1,3) = vir(1,3) + tmp2*mhat1*mhat3 vir(2,1) = vir(2,1) + tmp2*mhat2*mhat1 vir(2,2) = vir(2,2) + tmp1 * (vterm*mhat2*mhat2 - 1.d0) vir(2,3) = vir(2,3) + tmp2*mhat2*mhat3 vir(3,1) = vir(3,1) + tmp2*mhat3*mhat1 vir(3,2) = vir(3,2) + tmp2*mhat3*mhat2 vir(3,3) = vir(3,3) + tmp1 * (vterm*mhat3*mhat3 - 1.d0) #endif if( k1 > 1 .and. k1 <= nfft1 )then m1s = k1s - 1 if ( k1s > nf1 )m1s = k1s - 1 - nfft1 m2s = k2s - 1 if ( k2s > nf2 )m2s = k2s - 1 - nfft2 m3s = k3s - 1 if ( k3s > nf3 )m3s = k3s - 1 - nfft3 mhat1s = recip(1,1)*m1s+recip(1,2)*m2s+recip(1,3)*m3s mhat2s = recip(2,1)*m1s+recip(2,2)*m2s+recip(2,3)*m3s mhat3s = recip(3,1)*m1s+recip(3,2)*m2s+recip(3,3)*m3s msqs = mhat1s*mhat1s+mhat2s*mhat2s+mhat3s*mhat3s msq_inv = 1.d0/msqs eterms = exp(-fac*msqs)*prefac1(k1s)*prefac2(k2s)* & prefac3(k3s)*piv_inv*msq_inv ! struc2s = Q(1,k3,k1,k2q)*Q(1,k3,k1,k2q) + ! $ Q(2,k3,k1,k2q)*Q(2,k3,k1,k2q) ! tmp1 = eterms*struc2s tmp1 = eterms*struc2 energy = energy + tmp1 #ifndef noVIRIAL vterms = 2.d0*(fac*msqs + 1.d0)*msq_inv tmp2 = tmp1*vterms vir(1,1) = vir(1,1) + tmp1*(vterms*mhat1s*mhat1s - 1.d0) vir(1,2) = vir(1,2) + tmp2*mhat1s*mhat2s vir(1,3) = vir(1,3) + tmp2*mhat1s*mhat3s vir(2,1) = vir(2,1) + tmp2*mhat2s*mhat1s vir(2,2) = vir(2,2) + tmp1*(vterms*mhat2s*mhat2s - 1.d0) vir(2,3) = vir(2,3) + tmp2*mhat2s*mhat3s vir(3,1) = vir(3,1) + tmp2*mhat3s*mhat1s vir(3,2) = vir(3,2) + tmp2*mhat3s*mhat2s vir(3,3) = vir(3,3) + tmp1*(vterms*mhat3s*mhat3s - 1.d0) #endif end if q(1,k3,k1,k2q) = eterm * q(1,k3,k1,k2q) q(2,k3,k1,k2q) = eterm * q(2,k3,k1,k2q) end do ! k1 = k10, nf1+1 end do ! k3 = 1,nfft3 end do ! k2q = 1, mxzslabs eer = 0.5d0 * energy #ifndef noVIRIAL do m2 = 1,3 do m1 = 1,3 vir(m1,m2) = 0.5d0*vir(m1,m2) end do end do #endif return end subroutine scalar_sumrc !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine scalar_sumrc_orthog here] subroutine scalar_sumrc_orthog( & q, & ewaldcof,volume,recip,prefac1,prefac2,prefac3, & nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,eer,vir) implicit none # include "extra.h" # include "box.h" # include "md.h" #ifdef MPI # include "ew_parallel.h" # include "parallel.h" #endif integer ier,i,itot integer nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3 _REAL_ prefac1(nfft1),prefac2(nfft2), & prefac3(nfft3),ewaldcof,volume _REAL_ eer,vir(3,3) _REAL_ recip(3,3) integer k2q _REAL_ q(2,nfft3,nfftdim1,nfft2) _REAL_ pi,fac,denom,eterm,vterm,energy integer k1,k2,k3,m1,m2,m3,nff,indtop integer nf1,nf2,nf3 integer k10 _REAL_ mhat1,mhat2,mhat3,msq,struc2,msq_inv,piv_inv integer k1s,k2s,k3s,m1s,m2s,m3s _REAL_ mhat1s,mhat2s,mhat3s,msqs _REAL_ struc2s,eterms,vterms,denoms _REAL_ tmp1,tmp2,m2_m3_tbl,m2_m3_tbls _REAL_ recip11,recip22,recip33 _REAL_ recip11sq,recip22sq,recip33sq _REAL_, allocatable, save :: m1_tbl(:) _REAL_, allocatable, save :: m2_tbl(:) _REAL_, allocatable, save :: m3_tbl(:) logical, save:: first data first /.true./ indtop = nfft1*nfft2*nfft3 pi = 3.14159265358979323846d0 piv_inv = 1.d0/(pi*volume) fac = (pi*pi)/(ewaldcof*ewaldcof) nff = nfft1*nfft2 nf1 = nfft1/2 if ( 2*nf1 < nfft1 )nf1 = nf1+1 nf2 = nfft2/2 if ( 2*nf2 < nfft2 )nf2 = nf2+1 nf3 = nfft3/2 if ( 2*nf3 < nfft3 )nf3 = nf3+1 recip11 = recip(1,1) recip22 = recip(2,2) recip33 = recip(3,3) if( first) then allocate(m1_tbl(-(nfft1/2 + 1) : nfft1/2 + 1), & m2_tbl(-(nfft2/2 + 1) : nfft2/2 + 1), & m3_tbl(-(nfft3/2 + 1) : nfft3/2 + 1), & stat = ier) REQUIRE( ier==0 ) end if if( ntp > 0 .or. first) then recip11sq = recip11*recip11 recip22sq = recip22*recip22 recip33sq = recip33*recip33 itot = 0 do i = -(nfft1/2 + 1), nfft1/2 + 1 itot = itot + 1 m1_tbl(i) = -fac * dble(i) * dble(i) * recip11sq end do call vdexp(itot,m1_tbl(-nfft1/2-1),m1_tbl(-nfft1/2-1)) itot = 0 do i = -(nfft2/2 + 1), nfft2/2 + 1 itot = itot + 1 m2_tbl(i) = -fac * dble(i) * dble(i) * recip22sq end do call vdexp(itot,m2_tbl(-nfft2/2-1),m2_tbl(-nfft2/2-1)) itot = 0 do i = -(nfft3/2 + 1), nfft3/2 + 1 itot = itot + 1 m3_tbl(i) = -fac * dble(i) * dble(i) * recip33sq end do call vdexp(itot,m3_tbl(-nfft3/2-1),m3_tbl(-nfft3/2-1)) first = .false. end if energy = 0.d0 k10 = 1 do m2 = 1,3 do m1 = 1,3 vir(m1,m2) = 0.d0 end do end do !........Insist that Q(1,1,1,1) is set to 0 (true already for neutral) if(master)then q(1,1,1,1) = 0.d0 q(2,1,1,1) = 0.d0 end if !====================================================================== ! BIG LOOP !====================================================================== #ifdef MPI do k2q = 1, mxzslabs if(master)then k2=k2q k2s=mod(nfft2-k2+1,nfft2)+1 else k2 = k2q + mxzstart(mytaskid) k2s=mod(nfft2-k2+1,nfft2)+1 end if #else do k2q = 1, nfft2 k2s=mod(nfft2-k2q+1,nfft2)+1 k2=k2q #endif m2 = k2 - 1 if ( k2 > nf2 )m2 = k2 - 1 - nfft2 mhat2 = recip22*m2 m2s = k2s - 1 if ( k2s > nf2 )m2s = k2s - 1 - nfft2 mhat2s = recip22*m2s do k3 = 1,nfft3 k3s=mod(nfft3-k3+1,nfft3)+1 m3 = k3 - 1 if ( k3 > nf3 )m3 = k3 - 1 - nfft3 mhat3 = recip33*m3 m3s = k3s - 1 if ( k3s > nf3 )m3s = k3s - 1 - nfft3 mhat3s = recip33*m3s k10 = 1 if(master)then if(k3+k2 == 2) k10 = 2 end if m2_m3_tbl = piv_inv*m2_tbl(m2)*m3_tbl(m3) & *prefac2(k2)*prefac3(k3) m2_m3_tbls = piv_inv*m2_tbl(m2s)*m3_tbl(m3s) & *prefac2(k2s)*prefac3(k3s) do k1 = k10, nf1+1 k1s=nfft1-k1+2 m1 = k1 - 1 if ( k1 > nf1 ) m1 = k1 - 1 - nfft1 mhat1 = recip11*m1 msq = mhat1*mhat1+mhat2*mhat2+mhat3*mhat3 msq_inv = 1.d0/msq eterm = m1_tbl(m1) * m2_m3_tbl * & prefac1(k1) * msq_inv struc2 = q(1,k3,k1,k2q)*q(1,k3,k1,k2q) + & q(2,k3,k1,k2q)*q(2,k3,k1,k2q) tmp1 = eterm*struc2 energy = energy + tmp1 #ifndef noVIRIAL vterm = 2.d0*(fac*msq + 1.d0)*msq_inv tmp2 = tmp1*vterm vir(1,1) = vir(1,1) + tmp1 * (vterm*mhat1*mhat1 - 1.d0) vir(1,2) = vir(1,2) + tmp2*mhat1*mhat2 vir(1,3) = vir(1,3) + tmp2*mhat1*mhat3 vir(2,1) = vir(2,1) + tmp2*mhat2*mhat1 vir(2,2) = vir(2,2) + tmp1 * (vterm*mhat2*mhat2 - 1.d0) vir(2,3) = vir(2,3) + tmp2*mhat2*mhat3 vir(3,1) = vir(3,1) + tmp2*mhat3*mhat1 vir(3,2) = vir(3,2) + tmp2*mhat3*mhat2 vir(3,3) = vir(3,3) + tmp1 * (vterm*mhat3*mhat3 - 1.d0) #endif if( k1 > 1 .and. k1 <= nfft1 )then m1s = k1s - 1 if ( k1s > nf1 )m1s = k1s - 1 - nfft1 mhat1s = recip11*m1s msqs = mhat1s*mhat1s+mhat2s*mhat2s+mhat3s*mhat3s msq_inv = 1.d0/msqs eterms = m1_tbl(m1s) * m2_m3_tbls * & prefac1(k1s) * msq_inv tmp1 = eterms*struc2 energy = energy + tmp1 #ifndef noVIRIAL vterms = 2.d0*(fac*msqs + 1.d0)*msq_inv tmp2 = tmp1*vterms vir(1,1) = vir(1,1) + tmp1*(vterms*mhat1s*mhat1s - 1.d0) vir(1,2) = vir(1,2) + tmp2*mhat1s*mhat2s vir(1,3) = vir(1,3) + tmp2*mhat1s*mhat3s vir(2,1) = vir(2,1) + tmp2*mhat2s*mhat1s vir(2,2) = vir(2,2) + tmp1*(vterms*mhat2s*mhat2s - 1.d0) vir(2,3) = vir(2,3) + tmp2*mhat2s*mhat3s vir(3,1) = vir(3,1) + tmp2*mhat3s*mhat1s vir(3,2) = vir(3,2) + tmp2*mhat3s*mhat2s vir(3,3) = vir(3,3) + tmp1*(vterms*mhat3s*mhat3s - 1.d0) #endif end if q(1,k3,k1,k2q) = eterm * q(1,k3,k1,k2q) q(2,k3,k1,k2q) = eterm * q(2,k3,k1,k2q) end do ! k1 = k10, nf1+1 end do ! k3 = 1,nfft3 end do ! k2q = 1, mxzslabs eer = 0.5d0 * energy #ifndef noVIRIAL do m2 = 1,3 do m1 = 1,3 vir(m1,m2) = 0.5d0*vir(m1,m2) end do end do #endif return end subroutine scalar_sumrc_orthog