Skip to content

Commit

Permalink
PINT: corrected helium RDFs
Browse files Browse the repository at this point in the history
  • Loading branch information
@cschran @hforbert
cschran authored and hforbert committed Sep 2, 2019
1 parent 87d91fc commit 98a0cf8
Show file tree
Hide file tree
Showing 9 changed files with 232 additions and 182 deletions.
174 changes: 99 additions & 75 deletions src/motion/helium_common.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,12 @@ MODULE helium_common
helium_cell_shape_octahedron
USE kinds, ONLY: default_string_length,&
dp
USE mathconstants, ONLY: pi
USE memory_utilities, ONLY: reallocate
USE parallel_rng_types, ONLY: next_random_number
USE physcon, ONLY: angstrom,&
bohr
USE pint_public, ONLY: pint_calc_centroid,&
pint_com_pos
USE pint_public, ONLY: pint_com_pos
USE pint_types, ONLY: pint_env_type
USE splines_methods, ONLY: spline_value
USE splines_types, ONLY: spline_data_p_type,&
Expand Down Expand Up @@ -713,108 +713,130 @@ END SUBROUTINE helium_norm_rho
!> \brief Calculate helium radial distribution functions wrt positions given
!> by <centers> using the atomic weights given by <weights>.
!> \param helium ...
!> \param centers ...
!> \date 2009-07-22
!> \par History
!> 2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran]
!> \author Lukasz Walewski
! **************************************************************************************************
SUBROUTINE helium_calc_rdf(helium, centers)
SUBROUTINE helium_calc_rdf(helium)

TYPE(helium_solvent_type), POINTER :: helium
REAL(KIND=dp), DIMENSION(:), POINTER :: centers

CHARACTER(len=*), PARAMETER :: routineN = 'helium_calc_rdf', &
routineP = moduleN//':'//routineN

CHARACTER(len=default_string_length) :: msgstr
INTEGER :: bin, handle, ia, ib, ic, id, &
INTEGER :: bin, handle, ia, ib, ic, ind_hehe, &
n_out_of_range, nbin
REAL(KIND=dp) :: const, invdr, invp, nideal, pi, ri, &
rlower, rupper, sdens
REAL(KIND=dp) :: invdr, nideal, ri, rlower, rupper
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pref
REAL(KIND=dp), DIMENSION(3) :: r, r0
REAL(KIND=dp), DIMENSION(:), POINTER :: incr

CALL timeset(routineN, handle)

pi = 4.0_dp*ATAN(1.0_dp)
invdr = 1.0_dp/helium%rdf_delr
invp = 1.0_dp/helium%beads
nbin = helium%rdf_nbin
sdens = helium%wpref*invp*helium%atoms
ALLOCATE (incr(helium%rdf_num))
incr(:) = 0.0d0

! calculate the histogram of distances
n_out_of_range = 0
helium%rdf_inst(:, :, :) = 0.0_dp
DO ic = 1, SIZE(centers)/3
helium%rdf_inst(:, :) = 0.0_dp

r0(1) = centers(3*(ic-1)+1)
r0(2) = centers(3*(ic-1)+2)
r0(3) = centers(3*(ic-1)+3)
DO ia = 1, helium%atoms
ind_hehe = 0
! Calculate He-He RDF
IF (helium%rdf_he_he) THEN
ind_hehe = 1
DO ib = 1, helium%beads
DO ia = 1, helium%atoms
r0(:) = helium%pos(:, ia, ib)

DO ic = 1, helium%atoms
IF (ia == ic) CYCLE
r(:) = helium%pos(:, ic, ib)-r0(:)
CALL helium_pbc(helium, r)
ri = SQRT(r(1)*r(1)+r(2)*r(2)+r(3)*r(3))
bin = INT(ri*invdr)+1
IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN
! increment the RDF value for He atoms inside the r_6 sphere
helium%rdf_inst(ind_hehe, bin) = helium%rdf_inst(ind_hehe, bin)+1.0_dp
ELSE
n_out_of_range = n_out_of_range+1
END IF
END DO
END DO
END DO
END IF

! set the increments for this atom
incr(1) = invp

