C C $Id: atgrid.F,v 1.9 1998/07/16 16:39:19 jjv5 Exp arjan $ C C------------------------------------------------------------------------ subroutine atgrid(ipoint, natot, nbig, ibig, ierror) C C grid based subsetting C C written by Arjan van der Vaart, July-August 1997 C . January 1998 (pbc) C C C subsetting is atom-based for core and buffers. C implicit double precision (a-h,o-z) #include "divcon.dim" #include "divcon.h" C locals: logical nottaken dimension igsub1(maxatm), igsub2(maxatm), igsub3(maxatm), & igsub4(maxatm), nottaken(maxatm) CC-------------------------------------------------------------------CC C initialize ncell etc. This needs to be done everytime C this routine is called because it is dependent on the C maximum and minimum coordinates (can change during a NPT-MC-run) call setncell(ierror) if (ierror.ne.0) return C set up the grids: only atom-based is needed call mkagrid C initialize some stuff before subsetting maxatsub = 0 imaxatsub = 0 n = 0 ngsub = 0 ngatot = 0 if (prtsub) write(iout,'(" SUBSYSTEM DEFINITIONS:")') C do the actual subsetting do 90 iz=0,ngstep(3) do 10 i=1,3 ibg(3,i) = iextl(3,i) + igstep(3)*iz ied(3,i) = ibg(3,i) + iextr(3,i) gridmin(3,i) = iz*gstep(3) + extl(3,i) gridmax(3,i) = gridmin(3,i) + extr(3,i) 10 continue do 80 iy=0,ngstep(2) do 20 i=1,3 ibg(2,i) = iextl(2,i) + igstep(2)*iy ied(2,i) = ibg(2,i) + iextr(2,i) gridmin(2,i) = iy*gstep(2) + extl(2,i) gridmax(2,i) = gridmin(2,i) + extr(2,i) 20 continue do 70 ix=0,ngstep(1) do 30 i=1,3 ibg(1,i) = iextl(1,i) + igstep(1)*ix ied(1,i) = ibg(1,i) + iextr(1,i) gridmin(1,i) = ix*gstep(1) + extl(1,i) gridmax(1,i) = gridmin(1,i) + extr(1,i) 30 continue k1 = 0 k2 = 0 k3 = 0 k4 = 0 C . core is atom-based, based on *real* distances call getacorners(1,k1,igsub1,k2,igsub2) call getaribs(1,k1,igsub1,k2,igsub2) call getaplanes(1,k1,igsub1,k2,igsub2) call getacore(k1,igsub1) if (combsub.and.(k1.ne.0)) then call cmbcore(k1,igsub1,k2,igsub2) endif if (k1.ne.0) then if (ldgbuff1) then C . buffers are atom-based if (pbc) then do i=1,natoms nottaken(i) = .true. enddo do i=1,k1 nottaken(igsub1(i)) = .false. enddo call pgetacorners(2, k2, igsub2, k3, igsub3) call pgetaribs(2, k2, igsub2, k3, igsub3) call pgetaplanes(2, k2, igsub2, k3, igsub3) call pgetainner(2, k2, igsub2) call busort(k2,igsub2,nottaken) else call getacorners(2, k2, igsub2, k3, igsub3) call getaribs(2, k2, igsub2, k3, igsub3) call getaplanes(2, k2, igsub2, k3, igsub3) call getainner(2, k2, igsub2) endif if (ldgbuff2) then if (pbc) then do i=1,k2 nottaken(igsub2(i)) = .false. enddo call pgetacorners(3, k3, igsub3, k4, igsub4) call pgetaribs(3, k3, igsub3, k4, igsub4) call pgetaplanes(3, k3, igsub3, k4, igsub4) call pgetainner(3, k3, igsub3) call busort(k3,igsub3,nottaken) else call getacorners(3, k3, igsub3, k4, igsub4) call getaribs(3, k3, igsub3, k4, igsub4) call getaplanes(3, k3, igsub3, k4, igsub4) call getainner(3, k3, igsub3) endif else k3 = 0 endif else k2 = 0 endif nsub = nsub + 1 ngsub = ngsub + 1 if (nsub .gt. maxsub) then write(iout,'(" MAXIMUM NUMBER OF ", & "SUBSYSTEMS (",I5,") REACHED.", & /," CHANGE MAXSUB IN divcon.dim AND ", & "RECOMPILE.")') maxsub ierror = 1 return endif n = k1 + k2 + k3 natot = natot + n ngatot = ngatot + n if (natot .gt. mslist) then write(iout,'(/" ATOM LIST STORAGE LIMIT ", & "REACHED.",/" INCREASE MSLIST PARAMETER ", & "IN divcon.dim AND RECOMPILE")') ierror = 1 return endif if (n.gt.maxatsub) then maxatsub = n imaxatsub = nsub if (n.gt.nbig) then nbig = n ibig = nsub endif endif if (prtsub) call printsub(nsub, k1, igsub1, & k2, igsub2, k3, igsub3) iatom1(nsub) = ipoint do 40 iat=1,k1 iatoms(ipoint) = igsub1(iat) iabuff(ipoint) = 0 ipoint = ipoint + 1 40 continue do 50 iat=1,k2 iatoms(ipoint) = igsub2(iat) iabuff(ipoint) = 1 ipoint = ipoint + 1 50 continue do 60 iat=1,k3 iatoms(ipoint) = igsub3(iat) iabuff(ipoint) = 2 ipoint = ipoint + 1 60 continue C . sort atom list and carry along buffer status istart = iatom1(nsub) call bsort(n,iatoms(istart),iabuff(istart)) endif 70 continue 80 continue 90 continue iatom1(nsub+1) = ipoint natavg = ngatot/ngsub C-RDC write(iout,'(//" SUMMARY OF ATOM WISE, GRID BASED SUBSETTING:", C-RDC & /" -------------------------------------------", C-RDC & //" EFFECTIVE BOXDIMENSIONS FOR SUBSETTING:", C-RDC & /" X = ",F10.5," (",F10.5," - ",F10.5,")", C-RDC & /" Y = ",F10.5," (",F10.5," - ",F10.5,")", C-RDC & /" Z = ",F10.5," (",F10.5," - ",F10.5,")")') C-RDC & gdxyz(1),gxyzmin(1),gxyzmax(1), C-RDC & gdxyz(2),gxyzmin(2),gxyzmax(2), C-RDC & gdxyz(3),gxyzmin(3),gxyzmax(3) C-RDC write(iout,'(/" NUMBER OF SUBSYSTEMS = ",I4, C-RDC & /" AVERAGE SUBSYSTEM SIZE = ",I4, C-RDC & /" MAXIMUM SUBSYSTEM SIZE = ",I4, C-RDC & /" LARGEST SUBSYSTEM = #",I4/)') C-RDC & ngsub,natavg,maxatsub,imaxatsub if (screen) then C-RDC write(iscr,'(/"atgrid subsetting: nsub = ",i4,", av. size = ", C-RDC & i4," max. size = ",i4)') ngsub,natavg,maxatsub endif C in case eigenvectors will be written if (prtvec) then ncores = 1 isubend(0) = 0 isubend(1) = nsub endif end