*******> update.10 Author: Charles Lin and Ido Ben-Shalom Date: Jan, 3rd 2019 Programs: pmemd.cuda and pmemd Description: This adds a new implementation of the mcwat code that is speed optimized to take advantage of using steric grids and coarse energy calculations to avoid spending expensive computation time on intrinsically obvious unfavorable states. ------------------------------------------------------------------------------- src/pmemd/src/mcres.F90 | 1295 ++++++++++++++++++++++------------ src/pmemd/src/mdin_ctrl_dat.F90 | 24 +- src/pmemd/src/runmd.F90 | 21 +- test/cuda/mcwat/mdout.mcwat.GPU_DPFP | 182 +++-- test/cuda/mcwat/mdout.mcwat.GPU_SPFP | 180 +++-- test/cuda/mcwat/Run.mcwat | 15 +- 3 files changed, 1032 insertions(+), 785 deletions(-) diff --git src/pmemd/src/mcres.F90 src/pmemd/src/mcres.F90 index 510ed4e..ad4ab7f 100644 --- src/pmemd/src/mcres.F90 +++ src/pmemd/src/mcres.F90 @@ -36,20 +36,37 @@ module mcres_mod ! contains - real, parameter :: GSP = 0.2 ! GRID SPACING - real, parameter :: WATERRADIUS = 0.9 - double precision :: SHIFT + real, parameter :: GSP = 0.2 ! GRID SPACING + real, parameter :: WATERRADIUS = 0.9 + double precision :: SHIFT + double precision :: mnX,mnY,mnZ, mxX,mxY,mxZ, X,Y,Z, r + integer :: atom,delta,cnt,numberofwaters + integer :: error,stat, mxXVxl,mxYVxl,mxZVxl,nmEmpVxls,rndmWtr - double precision :: mnX,mnY,mnZ, mxX,mxY,mxZ, X,Y,Z, r - integer :: atom,reject,delta,cnt - integer :: error,stat, mxXVxl,mxYVxl,mxZVxl,nmEmpVxls,rndmWtr + integer, allocatable :: g(:,:,:) - integer, allocatable :: g(:,:,:) - integer :: rndmVxl(4) - double precision :: tempX, tempY, tempZ + integer, allocatable :: listEmptyVoxels(:,:) - logical :: empty + integer, allocatable :: waterList(:) + + integer :: rndmVxl(3) + double precision :: tempX, tempY, tempZ + double precision :: deltacrd(3) + double precision :: old_crd(3,3) + double precision :: new_crd(3,3) + logical :: empty,revert + integer :: numberOfAtoms + + integer :: maxCoarseXVxl, maxCoarseYVxl, maxCoarseZVxl + double precision :: coarseGSP + + integer, parameter :: UNIT = 30 + integer, allocatable :: coarseGrid(:,:,:,:) + + logical :: up + + character(4) :: mol_str contains @@ -91,7 +108,9 @@ subroutine setup_mcres(mol_str) ! Can be optimized later by only being used on first call. mc_err = .true. mcres_numatms = 0 - SHIFT = mcboxshift + SHIFT = max(GSP,mcboxshift) + + mcrescyc = nmc do i= 1, nres if(gbl_labres(i) .eq. mol_str) then @@ -100,6 +119,9 @@ subroutine setup_mcres(mol_str) end if end do + numberOfAtoms = mcres_numatms + +!print *,'numberOfAtoms = ',numberOfAtoms #ifdef MPI ! For legacy reasons we're going to do this backwards where false = 1 and ! true = 0 (to not break other versions of this code) @@ -138,47 +160,13 @@ end subroutine setup_mcres !******************************************************************************* ! -! Subroutine: wrap_for_MC -! -! Description: This function wraps the atomic coordinates. -! -!******************************************************************************* - -subroutine wrap_for_MC(atm_cnt, crd, crd_wrap) - - use mdin_ctrl_dat_mod - use pbc_mod - use prmtop_dat_mod - - implicit none - - integer, intent(in) :: atm_cnt - double precision, intent(inout) :: crd(3, atm_cnt) - double precision :: crd_wrap(3,atm_cnt) - -! local variables - double precision :: crd_copy(3, atm_cnt) - - crd_copy(:,:) = crd(:,:) - - if (ntb .gt. 0) then - call wrap_molecules(crd_copy) - if (ifbox .eq. 2) call wrap_to(crd_copy) - end if - - crd_wrap(:,:) = crd_copy(:,:) - return -end subroutine wrap_for_MC - -!******************************************************************************* -! ! Subroutine: findMnMxnPBox ! ! Description: This function finds the minimum and maximum of X, Y, and Z of the periodic box. ! !******************************************************************************* -subroutine findMnMxnPBox(crd,atm_cnt,mnX,mxX,mnY,mxY,mnZ,mxZ) +subroutine findMnMxnPBox(crd,atm_cnt,mxX,mxY,mxZ) use pbc_mod @@ -186,18 +174,13 @@ subroutine findMnMxnPBox(crd,atm_cnt,mnX,mxX,mnY,mxY,mnZ,mxZ) integer, intent(in) :: atm_cnt double precision, intent(in) :: crd(3, atm_cnt) - double precision, intent(out) :: mnX, mxX, mnY, mxY, mnZ, mxZ + double precision, intent(out) :: mxX, mxY, mxZ - mnX = 0 mxX = pbc_box(1) - mnY = 0 mxY = pbc_box(2) - mnZ = 0 mxZ = pbc_box(3) - - ! though the minimum is 0 for all axes, they are defined explicitly - - print *, 'mnX, mxX, mnY, mxY, mnZ, mxZ =', mnX, mxX, mnY, mxY, mnZ, mxZ + +! print *, 'mnX, mxX, mnY, mxY, mnZ, mxZ =', mnX, mxX, mnY, mxY, mnZ, mxZ end subroutine findMnMxnPBox @@ -209,21 +192,24 @@ end subroutine findMnMxnPBox ! !******************************************************************************* -subroutine findMxVxls(mnX,mxX,mnY,mxY,mnZ,mxZ,mxXVxl,mxYVxl,mxZVxl) +subroutine findMxVxls(mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl) + implicit none - double precision, intent(in) :: mnX,mxX,mnY,mxY,mnZ,mxZ - integer, intent(out) :: mxXVxl, mxYVxl, mxZVxl + double precision, intent(in) :: mxX,mxY,mxZ - mxXVxl = ceiling((mxX - mnX - SHIFT) / GSP) - mxYVxl = ceiling((mxY - mnY - SHIFT) / GSP) - mxZVxl = ceiling((mxZ - mnZ - 3) / GSP) ! exceeding the box into the solvent in one axis + integer, intent(out) :: mxXVxl, mxYVxl, mxZVxl -! mxZVxl = ceiling((mxZ - mnZ - SHIFT) / GSP) + mxXVxl = ceiling((mxX - SHIFT) / GSP)+1 + mxYVxl = ceiling((mxY - SHIFT) / GSP)+1 + mxZVxl = ceiling((mxZ - 3) / GSP)+1 ! exceeding the box into the solvent in one axis +! mxZVxl = ceiling((mxZ - SHIFT) / GSP)+1 ! print *, 'findMnMxVxls: mxXVxl, mxYVxl, mxZVxl = ', mxXVxl, mxYVxl, mxZVxl + end subroutine findMxVxls + !******************************************************************************* ! ! Subroutine: convertAtomicToGridUnits @@ -233,19 +219,47 @@ end subroutine findMxVxls ! !******************************************************************************* -subroutine convertAtomicToGridUnits(X,Y,Z,mnX,mnY,mnZ,gridX,gridY,gridZ,xVxl,yVxl,zVxl) +subroutine convertAtomicToGridUnits(X,Y,Z, mxX,mxY,mxZ, xVxl,yVxl,zVxl) + + use pbc_mod implicit none - double precision, intent(in) :: X,Y,Z, mnX,mnY,mnZ - double precision, intent(out) :: gridX,gridY,gridZ, xVxl,yVxl,zVxl + double precision, intent(in) :: mxX,mxY,mxZ + double precision, intent(inout) :: X,Y,Z + integer, intent(out) :: xVxl,yVxl,zVxl + + call wrapCoordinates(X,Y,Z,mxX,mxY,mxZ) - xVxl = ((X - mnX) / GSP) - yVxl = ((Y - mnY) / GSP) - zVxl = ((Z - mnZ) / GSP) + xVxl = ceiling(X/GSP) + yVxl = ceiling(Y/GSP) + zVxl = ceiling(Z/GSP) end subroutine convertAtomicToGridUnits + +subroutine wrapCoordinates(X,Y,Z,mxX,mxY,mxZ) + + use pbc_mod + + implicit none + + double precision, intent(inout) :: X,Y,Z + double precision, intent(in) :: mxX,mxY,mxZ + + double precision :: tempX,tempY,tempZ + + tempX = modulo(X,mxX) + tempY = modulo(Y,mxY) + tempZ = modulo(Z,mxZ) + + X = tempX + Y = tempY + Z = tempZ + +end subroutine wrapCoordinates + + !******************************************************************************* ! ! Subroutine: convertGridToAtomicUnits @@ -255,19 +269,20 @@ end subroutine convertAtomicToGridUnits ! !******************************************************************************* -subroutine convertGridToAtomicUnits(mnX,mnY,mnZ,xVxl,yVxl,zVxl,X,Y,Z) +subroutine convertGridToAtomicUnits(xVxl,yVxl,zVxl,X,Y,Z) implicit none - double precision, intent(in) :: mnX,mnY,mnZ,xVxl,yVxl,zVxl + double precision, intent(in) :: xVxl,yVxl,zVxl double precision, intent(out) :: X,Y,Z - X = ((xVxl*GSP) + mnX) - Y = ((yVxl*GSP) + mnY) - Z = ((zVxl*GSP) + mnZ) + X = (xVxl*GSP) + Y = (yVxl*GSP) + Z = (zVxl*GSP) end subroutine convertGridToAtomicUnits + !******************************************************************************* ! ! Subroutine: initializeGrid @@ -281,14 +296,14 @@ subroutine initializeGrid(g,mxXVxl,mxYVxl,mxZVxl) implicit none integer, intent(in) :: mxXVxl, mxYVxl, mxZVxl - integer, intent(inout) :: g(0:mxXVxl,0:mxYVxl,0:mxZVxl) + integer, intent(inout) :: g(mxXVxl,mxYVxl,mxZVxl) ! local variables integer :: l,m,n, mnXVxl,mnYVxl,mnZVxl - mnXVxl = (SHIFT / GSP) - mnYVxl = (SHIFT / GSP) - mnZVxl = (SHIFT / GSP) + mnXVxl = ceiling(SHIFT / GSP) + mnYVxl = ceiling(SHIFT / GSP) + mnZVxl = ceiling(SHIFT / GSP) do n = mnZVxl, mxZVxl do m = mnYVxl, mxYVxl @@ -297,10 +312,10 @@ subroutine initializeGrid(g,mxXVxl,mxYVxl,mxZVxl) end do end do end do -! print *,'END initializeGrid' end subroutine initializeGrid + !******************************************************************************* ! ! Subroutine: determineAtomRadius @@ -396,6 +411,7 @@ subroutine determineAtomRadius(crd,cnt,totalRadius,atm_cnt) end subroutine determineAtomRadius + !******************************************************************************* ! ! Subroutine: buildGrid @@ -407,33 +423,41 @@ end subroutine determineAtomRadius ! !******************************************************************************* -subroutine buildGrid(crd,g,atm_cnt,mnX,mnY,mnZ,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,cnt) +subroutine buildGrid(crd,g,atm_cnt,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,cnt) implicit none integer, intent(in) :: atm_cnt, mxXVxl,mxYVxl,mxZVxl double precision, intent(in) :: crd(3, atm_cnt) integer, intent(inout) :: g(:,:,:) - double precision, intent(in) :: mnX, mnY, mnZ, mxX, mxY, mxZ + double precision, intent(in) :: mxX, mxY, mxZ + integer, intent(out) :: delta,cnt !local variables integer :: counter,round, l + integer :: i,ic double precision :: X,Y,Z delta = 1 - do cnt = 1, atm_cnt - round = nint(atm_mass(cnt)) + do i = 1, atm_cnt + round = nint(atm_mass(i)) if (round > 3) then - X = crd(1, cnt) ! atomic units - Y = crd(2, cnt) - Z = crd(3, cnt) - call placeInGrid(crd,g,atm_cnt,mnX,mnY,mnZ,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,X,Y,Z,nmEmpVxls,cnt) + X = crd(1, i) ! atomic units + Y = crd(2, i) + Z = crd(3, i) + +! the wrapping is done inside "placeInGrid" + + call placeInGrid(crd,g,atm_cnt,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,X,Y,Z,nmEmpVxls,i) + ! at the buildGrid part, the variable nmEmpVxls is 0 !!! endif end do + ! print *,'END buildGrid' +!call flush() end subroutine buildGrid @@ -447,192 +471,151 @@ end subroutine buildGrid ! !******************************************************************************* -subroutine placeInGrid(crd,g,atm_cnt,mnX,mnY,mnZ,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,X,Y,Z,nmEmpVxls,cnt) +subroutine placeInGrid(crd,g,atm_cnt,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,X,Y,Z,nmEmpVxls,cnt) implicit none integer, intent(in) :: atm_cnt double precision, intent(in) :: crd(3, atm_cnt) integer, intent(inout) :: g(:,:,:) - double precision, intent(in) :: mnX,mnY,mnZ, mxX,mxY,mxZ, X,Y,Z + double precision, intent(in) :: mxX,mxY,mxZ + double precision, intent(inout) :: X,Y,Z integer, intent(in) :: mxXVxl, mxYVxl, mxZVxl, delta,cnt integer, intent(inout) :: nmEmpVxls ! local variables : integer :: p, round, l, m, n - double precision :: xVxl, yVxl, zVxl, r, totalRadius - double precision :: atmMnXVxl,atmMxXVxl,atmMnYVxl, atmMxYVxl,atmMnZVxl,atmMxZVxl, xRange,yRange,Zrange - double precision :: mnXVxl,mnYVxl,mnZVxl, gridX,gridY,gridZ + double precision :: r, totalRadius + integer :: xVxl, yVxl, zVxl + integer :: atmMnXVxl,atmMxXVxl,atmMnYVxl, atmMxYVxl,atmMnZVxl,atmMxZVxl, mnXVxl,mnYVxl,mnZVxl + double precision :: xRange,yRange,Zrange - mnXVxl = (SHIFT / GSP) - mnYVxl = (SHIFT / GSP) - mnZVxl = (SHIFT / GSP) + mnXVxl = ceiling(SHIFT / GSP) + mnYVxl = ceiling(SHIFT / GSP) + mnZVxl = ceiling(SHIFT / GSP) call determineAtomRadius(crd,cnt,totalRadius,atm_cnt) - call convertAtomicToGridUnits(X,Y,Z,mnX,mnY,mnZ,gridX,gridY,gridZ,xVxl,yVxl,zVxl) - - zRange = ((totalRadius)/GSP) - atmMnZVxl = ceiling(zVxl - zRange) - atmMxZVxl = floor(zVxl + zRange) + call convertAtomicToGridUnits(X,Y,Z, mxX,mxY,mxZ, xVxl,yVxl,zVxl) + + zRange = totalRadius/GSP + atmMnZVxl = int(zVxl - zRange) + 1 + atmMxZVxl = int(zVxl + zRange) - do n = max(int(atmMnZVxl), int(mnZVxl)), min(int(atmMxZVxl), int(mxZVxl)) + do n = max(atmMnZVxl,mnZVxl,1), min(atmMxZVxl, mxZVxl) yRange = sqrt(zRange*zRange - (zVxl - n)**2) - - atmMnYVxl = ceiling(yVxl - yRange) - atmMxYVxl = floor(yVxl + yRange) - do m = max(int(atmMnYVxl), int(mnYVxl)), min(int(atmMxYVxl), int(mxYVxl)) - xRange = sqrt(yRange*yRange - (yVxl - m)**2) + atmMnYVxl = int(yVxl - yRange) +1 + atmMxYVxl = int(yVxl + yRange) - atmMnXVxl = ceiling(xVxl - xRange) - atmMxXVxl = floor(xVxl + xRange) + do m = max(atmMnYVxl,mnYVxl,1), min(atmMxYVxl, mxYVxl) - do l = max(int(atmMnXVxl), int(mnXVxl)), min(int(atmMxXVxl), int(mxXVxl)) - if (g(l,m,n) .eq. 0 .and. delta .eq. -1) then - g(l,m,n) = 0 - else - g(l,m,n) = g(l,m,n) + delta - end if + xRange = sqrt(yRange*yRange - (yVxl - m)**2) - end do - end do - end do + if (xRange .gt. 0.0d0) then + + atmMnXVxl = int(xVxl - xRange) + 1 + atmMxXVxl = int(xVxl + xRange) -! print *,'END placeInGrid' + do l = max(atmMnXVxl,mnXVxl), min(atmMxXVxl, mxXVxl) -end subroutine placeInGrid + if (g(l,m,n) .eq. 0 .and. delta .eq. -1) then + g(l,m,n) = 0 + else + g(l,m,n) = g(l,m,n) + delta + end if -!******************************************************************************* -! -! Subroutine: makeProteinPdbFormat -! -! Description: This is a visualization function for debugging to generate -! the protein in a pdb format. In this case the number of atoms is hard coded. -! -!******************************************************************************* + end do -subroutine makeProteinPdbFormat(crd,mnX,mnY,mnZ,atm_cnt) + end if - implicit none + end do - integer, intent(in) :: atm_cnt - double precision, intent(in) :: crd(3, atm_cnt) - double precision, intent(in) :: mnX, mnY, mnZ + end do -! local variables - integer :: l, counter - double precision :: X, Y, Z +! print *,'END placeInGrid' - open(unit=7,FILE="protein.pdb", action="write") - counter = 0 - do l = 1, 1836 - - X = crd(1, l) - Y = crd(2, l) - Z = crd(3, l) - counter = counter + 1 - write (7,20) counter,X,Y,Z - 20 format ('ATOM ',i5,1x, ' RES',14x, 3f8.3) - end do - close(unit=7, status = 'keep') +end subroutine placeInGrid -end subroutine makeProteinPdbFormat !******************************************************************************* ! -! Subroutine: makeGridPdbFormat +! Subroutine: determineNumberOfEmptyVoxels ! -! Description: Similar to the previous function, this is a visualization function for debugging. -! This functions generates the grid in a pdb format. +! Description: This function is actually not necessary. it it just if you want +! to know the number of empty voxels ! !******************************************************************************* -subroutine makeGridPdbFormat(g,mxXVxl,mxYVxl,mxZVxl,mnX,mnY,mnZ) +subroutine determineNumberOfEmptyVoxels(g,mxXVxl,mxYVxl,mxZVxl,nmEmpVxls) implicit none - integer, intent(in) :: g(:,:,:) integer, intent(in) :: mxXVxl, mxYVxl, mxZVxl - double precision, intent(in) :: mnX, mnY, mnZ ! , xVxl, yVxl, zVxl + integer, intent(out) :: nmEmpVxls -! local variables - integer :: l,m,n, counter, mnXVxl,mnYVxl,mnZVxl - double precision :: X,Y,Z, xVxl,yVxl,zVxl + ! local variables + integer :: l,m,n,mnXVxl, mnYVxl, mnZVxl,counter - mnXVxl = (SHIFT / GSP) - mnYVxl = (SHIFT / GSP) - mnZVxl = (SHIFT / GSP) + mnXVxl = ceiling(SHIFT / GSP) + mnYVxl = ceiling(SHIFT / GSP) + mnZVxl = ceiling(SHIFT / GSP) - open(unit=7,FILE="grid.pdb", action="write") counter = 0 - do n = mnZVxl, mxZVxl - do m = mnYVxl, mxYVxl - do l = mnXVxl, mxXVxl - - if (g(l,m,n) .eq. 0) then - - xVxl = l - yVxl = m - zVxl = n - - call convertGridToAtomicUnits(mnX,mnY,mnZ,xVxl,yVxl,zVxl,X,Y,Z) - counter = counter + 1 - write (7,20) counter,X, Y, Z - 20 format ('ATOM ',i5,1x, ' WAT',14x, 3f8.3) - end if + do n = mnZVxl, mxZVxl-1 + do m = mnYVxl, mxYVxl-1 + do l = mnXVxl, mxXVxl-1 + + if (g(l,m,n) .eq. 0) then + counter = counter + 1 + + end if end do end do end do - close(unit=7, status = 'keep') - -!print *,'makeGridPdbFormat ' + nmEmpVxls = counter -end subroutine makeGridPdbFormat +end subroutine determineNumberOfEmptyVoxels -!******************************************************************************* -! -! Subroutine: determineNumberOfemptyVoxels -! -! Description: This function is actually not necessary. it it just if you want -! to know the number of empty voxels -! -!******************************************************************************* -subroutine determineNumberOfEmptyVoxels(g,mxXVxl,mxYVxl,mxZVxl,nmEmpVxls) +subroutine updateEmptyVoxelsList(g,mxXVxl,mxYVxl,mxZVxl,nmEmpVxls,listEmptyVoxels) implicit none - integer, intent(in) :: g(:,:,:) - integer, intent(in) :: mxXVxl, mxYVxl, mxZVxl + integer, intent(in) :: mxXVxl, mxYVxl, mxZVxl integer, intent(out) :: nmEmpVxls + integer, intent(inout) :: listEmptyVoxels(:,:) ! local variables - integer :: l,m,n, counter,mnXVxl, mnYVxl, mnZVxl + integer :: l,m,n, counter,mnXVxl, mnYVxl, mnZVxl - mnXVxl = (SHIFT / GSP) - mnYVxl = (SHIFT / GSP) - mnZVxl = (SHIFT / GSP) + mnXVxl = ceiling(SHIFT / GSP) + mnYVxl = ceiling(SHIFT / GSP) + mnZVxl = ceiling(SHIFT / GSP) counter = 0 - do n = mnZVxl, mxZVxl - do m = mnYVxl, mxYVxl - do l = mnXVxl, mxXVxl - if (g(l,m,n) .eq. 0) then + do n = mnZVxl, mxZVxl-1 + do m = mnYVxl, mxYVxl-1 + do l = mnXVxl, mxXVxl-1 + + if (g(l,m,n) .eq. 0) then counter = counter + 1 - end if + + listEmptyVoxels(counter,1) = l + listEmptyVoxels(counter,2) = m + listEmptyVoxels(counter,3) = n + + end if end do end do end do nmEmpVxls = counter - print *,'nmEmpVxls = ', nmEmpVxls -! print *,'END determineNumberOfEmptyVoxels' +end subroutine updateEmptyVoxelsList -end subroutine determineNumberOfEmptyVoxels !******************************************************************************* ! @@ -643,7 +626,7 @@ end subroutine determineNumberOfEmptyVoxels ! !******************************************************************************* -subroutine returnRandomVoxel(mxXVxl,mxYVxl,mxZVxl,tempX,tempY,tempZ) +subroutine returnRandomVoxel(nmEmpVxls,listEmptyVoxels,rndmVxl) #ifdef MPI use parallel_dat_mod @@ -653,54 +636,22 @@ subroutine returnRandomVoxel(mxXVxl,mxYVxl,mxZVxl,tempX,tempY,tempZ) implicit none - integer, intent(in) :: mxXVxl,mxYVxl,mxZVxl - double precision, intent(out) :: tempX,tempY,tempZ + integer, intent(in) :: nmEmpVxls + integer, intent(in) :: listEmptyVoxels(:,:) + integer, intent(out) :: rndmVxl(:) !local variable - double precision :: temp - integer :: mnXVxl, mnYVxl, mnZVxl - - mnXVxl = (SHIFT / GSP) - mnYVxl = (SHIFT / GSP) - mnZVxl = (SHIFT / GSP) + double precision :: temp + integer :: tempVxl call amrand_gen(mcres_randgen, temp) - tempX = (ceiling(temp * (mxXVxl - mnXVxl) ) + mnXVxl) - call amrand_gen(mcres_randgen, temp) - tempY = (ceiling(temp * (mxYVxl - mnYVxl) ) + mnYVxl) - call amrand_gen(mcres_randgen, temp) - tempZ = (ceiling(temp * (mxZVxl - mnZVxl) ) + mnZVxl) - -end subroutine returnRandomVoxel - -!******************************************************************************* -! -! Subroutine: yesNoEmptymVoxel -! -! Description: The function checks if the voxel is empty or occupied. -! -!******************************************************************************* - -subroutine yesNoEmptymVoxel(g,tempX,tempY,tempZ,empty,rndmVxl) - - implicit none - - integer, intent(in) :: g(:,:,:) - double precision, intent(in) :: tempX, tempY, tempZ - integer, intent(out) :: rndmVxl(:) - logical, intent(out) :: empty - - empty = .false. + tempVxl = ceiling(temp * nmEmpVxls) - if ( g(int(tempX),int(tempY),int(tempZ)) .eq. 0) then - empty = .true. + rndmVxl(1) = listEmptyVoxels(tempVxl,1) + rndmVxl(2) = listEmptyVoxels(tempVxl,2) + rndmVxl(3) = listEmptyVoxels(tempVxl,3) - rndmVxl(1) = tempX - rndmVxl(2) = tempY - rndmVxl(3) = tempZ - end if - -end subroutine yesNoEmptymVoxel +end subroutine returnRandomVoxel !******************************************************************************* ! @@ -710,55 +661,81 @@ end subroutine yesNoEmptymVoxel ! !******************************************************************************* -subroutine returnRandomWaterMolecule(crd,rndmWtr,mol_str,atm_cnt,mxXVxl,mxYVxl,mxZVxl) +subroutine returnRandomWaterMolecule(crd,rndmWtr,mol_str,atm_cnt,mxXVxl,mxYVxl,mxZVxl, mxX,mxY,mxZ) #ifdef MPI use parallel_dat_mod #endif use pmemd_lib_mod + use mdin_ctrl_dat_mod + use file_io_dat_mod + use prmtop_dat_mod + use mol_list_mod + implicit none integer, intent(in) :: atm_cnt, mxXVxl,mxYVxl,mxZVxl double precision, intent(in) :: crd(3, atm_cnt) integer, intent(out) :: rndmWtr character(4), intent(in) :: mol_str - + double precision, intent(in) :: mxX,mxY,mxZ ! local variables logical :: yesWat double precision :: temp - integer :: mnXVxl, mnYVxl, mnZVxl + integer :: mnXVxl, mnYVxl, mnZVxl + integer :: atomNumber + + double precision :: watX,watY,watZ + integer :: xVxl,yVxl,zVxl - mnXVxl = (SHIFT / GSP) - mnYVxl = (SHIFT / GSP) - mnZVxl = (SHIFT / GSP) + mnXVxl = ceiling(SHIFT / GSP) + mnYVxl = ceiling(SHIFT / GSP) + mnZVxl = ceiling(SHIFT / GSP) yesWat = .false. + do while (yesWat .eqv. .false.) + call amrand_gen(mcres_randgen, temp) - rndmWtr = ceiling(temp*nres) + rndmWtr = ceiling(temp*nres) + + atomNumber = gbl_res_atms(rndmWtr) + if(gbl_labres(rndmWtr) .eq. mol_str) then - if((crd(1, rndmWtr)/GSP) .gt. mxXVxl) then + + watX = crd(1, atomNumber) + watY = crd(2, atomNumber) + watZ = crd(3, atomNumber) + + call convertAtomicToGridUnits(watX,watY,watZ, mxX,mxY,mxZ, xVxl,yVxl,zVxl) + + if(xVxl .gt. (mxXVxl + 0.5) ) then yesWat = .false. - else if ((crd(1, rndmWtr)/GSP) .lt. mnXVxl) then + else if (xVxl .lt. (mnXVxl - 0.5) ) then yesWat = .false. - else if ((crd(2, rndmWtr)/GSP) .gt. mxYVxl) then + else if (yVxl .gt. (mxYVxl + 0.5) ) then yesWat = .false. - else if ((crd(2, rndmWtr)/GSP) .lt. mnYVxl) then + else if (yVxl .lt. (mnYVxl - 0.5) ) then yesWat = .false. - else if ((crd(3, rndmWtr)/GSP) .gt. mxZVxl) then + else if (zVxl .gt. (mxZVxl + 0.5) ) then yesWat = .false. - else if ((crd(3, rndmWtr)/GSP) .lt. mnZVxl) then + else if (zVxl .lt. (mnZVxl - 0.5) ) then yesWat = .false. else yesWat = .true. + end if - end if + + + + end if + end do -!! print *,'END returnRandomWaterMolecule' end subroutine returnRandomWaterMolecule + !******************************************************************************* ! ! Subroutine: plcWtrInVxl @@ -772,8 +749,8 @@ end subroutine returnRandomWaterMolecule ! !******************************************************************************* -subroutine plcWtrInVxl(frc,crd,rndmVxl,rndmWtr,X,Y,Z,mnX,mnY,mnZ,atm_cnt,atm_img_map,img_atm_map,& - mol_str,old_ene,nstep,mxXVxl,mxYVxl,mxZVxl,g,cnt,virial,ekcmt,reject,pme_err_est) +subroutine plcWtrInVxl(atm_cnt,crd,rndmVxl,rndmWtr,X,Y,Z,& + mol_str,mxXVxl,mxYVxl,mxZVxl,g,cnt, update_crd) use mol_list_mod use mdin_ctrl_dat_mod @@ -792,32 +769,24 @@ subroutine plcWtrInVxl(frc,crd,rndmVxl,rndmWtr,X,Y,Z,mnX,mnY,mnZ,atm_cnt,atm_img use cit_mod implicit none - - integer, intent(in) :: atm_cnt, mxXVxl,mxYVxl,mxZVxl, nstep - double precision, intent(in) :: frc(3, atm_cnt) + + integer :: atm_cnt double precision, intent(inout) :: crd(3, atm_cnt) + integer, intent(in) :: mxXVxl,mxYVxl,mxZVxl integer, intent(inout) :: rndmVxl(:) double precision, intent(out) :: X,Y,Z - double precision, intent(in) :: mnX, mnY, mnZ - integer,intent(inout) :: atm_img_map(atm_cnt) - integer,intent(inout) :: img_atm_map(atm_cnt) character(4), intent(in) :: mol_str - double precision, intent(inout) :: old_ene integer, intent(inout) :: g(:,:,:) - integer, intent(inout) :: reject,rndmWtr - double precision, intent(inout) :: virial(3) ! Only used for MD - double precision, intent(inout) :: ekcmt(3) ! Only used for MD - double precision :: pme_err_est ! Only used for MD + integer, intent(inout) :: rndmWtr + logical, intent(in) :: update_crd - ! LAPACK disnan - logical :: disnan !local variables - integer :: atom,l,m,n,k,i,j,f,b, numberOfAtoms + integer :: atom,l double precision :: temp, tempX, tempY, tempZ, oldX, oldY, oldZ - double precision :: deltaX,deltaY,deltaZ, watXVxl, watYVxl, watZVxl - double precision :: watX,watY,watZ, xVxl,yVxl,zVxl, new_ene - + double precision :: watXVxl, watYVxl, watZVxl !!! TODO - why double precision + double precision :: watX,watY,watZ + integer :: xVxl,yVxl,zVxl logical accept type(pme_pot_ene_rec) :: pot_ene @@ -826,27 +795,13 @@ subroutine plcWtrInVxl(frc,crd,rndmVxl,rndmWtr,X,Y,Z,mnX,mnY,mnZ,atm_cnt,atm_img logical :: need_virials logical :: new_list double precision :: random - double precision :: recip_11, recip_22, recip_33, scale_fac_x, scale_fac_y, scale_fac_z - - double precision :: f1, f2, f3, fractx, fracty, fractz - integer :: bktx, bkty, bktz, cur_bkt - double precision :: gridX,gridY,gridZ - integer :: numberofmoleculesinatom - #ifdef MPI integer :: my_atm_lst(my_atm_cnt) #endif integer, intent(inout) :: cnt ! inout or only out - -! this might be used for the coarse grid - - type(cit_tbl_rec) :: flat_cit(0 : cit_tbl_x_dim * & - cit_tbl_y_dim * & - cit_tbl_z_dim - 1) - new_list = .true. #ifdef CUDA @@ -871,149 +826,198 @@ subroutine plcWtrInVxl(frc,crd,rndmVxl,rndmWtr,X,Y,Z,mnX,mnY,mnZ,atm_cnt,atm_img call amrand_gen(mcres_randgen, temp) tempZ = temp - oldX = X !saving the old coordinates in case the move is rejected - oldY = Y - oldZ = Z + call convertAtomicToGridUnits(X,Y,Z, mxX,mxY,mxZ, xVxl,yVxl,zVxl) - call convertAtomicToGridUnits(X,Y,Z,mnX,mnY,mnZ,gridX,gridY,gridZ,xVxl,yVxl,zVxl) + watXVxl = rndmVxl(1) - 0.5 + watXVxl = watXVxl + tempX -!!! 3 scenarios -! X coordinates - if (rndmVxl(1) .eq. 1) then - watXVxl = rndmVxl(1) ! the coordinates of the grid map - watXVxl = watXVxl + (tempX*0.5) ! placing the atom within 0.5 grid spacing in all directions - else if (rndmVxl(1) .eq. mxXVxl) then - watXVxl = (rndmVxl(1) - 0.5) ! the coordinates of the grid map - watXVxl = watXVxl + (tempX*0.5) ! placing the atom within 0.5 grid spacing in all directions - else - watXVxl = (rndmVxl(1) - 0.5) ! the coordinates of the grid map - watXVxl = watXVxl + tempX ! placing the atom within 0.5 grid spacing in all directions - end if - -! Y coordinates - if (rndmVxl(2) .eq. 1) then - watYVxl = rndmVxl(2) ! the coordinates of the grid map - watYVxl = watYVxl + (tempY*0.5) ! placing the atom within 0.5 grid spacing in all directions - else if (rndmVxl(2) .eq. mxYVxl) then - watYVxl = (rndmVxl(2) - 0.5) ! the coordinates of the grid map - watYVxl = watYVxl + (tempY*0.5) ! placing the atom within 0.5 grid spacing in all directions - else - watYVxl = (rndmVxl(2) - 0.5) ! the coordinates of the grid map - watYVxl = watYVxl + tempY ! placing the atom within 0.5 grid spacing in all directions - end if - -! Z coordinates - if (rndmVxl(3) .eq. 1) then - watZVxl = rndmVxl(3) ! the coordinates of the grid map - watZVxl = watZVxl + (tempZ*0.5) ! placing the atom within 0.5 grid spacing in all directions - else if (rndmVxl(3) .eq. mxZVxl) then - watZVxl = (rndmVxl(3) - 0.5) ! the coordinates of the grid map - watZVxl = watZVxl + (tempZ*0.5) ! placing the atom within 0.5 grid spacing in all directions - else - watZVxl = (rndmVxl(3) - 0.5) ! the coordinates of the grid map - watZVxl = watZVxl + tempZ ! placing the atom within 0.5 grid spacing in all directions - end if + watYVxl = rndmVxl(2) - 0.5 + watYVxl = watYVxl + tempY - deltaX = (xVxl - watXVxl) ! calculating the shift between the old and new location in GRID units - deltaY = (yVxl - watYVxl) - deltaZ = (zVxl - watZVxl) + watZVxl = rndmVxl(3) - 0.5 + watZVxl = watZVxl + tempZ - numberOfAtoms = 3 -! numberOfAtoms = gbl_mol_atms_listdata(atom)%cnt !!!!! TODO correct this for the future + deltacrd(1) = X - (watXVxl*GSP) ! calculating the shift between the old and new location in ATOMIC units + deltacrd(2) = Y - (watYVxl*GSP) + deltacrd(3) = Z - (watZVxl*GSP) + + numberOfAtoms = mcres_numatms do l = 0, numberOfAtoms-1 - crd(1,cnt+l) = (crd(1,cnt+l) - (deltaX*GSP)) - crd(2,cnt+l) = (crd(2,cnt+l) - (deltaY*GSP)) - crd(3,cnt+l) = (crd(3,cnt+l) - (deltaZ*GSP)) + crd(1,cnt+l) = (crd(1,cnt+l) - deltacrd(1)) + crd(2,cnt+l) = (crd(2,cnt+l) - deltacrd(2)) + crd(3,cnt+l) = (crd(3,cnt+l) - deltacrd(3)) end do accept = .false. new_list = .true. -#ifdef CUDA - ! Upload new coordinates onto GPU for pme force calculation. CUDA - ! implementation could be faster by putting the mcres code on a cuda kernel - ! to avoid the need for the memory transfer. We can do this if testing shows - ! this approach to be useful - and we converge on a moveset. +end subroutine plcWtrInVxl - ! Note this is very inefficient right now - we need a gpu function that just uploads - ! the coordinates of the molecule that was moved. - call gpu_force_new_neighborlist() - call gpu_upload_crd(crd) + +!******************************************************************************* +! +! Subroutine: mcresenergy +! +! Description: +! computes the energy +! difference for the trial move, and uses the Metropolis criterion to accept or +! reject the move. If the move is accepted, then the grid counts at the new +! location are incremented by 1 and grid counts at the old location +! are decremented by 1 (add -1). +! +!******************************************************************************* + +subroutine mcresenergy(frc,crd,rndmVxl,rndmWtr,X,Y,Z,atm_cnt,atm_img_map,img_atm_map,& + mol_str,old_ene,nstep,mxXVxl,mxYVxl,mxZVxl,g,cnt,virial,ekcmt,pme_err_est,& + revert,nmEmpVxls,listEmptyVoxels) + + + use mol_list_mod + use mdin_ctrl_dat_mod + use pbc_mod + use gbl_constants_mod, only : KB + use energy_records_mod, only : pme_pot_ene_rec + use file_io_dat_mod + use pme_force_mod + use random_mod, only : random_state, amrand_gen, amrset_gen + use prmtop_dat_mod +#ifdef MPI + use parallel_dat_mod #endif + use pmemd_lib_mod + + use cit_mod + + implicit none + + integer, intent(in) :: atm_cnt, mxXVxl,mxYVxl,mxZVxl, nstep + double precision, intent(in) :: frc(3, atm_cnt) + double precision, intent(inout) :: crd(3, atm_cnt) + integer, intent(inout) :: rndmVxl(:) + double precision, intent(out) :: X,Y,Z + integer,intent(inout) :: atm_img_map(atm_cnt) + integer,intent(inout) :: img_atm_map(atm_cnt) + character(4), intent(in) :: mol_str + double precision, intent(inout) :: old_ene + integer, intent(inout) :: g(:,:,:) + double precision, intent(inout) :: virial(3) ! Only used for MD + double precision, intent(inout) :: ekcmt(3) ! Only used for MD + double precision :: pme_err_est ! Only used for MD + logical, intent (inout) :: revert + integer, intent(inout) :: rndmWtr,nmEmpVxls + integer, intent(inout) :: listEmptyVoxels(:,:) + + logical :: disnan + +!local variables + integer :: l, numberOfAtoms + double precision :: new_ene + logical :: accept, full_eval + type(pme_pot_ene_rec) :: pot_ene + double precision :: acceptance_criteria + logical :: need_virials + logical :: new_list + double precision :: random - ! Get new energy #ifdef MPI - call pme_force(atm_cnt, crd, frc, img_atm_map, atm_img_map, & - my_atm_lst, new_list, .true., need_virials, & - pot_ene, nstep, virial, ekcmt, pme_err_est) + integer :: my_atm_lst(my_atm_cnt) +#endif + + integer, intent(inout) :: cnt ! inout or only out + double precision :: delta_full + + + new_list = .true. + +#ifdef CUDA + ! Upload new coordinates onto GPU for pme force calculation. CUDA + ! implementation could be faster by putting the mcres code on a cuda kernel + ! to avoid the need for the memory transfer. We can do this if testing shows + ! this approach to be useful - and we converge on a moveset. + + ! Note this is very inefficient right now - we need a gpu function that just uploads + ! the coordinates of the molecule that was moved. + call gpu_force_new_neighborlist() + call gpu_upload_crd(crd) +#endif + + ! Get new energy +#ifdef MPI + call pme_force(atm_cnt, crd, frc, img_atm_map, atm_img_map, & + my_atm_lst, new_list, .true., need_virials, & + pot_ene, nstep, virial, ekcmt, pme_err_est) #else - call pme_force(atm_cnt, crd, frc, img_atm_map, atm_img_map, & - .true., .true., .false., & - pot_ene, nstep, virial, ekcmt, pme_err_est) + call pme_force(atm_cnt, crd, frc, img_atm_map, atm_img_map, & + .true., .true., .false., & + pot_ene, nstep, virial, ekcmt, pme_err_est) #endif - new_ene = pot_ene%total + new_ene = pot_ene%total - if(disnan(old_ene) .or. disnan(new_ene)) then - accept = .false. - else - ! Check acceptance criteria + numberOfAtoms = mcres_numatms ! mcres_numatms looks for the wat residue to get the numberof atoms in water - if(old_ene .gt. new_ene) then - ! Temp hack to prevent overflow. - if(abs(new_ene - old_ene) .lt. mcwatmaxdif) then - accept = .true. - end if - else + if(disnan(old_ene) .or. disnan(new_ene)) then + accept = .false. + else + ! Check acceptance criteria - call amrand_gen(mcres_randgen, random) + delta_full = (new_ene - old_ene) - if (mcverbosity .ge. 5) then - write(mdout, *) "call amrand_gen got: ", random - endif +! print *, 'delta_full, new_ene, old_ene = ', delta_full, new_ene, old_ene - acceptance_criteria = exp(-(new_ene-old_ene)/(temp0*KB)) + if(old_ene .gt. new_ene) then + ! Temp hack to prevent overflow. + if(abs(new_ene - old_ene) .lt. mcwatmaxdif) then + accept = .true. + end if + else -!print *, 'new_ene, old_ene = ', new_ene, old_ene + call amrand_gen(mcres_randgen, random) - if (mcverbosity .ge. 5) write(mdout, *) "Acceptance: ", acceptance_criteria + if (mcverbosity .ge. 5) then + write(mdout, *) "call amrand_gen got: ", random + endif - if(random .lt. acceptance_criteria) then - accept = .true. - end if - end if - end if + acceptance_criteria = exp(-(new_ene-old_ene)/(temp0*KB)) + + + if (mcverbosity .ge. 5) write(mdout, *) "Acceptance: ", acceptance_criteria + + if(random .lt. acceptance_criteria) then + accept = .true. + end if + end if + end if - ! We do a force check after we have accepted because force download from - ! gpu is fairly expensive we only want to do the download after we - ! accept. This check exists in case we mistakenly accept due to an - ! energy overflow. + ! We do a force check after we have accepted because force download from + ! gpu is fairly expensive we only want to do the download after we + ! accept. This check exists in case we mistakenly accept due to an + ! energy overflow. - if(accept .eqv. .true.) then + if(accept) then #ifdef CUDA call gpu_download_frc(frc) #endif - do l = 0, numberOfAtoms-1 + do l = 0, numberOfAtoms-1 - if(abs(frc(1,cnt+l)) .gt. mcwatmaxdif) then + if(abs(frc(1,cnt+l)) .gt. mcwatmaxdif) then accept = .false. - else if(abs(frc(2,cnt+l)) .gt. mcwatmaxdif) then + else if(abs(frc(2,cnt+l)) .gt. mcwatmaxdif) then accept = .false. - else if(abs(frc(3,cnt+l)) .gt. mcwatmaxdif) then + else if(abs(frc(3,cnt+l)) .gt. mcwatmaxdif) then accept = .false. - end if - end do - -! if(abs(new_ene - old_ene) .lt. mcwatmaxdif) then + end if + end do + end if -! an internal if condition, if it is true - if (accept .eqv. .true.) then ! needs to update grid + if (accept) then ! needs to update grid delta = -1 ! -1 meaning deleting the water - call placeInGrid(crd,g,atm_cnt,mnX,mnY,mnZ,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,X,Y,Z,nmEmpVxls,cnt) + + call placeInGrid(crd,g,atm_cnt,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,X,Y,Z,nmEmpVxls,cnt) X = crd(1,cnt) Y = crd(2,cnt) @@ -1021,53 +1025,361 @@ subroutine plcWtrInVxl(frc,crd,rndmVxl,rndmWtr,X,Y,Z,mnX,mnY,mnZ,atm_cnt,atm_img delta = 1 ! +1 meaning inserting the water mol - call placeInGrid(crd,g,atm_cnt,mnX,mnY,mnZ,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,X,Y,Z,nmEmpVxls,cnt) + call placeInGrid(crd,g,atm_cnt,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,X,Y,Z,nmEmpVxls,cnt) write(mdout, '(a,F12.4,a,F12.4)')"Old Energy: ", old_ene, " accepted New Energy: ", new_ene write(mdout, '(a,i9,a,a,a,i6,a)')"MCRES ACCEPT. Step: ", nstep,& " Residue: ",mol_str, ", Residue #: ", rndmWtr,"." -! write(mdout, '(a)')"" -!!! for testing purposes can also print the forces -! do l = 0, numberOfAtoms-1 ! return the molecule , the calculation is correct -! write(mdout,*)"Force: ",frc(1,atom+l),frc(2,atom+l),frc(3,atom+l) -!! write(mdout, '(a)')"" -! end do + revert = .false. - old_ene = new_ene ! only if it is true the energy is updated + call determineNumberOfEmptyVoxels(g,mxXVxl,mxYVxl,mxZVxl,nmEmpVxls) + + if (stat.ne.0) then + print *,'error: couldnt allocate memory for array, listEmptyVoxels=',stat + stop 1 + endif - end if + call updateEmptyVoxelsList(g,mxXVxl,mxYVxl,mxZVxl,nmEmpVxls,listEmptyVoxels) + + old_ene = new_ene ! only if it is true the energy is updated - else if (accept .eqv. .false.) then + else ! if (accept .eqv. .false.) then do l = 0, numberOfAtoms-1 ! return the molecule , the calculation is correct - crd(1,cnt+l) = crd(1,cnt+l) + (deltaX*GSP) - crd(2,cnt+l) = crd(2,cnt+l) + (deltaY*GSP) - crd(3,cnt+l) = crd(3,cnt+l) + (deltaZ*GSP) + + crd(1,cnt+l) = crd(1,cnt+l) + deltacrd(1) + crd(2,cnt+l) = crd(2,cnt+l) + deltacrd(2) + crd(3,cnt+l) = crd(3,cnt+l) + deltacrd(3) + end do - new_list = .true. + new_list = .true. + revert = .true. #ifdef CUDA - ! Uploads old coordinates back on to GPU if rejected. - ! Optimization needed here for only uploading the coordinates of the - ! atom that was moved. - call gpu_force_new_neighborlist() - call gpu_upload_crd(crd) + ! Uploads old coordinates back on to GPU if rejected. + ! Optimization needed here for only uploading the coordinates of the + ! atom that was moved. + call gpu_force_new_neighborlist() + call gpu_upload_crd(crd) #endif -!!! for testing purposes can also print the rejected steps -! write(mdout, '(a,E12.4,a,E12.4)')"Old Energy: ", old_ene, " rejected New Energy: ", new_ene -! write(mdout, '(a,i9,a,a,a,i6,a)')"MCRES REJECT. Step: ", nstep,& -! " Residue: ",mol_str, ", Residue #: ", rndmWtr,"." -! write(mdout, '(a)')"" + end if + +end subroutine mcresenergy + +subroutine findMxCoarseGrid(mxX, mxY, mxZ,coarseGSP, maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl) - reject = reject + 1 + implicit none -end if + double precision, intent(in) :: mxX, mxY, mxZ + integer, intent(out) :: maxCoarseXVxl, maxCoarseYVxl, maxCoarseZVxl + double precision, intent(out) :: coarseGSP -! print *,'END plcWtrInVxl' + coarseGSP = (es_cutoff/2) -end subroutine plcWtrInVxl + maxCoarseXVxl = ceiling(mxX /coarseGSP) + maxCoarseYVxl = ceiling(mxY /coarseGSP) + maxCoarseZVxl = ceiling(mxZ /coarseGSP) + +end subroutine findMxCoarseGrid + + +subroutine initializeCoarseGrid(coarseGrid,maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl) + + integer, intent(in) :: maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl + integer, intent(inout) :: coarseGrid(:,:,:,:) + +!local variables + integer :: l, m, n, o, numAtoms + + numAtoms = coarseGSP*5 + + do o = 1, numAtoms + do n = 1, maxCoarseZVxl + do m = 1, maxCoarseYVxl + do l = 1, maxCoarseXVxl + coarseGrid(l, m, n, o) = 0 + end do + end do + end do + end do + +end subroutine initializeCoarseGrid + + +subroutine buildCoarseGrid(crd,coarseGrid,atm_cnt,coarseGSP, maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl,mxX,mxY,mxZ) + + use mdin_ctrl_dat_mod + use pbc_mod, only : pbc_box + + implicit none + + double precision, intent(in) :: crd(3, atm_cnt) + double precision, intent(in) :: mxX,mxY,mxZ, coarseGSP + integer, intent(in) :: atm_cnt + integer, intent(inout) :: coarseGrid(:,:,:,:) + integer, intent(in) :: maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl + +!local variables + double precision :: X,Y,Z + integer :: tempX,tempY,tempZ, counter + + do cnt = 1, atm_cnt + X = crd(1, cnt) ! atomic units + Y = crd(2, cnt) + Z = crd(3, cnt) + + call wrapCoordinates(X,Y,Z,mxX,mxY,mxZ) + + tempX = ceiling(X / coarseGSP) + tempY = ceiling(Y / coarseGSP) + tempZ = ceiling(Z / coarseGSP) + + coarseGrid(tempX,tempY,tempZ,1) = (coarseGrid(tempX,tempY,tempZ,1) + 1) + coarseGrid(tempX,tempY,tempZ,(coarseGrid(tempX,tempY,tempZ,1)+1)) = cnt + end do + +end subroutine buildCoarseGrid + + +subroutine updateCoarseGrid(atm_cnt,mol_id,crd,coarseGrid,maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl,coarseGSP,revert, & +& mxX,mxY,mxZ,rndmWtr) + + use mol_list_mod + use pbc_mod + + implicit none + + integer, intent(in) :: atm_cnt + integer, intent(in) :: mol_id + double precision, intent(in) :: crd(3,atm_cnt) + integer, intent(inout) :: coarseGrid(:,:,:,:) + double precision, intent(in) :: coarseGSP + integer, intent(in) :: maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl + logical, intent(in) :: revert ! Reverts the grid backwards if true + double precision, intent(in) :: mxX,mxY,mxZ + integer, intent(in) :: rndmWtr + + + double precision :: new_crd(3) + integer :: gridold(3) + integer :: gridnew(3) + integer :: tempgrid(3) + double precision :: X,Y,Z + logical :: found + integer :: i, j, mol_atm_cnt, mol_offset, gridAtms,numberOfAtoms + + mol_offset = gbl_res_atms(rndmWtr) + + numberOfAtoms = mcres_numatms + + do j = mol_offset, mol_offset + numberOfAtoms - 1 + if(.not. revert) then + + !Update old + + X = modulo((crd(1,j)+(deltacrd(1))),mxX) + Y = modulo((crd(2,j)+(deltacrd(2))),mxY) + Z = modulo((crd(3,j)+(deltacrd(3))),mxZ) + + gridold(1) = ceiling(X / coarseGSP) + gridold(2) = ceiling(Y / coarseGSP) + gridold(3) = ceiling(Z / coarseGSP) + + !Update new + X = crd(1,j) + Y = crd(2,j) + Z = crd(3,j) + + call wrapCoordinates(X,Y,Z,mxX,mxY,mxZ) + + gridnew(1) = ceiling(X / coarseGSP) + gridnew(2) = ceiling(Y / coarseGSP) + gridnew(3) = ceiling(Z / coarseGSP) + + else + + X = modulo((crd(1,j)-(deltacrd(1))),mxX) + Y = modulo((crd(2,j)-(deltacrd(2))),mxY) + Z = modulo((crd(3,j)-(deltacrd(3))),mxZ) + + gridold(1) = ceiling(X / coarseGSP) + gridold(2) = ceiling(Y / coarseGSP) + gridold(3) = ceiling(Z / coarseGSP) + + X = crd(1,j) + Y = crd(2,j) + Z = crd(3,j) + + call wrapCoordinates(X,Y,Z,mxX,mxY,mxZ) + + gridnew(1) = ceiling(X / coarseGSP) + gridnew(2) = ceiling(Y / coarseGSP) + gridnew(3) = ceiling(Z / coarseGSP) + + end if + + found = .false. + + ! Remove ID from grid + gridAtms = coarseGrid(gridold(1),gridold(2),gridold(3),1)+1 + do i = 2, gridAtms + if(found) then + coarseGrid(gridold(1),gridold(2),gridold(3),i-1) = coarseGrid(gridold(1),gridold(2),gridold(3),i) + if(i .eq. coarseGrid(gridold(1),gridold(2),gridold(3),1)+1) then + coarseGrid(gridold(1),gridold(2),gridold(3),i) = 0 + end if + else + if(coarseGrid(gridold(1),gridold(2),gridold(3),i) .eq. j) then + found = .true. + coarseGrid(gridold(1),gridold(2),gridold(3),i) = 0 + end if + end if + end do + if(found) then + coarseGrid(gridold(1),gridold(2),gridold(3),1) = (coarseGrid(gridold(1),gridold(2),gridold(3),1) - 1) + + ! Add ID to grid + coarseGrid(gridnew(1),gridnew(2),gridnew(3),1) = (coarseGrid(gridnew(1),gridnew(2),gridnew(3),1) + 1) + i = coarseGrid(gridnew(1),gridnew(2),gridnew(3),1) + coarseGrid(gridnew(1),gridnew(2),gridnew(3),i+1)=j + else + write(0,*)"No match" + end if + end do + +end subroutine updateCoarseGrid + + +subroutine calculateCoarseEnergy(atm_cnt, mol_id, crd, CoarseGrid, coarseGSP, mX, mY, mZ, coarseene,rndmWtr) + + use mol_list_mod + use pbc_mod + use mdin_ctrl_dat_mod + use prmtop_dat_mod, only : ntypes, typ_ico, gbl_cn1, gbl_cn2, atm_qterm + use ene_frc_splines_mod + + implicit none + + integer, intent(in) :: atm_cnt, mol_id, rndmWtr + double precision, intent(in) :: crd(3,atm_cnt) + integer, intent(inout) :: coarseGrid(:,:,:,:) + double precision, intent(in) :: coarseGSP + integer, intent(in) :: mX, mY, mZ ! max of coarse grid + double precision :: coarseene + + double precision :: X,Y,Z + integer :: grid(3) + integer :: i,j,k, mol_atm_cnt, res_offset + integer :: bkt_x,bkt_y,bkt_z, bkt_offset,bkt_atms + integer :: wrap_bkt_x, wrap_bkt_y, wrap_bkt_z + integer :: iaci, ic + double precision :: atm_i(3) + double precision :: atm_j(3) + double precision :: cn1, cn2, delx, dely, delz, delr2 + double precision :: charge_i, charge_j + double precision :: delr2inv, r6, f12, f6 + integer :: ecnt, ljcnt + double precision :: elec_ene, lj_ene + double precision :: dens_efs, del_efs, du, du2, du3, b0 + integer :: ind + + ecnt = 0 + ljcnt = 0 + elec_ene = 0.0d0 + lj_ene = 0.0d0 + + bkt_offset = 3 + + res_offset = gbl_res_atms(rndmWtr) + coarseene = 0.0 + do j = res_offset, res_offset + numberOfAtoms-1 + + atm_i(1) = crd(1,j) + atm_i(2) = crd(2,j) + atm_i(3) = crd(3,j) + + grid(1) = ceiling(atm_i(1) / coarseGSP) + grid(2) = ceiling(atm_i(2) / coarseGSP) + grid(3) = ceiling(atm_i(3) / coarseGSP) + + iaci = ntypes * (atm_iac(j) - 1) + + do bkt_x=-bkt_offset,bkt_offset + wrap_bkt_x = modulo(grid(1)+bkt_x-1,mX)+1 + + do bkt_y=-bkt_offset,bkt_offset + wrap_bkt_y = modulo(grid(2)+bkt_y-1,mY)+1 + + do bkt_z=-bkt_offset,bkt_offset + wrap_bkt_z = modulo(grid(3)+bkt_z-1,mZ)+1 + bkt_atms = coarseGrid(wrap_bkt_x,wrap_bkt_y,wrap_bkt_z,1) + + do i=1,bkt_atms ! We ignore a distance check here and do all calcs + + k = coarseGrid(wrap_bkt_x, wrap_bkt_y, wrap_bkt_z, i+1) + + atm_j(1) = crd(1,k) + atm_j(2) = crd(2,k) + atm_j(3) = crd(3,k) + + delx = (atm_i(1) - atm_j(1)) - pbc_box(1)* (nint((atm_i(1)-atm_j(1))/ pbc_box(1))) + dely = (atm_i(2) - atm_j(2)) - pbc_box(2)* (nint((atm_i(2)-atm_j(2))/ pbc_box(2))) + delz = (atm_i(3) - atm_j(3)) - pbc_box(3)* (nint((atm_i(3)-atm_j(3))/ pbc_box(3))) + + delr2 = delx * delx + dely * dely + delz * delz + ic = typ_ico(iaci + atm_iac(k)) + + if (delr2 .le. 0.5001d0 .and. delr2 .gt. 0.0001d0 ) then + + if ( abs(j - k) .ge. mcres_numatms ) then ! this lines checks that the difference between the atom numbering is less than the atoms in that water model. + ! theoretically, it could be between atoms that are not from the same molecule, but we guess that if it happens + ! it is very rare, and it will just result in an additional full energy call. + coarseene = 50000000 + lj_ene = 25000000 + elec_ene = 25000000 + exit + + end if + + else if(delr2 .gt. 0.5001d0 .and. delr2 .le. vdw_cutoff*vdw_cutoff) then + + if(ic .gt. 0) then + cn1 = gbl_cn1(ic) + cn2 = gbl_cn2(ic) + delr2inv = 1.0d0/delr2 + r6 = delr2inv * delr2inv *delr2inv + f6 = cn2 * r6 + f12 = cn1 * r6 * r6 + lj_ene = lj_ene + f12 - f6 + ljcnt = ljcnt + 1 + end if + + charge_i = atm_qterm(j) + charge_j = atm_qterm(k) + dens_efs = efs_tbl_dens + del_efs = 1.0d0/dens_efs + ind = int(dens_efs * delr2) + du = delr2 - dble(ind) * del_efs + du2 = du * du + du3 = du * du2 + ind = ishft(ind, 3) ! 8 * ind + + elec_ene = elec_ene + charge_i*charge_j* (efs_tbl(1 + ind) + du * efs_tbl(2 + ind) + & + du2 * efs_tbl(3 + ind) + du3 * efs_tbl(4 + ind)) + ecnt = ecnt + 1 + end if + + end do + end do + end do + end do + end do + + coarseene = lj_ene + elec_ene + +end subroutine calculateCoarseEnergy !******************************************************************************* @@ -1127,32 +1439,24 @@ subroutine mcres(accept,mol_str, atm_cnt, crd, frc, & use parallel_dat_mod #endif use pmemd_lib_mod - use barostats_mod, only : mcbar_trial, mcbar_summary - use degcnt_mod - use dynamics_mod - use dynamics_dat_mod -#ifdef EMIL - use emil_mod -#endif use pbc_mod use pme_force_mod - use file_io_mod use file_io_dat_mod use img_mod use nb_pairlist_mod - use parallel_dat_mod - use parallel_mod use pmemd_lib_mod #ifdef MPI - use multipmemd_mod, only : free_comms -! use remd_exchg_mod + use parallel_dat_mod + use parallel_mod #endif use runfiles_mod use timers_mod use state_info_mod -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +use bintraj_mod + +implicit none character(4) :: mol_str integer :: atm_cnt @@ -1168,6 +1472,7 @@ subroutine mcres(accept,mol_str, atm_cnt, crd, frc, & double precision :: ekcmt(3) ! Only used for MD double precision :: pme_err_est ! Only used for MD integer, intent(in) :: nstep + logical :: revert #ifdef MPI integer :: my_atm_lst(my_atm_cnt) #endif @@ -1182,17 +1487,13 @@ subroutine mcres(accept,mol_str, atm_cnt, crd, frc, & integer cyc_cnt integer mol_atm_cnt real current_time - double precision :: crd_wrap(3,atm_cnt) - + double precision :: coarseene integer total_nstep + integer l + +!!! for the coarse energy + double precision :: coarse_before_ene, coarse_after_ene, delta_coarse - ! This is hardcoded should be old_crd(3,mol_atm_cnt) but since we don't - ! know what mol we'll be dealing with until mol_str finds the molecule, this - ! is hardcoded to only work with water at the moment. If fortran allowed - ! variable declaration in the middle of the code this would work a lot - ! better. - double precision old_crd(3,mcres_numatms) - double precision new_crd(3,mcres_numatms) double precision :: vel(3, atm_cnt) @@ -1222,52 +1523,106 @@ subroutine mcres(accept,mol_str, atm_cnt, crd, frc, & ! Store old energy old_ene = pot_ene%total - call wrap_for_MC(atm_cnt,crd,crd_wrap) -! call amrset_gen(mcres_randgen,ig) - - call findMnMxnPBox(crd,atm_cnt,mnX,mxX,mnY,mxY,mnZ,mxZ) - call findMxVxls(mnX,mxX,mnY,mxY,mnZ,mxZ,mxXVxl,mxYVxl,mxZVxl) + call findMnMxnPBox(crd,atm_cnt,mxX,mxY,mxZ) + call findMxVxls(mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl) ! Allocating the grid - not in a separate function - allocate(g(0:mxXVxl+50, 0:mxYVxl+50, 0:mxZVxl+50),stat=error) + allocate(g(mxXVxl+5, mxYVxl+5, mxZVxl+5),stat=error) if (stat.ne.0) then - print *,'error: couldnt allocate memory for array, g=',stat - stop + write(0,*)"error: couldnt allocate memory for array, g=",stat + stop 1 endif + coarseGSP = (es_cutoff/2) + call initializeGrid(g,mxXVxl,mxYVxl,mxZVxl) - call buildGrid(crd,g,atm_cnt,mnX,mnY,mnZ,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,cnt) + call buildGrid(crd,g,atm_cnt,mxX,mxY,mxZ,mxXVxl,mxYVxl,mxZVxl,delta,cnt) -! call makeProteinPdbFormat(crd,mnX,mnY,mnZ,atm_cnt) -! call makeGridPdbFormat(g,mxXVxl,mxYVxl,mxZVxl,mnX,mnY,mnZ) + call determineNumberOfEmptyVoxels(g,mxXVxl,mxYVxl,mxZVxl,nmEmpVxls) -reject = 0 -do cyc_cnt= 1, mcrescyc + allocate(listEmptyVoxels(int(nmEmpVxls*1.2),3),stat=error) - accept = .false. - call returnRandomVoxel(mxXVxl,mxYVxl,mxZVxl,tempX,tempY,tempZ) - call yesNoEmptymVoxel(g,tempX,tempY,tempZ,empty,rndmVxl) + if (stat.ne.0) then + write(0,*)"error: couldnt allocate memory for array, listEmptyVoxels=",stat + stop 1 + endif -if (empty .eqv. .true.) then + call updateEmptyVoxelsList(g,mxXVxl,mxYVxl,mxZVxl,nmEmpVxls,listEmptyVoxels) - call returnRandomWaterMolecule(crd,rndmWtr,mol_str,atm_cnt, mxXVxl, mxYVxl, mxZVxl) + call findMxCoarseGrid(mxX,mxY,mxZ,coarseGSP, maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl) - call plcWtrInVxl(frc,crd,rndmVxl,rndmWtr,X,Y,Z,mnX,mnY,mnZ,atm_cnt,atm_img_map,img_atm_map,& - mol_str,old_ene,nstep,mxXVxl,mxYVxl,mxZVxl,g,cnt,virial,ekcmt,reject,pme_err_est) + allocate(coarseGrid(maxCoarseXVxl+3, maxCoarseYVxl+3, maxCoarseZVxl+3, nint(coarseGSP*10)),stat=error) -end if + if (stat.ne.0) then + write(0,*)"error: couldnt allocate memory for array, coarseGrid=",stat + stop 1 + endif -end do + call initializeCoarseGrid(coarseGrid,maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl) + call buildCoarseGrid(crd,coarseGrid,atm_cnt,coarseGSP, maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl,mxX,mxY,mxZ) + + do cyc_cnt= 1, mcrescyc + + accept = .false. + call returnRandomVoxel(nmEmpVxls,listEmptyVoxels,rndmVxl) + coarseene=0.0 + call returnRandomWaterMolecule(crd,rndmWtr,mol_str,atm_cnt,mxXVxl,mxYVxl,mxZVxl, mxX,mxY,mxZ) + call calculateCoarseEnergy(atm_cnt, rndmWtr, crd, coarseGrid, coarseGSP,& + maxCoarseXVxl, maxCoarseYVxl, maxCoarseZVxl, coarseene,rndmWtr) + coarse_before_ene = coarseene + coarseene=0.0 + + ! This updates the water position to the new position. + call plcWtrInVxl(atm_cnt,crd,rndmVxl,rndmWtr,X,Y,Z,mol_str,mxXVxl,mxYVxl,mxZVxl,g,cnt, .true.) + + call updateCoarseGrid(atm_cnt,rndmWtr,crd,coarseGrid,maxCoarseXVxl, maxCoarseYVxl, maxCoarseZVxl, coarseGSP,& + & .false., mxX,mxY,mxZ,rndmWtr) + call calculateCoarseEnergy(atm_cnt, rndmWtr, crd, coarseGrid, coarseGSP,& + maxCoarseXVxl, maxCoarseYVxl, maxCoarseZVxl, coarseene,rndmWtr) + coarse_after_ene = coarseene + delta_coarse = coarse_after_ene - coarse_before_ene + + if(delta_coarse .le. 15) then + call mcresenergy(frc,crd,rndmVxl,rndmWtr,X,Y,Z,atm_cnt,atm_img_map,img_atm_map,& + mol_str,old_ene,nstep,mxXVxl,mxYVxl,mxZVxl,g,cnt,virial,ekcmt,pme_err_est,& + revert,nmEmpVxls,listEmptyVoxels) + else + do l = 0, numberOfAtoms-1 ! return the molecule + crd(1,cnt+l) = crd(1,cnt+l) + deltacrd(1) + crd(2,cnt+l) = crd(2,cnt+l) + deltacrd(2) + crd(3,cnt+l) = crd(3,cnt+l) + deltacrd(3) + end do -print *, cyc_cnt, 'energy cycle = ',old_ene + new_list = .true. + revert = .true. + end if + + if(revert) then + call updateCoarseGrid(atm_cnt, rndmWtr, crd, coarseGrid, maxCoarseXVxl,maxCoarseYVxl,maxCoarseZVxl,coarseGSP,& + & revert,mxX,mxY,mxZ,rndmWtr) + end if -print *,'the number of rejected energy calls is = ',reject + end do + + deallocate(listEmptyVoxels,stat=error) + + if (stat.ne.0) then + write(0,*)"error: couldnt allocate memory for array, listEmptyVoxels=",stat + stop 1 + endif deallocate(g,stat=error) if (error.ne.0) then - print *,'error in deallocating array g' + write(0,*)"error in deallocating array g" + stop 1 + endif + + deallocate(coarseGrid,stat=error) + if (error.ne.0) then + write(0,*)"error in deallocating array coarseGrid" + stop 1 endif end subroutine mcres diff --git src/pmemd/src/mdin_ctrl_dat.F90 src/pmemd/src/mdin_ctrl_dat.F90 index dd72c40..838025b 100644 --- src/pmemd/src/mdin_ctrl_dat.F90 +++ src/pmemd/src/mdin_ctrl_dat.F90 @@ -40,7 +40,7 @@ use file_io_dat_mod w_amd, emil_do_calc,isgld,isgsta,isgend,fixcom, & tishake, emil_sc,iemap, lj1264, efn,mcwat,mcint,& nucat, & ! gbneck2nu: check if atom belong to nucleic or not - mcrescyc,mcverbosity, & + mcrescyc,nmd,nmc,mcverbosity, & infe, & ! added by FENG PAN ineb, skmin, skmax, tmode, vv, & irandom, mask_from_ref, mbar_states, & @@ -75,11 +75,11 @@ use file_io_dat_mod mask_from_ref, mbar_states, & ! 96 geq_nstep, ginit_vel, gsyn_mass, & !99 gremd_acyc, gti_cpu_output, & !101 - mcint,mcrescyc,mcverbosity, & !104 - ineb,skmin,skmax,tmode,vv, & !109 - iphmd, midpoint !111 + mcint,nmd,mcrescyc,nmc,mcverbosity, & !106 + ineb,skmin,skmax,tmode,vv, & !111 + iphmd, midpoint !113 - integer, parameter :: mdin_ctrl_int_cnt = 111 + integer, parameter :: mdin_ctrl_int_cnt = 113 save :: / mdin_ctrl_int / @@ -302,7 +302,7 @@ use file_io_dat_mod efx,efy,efz,efn,effreq,efphase,infe,irandom, & mbar_lambda, mbar_states, mask_from_ref, & geq_nstep, ginit_vel, gsyn_mass, gremd_acyc, & - gti_cpu_output,mcwat,mcint,mcrescyc,mcverbosity,& + gti_cpu_output,mcwat,nmd,mcint,nmc,mcrescyc,mcverbosity,& mcresstr,mcwatmaxdif,ineb,skmin,skmax,tmode,vv,vfac,& mcboxshift, midpoint @@ -455,7 +455,9 @@ subroutine init_mdin_ctrl_dat(remd_method) mcverbosity=0 ! mcverbosity is a debug flag for mcwat mcwat=0 ! mcwat controls the Monte Carlo based water movement function. (set 1 to run) + nmd=0 mcint=0 ! mcint - interval for mcwat monte carlo moves. + nmc=0 mcrescyc=1 ! mcrescyc - controls number of mc moves per mcint - this is max attempts. If t gets ! an acceptance in less than mcrescyc it will break the loop and use that ! acceptance. @@ -1067,7 +1069,13 @@ subroutine validate_mdin_ctrl_dat(remd_method) write(mdout, '(a,a)') error_hdr, 'mcwat does not currently support MPI!' inerr = 1 #endif - if(mcint .lt. 1) then + if (mcrescyc .gt. 0) then + nmc=mcrescyc + end if + if (mcint .gt. 0) then + nmd=mcint + end if + if(nmd .lt. 1) then write(mdout, '(a,a)') error_hdr, 'mcint must be greater than 0!' inerr = 1 end if @@ -1075,7 +1083,7 @@ subroutine validate_mdin_ctrl_dat(remd_method) write(mdout, '(a,a)') error_hdr, 'mcwat=1 implies igb = 0!' inerr = 1 end if - if(mcrescyc .lt. 1) then + if(nmc .lt. 1) then write(mdout, '(a,a)') error_hdr, 'mcrescyc must be greater than 0!' inerr = 1 end if diff --git src/pmemd/src/runmd.F90 src/pmemd/src/runmd.F90 index e8e64d5..4be6b87 100644 --- src/pmemd/src/runmd.F90 +++ src/pmemd/src/runmd.F90 @@ -1357,18 +1357,21 @@ DBG_ARRAYS_TIME_STEP(nstep) ! Call mcres which will calculate energy with current coord set. ! attempts random residue move (solvent only) and then calculates ! new energy and applies metropolis criteria. - if (mcwat .gt. 0 .and. mod(nstep,mcint) .eq. 0) then + if (mcwat .gt. 0) then + mcint = nmd + if(mod(nstep,mcint) .eq. 0) then #ifdef MPI - call mcres(mcresaccept,mcresstr, atm_cnt, crd, frc, gbl_img_atm_map, & - gbl_atm_img_map, my_atm_lst, new_list, need_pot_enes, & - need_virials, pme_pot_ene, nstep, virial, ekcmt, & - pme_err_est) + call mcres(mcresaccept,mcresstr, atm_cnt, crd, frc, gbl_img_atm_map, & + gbl_atm_img_map, my_atm_lst, new_list, need_pot_enes, & + need_virials, pme_pot_ene, nstep, virial, ekcmt, & + pme_err_est) #else - call mcres(mcresaccept,mcresstr, atm_cnt, crd, frc, gbl_img_atm_map, & - gbl_atm_img_map, new_list, need_pot_enes, & - need_virials, pme_pot_ene, nstep, virial, ekcmt, & - pme_err_est) + call mcres(mcresaccept,mcresstr, atm_cnt, crd, frc, gbl_img_atm_map, & + gbl_atm_img_map, new_list, need_pot_enes, & + need_virials, pme_pot_ene, nstep, virial, ekcmt, & + pme_err_est) #endif + end if end if #ifdef MPI ! If we are doing REMD, there's no need to run pme_force on the first diff --git test/cuda/mcwat/Run.mcwat test/cuda/mcwat/Run.mcwat index b9c075f..d41c40a 100755 --- test/cuda/mcwat/Run.mcwat +++ test/cuda/mcwat/Run.mcwat @@ -66,17 +66,17 @@ $DO_PARALLEL $sander -O -i mdin -c eq1.x -o $output