DO ib = 1, helium%beads
r(:) = helium%pos(:, ia, ib)-r0(:)
CALL helium_pbc(helium, r)
ri = SQRT(r(1)*r(1)+r(2)*r(2)+r(3)*r(3))
bin = INT(ri*invdr)+1
IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN
! increment the RDF value for He atoms inside the r_6 sphere
DO id = 1, helium%rdf_num
helium%rdf_inst(id, bin, ic) = helium%rdf_inst(id, bin, ic)+incr(id)
END DO
ELSE
!WRITE(msgstr,*) helium%center * angstrom
!msgstr = "center = "//TRIM(ADJUSTL(msgstr))
!CALL cp_error_message(cp_warning_level, routineP, msgstr, error)
!WRITE(msgstr,*) r0 * angstrom
!msgstr = "r0 = "//TRIM(ADJUSTL(msgstr))
!CALL cp_error_message(cp_warning_level, routineP, msgstr, error)
!WRITE(msgstr,*) helium%pos(:,ia,ib) * angstrom
!msgstr = "pos = "//TRIM(ADJUSTL(msgstr))
!CALL cp_error_message(cp_warning_level, routineP, msgstr, error)
!WRITE(msgstr,*) ri * angstrom
!msgstr = "ri = "//TRIM(ADJUSTL(msgstr))
!CALL cp_error_message(cp_warning_level, routineP, msgstr, error)
!WRITE(msgstr,*) ri * invdr
!msgstr = "ri/dr = "//TRIM(ADJUSTL(msgstr))
!CALL cp_error_message(cp_warning_level, routineP, msgstr, error)
!WRITE(msgstr,*) helium%current_step
!msgstr = "step = "//TRIM(ADJUSTL(msgstr))
!CALL cp_error_message(cp_warning_level, routineP, msgstr, error)
n_out_of_range = n_out_of_range+1
END IF
! Calculate Solute-He RDF
IF (helium%solute_present .AND. helium%rdf_sol_he) THEN
DO ib = 1, helium%beads
DO ia = 1, helium%solute_atoms
r0(:) = helium%rdf_centers(ib, 3*(ia-1)+1:3*(ia-1)+3)

DO ic = 1, helium%atoms
r(:) = helium%pos(:, ic, ib)-r0(:)
CALL helium_pbc(helium, r)
ri = SQRT(r(1)*r(1)+r(2)*r(2)+r(3)*r(3))
bin = INT(ri*invdr)+1
IF ((0 .LT. bin) .AND. (bin .LE. nbin)) THEN
! increment the RDF value for He atoms inside the r_6 sphere
helium%rdf_inst(ind_hehe+ia, bin) = helium%rdf_inst(ind_hehe+ia, bin)+1.0_dp
ELSE
n_out_of_range = n_out_of_range+1
END IF
END DO
END DO
END DO
END DO

END IF

! for periodic case we intentionally truncate RDFs to spherical volume
! so we skip atoms in the corners
IF (.NOT. helium%periodic) THEN
IF (n_out_of_range .GT. 0) THEN
WRITE (msgstr, *) n_out_of_range
msgstr = "Number of bead positions out of range: "//TRIM(ADJUSTL(msgstr))
CPABORT(msgstr)
CPWARN(msgstr)
END IF
END IF
! for periodic case we intentionally truncate RDFs to spherical volume
! so we skip atoms in the corners

! normalize the histogram to get g(r)
! note: helium%density refers to the number of atoms, not the beads
const = 4.0_dp*pi*helium%density/3.0_dp
ALLOCATE (pref(helium%rdf_num))
pref(:) = 0.0_dp
IF (helium%periodic) THEN
! Use helium density for normalization
pref(:) = helium%density*helium%beads*helium%atoms
! Correct for He-He-RDF
IF (helium%rdf_he_he) THEN
pref(1) = pref(1)/helium%atoms*(helium%atoms-1)
END IF
ELSE
! Non-periodic case has density of 0, use integral for normalzation
! This leads to a unit of 1/volume of the RDF
pref(:) = 0.5_dp*helium%rdf_inst(:, 1)
DO bin = 2, helium%rdf_nbin-1
pref(:) = pref(:)+helium%rdf_inst(:, bin)
END DO
pref(:) = pref(:)+0.5_dp*helium%rdf_inst(:, helium%rdf_nbin)

