! #include "copyright.h" #include "assert.h" #include "dprec.h" !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Determine if the list needs rebuilding. !----------------------------------------------------------------------- ! --- CHECK_SKIN --- !--------------------------------------------------------------------- ! Check if any atom has moved more than half the skin distance, ! which is half the nbskin added to the cutoff in generating the ! verlet list; in which case a list build is flagged. ! An obvious parallel implementation, in which each processor checks its ! atoms and communicates the results to all other processors, produced ! on Linux clusters large losses with large numbers of processors. ! The separate sections based on nbtell exist for improved performance; ! computing the list update info has a small performance cost. subroutine check_skin(crd,savcrd,do_list_update) use trace implicit none _REAL_ crd(3,*) ! current atom coordinates, intent(in) _REAL_ savcrd(3,*) ! atom coordinates used to make list, intent(in) logical do_list_update ! true if a list update is needed, intent(out) # include "constants.h" # include "ew_unitcell.h" # include "extra.h" # include "memory.h" integer first_atom integer last_atom integer i integer nmoved_atoms ! total number of atoms triggering a list update _REAL_ dx,dy,dz,dis2 _REAL_ maxdis2 call trace_enter( 'check_skin' ) steps_since_list_build = steps_since_list_build + 1 first_atom = 1 last_atom = natom maxdis2 = zero if ( nbtell == 0 ) then ! list update info not requested do i = first_atom, last_atom dx = crd(1,i) - savcrd(1,i) dy = crd(2,i) - savcrd(2,i) dz = crd(3,i) - savcrd(3,i) dis2 = dx**2 + dy**2 + dz**2 maxdis2 = max(dis2,maxdis2) end do do_list_update = maxdis2 > fourth * skinnb * skinnb else ! list update info requested nmoved_atoms = 0 do i = first_atom, last_atom dx = crd(1,i) - savcrd(1,i) dy = crd(2,i) - savcrd(2,i) dz = crd(3,i) - savcrd(3,i) dis2 = dx**2 + dy**2 + dz**2 if ( dis2 > fourth * skinnb * skinnb) then maxdis2 = max(dis2,maxdis2) nmoved_atoms = nmoved_atoms + 1 end if end do do_list_update = nmoved_atoms > 0 if ( do_list_update ) then if ( master ) then write(6, '(1x,A,I7,/,1x,A,F8.4,I7)') & 'List Build Triggered: Number of atoms triggering = ', & nmoved_atoms, ' Maximum distance moved = ',sqrt(maxdis2) end if end if end if ! ( nbtell == 0 ) call trace_logical('do_list_update = ', do_list_update) call trace_exit( 'check_skin' ) return end subroutine check_skin !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Nonbond list management. !----------------------------------------------------------------------- ! --- NONBOND_LIST --- !----------------------------------------------------------------------- ! Handles set-up and error checking for calling of ! get_nb_list which creates the nonbond list. subroutine nonbond_list(crd,iac,ico,iblo,inb,ntypes, & natom,x,ix,r_stack,i_stack, & ipairs,ntnb,ibelly,belly,newbalance,cn1, & #ifdef PSANDER ,v,vold,ntp,xr, & #endif qsetup ) use trace implicit none integer ntnb integer ipairs(*) _REAL_ crd(3,*), r_stack(*),cn1(*) integer i_stack(*) integer natom,iac(*),ico(*),iblo(*),inb(*),ntypes _REAL_ x(*) integer ix(*),ibelly(*) integer newbalance integer ngrd1,ngrd2,ngrd3 integer igridpairs,sizgrdprs integer last_numlist,listdiff,listdiffmax,isiz logical belly,trial,balance # include "extra.h" # include "ew_cntrl.h" # include "ew_unitcell.h" # include "ew_localnb.h" # include "ew_mpole.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 integer tmplist(0:MPI_MAX_PROCESSORS), & alllist(0:MPI_MAX_PROCESSORS) #else /* not parallel needs numtasks and mytaskid */ integer numtasks,mytaskid #endif # include "ew_numtasks.h" # include "ew_tran.h" # include "def_time.h" integer ier integer, dimension(:), allocatable :: nghbptr integer, dimension(:), allocatable :: nghtran integer, dimension(:), allocatable :: gridpairs ! ARGUMENTS: (all are input) ! CRD is the array holding atomic coords. ! IAC and ICO are used to look up vdw and hbond interactions. ! ICO is an NTYPES by NTYPES array giving lookup into VDW and HBOND ! coefficients. ! IAC(i) is the atom type of atom i. The coefficients ! for 6-12 interactions for atoms i and j are indexed by ! ICO(iac(i),iac(j)). In practice ICO is unrolled to a 1 ! dimensional array. They are needed here to split nonbond ! list into vdw,hbonds (see pack_nb_list()) ! IBLO and INB are used to mask out some nonbond interactions ! (bonded pairs,1-3,1-4 pairs) etc. IBLO(i) stores the number ! of atoms masked by atom i. INB stores their indices. The ! masked pairs are stored back to back in INB. ! For our purposes we need to produce our own variants. ! See "load_mask". ! NTYPES is used with IAC and ICO ! X and IX are real and integer arrays which comprise the total "dynamic" ! memory for amber; coords,forces,bond lists etc are accessed as ! offsets in them. integer ifail,listtot,listtotall,nucgrd,issetup integer inddel, i logical do_list_update logical qsetup logical first data first /.true./ #ifdef PSANDER integer l_tmpcrd,l_tmpvel,l_tmpxr,l_tmpdp,l_tmpvel2 integer istart,inext,ilast,inwat,itmpro _REAL_ v(*),xr(*),vold(*) integer ntp #endif #ifdef PSANDER ! include "memory.h" integer ictmplist,iftmplist #endif save trial,balance save last_numlist,listdiffmax save issetup data issetup/0/ !---------------------------- code starts here --------------------------- call trace_enter( 'nonbond_list' ) #ifndef MPI mytaskid=0 numtasks=1 #endif ! --------------------------------------------------------------------- ! ------ FIRST TIME SETUP STUFF -------------------------------- ! -------- If qsetup is true then this is the ! ------- first time through, do the setup stuff ! ------ do not set the qsetup flag false till after the ! ----- if(.not. qsetup ) block if ( qsetup ) then #ifdef PSANDER ! ----- for spatial decomp, ew_startup is called from runmd ! where some other things are also done as a ! one-time setup. #else call ew_startup(natom,iblo,inb,x,ix,ix(ixtran),r_stack) nucgrd1_0=-1 nucgrd2_0=-1 nucgrd3_0=-1 call save_crds(natom,crd,x(lsavcrd)) #endif balance=.false. newbalance=0 if(periodic == 0) balance= (numtasks > 1) if(balance)newbalance=2 last_numlist=0 steps_since_list_build = 0 end if ! ------ DONE FIRST TIME SETUP STUFF -------------------------------- trial=newbalance > 0 if ( .not. qsetup )then #ifdef PSANDER ! The skin check is done in runmd for spatial #else !--------------------------------------------------------------------- ! ------- SKIN (BUFFER) CHECK see if new list is required ! ------- Do this block except on first pass. ----------- ! ------ For nocutoff, the list will build once ! ----- by skipping over this block to the build list part ! ------- if the nocutoff flag is set, then there is ! ------ no need to update the list, all pairs are ! ----- already determined if (nocutoff) then call trace_exit( 'nonbond_list' ) return end if ! ----- the skin check is done every step. ! ---- in old amber method (pre skin), update only when ntnb .ne. 0 ! --- this logic is retained if nbflag = 0 ! -- do not do any of this the first step, i.e. until qsetup = .false. ! if nbflag=0, use old logic; i.e. only update if ntnb .ne. 0 ! if nbflag=1, the skin check determines it if ( nbflag == 0 )then if ( nbtell /= 0 ) then ! list update info requested if(master)write(6,'(1x,A,I2)') 'OLD LIST LOGIC; ntnb = ',ntnb end if if ( ntnb == 0 ) then call trace_exit( 'nonbond_list' ) return end if else call check_skin(crd,x(lsavcrd),do_list_update) if ( do_list_update ) then call save_crds(natom,crd,x(lsavcrd)) else call trace_exit( 'nonbond_list' ) return end if end if #endif /* PSANDER */ end if ! ( .not. qsetup ) qsetup = .false. #ifdef PSANDER ! Good candidate for removal, may be done elsewhere.... call save_crds(natom,crd,x(lsavcrd)) #endif !------------------------------------------------------------------------ ! START THE LIST BUILD, FIRST THE LIST GRID SETUP !--------------------------------------------------------------------- call timer_start(TIME_BLDLST) if ( master .and. nbtell /= 0 ) then ! list update info requested write(6, '(1x,A,I7,A)') 'Building list: ', & steps_since_list_build, ' steps since previous list build.' end if steps_since_list_build = 0 #ifdef MPI if(i_do_direct)then call mpi_comm_size(direct_comm,numtasks,ierr) call mpi_comm_rank(direct_comm,mytaskid,ierr) #endif ! ---- Get nb grid information ! This is necessary in the case where the system dimensions are ! changing: nonperiodic and const pressure call map_coords(crd,natom,x(lfrction),recip) call save_frac_crds(natom,x(lfrction),x(lsavfrac)) call setup_grids(nghb,cutlist,reclng,dirlng, & nucgrd1,nucgrd2,nucgrd3,nghb1,nghb2,nghb3, & #ifdef PSANDER ix(iindexlo),ix(iindexhi),myindexlo,myindexhi,inddelta,1, & ! task=1 full setup #endif periodic,nogrdptrs,verbose) if ( nucgrd1*nucgrd2*nucgrd3 > nucgmax ) & call sander_bomb('nonbond_list', & ' volume of ucell too big, too many subcells', & ' list grid memory needs to be reallocated, restart sander') #ifdef PSANDER if(nucgrd1 /= nucgrd1_0 & .or. nucgrd2 /= nucgrd2_0 & .or. nucgrd3 /= nucgrd3_0)then ! ---------------------------------------------------------------- ! ------ machinery for redoing the spatial decomp when the grid ! changes still needs to be engineered. call sander_bomb('nonbond_list', & ' nonbond grid changed dimensions', & ' Spatial cannot continue. ') endif #endif ! --------- Nonperiodic systems: -------------------------------------- if(periodic == 0)then ! ----- If system has too few cells for the pointer method ! ----- to be efficient, use the no-grid-pointer system ! ----- where all forward cells are checked for pair atoms ! ----- rather than just the forward cells on the grid-pointer-list nogrdptrs=nogrdptrs .or.( (nucgrd1 <= 2*nghb+2) .or. & (nucgrd1 <= 2*nghb+2) .or. & (nucgrd1 <= 2*nghb+2) ) if(verbose >= 3)write(6,*)" List using nogrdptrs ",nogrdptrs ! ----- Make the subcells about 3 A in size now so that there ! ----- are more subcells, easier to load balance. if(nogrdptrs)then if(dirlng(1)/nucgrd1 > 3)then ngrd1 = dirlng(1)/3 ngrd2 = dirlng(2)/3 ngrd3 = dirlng(3)/3 if ( verbose >= 1)then write(6,'("| New Grids set up for nogrdptrs ")') write(6, '(5X,a,/,5X,i9,1X,i9,1X,i9)') & 'Number of grids per unit cell in x,y,z:', & ngrd1, ngrd2, ngrd3 end if nucgrd1=ngrd1 nucgrd2=ngrd2 nucgrd3=ngrd3 end if end if end if ! ---- end of nonperiodic grid stuff ------------------------ call fill_tranvec() nucgrd = nucgrd1*nucgrd2*nucgrd3 #ifdef MPI !--------------------------------------------------------------------- ! For trial run of a parallel run, an initial guess at ! the distribution of subcells must be made. (trial=.true.) ! For a periodic system, the balance is assumed to be ! good for evenly dividing the subcells. (balance = .false.) if(trial .or. .not.balance)then # ifdef PSANDER ! ------ THese are set in spatial_startup(), ! if the grid has changed so these values change, then ! abort the simulation until new machinery is in place. ! The ultimate goal is to get rid of indexlo-indexhi ! and replace with a list of cells for more flexibility ! in shape of spatial regions. # else inddel = (nucgrd-1) / numtasks + 1 if ( inddel == 0 )inddel = 1 myindexlo = 1+(mytaskid)*inddel myindexhi = myindexlo+inddel-1 if ( mytaskid == numtasks-1 ) myindexhi = nucgrd # endif last_numlist = 0 end if #else !--------------------------------------------------------------------- ! ----- for non parallel runs, do all the subcells myindexlo=1 myindexhi=nucgrd #endif !--------------------------------------------------------------------- ! ----- assign atoms to cells (atmcell) and make cell atom list (indatg) call grid_ucell(natom,nucgrd1,nucgrd2,nucgrd3,nucgrd, & ix(iatmcell),ix(inumatg),ix(iindoff),ix(iindatg),x(lfrction) & ,periodic & #ifdef PSANDER ,ix(iknown_atm) & #endif ) !--------------------------------------------------------------------- if(periodic == 0)then !----- Nonperiodic systems: isiz=max(2,(myindexhi-myindexlo+1)*(maxnptrs+1)) if(nogrdptrs)isiz=1 allocate( nghbptr(isiz), stat=ier ) REQUIRE( ier==0 ) isiz=max(2,(myindexhi-myindexlo+1)*(maxnptrs)) if(nogrdptrs)isiz=1 allocate( nghtran(isiz), stat=ier ) REQUIRE( ier==0 ) call grid_pointers( & nucgrd1,nucgrd2,nucgrd3,nghb1,nghb2,nghb3, & nghbptr,maxnptrs,numnptrs, & nghtran,ix(imy_grids), & nucgrd,myindexlo,myindexhi,ix(ixtran), & periodic,nogrdptrs) else !----- Periodic systems: ! only call grid_pointers if the unit cell has changed ! so much that nucgrd[123] have changed. Otherwise the ! grid pointers do not change. if(nucgrd1 /= nucgrd1_0 & .or. nucgrd2 /= nucgrd2_0 & .or. nucgrd3 /= nucgrd3_0)then nucgrd1_0=nucgrd1 nucgrd2_0=nucgrd2 nucgrd3_0=nucgrd3 call grid_pointers( & nucgrd1,nucgrd2,nucgrd3,nghb1,nghb2,nghb3, & ix(inghbptr),maxnptrs,numnptrs, & ix(inghtran),ix(imy_grids), & nucgrd,myindexlo,myindexhi,ix(ixtran), & periodic,nogrdptrs) end if end if ! (periodic == 0) ! --------------------------------------------------------------------- call grid_image(nucgrd1,nucgrd2,nucgrd3, & maximage,numimage,ix(inumatg),ix(iindoff), & ix(iindatg),ix(inumimg), & ix(inlogrid),ix(inhigrid),x(lfrction), & x(limgcrds),ucell,ix(ibckptr),verbose,periodic & ,ix(imy_grids) & ) if(balance)then sizgrdprs=nucgrd*2 else sizgrdprs=2 end if allocate( gridpairs(sizgrdprs), stat=ier ) REQUIRE( ier==0 ) if(periodic == 0)then call get_nb_list(iac,ico,ntypes,ifail, & listtot,nucgrd,maxnblst,natom, & ix(iatmlist),ix(ilist),ix(iscratch), & ix(iiwa),ix(iiwh), & ipairs, & ix(inumimg),ix(inlogrid),ix(inhigrid), & numnptrs,maxnptrs,nghbptr,nghb1,ix(ibckptr), & ix(inumvdw),ix(inumhbnd),cutlist, & x(limgcrds),ix(imask),ix(inummask),ix(imaskptr),verbose, & ix(iitran),ix(iktran),nghtran,tranvec, & ix(ixtran),ix(imygrdlist),numimage,nucgrd1, & nucgrd2*nucgrd3, & myindexlo,myindexhi, & belly,ibelly,balance,gridpairs, & periodic,nogrdptrs,cn1) else call get_nb_list(iac,ico,ntypes,ifail, & listtot,nucgrd,maxnblst,natom, & ix(iatmlist),ix(ilist),ix(iscratch), & ix(iiwa),ix(iiwh), & ipairs, & ix(inumimg),ix(inlogrid),ix(inhigrid), & numnptrs,maxnptrs,ix(inghbptr),nghb1,ix(ibckptr), & ix(inumvdw),ix(inumhbnd),cutlist, & x(limgcrds),ix(imask),ix(inummask),ix(imaskptr),verbose, & ix(iitran),ix(iktran),ix(inghtran),tranvec, & ix(ixtran),ix(imygrdlist),numimage,nucgrd1, & nucgrd2*nucgrd3, & myindexlo,myindexhi, & belly,ibelly,balance,gridpairs, & periodic,nogrdptrs,cn1) end if ! (periodic == 0) if ( ifail == 1) then write(6, '(5x,a,i10)') 'SIZE OF NONBOND LIST = ', listtot call sander_bomb('nonbond_list','Non bond list overflow!', & 'check MAXPR in locmem.f') end if #ifdef MPI if(trial)then ASSERT ( periodic == 0 ) call trace_mpi('mpi_allreduce', nucgrd, 'MPI_INTEGER', mpi_sum) call mpi_allreduce(gridpairs, & gridpairs(1+nucgrd),nucgrd, & mpi_integer, mpi_sum,commsander,ierr) call fix_grid_balance(myindexlo,myindexhi, & gridpairs(1+nucgrd), & nucgrd,numtasks,mytaskid, & listdiffmax) ! call free_istack(igridpairs) ! call free_istack(inghtran) ! call free_istack(inghbptr) deallocate( gridpairs, nghtran, nghbptr ) isiz=max(2,(myindexhi-myindexlo+1)*(maxnptrs+1)) if(nogrdptrs)isiz=1 allocate( nghbptr(isiz), stat=ier ) REQUIRE( ier==0 ) isiz=max(2,(myindexhi-myindexlo+1)*(maxnptrs)) if(nogrdptrs)isiz=1 allocate( nghtran(isiz), stat=ier ) REQUIRE( ier==0 ) allocate( gridpairs(2*nucgrd1*nucgrd2*nucgrd3), stat=ier ) REQUIRE( ier==0 ) call grid_pointers( & nucgrd1,nucgrd2,nucgrd3,nghb1,nghb2,nghb3, & nghbptr,maxnptrs,numnptrs, & nghtran,ix(imy_grids), & nucgrd,myindexlo,myindexhi,ix(ixtran), & periodic,nogrdptrs) call grid_image(nucgrd1,nucgrd2,nucgrd3, & maximage,numimage,ix(inumatg),ix(iindoff), & ix(iindatg),ix(inumimg), & ix(inlogrid),ix(inhigrid),x(lfrction), & x(limgcrds),ucell,ix(ibckptr), & verbose,periodic & ,ix(imy_grids) & ) call get_nb_list(iac,ico,ntypes,ifail, & listtot,nucgrd,maxnblst,natom, & ix(iatmlist),ix(ilist),ix(iscratch), & ix(iiwa),ix(iiwh), & ipairs, & ix(inumimg),ix(inlogrid),ix(inhigrid), & numnptrs,maxnptrs,nghbptr,nghb1,ix(ibckptr), & ix(inumvdw),ix(inumhbnd),cutlist, & x(limgcrds),ix(imask),ix(inummask),ix(imaskptr),verbose, & ix(iitran),ix(iktran),nghtran,tranvec, & ix(ixtran),ix(imygrdlist),numimage,nucgrd1, & nucgrd2*nucgrd3, & myindexlo,myindexhi, & belly,ibelly,balance,gridpairs, & periodic,nogrdptrs,cn1) trial=.false. last_numlist=listtot newbalance=0 else if (balance) then ! --- If this was not a rebalance step, then check ! whether the list has gotten more than a tolerance ! away from the last balance. If so, trigger ! a new balance by setting trial and newbalance ! newbalance will need to be communicated to ! all processors before the next list build. listdiff=iabs(last_numlist-listtot) if(listdiff > listdiffmax)then newbalance=2 end if end if ! (trial) if( (nbtell > 1 ) .or. ( balance .and. (verbose > 0)) )then if(master) write(6,105) if(master) write(6,106) do i=0,numtasks-1 tmplist(i)=0 end do tmplist(mytaskid)=listtot call trace_mpi('mpi_allreduce', numtasks, 'MPI_INTEGER', mpi_sum) call mpi_allreduce(tmplist(0),alllist(0),numtasks, & mpi_integer, mpi_sum,commsander,ierr) if(master) then write(6,110)(i,alllist(i),i=0,numtasks-1) do i=1,numtasks-1 tmplist(0)=tmplist(0)+alllist(i) end do write(6,103)tmplist(0) end if 110 format("| ",2i12) 103 format("| Total: ",i12) 105 format("| ",'----------------List Breakdown----------------') 106 format("| ",'list processor listtot') end if #else if(master)then if ( verbose > 0)write(6,*)'listtot = ',listtot end if #endif deallocate( gridpairs ) if( periodic == 0 )then deallocate( nghtran, nghbptr ) end if listtotall=listtot #ifdef MPI if(first .or. (nbtell >= 1))then call trace_mpi('mpi_allreduce', 1, 'MPI_INTEGER', mpi_sum) call mpi_allreduce(listtot,listtotall,1, & mpi_integer, mpi_sum,commsander,ierr) end if call mpi_comm_rank(commsander,mytaskid,ierr) call mpi_comm_size(commsander,numtasks,ierr) end if ! (i_do_direct) call mpi_barrier(commsander,ierr) #endif ! if( (nbtell > 1) .and. .not. master ) write(6, '(a,i10,a,i4)') & ! '| Local SIZE OF NONBOND LIST = ', listtot,' PE',mytaskid if( (first .or. (nbtell == 1)) .and. master ) then write(6, '(a,i10)') & '| Local SIZE OF NONBOND LIST = ', listtot write(6, '(a,i10)') & '| TOTAL SIZE OF NONBOND LIST = ', listtotall end if call amflsh(6) first = .false. call timer_stop(TIME_BLDLST) call trace_exit( 'nonbond_list' ) return end subroutine nonbond_list !----------------------------------------------------- ! GET_EE_FUNC !----------------------------------------------------- !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine get_ee_func here] subroutine get_ee_func(x,switch,d_switch_dx,ee_type) implicit none _REAL_ x,switch,d_switch_dx integer ee_type # include "constants.h" ! ---get switch function multiplying the Coulomb interaction 1/r ! r has been converted to x by x = dxdr*r for convenience if ( ee_type == 1 )then ! ---erfc function for ewald call erfcfun(x,switch) d_switch_dx = -two*exp(-x*x)/sqrt(pi) else if ( ee_type == 2 )then ! ---force shift cutoff if ( x < one )then switch = (one - x)**2 d_switch_dx = -two*(one - x) else switch = zero d_switch_dx = zero end if end if return end subroutine get_ee_func !---------------------------------------------------------------------------- ! --- GET_NB_ENERGY -- !---------------------------------------------------------------------------- ! ...the main routine for non bond energy (vdw and hbond) ! as well as direct part of ewald sum. It is structured for parallelism. !---------------------------------------------------------------------------- subroutine get_nb_energy(iac,ico,ntypes,charge, & cn1,cn2,asol,bsol,force,numatoms, & nucgrd,numimg,nlogrid,nhigrid,bckptr, & numvdw,numhbnd,ipairs,imagcrds, & ewaldcof,eedtbdns,eed_cub,eed_lin, & maxnblst,eelt,evdw,ehb,virial,eedvir, & filter_cut,ee_type,eedmeth,dxdr, & epol,dipole,field,mpoltype,r_stack,i_stack, & myindexlo,myindexhi) implicit none #ifdef MPI # include "ew_parallel.h" # include "parallel.h" #endif # include "flocntrl.h" # include "constants.h" integer l_real_df,l_real_x,l_real_y,l_real_z,l_real_r2,l_int integer numatoms,maxnblst,mpoltype integer iac(*),ico(*),ntypes,nucgrd,ee_type,eedmeth integer numimg(*),nlogrid(*),nhigrid(*) integer bckptr(*),numvdw(*),numhbnd(*) _REAL_ charge(*),cn1(*),cn2(*),asol(*),bsol(*) _REAL_ imagcrds(3,*),ewaldcof,eedtbdns,dxdr, & eed_cub(4,*),eed_lin(2,*),virial(3,3) integer ipairs(maxnblst) _REAL_ force(3,numatoms),eelt,epol,evdw,ehb _REAL_ eedvir,filter_cut, & dipole(3,*),field(3,*) integer myindexlo,myindexhi integer index,numpack,i,k,l,ncell_lo,ncell_hi,ntot,nvdw _REAL_ xk,yk,zk integer i_stack(*) _REAL_ r_stack(*) integer nn,ncache ! ---FLOW CONTROL FLAG (debug and future respa) if ( do_dir == 0 )return eelt = zero epol = zero evdw = zero ehb = zero eedvir = zero do l = 1,3 do k = 1,3 virial(k,l) = zero end do end do numpack = 1 call timer_start(TIME_SHORT_ENE) do index = myindexlo,myindexhi k = index if ( numimg(k) > 0 )then ncell_lo = nlogrid(k) ncell_hi = nhigrid(k) do k = ncell_lo,ncell_hi i = bckptr(k) xk = imagcrds(1,k) yk = imagcrds(2,k) zk = imagcrds(3,k) ntot = numvdw(i) + numhbnd(i) nvdw = numvdw(i) if ( ntot > 0 )then if ( mpoltype == 0 )then ! allocate 6 temporary caches for performance optimizations ncache = max( nvdw, numhbnd(i) ) call get_stack(l_real_df,ncache) call get_stack(l_real_x,ncache) call get_stack(l_real_y,ncache) call get_stack(l_real_z,ncache) call get_stack(l_real_r2,ncache) call get_istack(l_int,ncache) call short_ene(i,xk,yk,zk,ipairs(numpack),ntot,nvdw, & bckptr,imagcrds,ewaldcof,eedtbdns, & eed_cub,eed_lin,charge, & ntypes,iac,ico,cn1,cn2,asol,bsol,filter_cut, & eelt,evdw,ehb,force,virial,ee_type,eedmeth,dxdr, & eedvir,r_stack(l_real_df),r_stack(l_real_x) & ,r_stack(l_real_y),r_stack(l_real_z) & ,r_stack(l_real_r2),i_stack(l_int) ) call free_stack(l_real_r2) call free_stack(l_real_z) call free_stack(l_real_y) call free_stack(l_real_x) call free_stack(l_real_df) call free_istack(l_int) else if ( mpoltype == 1 )then call short_ene_dip(i,xk,yk,zk,ipairs(numpack),ntot,nvdw, & bckptr,imagcrds,ewaldcof,eedtbdns, & eed_cub,eed_lin,charge,dipole, & ntypes,iac,ico,cn1,cn2,asol,bsol,filter_cut, & eelt,epol,evdw,ehb,force,field,virial, & ee_type,eedmeth,dxdr,eedvir) end if numpack = numpack + ntot end if ! ( ntot > 0 ) end do ! k = ncell_lo,ncell_hi end if ! ( numimg(k) > 0 ) end do ! index = myindexlo,myindexhi call timer_stop(TIME_SHORT_ENE) return end subroutine get_nb_energy !---------------------------------------------------------------------------- ! --- GET_NB_LIST --- ********* ! ...this routine carries out the short range part of nonbond ! calculation. Note that I am using a gridding routine, or ! geometric hashing with renumbering of atoms to speed list ! generation and to promote locality of reference, important ! for cache memory usage. The pre-imaging is to avoid expensive ! coordinate transforms in the minimum image pair calculations. ! I set up pointers between the subcells to speed this process. ! NOTE: this comment also refers to routine pack_nb_list(). !---------------------------------------------------------------------------- subroutine get_nb_list(iac,ico,ntypes,ifail, & listtot,nucgrd,maxnblst,numatoms, & atmlist,list,scratch,iwa,iwh,ipairs, & numimg,nlogrid,nhigrid,numnptrs,maxnptrs, & nghbptr,nghb1, bckptr,numvdw,numhbnd,cutoff_list, & imagcrds,mask,nummask,maskptr,verbose, & itran,ktran,nghtran, & tranvec,xtran,mygrdlist,numimage,nucgrd1,nucgrd23, & myindexlo,myindexhi, & belly,ibelly,balance,gridpairs,periodic,nogrdptrs,cn1) ! ---main routine for list building ! note it is structured for parallelism implicit none # include "extra.h" #ifdef MPI # include "ew_parallel.h" # include "parallel.h" # include "mpif.h" #else /* not parallel needs numtasks and mytaskid */ integer numtasks,mytaskid #endif _REAL_ cutoff_list,imagcrds(3,*),cn1(*) _REAL_ tranvec(*) integer iac(*),ico(*),ntypes,ifail,listtot integer nucgrd,maxnblst,numatoms,maxnptrs integer atmlist(numatoms), & list(numatoms),scratch(numatoms), & iwa(numatoms),iwh(numatoms) & ,mygrdlist(*),numimage,nucgrd1,nucgrd23 integer xtran(7,10) integer ipairs(maxnblst) integer myindexlo,myindexhi integer ibelly(*) logical belly,balance integer np0,gridpairs(nucgrd,2) integer numimg(*),nlogrid(*),nhigrid(*) integer bckptr(*),numvdw(*),numhbnd(*),numnptrs,nghb1, & nghbptr(0:maxnptrs,*),mask(*),maskptr(*),nummask(*), & verbose integer itran(*),ktran(*),nghtran(maxnptrs,*) integer numlist,num integer index,index0,index1 integer j,jj,j0,k,ncell_lo,ncell_hi,m,m1,m2 integer i,kk,numpack integer nstart,tranyz,xtindex,jtran,indi,nstop,numindex _REAL_ xk,yk,zk _REAL_ cutoffsq integer periodic logical nogrdptrs #ifndef MPI numtasks=1 mytaskid=0 #endif do i=1,numimage mygrdlist(i)=0 end do ifail = 0 cutoffsq = cutoff_list*cutoff_list do j = 1,numatoms scratch(j) = 0 end do numpack = 0 # ifdef MPI if(balance)then np0=0 do i=1,nucgrd gridpairs(i,1)=0 end do end if # endif do index = myindexlo,myindexhi index0=index-myindexlo+1 if ( numimg(index) > 0 )then ncell_lo = nlogrid(index) ncell_hi = nhigrid(index) numlist = 0 !========================================================== if(nogrdptrs)then !========================================================== if(periodic /= 0) & call sander_bomb('get_nb_list', & 'Cannot run nogrdptrs with ntb=1', & 'turn off nogrdptrs ') do j0 = index,nucgrd m1 = nlogrid(j0) m2 = nhigrid(j0) if ( m2 >= m1 )then do m = m1,m2 numlist = numlist+1 atmlist(numlist) = m end do end if end do if ( numlist > 0 )then do k = ncell_lo,ncell_hi kk = k - ncell_lo + 1 i = bckptr(k) xk = imagcrds(1,k) yk = imagcrds(2,k) zk = imagcrds(3,k) mygrdlist(k)=1 call pack_nb_nogrdptrs(kk,i,xk,yk,zk,imagcrds,cutoffsq, & atmlist,numlist,list,num,numpack,bckptr, & iac,ico,ntypes,iwa,iwh,ipairs, & numvdw,numhbnd,maxnblst,ifail, & mask,nummask,maskptr,scratch, & mygrdlist, & belly,ibelly) if ( ifail == 1)then listtot = maxnblst return end if end do end if !========================================================== else ! nogrdptrs if statement !========================================================== ! -- get list of atoms in the neighborhood of cell(index0) numindex= numnptrs if(periodic == 0) numindex = nghbptr(0,index0) do j0 = 1,numindex index1 = nghbptr(j0,index0) if ( j0 == 1 )then nstart=0 else nstart=-nghb1 end if nstop=nghb1 if(periodic == 0)then indi=mod(index1-1,nucgrd1)+1 if(indi > nucgrd1-nghb1)nstop=nucgrd1-indi if(indi <= nghb1)nstart=1-indi if ( j0 == 1 ) nstart=0 end if xtindex=ishft(nghtran(j0,index0),-8) REQUIRE(xtindex <= 10 .and. xtindex >= 0) tranyz = nghtran(j0,index0)-ishft(xtindex,8) do j = nstart,nstop jtran = tranyz+xtran(j-nstart+1,xtindex) jj = index1+j-xtran(j-nstart+1,xtindex)*nucgrd1 m1 = nlogrid(jj) m2 = nhigrid(jj) if ( m2 >= m1 )then do m = m1,m2 numlist = numlist+1 atmlist(numlist) = m itran(numlist)=jtran end do end if end do end do ! j0 = 1,numindex if ( numlist > 0 )then do k = ncell_lo,ncell_hi kk = k - ncell_lo + 1 i = bckptr(k) xk = imagcrds(1,k) yk = imagcrds(2,k) zk = imagcrds(3,k) mygrdlist(k)=1 call pack_nb_list(kk,i,xk,yk,zk,imagcrds,cutoffsq, & atmlist,numlist,list,num,numpack,bckptr, & iac,ico,ntypes,iwa,iwh,ipairs, & numvdw,numhbnd,maxnblst,ifail, & mask,nummask,maskptr,scratch, & itran,ktran,tranvec,mygrdlist, & belly,ibelly,cn1) if ( ifail == 1)then listtot = maxnblst return end if end do end if !========================================================== end if ! nogrdptrs if statement !========================================================== end if ! ( numimg(index) > 0 ) # ifdef MPI if(balance)then gridpairs(index,1)=numpack-np0 np0=numpack end if # endif end do ! index = myindexlo,myindexhi listtot = numpack return end subroutine get_nb_list !------------------------------------------------------------------- ! --- GRID_IMAGE --- ! ...this routine grids the imaged atoms in the unit cell plus ! neighbor cells while building the array of sorted image atoms. ! The sorting should improve locality of reference... !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine grid_image(nucgrd1,nucgrd2,nucgrd3, & maximage,numimage,numatg,indoff,indatg, & numimg,nlogrid,nhigrid,fraction,imagcrds,ucell,bckptr, & verbose,periodic & ,my_grids & ) implicit none # include "constants.h" # include "extra.h" integer nucgrd1,nucgrd2,nucgrd3 integer maximage,numimage integer numatg(*),indatg(*),indoff(*) integer nlogrid(*),nhigrid(*) integer numimg(*),bckptr(maximage),verbose,periodic _REAL_ imagcrds(3,maximage) _REAL_ fraction(3,*),ucell(3,3) _REAL_ f1,f2,f3,shft integer index,nucgrd,i,j,n & ,my_grids(*) integer i0,i1,i2 shft=half if(periodic == 0)shft=zero nucgrd = nucgrd1*nucgrd2*nucgrd3 numimage = 0 do index = 1,nucgrd if(my_grids(index) == 1)then n = numatg(index) if (numimage + n > maximage ) then ! Maximage is only estimated via a density plus fudge factor ! see local_nb_get_sizes.f. So check for overflow. ! (happens to us with bad ! crystal set ups i.e. inappropriate periodic bdries, ! or with big cell fluctuations write(6,*)'maximage = ',maximage write(6,*)'numimage = ',numimage call sander_bomb('grid_image', & 'num image atoms exceeds maximage!!', & 'check ewaldlocmem ') end if i0 = indoff(index) i1 = i0 + 1 i2 = i0 + n nlogrid(index) = numimage + 1 nhigrid(index) = numimage + n numimg(index) = n do i = i1,i2 j = indatg(i) f1 = fraction(1,j)+shft f2 = fraction(2,j)+shft f3 = fraction(3,j)+shft numimage = numimage + 1 imagcrds(1,numimage) = f1*ucell(1,1) + f2*ucell(1,2) + & f3*ucell(1,3) imagcrds(2,numimage) = f1*ucell(2,1) + f2*ucell(2,2) + & f3*ucell(2,3) imagcrds(3,numimage) = f1*ucell(3,1) + f2*ucell(3,2) + & f3*ucell(3,3) bckptr(numimage) = j end do end if ! (my_grids(index) == 1) end do ! index = 1,nucgrd if(master) then if (verbose == 1)write(6,*)'grid_image: numimage = ',numimage end if return end subroutine grid_image !------------------------------------------------------------------- ! --- GRID_POINTERS --- !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine grid_pointers here] subroutine grid_pointers( & nucgrd1,nucgrd2,nucgrd3,nghb1,nghb2,nghb3, & nghbptr,maxnptrs,numnptrs,nghtran,my_grids, & nucgrd,myindexlo,myindexhi,xtran, & periodic,nogrdptrs) implicit none integer nucgrd1,nucgrd2,nucgrd3,nghb1,nghb2,nghb3 integer xtran(7,10) integer maxnptrs,nghbptr(0:maxnptrs,*),numnptrs & ,nghtran(maxnptrs,*) integer my_grids(*) integer nucgrd,myindexlo,myindexhi,index0 integer periodic integer i1,i2,i3,j2,j3,i2grd,i3grd,j3grd,j2grd,jj2 integer j2strt,j2stop,jj2strt,jj2stop,j3stop integer index,num integer k3off integer xtindex_1,xtindex_2 integer sizgrd12 integer index0_min,index0_max logical nogrdptrs index0_min = maxnptrs index0_max = 0 sizgrd12=nucgrd1*nucgrd2 do i1=1,nucgrd my_grids(i1)=0 end do ! ---get forward pointers for unit cell grid ! neighbor cells must be later in the subcell list ! to get unique atom - imageatom pairing !------------------------------------------------------------------------- ! Non-Periodic cell list pointers ! Will get ri if(periodic == 0)then if(nogrdptrs)then do i1=myindexlo,nucgrd1*nucgrd2*nucgrd3 my_grids(i1)=1 end do else do i3 = 1,nucgrd3 i3grd=sizgrd12*(i3-1) do i2 = 1,nucgrd2 i2grd=i3grd+nucgrd1*(i2-1) do i1 = 1,nucgrd1 index = i2grd+i1 if(index < myindexlo.or.index > myindexhi)goto 180 my_grids(index)=1 index0=index-myindexlo+1 index0_max = max (index0_max,index0) index0_min = min (index0_min,index0) ! ! Set up xtran index for this cell ! ! UNSTABLE till fixed: ! Need to make sure that atoms do not drift to edges ! Otherwise this simplified xtindex will not work. xtindex_1 = 1 xtindex_2 = 1 j2strt=-nghb2 if(i2+j2strt < 1)j2strt=1-i2 j2stop=nghb2 if(i2+j2stop > nucgrd2)j2stop=nucgrd2-i2 jj2strt=-nghb1 if(i1+jj2strt < 1)jj2strt=1-i1 jj2stop=nghb1 if(i1+jj2stop > nucgrd1)jj2stop=nucgrd1-i1 ! ! first do the plane that this cell is in. ! first row of neighbors starts with the cell itself (Cell0) ! num = 1 nghbptr(num,index0) = index ! ! no y or z translate, this cell HAS to be in uc ! nghtran(num,index0) = 5+ishft(1,8) do jj2=1,jj2stop my_grids(index+jj2)=1 end do ! !next rows of neighbors are the nghb2 rows of 2*nghb1+1 cells ! above. Use the center cell of the row as the pointer ! (with the same x position as Cell0) !This way the pointer cell must be in the uc wrt the x ! direction. do j2 = 1,j2stop j2grd=index+j2*nucgrd1 if(j2+i2 <= nucgrd2)then num=num+1 nghtran(num,index0)=5+ishft(xtindex_2,8) nghbptr(num,index0) = j2grd do jj2=jj2strt,jj2stop my_grids(j2grd+jj2 & -nucgrd1*xtran(jj2+1+nghb1,xtindex_2))=1 end do end if end do ! ! Now there are nghb3 planes of (2*nghb1+1)(2*nghb2+1) cells ! ahead that are good neighbors. ! j3stop=nghb3 if(i3+j3stop > nucgrd3)j3stop=nucgrd3-i3 do j3 = 1,j3stop j3grd=j3*sizgrd12 do j2 = j2strt,j2stop num=num+1 j2grd=index+j3grd+j2*nucgrd1 nghtran(num,index0)= & 5+ishft(xtindex_2,8) nghbptr(num,index0) = j2grd do jj2=jj2strt,jj2stop my_grids(j2grd+jj2)=1 end do end do end do nghbptr(0,index0) = num 180 continue end do ! i1 = 1,nucgrd1 end do ! i2 = 1,nucgrd2 end do ! i3 = 1,nucgrd3 end if ! (nogrdptrs) numnptrs = maxnptrs return end if ! (periodic == 0) !------------------------------------------------------------------------- do i3 = 1,nucgrd3 i3grd=sizgrd12*(i3-1) do i2 = 1,nucgrd2 i2grd=i3grd+nucgrd1*(i2-1) do i1 = 1,nucgrd1 index = i2grd+i1 if( (index >= myindexlo) .and. (index <= myindexhi) ) then my_grids(index)=1 index0=index-myindexlo+1 index0_max = max (index0_max,index0) index0_min = min (index0_min,index0) ! ! Set up xtran index for this cell ! if(i1+nghb1 > nucgrd1)then xtindex_1 = 1+i1+nghb1-nucgrd1 xtindex_2 = xtindex_1 + 2*nghb1 else if(i1-nghb1 < 1)then xtindex_1 = 1 xtindex_2 = 2*nghb1+2-i1 else xtindex_1 = 1 xtindex_2 = 1 end if ! ! first do the plane that this cell is in. ! first row of neighbors starts with the cell itself (Cell0) ! num = 1 nghbptr(num,index0) = index ! ! no y or z translate, this cell HAS to be in uc ! nghtran(num,index0) = 5+ishft(xtindex_1,8) do jj2=1,nghb1 my_grids(index+jj2 & -nucgrd1*xtran(jj2+1,xtindex_1))=1 end do ! ! next rows of neighbors are the nghb2 rows of 2*nghb1+1 cells ! above. Use the center cell of the row as the pointer ! (with the same x position as Cell0) ! This way the pointer cell must be in the uc wrt the x direction. ! do j2 = 1,nghb2 num=num+1 j2grd=index+j2*nucgrd1 if(j2+i2 <= nucgrd2)then nghtran(num,index0)=5+ishft(xtindex_2,8) else j2grd=j2grd-sizgrd12 nghtran(num,index0)=8+ishft(xtindex_2,8) end if nghbptr(num,index0) = j2grd do jj2=-nghb1,nghb1 my_grids(j2grd+jj2 & -nucgrd1*xtran(jj2+1+nghb1,xtindex_2))=1 end do end do ! ! Now there are nghb3 planes of (2*nghb1+1)(2*nghb2+1) cells ! ahead that are good neighbors. ! do j3 = 1,nghb3 j3grd=j3*sizgrd12 if(j3+i3 > nucgrd3)then j3grd=j3grd-sizgrd12*nucgrd3 k3off=9 else k3off=0 end if do j2 = -nghb2,nghb2 num=num+1 j2grd=index+j3grd+j2*nucgrd1 if(j2+i2 > nucgrd2)then nghtran(num,index0)=8+k3off+ishft(xtindex_2,8) j2grd=j2grd-sizgrd12 else if(j2+i2 < 1)then nghtran(num,index0)=2+k3off+ishft(xtindex_2,8) j2grd=j2grd+sizgrd12 else nghtran(num,index0)=5+k3off+ishft(xtindex_2,8) end if nghbptr(num,index0) = j2grd do jj2=-nghb1,nghb1 my_grids(j2grd+jj2 & -nucgrd1*xtran(jj2+1+nghb1,xtindex_2))=1 end do end do end do endif ! index within lo and hi end do ! 280 i1 = 1,nucgrd1 end do ! i2 = 1,nucgrd2 end do ! i3 = 1,nucgrd3 numnptrs = num return end subroutine grid_pointers !------------------------------------------------------------------- ! --- GRID_UCELL --- ! ...this routine grids the mapped atoms in the unit cell ! into the NUCGRD1 x NUCGRD2 x NUCGRD3 subcells according ! to their fractional coordinates. !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ subroutine grid_ucell(natom,nucgrd1,nucgrd2,nucgrd3,nucgrd, & atmcell,numatg,indoff,indatg,fraction,periodic & #ifdef PSANDER ,known_atm & #endif ) implicit none integer natom,nucgrd1,nucgrd2,nucgrd3,nucgrd,periodic integer atmcell(natom),numatg(nucgrd),indoff(nucgrd),indatg(natom) #ifdef PSANDER integer known_atm(natom) #endif _REAL_ fraction(3,natom) # include "constants.h" integer i,j,i1,i2,i3,index integer ichk _REAL_ shft shft=half if(periodic == 0)shft=0.0 do index = 1,nucgrd numatg(index) = 0 end do !----------------------------------------------------------------- ! ---find out which ucgrd subcell each atom is in. ! atmcell(i) is the ucgrd subcell that contains atom i ! numatg(index) is the number of atoms in the index ucgrd do i = 1,natom #ifdef PSANDER if(known_atm(i)==1)then #endif i1 = (fraction(1,i) + shft) * nucgrd1 + 1 i2 = (fraction(2,i) + shft) * nucgrd2 + 1 i3 = (fraction(3,i) + shft) * nucgrd3 + 1 index = nucgrd1*nucgrd2*(i3-1)+nucgrd1*(i2-1)+i1 atmcell(i) = index numatg(index) = numatg(index) + 1 #ifdef PSANDER endif #endif end do ! ! ---find the offset of the starting atoms for each ucgrd subcell. ! zero the numatg()entries as you go indoff(1) = 0 do i = 2,nucgrd indoff(i) = indoff(i-1) + numatg(i-1) numatg(i-1) = 0 end do ! ichk = natom - indoff(nucgrd) - numatg(nucgrd) ! this crash also only seems to occur when atoms can not be gridded due to ! truly awful coordinates e.g. NaN (We get with really bad input data ! sometimes) ! if ( ichk .ne. 0 )then ! write(6,*)'BIG PROBLEM in grid_ucell!!!' ! call mexit(6,1) ! endif !----------------------------------------------------------------- ! Fill indatg() as a list of atoms in each subcell such that ! the list of atoms in subcell 1 are at the beginning, ! subcell 2 list is right after that (starting at indoff(2)+1) ! numatg(nucgrd) = 0 do i = 1,natom #ifdef PSANDER if(known_atm(i)==1)then #endif index = atmcell(i) numatg(index) = numatg(index) + 1 j = numatg(index)+indoff(index) indatg(j) = i #ifdef PSANDER endif #endif end do return end subroutine grid_ucell !*********************************************************************** ! ------------------------------------------------------------------- ! --- PACK_NB_LIST --- *** ! ------------------------------------------------------------------- !*********************************************************************** !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine pack_nb_list here] #ifndef FUJITSU_VPP subroutine pack_nb_list(kk,i,xk,yk,zk,imagcrds,cutoffsq, & atmlist,numlist,list,num,numpack,bckptr, & iac,ico,ntypes,iwa,iwh, & ipairs, & numvdw,numhbnd,maxnblst,ifail, & mask,nummask,maskptr,scratch, & itran,ktran,tranvec,mygrdlist, & belly,ibelly,cn1) implicit none integer numpack integer kk,i,atmlist(*),numlist,list(*),num,bckptr(*), & iac(*),ico(*),ntypes,iwa(*),iwh(*), & numvdw(*),numhbnd(*), & maxnblst,ifail & ,mygrdlist(*) _REAL_ xk,yk,zk,imagcrds(3,*),cutoffsq,cn1(*) integer mask(*),nummask(*),maskptr(*),scratch(*) integer ipairs(*),ibelly(*) logical belly,deadi,deadik # include "parallel.h" # include "constants.h" _REAL_ dx,dy,dz,r2 integer iaci,ic,index,k,lpr,lhb,m,n,npr integer itran(*),ktran(*),jtran _REAL_ x_tran,y_tran,z_tran,tranvec(3,*) num = 0 m = maskptr(i) lpr = 0 lhb = 0 iaci = ntypes*(iac(i)-1) do n = 1,nummask(i) k = mask(m+n) scratch(k) = i end do deadi=.false. deadi=(ibelly(i) == 0).and.belly do m = kk+1,numlist jtran=itran(m) ! x_tran=tranvec(1,jtran) ! y_tran=tranvec(2,jtran) ! z_tran=tranvec(3,jtran) x_tran=tranvec(1,itran(m)) y_tran=tranvec(2,itran(m)) z_tran=tranvec(3,itran(m)) n = atmlist(m) k = bckptr(n) deadik=(ibelly(k) == 0).and.deadi if ( (scratch(k) /= i) .and. .not.deadik) then dx = imagcrds(1,n) - xk + x_tran dy = imagcrds(2,n) - yk + y_tran dz = imagcrds(3,n) - zk + z_tran r2 = dx*dx + dy*dy + dz*dz if ( r2 < cutoffsq ) then #ifdef TEST_CLOSEATOMS if(r2 < .5) & write(6,190)i,k,r2 190 format(" Atoms close:",2i8, & " r**2 ",f10.6) #endif if (mygrdlist(n) == 0) then mygrdlist(n)=1 end if num = num + 1 list(num) = n index = iaci+iac(k) ic = ico(index) if (ic >= 0) then if (cn1(ic) /= zero) then lpr = lpr+1 iwa(lpr) = ior(n,ishft(itran(m),27)) ! bitwise optimization else lhb = lhb+1 iwh(lhb) = ior(n,ishft(itran(m),27)) ! bitwise optimization end if else lhb = lhb+1 iwh(lhb) = ior(n,ishft(itran(m),27)) ! bitwise optimization end if end if end if end do ! m = kk+1,numlist if ( num + numpack > maxnblst )then write(6, '(/,a,i12,i12,a,i12,a,i4)') & ' * NB pairs ', num,numpack, & ' exceeds capacity (',maxnblst,')', mytaskid ifail=1 return end if ! ---now put all pairs into iwa do k = 1,lhb iwa(lpr+k) = iwh(k) end do npr = lpr+lhb numvdw(i) = lpr numhbnd(i) = lhb do k = 1,npr numpack = numpack+1 ipairs(numpack) = iwa(k) ! write(21,*)numpack,ishft(iwa(k),-27),iwa(k)-ishft(jtran,27) end do return end subroutine pack_nb_list #else ! Code optimized for vector processing subroutine pack_nb_list(kk,i,xk,yk,zk,imagcrds,cutoffsq, & atmlist,numlist,list,num,numpack,bckptr, & iac,ico,ntypes,iwa,iwh, & ipairs, & numvdw,numhbnd,maxnblst,ifail, & mask,nummask,maskptr,scratch, & itran,ktran,tranvec,mygrdlist, & belly,ibelly,cn1) implicit none integer numpack integer kk,i,atmlist(*),numlist,list(*),num,bckptr(*), & iac(*),ico(*),ntypes,iwa(*),iwh(*), & numvdw(*),numhbnd(*), & maxnblst,ifail & ,mygrdlist(*) _REAL_ xk,yk,zk,imagcrds(3,*),cutoffsq,cn1(*) integer mask(*),nummask(*),maskptr(*),scratch(*) integer ipairs(*),ibelly(*) logical belly,deadi,deadik # include "parallel.h" # include "constants.h" _REAL_ dx,dy,dz,r2 integer iaci,ic,index,k,lpr,lhb,m,n,npr integer itran(*),ktran(*),jtran _REAL_ x_tran,y_tran,z_tran,tranvec(3,*) integer max_vec parameter (max_vec = 50000) integer mm, ic_vec(max_vec), ior_vec(max_vec), numtmp integer n_vec(max_vec) _REAL_ r2_vec(max_vec) if (max_vec .lt. numlist) then write(6,*) "disaster in VPP version of pack_nb_list:" write(6,*) "max_vec too small" call mexit(6, 1) endif num = 0 m = maskptr(i) lpr = 0 lhb = 0 iaci = ntypes*(iac(i)-1) do n = 1,nummask(i) k = mask(m+n) scratch(k) = i end do deadi=.false. deadi=(ibelly(i) == 0).and.belly mm=0 do m = kk+1,numlist mm = mm + 1 jtran=itran(m) x_tran=tranvec(1,jtran) y_tran=tranvec(2,jtran) z_tran=tranvec(3,jtran) !x_tran=tranvec(1,itran(m)) !y_tran=tranvec(2,itran(m)) !z_tran=tranvec(3,itran(m)) n = atmlist(m) k = bckptr(n) dx = imagcrds(1,n) - xk + x_tran dy = imagcrds(2,n) - yk + y_tran dz = imagcrds(3,n) - zk + z_tran r2_vec(mm) = dx*dx + dy*dy + dz*dz index = iaci+iac(k) ic_vec(mm) = ico(index) ior_vec(mm) = ior(n,ishft(jtran,27)) enddo numtmp = num mm = 0 do m = kk+1,numlist mm = mm + 1 n = atmlist(m) k = bckptr(n) deadik=(ibelly(k) == 0).and.deadi if ( (scratch(k) /= i) .and. .not.deadik) then if ( r2_vec(mm) < cutoffsq ) then #ifdef TEST_CLOSEATOMS if(r2 < .5) & write(6,190)i,k,r2 190 format(" Atoms close:",2i8, & " r**2 ",f10.6) #endif numtmp = numtmp + 1 list(numtmp) = n n_vec(numtmp) = n if(ic_vec(mm).ge.0) then lpr = lpr+1 iwa(lpr) = ior_vec(mm) else lhb = lhb+1 iwh(lhb) = ior_vec(mm) endif endif endif end do do m = num+1, numtmp n = n_vec(m) if (mygrdlist(n).eq.0) then mygrdlist(n)=1 endif enddo num = numtmp if ( num + numpack .gt. maxnblst )then write(6, '(/,a,i12,i12,a,i12,a,i4)') & ' * NB pairs ', num,numpack, & ' exceeds capacity (',maxnblst,')', mytaskid ifail=1 return endif ! ---now put all pairs into iwa do k = 1,lhb iwa(lpr+k) = iwh(k) end do npr = lpr+lhb numvdw(i) = lpr numhbnd(i) = lhb do k = 1,npr numpack = numpack+1 ipairs(numpack) = iwa(k) ! write(21,*)numpack,ishft(iwa(k),-27),iwa(k)-ishft(jtran,27) end do return end subroutine pack_nb_list #endif !*********************************************************************** ! ------------------------------------------------------------------- ! --- PACK_NB_NOGRDPTRS --- *** ! ------------------------------------------------------------------- !*********************************************************************** !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine pack_nb_nogrdptrs here] subroutine pack_nb_nogrdptrs(kk,i,xk,yk,zk,imagcrds,cutoffsq, & atmlist,numlist,list,num,numpack,bckptr, & iac,ico,ntypes,iwa,iwh, & ipairs, & numvdw,numhbnd,maxnblst,ifail, & mask,nummask,maskptr,scratch, & mygrdlist, & belly,ibelly) implicit none integer numpack integer kk,i,atmlist(*),numlist,list(*),num,bckptr(*), & iac(*),ico(*),ntypes,iwa(*),iwh(*), & numvdw(*),numhbnd(*), & maxnblst,ifail & ,mygrdlist(*) _REAL_ xk,yk,zk,imagcrds(3,*),cutoffsq integer mask(*),nummask(*),maskptr(*),scratch(*) integer ipairs(*),ibelly(*) logical belly,deadi,deadik integer jtran # include "parallel.h" _REAL_ dx,dy,dz,r2 integer iaci,ic,index,k,lpr,lhb,m,n,npr num = 0 m = maskptr(i) lpr = 0 lhb = 0 iaci = ntypes*(iac(i)-1) do n = 1,nummask(i) k = mask(m+n) scratch(k) = i end do deadi=.false. deadi=(ibelly(i) == 0).and. belly do m = kk+1,numlist jtran=5 n = atmlist(m) k = bckptr(n) deadik=(ibelly(k) == 0).and.deadi if ( (scratch(k) /= i) .and. .not.deadik) then dx = imagcrds(1,n) - xk dy = imagcrds(2,n) - yk dz = imagcrds(3,n) - zk r2 = dx*dx + dy*dy + dz*dz if ( r2 < cutoffsq ) then #ifdef TEST_CLOSEATOMS if(r2 < .5) & write(6,190)i,k,r2 190 format(" Atoms close:",2i8, & " r**2 ",f10.6) #endif if (mygrdlist(n) == 0) then mygrdlist(n)=1 end if num = num + 1 list(num) = n index = iaci+iac(k) ic = ico(index) if(ic >= 0) then lpr = lpr+1 iwa(lpr) = ior(n,ishft(5,27)) ! bitwise optimization else lhb = lhb+1 iwh(lhb) = ior(n,ishft(5,27)) ! bitwise optimization end if end if end if end do ! m = kk+1,numlist if ( num + numpack > maxnblst )then write(6, '(/,a,i12,i12,a,i12,a,i4)') & ' * NB pairs ', num,numpack, & ' exceeds capacity (',maxnblst,')', mytaskid ifail=1 return end if ! ---now put all pairs into iwa do k = 1,lhb iwa(lpr+k) = iwh(k) end do npr = lpr+lhb numvdw(i) = lpr numhbnd(i) = lhb do k = 1,npr numpack = numpack+1 ipairs(numpack) = iwa(k) end do return end subroutine pack_nb_nogrdptrs !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Store the atomic coordinates. ! !----------------------------------------------------------------------- ! --- SAVE_CRDS --- !----------------------------------------------------------------------- ! Save the atomic coordinates that were used to build the list. ! This is needed for the skin test for buffered pairlists, ! where the current coordinates are compared with the saved ! coordinates relative to the skin criterion; see check_skin. !-------------------------------------------------------------------- subroutine save_crds(numatoms,crd,savcrd) implicit none integer numatoms _REAL_ crd(3*numatoms),savcrd(3*numatoms) ! These are actually two-dimensional (3,numatoms), but to enable ! vectorization on IA32 SSE platforms they are treated as ! one-dimensional; this may also improve software pipelining ! integer i do i = 1,3*numatoms savcrd(i) = crd(i) end do return end subroutine save_crds !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Store the fractional coordinates. ! !----------------------------------------------------------------------- ! --- SAVE_FRAC_CRDS --- !----------------------------------------------------------------------- ! This is needed in case you are wrapping coords in the box; ! note, td actually runs ewald without wrapping since it is easier ! to analyze but it should work equally well, either way. subroutine save_frac_crds(numatoms,fraction,savfrac) implicit none integer numatoms _REAL_ fraction(3*numatoms),savfrac(3*numatoms) ! These are actually two-dimensional (3,numatoms), but to enable ! vectorization on IA32 SSE platforms they are treated as ! one-dimensional; this may also improve software pipelining ! integer i do i = 1,3*numatoms savfrac(i) = fraction(i) end do return end subroutine save_frac_crds !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Calculate the direct Ewald component of the potentials !+ with polarizabilities. ! !----------------------------------------------------------------------- ! --- SHORT_ENE_DIP --- !----------------------------------------------------------------------- subroutine short_ene_dip(i,xk,yk,zk,ipairs,numtot,numvdw, & bckptr,imagcrds,ewaldcof,eedtbdns, & eed_cub,eed_lin,charge,dipole, & ntypes,iac,ico,cn1,cn2,asol,bsol,filter_cut, & eelt,epol,evdw,ehb,frc,field,virial, & ee_type,eedmeth,dxdr,eedvir) implicit none _REAL_ xk,yk,zk,imagcrds(3,*) integer i,numvdw,numtot,bckptr(*) integer ipairs(*),ee_type,eedmeth _REAL_ ewaldcof,eed_cub(4,*),eed_lin(2,*), & charge(*),dipole(3,*),virial(3,3),eedvir _REAL_ eedtbdns,filter_cut,dxdr integer ntypes,iac(*),ico(*) _REAL_ cn1(*),cn2(*),asol(*),bsol(*), & eelt,epol,evdw,ehb,frc(3,*),field(3,*) integer ic,j,m,n,ind,iaci _REAL_ del _REAL_ switch,d_switch_dx _REAL_ ee_vir_iso _REAL_ edx,edy,edz _REAL_ b0,b1,b2,b3,fac,dotir,dotjr,dotij,fact _REAL_ dphii_dx,dphii_dy,dphii_dz, & dphij_dx,dphij_dy,dphij_dz _REAL_ term,term0,term1,termi,termj,cgj _REAL_ filter_cut2,xx _REAL_ xktran(3,18) integer mask27 #ifdef LES # include "les.h" #endif # include "ew_nbe.h" # include "ew_tran.h" # include "constants.h" integer itran fac = two*ewaldcof*ewaldcof ee_vir_iso = zero del = one / eedtbdns dumx = zero dumy = zero dumz = zero edx = zero edy = zero edz = zero filter_cut2 = filter_cut*filter_cut cgi = charge(i) iaci = ntypes * (iac(i) - 1) mask27 = 2**27 - 1 do m=1,18 xktran(1,m) = tranvec(1,m) - xk xktran(2,m) = tranvec(2,m) - yk xktran(3,m) = tranvec(3,m) - zk end do #ifdef LES lestmp=nlesty*(lestyp(i)-1) #endif do m = 1,numvdw n = ipairs(m) itran=ishft(n,-27) n = iand(n,mask27) j = bckptr(n) delx = imagcrds(1,n) + xktran(1,itran) dely = imagcrds(2,n) + xktran(2,itran) delz = imagcrds(3,n) + xktran(3,itran) delr2 = delx*delx + dely*dely+delz*delz if ( delr2 < filter_cut2 )then delr = sqrt(delr2) delr2inv = one/delr2 x = dxdr*delr cgj = charge(j) if ( eedmeth == 1 )then ! -- cubic spline on switch ind = eedtbdns*x + 1 dx = x - (ind-one)*del switch = eed_cub(1,ind)+dx*(eed_cub(2,ind)+ & dx*(eed_cub(3,ind)+dx*eed_cub(4,ind)*third)*half) d_switch_dx = eed_cub(2,ind)+dx*(eed_cub(3,ind)+ & dx*eed_cub(4,ind)*half) else if ( eedmeth == 2 )then ! ---linear lookup on switch, deriv xx = eedtbdns*x + 1 ind = xx dx = xx - ind switch = (one - dx)*eed_lin(1,ind) + & dx*eed_lin(1,ind+1) d_switch_dx = (one - dx)*eed_lin(2,ind) + & dx*eed_lin(2,ind+1) else if ( eedmeth == 3 )then ! ---direct function call: call get_ee_func(x,switch,d_switch_dx,ee_type) else if ( eedmeth == 4 ) then ! ---use un-modified Coulomb interaction, no switch switch = one d_switch_dx = zero else write(6,*) 'bad eedmeth in ew_short_dip: ',eedmeth call mexit( 6,1 ) end if ! ( eedmeth == 1 ) ! TD Got the idea for B_l from Walter Smith's CCP5 article 1982 ! Ewald for point multipoles ! B_l satisfies grad_i B_l(|r_j - r_i|) = (r_j - r_i)B_{l+1}(|r_j-r_i|) ! grad_j B_l(|r_j - r_i|) = -grad_i B_l(|r_j - r_i|) b0 = switch*delr*delr2inv fact = d_switch_dx*dxdr b1 = (b0 - fact)*delr2inv fact = fac*fact b2 = (three*b1 - fact)*delr2inv fact = fac*fact b3 = (five*b2 - fact)*delr2inv ! B1 = (B0 - d_switch_dx*dxdr)*delr2inv ! B2 = (Three*B1 - fac*ewaldcof*d_switch_dx)*delr2inv ! B3 = (Five*B2 - fac*fac*ewaldcof*d_switch_dx)*delr2inv ! epol = dip_i dot grad_i of phii ! phii is direct sum electrostatic potential at i due to j ! so phii = cgj*B0 + dipj dot gradj of B0 = cgj*B0 - dotjr*B1 ! dphii_dx etc are derivatives with respect to r_i ! phij is direct sum electrostatic potential at j due to i ! dphij_dx etc are derivatives with respect to r_j dotjr = dipole(1,j)*delx+dipole(2,j)*dely+dipole(3,j)*delz dotir = dipole(1,i)*delx+dipole(2,i)*dely+dipole(3,i)*delz dotij = dipole(1,i)*dipole(1,j)+dipole(2,i)*dipole(2,j)+ & dipole(3,i)*dipole(3,j) ! gradi phii = cgj*rij*B1 + dipj*B1 - dotjr*rij*B2 ! so epol = -cgi*dotjr*B1 + (cgj*B1 - dotjr*B2)*dotir + dotij*B1 eelt = eelt + cgi*cgj*b0 term = cgj*dotir-cgi*dotjr+dotij epol = epol + term*b1 - dotir*dotjr*b2 term0 = cgi*cgj*b0 + term*b1 - dotir*dotjr*b2 ! so ene = ene + term0; dfx = dterm0_dx etc ! grad term0 = term1*rij + B1*grad_i term - grad_i dotir*dotjr*B2 ! grad_i term = -cgj*dip_i + cgi*dip_j ! grad_i dotir = -dip_i; similar for dotjr ! grad_i term0 = term1*rij + (-cgj*B1+dotjr*B2)*dip_i + ! (cgi*B1+dotir*B2)*dip_j term1 = cgi*cgj*b1 + term*b2 - dotir*dotjr*b3 termi = cgi*b1+dotir*b2 termj = cgj*b1-dotjr*b2 dfx = term1*delx + termi*dipole(1,j) - termj*dipole(1,i) dfy = term1*dely + termi*dipole(2,j) - termj*dipole(2,i) dfz = term1*delz + termi*dipole(3,j) - termj*dipole(3,i) ee_vir_iso = ee_vir_iso - dfx*delx - dfy*dely - dfz*delz ic = ico(iaci+iac(j)) r6 = delr2inv*delr2inv*delr2inv #ifdef LES lfac=lesfac(lestmp+lestyp(j)) f6 = cn2(ic)*r6*lfac f12 = cn1(ic)*(r6*r6)*lfac #else f6 = cn2(ic)*r6 f12 = cn1(ic)*(r6*r6) #endif evdw = evdw + f12 - f6 ! ---force related quantities df = (twelve*f12 - six*f6)*delr2inv dfx = dfx + df*delx dfy = dfy + df*dely dfz = dfz + df*delz virial(1,1) = virial(1,1) - dfx*delx virial(1,2) = virial(1,2) - dfx*dely virial(1,3) = virial(1,3) - dfx*delz virial(2,1) = virial(2,1) - dfy*delx virial(2,2) = virial(2,2) - dfy*dely virial(2,3) = virial(2,3) - dfy*delz virial(3,1) = virial(3,1) - dfz*delx virial(3,2) = virial(3,2) - dfz*dely virial(3,3) = virial(3,3) - dfz*delz frc(1,j) = frc(1,j) + dfx frc(2,j) = frc(2,j) + dfy frc(3,j) = frc(3,j) + dfz dumx = dumx + dfx dumy = dumy + dfy dumz = dumz + dfz ! ---field related quantities dphii_dx = termj*delx + b1*dipole(1,j) dphii_dy = termj*dely + b1*dipole(2,j) dphii_dz = termj*delz + b1*dipole(3,j) dphij_dx = -termi*delx + b1*dipole(1,i) dphij_dy = -termi*dely + b1*dipole(2,i) dphij_dz = -termi*delz + b1*dipole(3,i) edx = edx + dphii_dx edy = edy + dphii_dy edz = edz + dphii_dz field(1,j) = field(1,j) - dphij_dx field(2,j) = field(2,j) - dphij_dy field(3,j) = field(3,j) - dphij_dz end if ! ( delr2 < filter_cut2 ) end do ! m = 1,numvdw do m = numvdw+1,numtot n = ipairs(m) itran=ishft(n,-27) n = iand(n,mask27) j = bckptr(n) delx = imagcrds(1,n) + xktran(1,itran) dely = imagcrds(2,n) + xktran(2,itran) delz = imagcrds(3,n) + xktran(3,itran) delr2 = delx*delx + dely*dely+delz*delz if ( delr2 < filter_cut2 )then delr = sqrt(delr2) delr2inv = one/delr2 x = dxdr*delr cgj = charge(j) if ( eedmeth == 1 )then ! -- cubic spline on switch ind = eedtbdns*x + 1 dx = x - (ind-one)*del switch = eed_cub(1,ind)+dx*(eed_cub(2,ind)+ & dx*(eed_cub(3,ind)+dx*eed_cub(4,ind)*third)*half) d_switch_dx = eed_cub(2,ind)+dx*(eed_cub(3,ind)+ & dx*eed_cub(4,ind)*half) else if ( eedmeth == 2 )then ! ---linear lookup on switch, deriv xx = eedtbdns*x + 1 ind = xx dx = xx - ind switch = (one - dx)*eed_lin(1,ind) + & dx*eed_lin(1,ind+1) d_switch_dx = (one - dx)*eed_lin(2,ind) + & dx*eed_lin(2,ind+1) else if ( eedmeth == 3 )then ! ---direct function call: call get_ee_func(x,switch,d_switch_dx,ee_type) else if ( eedmeth == 4 ) then ! ---use un-modified Coulomb interaction, no switch switch = one d_switch_dx = zero else write(6,*) 'bad eedmeth in ew_short_dip: ',eedmeth call mexit( 6,1 ) end if ! ( eedmeth == 1 ) b0 = switch*delr*delr2inv fact = d_switch_dx*dxdr b1 = (b0 - fact)*delr2inv fact = fac*fact b2 = (three*b1 - fact)*delr2inv fact = fac*fact b3 = (five*b2 - fact)*delr2inv ! B0 = switch*delr*delr2inv ! B1 = (B0 - d_switch_dx*dxdr)*delr2inv ! B2 = (Three*B1 - fac*ewaldcof*d_switch_dx)*delr2inv ! B3 = (Five*B2 - fac*fac*ewaldcof*d_switch_dx)*delr2inv dotjr = dipole(1,j)*delx+dipole(2,j)*dely+dipole(3,j)*delz dotir = dipole(1,i)*delx+dipole(2,i)*dely+dipole(3,i)*delz dotij = dipole(1,i)*dipole(1,j)+dipole(2,i)*dipole(2,j)+ & dipole(3,i)*dipole(3,j) eelt = eelt + cgi*cgj*b0 term = cgj*dotir-cgi*dotjr+dotij epol = epol + term*b1 - dotir*dotjr*b2 term0 = cgi*cgj*b0 + term*b1 - dotir*dotjr*b2 term1 = cgi*cgj*b1 + term*b2 - dotir*dotjr*b3 termi = cgi*b1+dotir*b2 termj = cgj*b1-dotjr*b2 dfx = (term1)*delx + termi*dipole(1,j) - termj*dipole(1,i) dfy = (term1)*dely + termi*dipole(2,j) - termj*dipole(2,i) dfz = (term1)*delz + termi*dipole(3,j) - termj*dipole(3,i) ee_vir_iso = ee_vir_iso - dfx*delx - dfy*dely - dfz*delz #ifdef HAS_10_12 ! --- this code allows 10-12 terms; in many (most?) (all?) cases, the ! only "nominal" 10-12 terms are on waters, where the asol and bsol ! parameters are always zero; hence we can skip the L-J part; note ! that we still have to compute the electrostatic interactions ic = -ico(iaci+iac(j)) r10 = delr2inv*delr2inv*delr2inv*delr2inv*delr2inv f10 = bsol(ic)*r10 f12 = asol(ic)*(r10*delr2inv) ehb = ehb + f12 - f10 df = (twelve*f12 - ten*f10)*delr2inv #else df = zero #endif ! ---force related quantities dfx = dfx + df*delx dfy = dfy + df*dely dfz = dfz + df*delz virial(1,1) = virial(1,1) - dfx*delx virial(1,2) = virial(1,2) - dfx*dely virial(1,3) = virial(1,3) - dfx*delz virial(2,1) = virial(2,1) - dfy*delx virial(2,2) = virial(2,2) - dfy*dely virial(2,3) = virial(2,3) - dfy*delz virial(3,1) = virial(3,1) - dfz*delx virial(3,2) = virial(3,2) - dfz*dely virial(3,3) = virial(3,3) - dfz*delz frc(1,j) = frc(1,j) + dfx frc(2,j) = frc(2,j) + dfy frc(3,j) = frc(3,j) + dfz dumx = dumx + dfx dumy = dumy + dfy dumz = dumz + dfz ! ---field related quantities dphii_dx = termj*delx + b1*dipole(1,j) dphii_dy = termj*dely + b1*dipole(2,j) dphii_dz = termj*delz + b1*dipole(3,j) dphij_dx = -termi*delx + b1*dipole(1,i) dphij_dy = -termi*dely + b1*dipole(2,i) dphij_dz = -termi*delz + b1*dipole(3,i) edx = edx + dphii_dx edy = edy + dphii_dy edz = edz + dphii_dz field(1,j) = field(1,j) - dphij_dx field(2,j) = field(2,j) - dphij_dy field(3,j) = field(3,j) - dphij_dz end if ! ( delr2 < filter_cut2 ) end do ! m = numvdw+1,numtot frc(1,i) = frc(1,i) - dumx frc(2,i) = frc(2,i) - dumy frc(3,i) = frc(3,i) - dumz field(1,i) = field(1,i) - edx field(2,i) = field(2,i) - edy field(3,i) = field(3,i) - edz eedvir = eedvir + ee_vir_iso return end subroutine short_ene_dip !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Re-balance the grid by divvying cells among processors. ! !----------------------------------------------------------------------- ! --- FIX_GRID_BALANCE --- !----------------------------------------------------------------------- subroutine fix_grid_balance(myindexlo,myindexhi,gridpairs, & nucgrd,numtasks,mytaskid,listdiffmax) implicit none integer myindexlo,myindexhi,gridpairs(*) integer nucgrd,listtot,numtasks,mytaskid,listdiffmax integer i,j,n,lsum,i0,ishare n=0 lsum=0 myindexlo=0 myindexhi=0 listtot=0 do i=1,nucgrd listtot=listtot+gridpairs(i) end do i0=1 ishare=listtot/numtasks ! ------ set trigger for new balance as 10 percent of list size listdiffmax=max(1000,ishare/10) ! ------ give cells to each processor .lt.mytaskid till ! they have ishare, then give the next ! cells to this pe till it has its share, then return do i=1,nucgrd lsum=lsum+gridpairs(i) if(lsum >= ishare)then if(n == mytaskid)then myindexlo=i0 myindexhi=i return end if lsum=0 i0=i+1 n=n+1 end if end do if(n == mytaskid)then myindexlo=i0 myindexhi=nucgrd return end if myindexlo=nucgrd return end subroutine fix_grid_balance