C C $Id: mixgrid.F,v 1.9 1998/07/16 16:40:05 jjv5 Exp arjan $ C C------------------------------------------------------------------------ subroutine mixgrid(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 residue-based for core and atom based for C 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: since this is a mixed subsetting C scheme both atom-wise and residue-wise grids are needed call mkagrid call mkrgrid C initialize some stuff before subsetting maxatsub = 0 imaxatsub = 0 n = 0 ngsub = 0 ngatot = 0 C-RDC if (prtsub) write(iout,'(//" SUBSYSTEM DEFINITIONS:"/)') C do the actual subsetting do 130 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 120 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 110 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 residue-based call getrcorners(1,k1,igsub1,k2,igsub2) call getrribs(1,k1,igsub1,k2,igsub2) call getrplanes(1,k1,igsub1,k2,igsub2) call getrcore(k1,igsub1) if (combsub.and.(k1.ne.0)) then call cmbcore(k1,igsub1,k2,igsub2) endif if (k1.ne.0) then if (ldgbuff1) then do i=1,natoms nottaken(i) = .true. enddo do i=1,k1 nottaken(igsub1(i)) = .false. enddo C . buffers are atom-based if (pbc) then 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) 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 C . because part of the inner buffer was generated C . by the rejected atoms of the core (which is C . residue based), the same atom might occur more C . than once in igsub2: so make igsub2 unique. C . make sure none of the atoms of the inner buffer C . are already included in the core. call busort(k2, igsub2, nottaken) if (ldgbuff2) then do i=1,k2 nottaken(igsub2(i)) = .false. enddo if (pbc) then 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) 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 C . make sure none of the outer buffer atoms are C . already included in the core or inner buffer call busort(k3,igsub3,nottaken) endif else k2 = 0 endif nsub = nsub + 1 ngsub = ngsub + 1 if (nsub .gt. maxsub) then C-RDC write(iout,'(" MAXIMUM NUMBER OF ", C-RDC & "SUBSYSTEMS (",I5,") REACHED." C-RDC & /" CHANGE MAXSUB IN divcon.dim AND ", C-RDC & "RECOMPILE.")') maxsub ierror = 1 return endif n = k1 + k2 + k3 natot = natot + n ngatot = ngatot + n if (natot .gt. mslist) then C-RDC write(iout,'(/" ATOM LIST STORAGE LIMIT ", C-RDC & "REACHED.",/" INCREASE MSLIST PARAMETER ", C-RDC & "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 80 iat=1,k1 iatoms(ipoint) = igsub1(iat) iabuff(ipoint) = 0 ipoint = ipoint + 1 80 continue do 90 iat=1,k2 iatoms(ipoint) = igsub2(iat) iabuff(ipoint) = 1 ipoint = ipoint + 1 90 continue do 100 iat=1,k3 iatoms(ipoint) = igsub3(iat) iabuff(ipoint) = 2 ipoint = ipoint + 1 100 continue C . sort atom list and carry along buffer status istart = iatom1(nsub) call bsort(n,iatoms(istart),iabuff(istart)) endif 110 continue 120 continue 130 continue iatom1(nsub+1) = ipoint natavg = ngatot/ngsub C-RDC write(iout,'(//" SUMMARY OF RESIDUE WISE FOR CORE, ", C-RDC & "ATOM WISE FOR BUFFERS,", C-RDC & /" 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,'(/"mixgrid 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