!set integral of histogram to number of atoms:
pref(:) = pref(:)/helium%atoms
! Correct for He-He-RDF
IF (helium%rdf_he_he) THEN
pref(1) = pref(1)*helium%atoms/(helium%atoms-1)
END IF
END IF
! Volume integral first:
DO bin = 1, helium%rdf_nbin
rlower = REAL(bin-1, dp)*helium%rdf_delr
rupper = rlower+helium%rdf_delr
nideal = const*(rupper**3-rlower**3)
DO id = 1, helium%rdf_num
helium%rdf_inst(id, bin, :) = helium%rdf_inst(id, bin, :)/nideal
END DO
nideal = (rupper**3-rlower**3)*4.0_dp*pi/3.0_dp
helium%rdf_inst(:, bin) = helium%rdf_inst(:, bin)/nideal
END DO
! No normalization for density
pref(:) = 1.0_dp/pref(:)
DO ia = 1, helium%rdf_num
helium%rdf_inst(ia, :) = helium%rdf_inst(ia, :)*pref(ia)
END DO

DEALLOCATE (incr)
NULLIFY (incr)
DEALLOCATE (pref)

CALL timestop(handle)

Expand Down Expand Up @@ -2349,6 +2371,8 @@ END SUBROUTINE helium_update_coord_system
!> \param helium ...
!> \param pint_env ...
!> \date 2014-04-25
!> \par History
!> 2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran]
!> \author Lukasz Walewski
!> \note Sets the centers wrt which radial distribution functions are
!> calculated.
Expand All @@ -2361,14 +2385,14 @@ SUBROUTINE helium_set_rdf_coord_system(helium, pint_env)
CHARACTER(len=*), PARAMETER :: routineN = 'helium_set_rdf_coord_system', &
routineP = moduleN//':'//routineN
INTEGER :: i
INTEGER :: i, j
IF (helium%solute_present) THEN
CALL pint_calc_centroid(pint_env)
i = 3*helium%solute_atoms
helium%rdf_centers(1:i) = pint_env%centroid(:)
ELSE
helium%rdf_centers(:) = helium%center(:)
IF (helium%solute_present .AND. helium%rdf_sol_he) THEN
! Account for unequal number of beads for solute and helium
DO i = 1, helium%beads
j = ((i-1)*helium%solute_beads)/helium%beads+1
helium%rdf_centers(i, :) = pint_env%x(j, :)
END DO
END IF
RETURN
Expand Down
71 changes: 48 additions & 23 deletions src/motion/helium_io.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -413,7 +413,7 @@ SUBROUTINE helium_write_setup(helium)
! RDF related settings
IF (helium%rdf_present) THEN
WRITE (stmp, '(I6)') helium%rdf_num_ctr
WRITE (stmp, '(I6)') helium%rdf_num
CALL helium_write_line("RDF| number of centers: "//TRIM(stmp))
rtmp = cp_unit_from_cp2k(helium%rdf_delr, "angstrom")
WRITE (stmp, '(1X,F12.6)') rtmp
Expand Down Expand Up @@ -1378,12 +1378,12 @@ SUBROUTINE helium_print_rdf(helium_env)
CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_rdf', &
routineP = moduleN//':'//routineN

CHARACTER(len=default_string_length) :: comment, stmp
CHARACTER(len=default_string_length) :: stmp
INTEGER :: handle, ia, ic, id, itmp, iweight, k, &
nsteps, unit_nr
LOGICAL :: should_output
REAL(KIND=dp) :: inv_norm, rtmp
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: message
REAL(KIND=dp) :: inv_norm, rtmp, rtmp2
REAL(KIND=dp), DIMENSION(:, :), POINTER :: message
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: print_key

Expand All @@ -1404,38 +1404,46 @@ SUBROUTINE helium_print_rdf(helium_env)
IF (should_output) THEN
! work on the temporary array so that accumulated data remains intact
! save accumulated data of different env on same core in first temp
helium_env(1)%helium%rdf_inst(:, :, :) = 0.0_dp
helium_env(1)%helium%rdf_inst(:, :) = 0.0_dp
DO k = 1, SIZE(helium_env)
helium_env(1)%helium%rdf_inst(:, :, :) = helium_env(1)%helium%rdf_inst(:, :, :)+ &
helium_env(k)%helium%rdf_accu(:, :, :)
helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)+ &
helium_env(k)%helium%rdf_accu(:, :)
END DO

