*******> update.11
Author: Tyler Luchko
Date: December 17, 2019
Programs: sander, rism3d.snglpnt, nab
Description: Add a cluster size check for cluster-particle interaction in RISM treecode
Additional changes include
* Modifying cluster size check for TCF
* Adding more robust check for TCF and DCF treecode asymptotics
* modified syntax for passing writeVolume subroutine
* workarounds for PGI compiler
* fix for corner case in determining the k-space cutoff
--------------------------------------------------------------------------------
AmberTools/src/rism/amber_rism_interface.F90 | 114 ++++++-------------
AmberTools/src/rism/depend | 6 +-
AmberTools/src/rism/rism3d_c.F90 | 2 +-
AmberTools/src/rism/rism3d_potential_c.F90 | 18 ++-
AmberTools/src/rism/rism3d_tree_cp.F90 | 164 +++++++++++----------------
AmberTools/src/rism/rism_io.F90 | 20 ++++
6 files changed, 141 insertions(+), 183 deletions(-)
diff --git AmberTools/src/rism/amber_rism_interface.F90 AmberTools/src/rism/amber_rism_interface.F90
index ca22f9bc62..dcfffdd083 100644
--- AmberTools/src/rism/amber_rism_interface.F90
+++ AmberTools/src/rism/amber_rism_interface.F90
@@ -637,6 +637,9 @@ module amber_rism_interface
#endif /*RISM_CRDINTER*/
private :: rism_mpi_bcast
+
+
+
end module amber_rism_interface
#ifdef SANDER
@@ -2665,7 +2668,6 @@ contains
call rism_writeVolumetricData(rism_3d, step)
end subroutine rism_writeVolumetricDataC
-
!> Outputs volumetric data, such as solvent and electric potential
!! distributions, to their respective files. Each distribution
!! is written in a separate file with the step number before the
@@ -2676,6 +2678,7 @@ contains
subroutine rism_writeVolumetricData(this, step)
use constants, only : COULOMB_CONST_E, KB, PI
use amber_rism_interface
+ use rism_io
use rism3d_ccp4
use rism3d_opendx
use rism3d_xyzv
@@ -2683,18 +2686,6 @@ contains
use rism_util, only : checksum
implicit none
- !TODO: Move me to an I/O module.
- abstract interface
- subroutine writeVolumeInterface(file, data, grid, solute, o_rank, o_nproc, o_comm)
- use rism3d_grid_c
- use rism3d_solute_c
- character(len=*), intent(in) :: file
- type(rism3d_grid), intent(in) :: grid
- _REAL_, target, intent(in) :: data(grid%localDimsR(1), grid%localDimsR(2), grid%localDimsR(3))
- type(rism3d_solute), intent(in) :: solute
- integer, optional :: o_rank, o_nproc, o_comm
- end subroutine writeVolumeInterface
- end interface
#if defined(MPI)
include 'mpif.h'
@@ -2747,6 +2738,7 @@ contains
subroutine rism_writePdfTcfDcf(this, writeVolume, step, extension)
use constants, only : COULOMB_CONST_E, KB, PI
use amber_rism_interface
+ use rism_io
use rism3d_ccp4
use rism3d_opendx
use rism3d_xyzv
@@ -2754,20 +2746,10 @@ contains
use safemem
implicit none
- interface
- subroutine writeVolume(file, data, grid, solute, o_rank, o_nproc, o_comm)
- use rism3d_grid_c
- use rism3d_solute_c
- character(len=*), intent(in) :: file
- type(rism3d_grid), intent(in) :: grid
- _REAL_, target, intent(in) :: data(grid%localDimsR(1), grid%localDimsR(2), grid%localDimsR(3))
- type(rism3d_solute), intent(in) :: solute
- integer, optional :: o_rank, o_nproc, o_comm
- end subroutine writeVolume
- end interface
type(rism3d), intent(inout) :: this
character(len=*), intent(in) :: extension
integer, intent(in) :: step
+ procedure (writeVolumeInterface) :: writeVolume
character(len=64) :: suffix
character(len=16) :: cstep
@@ -2879,6 +2861,7 @@ contains
subroutine rism_writeSmearElectronMap(this, writeVolume, step, extension)
use constants, only : COULOMB_CONST_E, KB, PI
use amber_rism_interface
+ use rism_io
use rism3d_ccp4
use rism3d_opendx
use rism3d_xyzv
@@ -2887,20 +2870,10 @@ contains
use safemem
implicit none
- interface
- subroutine writeVolume(file, data, grid, solute, o_rank, o_nproc, o_comm)
- use rism3d_grid_c
- use rism3d_solute_c
- character(len=*), intent(in) :: file
- type(rism3d_grid), intent(in) :: grid
- _REAL_, target, intent(in) :: data(grid%localDimsR(1), grid%localDimsR(2), grid%localDimsR(3))
- type(rism3d_solute), intent(in) :: solute
- integer, optional :: o_rank, o_nproc, o_comm
- end subroutine writeVolume
- end interface
type(rism3d), intent(inout) :: this
character(len=*), intent(in) :: extension
integer, intent(in) :: step
+ procedure (writeVolumeInterface) :: writeVolume
character(len=64) :: suffix
character(len=16) :: cstep
@@ -2961,6 +2934,7 @@ contains
subroutine rism_writePotentialAsymptotics(this, writeVolume, step, extension)
use constants, only : COULOMB_CONST_E, KB, PI
use amber_rism_interface
+ use rism_io
use rism3d_ccp4
use rism3d_opendx
use rism3d_xyzv
@@ -2968,20 +2942,10 @@ contains
use safemem
implicit none
- interface
- subroutine writeVolume(file, data, grid, solute, o_rank, o_nproc, o_comm)
- use rism3d_grid_c
- use rism3d_solute_c
- character(len=*), intent(in) :: file
- type(rism3d_grid), intent(in) :: grid
- _REAL_, target, intent(in) :: data(grid%localDimsR(1), grid%localDimsR(2), grid%localDimsR(3))
- type(rism3d_solute), intent(in) :: solute
- integer, optional :: o_rank, o_nproc, o_comm
- end subroutine writeVolume
- end interface
type(rism3d), intent(inout) :: this
integer, intent(in) :: step
character(len=*), intent(in) :: extension
+ procedure (writeVolumeInterface) :: writeVolume
character(len=64) :: suffix
character(len=16) :: cstep
@@ -3038,6 +3002,7 @@ contains
subroutine rism_writeThermodynamicsDistributions(this, writeVolume, step, extension)
use constants, only : COULOMB_CONST_E, KB, PI
use amber_rism_interface
+ use rism_io
use rism3d_ccp4
use rism3d_opendx
use rism3d_xyzv
@@ -3045,20 +3010,10 @@ contains
use safemem
implicit none
- interface
- subroutine writeVolume(file, data, grid, solute, o_rank, o_nproc, o_comm)
- use rism3d_grid_c
- use rism3d_solute_c
- character(len=*), intent(in) :: file
- type(rism3d_grid), intent(in) :: grid
- _REAL_, target, intent(in) :: data(grid%localDimsR(1), grid%localDimsR(2), grid%localDimsR(3))
- type(rism3d_solute), intent(in) :: solute
- integer, optional :: o_rank, o_nproc, o_comm
- end subroutine writeVolume
- end interface
type(rism3d), intent(inout) :: this
character(len=*), intent(in) :: extension
integer, intent(in) :: step
+ procedure (writeVolumeInterface) :: writeVolume
character(len=64) :: suffix
character(len=16) :: cstep
@@ -3170,16 +3125,19 @@ contains
end if
! Outputting excess chemical potential map.
if (len_trim(excessChemicalPotentialfile) /= 0) then
- call writeThermo(excessChemicalPotentialfile, suffix, excessChemicalPotential_map, excessChemicalPotential_V_map)
+ call writeThermo(writeVolume,excessChemicalPotentialfile, suffix, excessChemicalPotential_map,&
+ excessChemicalPotential_V_map)
end if
! Outputting solvation energy map.
if (len_trim(solvationEnergyfile) /= 0 .and. rismprm%entropicDecomp == 1 .and. &
rism3d_canCalc_DT(rism_3d)) then
- call writeThermo(solvationEnergyfile, suffix, solvationEnergy_map, solvationEnergy_V_map)
+ call writeThermo(writeVolume,solvationEnergyfile, suffix, solvationEnergy_map,&
+ solvationEnergy_V_map)
end if
! Outputting solvent-solute energy map.
if (len_trim(solventPotentialEnergyfile) /= 0) then
- call writeThermo(solventPotentialEnergyfile, suffix, solventPotentialEnergy_map, solventPotentialEnergy_V_map)
+ call writeThermo(writeVolume,solventPotentialEnergyfile, suffix, solventPotentialEnergy_map, &
+ solventPotentialEnergy_V_map)
end if
! Outputting entropy map.
if (len_trim(entropyfile) /= 0 .and. rismprm%entropicDecomp == 1 .and. &
@@ -3187,7 +3145,8 @@ contains
call DAXPY(this%grid%totalLocalPointsR, -1d0, solvationEnergy_map, 1, excessChemicalPotential_map, 1)
call DAXPY(this%grid%totalLocalPointsR*this%solvent%numAtomTypes, -1d0, &
solvationEnergy_V_map, 1, excessChemicalPotential_V_map, 1)
- call writeThermo(entropyfile, suffix, excessChemicalPotential_map, excessChemicalPotential_V_map)
+ call writeThermo(writeVolume,entropyfile, suffix, excessChemicalPotential_map,&
+ excessChemicalPotential_V_map)
end if
@@ -3311,12 +3270,14 @@ contains
end if
! Outputting excess chemical potential map.
if (len_trim(excessChemicalPotentialGFfile) /= 0) then
- call writeThermo(excessChemicalPotentialGFfile, suffix, excessChemicalPotential_map, excessChemicalPotential_V_map)
+ call writeThermo(writeVolume,excessChemicalPotentialGFfile, suffix, excessChemicalPotential_map, &
+ excessChemicalPotential_V_map)
end if
! Outputting solvation energy map.
if (len_trim(solvationEnergyGFfile) /= 0 .and. rismprm%entropicDecomp == 1 .and. &
rism3d_canCalc_DT(rism_3d)) then
- call writeThermo(solvationEnergyGFfile, suffix, solvationEnergy_map, solvationEnergy_V_map)
+ call writeThermo(writeVolume,solvationEnergyGFfile, suffix, solvationEnergy_map, &
+ solvationEnergy_V_map)
end if
! Outputting entropy map.
if (len_trim(entropyGFfile) /= 0 .and. rismprm%entropicDecomp == 1 .and. &
@@ -3324,7 +3285,8 @@ contains
call DAXPY(this%grid%totalLocalPointsR, -1d0, solvationEnergy_map, 1, excessChemicalPotential_map, 1)
call DAXPY(this%grid%totalLocalPointsR*this%solvent%numAtomTypes, -1d0, &
solvationEnergy_V_map, 1, excessChemicalPotential_V_map, 1)
- call writeThermo(entropyGFfile, suffix, excessChemicalPotential_map, excessChemicalPotential_V_map)
+ call writeThermo(writeVolume,entropyGFfile, suffix, excessChemicalPotential_map, &
+ excessChemicalPotential_V_map)
end if
if (rismprm%molReconstruct) then
call rism_writeMolReconstruct(this, writeVolume, suffix)
@@ -3345,14 +3307,19 @@ contains
if (safemem_dealloc(solventPotentialEnergy_V_map)/= 0) &
call rism_report_error("RISM_WRITESOLVEDIST: failed to deallocate SOLVENTPOTENTIALENERGY_V_MAP")
contains
- subroutine writeThermo(root, suffix, data, data_sites)
+ subroutine writeThermo(wv, root, suffix, data, data_sites)
implicit none
+ ! we are passing writeVolume to work around what appears to be a
+ ! PGI compiler bug. writeThermo has scope to see writeVolume so
+ ! we should be able to use it directly. Doing so works with GNU
+ ! and Intel but not PGI
+ procedure(writeVolumeInterface) :: wv
character(len=*), intent(in) :: root, suffix
_REAL_, intent(in) :: data(:,:,:), data_sites(:,:,:,:)
- call writeVolume(trim(root)//trim(suffix), data, this%grid, &
+ call wv(trim(root)//trim(suffix), data, this%grid, &
this%solute, mpirank, mpisize, mpicomm)
do iv = 1, this%solvent%numAtomTypes
- call writeVolume(trim(root)//'.'//trim(rism_3d%solvent%atomName(iv))//trim(suffix), &
+ call wv(trim(root)//'.'//trim(rism_3d%solvent%atomName(iv))//trim(suffix), &
data_sites(:, :, :, iv), this%grid, &
this%solute, mpirank, mpisize, mpicomm)
end do
@@ -3367,6 +3334,7 @@ contains
subroutine rism_writeMolReconstruct(this, writeVolume, suffix)
use constants, only : COULOMB_CONST_E, KB, PI
use amber_rism_interface
+ use rism_io
use rism3d_ccp4
use rism3d_opendx
use rism3d_xyzv
@@ -3374,19 +3342,9 @@ contains
use safemem
implicit none
- interface
- subroutine writeVolume(file, data, grid, solute, o_rank, o_nproc, o_comm)
- use rism3d_grid_c
- use rism3d_solute_c
- character(len=*), intent(in) :: file
- type(rism3d_grid), intent(in) :: grid
- _REAL_, target, intent(in) :: data(grid%localDimsR(1), grid%localDimsR(2), grid%localDimsR(3))
- type(rism3d_solute), intent(in) :: solute
- integer, optional :: o_rank, o_nproc, o_comm
- end subroutine writeVolume
- end interface
type(rism3d), intent(inout) :: this
character(len=64), intent(in) :: suffix
+ procedure (writeVolumeInterface) :: writeVolume
_REAL_, pointer::excessChemicalPotential_V_map(:, :, :, :) => NULL(), &
solvationEnergy_V_map(:, :, :, :) => NULL(), &
solventPotentialEnergy_V_map(:, :, :, :) => NULL(),&
diff --git AmberTools/src/rism/depend AmberTools/src/rism/depend
index c7e18a43ee..bdffb2edb8 100644
--- AmberTools/src/rism/depend
+++ AmberTools/src/rism/depend
@@ -18,7 +18,6 @@ amber_rism_interface.o: \
rism3d_ccp4.o\
rism3d_opendx.o\
rism3d_xyzv.o\
- rism3d_grid_c.o\
files.h\
array_util.o
@@ -46,7 +45,6 @@ amber_rism_interface.NAB.o: \
rism3d_ccp4.o\
rism3d_opendx.o\
rism3d_xyzv.o\
- rism3d_grid_c.o\
files.h\
array_util.o
@@ -74,7 +72,6 @@ amber_rism_interface.SANDER.o: \
rism3d_ccp4.o\
rism3d_opendx.o\
rism3d_xyzv.o\
- rism3d_grid_c.o\
files.h\
array_util.o
@@ -623,7 +620,6 @@ rism3d_closure_c.o: \
safemem.o\
rism_util.o\
constants.o\
- rism3d_opendx.o\
bspline.o\
fftw3.o
@@ -834,6 +830,8 @@ rism_io.o: \
../include/dprec.fh\
safemem.o\
rism_report_c.o\
+ rism3d_grid_c.o\
+ rism3d_solute_c.o\
rism_util.o
diff --git AmberTools/src/rism/rism3d_c.F90 AmberTools/src/rism/rism3d_c.F90
index 24e2d30d44..6b02df8d6f 100644
--- AmberTools/src/rism/rism3d_c.F90
+++ AmberTools/src/rism/rism3d_c.F90
@@ -519,7 +519,7 @@ contains
this%centering = centering
this%ncuvsteps = ncuvsteps
nclosure = size(closure)
- t_closure => safemem_realloc(t_closure, len(t_closure), nclosure)
+ t_closure => safemem_realloc(t_closure, len(closure), nclosure)
t_closure = closure
t_cut = cut
t_mdiis_nvec = mdiis_nvec
diff --git AmberTools/src/rism/rism3d_potential_c.F90 AmberTools/src/rism/rism3d_potential_c.F90
index bb3ffa5f92..a52798b769 100644
--- AmberTools/src/rism/rism3d_potential_c.F90
+++ AmberTools/src/rism/rism3d_potential_c.F90
@@ -77,7 +77,11 @@ module asympk_cut
! maxval(abs(soluteCharges)), FOURPI, sqrt(2d0), maxval(abs(solventCharges)), &
! sum(abs(soluteCharges))/size(soluteCharges)
smearSq = smear**2
- cutoff = root_newton(asympck_cut, asympck_cut_deriv, 10d0, 1d-16)
+ ! We need to that we approach the
+ ! root from the left as approaching the root from the right can
+ ! cause the Newton to overshoot and end-up with NaNs. In the
+ ! case of a NaN cutoff, the effect is to have no cutoff.
+ cutoff = root_newton(asympck_cut, asympck_cut_deriv, 1d0, 1d-16)
end function asympck_cut_calc
!> Expression to find the k-space cut off for the DCF asymptotics.
@@ -514,6 +518,18 @@ contains
else
this%cut2_chlk = asympck_cut_calc(this%solute%charge, this%solvent%charge,&
this%chargeSmear, boxVolume, asympKSpaceTolerance)
+ ! isnan is not standard, but GNU doesn't support Fortran 2003's
+ ! ieee_arithmetic until version 5. Intel 15+ supports both; use
+ ! the PGI define to enable the most support for old compilers.
+#ifdef __PGI
+ if (ieee_is_nan(this%cut2_chlk)) then
+#else
+ if (isnan(this%cut2_chlk)) then
+#endif
+ call rism_report_error('Could not converge k-space asymptotics cutoff. '&
+ //'Try using a smaller error tolerance or no cutoff.')
+
+ end if
end if
end subroutine rism3d_potential_setcut_asympktolerance
diff --git AmberTools/src/rism/rism3d_tree_cp.F90 AmberTools/src/rism/rism3d_tree_cp.F90
index ea76604a82..781687f951 100644
--- AmberTools/src/rism/rism3d_tree_cp.F90
+++ AmberTools/src/rism/rism3d_tree_cp.F90
@@ -1,36 +1,55 @@
!
-! Authors: Leighton W. Wilson (lwwilson@umich.edu)
-! Henry A. Boateng (boateng@umich.edu)
+! Author: Leighton W. Wilson (lwwilson@umich.edu)
+! Department of Mathematics
+! University of Michigan, Ann Arbor
!
-! Department of Mathematics
-! University of Michigan, Ann Arbor
+! The cluster-particle treecode software found here is copyright
+! (c) 2013-2019 by the Regents of the University of Michigan.
!
-! Copyright (c) 2013-2017. The Regents of the University of Michigan.
-! All Rights Reserved.
+! The 3D-RISM-KH software found here is copyright (c) 2011-2012 by
+! Andriy Kovalenko, Tyler Luchko and David A. Case.
!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 3 of the License, or
-! (at your option) any later version.
+! This program is free software: you can redistribute it and/or modify it
+! under the terms of the GNU General Public License as published by the Free
+! Software Foundation, either version 3 of the License, or (at your option)
+! any later version.
!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
+! This program is distributed in the hope that it will be useful, but
+! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+! for more details.
!
-! You should have received a copy of the GNU General Public License
-! along with this program; if not, .
+! You should have received a copy of the GNU General Public License in the
+! ../../LICENSE file. If not, see .
!
-!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+! Users of the 3D-RISM capability found here are requested to acknowledge
+! use of the software in reports and publications. Such acknowledgement
+! should include the following citations:
+!
+! 1) A. Kovalenko and F. Hirata. J. Chem. Phys., 110:10095-10112 (1999);
+! ibid. 112:10391-10417 (2000).
+!
+! 2) A. Kovalenko, in: Molecular Theory of Solvation, edited by
+! F. Hirata (Kluwer Academic Publishers, Dordrecht, 2003), pp.169-275.
+!
+! 3) T. Luchko, S. Gusarov, D.R. Roe, C. Simmerling, D.A. Case, J. Tuszynski,
+! and A. Kovalenko, J. Chem. Theory Comput., 6:607-624 (2010).
+!
+! 4) I. Omelyan and A. Kovalenko, Mol. Simul. 39:25-48 (2013).
+
#include "../include/dprec.fh"
+!> Module for treecode calculations in 3D-RISM. Used to calculate
+!! the Coulomb potential and the real space long-range asymptotics
+!! for the TCF and DCF.
+!!
+!! Additionally contains routines for treecode-accelerated 3D-RISM
+!! force calculations, if the interest arises in the future.
MODULE treecode_procedures
USE safemem
IMPLICIT NONE
-! r8 is 8-byte (double precision) real
-
INTEGER,PARAMETER :: r8=SELECTED_REAL_KIND(15)
_REAL_,PARAMETER :: pi = 4_r8*ATAN(1.0_r8)
@@ -59,8 +78,8 @@
INTEGER :: numpar
_REAL_,DIMENSION(3) :: xyz_min, xyz_max, xyz_mid
- _REAL_ :: radius,sqradius,aspect
- INTEGER :: level,num_children,exist_ms
+ _REAL_ :: radius, sqradius, aspect, min_len
+ INTEGER :: level, num_children, exist_ms
_REAL_,POINTER :: ms(:) => NULL()
TYPE(tnode_pointer) :: child(8)
@@ -160,7 +179,7 @@
! local variables
_REAL_, DIMENSION(3) :: xyz_len
- _REAL_ :: lmax,t2
+ _REAL_ :: lmax
INTEGER :: i, err, loclev, numposchild
_REAL_, DIMENSION(6,8) :: xyzmms
@@ -205,10 +224,10 @@
xyz_len=p%xyz_max-p%xyz_min
lmax=MAXVAL(xyz_len)
- t2=MINVAL(xyz_len)
+ p%min_len=MINVAL(xyz_len)
- IF (t2 .NE. 0.0_r8) THEN
- p%aspect=lmax/t2
+ IF (p%min_len .GT. 0.0_r8) THEN
+ p%aspect=lmax/p%min_len
ELSE
p%aspect=0.0_r8
END IF
@@ -399,9 +418,6 @@
! source particle, setting the global variable TARPOS before the call. After
! the calls to COMPUTE_CP1, CP_TREECODE calls COMPUTE_CP2 to compute
! the energy using power series.
-! P is the root node of the target tree, xT,yT,zT are the coordinates of
-! the target particles while xS,yS,zS,qS are the coordinates and charges of the
-! source particles. The energy at target 'i', is stored in EnP(i).
!
INTEGER,INTENT(IN) :: numparsS,numparsT
@@ -447,9 +463,6 @@
! source particle, setting the global variable TARPOS before the call. After
! the calls to COMPUTE_CP1, CP_TREECODE calls COMPUTE_CP2 to compute
! the energy using power series.
-! P is the root node of the target tree, xT,yT,zT are the coordinates of
-! the target particles while xS,yS,zS,qS are the coordinates and charges of the
-! source particles. The energy at target 'i', is stored in EnP(i).
!
INTEGER,INTENT(IN) :: numparsS,numparsT
@@ -493,9 +506,6 @@
! source particle, setting the global variable TARPOS before the call. After
! the calls to COMPUTE_CP1, CP_TREECODE calls COMPUTE_CP2 to compute
! the energy using power series.
-! P is the root node of the target tree, xT,yT,zT are the coordinates of
-! the target particles while xS,yS,zS,qS are the coordinates and charges of the
-! source particles. The energy at target 'i', is stored in EnP(i).
!
INTEGER,INTENT(IN) :: numparsS,numparsT
@@ -553,20 +563,23 @@
! local variables
_REAL_,DIMENSION(3) :: xyz_t
- _REAL_ :: distsq
+ _REAL_ :: distsq, dist
INTEGER :: i, err
! determine DISTSQ for MAC test
xyz_t=tarpos-p%xyz_mid
distsq=SUM(xyz_t**2)
+ dist=SQRT(distsq)
! intialize potential energy and force
! If MAC is accepted and there is more than 1 particle in the
! box use the expansion for the approximation.
- IF ((p%sqradius .LT. distsq*thetasq) .AND. &
- (p%sqradius .NE. 0.0_r8)) THEN
+ IF (( p%sqradius .LT. distsq * thetasq) .AND. &
+ ( p%sqradius .GT. 0.0_r8 ) .AND. &
+ (dist-p%min_len .GE. 5.0_r8 * eta ) .AND. &
+ ( torderflat .LE. 4 * p%numpar )) THEN
CALL COMP_TCOEFF_TCF(xyz_t(1),xyz_t(2),xyz_t(3), kappa)
@@ -619,20 +632,23 @@
! local variables
_REAL_,DIMENSION(3) :: xyz_t
- _REAL_ :: distsq
+ _REAL_ :: distsq, dist
INTEGER :: i, err
! determine DISTSQ for MAC test
xyz_t=tarpos-p%xyz_mid
distsq=SUM(xyz_t**2)
+ dist=SQRT(distsq)
! intialize potential energy and force
! If MAC is accepted and there is more than 1 particle in the
! box use the expansion for the approximation.
- IF ((p%sqradius .LT. distsq*thetasq) .AND. &
- (p%sqradius .NE. 0.0_r8)) THEN
+ IF (( p%sqradius .LT. distsq * thetasq) .AND. &
+ ( p%sqradius .GT. 0.0_r8 ) .AND. &
+ (dist-p%min_len .GE. 5.0_r8 * eta ) .AND. &
+ ( torderflat .LE. p%numpar )) THEN
CALL COMP_TCOEFF_DCF(xyz_t(1),xyz_t(2),xyz_t(3), eta)
@@ -697,7 +713,7 @@
! If MAC is accepted and there is more than 1 particle in the
! box use the expansion for the approximation.
IF ((p%sqradius .LT. distsq*thetasq) .AND. &
- (p%sqradius .NE. 0.0_r8)) THEN
+ (p%sqradius .GT. 0.0_r8)) THEN
CALL COMP_TCOEFF_COULOMB(xyz_t(1),xyz_t(2),xyz_t(3))
@@ -1570,13 +1586,10 @@
forcesT) !output
IMPLICIT NONE
!
-! CP_TREECODE is the driver routine which calls COMPUTE_CP1 for each
-! source particle, setting the global variable TARPOS before the call. After
-! the calls to COMPUTE_CP1, CP_TREECODE calls COMPUTE_CP2 to compute
-! the energy using power series.
-! P is the root node of the target tree, xT,yT,zT are the coordinates of
-! the target particles while xS,yS,zS,qS are the coordinates and charges of the
-! source particles. The energy at target 'i', is stored in EnP(i).
+! PC_TREECODE_FORCES is the driver routine which calls COMPUTE_PC_FORCES for each
+! source particle, setting the global variable TARPOS before the call.
+! P is the root node of the source tree, solutePosition are the coordinates of
+! the target particles, and guv contains source grid information.
!
INTEGER,INTENT(IN) :: numparsS,numparsT
@@ -1602,13 +1615,9 @@
tarposq = soluteCharge(i)
CALL COMPUTE_PC_FORCES(p, guv, numparsS, force)
- !DO j=1,p%num_children
- ! CALL COMPUTE_PC(p%child(j)%p_to_tnode,guv,numparsS,forcetemp)
- ! force = force + forcetemp
- !END DO
! We have to do this backwards because we invert our grid directions to
-! agree with RISM's directions. A kludge indeed.
+! agree with RISM's directions
forcesT(1,i) = -tarposq * force(1)
forcesT(2,i) = -tarposq * force(2)
forcesT(3,i) = -tarposq * force(3)
@@ -1625,7 +1634,7 @@
RECURSIVE SUBROUTINE COMPUTE_PC_FORCES(p,guv,numpars,force)
IMPLICIT NONE
-! COMPUTE_CP1 is the recursive routine for computing the interaction
+! COMPUTE_PC_FORCES is the recursive routine for computing the interaction
! between a target particle and a source cluster. If the MAC is
! satisfied the power series coefficients for the current target
! cluster are updated. If the MAC is not satisfied then the algorithm
@@ -1778,8 +1787,8 @@
IMPLICIT NONE
!
! COMP_DIRECT directly computes the potential on the targets
-! in the current cluster due to the
-! current source (determined by the global variable TARPOS).
+! in the current cluster due to the current source
+! (determined by the global variable TARPOS).
!
INTEGER,INTENT(IN) :: numpars
TYPE(tnode),POINTER :: p
@@ -1969,28 +1978,6 @@
targetGridLimits, targetGridDims, targetGridIndices, &
treeLevel)
-
-! print tree information to stdout
-
-! OPEN(unit = 11, file ="test.txt",action='write',position='append')
-! WRITE(11,*) ' '
-! WRITE(11,*) 'Tree created. '
-! WRITE(11,*) 'Tree parameters: '
-! WRITE(11,*) ' '
-! WRITE(11,*) ' numpar: ',trootT%numpar
-! WRITE(11,*) ' x_mid: ',trootT%xyz_mid(1)
-! WRITE(11,*) ' y_mid: ',trootT%xyz_mid(2)
-! WRITE(11,*) ' z_mid: ',trootT%xyz_mid(3)
-! WRITE(11,*) ' radius: ',trootT%radius
-! WRITE(11,*) ' torder: ',torder
-! WRITE(11,*) ' theta: ',theta
-! WRITE(11,*) ' maxparnode: ',maxparnodeT
-! WRITE(11,*) ' '
-! WRITE(11,*) 'Tree will be run on ', numparsS, ' sources.'
-! WRITE(11,*) ' '
-! CLOSE(11)
-
-
!Call driver routine for cluster-particle
IF (potentialType == 0) THEN
CALL CP_TREECODE_TCF(treeRootNode, solutePosition, soluteCharge, potentialGrid, &
@@ -2003,7 +1990,6 @@
numSources, numTargets)
END IF
-
CALL CLEANUP(treeRootNode)
END SUBROUTINE TREECODE
@@ -2072,26 +2058,6 @@
sourceGridLimits, sourceGridDims, sourceGridIndices, &
treeLevel)
-! print tree information to stdout
-! OPEN(unit = 10, file ="test.txt",action='write',position='append')
-! WRITE(10,*) ' '
-! WRITE(10,*) 'Tree created. '
-! WRITE(10,*) 'Tree parameters: '
-! WRITE(10,*) ' '
-! WRITE(10,*) ' numpar: ',trootS%numpar
-! WRITE(10,*) ' x_mid: ',trootS%xyz_mid(1)
-! WRITE(10,*) ' y_mid: ',trootS%xyz_mid(2)
-! WRITE(10,*) ' z_mid: ',trootS%xyz_mid(3)
-! WRITE(10,*) ' radius: ',trootS%radius
-! WRITE(10,*) ' torder: ',torder
-! WRITE(10,*) ' theta: ',theta
-! WRITE(10,*) ' maxparnode: ',maxparnodeS
-! WRITE(10,*) ' levels: ',level
-! WRITE(10,*) ' '
-! WRITE(10,*) ' Tree will be run on ', numparsT, ' targets.'
-! WRITE(10,*)
-! CLOSE(10)
-
!Call driver routine for particle-cluster forces
CALL PC_TREECODE_FORCES(treeRootNode, guv, & !source info
solutePosition, soluteCharge, & !target info
diff --git AmberTools/src/rism/rism_io.F90 AmberTools/src/rism/rism_io.F90
index ae73e8d67f..02e1a17887 100644
--- AmberTools/src/rism/rism_io.F90
+++ AmberTools/src/rism/rism_io.F90
@@ -39,6 +39,26 @@ module rism_io
use rism_report_c
implicit none
+ !> Abstract interace for functions that write volumetric data
+ !! @param[in] file Name of output file
+ !! @param[in] data 3D array of _REAL_ data to write
+ !! @param[in] grid rism3d_grid object for data
+ !! @param[in] solute rism3d_solute object for data
+ !! @param[in] o_rank (optional) MPI rank
+ !! @param[in] o_nproc (optional) Number of MPI proceses
+ !! @param[in] o_comm (optional) MPI communicator
+ abstract interface
+ subroutine writeVolumeInterface(file, data, grid, solute, o_rank, o_nproc, o_comm)
+ use rism3d_grid_c
+ use rism3d_solute_c
+ character(len=*), intent(in) :: file
+ type(rism3d_grid), intent(in) :: grid
+ _REAL_, target, intent(in) :: data(grid%localDimsR(1), grid%localDimsR(2), grid%localDimsR(3))
+ type(rism3d_solute), intent(in) :: solute
+ integer, optional :: o_rank, o_nproc, o_comm
+ end subroutine writeVolumeInterface
+ end interface
+
contains
!> Read a 1D RDF from a space separated variable file.