C C $Id: setncell.F,v 1.6 1998/07/16 16:40:42 jjv5 Exp arjan $ C C------------------------------------------------------------------------ subroutine setncell(ierror) C sets parameters for grid-subsetting. C called everytime a gridsubsetting is done, C because the minimum and maximum coordinates C can change in a NPT-MC-simulation inducing a C change in the grid-subsetting parameters. C C conditions for these parameters: C C 1. the core must be build from 2 or more grid-cells C . (in 1 dimension) in order to maximize the use of C . symmetry (a lot of if-statements can be avoided in C . the subsetting in this way) C 2. the cell-length should not exceed the value of (core-overlap) C . (resolution error) C 3. the total number of cells shouldn't exceed the preset C . number in divcon.dim. An adaptive correction scheme C . is used to prevent this error. C C if dgbuff1 or dgbuff2 is smaller than rcell, more atoms C might be included in the buffer than would be expected (the actual C buffer region is <= rcell). A warning will be given in this C case. C C written by Arjan van der Vaart, August 1997 implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" parameter (eps=1.0E-6) C locals character dir dimension icore(3), idgbuff1(3), idgbuff2(3), igt(3), dir(3) data dir /'X','Y','Z'/ save dir CC---------------------------------------------------------CC if (resgr.and.pbc) then C . since geometric centers are confined in box and C . individual atoms are allowed outside of box, add C . gyration radius to boxlength for subsetting do i=1,3 gdxyz(i) = dbox(i)+2.0*gyr gxyzmin(i) = gcmin(i)-gyr gxyzmax(i) = gcmax(i)+gyr enddo elseif (pbc) then C . do subsetting over the entire box in order to include C . the proper images in the buffers later on do i=1,3 gdxyz(i) = dbox(i) gxyzmin(i) = boxmin(i) gxyzmax(i) = boxmax(i) enddo else C . do subsetting over the range of coordinates (which is C . usually smaller than the box) do i=1,3 gdxyz(i) = dxyz(i) gxyzmin(i) = xyzmin(i) gxyzmax(i) = xyzmax(i) enddo endif ldgbuff1 = dgbuff1.lt.eps ldgbuff2 = dgbuff2.lt.eps iigt = 0 do 10 i=1,3 if (ldgbuff1.and.ldgbuff2) then rcell(i) = min(0.5*core(i), (core(i)-overlap)) elseif (ldgbuff1) then rcell(i) = min(0.5*core(i), (core(i)-overlap), & dgbuff2) elseif (ldgbuff2) then rcell(i) = min(0.5*core(i), (core(i)-overlap), & dgbuff1) else rcell(i) = min(0.5*core(i), (core(i)-overlap), & dgbuff1, dgbuff2) endif ncell(i) = gdxyz(i)/rcell(i) d = gdxyz(i) - ncell(i)*rcell(i) if (d.gt.eps) ncell(i) = ncell(i) + 1 if (ncell(i).gt.mcell) then iigt = iigt + 1 igt(iigt) = i endif 10 continue C adaptive correction scheme n3o = ncell(1)*ncell(2)*ncell(3) ncell3 = n3o if (n3o.gt.mcell3) then x3 = n3o**0.3333333333333 n3o = x3 + 1 if (iigt.eq.3) then C-RDC write(iout,'(/" ERROR IN GRID-BASED SUBSETTING:", C-RDC & /" INCREASE MCELL IN divcon.dim TO ",I5)') n3o ierror = 1 return elseif (iigt.eq.2) then if (ncell(igt(1)).lt.ncell(igt(2))) then ncell(igt(1)) = mcell rcell(igt(1)) = gdxyz(igt(1))/ncell(igt(1)) else ncell(igt(2)) = mcell rcell(igt(2)) = gdxyz(igt(2))/ncell(igt(2)) endif ncell3 = ncell(1)*ncell(2)*ncell(3) if (ncell3.gt.mcell3) then ncell(igt(1)) = mcell rcell(igt(1)) = gdxyz(igt(1))/ncell(igt(1)) ncell(igt(2)) = mcell rcell(igt(2)) = gdxyz(igt(2))/ncell(igt(2)) ncell3 = ncell(1)*ncell(2)*ncell(3) if (ncell3.gt.mcell3) then C-RDC write(iout,'(/" ERROR IN GRID-BASED SUBSETTING:", C-RDC & /" INCREASE MCELL IN divcon.dim TO ",I5)') n3o ierror = 1 return endif endif else ncell(igt(1)) = mcell rcell(igt(1)) = gdxyz(igt(1))/ncell(igt(1)) ncell3 = ncell(1)*ncell(2)*ncell(3) if (ncell3.gt.mcell3) then C-RDC write(iout,'(/" ERROR IN GRID-BASED SUBSETTING:", C-RDC & /" INCREASE MCELL IN divcon.dim TO ",I5)') n3o ierror = 1 return endif endif endif ncell2 = ncell(1)*ncell(2) do 20 i=1,3 ncellm1(i) = ncell(i)-1 20 continue C check stuff do 30 i=1,3 icore(i) = core(i) / rcell(i) d = core(i) - icore(i)*rcell(i) if (d.gt.eps) icore(i) = icore(i) + 1 if (icore(i).lt.2) then C-RDC write(iout,'(/" ERROR IN GRID-BASED SUBSETTING:", C-RDC & /" CORE RESOLUTION IS TOO LOW. INCREASE CORE OR", C-RDC & " MCELL IN divcon.dim")') ierror = 1 return endif gstep(i) = core(i) - overlap igstep(i) = gstep(i)/rcell(i) if (igstep(i).lt.1) then C-RDC write(iout,'(/" ERROR IN GRID-BASED SUBSETTING:", C-RDC & /" CORE-OVERLAP RESOLUTION IS TOO LOW. DECREASE", C-RDC & " OVERLAP OR INCREASE CORE OR", C-RDC & /" MCELL IN divcon.dim")') ierror = 1 return endif C . note that ngstep = ngstep - 1 ngstep(i) = (ncell(i)-1)/igstep(i) if (ldgbuff1) then idgbuff1(i) = 0 else idgbuff1(i) = dgbuff1/rcell(i) d = dgbuff1 - idgbuff1(i)*rcell(i) if (d.gt.eps) idgbuff1(i) = idgbuff1(i) + 1 if (dgbuff1.lt.rcell(i)) then C-RDC write(iout,'(/" WARNING IN GRID-SUBSETTING:", C-RDC & /" THE GRIDCELL-LENGTH IN THE ",1A, C-RDC & " DIRECTION EXCEEDS DBUFF1.", C-RDC & /" MORE ATOMS MIGHT BE INCLUDED IN THE INNER", C-RDC & " BUFFER THAN EXPECTED.")') dir(i) endif endif if (ldgbuff2) then idgbuff2(i) = 0 else idgbuff2(i) = dgbuff2/rcell(i) d = dgbuff2 - idgbuff2(i)*rcell(i) if(d.gt.eps) idgbuff2(i) = idgbuff2(i) + 1 if (dgbuff2.lt.rcell(i)) then C-RDC write(iout,'(/" WARNING IN GRID-SUBSETTING:", C-RDC & /" THE GRIDCELL-LENGTH IN THE ",1A, C-RDC & " DIRECTION EXCEEDS DBUFF2.", C-RDC & /" MORE ATOMS MIGHT BE INCLUDED IN THE OUTER", C-RDC & " BUFFER THAN EXPECTED.")') dir(i) endif endif 30 continue ldgbuff1 = .not.ldgbuff1 ldgbuff2 = .not.ldgbuff2 do 40 i=1,3 iextl(i,1) = 0 iextl(i,2) = -idgbuff1(i) iextl(i,3) = iextl(i,2) - idgbuff2(i) iextr(i,1) = icore(i) - 1 iextr(i,2) = icore(i) + 2*idgbuff1(i) - 1 iextr(i,3) = iextr(i,2) + 2*idgbuff2(i) extl(i,1) = gxyzmin(i) extl(i,2) = gxyzmin(i) - dgbuff1 extl(i,3) = extl(i,2) - dgbuff2 extr(i,1) = core(i) extr(i,2) = core(i) + 2.0*dgbuff1 extr(i,3) = extr(i,2) + 2.0*dgbuff2 40 continue end