! average over processors
NULLIFY (message)
message => helium_env(1)%helium%rdf_inst(:, :, :)
message => helium_env(1)%helium%rdf_inst(:, :)
CALL mp_sum(message, helium_env(1)%comm)
itmp = helium_env(1)%helium%num_env
inv_norm = 1.0_dp/REAL(itmp, dp)
helium_env(1)%helium%rdf_inst(:, :, :) = helium_env(1)%helium%rdf_inst(:, :, :)*inv_norm
helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*inv_norm

! average over steps performed so far in this run
nsteps = helium_env(1)%helium%current_step-helium_env(1)%helium%first_step
inv_norm = 1.0_dp/REAL(nsteps, dp)
helium_env(1)%helium%rdf_inst(:, :, :) = helium_env(1)%helium%rdf_inst(:, :, :)*inv_norm
helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*inv_norm

iweight = helium_env(1)%helium%rdf_iweight
! average over the old and the current density (observe the weights!)
helium_env(1)%helium%rdf_inst(:, :, :) = nsteps*helium_env(1)%helium%rdf_inst(:, :, :)+ &
iweight*helium_env(1)%helium%rdf_rstr(:, :, :)
helium_env(1)%helium%rdf_inst(:, :, :) = helium_env(1)%helium%rdf_inst(:, :, :)/REAL(nsteps+iweight, dp)
helium_env(1)%helium%rdf_inst(:, :) = nsteps*helium_env(1)%helium%rdf_inst(:, :)+ &
iweight*helium_env(1)%helium%rdf_rstr(:, :)
helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)/REAL(nsteps+iweight, dp)

IF (logger%para_env%ionode) THEN

DO id = 1, helium_env(1)%helium%rdf_num ! loop over estimators ---
ia = 0
rtmp = cp_unit_from_cp2k(1.0_dp, "angstrom")
rtmp2 = 1.0_dp
IF (.NOT. helium_env(1)%helium%periodic) THEN
! RDF in non-periodic case has unit 1/bohr**3, convert to Angstrom:
rtmp2 = rtmp**(-3)
END IF

IF (helium_env(1)%helium%rdf_he_he) THEN
! overwrite RDF file each time it is written
ia = 1
stmp = ""
WRITE (stmp, *) id
WRITE (stmp, *) "He-He"
unit_nr = cp_print_key_unit_nr( &
logger, &
print_key, &
Expand All @@ -1444,22 +1452,39 @@ SUBROUTINE helium_print_rdf(helium_env)
file_position="REWIND", &
do_backup=.FALSE.)

comment = "Radial distribution function"
DO ic = 1, helium_env(1)%helium%rdf_nbin
WRITE (unit_nr, '(F20.10)', ADVANCE='NO') (REAL(ic, dp)-0.5_dp)*helium_env(1)%helium%rdf_delr*rtmp
WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%rdf_inst(ia, ic)*rtmp2
WRITE (unit_nr, *)
END DO

CALL m_flush(unit_nr)
CALL cp_print_key_finished_output(unit_nr, logger, print_key)
END IF

IF (helium_env(1)%helium%rdf_sol_he) THEN
! overwrite RDF file each time it is written
stmp = ""
WRITE (stmp, *) "Solute-He"
unit_nr = cp_print_key_unit_nr( &
logger, &
print_key, &
middle_name="helium-rdf-"//TRIM(ADJUSTL(stmp)), &
extension=".dat", &
file_position="REWIND", &
do_backup=.FALSE.)

DO ic = 1, helium_env(1)%helium%rdf_nbin
rtmp = (REAL(ic, dp)-0.5_dp)*helium_env(1)%helium%rdf_delr
rtmp = cp_unit_from_cp2k(rtmp, "angstrom")
WRITE (unit_nr, '(F20.10)', ADVANCE='NO') rtmp
DO ia = 1, helium_env(1)%helium%rdf_num_ctr
WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%rdf_inst(id, ic, ia)
WRITE (unit_nr, '(F20.10)', ADVANCE='NO') (REAL(ic, dp)-0.5_dp)*helium_env(1)%helium%rdf_delr*rtmp
DO id = 1+ia, helium_env(1)%helium%rdf_num
WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%rdf_inst(id, ic)*rtmp2
END DO
WRITE (unit_nr, *)
END DO

CALL m_flush(unit_nr)
CALL cp_print_key_finished_output(unit_nr, logger, print_key)

END DO ! loop over estimators ---
END IF

END IF
END IF
Expand Down

0 comments on commit 98a0cf8

Please sign in to comment.